Next Article in Journal
Machine Learning in P&C Insurance: A Review for Pricing and Reserving
Previous Article in Journal
Are Investors’ Attention and Uncertainty Aversion the Risk Factors for Stock Markets? International Evidence from the COVID-19 Crisis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Actuarial Approach for Modeling Pandemic Risk

LIDAM Institute of Statistics, Biostatistics and Actuarial Sciences, Université Catholique de Louvain, 5032 Isnes, Belgium
Current address: 20 voie du Roman Pays, 1348 Louvain-La-Neuve, Belgium.
Submission received: 23 November 2020 / Revised: 7 December 2020 / Accepted: 11 December 2020 / Published: 23 December 2020

Abstract

:
In this article, a model for pandemic risk and two stochastic extensions is proposed. It is designed for actuarial valuation of insurance plans providing healthcare and death benefits. The core of our approach relies on a deterministic model that is an efficient alternative to the susceptible-infected-recovered (SIR) method. This model explains the evolution of the first waves of COVID-19 in Belgium, Germany, Italy and Spain. Furthermore, it is analytically tractable for fair pure premium calculation. In a first extension, we replace the time by a gamma stochastic clock. This approach randomizes the timing of the epidemic peak. A second extension consists of adding a Brownian noise and a jump process to explain the erratic evolution of the population of confirmed cases. The jump component allows for local resurgences of the epidemic.

1. Introduction

The recent outbreak of COVID-19 reminds us that epidemics raise not only sanitary but also financial issues. There is a clear and growing need for covering the epidemic risk and for developing analytically tractable models. In this article, we propose a new deterministic model in which the contagion rate is inversely proportional to time instead of to the susceptible population. This model presents a great analytical tractability and replicates the first 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 payment in case of death.
The deterministic model is next extended in two different directions. In the first one, the time scale is random and ruled by a process called subordinator. Using a gamma process as stochastic clock preserves the analytical tractability and randomizes the dynamic of the pandemic. In the second extension, the evolution of the infectious population is noised by a Brownian diffusion and a jump process.
Multiple contributions are presented in this article. The literature on mathematical modeling of epidemic is abundant, but most of the existing solutions, as compartment models, do not provide any analytical expression for the evolution of populations in each compartment. The valuation of an epidemic-linked insurance requires one to calculate integrals of infected and susceptible population sizes and is therefore computationally intensive. The three models proposed in this article do suffer from this problem and have a high level of analytical tractability for actuarial applications. Furthermore, estimating their parameters does not pose any problem, and their empirical explanatory power is comparable to the one of the susceptible-infected-recovered (SIR) approach. Finally, the jump diffusion extension allows one to simulate realistic random epidemic scenarios and to value reinsurance treaties.
We first present an overview of previous research. This is followed by the introduction of the deterministic model that is compared to the susceptible-infected-recovered (SIR) approach. Section 4 studies the valuation of an insurance plan with healthcare and death benefits. The model is fitted to COVID-19 datasets for Belgium, Germany, Italy and Spain in the next section. Section 6 introduces the first 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 8. Section 9 and Section 10 explore the features of the jump diffusion model. Finally, we propose an estimation method and fit the model to COVID-19 datasets.

2. Overview of Previous Research

Communicable diseases have always been an important part of human history as underlined by Smith (2017), who reviews their nature and proposes a brief history of pandemics. He also introduces the mathematical models used to evaluate the risk that pandemics pose to human populations. Another detailed survey of these quantitative models, including future perspective, is available in Brauer (2017).
The reference model in epidemiology is the susceptible-infected-recovered (SIR) model proposed by Kermack and McKendrick (1927). Various disease outbreaks, including the SARS epidemic of 2002–2003, the concern about a possible H5N1 influenza epidemic in 2005, the H1N1 influenza pandemic of 2009 and the Ebola outbreak of 2014 have reignited interest in epidemic models, beginning with the reformulation of the Kermack–McKendrick model by Diekmann et al. (1995). In the SIR, the population is homogeneous, whereas in reality, an epidemic spreads differently according to factors like the age or susceptibility to infection. It is then necessary to follow the secondary infections in the subpopulations separately; this is done through the next-generation matrix as explained by Diekmann et al. (1990) and Van den Driessche and Watmough (2002). Anderson and May (1979) extended the SIR model by considering the host population as a dynamic variable rather than a constant. Tchuenche et al. (2007) studied the stability of a SIR model with a time delay in the contagion dynamic. Under the assumption that all individuals are susceptible, they showed that the endemic equilibrium is stable. Similarly, Zhang and Wang (2013) studied a nonautonomous SIRS epidemic model with time delay.
The SIR model has been extended to multiple compartments with labels such as M, S, E, I and R that are often used for the epidemiological classes. The class M contains infants with passive immunity inherited at birth. After that, the infant moves to the susceptible class S. When a susceptible individual has adequate contact with an infected individual such that transmission occurs, then the susceptible individual enters the exposed class E of those in the latent period, who are infected but not yet infectious. After the latent period ends, the individual enters the class I of infectious. We refer the reader to Hethcote (2000) for a comparison of MSEIR and SEIR models for various diseases, including measles in Niger and pertussis in the United States.
Drawing conclusions from mathematical models raises the question of the origin of data and of methodological best practices. Walters et al. (2018) reviewed the literature, highlighting common approaches and good practice and identifying research gaps. They extracted information from 78 records and found that most epidemiological data come from published journal articles, population data come from a wide range of sources and travel data mainly come from statistics or surveys. Rhodes et al. (2020) traced how models can be investigated as matters of correspondence and enactment in relation to their social and policy contexts.
At the beginning of a disease outbreak, there is a small number of infectious individuals and the transmission of infection is a stochastic event depending on the pattern of contacts between members of the population. The Watson and Galton (1874) process was one of the first approaches to successfully described this pattern. An alternative consists of introducing random noise into differential equations defining each compartment. Zhang and Wang (2013) explored this alternative and studied the asymptotic behavior of SIR model with Brownian noise and a jump process. The stochastic model containing a standard Brownian motion was studied by Caraballo and Colucci (2017). Caraball and Keraani (2018) explored the features of a stochastic SIR model with a fractional Brownian motion. The book by Daley and Gani (1999) contains an account of some of the more recent extensions.
In actuarial sciences, the literature on epidemics model is rather scarce. Jia and Tsui (2005) proposed and estimated a compartment model for severe acute respiratory syndrome (SARS) data. Chen and Cox (2009) employed the theory of real options and considered a regime-switching process for modeling the number of infected individuals. Feng and Garrido (2011) quantified the risk of infection with a classical epidemiological compartment model. They formulated financial arrangements between an insurer and insured using actuarial methodology and applied their framework to the SARS epidemic in 2003. Gathy and Lefèvre (2009) and Lefevre and Utev (1999) proposed extensions of deterministic compartment models and provided additional tools to account for randomness in epidemiological dynamics. Based on a Markov chain formulation of the susceptible-infected-recovered (SIR model), Lefèvre et al. (2017) developed a recursive method to calculate the cost components and the corresponding premium levels. More recently, Clara-Rahola (2020) proposed two distinct exponential models for the infection rate before and after lockdown. The fit to data from China, Spain, South Korea and Italy revealed that a crossover point between pre- and postlockdown infection rates is found one week after lockdown, which, in turn, is the average COVID-19 incubation period.

3. A Deterministic Epidemic Model

We consider a population of size S 0 hit 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:
I t = S 0 e α + μ t β t γ , t t 0 ,
where α , β , μ , γ R + . This function is the product of two terms. The first one is an exponential decreasing function, e α + μ t , whereas the second one is an increasing function. The number of infectious 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 differentiate Equation (1):
d I t = α + μ I t d t + S 0 e α + μ t β γ γ t γ 1 d t = α + μ I t d t + I t γ t d t .
This differential 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 first wave of COVID-19 in Belgium, Germany, Spain and Italy.
This model slightly differs from the susceptible-infected-recovered (SIR) model developed by Kermack and McKendrick (1927), which is a standard in the literature. In the SIR model, the evolution of the epidemics is described by the ordinary differential equations (ODEs):
d I t S I R = α S I R + μ S I R I t S I R d t + β S I R S t S I R S 0 I t S I R d t , d S t S I R = β S I R S t S I R S 0 I t S I R d t ,
where β S I R R + and S t S I R is the number of persons that are susceptible to be infected at time t. As in our model, α S I R and μ S I R are the recovery and mortality rates, respectively. The contagion rate per capita in the SIR is proportional to the population of susceptible, β S I R S t S I R S 0 , 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 benefit 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.
The basic reproduction number, R 0 , is defined 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 t S I R 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 us to take into account the impact of preventive measures to curb the epidemic, such as a lockdown or the wearing of masks. The empirical analysis of the next section confirms that this approach offers a better fit than the SIR model. Our model presents other interesting features. In particular, the peak of the epidemic is known and obtained by canceling the first order derivative of Equation (1):
t m a x = γ α + μ .
Combining Equations (1) and (4) allows us to evaluate the size of the infected population when the peak is reached:
I t m a x = S 0 e γ γ α + μ γ β γ .
Since the population of infectious persons may not exceed S 0 , we infer the necessary conditions β γ e γ t m a x γ . 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 differential equation:
d D t = μ I t d t .
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:
d S t = α I t d t I t γ t d t , t > t 0 .
By construction, the population of susceptible growths when infected individuals 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 individuals, infected individuals and deaths remains constant and equal to S 0 :
S t + I t + D t = S 0 , t 0 .

4. 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 susceptibles, 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 benefit rate is noted b. The plan terminates when the individual recovers or dies from the disease. In case of death, a lump sum benefit, c, is paid.
The risk-free rate is constant and denoted by r. If the insurance plan covers the whole population, premium, benefit rates and lump sum death capital must ensure the financial equilibrium of the plan. Under the assumption that the plan starts at time 0 and finishes at time T, discounted premiums have to cover all discounted benefits:
p 0 T e r s S s d s = b 0 T e r s I s d s + c 0 T e r s d D s .
As stated in the previous section, the discounted integral of I t admits a closed-form expression.
Proposition 1.
For r 0 , we have that:
0 T e r s I s d s = S 0 β γ θ γ + 1 Γ l γ + 1 , T θ ,
where θ = r + α + μ and Γ l ( γ + 1 , x ) = 0 x u γ e u d u is the gamma lower incomplete function.
Proof. 
Using Equation (1), we develop the integral as follows:
0 T e r s I s d s = S 0 β γ 0 T e r + α + μ s s γ d s ,
and we perform a change of variable u = θ s in order to rewrite the integral:
0 T e r s I s d s = S 0 β γ θ γ 1 0 θ T e u u γ d u .
From the definition of the gamma incomplete function, we directly obtain Equation (10). □
From this last proposition, we immediately infer the expressions of D t , 0 T e r s d D s and S t .
Corollary 1.
The cumulated number of deceases caused by the epidemic at time t 0 is equal to
D t = μ S 0 β γ α + μ γ 1 Γ l γ + 1 , t α + μ .
If θ = r + α + μ , the second term on the right-hand side of Equation (9) is:
0 T e r s d D s = μ S 0 β γ θ γ 1 Γ l γ + 1 , T θ .
The size of the population of susceptible at time t 0 is deduced from the relation S t + I t + D t = S 0 :
S t = S 0 S 0 e α + μ t β t γ μ S 0 β γ α + μ γ 1 Γ l γ + 1 , t α + μ .
The next proposition reports the closed-form expression for the premium rate solution of Equation (9).
Proposition 2.
For the benefit rate ( b , c ) , the fair premium rate that ensures the actuarial equilibrium of the plan is given by
p = b + c μ S 0 β γ θ γ 1 Γ l γ + 1 , T θ 0 T e r s S s d s ,
where the denominator is equal to:
0 T e r s S s d s = S 0 r 1 e r T S 0 β γ θ γ + 1 Γ l γ + 1 , T θ 1 + μ r + μ S 0 β γ r α + μ γ + 1 e r T Γ l γ + 1 , T α + μ .
Proof. 
From Equation (13), we infer that
0 T e r s S s d s = S 0 0 T e r s d s 0 T e r s I s d s 0 T e r s D s d s
The first integral on the right-hand side is equal to S 0 r 1 e r T whereas 0 T e r s I s d s is provided by Equation (10). From Equation (11), we infer that the integral of the discounted number of deaths is:
0 T e r s D s d s = μ S 0 β γ α + μ γ + 1 0 T e r s Γ l γ + 1 , s α + μ d s .
Integrating by parts leads to the following result:
0 T e r s α + μ γ + 1 Γ l γ + 1 , s α + μ d s = e r s r α + μ γ + 1 Γ l γ + 1 , s α + μ s = 0 s = T 0 T e r s r d d s 0 s e α + μ v v γ d v d s .
The integral in this last expression may also be reformulated in terms of an incomplete lower gamma function:
0 T e r s r d d s 0 s e α + μ v v γ d v d s = 1 r 0 T e r + α + μ s s γ d s = 1 r θ γ + 1 Γ l γ + 1 , T θ .
Combining these results leads to Equation (15) and the fair premium comes from the actuarial equilibrium equation. □

5. Empirical Illustration

We fit the model to data about the COVID-19 outbreak in Belgium, Germany, Italy and Spain. The first three countries are selected because they have reported the highest death rates in Europe during the 2020 first wave of COVID-19. In comparison, Germany has better managed the spread of the virus but the distribution of infected individuals over time has decreased at a lower pace than other countries considered in this study. We use the datasets from the library “coronavirus1” in R which provides daily time-series of the number of deaths and detected cases of COVID from the beginning of 2020 up to the end of July. We choose as starting date the day when the number of confirmed 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 the considered countries. Figure 1 shows these time series and Table 1 reports some statistics of the datasets. For Spain, the number of confirmed 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, I t . As the data sets only report new confirmed cases, we assume that contaminated individuals remain infected, on average, for 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 t o b s t = 1 , , n o b s is the time series of observations and n o b s is the number of days, parameters are obtained by a weighted least-square minimization:
( α + μ ^ , γ ^ , β ^ ) = arg min k = 1 n o b s ω k I k I k o b s 2 k = 1 n o b s ω k ,
where I k is the value of I t at time t k . As the impacts of α and μ on I t are indistinguishable, we first 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:
ω k = 0.05 k { 1 , , 25 } , ω k = 0.5 k { 36 , , 73 } , ω k = 1.5 k 74 .
These weights are chosen in order to obtain the best compromise fitting both the tail and the peaks of the I t curve. 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, fitted by least-square minimization. Our numerical experiments reveals that the SIR model fails to replicate the curve of I t . The only way to fit 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 offers at a first sight an excellent fit 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 confirms that our approach is a reliable alternative compared to 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 4. Two cases are considered. In the first one, collected premiums cover exclusively medical expenses. An allowance of 1000 EUR per day is paid during the treatment ( c = 365,000 EUR on a yearly basis). The second plan covers exclusively the death risk: a lump sum capital of 200,000 EUR is paid at the decease of an infected patient. The duration of both plans is six months and the risk-free rate is set to 2%. Table 4 and Table 5 report the fair premium rates calculated with our approach and the SIR model, respectively, for Belgium, Germany, Italy and Spain. We also test the sensitivity of these rates to variations of parameters. Per country, premium rates computed with our approach or the SIR model are similar. The premium rates for both benefits (62.38 EUR and 42.35 EUR per year) for Germany are the lowest due to the low number of confirmed cases and deaths reported by this country. The death coverage is the most expensive for Belgium (338 EUR/year), whereas Italy and Spain are in the same range (230.52 EUR and 232.57 EUR). The healthcare benefit is most expensive in Spain and Belgium (143.1 EUR/year and 138.54 EUR/year).

6. Time-Changed Extension

The model introduced in Section 3 is fully deterministic. In practice, we observe random fluctuations of the number of infected persons. In order to replicate such random variations, we propose two stochastic extensions of Equation (1). The first one developed in this section consists of replacing the time t by a stochastic clock, also called a subordinator. This clock is an increasing positive process denoted by τ t t 0 , defined on a probability space Ω endowed with a probability measure P and its natural filtration 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
f τ t ( x ) = 1 x > 0 x λ t 1 e x Γ ( λ t ) ,
where Γ ( a ) = 0 x a 1 e x d x is the standard gamma function such that Γ ( a + 1 ) = a Γ ( a ) for a > 0 . A straightforward calculation shows that the characteristic function for the gamma process is given by
E e i u τ t = 1 i u λ t = e t ψ ( u ) ,
where ψ ( u ) = λ ln 1 i u is the characteristic exponent and u R . 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:
ψ ( u ) = 0 λ x 1 e x ( 1 e i u x ) d x ,
from which we infer that the Lévy measure of τ t is ν ( d x ) = λ x 1 e x 1 ( x 0 ) d x . τ t is a process with finite variations. Therefore, for any function f ( t , τ t ) of time and of the subordinator, Itô’s lemma for semimartingales states that
d f ( t , τ t ) = f ( t , τ t ) t d t + f ( t , τ t + Δ τ t ) f ( τ t ) ) .
Whereas the F t -expectation of its infinitesimal variation is
E d f ( t , τ t ) | F t = f ( t , τ t ) t d t + 0 f ( t , τ t + x ) f ( τ t ) λ x 1 e x d x d t .
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:
I t = S 0 e α + μ τ t β τ t γ t 0 ,
where α , β , μ , γ R + . Using Itô’s lemma and first-order Taylor developments, we infer that:
d I t = S 0 e α + μ τ t + Δ τ t β τ t + Δ τ t γ S 0 e α + μ τ t β τ t γ = I t α + μ Δ τ t + γ τ t Δ τ t + O Δ τ t 2 .
This first-order approximation emphasizes the strong relation between the deterministic and the time-changed dynamics. Parameters α 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 gives the first two moments of I t t 0 .
Proposition 3.
The expected number of infected persons at time t is equal to:
E I t | F 0 = S 0 β γ Γ γ + λ t α + μ + 1 γ + λ t Γ ( λ t ) ,
whereas its variance is given by the following relation:
V I t | F 0 = S 0 2 β 2 γ Γ 2 γ + λ t 2 α + 2 μ + 1 2 γ + λ t Γ ( λ t ) Γ γ + λ t α + μ + 1 γ + λ t Γ ( λ t ) 2 .
Proof. 
The expectation of I t is rewritten in its integral form:
E I t | F 0 = S 0 0 u λ t 1 e u Γ ( λ t ) e α + μ u β u γ d u = S 0 β γ Γ ( λ t ) 0 u γ + λ t 1 e α + μ + 1 u d u .
Next, we do a change of variable v = α + μ + 1 u in order to obtain Equation (21). We obtain the moments of second order in a similar manner:
E I t 2 | F 0 = S 0 2 0 u λ t 1 e u Γ ( λ t ) e 2 α + μ u β u 2 γ d u = S 0 2 β 2 γ Γ ( λ t ) 0 u 2 γ + λ t 1 e 2 α + 2 μ + 1 u d u = S 0 2 β 2 γ Γ 2 γ + λ t 2 α + 2 μ + 1 2 γ + λ t Γ ( λ t ) .
Equation (22) is the difference of this second moment and of the square of the expectation. □
Notice that the maximum of the epidemics is reached at t m a x with τ t m a x = γ α + u . The cumulated number of deaths is the time-changed version of Equation (11):
D t = μ S 0 β γ α + μ γ 1 Γ l γ + 1 , τ t α + μ = μ S 0 β γ α + μ γ 1 0 τ t α + μ u γ e u d u = μ S 0 β γ 0 τ t v γ e α + μ v d v .
We use Itô’s lemma and first-order Taylor developments to check that the dynamics of D t is compliant with the one of I t :
d D t = μ S 0 β γ τ t τ t + Δ τ t v γ e α + μ v d v μ S 0 β τ t γ e α + μ τ t I t Δ τ t + O Δ τ t 2 .
which is the stochastic equivalent equation to (6). This equation confirms that the infinitesimal 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 semiclosed form expression and their valuation requires numerical integration:
E D t | F 0 = μ S 0 β γ 0 x λ t 1 e x Γ ( λ t ) 0 x v γ e α + μ v d v d x = μ S 0 β γ Γ ( λ t ) α + μ γ + 1 0 x λ t 1 e x 0 α + μ x u γ e u d u d x = μ S 0 β γ Γ ( λ t ) α + μ γ + 1 0 x λ t 1 e x Γ l γ + 1 , x α + μ d x .
and
E D t 2 | F 0 = μ S 0 β γ 2 0 x λ t 1 e x Γ ( λ t ) 0 x v γ e α + μ v d v 2 d x = μ S 0 β γ 2 Γ ( λ t ) α + μ 2 γ + 2 0 x λ t 1 e x 0 α + μ x u γ e u d u 2 d x = μ S 0 β γ 2 Γ ( λ t ) α + μ 2 γ + 2 0 x λ t 1 e x Γ l γ + 1 , x α + μ 2 d x .
The variance of the cumulated number of deaths up to time t is therefore:
V D t | F 0 = μ S 0 β γ 2 Γ ( λ t ) α + μ 2 γ + 2 0 x λ t 1 e x Γ l γ + 1 , x α + μ 2 d x 0 x λ t 1 e x Γ l γ + 1 , x α + μ d x 2 .
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:
S t = S 0 S 0 e α + μ τ t β τ t γ μ S 0 β γ α + μ γ 1 Γ l γ + 1 , τ t α + μ .
The expected size of the population of susceptible is simply equal to E ( S t | F 0 ) = S 0 E D t | F 0 E I t | F 0 where E I t | F 0 and E D t | F 0 are respectively provided by Equations (21) and (25).
If we consider the insurance plan introduced in Section 4, the fair premium rate that finances expected benefits is such that
p = b 0 T e r s E I s | F 0 d s + c 0 T e r s E d D s | F 0 0 T e r s E S s | F 0 d s .
Contrary to the deterministic model, the integrals present in this last expression do not admit a 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 equispaced times and if we note by Δ m the length of interarrival times, the integrals are computed as:
0 T e r s E S s | F 0 d s i = 1 m e r s i E S s i | F 0 Δ m , 0 T e r s E d D s | F 0 i = 1 m e r s i E D s i | F 0 E D s i 1 | F 0 , 0 T e r s E I s | F 0 d s i = 1 m e r s i E I s i | F 0 Δ m .
In the next Section, we present a different stochastic extension of our deterministic model that leads to analytical expressions for these integrals.

7. 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 financial option with the underlying being 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 measures2.
Proposition 4.
The value of an excess-of-loss reinsurance covering an excessive number of infectious is equal to:
C e r t E I t K + | F 0 = C e r t S 0 β γ Γ u λ t + γ , α + μ + 1 u K Γ ( λ t ) α + μ + 1 λ t + γ K Γ u λ t , u K Γ ( λ t ) ,
where u K is the positive solution of the following equation
ln S 0 ( α + μ ) u + γ ln β + γ ln u = ln K ,
and Γ u ( γ + 1 , x ) = x u γ e u d u is the upper gamma incomplete function.
Proof. 
We insert the definition of I t and rewrite the expectation as an integral with respect to the density of τ t :
C e r t E I t K + | F 0 = C e r t 0 u λ t 1 e u Γ ( λ t ) S 0 e α + μ u β u γ K + d u .
The integrand is positive if and only if u is above u K , such as it is defined in Equation (29). This allows us to rewrite the integral as a difference
0 u λ t 1 e u Γ ( λ t ) S 0 e α + μ u β u γ K + d u = S 0 β γ Γ ( λ t ) u K u λ t + γ 1 e α + μ + 1 u d u K u K u λ t 1 e u Γ ( λ t ) d u ,
where the second term is the ratio K Γ u λ t , u K Γ ( λ t ) . We replace the integrating variable by v that is α + μ + 1 u = v and infer that the first term is:
u K u λ t + γ 1 e α + μ + 1 u d u = α + μ + 1 u K v λ t + γ 1 e v d v α + μ + 1 λ t + γ = Γ u λ t + γ , α + μ + 1 u K α + μ + 1 λ t + γ ,
and we recover Equation (29). □
Notice that, in practice, the reinsurer will price reinsurance treaties with a set of parameters more conservative than the ones fitted 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:
C e r t E D t K + | F 0 = C e r t μ S 0 β γ α + μ 1 + γ u K u λ t 1 e u Γ ( λ t ) Γ l γ + 1 , u α + μ d u K K Γ u λ t , u K Γ ( λ t ) .
Unfortunately, the integral in this last equation does not admit a simple analytical form and must be computed numerically.

8. Estimation and Illustration

In order to illustrate the ability of the time-changed model to explain the evolution of a pandemic, we fit it to COVID-19 data sets for Belgium, Germany, Italy and Spain. As in Section 5, 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 to the number of infected persons adjusted by a time coefficient. More precisely, given that τ t has increments that are gamma-distributed, E Δ τ t = λ d t . From Equation (24), the expectation is E d D t | F 0 μ E I t | F 0 λ d t . Using a moment-matching approach, we then estimate the mortality rate as follows:
μ ^ = k = 1 n o b s d D k o b s λ ^ k = 1 n o b s I k d t .
Parameter estimates are presented in Table 6 and Figure 3 compares the expected number of infectious obtained with the time-changed model and its deterministic counterpart. Globally, we do not observe any similarities between parameters of time-changed and deterministic models. The goodness of fit, measured by the SSE, is also worse for the time-changed version than for the deterministic one. The time change does not fit the right tail of the I t curve better and overestimates, on average, the number of infected individuals at the beginning of the pandemic.
Table 7 shows the fair premium rates for the healthcare and death insurances valued in Section 5. Since the time-changed model predicts on average higher healthcare benefits during the growing phase of the outbreak, the healthcare insurance is slightly more expensive than the one valued with the deterministic model. The death benefits and death insurance premium rates are comparable in both models.
If we limit our analysis to a comparison of the expected number of infected individuals 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 different behavior of this model, we simulate 2000 samples paths with parameter estimates obtained for Italy. Figure 4 shows three of these paths and the average over the 2000 simulations. We see that the time-changed model generates curves of I t with a different 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 it 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 figures the same conclusions. The discontinuities of I t and the bimodal behavior of D t being unlikely in practice (at least for COVID-19), we investigate in the next section, an alternative stochastic extension of the deterministic model.

9. A Jump Diffusion 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 a Brownian motion W t t 0 and a compound jump process L t t 0 are defined. 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 defined on R + with a probability density function denoted by f J ( . ) . The expectation and variance of a jump are denoted by μ J = E ( J ) and V ( J ) = σ J 2 , respectively. The compound Poisson process L t t 0 is defined as the sum of jumps J k up to time t:
L t = k = 1 N t J k .
We assume that the dynamics of the population of infected individuals is ruled by the following geometric jump diffusion:
d I t = α + μ I t d t + I t γ t d t + σ I t d W t + I t d L t ,
where α , μ , γ t are the recovery, mortality and contagion rates, respectively. The term σ I t d W t , with σ R + , is a Gaussian noise, whereas I t d L 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:
I t = S 0 e α + μ + σ 2 2 t + σ W t β t γ k = 1 N t 1 + J k .
This result is a direct consequence of Itô’s lemma for jump diffusion. 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 μ d t of the population of sick persons. The differential equation ruling S t , the population of susceptible individuals, guarantees that the total size of the population remains equal to S 0 :
d D t = μ I t d t , d S t = I t α γ t d t σ I t d W t I t d L t .
The number of deaths is D t = μ 0 t I s d s , 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 first two moments of I t .
Proposition 5.
The expectation and variance of the number of infected persons at time t are respectively equal to
E I t | F 0 = S 0 e λ μ J α + μ t β t γ ,
and
V I t | F 0 = S 0 2 e 2 α + μ t β t 2 γ e 2 λ μ J t e σ 2 + λ E ( J 2 ) t 1 .
Proof. 
Given that W t is independent from N t , the expectation of I t is equal to a product of expectations:
E I t | F 0 = S 0 e α + μ t β t γ E e σ 2 2 t + σ W t | F 0 E k = 1 N t 1 + J k | F 0 .
As σ W t N ( 0 , σ 2 t ) is normal, E e σ 2 2 t + σ W t | F 0 = 1 . The moment-generating function of N t is E e ω N t = e λ t e ω 1 and N t is independent for J k k = 1 . . . , N t . Therefore, we conclude that the expectation of the product of jumps is:
E k = 1 N t 1 + J k | F 0 = E k = 1 N t 1 + μ J | F 0 = E e N t ln 1 + μ J | F 0 = e λ μ J t .
Combining these elements leads to Equation (35). In a similar manner, we calculate the second order moment of I t —that is,
E I t 2 | F 0 = S 0 2 e 2 α + μ t β t 2 γ E e σ 2 t + 2 σ W t | F 0 E k = 1 N t 1 + J k 2 | F 0 .
Since e 2 σ W t is log-normal, E e 2 σ W t | F 0 = e 2 σ 2 t and E e σ 2 t + 2 σ W t | F 0 = e σ 2 t . Furthermore, the expectation of the square of the product of jumps is equal to:
E k = 1 N t 1 + J k 2 | F 0 = E E k = 1 N t 1 + J k 2 | F 0 N t | F 0 = E k = 1 N t E 1 + J 2 | F 0 = E e N t ln ( 1 + 2 μ j + E ( J 2 ) ) | F 0 = e λ 2 μ J + E ( J 2 ) t .
We obtain, then, the second moment of I t and the variance in (36). □
The next proposition presents an analytical formula to evaluate expected healthcare benefits 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 6.
Let us consider a discount rate r R + . The integral of expected discounted numbers of infected is equal to:
0 t e r s E I s | F 0 d s = S 0 β γ θ J γ + 1 Γ l γ + 1 , θ J t ,
where θ J = r + α + μ λ μ J and Γ l ( γ + 1 , x ) is the gamma lower incomplete function.
The proof is similar to the one of Proposition 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 closed-form expressions in terms of incomplete lower gamma functions.
Corollary 2.
The expected cumulated number of deaths caused by the epidemic at time t 0 is equal to
E D t | F 0 = μ S 0 β γ α + μ λ μ J γ + 1 Γ l γ + 1 , t α + μ λ μ J .
If θ J = r + α + μ λ μ J , the expectation of the integral of discounted variation of D t is equal to:
E 0 T e r s d D s | F 0 = μ S 0 β γ θ J γ + 1 Γ l γ + 1 , θ J t .
Let us again consider the insurance plan introduced in Equation (4) with a maturity T. The premium rate is still denoted as p. The rate of healthcare expenses is b and the lump sum benefit 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 7.
For the benefit rates ( b , c ) , the fair premium rate that guarantees the actuarial equilibrium of the plan is given by:
p = b + c μ S 0 β γ θ J γ 1 Γ l γ + 1 , θ J T 0 T e r s E S s | F 0 d s ,
where the denominator is equal to:
0 T e r s E S s | F 0 d s = S 0 r 1 e r T S 0 β γ θ J γ + 1 Γ l γ + 1 , θ J T 1 + μ r + μ S 0 β γ e r T Γ l γ + 1 , α + μ λ μ J T r α + μ λ μ J γ + 1 .
Proof. 
The fair premium rate in a stochastic framework is given by Equation (28). The integral 0 T e r s E I s | F 0 d s is provided by Equation (37). As S t = S 0 I t D t , we infer that
0 T e r s E S s | F 0 d s = S 0 r 1 e r T S 0 β γ θ J γ + 1 Γ l γ + 1 , θ J T 0 T e r s E D s | F 0 d s ,
where 0 T e r s E D s | F 0 d s admits the following integral representation:
0 T e r s E D s | F 0 d s = μ 0 T e r s 0 s E I u | F 0 d u d s = μ S 0 β γ 0 T e r s 0 s e λ μ J α + μ u u γ d u d s .
Given that 0 s e λ μ J α + μ u u γ d u = Γ l γ + 1 , α + μ λ μ J s α + μ λ μ J γ + 1 , the double integral in this last expression becomes:
0 T e r s 0 s e λ μ J α + μ u u γ d u d s = 1 r e r s Γ l γ + 1 , α + μ λ μ J s α + μ λ μ J γ + 1 s = 0 s = T + 1 r 0 T e λ μ J r + α + μ s s γ d s = e r T Γ l γ + 1 , α + μ λ μ J T r α + μ λ μ J γ + 1 + Γ l γ + 1 , θ J T r θ J γ + 1
Therefore, the integral of the discounted expected number of deaths is equal to:
0 T e r s E D s | F 0 d s = μ S 0 β γ e r T Γ l γ + 1 , α + μ λ μ J T r α + μ λ μ J γ + 1 + μ S 0 β γ Γ l γ + 1 , θ J T r θ J γ + 1 .
Combining expressions (42) and (43) leads to Equation (41). □

10. Reinsurance in the Jump Diffusion 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 payoff 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:
I t | N t = n = S 0 β t γ e X t ( n ) ,
where X t ( n ) N μ X ( t , n ) ; σ X 2 ( t , n ) . The mean and variance of X t ( n ) are respectively equal to
μ X ( t , n ) = n ln ( 1 + κ ) α + μ + σ 2 2 t , σ X 2 ( 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 to be the same under the pricing and real measures.
Proposition 8.
Let us define
d 2 ( t , n ) = 1 σ X ( t , n ) ln K S 0 β t γ μ X ( t , n ) ,
d 1 ( t , n ) = d 2 ( t , n ) σ X ( t , n ) .
The value of an excess-of-loss reinsurance covering an excessive number of infected is equal to:
C e r t E I t K + | F 0 = C e r t n = 0 P ( N t = n ) S 0 β t γ e μ X + σ X 2 2 Φ d 1 ( t , n ) K Φ d 2 t , n ,
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 :
C e r t E I t K + | F 0 = C e r t n = 0 P ( N t = n ) E S 0 β t γ e X t ( n ) K + | F 0 .
The conditional expectations may be developed as the difference of two integrals:
E S 0 β t γ e X t ( n ) K + | F 0 = S 0 β t γ e μ X + σ X u K S 0 β t γ e μ X + σ X u e u 2 2 2 π d u K e μ X + σ X u K S 0 β t γ e u 2 2 2 π d u .
The second term, after a change of variable, is equal to:
K e μ X + σ X u K S 0 β t γ e u 2 2 2 π d u = K u 1 σ X ln K S 0 β t γ μ X e u 2 2 2 π d u = K Φ d 2 t , n .
where d 2 ( t , n ) = 1 σ X ( t , n ) ln K S 0 β t γ μ X ( t , n ) . In the same manner, the first term of Equation (47) becomes:
e μ X + σ X u K S 0 β t γ e μ X + σ X u e u 2 2 2 π d u = e μ X + σ X 2 2 u 1 σ X ln K S 0 β t γ μ X e u 2 2 σ X u + σ X 2 2 2 π d u = e μ X + σ X 2 2 s 1 σ X ln K S 0 β t γ μ X σ X e s 2 2 2 π d s = e μ X + σ X 2 2 Φ d 1 ( t , n ) ,
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.

11. Estimation of the Jump Diffusion Model

The jump diffusion model is fitted 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 o b s . The first step consists in estimating β , γ and η by minimizing the weighted sum of squares between the expected and observed numbers of infected persons:
( η ^ , γ ^ , β ^ ) = arg min k = 1 n o b s ω k E I k | F 0 I k o b s 2 k = 1 n o b s ω k .
Given that the expectation of I t in the jump diffusion approach coincides with the deterministic model, we obtain the same estimates β ^ , γ ^ as those in Table 2. In the second stage, we fit the jump process by the peak over threshold method. From Equations (33) and (35), we define Y t as the ratio of the process I t on its expectation:
Y t : = I t E I t | F 0 = e λ μ J σ 2 2 t + σ W t k = 1 N t 1 + J k .
Using Itô’s lemma, we infer that Y t is driven by the following infinitesimal dynamics:
d Y t Y t = + σ d W t + d L t λ μ J d t .
The value of Y t at time t k is noted Y k and we define Z k = Y k Y k 1 Y k the discrete approximation of d Y t Y t . If the time lag of one day between two successive observations is noted Δ , according to Equation (49), Z k is approximately the sum
Z k ϵ k + Δ L k λ μ J Δ ,
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 confidence level. To define the threshold, we fit a pure Gaussian process to the time series of Z k μ z Δ + σ z W Δ . The unbiased estimators of μ z and σ z are:
μ ^ z = 1 n o b s Δ j = 1 n o b s Z j σ ^ z 2 = 1 n o b s 1 Δ j = 1 n o b s Z j μ ^ z 2 .
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 t h 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:
N ( t k ) = max { j N | s j t k } .
Under the assumption that the diffusion 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 fit 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 confirmed 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 confidence level is set to 90%. Figure 6 shows those series and the threshold level. For Spain, we have not taken two abnormally high and low values into account due to retrospective adjustments of the number of confirmed cases. For Belgium and Spain, the Z k values have a mean close to zero and the variance seems more or less constant for different time windows. This confirms 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 values have a residual linear increasing trend and their variance seems to increase as it does 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.
By construction, the expectation of the jump diffusion model corresponds to the deterministic model of Section 3. Therefore, the fair premium rate of an insurance plan covering healthcare expenses and death benefits are the same with both approaches. We refer the reader to the previous Section 5 for numerical examples. However, the jump diffusion 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 approach, the jump diffusion model generates smoother sample paths which all appear as likely scenarios for the evolution of the infected population. These graphs also display the expected and the observed number of infected cases for each country. It reveals that the real sample path of I t for each country is a likely realization of the jump diffusion model, at least for t 0.07 . At the beginning of the pandemic, the real sample path of I t bounds from below simulated trajectories. This is a direct consequence of choices made in Section 5 for calibrating the average trend of the model. At the start of the epidemic, testing polices were in the process of being deployed and the number of infected cases was probably underestimated. This motivates our choice of underweighting observations collected in the early ascending phase. We also remark that the peak of infected cases is significantly higher than the observed one in some scenarios. To illustrate this, Table 10 reports the 90% and 95% percentiles of the simulated maximum of infected cases. The ratios of these 95% percentiles on real numbers of infected cases at the epidemic peak range from 121.23% for Italy up to 202.49% for Spain.
Table 11 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 12 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 set to t = 0.1 year and C = 1 , respectively. 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 .

12. Conclusions

The valuation of actuarial commitments requires us to integrate over time the size of infectious and susceptible populations. Since existing compartment models do not admit closed-form expressions for these quantities, actuarial calculations are in this framework computationally intensive and subject to numerical errors. The three models proposed in this article remedy to this issue and present a high degree of analytical tractability.
The first model is purely deterministic. 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 confirm that the model explains the first wave of such an epidemic. Furthermore, the insurance premium rate admits a closed-form expression within this framework. The main disadvantage of this approach is the absence of random effects that prevents to evaluate the incurred extreme costs.
The second model is a time-changed extension of the deterministic one. The time of the pandemic peak is randomized by observing the process on a stochastic time scale. The main advantage of this approach is that it preserves the main features of the deterministic model and leads to comparable premium rates. Nevertheless, the simulation study reveals that simulated scenarios display a different trend from what is observed for the COVID-19 outbreak. Furthermore, the stochastic clock modulates the speed at which the pandemic evolves but do not modify the sizes of infected and susceptible populations.
This article proposes a second stochastic extension of the deterministic model based on a jump diffusion process. In this approach, the rate at which patients cease to be considered as infected is noised by a Brownian motion. This allows us to randomize the duration of illness. The apparition of local clusters of infected causing a sudden increase of the number of contagious cases is replicated by the jump component. This model presents several interesting features. As it behaves, on average, as the deterministic approach does, it keeps a high analytical tractability for actuarial applications. On the other hand, the model is able to generate realistic noised sample paths of infected cases. This feature allows us to price reinsurance contracts, such as excess of loss treaties, that cannot be valued in a deterministic framework. Last but not least, it is remarkably easy to estimate its parameters with the proposed “peak-over-threshold” method.
Notice that, by construction, the contagion rate per capita decreases as 1 t , the size of the infected population, converges to zero after having reached the epidemic peak. In a similar way to the SIR model, the solutions proposed in this article are then designed for explaining a single epidemic wave with, eventually, random recovery duration and discovery of local clusters of infected individuals.
This observations opens the way to further research. Instead of a deterministic starting date for the epidemic, we can replace it by the jump time of a self-exciting point process, e.g., Hainaut and Moraux (2019). In its simplest version, the intensity of this process is persistent and suddenly increases as soon as a jump occurs. Within this approach, the starting date of the pandemic becomes random and the probability of observing a new epidemic wave raises after the first one but decay exponentially to its baseline level. Other possible extensions consist to randomize the mortality rates or to develop a compartmental version with subpopulations of infected individuals.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The author is grateful to the “Fonds de la Recherche Scientifique—FNRS” for financial support under Grant number 33658713.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Anderson, Roy M., and Robert M. May. 2008. Population biology of infectious diseases. Nature 280: 361–67. [Google Scholar] [CrossRef] [PubMed]
  2. Brauer, Fred. 2017. Mathematical epidemiology: Past, present, and future. Infectious Disease Modelling 2: 113–27. [Google Scholar] [CrossRef] [PubMed]
  3. Caraballo, Tomas, and Renato Colucci. 2017. A comparison between random and stochastic modeling for a SIR model. Communications on Pure and Applied Analysis 16: 151–62. [Google Scholar] [CrossRef] [Green Version]
  4. Caraballo, Tomas, and Sami Keraani. 2018. Analysis of a stochastic SIR model with fractional Brownian motion. Stochastic Analysis and Applications 36: 895–908. [Google Scholar] [CrossRef]
  5. Chen, Hua, and Samuel H. Cox. 2009. An option-based operational risk management model for pandemics. North American Actuarial Journal 13: 54–76. [Google Scholar] [CrossRef]
  6. Clara-Rahola, Joaquim. 2020. An empirical model for the spread and reduction of the COVID19 pandemic. Estudios de Economia Aplicada 38. [Google Scholar] [CrossRef]
  7. Daley, David J., and J. Gani. 1999. Epidemic Models: An Introduction. Cambridge Studies in Mathematical Biology 15. Cambridge: Cambridge University Press. [Google Scholar]
  8. Diekmann, O., J. A. P. Heesterbeek, and J. A. J. Metz. 1995. The legacy of Kermack and McKendrick. In Epidemic Models: Their Structure and Relation to Data. Edited by D. Mollison. Cambridge: Cambridge University Press, pp. 95–115. [Google Scholar]
  9. Diekmann, Odo, J. A. P. Heesterbeek, and Johan A. J. Metz. 1990. On the definition and the computation of the basic reproductive ratio R0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology 28: 365–82. [Google Scholar] [CrossRef] [Green Version]
  10. Feng, Runhuan, and Jose Garrido. 2011. Actuarial applications of epidemiological models. North American Actuarial Journal 15: 112–27. [Google Scholar] [CrossRef]
  11. Gathy, Maude, and Claude Lefèvre. 2009. From damage models to SIR epidemics and cascading failures. Advances in Applied Probability 41: 247–69. [Google Scholar] [CrossRef] [Green Version]
  12. Hainaut, Donatien, and Franck Moraux. 2019. A switching self-exciting jump diffusion process for stock prices. Annals of Finance Volume 15: 267–306. [Google Scholar] [CrossRef]
  13. Hethcote, Herbert W. 2000. The Mathematics of Infectious Diseases. SIAM Review 42: 599–653. [Google Scholar] [CrossRef] [Green Version]
  14. Jia, Na, and Lawrence Tsui. 2005. Epidemic Modelling Using SARS as a Case Study. North American Actuarial Journal 9: 28–42. [Google Scholar] [CrossRef]
  15. Kermack, William Ogilvy, and A. G. McKendrick. 1927. Contributions to the mathematical theory of epidemics—Part I. Proceedings of the Royal Society of London. Series A 115: 700–21. [Google Scholar]
  16. Lefèvre, Claude, and Sergey Utev. 1999. Branching Approximation for the Collective Epidemic Model. Methodology and Computing in Applied Probability Volume 1: 211–28. [Google Scholar] [CrossRef]
  17. Lefèvre, Claude, Picard Philippe, and Matthieu Simon. 2017. Epidemic risk and insurance coverage. Journal of Applied Probabilit 54: 286–303. [Google Scholar] [CrossRef] [Green Version]
  18. Rhodes, Tim, Kari Lancaster, Shelley Lees, and Melissa Parker. 2020. Modelling the pandemic: Attuning models to their contexts. BMJ Global Health 5: e002914. [Google Scholar] [CrossRef]
  19. Smith, Dominic. 2017. Pandemic Risk Modelling. In The Palgrave Handbook of Unconventional Risk Transfer. Edited by M. Pompella and N. A. Scordis. Cham: Palgrave Macmillan, pp. 463–495. [Google Scholar]
  20. Tchuenche, Jean M., Alexander Nwagwo, and Richard Levins. 2007. Global behaviour of an SIR epidemic model with time delay. Mathematical Methods in the Applied Sciences 30: 733–49. [Google Scholar] [CrossRef]
  21. Van den Driessche, P., and James Watmough. 2020. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences 180: 29–48. [Google Scholar] [CrossRef]
  22. Walters, Caroline E., Margaux M.I. Meslé, and Ian M. Hall. 2018. Modelling the global spread of diseases: A review of current practice and capability. Epidemics 25: 1–8. [Google Scholar] [CrossRef]
  23. Watson, H. W., and Francis Galton. 1874. On the probability of the extinction of families. The Journal of the Anthropological Institute of Great Britain and Ireland 4: 138–44. [Google Scholar] [CrossRef]
  24. Zhang, Xianghua, and Ke Wang. 2013. Stochastic SIR model with jumps. Applied Mathematics Letters 26: 867–74. [Google Scholar] [CrossRef]
1
See, https://github.com/RamiKrispin/coronavirus, package developed by Rami Krispin.
2
In practice, the reinsurer would use a set of parameters more conservative than the estimated one in order to include a safety margin.
Figure 1. Graphs of the number of COVID-19 reported cases and deaths from March to June, for Belgium, Germany, Italy and Spain.
Figure 1. Graphs of the number of COVID-19 reported cases and deaths from March to June, for Belgium, Germany, Italy and Spain.
Risks 09 00003 g001
Figure 2. Comparison of fitted I t with observed ones and those computed with the SIR model.
Figure 2. Comparison of fitted I t with observed ones and those computed with the SIR model.
Risks 09 00003 g002
Figure 3. Comparison of E I t | F 0 with the observed values and those computed with the deterministic model of Section 6.
Figure 3. Comparison of E I t | F 0 with the observed values and those computed with the deterministic model of Section 6.
Risks 09 00003 g003
Figure 4. Simulation of three sample paths for I t and comparison with the simulated average of 1000 sample paths for Italy.
Figure 4. Simulation of three sample paths for I t and comparison with the simulated average of 1000 sample paths for Italy.
Risks 09 00003 g004
Figure 5. Histograms of the cumulated number of deaths at time t = 0.1 and t = 0.25 for the Italian population. 2000 simulations.
Figure 5. Histograms of the cumulated number of deaths at time t = 0.1 and t = 0.25 for the Italian population. 2000 simulations.
Risks 09 00003 g005
Figure 6. Plot of Z k time series for Belgium, Germany, Italy and Spain. The threshold is calculated with q = 90 % .
Figure 6. Plot of Z k time series for Belgium, Germany, Italy and Spain. The threshold is calculated with q = 90 % .
Risks 09 00003 g006
Figure 7. Simulations of 1000 sample paths of I t , with parameter estimates in Table 9 and constant J. The thick dotted blue line is the expectation of I t , whereas the thick black line is the observed number of infected cases.
Figure 7. Simulations of 1000 sample paths of I t , with parameter estimates in Table 9 and constant J. The thick dotted blue line is the expectation of I t , whereas the thick black line is the observed number of infected cases.
Risks 09 00003 g007
Table 1. Statistics about the time series of deaths and confirmed COVID cases.
Table 1. Statistics about the time series of deaths and confirmed COVID cases.
CountryStarting DateNumber of DaysNumber of Infected IndividualsNumber of DeathsPopulation Size, S 0
Belgium7/3/202010159,991966111,589,623
Germany8/3/2020100186,883880783,770,952
Italy27/2/2020110236,83734,35960,461,826
Spain7/3/2020101243,70927,13146,934,632
Table 2. Parameter estimates for the model ruled by Equation (1) (unit of time t: year). The SSE is the value of the optimization criterion in Equation (16).
Table 2. Parameter estimates for the model ruled by Equation (1) (unit of time t: year). The SSE is the value of the optimization criterion in Equation (16).
α ^ γ ^ β ^ μ ^ t max (days) I max SSE
Belgium40.7184.746.6064.4573817,829535,988
Germany40.6333.1243.6931.2392765,66512,513,688
Italy30.8783.3823.7093.9313565,1033,093,842
Spain46.6313.9376.9792.9662989,47924,203,834
Table 3. Parameter estimates of the SIR model.
Table 3. Parameter estimates of the SIR model.
Adjusted S ^ 0 α ^ SIR β ^ SIR μ ^ SSE
Belgium38,20312.407107.764.45716,109
Germany128,25016.816131.6491.2391,372,666
Italy118,25010.616113.4253.9311,175,335
Spain171,75015.914146.5162.9662,743,222
Table 4. Fair premium rates with the model in (1). Duration T = 0.5 year and r = 2 % .
Table 4. Fair premium rates with the model in (1). Duration T = 0.5 year and r = 2 % .
bcFair p β + 1 % β 1 % α + 1 % α 1 % γ + 1 % γ 1 %
Be365,0000138.54145.23132.09131.58145.93136.83140.32
0200,000338.35354.7322.59321.36356.41334.19342.71
Ge365,000062.3864.3560.4559.9564.9460.2164.65
0200,00042.3543.6941.0440.7044.0940.8843.89
It365,0000107.02110.68103.44102.95111.28103.89110.27
0200,000230.52238.42222.81221.77239.7223.78237.53
Sp365,0000143.1148.82137.54136.64149.94140.5145.8
0200,000232.57241.87223.54222.07243.68228.35236.95
Table 5. Fair premium rates with the model (1). Duration T = 0.5 year and r = 2 % .
Table 5. Fair premium rates with the model (1). Duration T = 0.5 year and r = 2 % .
bcFair p, SIR β SIR + 1 % β SIR 1 % α SIR + 1 % α SIR 1 %
Be365,0000142.43142.45142.4141.38143.49
0200,000347.84347.89347.79345.29350.43
Ge365,000061.8661.8661.8561.2862.44
0200,00042.0042.0041.9941.6142.39
It365,000097.9597.9597.9497.2498.66
0200,000210.98211.00210.97209.47212.52
Sp365,0000141.51141.52141.51140.33142.72
0200,000229.99230229.98228.07231.95
Table 6. Parameter estimates for the model ruled by Equation (19) (unit of time t: year). The SSE is the value of the optimization criterion in Equation (16).
Table 6. Parameter estimates for the model ruled by Equation (19) (unit of time t: year). The SSE is the value of the optimization criterion in Equation (16).
α ^ γ ^ β ^ μ ^ λ ^ SSE
Belgium44.806105.0211.1130.14927.7761,005,094
Germany17.20829.8661.2950.04129.35318,592,728
Italy12.40725.0881.090.14925.3075,906,960
Spain67.31181.6022.1330.11223.42740,570,084
Table 7. Fair premium rates with the model in (1). Duration T = 0.5 year and r = 2 % .
Table 7. Fair premium rates with the model in (1). Duration T = 0.5 year and r = 2 % .
bcFair pModel 1, Fair p λ + 1 % λ 1 %
Be365,0000148.76138.54147.29150.26
0200,000336.76338.35336.77336.76
Ge365,000063.9862.3863.3464.62
0200,00041.942.3541.941.9
It365 0000110.99107.02109.9112.11
0200 000228.32230.52228.32228.31
Sp365,0000161.6143.1160163.23
0200,000227.73232.57227.73227.73
Table 8. Expectations, standard deviations, 5% and 95% percentiles of simulated D t . 2000 simulations.
Table 8. Expectations, standard deviations, 5% and 95% percentiles of simulated D t . 2000 simulations.
E ( D 0.25 ) V ( D 0.25 ) E ( D 0.10 ) V ( D 0.10 ) D 0.10 : 5% D 0.10 : 95%
PercentilePercentile
Belgium9574.0271005.1825138.0694504.64809718.298
Germany8708.296432.4376288.873403.71208741.932
Italy33,735.4733497.94419071.4114,798.839034,382.855
Spain26,533.641116.69519,959.17110,804.449026,602.104
Table 9. Parameter estimates for the model ruled by Equation (33) (unit of time t: year). The SSE is the value of the optimization criterion in Equation (16).
Table 9. Parameter estimates for the model ruled by Equation (33) (unit of time t: year). The SSE is the value of the optimization criterion in Equation (16).
α ^ γ ^ β ^ μ ^ σ ^ λ ^ μ J ^ σ J ^ SSE
Belgium41.9444.7406.6064.4570.58920.2780.060.007535,988
Germany43.7583.1253.6931.2390.72730.8450.1010.02212,513,688
Italy31.883.3823.7093.9310.37513.5190.0740.0243,093,842
Spain51.4653.9376.9792.9661.51225.3470.1910.02624,203,834
Table 10. Statistics about the observed and simulated epidemic peaks (1000 simulations).
Table 10. Statistics about the observed and simulated epidemic peaks (1000 simulations).
Maximum90% Percentile95% PercentileRatio 95%
Number ofMax. NumberMax. NumberPercentile on
Infectedof Infectedof InfectedMax. Cases
Belgium18,22523,70225,260138.6%
Germany71,21993,506100,764141.49%
Italy70,23379,82785,141121.23%
Spain97,400164707197,224202.49%
Table 11. Expectation, standard deviation and 5%–95% percentiles of simulated cumulated number of deaths, D t . 1000 simulations.
Table 11. Expectation, standard deviation and 5%–95% percentiles of simulated cumulated number of deaths, D t . 1000 simulations.
D 0.25 ExpectedStandard5%95%
DeviationPercentilePercentile
Belgium9379.4041850.3856664.90812,717.609
Germany8576.7882079.5895671.55412,163.033
Italy33,166.374704.3626,483.3141,795.193
Spain24,948.89512,595.32110,641.10746,904.273
Table 12. Prices of excess of loss reinsurance treaties. r = 2 % .
Table 12. Prices of excess of loss reinsurance treaties. r = 2 % .
K e r 0.10 E I 0.10 K + K e r 0.25 E D 0.25 K +
ine Belgium14,9983154.5369661610.114
Germany46,72111,411.3648807705.235
Italy59,2097088.72634,3591362.358
Spain60,92725,778.16427,1313739.249
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hainaut, D. An Actuarial Approach for Modeling Pandemic Risk. Risks 2021, 9, 3. https://0-doi-org.brum.beds.ac.uk/10.3390/risks9010003

AMA Style

Hainaut D. An Actuarial Approach for Modeling Pandemic Risk. Risks. 2021; 9(1):3. https://0-doi-org.brum.beds.ac.uk/10.3390/risks9010003

Chicago/Turabian Style

Hainaut, Donatien. 2021. "An Actuarial Approach for Modeling Pandemic Risk" Risks 9, no. 1: 3. https://0-doi-org.brum.beds.ac.uk/10.3390/risks9010003

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop