Laurent Lesage

Madalina Deaconu

Antoine Lejay

Jorge Augusto Meira

Geoffrey Nichil

Radu State

Hawkes processes framework with a Gamma density as excitation function: application to natural disasters for insurance

Keywords: Point processes, Hawkes processes, insurance, EM algorithm, natural disasters 1

Hawkes process are temporal self-exciting point processes. They are well established in earthquake modelling or finance and their application is spreading to diverse areas. Most models from the literature have two major drawbacks regarding their potential application to insurance. First, they use an exponentially-decaying form of excitation, which does not allow a delay between the occurrence of an event and its excitation effect on the process and does not fit well on insurance data consequently. Second, theoretical results developed from these models are valid only when time of observation tends to infinity, whereas the time horizon for an insurance use case is of several months or years. In this paper, we define a complete framework of Hawkes processes with a Gamma density excitation function (i.e. estimation, simulation, goodness-of-fit) instead of an exponential-decaying function and we demonstrate some mathematical properties (i.e. expectation, variance) about the transient regime of the process. We illustrate our results with real insurance data about natural disasters in Luxembourg.

The Gamma density excitation function allows us to deal with this delay and could be seen as a generalization of the exponential excitation function, since both are equivalent for a particular set of parameters. In this paper, we develop a framework for Hawkes processes with Gamma density excitation function regarding estimation and simulation. A new result in this context is the calculus of the mean and the variance of the process in the transient regime. At our best knowledge in the literature, there only exists limit results (see [START_REF] Gao | Functional central limit theorems for stationary Hawkes processes and application to infiniteserver queues[END_REF]). Such limit results established on the stationary regime are not suitable to deal with practical cases in the insurance context, where the time horizon is of several months. To illustrate our work, we develop the use case of natural disasters in Luxembourg.

Main contributions. The major contributions of this paper are to:

1. Define Hawkes processes and the elements needed for a complete use case study (estimation, simulation, goodness-of-fit, etc.);

2. Study mathematical expressions and properties of Hawkes processes (i.e. expectation, variance and central limit theorem) with a Gamma density as excitation function;

3. Model an insurance use case with real data by this specific form of Hawkes processes as an illustration of the mathematical tools.

Plan. The remainder of this paper is organized as follows. The second section defines a Hawkes process, introduces all useful notations and proposes methods to simulate the process. The third section contains the development of our results on mathematical properties of Hawkes processes. We also study a use case applied to insurance in the fourth section (natural disasters in Luxembourg).

Hawkes processes definition and simulation

In this section, we define the notion of Hawkes process, and we introduce useful notations. We also introduce a method to estimate and to simulate the Hawkes processes.

Point process

Hawkes processes are a particular class of point processes. We first propose a definition of a point process, inspired by [START_REF] Daley | An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods[END_REF].

Definition 1. (Point process).

A point process T is an increasing sequence of random variables T = {T 1 , T 2 , ...} called occurrences which takes values in [0, +∞[, such that P(0 ≤ T 1 ≤ T 2 ≤ ...) = 1.

Definition 2. (History of occurrences). We define (H(t), t ≥ 0) the filtration and the history of the occurrences until time t.

Another useful representation of a sequence of occurrences is the counting notation.

Definition 3. (Counting process). We say that N is a counting process if it is an almost surely finite stochastic process and a right-continuous step function with +1 increments after each step, taking values in N, and N (0) = 0.

The two notions of point process and counting process are interchangeable. Indeed, a counting process could be seen as the cumulative count of a point process. The link between T and N is:

∀t ∈ R + , N (t) = # T i ∈ T T i ≤ t , i ≥ 1,
where # denotes the cardinality of a set and T i is the time of the i th occurrence in the time interval [0, T ], where T is the final time of observation.

Intensity function

The most frequently used characterization of a

λ(t|H(t)) = lim h→0 + E[N (t + h) -N (t)|H(t)] h , (2.1)
then λ(t|H(t)) is the conditional intensity function of N (•). We denote λ(t|H(t)) by λ * (t) for the remainder of the paper.

The following proposition is inspired from Proposition 2.2 of [START_REF] Rasmussen | Temporal Point Processes and the Conditional Intensity Function[END_REF] and states that the conditional intensity function uniquely defines the point process. Thus in the following we use this characterization.

Proposition 1. Let us consider a conditional intensity function λ * (•), defined for any time period [u, t], 0 ≤ u ≤ t, such that:

• λ * (•) non-negative and integrable on [u, t],

• lim t→+∞ t u λ * (s)ds = +∞.

Then there exists an unique point process with λ * (t) conditional intensity function.

Hawkes process definition

We introduce now the definition of Hawkes processes.

λ * (t) = λ + t 0 µ(t -u)dN (u) = λ + N (t) k=1 µ(t -t k ), (2.2)
where N (•) is defined in Definition 3.

Hawkes processes are point processes who have a so-called "self-exciting" property. It means that each occurrence increases the intensity, according to the excitation function. A higher intensity leads to more occurrences, leading to a higher frequency, etc. A Hawkes process with a null excitation function is a homogeneous Poisson process of rate λ.

Definition 6. (Γ-Hawkes processes). We define the Γ-Hawkes processes, a Hawkes process whose excitation function is of form:

µ(t) = α t k1-1 exp(-t k2 ) k k1 2 Γ(k 1 ) , t ≥ 0, α ∈ R * + , k 1 , k 2 ∈ R * + , (2.3)
where Γ(•) is the Gamma function. k 1 is the scale parameter and k 2 is the shape parameter.

Remark 1. For k 1 = 1, we obtain the exponential excitation function:

µ(t) = α exp(-t k2 ) k 2 , t ≥ 0, α ∈ R * + , k 2 ∈ R * + , (2.4)
used for instance originally by Hawkes [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF]. In this situation:

λ * (t) = λ + N (t) k=1 α exp(-t-t k k2 ) k 2 .
(2.5)

For k 1 > 1, this excitation function allows to avoid a discontinuity in the intensity function when a event occurs, and introduce a delay between an event and its impact on the process. Indeed, in that case, µ(0) = 0. This property is very helpful when considering events for which an occurrence and the resulting self-excitation are delayed.

Simulation method

Being able to simulate Hawkes process is useful to check or estimate numerically mathematical properties of the process. We present here the method introduced by Ogata in [START_REF] Ogata | Simulation Method for Point Processes[END_REF].

The idea of Ogata's thinning algorithm is that randomly removing points from a "faster" simulated Poisson process (i.e. with a higher intensity function than targeted Hawkes process) allows to simulate a Hawkes process. This idea is based on the proposition below, which is inspired from Theorem 4.2 of [START_REF] Chen | Thinning Algorithms for Simulating Point Processes[END_REF].

Proposition 2. We consider a point process, with conditional intensity function λ * (•). We also consider an upper-bound for λ * (•), λ ≥ 0, and a Poisson process N of parameter λ of jumping times Ti , i ≥ 1. Then the point process defined by:

N (t) = i≥1 1 { Ti≤t,Ui≤ λ * ( Ti ) λ } , (2.6)
where (U i ) are i.i.d and follow an uniform distribution on [0, 1], is a process of conditional intensity function λ * (•).

Branching structure

In this subsection, we present an interpretation of Hawkes processes introduced in [START_REF] Halpin | Modelling Dyadic Interaction with Hawkes Processes[END_REF]. This approach is suitable for discussing mathematical properties and for parameters estimation. Events occurring in a Hawkes process could be seen as a population, which expands either by immigration or births. Immigrants correspond to occurrences due to the background intensity, which are spontaneous. Births correspond to occurrences due to the self-exciting part of the intensity, which are a consequence of the previous occurrences. We can consider that the Hawkes process is a sum of independent Poisson processes, where each occurrence would be a consequence of one of these Poisson processes:

λ * (t) = λ + N (t) k=1 µ(t -t k ) = N (t) i=0 λ * i (t), (2.7)
where λ * (•), N (•) and µ(•) were introduced previously.

• The number of spontaneous occurrences is a Poisson process of rate λ over [0, τ ], where τ ≥ 0 is the final time of observation. We denote λ * 0 (t) = λ;

• The number of responses to each event

t i , i ≥ 1 is a Poisson process of rate µ(t -t i ) over [t i , τ ]. We denote λ * i (t) = αµ(t -t i ).
Definition 7. (Process of response). For each event T i , we introduce the random variable Z i , which indicates the process to which the event belongs.

Figure 1: A realisation of a Hawkes process and its branching structure

As an example, we consider the realisation of a Hawkes process, represented in the Figure 1.

• Z 1 = 0: event T 1 is a spontaneous occurrence.

• Z 2 = 1: event T 2 is a response to event T 1 .

• Z 3 = 1: event T 3 is a response to event T 1 .

• Z 4 = 3: event T 4 is a response to event T 3 .

• Z 5 = 0: event T 5 is a spontaneous occurrence.

• And so on.

Estimation method: the EM algorithm

This subsection details a method to estimate parameters of Hawkes processes, which is developed in [START_REF] Halpin | Modelling Dyadic Interaction with Hawkes Processes[END_REF]. First, we define the likelihood of Hawkes processes.

Definition 8. (Observations likelihood).

Let us consider N (•) a point process of intensity function λ * (•), on time period [0, τ ], τ > 0. We denote T = {T 1 , ..., T N (τ ) } the jumping times until time τ . We observe N (τ ) = n and T = {t 1 , .., t n }. Observations log-likelihood ({t 1 , .., t n }, θ) of N (•) is:

({t 1 , .., t n }, θ) = n i=1 log λ * (t i ) - τ 0 λ * (s)ds, (2.8)
where θ is the set of parameters to estimate (i.e. λ and parameters of the excitation function µ).

Then we define the so-called complete likelihood, when assuming that {Z 1 , ..., Z n } introduced in Definition 7 are observed.

Definition 9. (Complete likelihood). Assuming we observe {Z 1 , ..., Z n } = {z 1 , ..., z n }, then complete log-likelihood is:

({t 1 , .., t n }, {z 1 , z 2 , ..., z n }, θ) = n z=1 n i=1 log λ * z (t i ) 1 zi=z -Λ z (τ ) , (2.9)
where Λ z (τ ) = τ 0 λ * z (s)ds. This complete likelihood allows us to build the Expectation-Maximization algorithm (or EM algorithm), which will estimate the parameters Hawkes processes by minimizing this likelihood. We first define the quantity Q(θ|θ ) below.

Definition 10. We define:

Q(θ|θ ) = E l({t 1 , .., t n }, {Z 1 , ..., Z n }, θ) {t 1 , .., t n }, θ , (2.10)
where θ and θ are two sets of parameters and whose computation is the E step of the EM algorithm (see Proposition 3).

We can show that:

Q(θ|θ ) = n z=1 n i=1 log λ * z (t i ) P Z i = z {t 1 , .., t n }, θ -Λ z (τ ) . (2.11)
This equation (2.11) will be useful to describe the EM algorithm. This algorithm is based on the following proposition, which states that the iterative maximization of Q(θ|θ ) allows to converge to a local maximum of the complete likelihood. In other words, maximizing Q(θ|θ ) allows to estimate parameters (see Section 1.5 of [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF] for the proof).

Proposition 3. The sequence defined by

θ (m+1) = argmax θ Q(θ|θ (m)
) converges to a local maximum of ({t 1 , .., t n }, θ), the observations likelihood.

Then we write the EM algorithm:

1. Initialize m = 0, θ (m) = θ 0 2.
While convergence:

•

[E step] Compute P Z i = z {t 1 , .., t n }, θ (m) • [M step] Compute θ (m+1) = argmax θ Q(θ|θ (m)
) → maximization by Nelder-Mead method (see [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF])

• m = m + 1.

Mathematical properties

In this section, we are going to explicit some mathematical properties of the Γ-Hawkes processes (defined by (2.3)): we calculate the expectation and the variance for any time t ≥ 0 and we state a central limit theorem.

Assumption 1. For the whole Section 3, we consider a Γ-Hawkes process of conditional intensity function:

λ * (t) = λ + t 0 µ(t -u)dN (u) = λ + α N (t) k=1 β k1 (t -t k ) k1-1 exp(-β(t -t k )) Γ(k 1 ) , (3.1)
where

α ∈ R * + , β, k 1 ∈ R * + , N (t)
≥ 0 (number of events occurred between 0 and t), t k , k ∈ {1, ..., N (t)} for which we are going to explicit results for k 1 = 1 and k 1 = 2.

This study aims to understand the behaviour of the Hawkes processes in function of parameters introduced previously.

Behaviour in function of values of α

In this section, we are going to see that the Hawkes process dynamics depend on values of α (defined in Assumption 1). To do so, let us focus on g(t) = E[λ * (t)], which will be used to calculate the expectation of the Hawkes process later.

We could show that g follows a renewal equation:

g(t) = λ + t 0 µ(t -u)g(u)du. (3.2)
According to [START_REF] Asmussen | Applied Probability and Queues[END_REF], the behaviour of g(t), when t goes to infinity, depends on the value of:

+∞ 0 µ(u)du = +∞ 0 α β k1 u k1-1 exp(-βu) Γ(k 1 ) du = α. (3.3)
We distinguish, in the classic literature about point processes:

• the defective case: α < 1;

• the proper case: α = 1;

• the excessive case: α > 1.

The Hawkes process intensity goes to infinity in the excessive case, leading to an explosion of the counting process. This situation does not suit to any real event we would like to model. The proper case is not studied here, because it does not correspond to a . We concentrate our attention in the defective case. In this case, g admits a finite limit in +∞, which will be calculated later.

Assumption 2. For the rest of the Section 3, we consider a Hawkes process of conditional intensity function described in Assumption 1, with 0 < α < 1 (defective case).

Expectation: E(N (t))

The first indicator that seems relevant for our study was the evaluation of the average number of events through time. We will develop the calculations for two values of the shape parameter: k 1 = 1 (exponential excitation function) and k 1 = 2. Since a Hawkes process is the sum of a background part of intensity λ and an self-exciting part, we expect that the expectation is higher than for a Poisson process of parameter λ, which is λt for all time t > 0.

Proposition 4. The expectation of the Γ-Hawkes process (i.e. conditional intensity function given by Equation (3.1)) is:

1. For k 1 = 1 (i.e. µ(t) = α exp(-t k 2 ) k2
, exponential decay):

E(N (t)) = λ 1 -α t - αλk 2 (1 -α) 2 1 -exp - 1 -α k 2 t . (3.4) 2. For k 1 = 2 (i.e. µ(t) = α t exp(-t k 2 ) k 2 2
):

E(N (t)) = λ 1 -α t + λ √ αk 2 2(1 + √ α) 2 1 -exp - 1 + √ α k 2 t - λ √ αk 2 2(1 - √ α) 2 1 -exp - 1 - √ α k 2 t .
(3.5)

See Appendix 1 for the proof.

Let us discuss the trend of the expectation for both values of k 1 . We first notice that for α → 0 (null excitation function), we obtain E(N (t)) -→ α→0 λt. This result is natural since we already mentioned that a

Hawkes process with a null excitation function is equivalent to a Poisson process of parameter λ. For the both values of the shape parameter, the expectation reveals three regimes for the Γ-Hawkes process. As an example, we represent the expectation for k 1 = 1, k 2 = 9, α = 0.9, λ = 0.3, which follows equation (3.4), in Figure 2. It contains four plots:

• Top-left: expectation from equation (3.4) in blue and the so-called stationary regime t → λ 1-α t in red, from t = 0 to 3000;

• Bottom-left: the relative difference between the expectation and the stationary regime on this time period;

• Top-right: focus on the expectation from t = 0 to 500;

• Bottom-right: plot of the exponential term from (3.4), t → exp -1-α k2 t . We distinguish three regimes in this example:

• The transient regime (from t = 0 to 500 approximately): it is the time period necessary for the exponential term in (3.4) to become null;

• The intermediary regime (from t = 500 to 2000 approximately): the exponential term is null but the constant term could not be neglected before the stationary regime. On this period, the relative difference between the expectation and the stationary regime is above 5%;

• The stationary regime: for any value of the parameter scale, and for any Hawkes process (not only for Γ-Hawkes processes), E(N (t)) ∼ t→+∞ λ 1-α t (see [START_REF] Gao | Functional central limit theorems for stationary Hawkes processes and application to infiniteserver queues[END_REF]). The constant term in (3.4) becomes negligible before the linear term. After the first two regimes, the average number of occurrences of the Hawkes process is the same as for a Poisson process of parameter λ 1-α . We conclude in this discussion that the expectation could not be approximated correctly by the stationary regime from t = 0 to 2000. We will see that this point will have an impact on our use case study.

Variance: V(N (t))

We study now the variance of the Hawkes process, in order to see how the occurrences are spread out from the expectation calculated previously.

Proposition 5. The variance of the Γ-Hawkes process (i.e. the conditional intensity function is given by equation (3.1)) is:

1. For k 1 = 1 (i.e. µ(t) = α exp(-t k 2 ) k2
):

V(N (t)) = λ (1 -α) 3 t - λα(2 -α)k 2 (1 -α) 4 1 -exp - 1 -α k 2 t . (3.6) 2. For k 1 = 2 (i.e. µ(t) = α t exp(-t k 2 ) k 2
2

):

V(N (t)) = C 1 + C 2 t + C 3 exp(-ω 1 t) + C 4 exp(-ω 2 t
), where:

)

ω 1 = 1 + √ α k 2 , (3.8
)

ω 2 = 1 - √ α k 2 , (3.9)
C 1 = -αk 2 λ[8 -5α + α 2 ] 2(1 -α) 4 , (3.10)
C 2 = λ (1 -α) 3 , (3.11)
C 3 = - √ αλk 2 [4 -3α - √ αα] 4(1 -α) 2 (1 + √ α) 2 , (3.12
)

C 4 = √ αλk 2 [4 -3α + √ αα] 4(1 -α) 2 (1 - √ α) 2 . (3.13)
See Appendix 2 for the proof.

Let us observe the behaviour of the variance. We first notice that for α → 0 (null excitation function), we obtain V(N (t)) -→ α→0 λt. This result is natural since we already mentioned that a Hawkes process with a null excitation function is equivalent to a Poisson process of parameter λ. Secondly, we remark that

V(N (t)) ∼ t→+∞ λ (1-α) 3 t.
We have seen previously that E(N (t)) ∼ t→+∞ λ 1-α t: it is natural that we obtain as a limit a variance greater than λ 1-α t since a Hawkes process is more volatile than a Poisson process. The limit of the variance is calculated in Proposition 3 from [START_REF] Gao | Functional central limit theorems for stationary Hawkes processes and application to infiniteserver queues[END_REF], which gives the same result. It is also possible to do an analysis in terms of regimes similar to the expectation.

Central limit theorem for Hawkes processes

From Theorem 2 of [START_REF] Bacry | Scaling limits for Hawkes processes and application to financial statistics[END_REF], we state a central limit theorem for Hawkes process, when t → +∞. The following proposition is true for any k 1 ≥ 1: Proposition 6. (Central limit theorem for Hawkes processes)

lim t→+∞ N (t) -λt √ t d = σB(t), (3.14)
where:

• λ = λ 1-α ; • σ 2 = λ (1-α) 3 ;
• B(•) is a standard Brownian motion;

• d = means equality in distribution.
We recognize in this central limit theorem the limits of expectation and variance calculated previously. The central limit theorem gives the distribution of the number of events at a large time horizon.

Use case: natural disasters in Luxembourg

In this section, we present an insurance use case. Natural disasters are one of the types of rare events who could occur in any insurance setting. We apply it in the case of Luxembourg, with data from insurance. They are mainly storms and gusts of wind, whose frequency seems to increase for these last five years. This trend may be explained by climate imbalance (see public data about climate change in Luxembourg [22]).

After introducing the data provided by Foyer Assurances in Section 4.1, we present the parameters estimated on these data in Section 4.2. Then we justify the choice of using a Γ-Hawkes process and its study on a transient regime by first comparing the goodness-of-fit of our approach with several models by a statistical test, then simulating data and observing the number of events from these simulations to check whether our Γ-Hawkes process fits with the data.

Presentation of the data

Foyer Assurances has provided data about natural disasters and home insurance. Figure 3 shows the number of home insurance claims per day processed by Foyer Assurances for which the natural disasters cover was used, from January 2015 to December 2019. This timeline shows the low frequency and the high severity of this kind of event. The peaks correspond to specific weather conditions:

• 31/03/15: "Niklas" storm;

• 16/09/15: "Henri" ex-tropical storm;

• 09/02/16: gusts of wind in East France and Luxembourg;

• 13/01/17: "Egon" cyclone: storm, snow;

• 03/01/18: "Eleanor" storm;

• 10/03/19: gusts of wind in East France and Luxembourg;

• 09/08/19: exceptional tornado in Luxembourg.

Parameters estimation

To model the type of event introduced previously, we propose a Γ-Hawkes process:

λ * N D (t) = λ N D + n k=1 α N D (t -t k ) k 1,N D -1 Γ(k 1,N D )k k 1,N D 2,N D exp - (t -t k ) k 2,N D , (4.1)
where ND is an abbreviation for natural disasters and with k N D 1 = 2. This choice of model could be justified a priori as follows:

• Self-exciting property: a unique background event which occurs rarely (e.g. a storm) and is taken into account by the background intensity in the model, triggers a sequence of numerous claims around the country. The more claims there is, the more serious the event is, the more likely new claims are: this dynamic is modelled by the self-exciting part of the intensity. Thus, we are expecting that the exciting part of the process has more importance than the background part, which should lead to a α N D value greater than λ N D .

• Gamma density intensity function of parameter k N D 1 = 2: there is a delay between the occurrence of claims on the different Luxembourgish cities. The storm reaches the Luxembourg regions at different times. This function allows us to model this delay.

EM algorithm estimation leads to the following parameters in Table 1. We used the R library hawkes, developed by the authors of [START_REF] Halpin | Modelling Dyadic Interaction with Hawkes Processes[END_REF], to implement the estimation.

Parameter Estimated value λ N D 0.26 α N D 0.92 k N D 2

Table 1: Estimated values of the HP parameters by the EM algorithm

We see that α N D > λ N D , as expected, and that α N D is close to 1. Since the expectation is proportional to Figure 4 presents the evolution of the cumulative number of events (in blue) and the estimated intensity of the Γ-Hawkes process (in red). We see that intensity peaks correspond to the meteorological events described previously.

Goodness-of-fit and comparison with other models

After the parameters estimation, we need to check if this model makes sense. The following proposition from [START_REF] Daley | An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods[END_REF], called residual analysis, allows us to test the quality of the model. Thanks to this proposition, we could test with many procedures whether our data fit with a Hawkes process. The previous proposition shows us that testing whether {t 1 , t 2 , ...} follows a Hawkes process with intensity function λ is equivalent to test whether {Λ(t 1 ), Λ(t 2 ), ...} form a Process process of parameter 1. It is also equivalent to test whether every interarrival time {Λ(t 1 ), Λ(t 2 ) -Λ(t 1 ), Λ(t 3 ) -Λ(t 2 ), ...} are independent and identically distributed and follow an Exponential law of parameter 1. This could be done by:

• A Quantile-Quantile plot (QQ-plot), which represents quantiles of both distributions;

• Performing a Kolmogorov-Smirnov test, whose statistic D is the maximum absolute difference between the two distributions.

Figure 5 is a QQ-plot which allows us to compare the distribution of our data and the distribution of the Γ-Hawkes process with the estimated parameters. We also perform a Kolmogorov-Smirnov test. The p-value of the statistical test is greater than 5%, which means that we cannot reject the hypothesis that the data follow the estimated Γ-Hawkes process.

We compare the results with three other models:

• An homogeneous Poisson process, since it is the simplest point process and is equivalent to a Hawkes process with no self-excitation, whose rate is the average number of events per day;

• A Hawkes process with an exponential excitation function, in order to check whether another excitation function would have better performances than the Gamma density, whose parameters are estimated by the EM algorithm:

Parameter Estimated value λ N D 0.29 α N D 0.91 k N D 2 9.03
Table 2: Estimated value of the HP parameters by the EM algorithm for an exponential excitation function • We also evaluated a time series approach, but it does not allow to catch the severe peaks which characterize our data. This kind of approach is more efficient on data with seasonality and a regular trend.

The QQ-plots for the two models are represented on Figure 6: As expected, an homogeneous Poisson process is not enough to model this type of event. We need to take into account the self-excitation property of the process. We could also see that the fitting is better with a Gamma excitation function than an exponential decay, which validates our study a posteriori.

Simulations

Moreover, we must make sure that our model is able to replicate severe events like those observed in our data. We verify that simulated trajectories present the same structure than the actual time series: few severe events between very calm periods. Figure 7 shows the simulation of 20 trajectories over five years according to our estimated Hawkes process. 1-α N D , we slighty change λ N D and α N D so that stationary expectation remains the same. Moreover, as we would like to have larger jumps in the trajectory, it means that the importance of the self-exciting part of the process should be higher. Thus, we will increase the value of α N D . We propose the same type of simulation over five years for different sets of (λ N D , α N D ) in Figure 8. Results show that increasing α N D allows indeed the process to be more volatile and with more severe events. Simulated trajectories are more alike actual data for the different sets of parameters. The higher is α N D , the larger are the jumps. In order to quantify the capacity to generate scenarios with severe events, we evaluate by a Monte-Carlo method the average number of times the process reaches a threshold of 500 events in a period of time of 1825 days (five years), by performing 1000 simulations. We present the results on Figure 9.

Parameters (λ N D , α N D ) Value (0.19,0.94)
1.75 (0.17,0.95) 2.08 (0.13,0.96) 2.28 (0.07,0.98) 2.37

Figure 9: Average number of times the process reaches a threshold of 500 events in a period of time of 1825 days

The Monte-Carlo estimation shows that the higher is α, the higher is the average value of severe events. Therefore it confirms that another set of parameters could be more appropriate in terms of ability to generate disasters observed on real data.

In terms of goodness-of-fit, Figure 10 shows QQ-plot for each set of parameters. Figure 10 illustrates that goodness-of-fit is worse for each set of parameters, when comparing to the first set of parameters we estimated by our approach (see Figure 5). Indeed, for each set of parameters we could reject the hypothesis that the Hawkes process fits the data. We eventually have a set of parameters which maximizes the goodness-of-fit (the one estimated by the EM algorithm), but another which seems to catch better the dynamics of the time series (a set with a higher α N D ). It is explained by the fact that the goodness-of-fit is evaluated on the entire distribution, by computing the distance on all the quantiles of the distribution. Whereas the dynamics of the natural disasters rely on the accumulation of grouped events, which are the smallest quantiles of the distribution. By the way we could see on Figure 10 that the distribution gets worse on medium quantiles. There is an imbalance between the quality of simulations and the statistical tests, which illustrates the difficulty to estimate the right set of parameters for a Hawkes process.

Distribution of the number of events

In this subsection, we observe the distribution of N (t), where t = five years. We perform 1000 simulations of the Hawkes processes, with the first set of parameters estimated, for a time period of five years. We observe the errors distribution, N (t)-λt √ t

, and we compare with a Normal distribution of variance σ 2 = λ (1-α) 3 t. We see that the mean of errors is different from zero. The distribution of errors does not follow the Normal distribution, as the process is not observed long enough to be close to its limit. This observation justifies our mathematical study a posteriori, since the study of limits is not enough. We need the transitory values of the expectation and the variance, even for a long five years period.

Numerically, from equations (3.5) 4.3), we could see that we are still in either the transient state or the intermediary state for t = 5 years = 1825 days: while exponential terms are close to zero for both expectation and variance, the difference between the actual value (third column of Table 3) and the limit value (fourth column) is explained by the constant terms. The order of magnitude of these constant terms shows that the limit values are not a good approximation for a period of several years. That justifies our study on the different regimes.

Conclusion and future work

In this paper, we introduced a complete framework for the study of Hawkes process: definition, estimation and simulation. We presented a specific form of Hawkes processes, the Γ-Hawkes process, and we developed mathematical properties for this category of process. We applied this model to an insurance use case, i.e. natural disasters, and we observed that the considered Hawkes process is well fitted with historical data.

For future work, we intend to apply this model to other insurance use cases, such that:

• Life events prediction: births, marriage, job change, etc. We could use a multi-variate version of Hawkes processes (see [START_REF] Laub | [END_REF]), since one type of event could trigger another (e.g. a marriage could lead to a birth);

• Epidemic: we already studied the dynamics of Covid-19 pandemic in [START_REF] Lesage | process to make aware people of the severity of COVID-19 outbreak: application to cases in France[END_REF]. We intend to adapt the study to a long class of diseases;

• Workload prediction: we would like to anticipate peaks of workload for Foyer Assurances employees in order to optimize staffing.

6 Appendix

Appendix 1

Proof of the Proposition 4.

1. Laplace transform of g

(t) = E[λ * (t)] Since E(N (t)) = t 0 E[λ * (t)
]dt, we first focus on E[λ * (t)], that we denote g(t). We saw that:

g(t) = λ + t 0 µ(t -u)g(u)du. (6.1)
We recognize a convolution product of g and µ, denoted g * µ. Therefore it is much easier to work into Laplace domain, since L g * µ(t) (s) = L g(t) (s) × L µ(t) (s). By transforming (6.1) into Laplace domain:

L g(t) (s) = L λ + t 0 µ(t -u)g(u)du (s) = L λ (s) + L µ(t) (s) × L g(t) (s).
By rearranging, we obtain:

L g(t) (s) = L λ (s) 1 -L µ(t) (s) . (6.2)

Inverse of g(t) Laplace transform

We turn the Laplace transform into the temporal domain. In our study, µ is a Gamma density function and is given by equation (2.3). We can show that the Laplace transform is (see [START_REF] D'azzo | Linear Control Systems Analysis and Design[END_REF]):

L µ(t) (s) = α (1 + k 2 s) k1 . ( 6.3)
By using (6.2) in (6.3), we obtain:

L g(t) (s) = L λ (s) 1 -L µ(t) (s) = λ s 1 - α (1+k2s) k 1 . • For k 1 = 1 (see (2.5)): L g(t) (s) = λ s 1 -α (1+k2s) = λ 1-α s - αλ 1-α 1-α k2 + s . ( 6.4)
By transforming (6.4) into the temporal domain, we obtain:

g(t) = E(λ(t)) = λ 1 -α - αλ 1 -α exp - 1 -α k 2 t . • For k 1 = 2 (see 2.4): L g(t) (s) = λ k 2 2 s 1+ √ α k2 + s 1- √ α k2 + s + 2λ k 2 1+ √ α k2 + s 1- √ α k2 + s + λs 1+ √ α k2 + s 1- √ α k2 + s . (6.5)
By turning (6.5) into the temporal domain, we obtain:

g(t) = E(λ(t)) = λ 1 -α + λ √ α 2(1 + √ α) exp - 1 + √ α k 2 t - λ √ α 2(1 - √ α) exp - 1 - √ α k 2 t .
3. Calculation of E(N (t)) from g(t)

From g(t), we calculate E(N (t)):

E(N (t)) = t 0 g(t)dt,
which leads us to (3.4) and (3.5).

Appendix 2

Proof of the Proposition 5, inspired from Appendix A.2 of [START_REF] Laub | [END_REF].

1. Link between the variance and the covariance density Φ(τ )

We calculate the variance from the so-called covariance density Φ(τ ), τ > 0, defined in [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF]. According to [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF], the covariance density is given by the following equation:

Φ(τ ) = λµ(τ ) + τ -∞ µ(τ -u)Φ(u)du = λµ(τ ) + +∞ 0 µ(τ + u)Φ(u)du + τ 0 µ(τ -u)Φ(u)du, (6.6)
where λ = λ 1-α . The variance V(N (t)) and the covariance density are linked by the following equation: )

V(N (
dτ = α k 2 2 (s + 1 k2 ) 1 (s + 1 k2 ) L Φ(τ ) ( 1 k 2 ) -L Φ(τ ) ( 1 k 2 ) .
By replacing l 1 (s) in (6.8) we obtain:

L Φ(τ ) (s) = λL µ(τ ) (s) + α k 2 2 (s + 1 k2 ) 1 (s + 1 k2 ) L Φ(τ ) ( 1 k 2 ) -L Φ(τ ) ( 1 k 2 )
+ L Φ(τ ) (s) × L µ(τ ) (s). (6.9)

To solve (6.9) and find Φ, we need to calculate L Φ(τ ) ( 1 k2 ) and L Φ(τ ) ( 1 k2 ).

First, by substituting s = 1 k2 in (6.9), rearranging and using (6.3) for L µ(τ ) (s), we obtain:

L Φ(τ ) ( 1 k 2 ) + (2 -α)k 2 α L Φ(τ ) ( 1 k 2 ) = λk 2 2 . (6.10)
Second, by differentiating (6.9), we obtain:

L Φ(τ ) (s) = λL µ(τ ) (s) - 2αL Φ(τ ) ( 1 k2 ) k 2 2 (s + 1 k2 ) 3 + αL Φ(τ ) ( 1 k2 ) k 2
2 (s + 1 k2 ) 2 + L µ(τ ) (s)L Φ(τ ) (s) + L µ(τ ) (s)L Φ(τ ) (s). (6.11) By substituting s = 1 k2 in (6.11), rearranging and using (6.3) for L µ(τ ) (s) and L µ(τ ) (s), we obtain:

L Φ(τ ) ( 1 k 2 ) + αk 2 (2 -α) L Φ(τ ) ( 1 k 2 ) = -λk 2 α 2(2 -α) .
(6.12) (6.10) and (6.12) form a system of two equations with two variables, L Φ(τ ) ( 1 k2 ) and L Φ(τ ) ( 1 k2 ). Solving this system leads to:

L Φ(τ ) ( 1 k 2 ) = λα 4(1 -α) , L Φ(τ ) ( 1 k 2 ) = -λαk 2 4(1 -α) .
We can now solve (6.9). After rearranging:

L Φ(τ ) (s) = α[ λ + L Φ(τ ) ( 1 k2 ) + √ α k2 L Φ(τ ) ( 1 k2 )] k 2 2 (s + 1+ √ α k2 )(s + 1- √ α k2 ) - αL Φ(τ ) ( 1 k2 ) k 2 2 (s + 1- √ α k2 ) = - √ α[ λ + L Φ(τ ) ( 1 k2 ) + √ α k2 L Φ(τ ) ( 1 k2 )] 2k 2 (s + 1+ √ α k2 ) + √ α[ λ + L Φ(τ ) ( 1 k2 ) - √ α k2 L Φ(τ ) ( 1 k2 )] 2k 2 (s + 1- √ α k2 )
.

Inverse Laplace transform of Φ(τ )

We turn the Laplace transform of Φ(τ ) into the temporal domain.

For k 1 = 2, we can show that: Φ(τ ) = D 1 exp(-ω 1 τ ) + D 2 exp(-ω 2 τ ), (6.13) where:

ω 1 = 1 + √ α k 2 , (6.14
)

ω 2 = 1 - √ α k 2 , (6.15)
D 1 = - √ α[ λ + L Φ(τ ) ( 1 k2 ) + √ α k2 L Φ(τ ) ( 1 k2 )] 2k 2 = - √ α λ[4 -3α - √ αα] 8k 2 (1 -α) 2
, (6.16)

D 2 = √ α[ λ + L Φ(τ ) ( 1 k2 ) - √ α k2 L Φ(τ ) ( 1 k2 )] 2k 2 = √ α λ[4 -3α + √ αα] 8k 2 (1 -α) 2
. (6.17 4. Calculation of V(N (t)) from Φ(τ )

From Equations (6.7) and (6.13), we calculate V(N (t)).

Figure 2 :

 2 Figure 2: Expectation for k 1 = 1, k 2 = 9, α = 0.9, λ = 0.3

Figure 3 :

 3 Figure 3: Number of claims per day processed by Foyer Assurances and labelled as natural disasters since 2015

1 1 -

 1 α N D (we recall that E(N (t)) ∼ t→+∞ λ N D 1-α N D t), this value of α N D should lead to a high average number of events.

Figure 4 :

 4 Figure 4: Counting process and intensity since 2015 / Cumulated number of natural disasters events in Luxembourg since 2015 (in blue) and the intensity of the estimated Hawkes process (in red)

Proposition 7 .

 7 (Residual analysis) We consider a sequence of occurrence times {t 1 , t 2 , ...} and a monotonic, continuous compensator Λ(t) = t 0 λ * (s)ds such that lim t→+∞ Λ(t) = +∞ almost surely. The sequence {Λ(t 1 ), Λ(t 2 ), ...} is a Poisson process with an unit rate if and only if {t 1 , t 2 , ...} is a realisation from the point process defined by Λ(•).

Figure 5 :

 5 Figure 5: QQ-plot and Kolmogorov-Smirnov test for set of parameters estimated by EM algorithm

Figure 6 :

 6 Figure 6: QQ-plot and Kolmogorov-Smirnov test for an homogeneous Poisson process (left) and a Hawkes process with an exponential excitation function (right)

Figure 7 :

 7 Figure 7: Simulation of 20 trajectories over five years for set of parameters estimated by EM algorithm

Figure 8 :

 8 Figure 8: Simulation of 20 trajectories over five years for four set of parameters

Figure 10 :

 10 Figure 10: QQ-plot and Kolmogorov-Smirnov test for four set of parameters

Figure 11 :

 11 Figure 11: Distribution of errors at the end of the five years

Φ(t 2 -t 1 )dt 1 dt 2 . ( 6 . 7 ) 2 . 2 Φ 2 Φ 2 Φ

 212672222 Laplace transform of Φ(τ )We first calculate Φ(τ ). We recognize a convolution product in equation(6.6). Thus we work in the Laplace domain:L Φ(τ ) (s) = L λµ(τ ) + +∞ 0 µ(τ + u)Φ(u)du + τ 0 µ(τ -u)Φ(u)du (s) = λL µ(τ ) (s) + L +∞ 0 µ(τ + u)Φ(u)du (s) =:l1(s) +L Φ(τ ) (s) × L µ(τ ) + u) exp(-τ +u k2 ) k 2 (u)du (s) = L +∞ 0 ατ exp(-τ k2 ) exp(-u k2 ) k 2 (u)du (s) + L +∞ 0 αu exp(-τ k2 ) exp(-u k2 ) k 2 (u)du (s).By using the definition of the Laplace transform:

  )

  Definition 5. (Hawkes process). Let us consider λ > 0, the background intensity, and µ : [0, +∞[→ [0, +∞[, the excitation function. We denote {t 1 , .., t N (t) } the sequence of past occurrences until time t. A point process is a Hawkes process if its conditional intensity function is of the form:

Table 3 :

 3 Observed, calculated and limit value of expectation and varianceThis table confirms that calculated expectation and variance give the correct values and that focusing on limits is not sufficient to study a process over several years, which is a quite long period for insurance use cases. Thanks to numeric values provided by equations (4.2) and (

	and (3.7) we obtain:
	E(N (t)) = -655.55 +3.25 t -0.28 exp(-0.22 t) + 655.84 exp(-0.004 t),	(4.2)
	V(N (t)) = -82 154.94 +383.66 t -6.01 exp(-0.22 t) + 82 160.95 exp(-0.004 t),	(4.3)
	with t in days.
	The following table presents the observed, calculated and limit value of expectation and variance at the final
	time of simulations.
	Quantity	Observed	Calculated	Limit
		value	value	value: lin-
				ear term
	Expectation	5265.2	5275.8	5931.2
	Variance	621988.3	618039.4	700177.6

Acknowledgements

We developed this work through a strong collaboration with Foyer Assurances, leader of individual and professional insurance in Luxembourg, which provided the domain specific knowledge and use cases. The first author also gratefully acknowledges the funding received towards our project (number 13659700) from the Luxembourgish National Research Fund (FNR).