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.
76,642 characters · 16 sections · 71 citation commands
Sparse HP Filter: Finding Kinks in the COVID-19 Contact Rate
\thispagestyle{empty}
\onehalfspacing
\setcounter{page}{1} \pagenumbering{arabic}
Since March 2020, there has been a meteoric rise in economic research on COVID-19. New research outputs have been appearing on the daily and weekly basis at an unprecedented level.\footnote{The major outlets for economists are: arXiv working papers, NBER working papers, and CEPR's new working paper series called “Covid Economics: Vetted and Real-Time Papers” among others.} To sample a few, Ng:NBER quantified the macroeconomic impact of COVID-19 by using data on costly and deadly disasters in recent US history; Manski:2020:JoE and Manski:NBER applied the principle of partial identification to the infection rate and antibody tests, respectively; CKS:2020 used the US state-level data to study determinants of social distancing behavior.
Across a wide spectrum of research, there is a rapidly emerging strand of literature based on a Susceptible-Infected-Recovered (SIR) model and its variants Hethcote:2000. Many economists have embraced the SIR-type models as new tools to study the COVID-19 pandemic. Avery:NBER provided a review of the SIR models for economists, calling for new research in economics. A variety of economic models and policy simulations have been built on the SIR-type models. See Acemoglu:NBER, Alvarez:NBER, Atkeson:NBER, Eichenbaum:NBER, Pindyck:NBER, Stock:NBER, kim2020estimating, and Toda among many others.
One of central parameters in the SIR-type models is the contact rate, typically denoted by $\beta$.\footnote{It is also called the transmission rate by Stock:NBER.} It measures “the average number of adequate contacts (i.e., contacts sufficient for transmission) of a person per unit time” Hethcote:2000. The contact number $\beta/\gamma$ is the product between $\beta$ and the average infectious period, denoted by $1/\gamma$; the contact number is interpreted as “the average number of adequate contacts of a typical infective during the infectious period” Hethcote:2000.
The goal of this paper is to estimate the time-varying COVID-19 contact rate, say $\beta_t$. In canonical SIR models, $\beta$ is a time-constant parameter. However, it may vary over time due to multiple factors. For example, as pointed by Stock:NBER, self-isolation, social distancing and lockdown may reduce $\beta$. To estimate a SIR-type model, FVC:ver2 allowed for a time-varying contact rate to reflect behavioral and policy-induced changes associated with social distancing. In particular, they estimated $\beta_t$ using data on deaths at city, state and country levels. Their main focus was to simulate future outcomes for many cities, states and countries.
Researchers have also adopted nonlinear time-series models from the econometric toolbox. For example, Li:Linton analyzed the daily data on the number of new cases and the number of new deaths with a quadratic time trend model in logs. Their main purpose was to estimate the peak of the pandemic. Liu2020 studied the density forecasts of the daily number of active infections for a panel of countries/regions. They modeled the growth rate of active infections as autoregressive fluctuations around a piecewise linear trend with a single break. Hartl:et:al used a linear trend model in logs with a trend break to fit German confirmed cases. Harvey:Kattuman used a Gompertz model with a time-varying trend to fit and forecast German and UK new cases and deaths.
In this paper, we aim to synthesize the time-varying contact rate with nonparametric time series modeling. Especially, we build a new nonparametric regression model for $\beta_t$ that allows for a piecewise linear trend with multiple kinks at unknown dates. We analyze daily data from Johns Hopkins University Center for Systems Science and Engineering JHU and suggest a particular transformation of data that can be regarded as a noisy measurement of time-varying $\beta_t$. Our measurement of $\beta_t$, which is constructed from daily data on confirmed, recovered and deceased cases, is different from that of FVC:ver2 who used only death data. We believe both measurements are complements to each other. However, the SIR model is at best a first-order approximation to the real world; a raw series of $\beta_t$ would be too noisy to draw on inferences regarding the underlying contact rate. In fact, the raw series exhibits high degrees of skewness and time-varying volatility even after the log transformation.
To extract the time-varying signal from the noisy measurements, we consider nonparametric trend filters that produce possibly multiple kinks in $\beta_t$ where the kinks are induced by government policies and changes in individual behavior. A natural candidate method that yields the kinks is $\ell_1$ trend filtering kim2009. However, $\ell_1$ trend filtering is akin to LASSO; hence, it may have a problem of producing too many kinks, just like LASSO selects too many covariates. In view of this concern, we propose a novel filtering method by adding a constraint on the maximum number of kinks to the popular HP (HP) filter. It turns out that this method produces a smaller number of the kink points than $\ell_1$ trend filtering when both methods fit data equally well. In view of that, we call our new method the sparse HP filter. We find that the estimated kinks are well aligned with actual events in each country. To document and monitor outbreaks of COVID-19, we propose to use piecewise constant contact growth rates using the piecewise linear trend estimates from the sparse HP filter. They provide not only an informative summary of past outbreaks but also a useful surveillance measure.
The remainder of the paper is organized as follows. In Section (ref), we describe a simple time series model of the time-varying contact rate. In Section (ref), we introduce two classes of filtering methods. In Section (ref), we have a first look at the US data, as a benchmark country. In Section (ref), we present empirical results for five countries: Canada, China, South Korea, the UK and the US. In Section (ref), we establish risk consistency of both the sparse HP and $\ell_1$ trend filters. Section (ref) concludes and appendices include additional materials. The replication R codes for the empirical results are available at \url{https://github.com/yshin12/sparseHP}. Finally, we add the caveat that the empirical analysis in the paper was carried out in mid-June using daily observations up to June 8th. As a result, some remarks and analysis might be out of sync with the COVID-19 pandemic in real time.
In this section, we develop a time-series model of the contact rate. Our model specification is inspired by the classical SIR model which has been adopted by many economists in the current coronavirus pandemic.
We start with a discrete version of the SIR model, augmented with deaths, adopted from Pindyck:NBER:
where the (initial) population size is normalized to be 1, $S_t$ is the proportion of the population that is susceptible, $I_t$ the fraction infected, $D_t$ the proportion that have died, and $R_t$ the fraction that have recovered. The parameter $\gamma = \gamma_r + \gamma_d$ governs the rate at which infectives transfer to the state of being deceased or recovered.
In the emerging economics literature on COVID-19, the contact rate $\beta$ is viewed as the parameter that can be affected by changes in individual behavior and government policies through social distancing and lockdown. We follow this literature and let $\beta = \beta_t$ be time-varying.
Let $C_t$ be the proportion of confirmed cases, that is $C_t = I_t + R_t + D_t$. In words, the confirmed cases consist of actively infected, recovered and deceased cases. Use the equations in (ref) to obtain
Assume that we have daily data on $\Delta C_t$, $\Delta R_t$ and $\Delta D_t$. From these, we can construct cumulative $C_t$, $R_t$ and $D_t$. Then $S_t = 1 - C_t$ and $I_t = C_t - R_t - D_t$. This means that we can obtain time series of $\beta_t$ from $Y_t$. We formally assume this in the following.
By Assumption (ref), we can construct $Y_t = \Delta C_{t}/(I_{t-1} S_{t-1})$. Assumption (ref) is a key assumption in the paper. We use daily data from JHU CSSE and they are subject to measurement errors, which could bias our estimates. In Appendix A, we show that the time series model given in this section is robust to some degree of under-reporting of confirmed cases. However, our estimates are likely to be biased if the underreporting is time-varying. For example, this could happen because testing capacity in many countries has expanded over the time period. Nonetheless, we believe that our measurement of $Y_t$ primarily captures the genuine underlying trend of $\beta_t$. Moreover, because the SIR model in (ref) is at best a first-order approximation, a raw series of $Y_t$ would be too noisy to be used as the actual series of the underlying contact rate $\beta_t$. In other words, $\beta_t \neq Y_t$ in actual data and it would be natural to include an error term in $Y_t$. Because $\beta_t$ has to be positive, we adopt a multiplicative error structure and make the following assumption.
Define
Under Assumption (ref), (ref) can be rewritten as
The time-varying parameter $\log \beta_t$ would not be identified without further restrictions. Because it is likely to be affected by government policies and cannot change too rapidly, we will assume that it follows a piecewise trend:
The main goal of this paper is to estimate $\log \beta_t$ and its kinks under Assumptions (ref), (ref) and (ref).
We consider two different classes of trend filtering methods to produce piecewise estimators of $f_{0,t} := \log \beta_t$. The first class is based on $\ell_1$ trend filtering, which has become popular recently. See, e.g., kim2009, tibshirani2014, and wang2016trend among others.
The starting point of the second class is the HP filter, which has been popular in macroeconomics and has been frequently used to separate trend from cycle. The standard convention in the literature is to set $\lambda = 1600$ for quarterly time series. For example, Ravn:02 suggested a method for adjusting the HP filter for the frequency of observations; deJong:16 and Explicit:HP established some representation results; Hamilton:2018 provided criticism on the HP filter; phillips2019boosting advocated a boosted version of the HP filter via $L_2$-boosting bHP that can detect multiple structural breaks. We view that the kinks might be more suitable than the breaks for modelling $\beta_t$ using daily data. It is unlikely that in a few days, the degree of contagion of COVID-19 would be diminished with an abrupt jump by social distancing and lockdown. The original HP filter cannot produce any kink just as ridge regression does not select any variable. We build the sparse HP filter by drawing on the recent literature that uses an $\ell_0$-constraint or -penalty bertsimas2016,chen2018,chen2018arXiv,Huang:2018.
In $\ell_1$ trend filtering, the trend estimate $f_t$ is a minimizer of
which is related to HP filtering; the latter is the minimizer of
In this paper, the main interest is to find the kinks in the trend. For that purpose, $\ell_1$ trend filtering is more suitable than the HP filtering. The main difficulty of using (ref) is the choice of $\lambda$. This is especially challenging since the time series behavior of $y_t$ is largely unknown.
The $\ell_1$ trend filter is akin to LASSO. In view of an analogy to square-root LASSO sqrt-lasso, it might be useful to consider a square-root variant of (ref):
We will call (ref) square-root $\ell_1$ trend filtering. Both (ref) and (ref) can be solved via convex optimization software, e.g., CVXR CVXR.
As an alternative to $\ell_1$ trend filtering, we may exploit Assumption (ref) and consider an $\ell_0$-constrained version of trend flitering:
The formulation in (ref) is related to the method called best subset selection bertsimas2016,chen2018. It requires only the input of $\kappa$. However, because of the nature of the $\ell_0$-(pseudo)norm, it would not work well if the signal-to-noise ratio (SNR) is low hastie2017extended,mazumder2017subset. This is likely to be a concern for our measurement of the log contact rate.
To regularize the best subset selection procedure, it has been suggested in the literature that (ref) can be combined with $\ell_1$ or $\ell_2$ penalization bertsimas2020,mazumder2017subset. We adopt bertsimas2020 and propose an $\ell_0$-constrained version of the Hodrick-Prescott filter:
As in (ref), the tuning parameter $\kappa$ controls how many kinks are allowed for. Thus, we have a direct control of the resulting segments of different slopes. The $\ell_2$ penalty term is useful to deal with the low SNR problem with the COVID 19 data. We will call (ref) sparse HP trend filtering.
Problem (ref) can be solved by mixed integer quadratic programming (MIQP). Rewrite the objective function in (ref) as
subject to $z_t \in \{0, 1\}, t=2,\ldots,T-1$, $\underline{f} \leq f_t \leq \overline{f}$, $\sum_{t=2}^{T-1} z_t \le \kappa$, and
This is called a big-M formulation that requires that \[ \max_t | f_{t-1} - 2 f_t + f_{t+1} | \leq M. \] We need to choose the auxiliary parameters $\underline{f}$, $\overline{f}$ and $M$. We set $\underline{f} = \min y_t$ and $\overline{f} = \max y_t$. One simple practical method for choosing $M$ is to set
To implement the proposed method, it is simpler to write the MIQP problem above in matrix notation. Let $\bm{y}$ denote the $(T \times 1)$ vector of $y_t$'s and $\bm{1}$ a vector of 1's whose dimension may vary. We solve
subject to $\bm{z} \in \{0, 1\}^{T-2}$, $\underline{f} \bm{1} \leq \bm{f} \leq \overline{f} \bm{1}$, $\bm{1}^\top \bm{z} \leq \kappa$, $-M \bm{z} \leq \bm{D} \bm{f} \leq M \bm{z}$, where $\bm{D}$ is the $(T-2) \times T$ second-order difference matrix such that $$ \bm{D} = \left[
\right] $$ with entries not shown above being zero. Let $\widehat{\bm{f}}$ and $\widehat{\bm{z}}$ denote the resulting maximizers. It is straightforward to see that $\widehat{\bm{f}}$ also solves \eqref{sparseHP}. Therefore, $\widehat{\bm{f}}$ is the $(T \times 1)$ vector of trend estimates and $\widehat{K} := \{ t =2,\ldots,T-1: \widehat{z}_t = 1 \}$ is the index set of estimated kinks. The MIQP problem can be solved via modern mixed integer programming software, e.g., \textbf{Gurobi}. Because the sample size for $y_t$ is typically less than 100, the computational speed of MIQP is fast enough to carry out cross-validation to select tuning parameters. We summarize the equivalence between the original and MIQP formulation in the following proposition.
We first consider the sparse HP filter. There are two tuning parameters: $\lambda$ and $\kappa$. It is likely that there will be an initial stage of coronavirus spread, followed by lockdown or social distancing. Even without any policy intervention, it will come down since many people will voluntarily select into self-isolation and there is a chance of herd immunity. Hence, the minimum $\kappa$ is at least 1. If $\kappa$ is too large, it is difficult to interpret the resulting kinks. In view of these, we set the possible values $\kappa \in \mathcal{K} = \{2,3,4\}$. For each pair of $(\kappa,\lambda)$, let $\widehat{\bm{f}}_{-s}(\kappa,\lambda)$ denote the leave-one-out estimator of $\bm{f}_s$. That is, it is the sparse HP filter estimate by solving:
The only departure from (ref) is that we replace the fidelity term $\sum_{t=1}^T (y_t - f_t)^2$ with $\sum_{t=1, t \neq s}^T (y_t - f_t)^2$. We choose the optimal $(\kappa,\lambda)$ by
where $\mathcal{L}$ is the set for possible values of $\lambda$. We view $\lambda$ as an auxiliary tuning parameter that mitigates the low SNR problem. Hence, we take $\mathcal{L}$ to be in the range of relatively smaller values than the typical values used for the HP filter. In the numerical work, we let $\Lambda$ to a grid of equi-spaced points in the $\log_2$-scale.
We now turn to the HP, $\ell_1$ and square-root $\ell_1$ trend filters. For each filter, we choose $\lambda$ such that the fidelity term $\sum_{t=1}^T (y_t - f_t)^2$ is the same as that of the sparse HP filter. In this way, we can compare different methods holding the same level of fitting the data. Alternatively, we may choose $\lambda$ by leave-one-out cross validation for each filtering method. However, in that case, it would be more difficult to make a comparison across different methods. Since our main focus is to find the kinks in the contact rate, we will fine-tune all the filters to have the same level of $\sum_{t=1}^T (y_t - f_t)^2$ based on the sparse HP filter's cross validation result.
As a benchmark, we have a first look at the US data. The dataset is obtained via R package coronavirus Krispin, which provides a daily summary of COVID-19 cases from Johns Hopkins University Center for Systems Science and Engineering JHU. Following Liu2020, we set the first date of the analysis to begin when the number of cumulative cases reaches 100 (that is, March 4 for the US). To smooth data minimally, we take $Y_t$ in (ref) to be a three-day simple moving average: that is, $Y_t = (\breve{Y}_{t}+\breve{Y}_{t-1} + \breve{Y}_{t-2})/3$, where $\breve{Y}_{t}$ is the daily observation of $Y_t$ constructed from the dataset.\footnote{Liu2020 used one-sided three-day rolling averages; FVC:ver2 took 5-day centered moving averages.} Then, we take the log to obtain $y_t = \log Y_t$.
Figure (ref) has four panels. The top-left panel shows the fraction of daily positives, the top-right panel the fraction of lagged cumulative infectives, the bottom-left panel the fraction of lagged cumulative susceptibles, and the bottom-right $Y_t = \Delta C_{t}/(I_{t-1} S_{t-1})$. In the US, statewide stay-at-home orders started in California on March 20 and extended to 30 states by March 30 NYtimes. The inserted vertical line in the figure corresponds to March 30, which we will call the “lockdown” date for simplicity, although there was no lockdown at the national level. As a noisy measurement of $\beta_t$, $Y_t$ shows enormous skewness and fluctuations especially in the beginning of the study period. This indicates that the signal-to-noise ratio is high and is time-varying as well. This pattern of the data has motivated Assumption (ref). Because $S_{t-1}$ is virtually one throughout the analysis period (0.994 on June 8, which is the last date of the sample), $Y_t \approx \Delta C_{t}/I_{t-1}$, which is daily positives divided by the lagged infectives.
Figure (ref) shows the raw data along with parametric fitting. The top-left panel shows the logarithm of $Y_t$, which still exhibits some degree of skewness and time-varying variance. The fitted regression line is based on the following parametric regression model:
where $t_0$ is March 30. The simple idea behind (ref) is that an initial, time-constant contact rate began to diminish over time after a majority of US states imposed stay-at-home orders.
In simple SIR models, the contact number $\beta/\gamma$ is identical to the basic reproduction number denoted by $R_0$, which is viewed as a key threshold quantity in the sense that “an infection can get started in a fully susceptible population if and only if $R_0 > 1$ for many deterministic epidemiology models” Hethcote:2000. Since $\beta_t$ is time-varying in our framework, we may define a time-varying basic reproduction number by $R_0(t) := \beta_t/\gamma$.
The top-right panel shows the estimates of time-varying $R_0(t)$:\footnote{The formula given in (ref) is valid if errors are homoskedastic, which is unlikely to be true in actual data. However, we present (ref) here because it is simpler. Our main analysis focuses on estimation of the kinks based on $y_t$, not on estimating $R_0(t)$. We use the latter mainly to appreciate the magnitude of the kinks.}
where $\gamma = 1/18$ is taken from Acemoglu:NBER. This corresponds to 18 days of the average infectious period. The parametric estimates of $R_0(t)$ started above 4 and reached $0.15$ at the end of the sample period.
The left-bottom panel shows the residual plot in terms of $y_t$ and the right-bottom panel the residual plot in terms of $R_0(t)$. In both panels, the estimated residuals seem to be biased and show autocorrelation. Especially, the positive values of residuals at the end of the sample period is worrisome because the resulting prediction would be too optimistic.
In this section, we present estimation results for five countries: Canada, China, South Korea, the UK and the US. These countries are not meant to be a random sample of the world; they are selected based on our familiarity with them so that we can interpret the estimated kinks with narratives. We look at the US as a benchmark country and provide a detailed analysis in Section (ref). A condensed version of the estimation results for other countries are provided in Section (ref).
Figure (ref) summarizes the results of leave-one-out cross validation (LOOCV) as described in Section (ref). The range of tuning parameters were: $\kappa\in\{2,3,4\}$ and $\lambda=\{2^{0}, 2^{1}, \ldots, 2^{5}\}$. We can see that the choice of $\kappa$ seems to matter more than that of $\lambda$. Clearly, $\kappa = 2$ provides the worst result and $\kappa = 3$ and $\kappa = 4$ are relatively similar. The LOOCV criterion function was minimized at $(\widehat{\kappa},\widehat{\lambda})=(4,1)$.
Based on the tuning parameter selection in Figure (ref), we show estimation results for the sparse HP filter in Figure (ref). The structure of Figure (ref) is similar to that of Figure (ref). The top-left panel shows estimates of the sparse HP filter along with the raw series of $y_t$ and the parametric estimates shown in Figure (ref). The top-right panel displays counterparts in terms of $R_0(t)$. The bottom panels exhibit residual plots for the $\log \beta_t$ and $R_0(t)$ scales. The trend estimates from the sparse HP filter fit the data much better than the simple parametric estimates. The estimated kink dates are: March 16, March 20, April 14, and May 13. There are five periods based on them.
To provide narratives on these dates, President Trump declared a national emergency on March 13; The Centers for Disease Control and Prevention (CDC) recommended no gatherings of 50 or more people on March 15; New York City's public schools system announced that it would close on March 16; and California started stay-at-home orders on March 20 NYtimes:timeline,NYtimes. These events indicate that the second period was indeed the peak of the COVID-19 epidemic in the US. The impact of social distancing and stay-at-home orders across a majority of states is clearly visible in the third period. The fourth and fifth periods include state reopening: for example, stay-at-home order expired in Georgia and Texas on April 30; in Florida on May 4; in Massachusetts on May 18; in New York on May 28 NYtimes:reopening. In short, unlike the parametric model with a single kink, the nonparametric trend estimates detect multiple changes in the slopes and provide kink dates, which are well aligned with the actual events.
We now turn to different filtering methods. In Figure (ref), we show selection of $\lambda$ for the HP, $\ell_1$ and square-root $\ell_1$ filters. As explained in Section (ref), the penalization parameter $\lambda$ is chosen to ensure that all different methods have the same level of fitting the data. Figure (ref) shows the estimation results for the HP filter. The HP trend estimates trace data pretty well after late March, as clear in residual plots. However, there is no kink in the estimates due to the nature of the $\ell_2$ penalty term in the HP filter. The tuning parameter was $\lambda = 30$, which is 30 times as large as the one used in the sparse HP filter. This is because for the HP filter, $\lambda$ is the main tuning parameter; however, for the sparse HP filter, $\lambda$ plays a minor role of regularizing the $\ell_0$ constrained method.
In Figure (ref), we plot estimation results using $\ell_1$ trend filtering. The results look similar to those in Figure (ref), but there are now 10 kink points: March 7, March 15, March 16, March 20, March 21, March 30, April 14, April 21, May 12, May 27. They are dates $t$ such that $|\Delta^2 \log \widehat{\beta}_t| > \eta$, where $\Delta^2$ is the double difference operator and $\eta=10^{-6}$ is an effective zero.\footnote{The results are robust to the size of the effective zero and do not change even if we set $\eta=10^{-3}$. Gurobi used for the Sparse HP filtering also imposes some effective zeros in various constraints. We use the default values of them. For example, the integer tolerance level and the general feasibility tolerance level are $10^{-5}$ and $10^{-6}$, respectively.} The tuning parameter $\lambda=0.9$ was chosen by minimizing the distance between the fidelity of $\ell_1$ and that of the Sparse HP. Recall that the sparse HP filter produces the kinks on March 16, March 20, April 14, and May 13. In other words, the $\ell_1$ filter estimates 6 more kinks than the sparse HP filter when both fit the data equally well. It is unlikely that two adjacent dates (March 15-16 and March 20-21) correspond to two different regimes in the time-varying contact rate. This suggests that the $\ell_1$ filter may over-estimate the number of kinks. Figure (ref) shows estimation results for the square-root $\ell_1$ trend filters. The chosen $\lambda=0.5$ was smaller than that of the $\ell_1$ trend filter due to the change in the scale of the fidelity term; however, the trend estimates look very similar and the estimated kinks are identical between the $\ell_1$ and square-root $\ell_1$ trend filters. In Figure (ref), we plot the sparse HP filter estimates along with $\ell_1$ filter estimates. Both methods have produced very similar trend estimates, but the number of kinks is substantially different: only 4 kinks for the sparse HP filter but 10 kinks for the $\ell_1$ filter.
In this section, we provide condensed estimation results for other countries. We focus on the sparse HP and $\ell_1$ filters whose tuning parameters are chosen as in the previous section. Appendices (ref) and (ref) contain the details of the selection of tuning parameters.
Figure (ref) shows the empirical results of Canada. The estimated kink dates are: March 18 and April 11. Based on them, we can classify observations into three periods:
Quebec and Ontario are the two provinces hardest hit by COVID-19. In Quebec, daycares, public schools, and universities are closed on March 13 followed by non-essential businesses and public gathering places on March 15. Montreal declared state of emergency on March 27 CTV-Quebec:timeline. Similarly, all public schools in Ontario are closed on March 12. The state of emergency was announced in Ontario on March 17 and ordered to close all non-essential businesses on March 23 Global-Ontario:timeline. We set the lockdown date in Canada on March 13 as other provincial governments as well as the federal government started to recommend the social distancing measures strongly along with the cancellation of various events on the date CBC:March13. These tight lockdown and social distancing measures seemed to contribute the sharp decline of the contact rate in the second period. Both governments started to announce the plans to lift the lockdown measures at the end of April, which corresponds to the third period. Lockdown fatigue would also cause the slower decrease of the contact rate. In sum, a series of social distancing measures have been effective to decrease the contact rate but with some lags. The sparse HP filtering separates these periods reasonably well. However, the $\ell_1$ filtering overfits the model with 5 kinks.
Figure (ref) shows the results for China. Since the pandemic is almost over in China, we use the data censored on April 26th when the 3-day-average of newly confirmed cases is less than 10. The estimated kink dates are: January 28, March 14, March 24, and April 18. Based on them, we can classify observations into five periods:
Figure (ref) shows the results for South Korea. For the same reason in China, we use the data censored on April 29. The estimated kink dates are: March 3, March 15, April 2, and April 21. Based on them, we can classify observations into five periods:
Figure (ref) shows the empirical results of the UK. The estimated kink dates are: March 12 and March 14. Based on them, we can classify observations into three periods:
Overall, the trend of the contact rate is quite similar to those of the US and Canada. The location of the kinks are around more in the initial periods but it shows the steady downward trending after the prime minister's lockdown announcement. This results in a smooth curve in the $R_0(t)$ scale. The trend estimates of the $\ell_1$ filter is almost identical to those of the sparse HP filter; however, it indicates 10 kinks, which seem overly excessive.
The sparse HP filter produces the kinks where the slope changes in the $\log \beta_t$ scale, thereby providing a good surveillance measure for monitoring the ongoing epidemic situation. The policy responses are based on various scenarios and the contact rate is one of the most important measures that determine different developments. As a summary statistic of the time-varying contact rate, we propose to consider the time-varying growth rate of the contact rate, which we call contact growth rates:
Recall that we have defined the time-varying basic reproduction number by $R_0(t) = \beta_t/\gamma$. Because $\gamma$ is fixed over time, we have that
Therefore, $\xi(t)$ can be interpreted as the time-varying growth rate of the basic reproduction number; it does not require the knowledge of $\gamma$ and solely depends on $\beta_t$. Furthermore, by simple algebra,
which implies that $\xi(t)$ will be piecewise constant if $\log \beta_t$ is piecewise linear. This simple algebraic relationship shows that a change in the slope at the kink in the $\log \beta_t$ scale is translated to a break in the time-varying contact growth rates and therefore in growth rates of the time-varying basic reproduction number. When $\xi (t)$ is a large positive number, that will be a warning signal for the policymakers. On the contrary, if $\xi(t)$ is a big negative number, that may suggest that policy measures imposed before are effective to reduce the contagion.
Table (ref) reports the time-varying contact growth rates in the five countries that we investigate, using the sparse HP trend estimates. For the US, the explosive growth rate of 7.5% in the second period is followed by the negative growth rates of $-7.7$%, $-3.4$%, and $-1$%, albeit at diminishing magnitudes. The trajectory of Canada is similar to that of the US. The growth rates of China fluctuated up and down: it started with a high positive 15% followed by $-12$%; a sharp V-turn at the end of the second period (March 14) with the resulting explosive growth rate of 30%, followed by moderate 4% and impressive $-23$%. It might be the case that the up-and-down pattern observed in China is in part due to data quality issues since China was the first country to experience the pandemic. For South Korea, we can see the stunning drop of the growth rates culminating on March 15 (the end of the second period). A modest positive growth rate during period 3 is offset by a larger magnitude of negative growth rate in period 4. The UK has experienced steady---but not spectacular---negative growths over the sample period following a sharp fluctuation in mid-March. This hints the degrees of effectiveness of the UK lockdown policy. As early pandemic epicenters, China and South Korea experienced V-turns in the time-varying growth rates of basic reproduction number. Canada, the UK and the US may face similar trajectories as they reopen their countries. Our surveillance statistic can be a useful indicator to monitor a new outbreak of COVID-19. However, it will be mainly useful for a short-term projection of the contact growth rate because it is not designed to make long-term trend predictions.
In this section, we examine theoretical properties of the sparse HP and $\ell_1$ filters in terms of risk consistency. Let $\| \cdot \|_0$ denote the usual $\ell_0$-(pseudo)norm, that is the number of nonzero elements, and let $\| \cdot \|_r$ and $\| \cdot \|_{\infty}$, respectively, denote the $\ell_r$ norm for $r=1,2$ and the sup-norm.
Define
where $M$ is defined in $\eqref{def:M}$. For each $\bm{f} \in \mathcal{F}$, define
Let $\bm{f^*}$ denote the ideal sparse filter in the sense that
Let $\bm{\widehat{f}}$ denote the sparse HP filter defined in Section (ref). Then,
is always nonnegative. Following the literature on empirical risk minimization, we bound the excess risk $R$ in (ref) and establish conditions under which it converges to zero.
Recall that the sparse HP filter minimizes
subject to $\bm{f} \in \mathcal{F}$.
Let $S_n( \bm{f} ) := T^{-1} (\bm{y}- \bm{f})^\top (\bm{y}- \bm{f})$. Write
Therefore, it suffices to bound two terms above. For the second term, we can use (ref) and (ref) to bound
We summarize discussions above in the following lemma.
To derive an asymptotic result, we introduce subscripts indexed by the sample size $ T $, when necessary for clarification. Let $ \mathcal{G}_{\kappa} $ denote the set of every continuous and piecewise linear function whose slopes and the function itself is bounded by $ C_1 $ and $ C_2 $, respectively, and the number of kinks is bounded by $ \kappa $.
Then, we have the following proposition.
Theorem (ref) establishes the consistency in terms of the excess risk $R$. Assumption (ref) provides sufficient conditions for (ref) and (ref). Condition (ref) is a uniform law of large numbers for the class $\mathcal{F}$ and condition (ref) imposes a weak condition on $\lambda$. Proposition (ref) follows immediately from Lemma (ref) once (ref) and (ref) are established.
The $\ell_1$ trend filtering ((ref)) can be expressed as $$ \widetilde{\bm{f}} :=\arg\min_{\bm{f}\in\mathbb R^T} \|\bm{y}-\bm{f}\|_2^2+\lambda \|\bm{D}\bm{f}\|_1. $$ We now derive the deviation bound for $\|\widetilde{\bm{f}} -\bm{f}^*\|_2.$ First, the problem is equivalent to a regular LASSO problem as stated in Lemma (ref) below.
Write $ \bm{D}= ( \bm{D}_3 , \bm{D}_2) $ where $\bm{D}_2$ has two columns. Additionally, write $$ \bm{G}_2:=
,\quad \bm{g}_1:=
, $$ where $\bm{0}$ is $2\times (T-2)$, $\bm{g}_1$ is $T\times 2$ and $\bm{G}_2$ is $T\times (T-2).$ Let $ \bm{P}_{\bm{g}_1} = \bm{g}_1(\bm{g}_1^\top \bm{g}_1)^{-1}\bm{g}_1^\top. $
Next, let $J$ denote the indices of $t$ so that $f_{0,t-1}-f_{0,t}\neq f_{0,t}-f_{0,t+1}$ when $t\in J$; let $J^c$ denote the indices of $t$ so that $f_{0,t-1}-f_{0,t}= f_{0,t}-f_{0,t+1}$ when $t\in J$. Here, $\{ f_{0,t}: t=1, \ldots,T \}$ denote the true elements of $\bm{f}$. For a generic vector $\theta\in\mathbb R^{T-2}$, let $\theta_J$ and $\theta_{J^c}$ respectively be its subvectors whose elements are in $J$ and $J^c.$ No we define the restricted eigenvalue constant $$ \zeta:=\inf_{\|\theta_{J^c}\|_1\leq 9\|\theta_J\|_1} \frac{ \|\frac{1}{\sqrt{T}} \widetilde{\bm{X}} \theta\|_2^2 }{\|\theta\|_2^2}. $$
To achieve risk consistency, $\lambda$ has to be chosen to make the second term on the right-hand side of (ref) asymptotically small and to ensure that the event $2.5 \| {\bm{u}}^\top \widetilde{\bm{X}}\|_{\infty}<\lambda$ holds with high probability. The first term on the right-hand side of (ref) will converge to zero under mild conditions on $\bm{u}$. It is reassuring that the $\ell_1$ trend filter fits COVID-19 data well in our empirical results.
In this subsection, we obtain risk consistency of $\exp(\widehat{f}_t)$ and $\exp(\widetilde{f}_t)$. To do so, we first rewrite the excess risk in (ref) as \[ R(\bm{f}, \bm{f^*})=\frac{1}{T}\sum_{t=1}^T\mathbb E (f_t-f_t^*)^2, \] for $\bm{f} = ( f_1,..., f_T)^\top$ and $\bm{f^*}=(f_1^*,...,f_T^*)^\top$. We have proved in the previous sections that $R(\bm{\widehat{f}}, \bm{f^*}) \rightarrow_P 0$ and $R(\bm{\widetilde{f}}, \bm{f^*}) \rightarrow_P 0$. Then, under the assumption that there exists a constant $C < \infty$ such that $\max_t|f_t|+\max_t|f_t^*|<C$ for all $\{f_t\}$ on the parameter space, we get, uniformly for all $f$ on the parameter space, \[ \frac{1}{T}\sum_{t=1}^T\mathbb E (\exp(f_t)-\exp(f_t^*))^2 = \frac{1}{T}\sum_{t=1}^T\mathbb E \exp(2\tilde f_t)(f_t-f_t^*)^2 \leq \exp(2C) \frac{1}{T}\sum_{t=1}^T\mathbb E (f_t-f_t^*)^2, \] where in the first equality we used the mean value theorem for some $\tilde f_t\in(f_t, f_t^*).$ Therefore, $$ R(\exp(\bm{\widehat{f}}), \exp({\bm{f^*}})) \leq \exp(2C) R(\bm{\widehat{f}}, \bm{f^*}) =o_P(1) $$ and the analogous result folds for $\exp(\bm{\widetilde{f}})$.
We have developed a novel method to estimate the time-varying COVID-19 contact rate using data on actively infected, recovered and deceased cases. Our preferred method called the sparse HP filter has produced the kinks that are well aligned with actual events in each of five countries we have examined. We have also proposed contact growth rates to document and monitor outbreaks. Theoretically, we have outlined the basic properties of the sparse HP and $\ell_1$ filters in terms of risk consistency. The next step might be to establish a theoretical result that may distinguish between the two methods by looking at kink selection consistency. It would also be important to develop a test for presence of kinks as well as an inference method on the location and magnitude of kinks and on contact growth rates. In the context of the nonparametric kernel regression of the trend function, Delgado2000 explored the distribution theory for the jump estimates but did not offer testing for the presence of a jump. Compared to the kernel smoothing approach, it is easier to determine the number of kinks using our approach, as we have demonstrated. Furthermore, the linear trend specification is more suitable for forecasting immediate future outcomes, at least until the next kink arises. The long-term prediction is more challenging and it is beyond the scope of this paper. Finally, it would be useful to develop a panel regression model for the contact rate at the level of city, state or country. These are interesting research topics for future research.