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-exposedinfectious-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 (1949 Bartlett ( , 1956 and Bailey (1950 Bailey ( , 1953 Bailey ( , 1957 . 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 Zhao, 2011) , which is an extension of the classical Hawkes process (Hawkes, 1971a,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.
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 N t with N 0 = 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 N t with two phases: Phase 1 (Full Contagion): For the first phase period t ∈ [0, ), N t follows a dynamic contagion process with stochastic intensity
• δ > 0 is the constant rate of exponential decay; 
Phase 2 (Self Contagion): For the second phase period t ∈ ( , ∞), N t is a pure self-exciting Hawkes process with stochastic intensity The point process N t 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 N t ≡ {T i } i=1,... with N 0 = 0 and stochastic intensity 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.
• The cumulated number of infected cases, N t , is a cluster point process, which consists of two types of points: outside-imported cases and inside-infected cases.
• The arrivals of outside-imported cases follows a Cox process with shot-noise intensity
where externally-exciting jumps arrive as a Poisson process N * t at time points {T * k } k=1,... with sizes (marks) {Z k } k=1,... , and they disappear after time point when the interventions took effect, and there will be no any increase of imported cases in a long run.
• Each imported case may infect other individuals inside and thereby causes new cases, and each of these new cases would further infect others inside, and so on. The infection of any new cases caused by the previous infected cases follows a Cox process with exponentially decaying intensity Y · e −δ(t−T·) , where Y · follows a two-phase distribution G(y; t) and T · is the infection time of the previous infected case. After the interventions took effect, the COVID-19 becomes less easy to spread on average. This is captured by our assumption of two-phase distribution (2.4) for Y · here.
• Overall, the superposition of all these infected cases form a point process N t , a two-phase dynamic contagion process with stochastic intensity (2.3).
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 N t , 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 , N t )
Let {F t } t≥0 be the natural filtration of the point process N t , i.e. F t = σ (N s , s ≤ t) and assume that the intensity process {λ t } t≥0 being F t -adapted. The joint Laplace transform and probability generating function for (λ t , N t ) 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 N t is of the form
θ Ns e −A(s)λs , < s ≤ t, where the boundary condition is
and c(t) is determined by
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 N t is provided in Theorem 3.1 of Dassios and Zhao (2011) . For < s ≤ t, given F s , the infinitesimal generator of the dynamic contagion process (λ s , N s , s) acting on a function f (λ, n, s) within its domain Ω(A) is given by
Consider a function f (λ, n, s) of form f (λ, n, s) = θ n e −A(s)λ , and substitute this into Af = 0, we then have the ODE
adding the boundary condition A(t) = v, gives the ODE in (3.2).
The moments of λ t and N t can be obtained by differentiating the joint Laplace transform and probability generating function of λ t and N t , 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 F s for s ≤ t is given by
and the conditional expectation of the point process N t given F s is of the form
Proof. The result in (3.5) immediately follows Theorem 3.6 in Dassios and Zhao (2011) . And
which directly implies the result in (3.6).
Proposition 3.2. The conditional second moment of the process λ t given F s for s ≤ t is given by
and the conditional joint expectation of the process λ t and the point process N t given F s for s > is of the form
and the second moment of point process N t given F s for s > is of the form
Proof. These results immediately follows Lemma 3.1, Lemma 3.2 and Theorem 3.9 in Dassios and Zhao (2011) .
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 ‹ T to be elimination time such that
(3.10)
The condition probability of the elimination time is provided in Proposition 3.3.
Proof. Given ‹ T being the timing of the last ever event, the event { ‹ T ≤ t} implies that N u = N t for any u ≥ t, which also lead to
Hence, we have
And according to Theorem 3.1, by setting θ = 1 in (3.1), the result follows immediately.
Given the last ever event { ‹ T < t}, one could obtain the expected size of the epidemic at time t.
The relevant details are presented in Corollary 3.1.
Proof. According to Theorem 3.1 and Proposition 3.3, we have
, the result immediately follows (3.13).
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 N t 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 N ∞ conditional on F s is given
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
This implies that, the final epidemic size N ∞ conditional on F follows a mixed-Poisson distribution with the probability mass function
where m(v) is the density function of the mixing distribution,
which is an inverse Gaussian distribution with parameters β δβ−1 λ s and β 2δ λ 2 s .
We provide a calibration scheme based on the daily increments of the two-phase dynamic contagion process N t . Let us first denote the observations of the daily confirmed COVID-19 cases as {C t } t=0,1,2,...,T . The mean square error (MSE) between the expected daily increments of N t and the actual reported daily confirmed COVID-19 cases is given as
(4.1)
We consider the calibration based on minimising the MSE (4.1), i.e., we choose parameters (α,β,δ,ˆ ,ˆ )
such that
with α, β, δ, ≥ 0 and ∈ N + .
Without loss of generality, for simplicity, we assume Z k = 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 Z k , 
The The branching ratio (BR), which demonstrates the average infection rate, is determined by Table 2 , we compare the estimated branching ratios before and after the government interventions came into effect, namely R b and R a , 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 R b and R a . The comparison of R b and R a 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, 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. Table 1 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.
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.
In Table 3 , we also presents the date of day 0, i.e. Date 0 , and the date of government interventions imposed, i.e. Date G . The delay period of government interventions came into effect therefore can be obtained given the estimatedˆ , with Date 0 and Date G . 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 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 R b = R a . 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. 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, R b , and the vertical axis represents the branching ratio after the government interventions came into effect, Ra. 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.
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 deviations 1 for some typical countries and territories. The black dash line in each graph of Figure Table 3 and actual daily confirmed COVID-19 cases for France, Germany, Switzerland, New York.
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 P Ä ‹ 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. 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.
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ˆ ,ˆ , R a against 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 Table 3 . 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 R a 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.
Beside the conditional probability for the elimination time ‹ T , the epidemic size N t given
{T ≤ t} can also be predicted according to the join expectation of N t and { ‹ T ≤ t} derived in Corollary 3.1. In Table 5 , we report the 95% confidence interval for elimination time ‹ for regions and countries with calibration parameters in 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 T , the conditional expectation of the elimination time T , the expected elimination date and the conditional expectation of the epidemic size Nt given T ≤ t with t = E[ T |Fˆ ] under the calibrated parameters (α,β,δ,ˆ ,ˆ ) for these regions and countries suggested in Table 3 .
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 selfexciting 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., P (N ∞ − N = k | F ), for some regions and countries under the calibrated parameters (α,β,δ,ˆ ,ˆ ) for these regions and countries suggested in Table 3 .
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 Table 3. 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.