1. Introduction. 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 (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, 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
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.
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, which partiallyprovides evidence that strict social distancing rules help slow down the virusgrowth and flatten the curve. Moreover, our analysis further indicates that compared to developed countries, most developing countries are still in the earlystages 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.
2. Methodology. 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 
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.
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,
t=i . 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 (2.4) G n = max k∈{h,··· ,n−h} T n,δ (k), T n,δ (k) = D n (1, k, n) ⊤ V n,δ (1, k, n) −1 D n (1, k, n).
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. 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). 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, lim n→∞ P(| κ − κ| < η) = 1.
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 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 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. 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 Table 3 .1 Size and size-adjusted power for SN test and supLM test. 3.2. Multiple change-point estimation. We examine the numerical performance of SN-NOT by considering the following DGP with n = 100:
.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
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 obtain the data from 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 1 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.
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 (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
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 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. 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/ 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. 
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 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 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. 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 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. 
Disease Control and Prevention (CDC) 4 , 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. 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
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.
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"). Table 5 .1 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. Let X n ∈ R d with dimension d > 0 be a set of random vector defined in a probability space
(Ω, P, F). For a corresponding set of constants a n , we say X n = O s p (a n ) if for any ε > 0, there exists a finite M > 0 and a finite N > 0 such that for n > N , P( X n /a n d > M ) + P( X n /a n d < 1/M ) < ε, where d denotes the L d norm.
Proof of Theorem 2.1 (i) It is a direct application of Theorem 3.1 in Rho and Shao (2015) and continuous mapping theorem. In particular, the result of (i) in Theorem 3.1 in Rho and Shao (2015) corresponds to the case of (i) in Assumption 2.1 for linear processes while the result of (ii) in Theorem 3.1 in Rho and Shao (2015) corresponds to the case of (ii) in Assumption 2.1 for nonlinear processes.
(ii) On one hand, note that the continuous mapping theorem indicates that L n,δ (1, τ τ τ , n) ⇒ Γ 2 L δ (κ), and R n,δ (1, τ τ τ , n) ⇒ Γ 2 R δ (κ). and it follows that V n,δ (1, τ τ τ , n) ⇒ Γ 2 V δ (κ).
On the other hand,
and it is clear that
Then the continuous mapping theorem indicates that (n b 2 2 ) −1 D n (1, τ τ τ , n) ⊤ V n,δ (1, τ τ τ , n) −1 D n (1, τ τ τ , n)
Here the last equality uses the fact that RHS of (A.1) is greater than 0 with probability 1, or equivalently, L δ (κ) and R δ (κ) is positive definite with probability 1, which will hold by similar arguments in Lemma A.1 using CauchySchwarz inequality.
Observe that max k T n,δ (k) ≥ D n (1, τ τ τ , n) ⊤ V n,δ (1, τ τ τ , n) −1 D n (1, τ τ τ , n) = O s p (n b 2 2 ). The result follows by noting n b 2 2 → L and L → ∞,
Proof of Theorem 2.2 Note that by (A.1), we have shown that with probability 
n,η }, we have
Let ν = lim 
Next, since k < τ τ τ − nη, we decompose R n,δ (1, k, n) by
:=R n,δ,1 (1, k, n) + R n,δ,2 (1, k, n).
It follows easily that V n,δ (1, k, n) −1 ≤ R n,δ (1, k, n) −1 ≤ R n,δ,2 (1, k, n) −1 where for semipositive definite matrices A and B, A ≤ B indicates B − A is semi-positive definite.
In addition, we have
where for r ∈ (κ, 1) uniformly, we can show
Therefore, if follows that
By Lemma A.1, when ν < κ, R δ,2 (ν) is invertible, hence (n b 2 2 ) −1 D n (1, k, n) ⊤ V n,δ (1, k, n) −1 D n (1, k, n)
Lemma A.1. R δ,2 (ν), defined in (A.3), is invertible for ν < κ and b 2 = 0.
Note that
Therefore we obtain
Then, since R δ,2 (ν) = 1−δ κ+δ (g 1 (r, ν, κ, b 1 , b 2 ), g 2 (r, ν, κ, b 1 , b 2 ) ⊤ (g 1 (r, ν, κ, b 1 , b 2 ), g 2 (r, ν, κ, b 1 , b 2 ) dr, the invertibility of R δ,2 (ν) is equivalent to that det(R δ,2 (ν)) > 0 (as R δ,2 (ν) is clearly semi-positive definite), i.e.
which is implied by CauchySchwarz inequality as long as
is not a constant for all r ≥ κ.
To see this, suppose R δ,2 (ν) is not invertible, then (A.5) is a constant for all r ≥ κ. Note that the numerator and the denominator of RHS of (A.5) can be written in a quadratic form of r as 2w 1 r 2 + (2νw 1 − 3w 2 )r + (2ν 2 w 1 − 3νw 2 ), (A.6) 0r 2 − 6w 1 r + (−6νw 1 + 12w 2 ), (A.7) respectively.
Therefore, comparing coefficients of the quadratic functions (A.6) and (A.7) w.r.t r, it must hold that w 1 = 0, and hence w 2 = 0, i.e.
Solving these equations for b 1 and b 2 we obtain that b 1 = b 2 = 0, contradiction.
Hence, R δ,2 (ν) is invertible.
In this section, we extend the piecewise linear structure in model (1.1) of the main text to a piecewise polynomial structure. We further apply a piecewise quadratic trend model to analyze the cumulative confirmed cases in 8 representative countries as in Section 4.2.
We extend the piecewise linear trend model (1.1) by allowing higher order polynomial terms. Specifically, let the time series
where F (p) t = (1, t/n, · · · , (t/n) p ) ⊤ and β β β t = (β 0,t , · · · , β p,t ) ⊤ are the coefficients at time t with fixed p ≥ 1. Same as in model (1.1), {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. Model (B.1) extends the piecewise linear model by allowing for polynomial trends and provides more flexibility of modeling observations in each segment.
The estimation procedure of model (B.1) is essentially the same as the one for model (1.1). Given the grid parameter ǫ, we let h = ⌊ǫn⌋. Define F (p) (s) = (1, s, · · · , s p ) ⊤ . For 1 ≤ i < j ≤ n, we denote β β β i,j = j t=i F (p) (t/n)F (p) (t/n) ⊤ −1 j t=i F (p) (t/n)Y t as the OLS estimator of β β β based on {Y t } j t=i . Let the trimming parameter satisfy 0 ≤ δ < ǫ/2. 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 (p) n (t 1 , k, t 2 ), and the self-normalizer V (p) n,δ (t 1 , k, t 2 ) = L (p) n,δ (t 1 , k, t 2 ) + R (p) n,δ (t 1 , k, t 2 ) in the same spirit as (2.1), (2.2) and (2.3) by:
n (t 1 , k, t 2 ) = (k − t 1 + 1)(t 2 − k) (t 2 − t 1 + 1) 3/2 ( β β β t 1 ,k − β β β k+1,t 2 ),
n,δ (t 1 , k, t 2 ) = k−p−1−⌊nδ⌋ i=t 1 +p+⌊nδ⌋ (i − t 1 + 1) 2 (k − i) 2 (k − t 1 + 1) 2 (t 2 − t 1 + 1) 2 ( β β β t 1 ,i − β β β i+1,k ) ⊗2 , R (p) n,δ (t 1 , k, t 2 ) = t 2 −p−⌊nδ⌋ i=k+2+p+⌊nδ⌋ (i − 1 − k) 2 (t 2 − i + 1) 2 (t 2 − t 1 + 1) 2 (t 2 − k) 2 ( β β β i,t 2 − β β β k+1,i−1 ) ⊗2 .
Then, the test statistic targeting against the one change-point alternative is defined as: The proof is a simple extension of Appendix A, hence omitted.
B.2. Analysis of cumulative confirmed cases in 8 representative countries.
We use the piecewise quadratic trend model, i.e. model (B.1) with p = 2 to re-analyze the cumulative confirmed cases in the 8 countries as in Section 4.2. Figure B .1 gives the estimated models for each country. As can be seen, compared to Figure 4 .1 in the main text, the estimated number of change-points decreases for every country, which is intuitive as more flexibility is brought into the model. For most countries, a piecewise quadratic model with one or two change-points fits the data reasonably well.
However, compared to the piecewise linear trend model, the quadratic model losses its interpretability as the parameters of each segment cannot be naturally linked to growth rate. Thus the meaning of "change-point" needs a more delicate definition. Moreover, within each segment, the growth rate of the virus still changes from day to day, making it difficult to interpret the behavior of the estimated segments. For example, we find that most estimated change-points can hardly be associated with the initiations of emergency public health measures, as the intervention effect may have been absorbed into the quadratic function. Therefore, we prefer the piecewise linear trend model for the analysis. 