key: cord-0957037-5rczft4r
authors: Jiang, Feiyu; Zhao, Zifeng; Shao, Xiaofeng
title: Time series analysis of COVID-19 infection curve: A change-point perspective
date: 2020-07-30
journal: J Econom
DOI: 10.1016/j.jeconom.2020.07.039
sha: a09daff9a7fd822b43710e1f076f86c382178389
doc_id: 957037
cord_uid: 5rczft4r

In this paper, we model the trajectory of the cumulative confirmed cases and deaths of COVID-19 (in log scale) via a piecewise linear trend model. The model naturally captures the phase transitions of the epidemic growth rate via change-points and further enjoys great interpretability due to its semiparametric nature. On the methodological front, we advance the nascent self-normalization (SN) technique (Shao, 2010) to testing and estimation of a single change-point in the linear trend of a nonstationary time series. We further combine the SN-based change-point test with the NOT algorithm (Baranowski et al., 2019) to achieve multiple change-point estimation. Using the proposed method, we analyze the trajectory of the cumulative COVID-19 cases and deaths for 30 major countries and discover interesting patterns with potentially relevant implications for effectiveness of the pandemic responses by different countries. Furthermore, based on the change-point detection algorithm and a flexible extrapolation function, we design a simple two-stage forecasting scheme for COVID-19 and demonstrate its promising performance in predicting cumulative deaths in the U.S.

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 (Fan and Mackey, 2017) , climate science (Gromenko et al., 2017) , economics (Bai, 1994 (Bai, , 1997 Cho and Fryzlewicz, 2015) , finance (Fryzlewicz, 2014) , medical science (Chen and Gupta, 2011) , and signal processing (Chen and Gu, 2018) ; see Perron (2006) , Aue and Horváth (2013) and Truong et al. (2020) for some recent reviews. However, 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

where (a t , b t ) is the linear trend (intercept and slope) of E(Y t ) at time t, {u t } is a weakly dependent stationary error process, τ τ τ = (τ 1 , · · · , τ m ) denotes the m ≥ 0 change-points with the convention that τ 0 = 0 and τ m+1 = n, and we require β β β (i) = β β β (i+1) , i = 1, · · · , m.

In this paper, we set {Y t } n t=1 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 changepoints 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 (1.1) 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.

J o u r n a l P r e -p r o o f Journal Pre-proof An important part in estimation of (1.1) is to recover the unknown number m and location τ τ τ 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 Baranowski et al. (2019) , and a model selection based procedure is derived in Maidstone and Letchford (2019) . However, both methods assume temporal independence of {u t }, which can be restrictive as serial dependence is commonly found in time series data. Although Baranowski et al. (2019) briefly discussed possible extensions to temporally dependent series, potentially important issues such as choice of tuning parameters seem not carefully addressed. Bai and Perron (1998) 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 Shao (2010) , we propose a novel SN-based change-point detection procedure for the estimation of (1.1) 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 Shao (2015) 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, J o u r n a l P r e -p r o o f Journal Pre-proof 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 changepoint detection in model (1.1) 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 changepoint estimation procedure is further developed by combining the proposed SN test with the NOT algorithm in Baranowski et al. (2019) .

2.1. Testing for a single change-point. We start with a change-point testing problem where for model (1.1) we want to test the null hypothesis H 0 of no change-point against the alternative H a of one change-point:

where β β β t = (a t , b t ), τ = κn is an unknown change-point satisfying < κ < 1 − for some 0 < < 1/2 and is the commonly used trimming parameter in the change-point analysis (see e.g. Andrews (1993) ).

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.

J o u r n a l P r e -p r o o f

Assumption 2.1. The error process {u t } is strictly stationary such that E(u t ) = 0, E(u 4 t ) < ∞ and the long-run variance satisfies Γ 2 = lim n→∞ Var(n −1/2 n t=1 u t ) ∈ (0, ∞). Denote {e t } as a sequence of i.i.d. random variables with zero mean and unit variance, we further assume that {u t } admits one of the following two representations:

for some measurable function G and F t = (e t , e t−1 , · · · ). For some

is an i.i.d. copy of e 0 and X 4 = (E(X 4 )) 1/4 for a random variable X.

Assumption 2.1(i) is popular in the linear process literature to ensure the central limit theorem and the invariance principle. Assumption 2.1(ii) is basically equivalent to the geometric moment contracting condition for the nonlinear causal process (Wu and Shao (2004) , Wu (2005) ), which implies invariance principle.

Earlier works on this testing problem include Andrews (1993) and Bai and Perron (1998) where Lagrangian multiplier, Wald, likelihood ratio and F statistics are considered. These tests typically require an estimator of the long-run variance (LRV) Γ due to the unknown temporal dependence of the error process {u t }. However, as pointed out in Shao and Zhang (2010) , 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 Crainiceanu and Vogelsang (2007) and Shao and Zhang (2010) . To avoid the bandwidth selection involved in the estimation of LRV, we instead adapt the idea of self-normalization in Shao (2010) , which was originally proposed for inference of stationary time series and was generalized to change-point testing for piecewise stationary time series in Shao and Zhang (2010) and Zhang and Lavitas (2018) . See Shao (2015) for a review of SN.

To proceed, we first introduce some notations. Given , denote h = n . For a vector x, denote the l 2 norm as x 2 and denote x ⊗2 = xx . Define F (s) = (1, s) , for 1 ≤ i < j ≤ n, 

For any 1 ≤ t 1 < k < t 2 ≤ n, given the subsample {Y t } t 2 t=t 1 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 β β β with pre-k samples {Y t } k t=t 1 and post-k samples {Y t } t 2 t=k+1 . Intuitively, a large max h≤k≤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 < δ < /2 as a local trimming parameter, we define the self-normalizer

The local trimming parameter δ is introduced to make sure all the subsample estimates of β β β in the self-normalizer V n,δ (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 ( , δ).

Based on the contrast statistic D n (1, k, n) and the self-normalizer V n,δ (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,δ (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 2.1.

Journal Pre-proof

is a standard Brownian motion. Theorem 2.1 states the asymptotic behavior of the SN test statistic G n under H 0 and H a respectively.

Theorem 2.1. Suppose Assumption 2.1 holds. Let G n be defined in (2.4), we have

Due to self-normalization, the limiting distribution G( , δ) in (2.5) is pivotal and invariant to the LRV. The corresponding critical values can be easily obtained via simulation. Table 2 .1 gives the 1 − α quantiles of G( , δ) for some combinations of ( , δ) (based on 10000 replications). Note that the limiting null distribution G( , δ) explicitly depends on the choice of ( , δ), thus the impact of trimming parameters ( , δ) is accounted for at the first order, in the same spirit of the fixed-b asymptotics (Kiefer and Vogelsang (2005) ). See also Zhou and Shao (2013) . Throughout the paper, we set ( , δ) = (0.1, 0.02).

J o u r n a l P r e -p r o o f Give that the null hypothesis H 0 is rejected, we estimate the change-point τ by τ = arg max k∈{h,··· ,n−h} T n,δ (k). The following theorem gives the consistency result of κ = n −1 τ .

Theorem 2.2. Under H a , suppose Assumption 2.1 holds, and n b 2 2 → ∞ as n → ∞. Then, we have that for any η > 0,

Theorem 2.2 allows a diminishing change size b 2 with the sample size n as long as n b 2 2 → ∞. Note that no consistency result is provided in Shao and Zhang (2010) 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 τ due to the complicated nature of the self-normalizer V n,δ 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 (Olshen et al., 2004) , which is common in real data. Several variants of BS have been proposed to address this drawback, such as wild binary segmentation (WBS) (Fryzlewicz (2014)) and Narrowest-Over-Threshold (NOT) (Baranowski et al. (2019) ). Since NOT is shown to J o u r n a l P r e -p r o o f 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 } n t=1 instead of the entire sample {Y t } n t=1 . 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.

M } as the set of M random intervals such that each pair of integers (s i , e i ) are drawn uniformly from {1, · · · , n} and satisfy 1 ≤ s i < e i ≤ n and e i − s i + 1 ≥ 2h. For each random interval (s, e) ∈ F M n , we calculate the SN test

SN-NOT finds the narrowest interval (s, e) ∈ F M n where the test statistic G n,δ (s, e) exceeds a given threshold ζ n and estimates the change-point as τ = arg max k∈{s+h−1,··· ,e−h} T n,δ (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 τ 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 1. We propose to select the threshold ζ n as follows. Generate B sequences of i.i.d N (0, 1) random variables {ε b t } n t=1 , b = 1, · · · , B; for the bth sample, we calculate

The threshold ζ n is set as the 95% sample quantile of {ζ b n } B b=1 . Since the SN test statistic is asymptotically pivotal, this threshold is expected to well approximate the 95% quantile of 3.1. Testing size and Power. We generate the data from model (1.1) with sample size n = 100, 500 and 1000 respectively. For the size performance, we let β β β = (3, 0.05n) while for the power performance, we let β β β (1) = (3, 0.06n) and β β β (2) = (3 + 0.015n, 0.03n) with the change-point τ = n/2. The error process {u t } is generated via an AR(1) model where

For comparison, we also implement the supLM test defined in Andrews (1993) (using J o u r n a l P r e -p r o o f 

The error process {u t } is generated via an AR(1) model where u t = ρu t−1 + e t , e t i.i.d.

∼ N (0, (1 − ρ 2 )σ 2 ) with ρ = 0, ±0.2, ±0.5 and σ = 0.15. For comparison, we also implement the multiple change-point detection procedure proposed in Bai and Perron (1998) (denoted as BP hereafter), which is the most widely used detection algorithm allowing for temporal dependence in the error term of model (1.1). BP is implemented using function 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 τ τ τ o and the set of estimated changepoints as τ τ τ , we define d 1 (τ τ τ o , τ τ τ ) = max τ 1 ∈ τ τ τ min τ 2 ∈τ τ τ o |τ 1 −τ 2 | and d 2 (τ τ τ o ,τ τ τ ) = max τ 1 ∈τ τ τ o min τ 2 ∈ τ τ τ |τ 1 − τ 2 |, where d 1 measures the over-segmentation error of τ τ τ and d 2 measures the undersegmentation error of τ τ τ . The Hausdorff distance is then defined as d H (τ τ τ o , τ τ τ ) = max(d 1 (τ τ τ o , τ τ τ ), d 2 (τ τ τ o , τ τ τ )).

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 Hubert and Arabie (1985) . We study the cumulative confirmed cases and deaths (in log scale) of each country via the piecewise linear trend model (1.1), where given {Y t }, the change-points (τ 1 , · · · , τ m )

is estimated by the SN-NOT algorithm. An OLS is then used to recover the linear model for the ith estimated segment {Y t } τ i t= τ i−1 +1 , i = 1, 2, · · · , m + 1. With a slight abuse of notation, denote b i as the estimated slope for the ith segment. We define the normalized slope S i = b i /n for each segment. As can be seen from (1.1), the normalized slope S i measures E[Y t+1 − Y t ] for the ith 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.

Detailed analysis of cumulative confirmed cases in 8 representative countries. 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 4 .1 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 m+1 ) are also presented. In addition, we report the lag-1 sample autocorrelation ρ of the error process. From the table, we can see all of these countries have been affected J o u r n a l P r e -p r o o f 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. 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 2 . 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 2 A member of the Shincheonji religious organization was diagnosed as 31st case in Daegu, see https:

//foreignpolicy.com/2020/02/27/coronavirus-south-korea-cults-conservatives-china/ J o u r n a l P r e -p r o o f 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 4 .1, 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.

Analysis of cumulative confirmed cases in 30 countries. 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≤i≤ m+1 n −1 b i and S cur = n −1 b m+1 respectively. Combined together, the two measures allow us to obtain an overall picture of the phase when J o u r n a l P r e -p r o o f 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 4 .2, we plot S max against S cur for each country. Note that by their relative positions in Figure 4 .2, 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 To take the time factor into consideration, in Figure 4 .3, we plot the ratio S cur /S max against the days in between (i.e. τ cur − τ max with τ max as the start date for the segment J o u r n a l P r e -p r o o f with the largest slope and τ cur = τ 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 τ cur − τ 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 4 .3, 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. to COVID-19 vary from nation to nation, thus comparative analysis across countries should be interpreted with caution. Table 4 .2 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 4 .1.

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 4 .1. 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 4 .1, 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 3 . Figure 4 .5 plots S max against S cur for each country. Similar to the results 2020-05-26/spanish-health-ministry-lowers-coronavirus-death-toll-by-nearly-2000.html. for confirmed cases in Figure 4 .2, 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. 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. 5.1. Method. 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, Pesaran and Timmermann (2002) and Bauwens et al. (2015) .

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 } n t=1 , a piecewise linear trend model is estimated via SN-NOT with change-points τ τ τ . In the second stage, a flexible function f (t) is fitted on the last segment {Y t } n t= τ m +1 with the assumption that E(Y t ) = f (t) and the k-day ahead forecast for cumulative deaths can be readily made via extrapolation of 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 } n t=1 exhibits relatively stable behavior and thus facili-tates 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) = L 1 + exp − α(t/n − t 0 ) .

Based on {Y t } n t= τ m +1 , 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 Death n+k = exp( 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 5 . 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 6 available on the CDC website: "LANL" by Los Alamos National Laboratory (2020), "Imperial" by Unwin et al. (2020) , "UT" by 5 https://www.cdc.gov/coronavirus/2019-ncov/covid-data/forecasting-us.html 6 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.

J o u r n a l P r e -p r o o f University of Texas (2020), "YYG" by Gu (2020) and "MOBS" by Laboratory for the Modeling of Biological and Socio-technical Systems (2020). 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"). (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 J o u r n a l P r e -p r o o f it can serve as a meaningful addition to the existing set of forecasting models for tracking the COVID-19 pandemic. J o u r n a l P r e -p r o o f

Tests for parameter instability and structural change with unknown change point

Structural breaks in time series

Least squares estimation of a shift in linear processes

Estimation of a change point in multiple regression models

Estimating and testing linear models with multiple structural changes

Narrowest-over-threshold detection of multiple change points and change-point-like features

The contribution of structural break models to forecasting macroeconomic series

Parametric statistical change point analysis: with applications to genetics, medicine, and finance

Subspace change-point detection: A new model and solution

Multiple change-point detection for high-dimensional time series via sparsified binary segmentation

Nonmonotonic power for tests of mean shift in a time series

An empirical bayesian analysis of simultaneous changepoints in multiple data sequences

Wild binary segmentation for multiple change-point detection

Detection of change in the spatiotemporal mean function

YYG Model

Comparing partitions

A new asymptotic theory for heteroskedasticity-autocorrelation robust tests

Laboratory for the Modeling of Biological and Socio-technical Systems (2020)

LANL Model

Detecting changes in slope with an L0 penalty

Circular binary segmentation for the analysis of array-based dna copy number data

Dealing with structural breaks

Market timing and return prediction under model instability

A self-normalized approach to confidence interval construction in time series

Self-normalization for time series: a review of recent developments

Testing for change points in time series

Selective review of offline change point detection methods

The University of Texas COVID-19 Modeling Consortium

Report 23 -state-level tracking of covid-19 in the united states

Nonlinear system theory: Another look at dependence

Limit theorems for iterated random functions

Unsupervised self-normalized change-point testing for time series

Inference for linear models with dependent errors