A Two-Phase Dynamic Contagion Model for COVID-19

Zezhun Chen
London School of Economics
Angelos Dassios
London School of Economics
Valerie Kuan
University College London
Jia Wei Lim
Brunel University London
Yan Qu*
University of Warwick
Budhi Surya
Victoria University of Wellington
Hongbiao Zhao
Shanghai University of Finance and Economics

10th June 2020

*Correspondence Email: y.qu3@lse.ac.uk

Abstract

In this paper, we propose a continuous-time stochastic intensity model, namely, two-phase dynamic contagion process (2P-DCP), for modelling the epidemic contagion of COVID-19 and investigating the lockdown effect based on the dynamic contagion model introduced by Dassios and Zhao (2011). It allows randomness to the infectivity of individuals rather than a constant reproduction number as assumed by standard models. Key epidemiological quantities, such as the distribution of final epidemic size and expected epidemic duration, are derived and estimated based on real data for various regions and countries. The associated time lag of the effect of intervention in each country or region is estimated. Our results are consistent with the incubation time of COVID-19 found by recent medical study. We demonstrate that our model could potentially be a valuable tool in the modeling of COVID-19. More importantly, the proposed model of 2P-DCP could also be used as an important tool in epidemiological modelling as this type of contagion models with very simple structures is adequate to describe the evolution of regional epidemic and worldwide pandemic.

Keywords: Stochastic intensity model; Stochastic epidemic model; Two-phase dynamic contagion process; COVID-19; Lockdown
Mathematics Subject Classification (2010): Primary: 60G55; Secondary: 60J75
JEL Classification: I1; H0; E6

1 Introduction

In the early stages of epidemic modelling, the spread of diseases was formulated as a deterministic process. The classical deterministic model of susceptible-infectious-recovered (SIR) was introduced in the seminal paper of Kermack and McKendrick (1927). It models the spread and ultimate containment of an infection in a setting where those who recover are immune to the disease and thus the susceptible population declines over time. Many epidemic models are variations of the SIR model, see Brauer et al. (2008); Keeling and Rohani (2008); Diekmann et al. (2013) and Martcheva (2015). For example, during the outbreak of COVID-19 since December 2019, a commonly adopted approach for predicting the number of infections is the susceptible-exposed-infectious-recovered (SEIR) model, which adds an exposed period to the SIR model for accounting the reported incubation period of COVID-19 during which individuals are not yet infectious, e.g. Berger et al. (2020); Liu et al. (2020) and Tian et al. (2020). More recently, Acemoglu et al. (2020) develop a multi-risk SIR model, which takes into account that different subpopulations have different risks and is applied to analysing optimal lockdown.

However, the random nature of the epidemics spread in our real world suggests that a stochastic model is needed. A continuous-time stochastic counterpart of SIR model was first proposed by McKendrick (1925), and then a variety of stochastic models were studied in the literature, e.g. Bartlett (19491956) and Bailey (195019531957). For more recent developments on stochastic epidemic models in general, see Daley and Gani (1999); Andersson and Britton (2000); Allen (2008) and Fuchs (2013). In particular, many researchers adopted branching processes. Ball (1983) used the birth-and-death process for constructing a sequence of general stochastic epidemics, and Ball and Donnelly (1995) used branching processes to approximate the early stages of epidemic dynamics, see also Britton (2010) and Ball et al. (2016).

In this paper, we propose a continuous-time stochastic epidemic model, namely, the two-phase dynamic contagion process (2P-DCP), for modelling the epidemic contagion. It is a branching process, and can be considered as a generalisation of dynamic contagion process (Dassios and Zhao2011), which is an extension of the classical Hawkes process (Hawkes1971a,b). In fact, Hawkes process and its various generalisations were originally used for modelling earthquakes in seismology, and recently become extremely popular for modelling financial contagion in economics, see Bowsher (2007); Large (2007); Embrechts et al. (2011); Bacry et al. (2013a,b); Aït-Sahalia et al. (2015); Dassios and Zhao (2017a,b) and Qu et al. (2019). Analogously, we advocate that they are also applicable to epidemiology. As not all individuals are equally infectious in reality, the main advantage of this Hawkes-based approach is that, it allows randomness to the infectivity of individuals, rather than a constant reproduction number (the average number of subsequent infections of an infected individual) in standard models. In this paper, we adopt the 2P-DCP as a more realistic and parsimonious example of Hawkes-based models for modelling the current progression of COVID-19 and investigating the lockdown effect. Key epidemiological quantities, such as the distribution of final epidemic size and expected epidemic duration, have been derived and estimated based on real data. Pandemics have largely shaped the history of human being as described in the popular book by McNeill (1976), and have made huge impacts to our society and and economy. However, mathematical models developed in epidemiology and economics don’t talk to each other much until the current outbreak of COVID-19, which needs urgent calls (e.g. from The Royal Society) for researchers across disciplines to work together and jointly support the scientific modelling for epidemics, see recent intensive interplays between the two fields, e.g. Acemoglu et al. (2020); Alvarez et al. (2020); Atkeson (2020); Eichenbaum et al. (2020) and Guerrieri et al. (2020). Our paper also responds to the calls by introducing the Hawkes-based approach as a potentially very valuable tool for epidemic modelling.

This paper is organised as follows: Section 2 offers the preliminaries including an introduction and formal mathematical definitions for our stochastic epidemic model, a two-phase dynamic contagion process. Key distributional properties, such as the distribution of final epidemic size and expected epidemic duration, are provided in Section 3. In Section 4, our model is implemented based on real data, and the associated time lag of the effect of intervention in each country or region is estimated. Finally, Section 5 draws a conclusion for this paper, and proposes some issues for possible further extensions and future research.

2 Two-Phase Dynamic Contagion Process

In this section, we introduce a two-phase dynamic contagion process (2P-DCP) for modelling the dynamics of COVID-19 contagion. The unobservable effective time that aggregated government interventions (e.g. lockdown of a city or country) came into effect is denoted by ℓ > 0, which divides the COVID-19 epidemic dynamics into two phases. Note that, the time point is different from the exact timing of intervention that can be observed. The cumulated number of infected individuals is described by a counting process Nt with N0 = 0, and it is modelled by a two-phase dynamic contagion process defined as below.

Definition 2.1 (Two-Phase Dynamic Contagion Process (2P-DCP)). A two-phase dynamic contagion process (2P-DCP) is a point process Nt with two phases:

Phase 1 (Full Contagion):
For the first phase period t [0,ℓ), Nt follows a dynamic contagion process with stochastic intensity
                *
         -δt  ∑Nt    -δ(t- T*)  N∑t    -δ(t- T)
λt =  λ0e   +     Zke      k +    Yie      i,    t ∈ [0,ℓ],
              k=1              i=1
(2.1)

where

  • λ0 > 0 is the initial intensity at time t = 0;
  • δ > 0 is the constant rate of exponential decay;
  • Nt*≡{Tk*}k=1,... is a Poisson process of constant rate ϱ > 0 arriving in time t ;
  • {Zk}k=1,...,N* are i.i.d. externally-exciting jump sizes, realised at times {Tk*}k=1,...,N*, with distribution H(y);
  • {Y i}i=1,...,N are i.i.d. self-exciting jump sizes of the first phase, realised at times {Ti}i=1,...,N , with distribution G1(y).
Phase 2 (Self Contagion):
For the second phase period t (ℓ,), Nt is a pure self-exciting Hawkes process with stochastic intensity
                  ∑Nt
λt =  λℓe-δ(t-ℓ) +      Yie-δ(t- Ti),    t ∈ (ℓ,∞ ),
                 i=Nℓ+1
(2.2)

where

  • λ is the initial intensity of the second phase starting at the cutoff time point , which is the terminal intensity of the first phase;
  • {Y i}i=N+1,... are i.i.d. self-exciting jump sizes of the second phase, realised at times {Ti}i=N+1,..., with distribution G2(y). Note that compared with the average mean of the self-exciting jump size for Phase 1, the mean of self-exciting jump in Phase 2 is smaller, which demonstrates that the COVID-19 becomes less contagious on average after time point .


PIC

Figure 1: Two-phase dynamic contagion process


The point process Nt and its intensity process λt are illustrated in Figure 1. Overall, we can more compactly define our new pandemic model, a two-phase dynamic contagion process, as a counting process Nt ≡{Ti}i=1,... with N0 = 0 and stochastic intensity

              N *t                   Nt
λ  =  λ e-δt + ∑ Z  1    e-δ(t- Tk*) + ∑ Y e-δ(t- Ti),    t ≥ 0,
 t     0           k {t≤ℓ}               i
              k=1                   i=1
(2.3)

where

This equivalent definition as a dynamic contagion process has an advantage: it can be viewed as a branching process and has a more intuitive interpretation with regard to a pandemic. The cluster-process presentation is provided as follows.

3 Distributional Properties

In this section, we outline key distributional properties for the two-phase dynamic contagion process. We derived the conditional joint Laplace transform of λt and probability generating function of Nt, which are the key results to further derive the elimination probability of the epidemic and the distribution of the final epidemic size.

Joint Distribution of (λt,Nt)

Let {Ft}t0 be the natural filtration of the point process Nt, i.e. Ft = σ(Ns,s ≤ t) and assume that the intensity process {λt}t0 being Ft-adapted. The joint Laplace transform and probability generating function for (λt,Nt) is provided in Theorem 3.1 as below.

Theorem 3.1. For time s t, the conditional joint Laplace transform and probability generating function for λt and the point process Nt is of the form

                 8
                 ><   Ns c(s)-A (s)λs
E   θNte-vλt|F    =   θ  e   e       ,   0 ≤ s ≤ t ≤ ℓ,            (3.1)
             s   >:  Ns - A(s)λs
                   θ  e       ,       ℓ < s ≤ t,
where A(s) is determined by the nonlinear ordinary differential equation (ODE)
A ′(s)- δA (s)+ 1 - θˆg(A(s);s) = 0,
(3.2)

where the boundary condition is A(t) = v with

ĝ(u;t) = 0e-uydG(y;t),

and c(t) is determined by

             ∫t
                     ˆ
c(t) = ϱ1{t≤ℓ}    1 - h(A(u)) du,
             0
(3.3)

with

ĥ(u) = 0e-uydH(y).

Proof. For 0 s t , λt is the intensity process of the dynamic contagion process introduced in Dassios and Zhao (2011). The corresponding conditional joint Laplace transform, probability generating function for the process λt and the point process Nt is provided in Theorem 3.1 of Dassios and Zhao (2011). For ℓ < s t, given Fs, the infinitesimal generator of the dynamic contagion process (λs,Ns,s) acting on a function f(λ,n,s) within its domain Ω(A) is given by


                              ∞∫
Af (λ,n,s) = ∂f-- δλ ∂f-+ λ     f(λ + y,n + 1,s)dG (y;s)- f (λ,n,s)  .
             ∂s      ∂λ
                              0
(3.4)

Consider a function f(λ,n,s) of form

f(λ,n,s) = θne-A(s)λ,

and substitute this into Af = 0, we then have the ODE

A(s) = δA(s) - 1 + θĝ(A(s),s),

adding the boundary condition A(t) = v, gives the ODE in (3.2). __

The moments of λt and Nt can be obtained by differentiating the joint Laplace transform and probability generating function of λt and Nt, and the results are provided in Proposition 3.1 and 3.2 as below.

Proposition 3.1. The conditional expectation of the process λt given Fs for s t is given by

           8> ϱμ 1              ϱμ  1
           <  --H-{κt≤-ℓ} +  λs - -H-κ{t≤ℓ} e-κ(t-s),    κ ⁄= 0,
E [λt|Fs] = >                                                        (3.5)
           : λs + ϱ μH1 {t≤ℓ}(t- s),                κ = 0,
and the conditional expectation of the point process Nt given Fs is of the form
           8
           ><       ϱμH1{t≤-ℓ}(t-s)          ϱμH1-{t≤ℓ}  1-e-κ(t-s)
E[N |F ] =   Ns +       κ      +  λs -    κ         κ    ,    κ ⁄= 0,    (3.6)
    t s    >:                 1               2
             Ns + λs(t- s)+  2ϱμH 1{t≤ℓ}(t- s) ,               κ = 0,
where
μH = 0ydH(y), μ G = 0ydG 1(y)1{t} + 0ydG 2(y)1{t>ℓ},

and κ = δ -μG.

Proof. The result in (3.5) immediately follows Theorem 3.6 in Dassios and Zhao (2011). And since Nt -Ns - stλ udu is a martingale, we have

E[Nt - Ns |Fs] = stE[λu|Fs ]du,

which directly implies the result in (3.6). __

Proposition 3.2. The conditional second moment of the process λt given Fs for s t is given by

           8
           >>>   λ2e-2κt + 2ϱμH1{t≤ℓ}+μ2G  λs - ϱμH   e-κ(t-s) - e-2κ(t-s)
           >>>>    s              κ              κ
           >><     +  (2ϱμH+-μ2G)ϱμH1-{t≤ℓ}+ ϱμ2H1{t≤ℓ}    1 - e-κ(t-s)  ,       κ ⁄= 0,
E   λ2t|Fs   =                 2κ2              2κ                                (3.7)
           >>>   λ2+ λ μ   t+ (2λ μ  + μ   )ϱ1     (t - s)
           >>>>    s    s 2G       0 H     2H    {t≤ℓ}
           >>:    +   ϱ2μ2 + 1ϱμ  μ     1    (t- s)2,                      κ = 0,
                      H   2   H 2G   {t≤ ℓ}
and the conditional joint expectation of the process λt and the point process Nt given Fs for s > ℓ is of the form
             8
             >>       - κ(t-s)            λsμ2G-         -κ(t-s)
             >>>  λsNse       +  λsμG +    κ   (t- s)e
             <        λ2s   λsμ2G     - κ(t-s)   -2κ(t- s)
E[λtNt|Fs] = >>    +  -κ - --κ--  (e       - e       ),         κ ⁄= 0,     (3.8)
             >>>            2                  1            2
             : λsNs +  λs + λsμG  (t - s)+ 2λsμ2G (t- s) ,    κ = 0
and the second moment of point process Nt given Fs for s > ℓ is of the form
            8
            >>>   λsμ2G - λ2s  (1- e-2κ(t- s))-   2λsμG-+  2λsμG2- (t- s)e-κ(t-s)
              >><     κ3     κ2                      κ      κ2
E N 2t |Fs =     +   λs+2λsNs-+ 2λ2s+2λsμ2G- (1-  e-κ(t-s)),                      κ ⁄= 0,(3.9)
            >>>           κ          κ2
            >>: λ (μ  + 2N  )(t-  s) +   λ2+ λ  μ   (t-  s)2 + 1λ μ  (t - s)3,      κ = 0,
               s  G      s           s    s G           3 s 2G
where
μ2H = 0y2dH(y), μ 2G = 0y2dG 1(y)1{t} + 0y2dG 2(y)1{t>ℓ},

and κ = δ -μG.

Proof. These results immediately follows Lemma 3.1, Lemma 3.2 and Theorem 3.9 in Dassios and Zhao (2011). __

Probability of Elimination Time eT

After the government interventions come into effect, the contagion rate will dramatically decline and new cases will drop abruptly almost to nothing in the near future. It is therefore of great interests to calculate the probability of elimination time, i.e. the time that the last ever case arrives, after government interventions come into effect. Let eT to be elimination time such that

Te:= inf {t > ℓ : ∀s ≥ t, Ns - Nt = 0}.
(3.10)

The condition probability of the elimination time is provided in Proposition 3.3.

Proposition 3.3. For s t, the elimination probability is given by

                  -A(s)λ
ℙ  eT ≤ t|Fs  = e     s,
(3.11)

where A(s) is determined by the ODE in (3.2) with boundary condition A(t) = 1
δ.

Proof. Given e
T being the timing of the last ever event, the event { e
T t} implies that Nu = Nt for any u t, which also lead to

λu = e-δ(u-t)λ t.

Hence, we have


ℙ Te≤  t|Fs   =   E 1 {Te≤t} | Fs
                   2         ∞∫                    3
                   4             -δ(u- t)         5
             =   E  exp   -   λte      du   | Fs
                              t
                          λt
             =   E  exp  -  δ  | Fs                          (3.12)
And according to Theorem 3.1, by setting θ = 1 in (3.1), the result follows immediately. __

Joint Expectation of Epidemic Size Nt and Elimination Time  e
T

Given the last ever event {Te < t}, one could obtain the expected size of the epidemic at time t. The relevant details are presented in Corollary 3.1.

Corollary 3.1. For s e
T t, the conditional joint expectation of Nt and 1{eTt}is of the following form

                      d    Ns - A(s)λs
E  Nt1 {Te≤t} | Fs = dθ- θ   e       |θ=1- ,                (3.13)
where A(s) satisfies the ODE in (3.2) with boundary condition A(t) = 1δ.

Proof. According to Theorem 3.1 and Proposition 3.3, we have

    Nt                   h Nt - λδt    i
E θ  1 {Te≤t} | Fs  =   E θ   e   | Fs
                  =   θNse-A(s)λs.                     (3.14)
Since E
Nt1      | Fs
    {Te≤t} = d-
dθE
θNt1      | Fs
    {Te≤t}

θ=1- , the result immediately follows (3.13). __

Distribution of Final Epidemic Size N

The final epidemic size is one of the most important epidemiological quantities to study. In fact, under the two-phase dynamic contagion model, the final epidemic size is the value of the point process Nt when time goes to infinity. Conditional on s > ℓ, since there are no externally-exciting jumps in the intensity, the distribution of N can be characterised by Proposition 3.4 as below.

Proposition 3.4. For ℓ < s, the probability generating function of Nconditional on Fs is given as

E   θN∞ | F   = e- v*λs,
          s
(3.15)

where

v* = 1-
δ      ∫∞
1 - θ   e-v*ydG (y)
               2
      0.

Proof. The result immediately follows Theorem 3.5 in Dassios and Zhao (2011). __

While the government interventions come into effect, if we assume i.i.d. self-exciting jump sizes Y i ~ Exp(β) for i = N + 1,...., then, we have an explicit expression for the probability generating function of N as

E[θN Fs] = exp   p ---------------------
- --(δβ---1)2 +-4δβ-(1---θ)--(δβ---1)λ
                  2δ                s, s > ℓ.

This implies that, the final epidemic size N conditional on F follows a mixed-Poisson distribution with the probability mass function

                  ∞∫  k - v
ℙ (N∞  = k | Fs) =  v-e---m (v)dv,     k = 0,1,....,
                      k!
                  0
(3.16)

where m(v) is the density function of the mixing distribution,

m(v) := exp
δβ---1       δβ---1 2 δ-   2βδλ2s-
  2δ   λs -    2δ     βv -  2v   ---
--2βδλs-
√ ---32
  2πv,

which is an inverse Gaussian distribution with parameters --β-
δβ-1λs and β-
2δλs2.

4 Empirical Study

We provide a calibration scheme based on the daily increments of the two-phase dynamic contagion process Nt. Let us first denote the observations of the daily confirmed COVID-19 cases as {Ct}t=0,1,2,...,T . The mean square error (MSE) between the expected daily increments of Nt and the actual reported daily confirmed COVID-19 cases is given as

                    T-1
                  1-∑                       2
MSE (α,β,δ,ϱ,ℓ) = T    (E [Nt+1 - Nt ]- Ct+1) .
                    t=0
(4.1)

We consider the calibration based on minimising the MSE (4.1), i.e., we choose parameters (αˆ,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) such that

MSE(ˆα,ˆ
β , ˆ
δ ,ϱˆ ,ˆ
ℓ )  :=  min(α,β,δ,ϱ,ℓ)MSE(α,β,δ,ϱ,ℓ),

with α,β,δ,ϱ 0 and +.

Without loss of generality, for simplicity, we assume Zk = 1 for any k, Y i ~ Exp(α) for i = 1,...,N and Y i ~ Exp(β) for i = N + 1,..., and λ0 = 0 in (2.3) for model calibration. Other assumptions for Zk,{Y i}i=1,...N,N+1,...0 can also be used if necessary. We provide two empirical examples. The first one concentrates on the COVID-19 pandemic in mainland China during the period from early January to late March 2020. The second one focuses on the worldwide COVID-19 pandemic during the period from mid February to early May 2020. The data we used are publicly available. Datasets are mostly cited from the associated official government health department websites for non-European countries and from the European Centre of Disease Control (ECDC) for European countries.

COVID-19 Pandemic in Mainland China

The daily confirmed COVID-19 cases for regions in China can be obtained from daily reports of the National Health Commission of the PRC. We use the reported daily confirmed cases for several regions in China from 2020-01-19 to 2020-03-31 as examples for model calibration. The corresponding estimation results are illustrated in Table 1. We can see that the estimator ˆϱ are quite different from each other. In particular, these regions that are close to Hubei, namely Henan, Hunan, Anhui, have relatively larger intensities for externally-exciting jumps, which means that these regions experienced more external shocks from Wuhan and these external shocks can be originally infected individuals from Wuhan. Naturally, without taking into account Wuhan, the rest of Hubei has the largest intensity ˆϱ for externally-exciting jumps. The main government intervention established by the Chinese authority is the announcement of the completely lockdown of Wuhan and later the whole Hubei Province on 23 January 2020. One or two days later, all other regions enforced the quarantine and raised the alert of public health emergency. Setting the date 2020-01-19 as the initial time t = 0, then the government intervention took place when t = 4 and came into effect when t = ˆℓ . Since the delay period of the government interventions is the difference of the date when the government interventions came into effect and the date when the government introduced the restriction measures, we can therefore observe from Table 1 that the estimated delays of the government interventions for different regions therefore are between 5 and 15 days, which are consistent to the incubation time of COVID-19 for most people, e.g. as found in a highly cited medical study of Lauer et al. (2020).



Table 1: Calibration parameters (ˆα,ˆ
β ,ˆ
δ ,ˆϱ ,ˆ
ℓ ) for total confirmed COVID-19 cases from 2020-01-19 to 2020-03-31 for various regions in China.

_________________________________________________________________________________








PICT                                                                           ˆα                                                                                   βˆ                                                                                            ˆδ                                                                                                    ˆϱ                                                                                                          ˆℓ                                                                                                                 ˆˆ ˆ
                                                                                                            MSE(αˆ,Nβ,δ,ˆϱ,ℓ)














Heilongjiang 2.4311238.0361460.2523240.93457416 0.036414







Sicuan 5.5092847.4028510.21608 3.91349611 0.023037







Shandong 5.3034777.8680650.3402735.29646519 0.033823







Jiangxi 2.9365127.9036770.2859223.78210415 0.070023







Anhui 3.6441567.7739380.3391837.16124418 0.037944







Hunan 4.4679527.6603230.2938608.35100815 0.053281







Zhejiang 1.3987418.0651250.2604491.62782209 0.144322







Henan 3.3709327.8677860.2863576.66910515 0.042004







Guangdong 2.3168797.9879040.2561432.36050712 0.073389







Hubei 4.46843528.467800.17387458.2400116 0.654813







_________________________________________________________________________________

The branching ratio (BR), which demonstrates the average infection rate, is determined by E[Y i]∕δ. In Table 2, we compare the estimated branching ratios before and after the government interventions came into effect, namely Rb and Ra, respectively. It is clear that the branching ratio for every region decreases significantly when the state changed, i.e. government interventions came into effect. One can also access the efficiency for when regions implemented the restriction packages introduced by the central government by comparing the corresponding branching ratios Rb and Ra. The comparison of Rb and Ra for regions in China are presented in Figure 2. We can see that the government restriction packages had been well-implemented for all regions in China. In particular, Hubei, where the strictest measure, i.e. the completely lockdown of Wuhan and Hubei, had been introduced, shows a dramatic drop of contagion rate after the interventions came into effect.



Table 2: The estimated branching ratio (BR) before and after the government interventions came into effect, namely Rb,Ra respectively, for regions in China.

_________________________________________________________________




PICT Rb Ra






Heilongjiang 1.6301760.493166



Sicuan 0.8400210.625153



Shandong 0.5541300.373512



Jiangxi 1.1910250.442510



Anhui 0.8090390.379250



Hunan 0.7616410.444234



Zhejiang 2.7449890.476066



Henan 1.0359570.443853



Guangdong 1.6850520.488747



Hubei 1.2870940.202028



_________________________________________________________________


PIC

Figure 2: Comparison of the branching ratios before and after the government interventions came into effect for regions in China. The horizontal axis represents the branching ratio before the government interventions came into effect, namely Rb and the vertical axis represents the branching ratio after the government interventions came into effect, i.e. Ra.


Figure 3 and 4 demonstrate comparisons between the expected daily/total confirmed cases for the two-phase dynamic contagion model under the calibrated parameter (ˆα,βˆ ,ˆδ ,ˆϱ ,ˆℓ ) in Table 1 and the actual daily/total confirmed COVID-19 cases for the period 2020-01-19 to 2020-03-31. We observe that the model allows different shapes of trend before the interventions came into effect. All these regions indicate relatively smooth exponential decay of daily new cases after the peak. In addition, the estimated cumulative confirmed cases are very close to the actual total confirmed COVID-19 cases, which further confirms that our new model is a good candidate for describing the propagation process.


PIC

Figure 3: Model calibration comparisons between the expected daily confirmed cases under the calibrated parameters (ˆα,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) in Table 1 and actual daily confirmed COVID-19 cases from 2020-01-19 to 2020-03-31 for Guangdong, Henan, Hunan, Anhui, Jiangxi, and Hubei of China.



PIC

Figure 4: Comparisons between total confirmed COVID-19 cases and total estimated cases under the calibrated parameters (αˆ,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) in Table 1 from 2020-01-19 to 2020-03-31 for Guangdong and Hubei, China.


COVID-19 Pandemic for the World

From mid-February 2020, the COVID-19 started to spread in other countries around the world. At beginning, only a small number of initial cases were reported for some countries in Europe, South/East Asia and North American. However, lately, several large outbreaks were reported in South Korea, Italy, Iran, Spain, Japan and the total number of cases outside China quickly passed the China’s. The WHO then recognized the spread of COVID-19 as pandemic on 2020-03-11. We could use this as a second example to confirm our observations from the last exercise. The calibration settings were the same as the previous one. We use the reported daily confirm cases for different regions and countries around world from mid-February to early May 2020. Note that, due to the fact that the pandemic reached each country or territory at different time and the corresponding government interventions also imposed and came into effect at different times, there is no sense to calibrate the model using the data within the same truncated time series. Table 3 presents the estimation results (ˆα,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) of (α,β,δ,ϱ,ℓ) for various countries and territories. We notice that regions and countries like Italy, China, New York have much larger ϱ compared with other areas. This phenomena is reasonable as these areas have specific outbreak area which created external shocks to other part of the regions and countries and hence the number of confirmed cases increased more rapidly than other regions and countries.



Table 3: Calibration parameters (ˆα,ˆ
β ,ˆ
δ ,ˆϱ ,ˆ
ℓ ) for total confirmed COVID-19 cases from mid-February, 2020 to early May, 2020.

_________________________________________________________________________________










PICT Date0 DateG                                                                                               ˆα                                                                                                       ˆβ                                                                                                               ˆδ                                                                                                                        ˆϱ                                                                                                                              ˆℓ                                                                                                                                     ˆ ˆ ˆ
                                                                                                                                 MSE(ˆα,βN,δ,ˆϱ,ℓ)-


















Australia 2020-02-272020-03-151.7038073.3160830.4015500.489762 28 0.124391









Austria 2020-03-012020-03-102.8602917.7819220.222730 5.788434 24 0.326694









China(Mainland) 2020-01-192020-01-232.8462366.5342690.24277383.16086517 1.003152









Czech 2020-03-042020-03-103.4646547.8428620.1749902.247588 25 0.359398









France 2020-02-252020-03-133.4534596.9124100.19189015.49253536 5.410160









Germany 2020-02-242020-03-123.2467986.4020970.20132722.86372732 3.464682









Greece 2020-02-262020-03-104.1818317.7781440.2008861.337306 35 0.172176









Hong Kong 2020-03-012020-03-232.7823698.0544020.2870480.512120 31 0.041714









Iceland 2020-02-282020-03-133.4059017.9086910.2753781.892734 33 0.163572









Italy 2020-02-212020-03-053.0681035.4456460.21104823.83043130 0.911767









Latvia 2020-03-072020-03-134.2273887.8192750.2082091.145524 22 0.099024









New York 2020-02-292020-03-122.4499063.2035630.35659270.93760433 2.926862









New Zealand 2020-03-122020-03-162.1229998.4386490.2046670.551023 14 0.064829









Norway 2020-02-212020-03-124.4965017.5195730.1958596.553447 29 0.303364









South Korea 2020-02-162020-02-201.3095078.7719680.2443165.641225 13 0.418018









Switzerland 2020-02-252020-03-133.1958297.0468640.20261514.08374928 1.571950


















_________________________________________________________________________________

In Table 3, we also presents the date of day 0, i.e. Date0, and the date of government interventions imposed, i.e. DateG. The delay period of government interventions came into effect therefore can be obtained given the estimated ˆ
ℓ , with Date0 and DateG. The details for the delay periods of regions and countries are illustrated in Figure 5. We can see that the delay of the interventions for different regions and countries is around 8 ~ 21 days. In fact, the delay period can be considered as a criteria for evaluating the effectiveness of restrictions imposed by the authorities to prevent further spread of COVID-19. In general, most regions and countries with short delay periods normally took tougher restrictions or more effective measures to stop the spread of virus. New Zealand and South Korea are two typical examples. The authority of New Zealand introduced a nationwide lockdown by closing all borders and entry ports to all nonresidents. On contrary, the South Korea authority introduced one of the largest and best-organised epidemic control programs to screen the mass population for the virus with isolation, tracing, quarantine took place simultaneously without further lockdown. Therefore, the estimated delay periods for these two countries are only 9 days for New Zealand and 8 days for South Korea, which are much shorter than the average incubation time of COVID-19. For most European countries, due to the containment restriction measures such as quarantines and curfews were not strictly put into effect, the associated delay periods are relatively longer than the incubation time of COVID-19.

The estimated branching ratios before and after the government interventions came into effect for different countries and territories are reported in Table 4. And Figure 6 demonstrates a comparison between the BRs before the government interventions took effect and the BRs after the interventions worked, with a blue dash line that represents Rb = Ra. We can immediately see that for most regions and countries, the BR dropped dramatically after government interventions came into effect, which suggests that the restriction measures imposed by the authority indeed reduce the contagion/infection rate significantly.


PIC

Figure 5: Comparison of the delay period for different regions and countries around the world. The horizontal axis represents the abbreviation of regions and countries listed in Table 4 and the vertical axis represents the number of days for the government interventions came into effect after the government announcement the relevant measures.



PIC

Figure 6: Comparison of the branching ratios before and after the government interventions came into effect for different regions and countries around the world. The horizontal axis represents the branching ratio before the government interventions came into effect, Rb, and the vertical axis represents the branching ratio after the government interventions came into effect, Ra.




Table 4: The estimated branching ratio (BR) before and after the government interventions came into effect, namely Rb,Ra respectively, for regions and countries around the would.

_________________________________________________




PICT Rb Ra






Australia 1.4616380.750991



Austria 1.5696810.576945



China(Mainland) 1.4472000.630380



Czechia 1.6493990.728637



France 1.5090160.753908



Germany 1.5298260.775845



Greece 1.1903770.639993



Hong Kong 1.2520770.432526



Iceland 1.0662000.459162



Italy 1.5443640.870102



Latvia 1.1367980.614084



New York 1.1446670.875377



New Zealand 2.3014550.579001



Norway 1.1354880.678991



South Korea 3.1256430.466606



Switzerland 1.5443490.700379






_________________________________________________

A comparison between the expected daily confirmed cases for the two-phase dynamic contagion model under the calibrated parameters (αˆ,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) is reported in Table 3, and the actual daily confirmed COVID-19 cases for different regions and countries over the period of mid-February to early May are presented in Figure 7. In general, we can see that the model can precisely catch the trend of infection, this further confirms that the two-phase dynamic contagion model is effective. Note that, we have smoothed the biggest jump of daily confirmed cases in China for better illustration and fitting purpose. This is due to a change in the confirmation standard established by the Chinese authority on 2020-02-12.


PIC

Figure 7: Model calibration comparisons between the expected daily confirmed cases under the calibrated parameters (ˆα,βˆ ,δˆ ,ϱˆ ,ℓˆ ) in Table 3 and actual daily confirmed COVID-19 cases for Australia, Austria, China(Mainland), Italy, New Zealand, and South Korea.


In Figure 8, we compare the estimated daily increment with the actual confirmed COVID-19 cases for France, Germany, Switzerland, and New York. The daily records of confirmed cases for these areas were not as smooth as those countries illustrated in Figure 7. The spikes within the graphs could be caused by many reasons such as the delay of reports, testing capacity, hospital capacity, diagnostic methods and etc. For instance, the daily confirmed cases for France, Germany, Switzerland and New York suddenly declined on a regular basis, which mostly happened during the weekends. Even so, we can see the model can still capture the trend of infectious evolution. In Figure 9, we also compare the actual total confirmed COVID-19 cases against the cumulative estimated cases with a confidence interval within two standard deviations1 for some typical countries and territories. The black dash line in each graph of Figure 9 represents the end of data collection period for calibration. The red curve on the left of the black dash line shows the historical data used for calibration and the blue curve is the corresponding estimated result. The red and blue curves on the right of the black dash line demonstrate a comparison between the predicted and actual confirmed COVID-19 cases for countries and regions from the end of their data collection period to the end of May, 2020. We can see that the estimated curves of the number of confirmed infections under the two-phase dynamic contagion model well fitted the actual propagation process of the COVID-19. In addition, the forecasted infection cases in the coming weeks after the end of data collection period also well suited the up to date actual total confirmed COVID-19 cases.


PIC

Figure 8: Model calibration comparisons between the expected daily confirmed cases under the calibrated parameters (ˆα,βˆ ,δˆ ,ϱˆ ,ℓˆ ) in Table 3 and actual daily confirmed COVID-19 cases for France, Germany, Switzerland, New York.



PIC

Figure 9: Comparisons between total confirmed COVID-19 cases and total estimated cases under the calibrated parameters (ˆα,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) in Table 3 from mid-February onward for New Zealand, Austria, Germany, France, Italy, and New York. The red curve represents total confirmed COVID-19 cases. The blue curve represents the total estimated cases and the left zone of the black dash line illustrates historical data used for calibration, and the right zone demonstrates the predicted estimated cases. The shadowed region plots the values within two standard deviations.


4.1 Elimination Probability and Final Epidemic Size

According to Proposition 3.3, one could obtain the elimination probability of the epidemic by numerically solving the ODE (3.2). Based on the calibration parameters (ˆα,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) provided in Table 3 for regions and countries, we could obtain the associated elimination probabilities. Figure 10 illustrates how  e
T ≤ t|Fˆℓ varies for different countries and territories. We can see that for regions and countries with effective restriction measures, the probability for a shorter period to observe the last ever event arrives after government interventions come into effect will be much higher. On contrary, for some regions and countries, longer periods are needed for elimination probabilities to be closed to 1. For instance, we can see that there is still a long way to go to end the COVID-19 pandemic for Italy.


PIC

Figure 10: Conditional elimination probability (eT t|Fℓˆ ) for China (Mainland), New Zealand, South Korea, Germany, Italy, New York, Iceland and Hong Kong under the associated calibration parameters (ˆα,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) for these regions and countries suggested in Table 3.


The elimination time of the pandemic depends on many decisive factors, such as the initial intensity of the externally-exciting jumps, the time needed for the government interventions to come into effect, the size of the branching ratio after the government interventions came into effect, etc. Figure 11, 12, and 13 illustrate comparisons between the estimated  ˆ
ℓ , ϱˆ , Ra against E[e
T|Fˆℓ ], respectively. From Figure 11, we can see that for most countries and territories, the quicker the government interventions come into effect, the faster the pandemic will end. However, some places like Hong Kong and Iceland still have relatively fast elimination time even though it takes longer for the government interventions to come into effect. This is probably because the restriction measures for these places were imposed so early that reporting procedures were not properly in place yet. Figure 12, and 13 clearly demonstrate that the externally-exciting jump intensity ϱ and the branching ratio after the government came into effect Ra are important factors that determine the extinction time of the pandemic. In general, to reduce the extinction time of the pandemic, the first priority for the authorities should be introducing restriction measures such as national/subnational lockdown to reduce the intensity of the external imported cases, and while the external imported cases are controlled and thereafter negligible, the governments should simultaneously introduce enforced restrictions to prevent further transmission. If the government intervention strategies were effectively implemented without a lack of civic spirit, the infection rate of the virus after the these intervention measures come into place will be reduced significantly, and therefore lead to a quicker elimination of the COVID-19 pandemic. Note that the prediction for expected elimination time for regions and countries is based on the assumption that the government intervention measures are still taking into place in some form and propagation of the disease continues as in Phase 2. Relaxation of the government intervention measures will inevitably delay the disease elimination for most regions and countries.


PIC

Figure 11: Comparison of the conditional expectation of the elimination time Te and the estimated government interventions came into effect time for different regions and countries around the world. The horizontal axis represents ˆℓ , and the vertical axis represents E[eT|Fˆℓ ] .



PIC

Figure 12: Comparison of the conditional expectation of the elimination time Te and the estimated intensity for externally-exciting jumps for different regions and countries around the world. The horizontal axis represents ˆϱ , and the vertical axis represents E[Te|Fˆℓ ] .



PIC

Figure 13: Comparison of the conditional expectation of the elimination time eT and the branching ratio after the government interventions came into effect for different regions and countries around the world. The horizontal axis represents Ra, and the vertical axis represents E[eT|Fˆℓ ].


Beside the conditional probability for the elimination time  e
T, the epidemic size Nt given {˜
T t} can also be predicted according to the join expectation of Nt and {Te t} derived in Corollary 3.1. In Table 5, we report the 95% confidence interval for elimination time eT the condition expectation of the elimination time e
T, E[ e
T|Fˆℓ ], the expected elimination date DateE, and the conditional expectation of the epidemic size Nt, E[Nt|Fˆℓ ∩{Te t}] with t = E[eT|Fˆℓ ], for regions and countries with calibration parameters in Table 3. Note that, the regions and countries with more confirmed COVID-19 cases before government interventions came into effect will experience longer time to reach elimination state, like France, Germany, Italy, and New York. And the corresponding expected epidemic size for these areas are also much larger. Note that since we have smoothed the biggest jump of daily confirmed cases, adding up with the cases which have been smoothed, the actual conditional expectation of the epidemic size is about 83113, which is very close to the current total confirmed cases 82,993 on 2020-05-27. In general, not only the expected epidemic size is very close to the actual total confirmed cases for the listed regions and countries, but also the estimated elimination date is very close to the actual eradicate date. New Zealand is one typical example that can be used to demonstrate the effectiveness of our model in predicting the key epidemiological quantities. New Zealand authority has officially declared that the country has completely eradicated COVID-19 for now with total of 1154 confirmed COVID-19 cases on 2020-06-08. This elimination date and the final epidemic size are very close to what we predicted for New Zealand under the two-phase dynamic contagion model, the predicted elimination date is around 2020-06-04 and the predicted epidemic size is about 1250. More remarkably, the historical data we used for model calibration for New Zealand is from 2020-03-12 to 2020-04-13, which is already one and half months ago. This clearly shows that the two-phase dynamic contagion model is pretty powerful in forecasting cumulative confirmed COVID-19 cases, predicting possible elimination duration for the pandemic, and evaluating effectiveness of relevant government intervention measures.



Table 5: The 95% confidence interval for elimination time e
T, the conditional expectation of the elimination time e
T, the expected elimination date and the conditional expectation of the epidemic size Nt given eT t with t = E[eT|Fˆℓ ] under the calibrated parameters (αˆ,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) for these regions and countries suggested in Table 3.

_________________________________________________________________________________






PICT(                                                                         eT (t1,t2)|F                                                                                    ˆℓ ) = 95%E[                                                                                              eT|F                                                                                                 ˆℓ ]      DateE E[Nt|F                                                                                                                        ˆℓ ∩{                                                                                                                           eT t}]










Australia (47,97) 66 2020-05-30 7059.46





Austria (69,122) 89 2020-06-21 15610.22





China (Mainland) (87,142) 108 2020-05-22 69675.10





Czechia (111,216) 150 2020-08-25 8932.448





France (166,272) 206 2020-10-23 154112.00





Germany (175,285) 216 2020-10-28 189116.50





Greece (59,127) 85 2020-06-24 2650.07





Hong Kong (23,54) 35 2020-05-05 945.81





Iceland (28,61) 41 2020-05-11 1787.79





Italy (264,445) 329 2021-02-13 276480.1





Latvia (40,102) 64 2020-05-31 764.80





New York (149,261) 192 2020-10-10 209311.80





New Zealand (49,107) 71 2020-06-04 1249.77





Norway (83,162) 113 2020-07-15 8024.73





Switzerland (121,203) 152 2020-08-22 33414.19





_________________________________________________________________________________

The conditional distribution of the final epidemic size N can be obtained by numerically inverse the probability generating function provided in Proposition 3.4. Since we assume the self-exciting jumps follows an exponential distribution after government interventions came into effect, the final epidemic size N conditional on F follows a mixed-Poisson distribution with probability mass function specified in (3.16). Figure 14 demonstrates the conditional probability mass function of the difference between the finial epidemic size N and the total number of confirm cases N when government interventions came into effect, i.e., (N ∞ - N ℓ = k | F ℓ), for some regions and countries under the calibrated parameters (ˆα,ˆβ ,ˆδ ,ˆϱ ,ˆℓ ) for these regions and countries suggested in Table 3.


PIC

Figure 14: Probability mass function 
N ∞ - N ˆℓ = k | F ˆℓ for Latvia, New Zealand, Hong Kong and Iceland under the calibrated parameters (ˆα,βˆ ,ˆδ ,ϱˆ ,ℓˆ ) for these regions suggested in Table 3.


5 Concluding Remarks

In this paper, we have introduced a two-phase dynamic contagion process for modelling the current spread of COVID-19. This model allows randomness to the infectivity of individuals rather than a constant reproduction number as commonly assumed by standard models. Key episdemiological quantities, such as the distribution of final epidemic size and expected epidemic duration, are derived and estimated based on real data for various regions and countries. In addition, the associated time lag of the effect of intervention in each country or region has been estimated, and our empirical results are consistent to the incubation time of COVID-19 for most people found by existing medical study such as Lauer et al. (2020). The aim of this paper is to demonstrate that our model, as a representative of Hawkes-based processes, could be a valuable tool for epidemic modelling. In fact, the vast literature of Hawkes-based processes would also be relevant and potentially be applicable. For example, multivariate extensions of Hawkes-based processes, such as Embrechts et al. (2011) for analysing financial high-frequency data, could be adopted for modelling the cross-region epidemic contagion. Lévy-driven extensions, such as Qu et al. (2019) for portfolio credit risk, may perform better in capturing the heavy tail property of epidemic distribution. In addition, easing of the government interventions will lead to change of parameters and delay extinction times. The model can be adjusted by introducing an additional phase with change of parameters. When countries cycle between periods of restrictions and relaxations to manage COVID-19, we can adjust the two-phase dynamic contagion model by replacing the constant parameters to piecewise time dependent parameters. These are all proposed as future research.

References

   Acemoglu, D., Chernozhukov, V., Werning, I., and Whinston, M. D. (2020). A multi-risk SIR model with optimally targeted lockdown. National Bureau of Economic Research Working Paper.

   Aït-Sahalia, Y., Cacho-Diaz, J., and Laeven, R. J. (2015). Modeling financial contagion using mutually exciting jump processes. Journal of Financial Economics, 117(3):585–606.

   Allen, L. J. (2008). An introduction to stochastic epidemic models. In Allen, L. J., Brauer, F., Van den Driessche, P., and Wu, J., editors, Mathematical Epidemiology, chapter 3, pages 81–130. Springer, Berlin.

   Alvarez, F. E., Argente, D., and Lippi, F. (2020). A simple planning problem for COVID-19 lockdown. National Bureau of Economic Research Working Paper.

   Andersson, H. and Britton, T. (2000). Stochastic Epidemic Models and Their Statistical Analysis. Springer, New York.

   Atkeson, A. (2020). What will be the economic impact of COVID-19 in the US? rough estimates of disease scenarios. National Bureau of Economic Research Working Paper.

   Bacry, E., Delattre, S., Hoffmann, M., and Muzy, J.-F. (2013a). Modelling microstructure noise with mutually exciting point processes. Quantitative Finance, 13(1):65–77.

   Bacry, E., Delattre, S., Hoffmann, M., and Muzy, J.-F. (2013b). Some limit theorems for Hawkes processes and application to financial statistics. Stochastic Processes and their Applications, 123(7):2475–2499.

   Bailey, N. T. J. (1950). A simple stochastic epidemic. Biometrika, 37(3/4):193–202.

   Bailey, N. T. J. (1953). The total size of a general stochastic epidemic. Biometrika, 40(1/2):177–185.

   Bailey, N. T. J. (1957). The Mathematical Theory of Epidemics. Griffin, London.

   Ball, F. (1983). The threshold behaviour of epidemic models. Journal of Applied Probability, 20(2):227–241.

   Ball, F., Britton, T., and Neal, P. (2016). On expected durations of birth–death processes, with applications to branching processes and SIS epidemics. Journal of Applied Probability, 53(1):203–215.

   Ball, F. and Donnelly, P. (1995). Strong approximations for epidemic models. Stochastic Processes and their Applications, 55(1):1–21.

   Bartlett, M. (1949). Some evolutionary stochastic processes. Journal of the Royal Statistical Society. Series B (Methodological), 11(2):211–229.

   Bartlett, M. S. (1956). Deterministic and stochastic models for recurrent epidemics. In Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, volume 4, pages 81–109.

   Berger, D. W., Herkenhoff, K. F., and Mongey, S. (2020). An SEIR infectious disease model with testing and conditional quarantine. National Bureau of Economic Research Working Paper.

   Bowsher, C. G. (2007). Modelling security market events in continuous time: Intensity based, multivariate point process models. Journal of Econometrics, 141(2):876–912.

   Brauer, F., Van den Driessche, P., and Wu, J. (2008). Mathematical Epidemiology. Springer, Berlin.

   Britton, T. (2010). Stochastic epidemic models: a survey. Mathematical Biosciences, 225(1):24–35.

   Daley, D. J. and Gani, J. (1999). Epidemic Modelling: An Introduction. Cambridge University Press, Cambridge.

   Dassios, A. and Zhao, H. (2011). A dynamic contagion process. Advances in Applied Probability, 43(3):814–846.

   Dassios, A. and Zhao, H. (2017a). Efficient simulation of clustering jumps with CIR intensity. Operations Research, 65(6):1494–1515.

   Dassios, A. and Zhao, H. (2017b). A generalised contagion process with an application to credit risk. International Journal of Theoretical and Applied Finance, 20(1):1–33.

   Diekmann, O., Heesterbeek, H., and Britton, T. (2013). Mathematical Tools for Understanding Infectious Disease Dynamics. Princeton University Press, New Jersey.

   Eichenbaum, M. S., Rebelo, S., and Trabandt, M. (2020). The macroeconomics of epidemics. National Bureau of Economic Research Working Paper.

   Embrechts, P., Liniger, T., and Lin, L. (2011). Multivariate Hawkes processes: an application to financial data. Journal of Applied Probability, 48A:367–378.

   Fuchs, C. (2013). Inference for Diffusion Processes: With Applications in Life Sciences. Springer-Verlag, Berlin.

   Guerrieri, V., Lorenzoni, G., Straub, L., and Werning, I. (2020). Macroeconomic implications of COVID-19: Can negative supply shocks cause demand shortages? National Bureau of Economic Research Working Paper.

   Hawkes, A. G. (1971a). Point spectra of some mutually exciting point processes. Journal of the Royal Statistical Society. Series B (Methodological), 33(3):438–443.

   Hawkes, A. G. (1971b). Spectra of some self-exciting and mutually exciting point processes. Biometrika, 58(1):83–90.

   Keeling, M. J. and Rohani, P. (2008). Modeling Infectious Diseases in Humans and Animals. Princeton University Press, New Jersey.

   Kermack, W. O. and McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society, 115(772):700–721.

   Large, J. (2007). Measuring the resiliency of an electronic limit order book. Journal of Financial Markets, 10(1):1–25.

   Lauer, S. A., Grantz, K. H., Bi, Q., Jones, F. K., Zheng, Q., Meredith, H. R., Azman, A. S., Reich, N. G., and Lessler, J. (2020). The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Annals of Internal Medicine, 172(9):577–582.

   Liu, Y., Gayle, A. A., Wilder-Smith, A., and Rocklöv, J. (2020). The reproductive number of COVID-19 is higher compared to SARS coronavirus. Journal of Travel Medicine, 27(2):1–4.

   Martcheva, M. (2015). An Introduction to Mathematical Epidemiology. Springer, New York.

   McKendrick, A. (1925). Applications of mathematics to medical problems. Proceedings of the Edinburgh Mathematical Society, 44:98–130.

   McNeill, W. H. (1976). Plagues and Peoples. Anchor, New York.

   Qu, Y., Dassios, A., and Zhao, H. (2019). Efficient simulation of Lévy-driven point processes. Advances in Applied Probability, 51(4):927–966.

   Tian, H., Liu, Y., Li, Y., Wu, C.-H., Chen, B., Kraemer, M. U. G., Li, B., Cai, J., Xu, B., Yang, Q., Wang, B., Yang, P., Cui, Y., Song, Y., Zheng, P., Wang, Q., Bjornstad, O. N., Yang, R., Grenfell, B. T., Pybus, O. G., and Dye, C. (2020). An investigation of transmission control measures during the first 50 days of the COVID-19 epidemic in China. Science, 368(6491):638–642.