Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
70,903 characters · 15 sections · 46 citation commands
Time series analysis of COVID-19 Infection Curve: a change-point perspective
Since the initial outbreak of the novel coronavirus in Wuhan, China in early January 2020, the COVID-19 pandemic has rapidly spread across the world. Due to the high infectivity of the virus and the lack of immunity in the human population, the epidemic grows exponentially without intervention, and thus can greatly stress the public health system and bring enormous disruption to economy and society. Thus, a crucial task facing every country is to reduce the transmission rate and flatten the (infection) curve. Various emergency measures, such as regional lockdown and mass testing, have been taken by different countries and a natural question is whether (and to what degree) these interventions are effective in slowing down the pandemic. Additionally, each country is at a different stage of the epidemic and it is essential for countries to understand its own pattern of virus growth, as such information is critical for important policy decisions such as extending lockdown or reopening. To (at least partially) answer these questions, a natural step is to analyze the trajectory of the infection curve of COVID-19 since the initial outbreak in each country.
In this paper, we propose to model the time series of cumulative confirmed cases and deaths (in log scale) of each country via a piecewise linear trend model (see formal definition later). In other words, we model the mean of the logarithm of cumulative infection as a linear trend with an unknown number of potential changes in the intercept and slope, as it is natural to expect that the spread of COVID-19 may experience several phases, where the initial growth is typically rapid due to absence of immunity and lack of preparation, and the spread may then evolve into phases with slower growth depending on government intervention and public health responses (i.e. flattening the curve). The estimation of such a model can be formulated as a change-point detection problem.
In recent years, change-point analysis has become an increasingly active research area in statistics and econometrics thanks to its applications across a wide range of fields, including bioinformatics Fan2017, climate science Gromenko2017, economics bai1994least,bai1997estimation,Cho2015, finance fryzlewicz2014wild, medical science chen2011parametric, and signal processing SP2018; see perron2006dealing, aue2013structural and truong2020selective for some recent reviews. However, most existing change-point literature operates under the piecewise stationarity assumption, where it is assumed that the time series of interest is (potentially) non-stationary but can be partitioned into piecewise stationary segments such that observations within each segment are stationary and share a common parameter of interest such as mean or variance. While the piecewise stationarity assumption is proven to be reasonable and fruitful for many applications, methods developed under this framework cannot handle time series with intrinsic non-stationarity, such as the cumulative infection curve of COVID-19.
A simple but important class of time series with intrinsic non-stationarity is the piecewise linear trend model, which has the following mathematical formulation. Let the time series $\{Y_t\}_{t=1}^n$ admit
where $(a_t, b_t)$ is the linear trend (intercept and slope) of $\mathbb{E}(Y_t)$ at time $t$, $\{u_t\}$ is a weakly dependent stationary error process, $\mathbb{\pmb{\tau}}=(\tau_1,\cdots,\tau_m)$ denotes the $m\geq 0$ change-points with the convention that $\tau_0=0$ and $\tau_{m+1}=n$, and we require $\mathbb{\pmb{\beta}}^{(i)}\neq \mathbb{\pmb{\beta}}^{(i+1)}, i=1,\cdots,m$. In this paper, we set $\{Y_t\}_{t=1}^n$ to be the time series of daily cumulative confirmed cases or deaths (in log scale) of COVID-19. Due to the log transformation, the slope $b_t$ naturally measures the growth rate of the virus at day $t$.
The piecewise linear trend model is intuitive, interpretable and is useful for tracking the dynamics of a pandemic as it naturally segments the spread process into phases with (approximately) the same growth rate. The slope of the last segment can shed light on the current status of the pandemic and provide short-term forecast, while the estimated change-points can be compared with dates when emergency measures such as lockdown were introduced to help assess the effectiveness of different policies. Also, the semiparametric nature of (ref) helps to achieve model flexibility while maintaining simplicity, which is advantageous for modeling the cumulative cases at the early stage of a pandemic as the time series is relatively short, curbing the use of sophisticated fully nonparametric methods.
An important part in estimation of (ref) is to recover the unknown number $m$ and location $\mathbb{\pmb{\tau}}$ of the change-points. As discussed above, such a problem has mostly been ignored in the change-point literature with only a few exceptions. A CUSUM based detection algorithm is proposed in baranowski2019narrowest, and a model selection based procedure is derived in Fearnhead2019. However, both methods assume temporal independence of $\{u_t\}$, which can be restrictive as serial dependence is commonly found in time series data. Although baranowski2019narrowest briefly discussed possible extensions to temporally dependent series, potentially important issues such as choice of tuning parameters seem not carefully addressed. bai1998estimating can detect structural breaks in the linear trend model under serial dependence. However, numerical study (see Section 4) suggests that their method is relatively sensitive to positive temporal dependence, which is indeed exhibited by the COVID-19 data, and may give less favorable estimation performance under small sample size.
Based on the self-normalization (SN) idea in shao2010self, we propose a novel SN-based change-point detection procedure for the estimation of (ref) that is robust to temporal dependence both in asymptotic theory and in finite sample. The essential idea of SN is using an inconsistent variance estimator to absorb the unknown serial dependence in the data. See a brief review of SN in Section 2.1 and shao2015self for a comprehensive overview of recent developments of SN for low dimensional time series.
Using the proposed SN method and the piecewise linear trend model, we analyze the time series of cumulative confirmed cases and deaths of COVID-19 (in log scale) in 30 major countries. We find that the spread of coronavirus in each country can typically be segmented into several phases with distinct growth rates and countries with geographical proximity share similar spread patterns, which is particularly evident for continental European countries and developing countries in Latin America. In addition, the transition date from rapid growth phases to moderate growth phases is typically associated with the initiation of emergency measures such as lockdown and mass testing with contact tracing, which partially provides evidence that strict social distancing rules help slow down the virus growth and flatten the curve. Moreover, our analysis further indicates that compared to developed countries, most developing countries are still in the early stages of the pandemic and are generally less efficient in terms of controlling the spread of coronavirus, thus may need more international aids to help contain the epidemic.
Combining the SN-based change-point detection algorithm with a flexible extrapolation function, we further design a simple two-stage forecasting scheme for COVID-19. The proposed method is used to forecast the cumulative deaths in the U.S. and is found to deliver accurate prediction valuable to data-driven public health decision-making.
In this section, we propose a novel SN-based method for change-point detection in model ((ref)) that is robust against a wide range of temporal dependence. Specifically, an SN-based test statistic is first proposed for testing a single change-point alternative and then modified to consistently estimate the change-point. A multiple change-point estimation procedure is further developed by combining the proposed SN test with the NOT algorithm in baranowski2019narrowest.
We start with a change-point testing problem where for model (ref) we want to test the null hypothesis $H_0$ of no change-point against the alternative $H_a$ of one change-point:
where $\mathbb{\pmb{\beta}}_t=(a_t,b_t)$, $\tau=\lfloor\kappa n\rfloor$ is an unknown change-point satisfying $\epsilon<\kappa<1-\epsilon$ for some $0<\epsilon<1/2$ and $\epsilon$ is the commonly used trimming parameter in the change-point analysis (see e.g. andrews1993tests).
Throughout this paper, we operate under the following mild assumption of $\{u_t\}$, which covers a wide range of weakly dependent error process and is weaker than most existing literature where independence of $\{u_t\}$ is assumed.
Assumption (ref)(i) is popular in the linear process literature to ensure the central limit theorem and the invariance principle. Assumption (ref)(ii) is basically equivalent to the geometric moment contracting condition for the nonlinear causal process (wushao2004, wu2005), which implies invariance principle.
Earlier works on this testing problem include andrews1993tests and bai1998estimating where Lagrangian multiplier, Wald, likelihood ratio and $F$ statistics are considered. These tests typically require an estimator of the long-run variance (LRV) $\Gamma$ due to the unknown temporal dependence of the error process $\{u_t\}$. However, as pointed out in shao2010testing, the size and power performance of these tests may depend crucially on the selection of various tuning parameters. In particular, if a data-driven bandwidth parameter is used for the estimation of LRV, an undesirable non-monotonic power phenomenon may occur; see cv2007 and shao2010testing. To avoid the bandwidth selection involved in the estimation of LRV, we instead adapt the idea of self-normalization in shao2010self, which was originally proposed for inference of stationary time series and was generalized to change-point testing for piecewise stationary time series in shao2010testing and zhang2018unsupervised. See shao2015self for a review of SN.
To proceed, we first introduce some notations. Given $\epsilon,$ denote $h=\lfloor \epsilon n\rfloor$. For a vector $x$, denote the $l_2$ norm as $\|x\|_2$ and denote $x^{\otimes2}=xx^{\top}$. Define $F(s)=(1,s)^{\top}$, for $1\leq i<j\leq n$, we denote $\widehat{\mathbb{\pmb{\beta}}}_{i,j}=\Big[\sum_{t=i}^{j}F(t/n)F(t/n)^{\top}\Big]^{-1}\sum_{t=i}^{j}F(t/n)Y_t$ as the OLS estimator of $\mathbb{\pmb{\beta}}$ based on $\{Y_{t}\}_{t=i}^j$. For any $1\leq t_1 <k <t_2\leq n$, given the subsample $\{Y_t\}_{t=t_1}^{t_2}$ and a potential change-point $k$, we define a contrast statistic $D_{n}$ where
Note that $D_{n}(t_1,k,t_2)$ is a normalized difference between the OLS estimates of $\mathbb{\pmb{\beta}}$ with pre-$k$ samples $\{Y_t\}_{t=t_1}^{k}$ and post-$k$ samples $\{Y_t\}_{t=k+1}^{t_2}$. Intuitively, a large $\max_{h\leq k\leq n-h}\|D_{n}(1,k,n)\|_2$ leads to the rejection of $H_0$. However, the asymptotic distribution of $D_{n}(1,k,n)$ depends on the unknown LRV of $\{u_t\}$, and as discussed before the accurate estimation of LRV is rather challenging and problematic in practice.
To bypass the problematic estimation of LRV, we utilize the self-normalization technique. Define $0< \delta<\epsilon/2$ as a local trimming parameter, we define the self-normalizer $V_{n,\delta}(t_1,k,t_2)=L_{n,\delta}(t_1,k,t_2)+ R_{n,\delta}(t_1,k,t_2)$ where
The local trimming parameter $\delta$ is introduced to make sure all the subsample estimates of $\mathbb{\pmb{\beta}}$ in the self-normalizer $V_{n,\delta}(t_1,k,t_2)$ are constructed with a subsample of size being a positive fraction of n, which is a technical condition necessary in our theoretical analysis. We later discuss the implication of the trimming parameters $(\epsilon, \delta)$.
Based on the contrast statistic $D_n(1,k,n)$ and the self-normalizer $V_{n,\delta}(1,k,n)$, we propose an SN-based test statistic $G_n$ for testing the single change-point alternative where
Intuitively, due to the presence of the self-normalizer, the LRVs in $D_{n}(1,k,n)$ and $V_{n,\delta}(1,k,n)$ cancel out with each other, leading to a test statistic $G_n$ that is invariant to LRV. This phenomenon is made formal in Theorem (ref).
Denote $\overset{\mathcal D}{\longrightarrow}$ as convergence in distribution and $\mathbf{b}=\mathbb{\pmb{\beta}}^{(2)}-\mathbb{\pmb{\beta}}^{(1)}$. Define $Q(r)=\int_{0}^{r}F(s)F(s)^{\top}ds$ and $B_F(r)=\int_{0}^{r}F(s)dB(s)$ where $B(\cdot)$ is a standard Brownian motion. Theorem (ref) states the asymptotic behavior of the SN test statistic $G_n$ under $H_0$ and $H_a$ respectively.
Due to self-normalization, the limiting distribution $G(\epsilon, \delta)$ in ((ref)) is pivotal and invariant to the LRV. The corresponding critical values can be easily obtained via simulation. Table (ref) gives the $1-\alpha$ quantiles of $G(\epsilon, \delta)$ for some combinations of $(\epsilon, \delta)$ (based on 10000 replications). Note that the limiting null distribution $G(\epsilon, \delta)$ explicitly depends on the choice of $(\epsilon,\delta)$, thus the impact of trimming parameters $(\epsilon,\delta)$ is accounted for at the first order, in the same spirit of the fixed-$b$ asymptotics (kv2005). See also zhou2013inference. Throughout the paper, we set $(\epsilon,\delta)=(0.1,0.02)$.
Give that the null hypothesis $H_0$ is rejected, we estimate the change-point ${\tau}$ by $\widehat{\tau}= \arg\max_{k\in\{h,\cdots,n-h\}}T_{n,\delta}(k).$ The following theorem gives the consistency result of $\widehat{\kappa}=n^{-1}\widehat{\tau}$.
Theorem (ref) allows a diminishing change size $\|\mathbf{b}\|_2$ with the sample size $n$ as long as $n\|\mathbf{b}\|^2_2\to\infty$. Note that no consistency result is provided in shao2010testing for the change-point location estimation, and our result seems to be the first formal attempt based on the SN technique. However, it is challenging to obtain an explicit rate of convergence for $\widehat{\tau}$ due to the complicated nature of the self-normalizer $V_{n,\delta}$ and we leave it for future investigation.
To extend single change-point testing to multiple change-point estimation, the classical idea is to combine the change-point test with binary segmentation (BS). Although conceptually and computationally simple, it is well known that BS can cause severe power loss for detecting non-monotonic changes Olshen2004, which is common in real data. Several variants of BS have been proposed to address this drawback, such as wild binary segmentation (WBS) (fryzlewicz2014wild) and Narrowest-Over-Threshold (NOT) (baranowski2019narrowest). Since NOT is shown to be superior to WBS, we combine the SN-based test with the NOT algorithm to estimate multiple change-points and name our algorithm SN-NOT.
The essential idea of SN-NOT is to compute the SN test on a large collection of random subsamples of $\{Y_t\}_{t=1}^n$ instead of the entire sample $\{Y_t\}_{t=1}^n$. With high probability, some subsamples will only contain a single change-point, where the SN test statistics are expected to exhibit large values, leading to the discovery of a change-point.
Denote $F_n^{M}=\{(s_i,e_i):i=1,\cdots,M\}$ as the set of $M$ random intervals such that each pair of integers $(s_i,e_i)$ are drawn uniformly from $\{1,\cdots, n\}$ and satisfy $1\leq s_i<e_i\leq n$ and $e_i-s_i+1\geq 2h$. For each random interval $(s,e) \in F_n^M$, we calculate the SN test $$ G_{n,\delta}(s,e)=\max_{k\in\{ s+h-1,\cdots,e-h\}}T_{n,\delta}(s,k,e),\quad T_{n,\delta}(s,k,e)= D_{n}(s,k,e)V_{n,\delta}(s,k,e)^{-1}D_{n}(s,k,e)^{\top}. $$ SN-NOT finds the narrowest interval $(s,e)\in F_n^M$ where the test statistic $G_{n,\delta}(s,e)$ exceeds a given threshold $\zeta_n$ and estimates the change-point as $\widehat{\tau}=\arg\max_{k\in\{ s+h-1,\cdots,e-h\}}T_{n,\delta}(s,k,e)$. Note that for large $M$, with high probability there is only one change-point in this narrowest interval, which thus remedies the drawback of BS in detecting non-monotonic changes. Once a change-point $\widehat{\tau}$ is identified, SN-NOT then divides the sample into two subsamples accordingly and apply the same procedure on each of them. The process is implemented recursively until no change-point is detected. In addition to the advantage of detecting non-monotonic changes, SN-NOT broadens the applicability of the NOT algorithm itself by allowing for temporal dependence in the error process thanks to the self normalization technique.
The detailed implementation of SN-NOT is given in Algorithm (ref). We propose to select the threshold $\zeta_n$ as follows. Generate $B$ sequences of i.i.d $\mathcal{N}(0,1)$ random variables $\{\varepsilon_t^b\}_{t=1}^n$, $b=1,\cdots, B$; for the $b$th sample, we calculate $$ \zeta_n^{b}=\arg\max_{i=1,\cdots,M}G_{n,\delta}(s_i,e_i),\quad b=1,\cdots,B. $$ The threshold $\zeta_n$ is set as the 95% sample quantile of $\{\zeta_n^{b}\}_{b=1}^{B}$. Since the SN test statistic is asymptotically pivotal, this threshold is expected to well approximate the $95\%$ quantile of the finite sample distribution of the maximum SN test statistic on the $M$ random intervals under null. Throughout this paper, we set $B=1000$, $M=300$.
In this section, we study the finite sample performance of the SN test in testing single change-point and the SN-NOT algorithm in detecting multiple change-points through numerical experiments. All results are reported based on 1000 replications.
We generate the data from model ((ref)) with sample size $n=100$, $500$ and $1000$ respectively. For the size performance, we let $\mathbb{\pmb{\beta}}=(3,0.05n)$ while for the power performance, we let $\mathbb{\pmb{\beta}}^{(1)}=(3,0.06n)$ and $\mathbb{\pmb{\beta}}^{(2)}=(3+0.015n,0.03n)$ with the change-point $\tau=n/2$. The error process $\{u_t\}$ is generated via an AR(1) model where $u_t=\rho u_{t-1}+e_t$, $e_t\overset{i.i.d.}{\sim}$ $\mathcal{N}(0,(1-\rho^2)\sigma^2)$ with $\rho=0,\pm0.2,\pm 0.5$ and $\sigma=0.15$.
For comparison, we also implement the supLM test defined in andrews1993tests (using function {\tt sctest} of the R package strucchange) with the same trimming parameter $\epsilon=0.1$. The results are summarized in Table (ref) at significance levels $\alpha=5\%$ and 10%. It can be seen that when $n$ is small, both methods have distorted sizes. In particular, SN is prone to be conservative when $\rho$ is negative and oversized when $\rho$ is positive while supLM is undersized in all cases. As $n$ increases, we find that both tests tend to have more accurate sizes. For $n=100$, supLM test has slightly higher power than SN test while for $n=500$ and $n=1000$, SN test beats supLM test under positive $\rho$. Note that both tests are more powerful under negative $\rho$.
We examine the numerical performance of SN-NOT by considering the following DGP with $n=100$:
The error process $\{u_t\}$ is generated via an AR(1) model where $u_t=\rho u_{t-1}+e_t$, $e_t\overset{i.i.d.}{\sim}\mathcal{N}(0,(1-\rho^2)\sigma^2)$ with $\rho=0,\pm0.2,\pm 0.5$ and $\sigma=0.15$. For comparison, we also implement the multiple change-point detection procedure proposed in bai1998estimating (denoted as BP hereafter), which is the most widely used detection algorithm allowing for temporal dependence in the error term of model ((ref)). BP is implemented using function {\tt breakpoints} of the R package strucchange.
To assess the accuracy of change-point estimation, we define the Hausdorff distance between two sets. Denote the set of true change-points as $\mathbb{\pmb{\tau}}_o$ and the set of estimated change-points as $\widehat{\mathbb{\pmb{\tau}}}$, we define $d_1(\mathbb{\pmb{\tau}}_o, \widehat{\mathbb{\pmb{\tau}}})=\max_{\tau_1 \in \widehat{\mathbb{\pmb{\tau}}}}\min_{\tau_2 \in \mathbb{\pmb{\tau}}_o}|\tau_1-\tau_2|$ and $d_2(\mathbb{\pmb{\tau}}_o, \hat{\mathbb{\pmb{\tau}}})=\max_{\tau_1 \in \mathbb{\pmb{\tau}}_o}\min_{\tau_2 \in \widehat{\mathbb{\pmb{\tau}}}}|\tau_1-\tau_2|$, where $d_1$ measures the over-segmentation error of $\widehat{\mathbb{\pmb{\tau}}}$ and $d_2$ measures the under-segmentation error of $\widehat{\mathbb{\pmb{\tau}}}$. The Hausdorff distance is then defined as $d_H(\mathbb{\pmb{\tau}}_o, \widehat{\mathbb{\pmb{\tau}}})=\max(d_1(\mathbb{\pmb{\tau}}_o, \widehat{\mathbb{\pmb{\tau}}}), d_2(\mathbb{\pmb{\tau}}_o, \widehat{\mathbb{\pmb{\tau}}})).$ In addition, we report the adjusted Rand index (ARI) which measures the similarity between two partitions of the same observations. Roughly speaking, a higher ARI (with the maximum value of 1) means more accurate change-point estimation. For the definition and detailed discussions of ARI, we refer to hubert1985comparing.
Table (ref) summarizes the numerical result where we report ARI, $d_1$, $d_2$, $d_H$ and the frequency of $|\widehat{m}-m_o|$ for SN-NOT and BP. It can be seen that SN-NOT is overall better than BP in terms of ARI, $d_H$ and the estimated number of change-points when $\rho\geq 0$. This finding suggests using SN-NOT could be more advantageous for analyzing COVID-19 data, which exhibit positive temporal dependence (see the last column of Table (ref)). For applications where negatively correlated error is expected, BP could be a better choice.
In this section, based on the proposed SN-NOT algorithm, we provide detailed in-sample analysis of the cumulative confirmed cases (Section 4.2-4.3) and deaths (Section 4.4) of COVID-19 (in log scale) in 30 major countries.
We focus on G20 (with 19 sovereign countries{\footnote{G20 is an international forum for the governments and central bank governors from 19 countries and the European Union. We will view members of the European Union as individual countries because the responses to COVID-19 usually come from the national level.}}) and 11 other countries leading the total infected cases as of May 27, 2020, including Australia (AUS), Argentina (ARG), Belgium (BEL), Brazil (BRA), Canada (CAN), Chile (CHI), China (CHN), France (FRA), Germany (GER), India (IND), Indonesia(INA), Iran (IRI), Italy (ITA), Japan (JPN), Mexico (MEX), Netherlands (NED), Pakistan (PAK), Peru (PER), Portugal (POR), Qatar(QAT), Russia (RUS), Saudi Arabia (KSA), Spain (ESP), South Africa (RSA), South Korea (ROK), Sweden (SWE), Switzerland (SUI), Turkey (TUR), United Kingdom (GBR), United States (USA).
We obtain the data from \url{https://ourworldindata.org/coronavirus-source-data} maintained by “Our World in Data", where cumulative measures such as confirmed cases and deaths are updated daily for each nation. For each country, the logarithm of cumulative confirmed cases (or deaths) $\{Y_t\}$ starts on the date when the cumulative cases (or deaths) exceeded 20 and ends on May 27.
We study the cumulative confirmed cases and deaths (in log scale) of each country via the piecewise linear trend model (ref), where given $\{Y_t\}$, the change-points $(\tau_1,\cdots, \tau_{\widehat{m}})$ is estimated by the SN-NOT algorithm. An OLS is then used to recover the linear model for the $i$th estimated segment $\{Y_t\}_{t=\widehat{\tau}_{i-1}+1}^{\widehat{\tau}_i}$, $i=1,2,\cdots,\widehat{m}+1$. With a slight abuse of notation, denote $\widehat{b}_{i}$ as the estimated slope for the $i$th segment. We define the normalized slope $S_i=\widehat{b}_{i}/n$ for each segment. As can be seen from ((ref)), the normalized slope $S_i$ measures $\mathbb{E}[Y_{t+1}-Y_t]$ for the $i$th segment, which can be interpreted as the “log-return" and measures the daily growth rate of the cumulative confirmed cases (or deaths) in the original scale.
Methodologically speaking, for cumulative confirmed cases, the piecewise linearity allows us to assess the growth rate of the coronavirus at any given time and further facilitates short-term forecast. In particular, the estimated slope $S_i$ of each segment indicates the pace of the growth rate during the corresponding period. Moreover, by comparing the slope before and after each change-point, we can quantitatively assess the changes in growth rate, which partially measure the effectiveness of policies taken by the government.
We first conduct a detailed case study for eight representative countries that either lead confirmed cases (the U.S., Brazil, Russia, and India) in the corresponding continent or receive most media attention (the U.K., Spain, Italy, and South Korea).
Table (ref) summarizes the detailed estimation result for each country (in descending order of the cumulative confirmed cases), where we report the starting date of the series, length of the series $n$, the estimated number of change-points, dates of the first, second and latest estimated change-point. The first ($S_1$), the second ($S_2$) and the current normalized slope ($S_{\widehat{m}+1}$) are also presented. In addition, we report the lag-1 sample autocorrelation $\widehat{\rho}$ of the error process. From the table, we can see all of these countries have been affected by the coronavirus for more than two months. The average length of segments between two adjacent change-points is around 13-20 days, indicating that the spread rate can be relatively steady for a window of 2-3 weeks. The latest change-point for most countries appeared in May except for Brazil. We also note that the current normalized slopes (i.e. growth rate) vary considerably across countries with comparably large values in Brazil and India. Meanwhile, the lag-1 sample autocorrelation $\widehat{\rho}$ are all positive, which suggests the use of SN-NOT instead of BP as discussed in Section 3.2. In Figures (ref) and (ref) of the supplementary material, we further plot the lag-1 to lag-30 ACF and PACF of the residuals, which rules out the scenario of long memory and supports the validity of Assumption 2.1.
Figure (ref) visualizes the estimated piecewise linear models for the eight countries, which gives a more direct perception of how the growth rate changes over time. Note that the U.S. and South Korea are the only two countries that witnessed an increase in the slope after the first change-point. For the U.S., the first change-point is March 4, one day after the first confirmed case appeared in New York. Since then, the pandemic underwent an outbreak in the New York state, which has been the leading state in the U.S. in terms of infected cases. The second change-point appeared on March 24, after which the slope began to drop. This is also noteworthy as on March 20, the U.S. began barring entry of foreign nationals who had traveled to 28 European countries within the past 14 days. While in South Korea, after February 18, the infected cases increased drastically, and the slope dropped after March 3. We find that the first change-point is the day when the first super-spreader in South Korea was diagnosed {\footnote{A member of the Shincheonji religious organization was diagnosed as 31st case in Daegu, see \url{https://foreignpolicy.com/2020/02/27/coronavirus-south-korea-cults-conservatives-china/}}}. The second change-point, March 3, is when the drive-through testing was made widely available to Korean citizens.
The growth rate decreased after the first change-point in other countries. For the U.K., the first and second change-points are quite close. In particular, we find the U.K. governments gradually increased the restrictions on freedom of movement for the general public between these two change-points (March 20 and March 29). This could help explain why both change-points are associated with significant drops in the virus growth rate. In addition, we find that Italy extended the quarantine lockdown from region-focused to nationwide on March 10, one day after the first estimated change-point. For Spain, the first change-point is estimated as March 14, which is one day after Spain declared the nationwide state of emergency. Similar to Italy, the slopes dropped drastically after the first change-point. Generally speaking, the first or second change-point of these countries are closely associated with the date when local or nationwide interventions from the governments were initiated. These countries typically transition from a rapid growth phase to a moderate growth phase after the first or second change-point. This may serve as evidence that government intervention such as lockdown and massive testing could effectively slow down the spread of the coronavirus.
From Figure (ref), we also find the situations in Brazil, Russia and India rather somber, as of May 27. Russia is still transitioning from the rapid growth phase to the moderate growth phase, while the fast growing trend in Brazil has not changed since April 12. Even though Brazil managed to bring down the slope by a significant amount at the first change-point on March 25, it seemed the right-wing government took few follow-up effective measures. The situation in India is also grim where the decreases of growth rate at the first and second change-points are quite small and the current growth rate is still high, suggesting that stricter measures to be taken. In summary, these three countries still have a long way to go in terms of slowing down the spread of COVID-19.
We further extend the scope of analysis to 30 countries to obtain a relatively complete picture of the pandemic situations around the world. Specifically, we conduct a comparative study based on two important quantities: the maximum normalized slope and the current normalized slope, which are estimated by $S_{max}={\max_{1\leq i\leq \widehat{m}+1}n^{-1}\widehat{b}_i}$ and $S_{cur}=n^{-1}\widehat{b}_{\widehat{m}+1}$ respectively. Combined together, the two measures allow us to obtain an overall picture of the phase when the virus transmitted fastest and the current situation in each country. In particular, $S_{max}$ provides information on the growth rate at the early stage of the pandemic for a particular country. In this phase, often no government regulations are imposed so it depicts the worst scenario if no emergency measure is taken. $S_{cur}$ gives the ongoing epidemic growth rate and could help make predictions in the short run.
In Figure (ref), we plot $S_{max}$ against $S_{cur}$ for each country. Note that by their relative positions in Figure (ref), the 30 countries can be roughly grouped into three clusters: East Asian countries and Australia, European and North American countries and Other developing countries. We find that countries within the same cluster tend to have similar current growth rate. China, South Korea, and Australia are among the best with $S_{cur}$ close to zero. Most European and North American countries are in the second tier while countries in continental Europe generally have slower ongoing virus growth than the U.K., the U.S. and Canada. The only exceptions are Sweden and Russia. In fact, Sweden adopted a different strategy than other countries in that no lockdown has been imposed by the government and large parts of its society remain open. Note that Figure (ref) does not take the time effect into account, thus the cluster along the horizontal direction may also be attributed to the cluster of similar eruption time of the virus. This could help explain why Russia is closer to developing countries and why Latin American countries have the largest $S_{cur}$.
To take the time factor into consideration, in Figure (ref), we plot the ratio $S_{cur}/S_{max}$ against the days in between (i.e. $\tau_{cur}-\tau_{max}$ with $\tau_{max}$ as the start date for the segment with the largest slope and $\tau_{cur}=\tau_{\widehat{m}}$ as the latest change-point), which allows us to further understand how the growth rate changes from its peak to the current status with time. Horizontally speaking, for the same ratio $S_{cur}/S_{max}$, if country A is to the left of country B, then A acts faster than B in bringing down the virus growth from its peak value. Vertically speaking, for the same time length $\tau_{cur}-\tau_{max}$, if A is below B, then A is more effective than B in reducing the growth rate.
We again find that most European and North American countries tend to share similar characteristics. The growth rates in the current phases for these countries are less than one-tenth of their peak value, and it took them about two to three months to achieve that. From the lower panel in Figure (ref), we find that South Korea, China and Australia outperform other countries as the ratios were brought to near zero in around 65 days. Again, we find that continental European countries (except Russia and Sweden) perform better than U.S, Canada and U.K.
Most developing countries are on the top-left of the plot, suggesting that they are still in the relatively early stage of the pandemic and the situation has not improved much since the beginning of the outbreak. In addition, we find Latin American countries, such as Mexico, Brazil, Chile, and Peru, tend to cluster. Given their geographical proximity, this is not a surprise. We note that developing countries tend to be less efficient in slowing the spread of COVID-19. For example, with roughly the same amount of time, the ratios in India and Argentina are three times larger than developed countries. In summary, more caution and attention should be given to the epidemic in developing countries as they may need more international aids compared to the developed countries.
Based on the same methodology, we analyze cumulative deaths in the 30 countries. Note that unlike confirmed cases, public health interventions naturally have a longer lagged effect on coronavirus-related deaths, as severe symptoms may not develop immediately upon infection. Thus, we believe a change-point analysis on cumulative confirmed cases should be preferred in terms of quantifying the effectiveness of emergency policies. Additionally, the criteria for certifying deaths due to COVID-19 vary from nation to nation, thus comparative analysis across countries should be interpreted with caution.
Table (ref) summarizes the detailed estimation result for cumulative deaths in the eight representative countries. Notably, for each country, the estimated number of change-points for deaths is smaller than or equal to that for cumulative confirmed cases in Table (ref). This is intuitive as the history of cumulative deaths is shorter and number of deaths largely depend on infections (with a lag). Note that the duration between the starting date and the first change-point for cumulative deaths is around 2-3 weeks, which is consistent with that for confirmed cases in Table (ref). The same phenomenon also applies to the duration between the first and second change-points. This consistency in part confirms the validity of the change-point estimation results and indicates a 2-3 weeks response lag between changes in growth rate of infections and changes in growth rate of deaths. We note that Italy and Spain have the highest growth rate of cumulative deaths before the first change-point, which highlights the extreme importance of “flattening the curve", as it is known that the exponential surge of coronavirus cases exhausted the public health system in the two countries at the early stage of the pandemic.
Figure (ref) further plots the estimated piecewise linear models for cumulative deaths in the eight countries. The pattern exhibited by each country is largely consistent with its pattern in Figure (ref), except for South Korea. Note that the start date of the cumulative death curve in South Korea is almost 30 days later than the start date of the cumulative confirmed cases, which partially explains the different pattern around its first change-point.
We further conduct a comparative analysis for cumulative deaths in 30 countries. We exclude China, Spain and Qatar in the analysis as the death tolls were either revised or unavailable{\footnote{China revised its death toll upwards on April 17, see {\url{https://www.nytimes.com/2020/04/17/world/asia/china-wuhan-coronavirus-death-toll.html}}. The death toll is not available for Qatar.}}. Figure (ref) plots $S_{max}$ against $S_{cur}$ for each country. Similar to the results for confirmed cases in Figure (ref), European and North American countries tend to cluster while developing countries generally have higher ongoing growth rates $S_{cur}$.
Note that South Korea and Australia deliver the best responses with small $S_{max}$ and near-zero $S_{cur}$ for cumulative deaths. However, it is unexpected to see that western developed countries, such as Italy and the U.K., experience the largest maximum growth rate. Since the maximum growth rate always takes place in the first segment of the cumulative death curve, it indicates that the coronavirus may take these countries by surprise and the health systems may not be well prepared for the flood of coronavirus patients in the early stage of the pandemic. Another notable pattern is that Latin American countries tend to have larger values in both maximum and current growth rates than other developing countries, signaling the possibility of Latin America becoming the next epicenter of the COVID-19 pandemic.
Figure (ref) plots $S_{cur}/S_{max}$ against $\tau_{cur}-\tau_{\max}$ for cumulative deaths in each country, where the observed patterns are similar to the ones for cumulative confirmed cases in Figure (ref). Specifically, developing countries again tend to be less efficient in slowing the spread of COVID-19, where with roughly the same amount of time, the ratios $S_{cur}/S_{max}$ in developing countries are noticeably larger than developed countries.
As stated by the Centers for Disease Control and Prevention (CDC)\footnote{\url{https://www.cdc.gov/coronavirus/2019-ncov/covid-data/forecasting-us.html#why-forecasting-critical}}, accurate forecast of COVID-19 deaths is critical for public health decision-making, as it projects the likely impact of coronavirus to health systems in coming weeks and helps government officials develop data-driven public health policies for controlling the pandemic.
In Section 5.1, we propose a simple and intuitive forecasting scheme for cumulative deaths due to COVID-19 by combining SN-NOT with a flexible extrapolation function. In Section 5.2, we further demonstrate its promising performance in predicting cumulative deaths in the U.S.
As suggested by the analysis in Section 4, the spread of coronavirus typically experiences several different stages due to external interventions. While a sophisticated epidemiology model based on differential equations may manage to take into account information about interventions and characterize the entire cumulative death curve, a more natural (and simpler) solution from the change-point aspect is to first segment the time series into periods with relatively stable behavior and then generate forecast based on observations in the last segment, see for example, pesaran2002market and Bauwens2015.
Following this idea, we propose an SN-NOT based two-stage approach for cumulative deaths prediction. Specifically, in the first stage, given the cumulative deaths (in log scale) $\{Y_t\}_{t=1}^n$, a piecewise linear trend model is estimated via SN-NOT with change-points $\widehat{\mathbb{\pmb{\tau}}}$. In the second stage, a flexible function $f(t)$ is fitted on the last segment $\{Y_t\}_{t=\widehat{\tau}_{\widehat{m}}+1}^n$ with the assumption that $\mathbb{E}(Y_t)=f(t)$ and the $k$-day ahead forecast for cumulative deaths can be readily made via extrapolation of $\widehat{f}(t)$.
Note that the purpose of the first stage (in-sample) change-point analysis is to identify the most recent segment where $\{Y_t\}_{t=1}^n$ exhibits relatively stable behavior and thus facilitates the second stage (out-of-sample) forecast. As demonstrated in Section 4, the piecewise linear trend model with SN-NOT is sufficient for this task. However, as for prediction in the second stage, any flexible extrapolation function $f(t)$ can be considered, as it is expected that a linear function may only provide a reasonable forecast for short horizons due to its limited flexibility.
In the following, we consider three commonly used extrapolation functions (in the order of increasing flexibility) in the literature, including the linear function $f(t)=a+b(t/n)$, the quadratic function $f(t)=c+ d(t/n)+e(t/n)^2$ and the logistic function $f(t)=\dfrac{L}{1+\exp\big(-\alpha(t/n-t_0)\big)}$.
Based on $\{Y_t\}_{t=\widehat{\tau}_{\widehat{m}}+1}^n$, a standard OLS can be used to estimate the linear and quadratic functions and a standard nonlinear least square can be used to estimate the logistic function. The $k$-day ahead forecast for $Y_{n+k}$ is formulated respectively as
The prediction for cumulative deaths on day $n+k$ is $\widehat{\mathrm{Death}}_{n+k}=\exp(\widehat{Y}_{n+k})$.
We apply the SN-NOT based prediction method to forecast cumulative deaths in the U.S. and compare its performance with other forecasting models listed on the CDC website\footnote{\url{https://www.cdc.gov/coronavirus/2019-ncov/covid-data/forecasting-us.html}}. Specifically, following the CDC website, the forecast is generated on five dates, April-27, May-04, May-11, May-18 and May-25, and the forecast horizon is 5-day (one-week) ahead and 12-day (two-week) ahead.
We compare with five forecasting models{\footnote{Other models can be found on the CDC website. The five models are chosen as their predictions are available on all the aforementioned dates while other models only report on some of the recent dates.}} available on the CDC website: “LANL" by LANL, “Imperial" by verity2020estimates, “UT" by UT, “YYG" by YYG and “MOBS" by MOBS. These forecasting methods are mainly ensembles of complex mechanistic models (such as SEIR and SEIS), known as compartmental models in epidemiology, which track the spread of infectious disease via a system of differential equations. To highlight the importance of the first-stage change-point analysis, we additionally report the forecast given by fitting a logistic function on the entire time series without segmentation (and name it “Logistic").
Table (ref) reports the prediction results and the findings can be summarized as follows.
(1) SNL gives comparable performance to other methods for the 5-day ahead forecast, while it considerably overestimates deaths at the 12-day horizon. In other words, linear extrapolation can only be used for short-term forecasts. This is not surprising as the linear function essentially assumes a constant growth rate for the cumulative deaths. While such an approximation is reasonable for short-term, it may not be able to track the growth rate for a long period to make accurate predictions. SNQ generally performs better than SNL due to its increased flexibility, though it tends to underestimate at the 12-day horizon as the quadratic function may pass its peak for long-horizon extrapolation.
(2) SNLG is consistently a top performer among all models thanks to the flexibility of the logistic function, which ensures the fitted curve is non-decreasing and is capable of tracking both increasing and decreasing growth rate. Note that there is a drastic performance difference between the two-stage SNLG forecast and the pure Logistic forecast, which indicates the value of the first-stage change-point estimation for identifying the most recent segment where cumulative deaths exhibit relatively stable behavior.
In summary, the SN-NOT based two-stage prediction, in particular SNLG, provides decent forecasts for the cumulative deaths in the U.S. Considering that SNLG is solely based on the time series of cumulative deaths, this result is rather promising and further confirms the value and validity of the change-point analysis. Though by no means SNLG can replace the complex mechanistic models built on epidemiology principles, we believe it can serve as a meaningful addition to the existing set of forecasting models for tracking the COVID-19 pandemic.