An actuarial approach for modeling pandemic risk

This article proposes a model for pandemic risk and two stochastic extensions. It is designed for actuarial valuation of insurance plans providing healthcare and death bene(cid:28)ts. The core of our approach relies on a deterministic model that is an e(cid:30)cient alternative to the susceptible-Infected-Recovered (SIR) method. This model explains the evolution of the (cid:28)rst waves of COVID-19 in Belgium, Germany, Italy and Spain. Furthermore, it is analytically tractable for fair pure premium calculation. In a (cid:28)rst extension, we replace the time by a Gamma stochastic clock. This approach randomizes the timing of the epidemic peak. A second extension consists in adding a Brownian noise and a jump process to explain the erratic evolution of the population of con(cid:28)rmed cases. The jump component allows for local resurgences of the epidemic.


Introduction
The recent outbreak of COVID-19 reminds us that epidemics are sanitary but also nancial issues. There is a clear growing need for covering the epidemic risk and developing analytically tractable models. A classical deterministic model is the Susceptible-Infected-Recovered (SIR) model proposed by Kermack and Mc Kendrick (1927). About the works on the deterministic SIR model and its generalization, the reader can refer to Anderson and May (1979) E-mail to: donatien.hainaut(at)uclouvain.be components and the corresponding premium levels.
In this article, we propose an alternative to the SIR model in which the contagion rate is inversely proportional to time instead to the population of susceptible. This model presents a greater analytical tractability than the SIR and replicates the rst wave of COVID-19 in Belgium, Germany, Italy and Spain. In this framework, we infer a closed form expression for the fair premium rate of an insurance plan covering healthcare expenses of infected persons and providing a lump sum capital in case of death.
Two stochastic extensions are next proposed. The rst one consists to replace the time by a stochastic clock, called subordinator. When the subordinator is a Gamma process, the analytical tractability of the model is preserved and the time of the pandemic peak becomes random.
In the second stochastic extension, we add a jump-diusion component. This model allows for local resurgence of the epidemic. The premiums are identical to these of the deterministic approach, but the model generates various realistic sample paths for the number of sicks and deaths. Finally, we show that reinsurance treaties can be valued with both the time-changed and jump-diusion models.
The outline of this article is the following. We rst present the deterministic model and compare its features with the SIR approach. Section 3 studies the valuation of an insurance plan with healthcare and death benets. The model is tted to COVID-19 datasets for Belgium, Germany, Italy and Spain in the next section. Section 5 introduces the rst stochastic variant of this model, based on a random Gamma clock. The valuation of reinsurance treaties is developed in the next paragraph and the time-changed model is estimated in Section 7. Sections 8 and 9 explore the features of the jump-diusion model. Finally, we propose an estimation method and t the model to COVID-19 datasets.

A deterministic epidemic model
We consider a population of size S 0 hits by an epidemic disease at time 0. We propose to model the number of infectious persons at time t ≥ 0 , denoted by I t , by the following relation: where α, β, µ, γ ∈ R + . This function is the product of two terms. The rst one is an exponential decreasing function, e −(α+µ)t , whereas the second one is an increasing function. The number of infecteds decreases exponentially with a rate (α + µ). The parameter α is the recovery rate from the disease whereas µ is assumed to be the death rate of infected persons. The average duration before recovery is then 1/α. At time 0, the initial number of infected individuals is equal to I 0 = 0. In order to understand the role played by the parameter β in the dynamics of infectious population, we derive Equation (1): This dierential equation reveals that the initial contagion rate per capita is equal to γ t . This is a decreasing function to 0 when t → ∞. Empirical tests in following sections emphasizes that Equation (2) explains the evolution of the rst wave of COVID-19 in Belgium, Germany, Spain and Italy.
This model slightly diers from the the Susceptible-Infected-Recovered (SIR) model developed by Kermack and Mc Kendrick (1927), which is a standard in the literature. In the SIR model, the evolution of the epidemics is described by the ordinary dierential equations (ODE's): where β SIR ∈ R + and S SIR t is the number of persons that are susceptible to be infected at time t. As in our model, α SIR and µ SIR are respectively the recovery and mortality rates. The contagion rate per capita in the SIR is proportional to the population of susceptible, whereas it is a function of time, γ t , in the approach proposed in this article. This assumption allows us to obtain the closed form expression (2) for the infectious populations. This is the greatest benet of this model compared to the SIR that does not admit analytical solutions for the system (3). In actuarial applications developed in following sections, we have to integrate I t . Therefore having an analytical formula allows us to avoid numerical integration of an ODE numerical solution and propagation of numerical errors. Furthermore, the model (2) is easily extended to stochastic frameworks, as detailed in following sections.
On the other hand, the basic reproduction number, R 0 is dened as the average number of secondary cases arising from a typical primary case. Under the assumption that the population of susceptible is large, we have S SIR t S 0 ≈ 1 and the basic reproduction number of the SIR model is R 0 = β α+µ . In our model the reproduction number is instead a function of time equal to R 0 (t) = γ t(α+µ) . The time-varying R 0 allows to take into account the impact of preventive measures to curb the epidemic, as a lockdown or the wear of masks. The empirical analysis of the next section conrms that this approach oers a better t than the SIR model. Our model presents other interesting features. In particular, the peak of the epidemics is known and obtained by canceling the rst order derivative of Equation (1): Combining Equations (4) and (2) allows us to evaluate the size of the infected population when the peak is reached: Since the population of infectious persons may not exceed S 0 , we infer the necessary conditions β γ ≤ e γ (tmax) γ . As µ is the mortality rate, the total number of deaths up to time t is a function, denoted by D t , that is solution of the ordinary dierential equation: Under the assumptions that recovery does not not provide a protective immunity and that there is no entry into or departure from the population, the size of the population of susceptible, denoted by S t , is then solution of the following ODE: By construction, the population of susceptible growths when infecteds recover from the disease and is decreased by the number of new contaminated. As we do not consider new entrants in the population, the sum of the number of susceptible, infecteds and deaths remains constant and equal to S 0 : 3 Actuarial valuation of an insurance plan In the same spirit as Feng and Garrido (2011), we consider an infectious disease insurance plan that collects premiums in the form of continuous annuities from susceptible, as long as they are healthy. The premium rate is assumed constant and noted p. Collected premiums cover medical expenses which are continuously reimbursed for each infected policyholder during the period of treatment. The benet rate is noted b. The plan terminates when the individual recovers or dies from the disease. In case of death, a lump sum benet, c, is paid.
The risk-free rate is constant and denoted by r. If the insurance plan covers the whole population, premium, benet rates and lump sum death capital must ensure the nancial equilibrium of the plan. Under the assumption that the plan starts at time 0 and nishes at time T , discounted premiums have to cover all discounted benets: As stated in the previous section, the discounted integral of I t admits a closed form-expression.
Proof Using Equation (1), we develop the integral as follows: If we remember the denition of the Gamma incomplete function we directly obtain Equation (10).
From this last proposition, we immediately infer the expressions of D t , T 0 e −rs dD s and S t .
Corollary 3.2. The cumulated number of deceases caused by the epidemic at time t ≥ 0 is equal to If θ = r + α + µ, the second term on the right hand side of Equation (9) is: T 0 e −rs dD s = µS 0 β γ θ −γ−1 Γ l (γ + 1, T θ) . (12) The size of the population of susceptible at time t ≥ 0 is deduced from the relation S t +I t +D t = S 0 : The next proposition reports the closed form expression for the premium rate solution of Equation (9). Proposition 3.3. For benets rate (b, c), the fair premium rate that ensures the actuarial equilibrium of the plan is given by T 0 e −rs S s ds , (14) where the denominator is equal to: Proof From Equation (13), we infer that  (11), we infer that the integral of the discounted number of deaths is: T 0 e −rs Γ l (γ + 1, s (α + µ)) ds.

Empirical Illustration
We t the model to data about the COVID-19 outbreak in Belgium, Germany, Italy and Spain.
We use the datasets from the library coronavirus 1 in R which provides daily time-series of the number of deaths and detected cases of COVID from the beginning of 2020 up to end July. We choose as starting date the day when the number of conrmed cases passes above a threshold set to 0.005% of the total population of the country. As the model is designed for modeling a single epidemic wave, the ending date is set to the 15th of June 2020, which corresponds to the end of the lockdown period in considered countries. Figure 1 shows these time-series and Table   1 See, https://github.com/RamiKrispin/coronavirus, package developed by Rami Krispin.
1 reports some statistics of the datasets. For Spain, the number of conrmed cases or deaths is negative for a few days. This is due to retrospective corrections.  Both the SIR and our model aim to describe the evolution of the number of infected persons, As the data sets only report new conrmed cases, we assume that contaminated individuals remain infected on average during 12 days, which is slightly less than the duration of the quarantine imposed e.g. in Belgium (14 days) after being in touch with a contaminated person. If I obs t t=1,...,n obs is the time-series of observations and n obs , the number of days, parameters are obtained by a weighted least square minimization: where I k is the value of I t at time t k . As the impacts α and µ on I t are indistinguishable, we rst estimate their sum. The annualized mortality rate is estimated as the ratio of the total number of deaths on the cumulated number of infecteds forecast by the model, multiplied by 365. Given that the COVID testing was far from being generalized in March, it is likely that the number of infected cases was higher than the one reported. In order to take this into account in the estimation procedure, more importance is granted to most recent daily observations as follows: The results of the calibration procedure are reported in Table 2.  We benchmark the capacity of Equation (1) to model the evolution of the infected population with the SIR model, tted by least square minimization. Our numerical experiments reveals that the SIR model fails to replicate the curve of I t . The only way to t this model consists to consider that S 0 is also a parameter. Parameter estimates are reported in Table 3. Figure 2 compares observed to forecast I t with our and SIR models. The SIR oers at a rst sight an excellent t but considering that S 0 is adjustable is hard to justify. Furthermore, the adjusted S 0 are considerably smaller than the real size of considered populations. This conrms that our approach is a reliable alternative compared to the SIR model.  Table 3: Parameter estimates of the SIR model.
Next, we use parameter estimates in order to evaluate the fair premium rates of an insurance plan, such as described in Section 3. Two cases are considered. In the rst one, collected premiums cover exclusively medical expenses. An allowance of 1000¿ per day is paid during the treatment (c = 365 000¿ on a yearly basis). The second plan covers exclusively the death risk: a lump sum capital of 200 000¿ is paid at the decease of an infected patient. The duration of both plans is 6 month and the risk-free rate is set to 2%.     Table 5: Fair premium rates with the model (1). Duration T = 0.5 year and r = 2%.

Time-changed extension
The model introduced in Section 2 is fully deterministic. In practice, we observe random uctuations of the number of infected persons. In order to replicate such random variations, we propose two stochastic extensions of Equation (1). The rst one developed in this section, consists to replace the time t by a stochastic clock, also called subordinator. This clock is an increasing positive process denoted by (τ t ) t≥0 , dened on a probability space Ω endowed with a probability measure P and its natural ltration (F t ) t≥0 . We consider that τ t is a Gamma process: i.e. τ t is Gamma distributed with expectation and variance equal to λt. The probability density function of τ t is given by where Γ(a) = ∞ 0 x a−1 e −x dx is the standard Gamma function that is such that Γ(a+1) = aΓ(a), for a > 0. A straightforward calculation shows that the characteristic function for the gamma process is given by is the characteristic exponent. This Lévy-Khinchine representation of the characteristic function reveals that τ t is also a Lévy process. Indeed ψ(u) may be rewritten as the following integral: from which we infer that the Lévy measure of τ t is ν(dx) = λx −1 e −x 1 (x≥0) dx. τ t is a process with nite variations. Therefore, for any function f (t, τ t ) of time and of the subordinator, the Itô's lemma for semi-martingales states that Whereas the F t −expectation of its innitesimal variation is The time changed version of the deterministic model is obtained by replacing t with the chronometer τ t . The dynamics of the population of infectious individuals is then: where α, β, µ, γ ∈ R + . Using the Itô's lemma and First-order Taylor developments, we infer that: This rst order approximation emphasizes the strong relation between the deterministic and the time changed dynamics. α and µ may still be interpreted as recovery and death rates but over a random time interval of size ∆ t . The contagion rate per capita at time t is equal to γ τ t− for a period ∆ t . The next proposition shows the rst two moments of (I t ) t≥0 .
Proposition 5.1. The expected number of infected persons at time t is equal to: whereas its variance is given by the following relation: Proof. The expectation of I t is rewritten in its integral form: Next we do a change of variable v = (α + µ + 1) u in order to obtain Eq. (21). We obtain the moments of second order in a similar manner: .
Equation (22) is the dierence of this second moment and of the square of the expectation.
Notice that the maximum of the epidemics is reached at t max with τ tmax = γ α+u . The cumulated number of deaths is the time-changed version of Equation (11): We use the Itô's lemma and First-order Taylor developments to check that the dynamics of D t is compliant with the one of I t : which is the stochastic equivalent equation to (6). This equation conrms that the innitesimal variation of D t is a fraction µ∆τ t of the population of infecteds. Unfortunately, the expectation and variance of the number of deaths only admit a semi-closed form expression and their valuation requires a numerical integration: The variance of the cumulated number of deaths up to time t is therefore: The size of the population of susceptible at time t ≥ 0 is deduced from the relation S t + I t + D t = S 0 and is given by the following expression: The expected size of the population of susceptible is simply equal to (21) and (25).
If we consider the insurance plan introduced in Section 3, the fair premium rate that nance expected benets is such that Contrary to the deterministic model, the integrals present in this last expression do not admit closed-form expression but can easily be approached by a sum over a partition of the interval [0, T ]. If we consider a partition {s 0 = 0, s 1 , ...., s m = T } of equi-spaced times and if we note by ∆ m the length of inter-arrival times, the integrals are computed as the following sums: In the next Section, we present a dierent stochastic extension of our deterministic model that leads to analytical expressions for these integrals.
6 Reinsurance in the time-changed model The introduction of randomness in the dynamic of the epidemic allows us to price various reinsurance coverages. As illustration, we consider excess-of-loss contracts providing a compensation if the number of infected or the number of deaths exceeds a certain threshold. The following proposition provides an analytical expression for a reinsurance treaty that plans the payment of an amount C (I t − K) at time t, where C, K ∈ R + , if I t > K. This contract is similar to a nancial option with an underlying that is the size of the infected populations. The risk-free rate is noted r and the model parameters for I t and τ t are assumed the same under the pricing and real measures 2 .
Proposition 6.1. The value of an excess-of-loss reinsurance covering an excessive number of infecteds is equal to: where u K is the positive solution of the following equation and Γ u (γ + 1, x) = ∞ x u γ e −u du is the Gamma upper incomplete function. Proof We insert the denition of I t and rewrite the expectation as an integral with respect to the density of τ t : The integrand is positive if and only if u is above u K , such as dened in Equation (29). This allows us to rewrite the integral as a dierence where the second term is the ratio . We replace the integrating variable by v that is (α + µ + 1) u = v and infer that the rst term is: and we retrieve Equation (29).
Notice that in practice, the reinsurer will price reinsurance treaties a set of parameters more conservative than the ones tted to data in order to include a safety loading. The time-changed model can also be used for the evaluation of contracts covering an excess of mortality caused by the epidemic. To illustrate this, we consider a treaty that plans the payment of an amount C (D t − K) at time t, where C, K ∈ R + , if D t > K. The value of this contract is similar to an option written on D t . If r is the risk-free rate, u K is the positive solution of the equation: Γ l (γ + 1, u (α + µ)) = K µS 0 β γ (α + µ) −γ−1 , then the value of the reinsurance treaty is: Unfortunately, the integral in this last equation does not admit simple analytical form and must be numerically computed.

Estimation and Illustration
In order to illustrate the ability of the time-changed model to explain the evolution of a pandemic, we t it to COVID-19 data sets for Belgium, Germany, Italy and Spain. As in Section 4, parameter estimates are found by a weighted least square minimization between expected and observed sizes of infected population. We use the same weights as those in Equation (17). Given that the remission and mortality rates have the same impact on I t , we estimate their sum. The force of mortality is next found by considering the ratio deaths on number of infected persons adjusted by a time coecient. More precisely, given that τ t has increments that are Gamma distributed, E (∆τ t ) = λ dt. From Equation (24), the expectation E (dD t |F 0 ) ≈ µE (I t− |F 0 ) λ dt.
Using a moment matching approach, we then estimate the mortality rate as follows: Parameter estimates are presented in Table 6 and Figure Table 7: Fair premium rates with the model (1). Duration T = 0.5 year and r = 2%.
If we limit our analysis to a comparison of the expected number of infecteds and premium rates, we have the impression that both deterministic and time-changed models are quite similar.
However, this is far from being the case given that the second one is stochastic. In order to emphasize the dierent behaviour of this model, we simulate 2000 samples paths with parameter estimates obtained for Italy. Figure 4 shows 3 of these paths and the average over the 2000 simulations. We see that the time-changed model generates curves of I t with a very dierent shape than the average sample path. The simulated peaks of the epidemic may be far above the observed one and the timing of this peak displays a high variance. The sample paths are also much more discontinuous than the real evolution of I t . It is also interesting to look at the distribution of the cumulated number of deaths. Figure 5 presents the histograms of D t for 2000 simulations at time t = 0.1 and t = 0.25. For t = 0.1, we have a bimodal distribution whereas the number of cumulated deaths after a quarter is nearly deterministic. This is explained by the type of randomness driving the model. The stochastic clock either delays or advances the time of the epidemic peak but does not modify the ultimate total number of deaths caused by the pandemic. Table 8 reports the expectations, standard deviations, 5% and 95% percentiles of D t , computed by simulations for the other countries. We draw from those gures the same conclusions. The discontinuities of I t and the bimodal behaviour of D t being unlikely in practice (at least for the COVID-19), we investigate in the next section, an alternative stochastic extension of the deterministic model.

A jump-diusion model
This section presents an alternative stochastic model for the dynamic of the infected population that includes random noise and local resurgence of the epidemic. This model also has an excellent analytical tractability and is estimated by a peak over threshold approach. We consider a probability space (Ω, F, P ) on which are dened a Brownian motion, (W t ) t≥0 and a compound jump process (L t ) t≥0 . We denote by (N t ) t≥0 , a Poisson process with intensity λ ∈ R + , and by (J k ) k∈N ∼ J, i.i.d. random variables dened on R + with a probability density function denoted by f J (.). The expectation and variance of a jump are respectively noted µ J = E(J) and V(J) = σ 2 J . The compound Poisson process (L t ) t≥0 is dened as the sum of jumps J k up to time t: We assume that the dynamics of the population of infecteds is ruled by the following geometric jump diusion: where α, µ, γ t are respectively the recovery, mortality and contagion rates. The term σI t dW t , with σ ∈ R + , is a Gaussian noise whereas I t dL t introduces random discontinuities caused by the discovery of clusters of infection. The next proposition presents the solution of Equation (32).
Under the assumption that I t is ruled by Equation (32), the size of the population of infected people is equal to: This result is a direct consequence of the Itô's lemma for jump-diusion. In order to evaluate an insurance plan, we need to adopt a dynamic for the cumulated number of deaths, D t . As in previous models, we assume that the instantaneous number of death at time t is a proportion µdt of the population of sick persons. The dierential equation ruling S t , the population of susceptible, guarantees that the total size of the population remains equal to S 0 : The number of deaths is D t = µ t 0 I s ds whereas S t = S 0 − I t − D t . Of course we could include a random noise in the dynamics of D t but this would not fundamentally modify the results developed in the remainder of this section. The next proposition presents the rst two moments of I t . and Proof Given that W t is independent from N t , the expectation of I t is equal to a product of expectations: On the other hand, the moment generating function of N t is E e ωNt = e λt(e ω −1) and N t is independent for (J k ) k=1...,Nt . Therefore, we conclude that the expectation of the product of jumps is: Combining these elements leads to Equation (35). In a similar manner, we calculate the second order moment of I t that is Since e 2σWt is log-normal, E e 2σWt |F 0 = e 2σ 2 t and E e −σ 2 t+2σWt |F 0 = e σ 2 t . Furthermore, the expectation of the square of the product of jumps is equal to: We obtain then the second moment of I t and the variance (36).
The next proposition presents an analytical formula for evaluating expected healthcare benets paid up to time t. This result is similar to the one for the deterministic model, except that it takes into account the frequency and the average size of jumps. Proposition 8.2. Let us consider a discount rate r ∈ R + . The integral of expected discounted numbers of infected is equal to: where θ J = r + α + µ − λµ J and Γ l (γ + 1, x) is the Gamma lower incomplete function.
The proof is similar to the one of Proposition 3.1. We insert the expression (34) of E(I s |F 0 ) in the integral and obtain Equation (37). Contrary to the time-changed model, the expected number of deaths and its discounted integral admit a closed form expression in terms of incomplete lower Gamma functions.
Corollary 8.3. The expected cumulated number of deceases caused by the epidemic at time t ≥ 0 is equal to If θ J = (r + α + µ) − λµ J , the expectation of the integral of discounted variation of D t is equal to: Let us again consider the insurance plan introduced in Equation (3) with a maturity T . The premium rate is still noted p. The rate of healthcare expenses is b and the lump sum benet in case of death is c. The next proposition presents the fair premium rate that ensures the equilibrium of this plan under the assumption that the pricing and real measures are identical. Proposition 8.4. For benets rates (b, c), the fair premium rate that guarantees the actuarial equilibrium of the plan is given by: where the denominator is equal to: Proof The fair premium rate in a stochastic framework is given by Equation (28). The integral T 0 e −rs E (I s | F 0 ) ds is provided by Equation (37).
An therefore, the integral of the discounted expected number of deaths is equal to: Combining expressions (42) and (43)  9 Reinsurance in the jump-diusion model As in the time-changed model, we can evaluate various reinsurance coverages by Monte-Carlo simulations. When jumps are constant and equal to J = κ ∈ R + , reinsurance treaties with a payo dependent on I t can be valued with a closed-form expression. In this case, conditionally to N t = n jumps, I t is log-normal: where X (n) t ∼ N µ X (t, n) ; σ 2 X (t, n) . The mean and variance of X (n) t are respectively equal to µ X (t, n) = n ln(1 + κ) − α + µ + σ 2 2 t , σ 2 X (t, n) = σ 2 t .
The following proposition provides an analytical expression for a reinsurance treaty that plans the payment of an amount C (I t − K) at time t, where C, K ∈ R + , if I t > K. Model parameters are assumed the same under the pricing and real measures.
Proposition 9.1. Let us dene The value of an excess-of-loss reinsurance covering an excessive number of infecteds is equal to: where P (N t = n) = (λt) n n! e −λt and Φ(.) is the cumulative probability distribution of a standard normal random variable. Proof We can rewrite the price of this treaty as a sum of conditional expectations with respect to N t : The conditional expectations may be developed as the dierence of two integrals: The second term, after a change of variable, is equal to: In the same manner, the rst term of Equation (47) becomes: and we retrieve the result.
By construction, the cumulated number of deaths up to time t is proportional to the integral of I t . Unfortunately, the statistical distribution of D t is unknown and therefore reinsurance treaties covering excess of mortality must be valued by simulations.

Estimation of the jump diusion model
The jump-diusion model is tted in two steps. Let us recall that E (I k |F 0 ) = S 0 e ηt k (βt k ) γ where η = λµ J − (α + µ) and t k are the observation times for k = 1, ..., n obs . The rst step consists to estimate β, γ and η by minimizing the weighted sum of squares between the expected and observed numbers of infected persons: ( η, γ, β) = arg min Given that the expectation of I t in the jump-diusion approach coincides with the deterministic model, we obtain the same estimates β, γ as those in Table 2. In the second stage, we t the jump process by the peak over threshold method. From Equations (33) and (35), we dene Y t as the ratio of the process I t on its expectation: (1 + J k ) .
Using the Itô's lemma, we infer that Y t is driven by the following innitesimal dynamics: The value of Y t at time t k is noted Y k and we dene Z k = Yt . If the time lag of one day between to successive observations is noted ∆, according to Equation (49), Z k is approximately the sum where k ∼ N 0, σ √ ∆ and ∆L k = L t k − L t k−1 . A jump is believed to occur when Z k is above a threshold, noted by g(q), where q is a condence level. To dene the threshold, we t a pure Gaussian process to the time-series of Z k ∼ µ z ∆ + σ z W ∆ . The unbiased estimators of µ z and σ z are: If Φ(.) denotes the cumulative distribution function of a standard normal, g(q) is set to the q percentile of the Gaussian process: g(q) = µ z ∆ + σ z √ ∆Φ −1 (q). The time s k of the k th jump is therefore: s k = min{t j ∈ {t 1 , ..., t n } | z j ≥ g(q) , j ≥ k} . and the sample path of (N t ) t≥0 is approached by the following time series: Under the assumption that the diusion is negligible with respect to jumps, we assimilate z k −µ z ∆ to J k when a jump is detected at time t k .
We have applied this estimation procedure to COVID-19 data sets for Belgium, Germany, Italy and Spain. The estimates β and γ are found by minimization of weighted least squares as shown in Equation (48). We use the same weights as those used to t previous models (see Equation (17)). The other parameters are estimated by the peak over threshold method and are reported in Table 9. Given that the number of conrmed cases in March is underestimated due to the penury of tests, we mainly focus on the period from mid-April to mid-June to calculate the time-series of Z k . The threshold condence level is set to 90%. Figure 6 shows those series and the threshold level. For Spain, we have not take into account two abnormally high and low values due to retrospective adjustments of the number of conrmed cases. For Belgium and Spain, Z k 's have a mean close to zero and the variance seems more or less constant for dierent time windows. This conrms that the postulated dynamic in Equation (50) for Z k is acceptable for those countries. For Italy, the variance of Z k seems to raise during the month of June, but their mean is close to zero. For Germany, the Z k 's have a residual linear increasing trend and their variance seems to increase as for Italy. The consequence is that the peak-over-threshold method tends to overestimate the size and frequency of jumps. Nevertheless, this is a conservative approach from the insurer point of view and therefore we accept these parameter estimates. For the same reason, we do not consider negative jumps.

Spain
Year # z(t) Threshold Figure 6: Plot of Z k time-series for Belgium, Germany, Italy and Spain. The threshold is calculated with q = 90%.
By construction, the expectation of the jump-diusion model corresponds to the deterministic model of Section 2. Therefore, the fair premium rate of an insurance plan covering healthcare expenses and death risk are the same with both approaches. We refer the reader to the previous Section 4 for numerical examples. However, the jump-diusion model can generate a wide variety of sample paths for I t . This point is illustrated in Figure 7 that shows 1000 simulated paths and the expectation of I t over a quarter. Notice that these simulations are performed with the assumption that jumps J are constant and equal to µ J . Contrary to the time-changed model, the jump-diusion model generates smoother sample paths which all appear as likely scenarios for the evolution of the infected population. Table 10 reports statistics about the simulated number of deceases over one quarter. Globally, the mean and standard deviation of D 0.25 seems credible and the distribution of D t does not display any bimodality as in the time-changed model.  Table 9 and constant J. The thick line is the expectation of I t .   Table 10: Expectation, standard deviation and 5%-95% percentiles of simulated cumulated number of deaths, D t . 1000 simulations. Table 11 presents the prices of a few excess-of-loss reinsurances with I t and D t as underlying risks. The treaties on I t have a threshold K equal to one fourth of the total number of reported cases over the considered period of time (see Table 1). The time horizon of the contract and the capital by unit in excess are respectively set to t = 0.1 year and C = 1¿. Prices are calculated with the closed form expression (46). We use Monte-Carlo simulations for the valuation of reinsurance contracts covering an excessive mortality. The threshold K is set to the observed number of deaths observed over the considered period (see Table 1). The time horizon and the capital by unit in excess are respectively set to t = 0.1 year and C = 1¿.

Conclusions
This article develops a model for actuarial valuation of insurance plans covering claims caused by a pandemic. The basic reproduction number, R 0 , decays with time in order to replicate the impact of preventive measures to curb the epidemic. Contrary to the SIR, the empirical tests performed on COVID-19 data conrm that our model can explain the rst wave of such an epidemic. Furthermore, this model presents a high level of analytical tractability and the fair premium rate admits a closed form expression.
The time-changed extension randomizes the time of the pandemic peak by observing the epidemic on a stochastic time scale. This preserves the main features of the deterministic model and leads to comparable premium rates. Nevertheless, the simulation study reveals that sample paths of I t display a very dierent trend from what is observed for the COVID-19 outbreak.
In a second extension, we add a Brownian noise and a jump process in order to replicate local resurgences of the virus. In this setting, fair insurance premiums and excess of loss reinsurance treaties still admit closed form expressions. The model is easy to estimate by a peak over threshold approach and empirical results show that it generates realistic sample paths for I t .
This model is designed to explain the rst wave of an epidemic, with eventual local resurgence.
Even if the rst results are encouraging, they remain several open questions. For example, further research has then to be conducted to allow for a second massive wave of contagion. It would also be interesting to consider a multivariate extension in order to model cross-borders transmission of a virus. Finally, this approach fails to explain epidemic with a fat right tail for I t , such as observed in USA or Brazil for the COVID-19.