Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
70,400 characters · 9 sections · 62 citation commands
Monitoring the pandemic: A fractional filter for the COVID-19 contact rate
\thispagestyle{empty} \setcounter{page}{0} \paragraph{\bf Abstract.}
\paragraph{\bf Keywords.} COVID-19, filtering, long memory, SIR model, unobserved components.
\paragraph{\bf JEL-Classification.} C22, C51, C52.
Since the outbreak of COVID-19 reducing social contacts is widely viewed as the key way to contain the spread of the virus. In terms of the Susceptible-Infected-Recovered (SIR) model\footnote{The SIR model -- in its various variants -- has recently become a popular tool to study the economic impact of the pandemic and for policy simulations, see AceCheWe2020, AveBosCl2020, Kor2021, LiuMooSc2021 among others.}, this relates to the contact rate, defined as the average number of contacts per person per time unit multiplied by the probability of disease transmission between a susceptible and an infected individual Het2000. The probability of disease transmission should only depend on characteristics that are specific to the virus. Therefore, the contact rate can be interpreted as a proxy for aggregate social behavior and is the key variable addressed by social distancing measures. Knowing the trajectory of the contact rate would allow to draw inference on the impact of policy measures on contact reduction, to real-time monitor the dynamics of virus dispersion, and to design policy rules based on the current pandemic situation. Since the contact rate itself is unobservable, appropriate methods to estimate the contact rate are required, and will be considered in this paper.
At the early stage of the pandemic, first estimates for the natural logarithm of the contact rate were obtained by fitting a deterministic, linear time trend with structural breaks to transformations of data on confirmed, recovered, and deceased cases HarWaeWe2020, LeeLiaSe2021, LiuMooSc2021. Modeling the log contact rate by a piece-wise linear time trend was a reasonable and pragmatic approximation given the short time series on case numbers available at that time. However, it implies that contact rate growth evolves deterministically as a straight line with jumps at the break dates. This assumption is likely to be violated by the behavior of individuals. While structural breaks may be suitable to identify turning points of the contact rate, they are inappropriate for monitoring the current pandemic situation, as breaks require at least some post-break observations to be well identified.
This paper aims to improve estimates for the contact rate of COVID-19 by taking into account key features of aggregate social behavior. In detail, the log contact rate, as denoted by $\log \beta_t$, is modeled as an unobserved, fractionally integrated process of (unknown) order $d \in \mathbb{R}_+$, generated by stochastic shocks $\{\eta_i\}_{i=1}^t$.\footnote{Fractional integration techniques have been found useful for describing the aggregate behavior of individuals in a variety of applications, e.g.\ for explaining the Deaton paradox DieRud1991 and for the estimation of the business cycle HarTscWe2020.} The stochastic specification of the contact rate is motivated by the consideration that social decisions, e.g.\ on whether to meet, are made conditional on the information available at that time, e.g.\ on current social distancing measures or the state of the pandemic. As information does not evolve deterministically but appears as stochastic shocks, this suggests to treat $\log \beta_t$ as a stochastic process generated by the information shocks $\{\eta_i\}_{i=1}^t$. Specifying $\log \beta_t$ as a fractionally integrated process accounts for strong persistence and nonstationarity (in short:\ long memory) of social behavior. In contrast to structural breaks but also to random walks, the fractional specification allows social behavior to gradually adjust to new information both at the individual and at the aggregate level. Individually, this reflects a gradual reduction or increase of contacts as new information becomes available (e.g.\ as new contact restrictions are imposed), while on aggregate it allows individuals to react heterogeneously both in terms of speed and intensity to novel information. As the persistence of the log contact rate is unknown, the integration order $d$ is treated as an unknown parameter to be estimated.
Methodologically, this paper contributes to the literature on time series filtering by setting up a novel unobserved components (UC) model that does not require prior knowledge about the integration order of the variable under study. Current UC models and related filtering techniques rely heavily on prior assumptions about the integration order $d$ and typically assume $d=1$ Har1985, MorNelZi2003, ChaMilPa2009 or $d=2$ Cla1987, HodPre1997, OhZivCr2008 to be known. In contrast, the novel UC model reflects that the degree of persistence of the log contact rate is unknown. It allows to decompose a noisy measurement for the log contact rate that is based on a transformation of data on confirmed, recovered, and deceased cases, into measurement errors, seasonal components, and the unobserved log contact rate itself. As the latter is modeled by a fractionally integrated process, the model is called the fractional UC model.
The second methodological contribution of this paper is to derive a computationally much simpler estimator for the model parameters and the unobserved components compared to current state space methods. Current methods typically rely on the Kalman filter to set up a conditional (quasi-)likelihood function for the estimation of the model parameters. Given the parameter estimates, a time-varying signal for the unobserved components is then obtained from the Kalman smoother. Both the Kalman filter and smoother become computationally infeasible when the dimension of the state vector of UC models is high, as for fractionally integrated processes. To address this problem, this paper proposes a computationally simple modification of the Kalman filter and smoother that is termed the fractional filter. While filtered and smoothed estimates from the fractional filter are identical to the Kalman filter and smoother, the fractional filter avoids the computationally intensive recursions for the conditional variance. The fractional filter provides a closed-form expression for the prediction error of UC models, based on which a conditional-sum-of-squares (CSS) estimator for the fractional integration order and other model parameters is set up. While the CSS estimator has been found useful for the estimation of ARFIMA models, see HuaRob2011 and Nie2015, it has not been considered in the UC literature so far. The CSS estimator minimizes the sum of squared prediction errors that is proportional to the exponent in the conditional (quasi-)likelihood function based on the Kalman filter. Due to the computational gains from the fractional filter, the CSS estimator allows to estimate UC models with richer long-run dynamics. The paper provides the asymptotic theory for the CSS estimator, showing it to be consistent and asymptotically normally distributed, while the finite sample properties are assessed by a Monte Carlo study.
Using data from the Johns Hopkins University Center for Systems Science and Engineering DonDuGa2020, estimates for contact and reproduction rate are presented for Canada, Germany, Italy, and the United States, where benefits from the new methods directly become apparent: First, estimation results are not only well in line with the chronology of the pandemic, but also allow to identify different contact regimes generated by the strengthening and easing of contact restrictions. Second, a recursive window evaluation shows contact rate estimates at the end of a truncated sample to largely overlap with those based on the full sample information. This makes the fractional filter a suitable candidate for monitoring outbreaks at the current frontier of the data. And third, the proposed estimation and filtering techniques are shown to be fairly robust to under-reporting of recovered cases, which is of particular importance for the US, as several states do not report data on recovered individuals. While under-reporting heavily downward-biases contact and reproduction rate estimates in LeeLiaSe2021, this is shown not to be the case for the fractional filter.
The remaining paper is organized as follows: Section (ref) motivates the specification of the contact rate and sets up the fractional UC model. Section (ref) introduces the fractional filter for $\log \beta_t$, covers parameter estimation via the CSS estimator and presents the asymptotic theory. Section (ref) contains empirical results for Canada, Germany, Italy, and the United States, while section (ref) concludes. The appendices include proofs for consistency and asymptotic normality of the CSS estimator as well as a Monte Carlo study on the finite sample properties.
To motivate the estimation of the contact rate, consider the discrete SIR model, augmented to include deaths, which also forms the starting point of Pin2020 and LeeLiaSe2021
In (ref), the (initial) population size, normalized to be one, is decomposed into $S_t$, the proportion of the population susceptible in $t$, $I_t$, the fraction of the population infected in $t$, $D_t$, the fraction that has died until $t$, and $R_t$, the proportion that has recovered until $t$. In (ref), $\gamma = \gamma_d + \gamma_r$ denotes the rate at which infected either die, see (ref), or recover, see (ref), and obviously $\gamma_d, \gamma_r \geq 0$. Thus, $\gamma I_{t-1}$ denotes the fraction of outflows of infected at $t$. The fraction of new infections at $t$ is captured by $\beta_t S_{t-1}I_{t-1}$, where $S_{t-1}I_{t-1}$ can be interpreted as the average probability of a contact being between a susceptible subject and an infected subject. $\beta_t > 0$ is called the contact (or transmission) rate. It equals the average number of contacts per person per time unit multiplied by the probability of disease transmission between a susceptible and an infectious person Het2000. As in LeeLiaSe2021, the contact rate is allowed to be time-varying. This reflects social behavior to change over time, e.g.\ in response to policy changes or to novel information on the pandemic. Since the contact rate determines inflows into infected, see (ref), it is the key variable tackled by social distancing policies.
Based on the contact rate, the reproduction rate $\mathcal{R}_{ t} = \beta_t/\gamma$ can be derived. It is the average number of infections caused by an infected subject during the infectious period $1/\gamma$ at the early stage of the pandemic (where $S_{t-1} \approx 1$). $\mathcal{R}_{t}$ is an indicator for the current dynamics of the pandemic, as for $\mathcal{R}_{t} < 1$ outflows from infected exceed inflows, causing $\Delta I_t$ to converge, see (ref) where $0 \leq S_{t-1} \leq 1$. Thus, if policy seeks to contain the spread of COVID-19, then it must control the contact rate, which controls the reproduction rate $\mathcal{R}_t$.
As shown by LeeLiaSe2021, from (ref) to (ref) a measurement for the contact rate $\beta_t$ can be obtained directly: Denote $C_t = I_t + R_t + D_t$ as the fraction of confirmed cases (consisting of infected, recovered, and deceased cases) and use $\Delta C_t = \Delta I_t + \Delta R_t + \Delta D_t$ together with (ref) to (ref) to obtain $\Delta C_t = \beta_t S_{t-1} I_{t-1} - \gamma I_{t-1} + (\gamma_d + \gamma_r) I_{t-1} = \beta_t S_{t-1} I_{t-1} $. Solving for $\beta_t$ yields
see LeeLiaSe2021. As argued there, if for each $t$ the data $(C_t, R_t, D_t)$ can be observed, then the time-varying contact rate can be calculated straightforwardly via (ref) using $S_t = 1 - C_t$, as well as $I_t = C_t - R_t - D_t$.
Unfortunately, reported case numbers for $C_t$, $R_t$, and $D_t$, such as the daily data from JHU CSSE used in the applications in section (ref), suffer from measurement errors, see e.g.\ HorLiuSc2021. In addition, they display a strong weekly seasonal pattern that is likely to be driven by a varying number of tests conducted over the different days of the week BerSelAg2020. Under the assumption that $Y_t$ is measured with a proportionally constant error variance resulting from seasonality and measurement errors, one has the following structure for the natural logarithm of the observable $\tilde Y_t$.
Assumption (ref) specifies an unobserved components (UC) model where the observable noisy measurement $\log \tilde Y_t$ is decomposed into an unobservable measurement error $u_t$, seasonal components $\sum_{i=1}^7 \alpha_i s_{i, t}$, and the log contact rate $\log \beta_t$. The log specification accounts for a proportional impact of measurement errors and seasonality, and forces the contact rate to be strictly positive.
As the different components are not separately identified, an additional assumption on the dynamic structure of the contact rate is required. Empirical models of COVID-19 case numbers have so far assumed $\log \beta_t$ to follow a piece-wise linear time trend with structural breaks, see HarWaeWe2020, LeeLiaSe2021, LiuMooSc2021. As an alternative, the UC literature suggests to model time-varying coefficients as random walks DurKoo2012. Both specifications assume contact rate growth $\Delta \log \beta_t$ only to be contemporaneously affected either by structural breaks or by stochastic shocks, an assumption that is likely to be violated. Reflecting that the persistence properties of social behavior, and thus of the contact rate, are unknown, assumption (ref) specifies the log contact rate as a fractionally integrated process of unknown order $d$.
Under assumption (ref), the log contact rate $\log \beta_t$ is a stochastic long memory process generated by the shocks $\{\eta_i\}_{i=1}^t$. The shock $\eta_t$ models the information new in $t$, such as news reports or policy announcements. Social decisions, reflected in $\log \beta_t$, however may additionally depend on past information $\eta_{t-1},...,\eta_1$. Together, $\{\eta_i\}_{i=1}^t$ forms the information available at $t$, conditional on which social decisions, e.g.\ on whether to meet, are made. The specification takes into account that new information does not evolve deterministically, but appears as stochastic shocks, which cannot be captured by a deterministic specification as e.g.\ in LeeLiaSe2021.
The degree of persistence of the log contact rate is determined by the integration order $d$, which controls for the persistent impact of past shocks via the fractional difference operator $\Delta^d_+$. The latter exhibits a polynomial expansion in the lag operator $L$ of order infinite
The $+$-subscript denotes a truncation of an operator at $t \leq 0$, $\Delta_+^d x_t= \sum_{i=0}^{t-1}\pi_i(d) x_{t-i}$, which reflects the type II definition of fractionally integrated processes MarRob1999. For $d=1$ the log contact rate is a random walk, which follows from plugging $d=1$ into (ref). Consequently, assumption (ref) encompasses the predominant specification in the UC literature. However, assumption (ref) allows for a far more general dynamic impact of past shocks $\eta_1,...,\eta_t$ on $\log \beta_t$, as can be seen by plugging $x_t = \Delta_+^{-d}\eta_t$ into $\log \beta_t = \mu + x_t$, which gives
While a random walk is an unweighted sum of past shocks $\eta_1,...,\eta_t$, so that $\pi_i(-1)=1$ for all $i=1,...,t-1$, allowing $d \neq 1$ yields non-uniform weights of past shocks in the impulse response function of $\log \beta_t$ and thus a gradual adjustment of the log contact rate to new information. This reflects that social behavior adjusts step-wise to new information both at the individual and the aggregate level. As processing new information on the Coronavirus and revising individual decisions (e.g.\ meeting friends, traveling, working from home) takes time and evolves gradually, individuals can be expected to step-wise adjust their contacts in response to new information. Overall, individuals will react heterogeneously both in terms of speed and intensity to novel information: Some will anticipate new information faster than others, and the extent of reaction will depend on individual characteristics such as risk awareness and attitudes. Such gradual adjustments are well captured by assumption (ref), in particular when $1 < d < 2$: In that case, contact rate growth $\Delta \log \beta_t \sim I(d-1)$ is strongly persistent and mean-reverting, as will become apparent in the applications in section (ref). Strong persistence reflects the gradual adjustment of social behavior to new information, while mean-reversion ensures an asymptotically declining impact of past information to today's contact rate growth.
The remaining assumptions are imposed mainly for technical reasons. The type II definition of fractional integration assumes zero starting values for the fractionally integrated process by truncating the polynomial expansion of the fractional difference operator, $\Delta_+^d x_t= \sum_{i=0}^{t-1}\pi_i(d) x_{t-i}$. It is required to treat the asymptotically stationary ($d<1/2$, from now on `stationary' for brevity) and the asymptotically nonstationary case ($d > 1/2$, from now on `nonstationary') alongside each other. While the type II definition may be a strong assumption for some time series, it is plausible for the contact rate, as we have data covering roughly the whole pandemic. Thus, the pre-sample shocks $\eta_i$, $i \leq 0$, should be zero. Independence of $u_t$ and $\eta_t$ follows from the characterization of $u_t$ as a measurement error that should not influence the contact rate. In general, the assumption can be relaxed to allow for $\operatorname{Corr}(\eta_t, u_t) \neq 0$, as for instance in correlated UC models MorNelZi2003, and will not affect the asymptotic results in section (ref). The distributional assumptions on $\eta_t$ and $u_t$ are somewhat weaker than the assumption of Gaussian white noise on which UC models typically rely MorNelZi2003. They will be shown to be largely satisfied in the applications of section (ref). Finally, $d > 0$ is required to separately identify $\log \beta_t$ and $u_t$.
In this section, the fractional filter is derived. It is a computationally simple modification of the Kalman filter that avoids the Kalman recursions for the conditional variance. The modification is necessary, as the Kalman filter becomes computationally infeasible for UC models when the dimension of the state vector is high, as for fractionally integrated processes. The fractional filter provides a closed-form expression for the prediction error of the UC model. Based on that, a conditional-sum-of-squares (CSS) estimator for the model parameters is set up. It minimizes the sum of squared prediction errors obtained from the fractional filter. The CSS estimator is shown to be consistent and asymptotically normally distributed. Given the CSS parameter estimates, the log contact rate can be estimated by the fractional filter given the full sample information. Finally, estimation of the mean and seasonal components is considered.
Under assumptions (ref) and (ref), the fractional UC model is given by
Denote $\mu_0, \alpha_{1, 0},...,\alpha_{7, 0}, d_0, \sigma_{\eta, 0}^2, \sigma_{u, 0}^2$ as the true parameters of the data-generating mechanism. Leaving aside the deterministic terms for the moment, by defining $y_t = \log \tilde Y_t - \mu - \sum_{i=1}^7 \alpha_i s_{i, t}$, the stochastic part of the fractional UC model (ref) is
In the following, let $\theta = (d, \sigma_\eta^2, \sigma_{u}^2)' \in \varTheta$ denote the vector holding the parameters of (ref), and let $\theta_0 = (d_0, \sigma_{\eta,0}^2, \sigma_{u, 0}^2)' \in \varTheta$, where $\varTheta = D \times \varOmega_\eta \times \varOmega_u$ denotes the parameter space with $D = \{d \in \mathbb{R} | 0 < d \leq d_{max}\}$ and $\varOmega_i = \{\sigma_{i}^2 \in \mathbb{R} | 0 < \sigma_i^2 < \infty\}$, $i = \eta, u$. Define $\mathcal{F}_t$ as the $\sigma$-algebra generated by $y_1,...,y_t$, and let the expected value operator $\operatorname{E}_\theta(z_t)$ of an arbitrary random variable $z_t$ denote that expectation is taken with respect to the distribution of $z_t$ given $\theta$, so that $\operatorname{E}_{\theta_0}(z_t) = \operatorname{E}(z_t)$. Furthermore, let $\varSigma^{(i, j)}$ denote the $(i, j)$-th entry of an arbitrary matrix $\varSigma$.
Estimation of the parameters $\theta_0$ is carried out by the CSS estimator that minimizes the sum of squared prediction errors of model (ref). The prediction error is defined as the one-step ahead forecast error of $y_{t+1}$ given $\mathcal{F}_t$
It depends on $\operatorname{E}_\theta(x_{t+1} | \mathcal{F}_{t})$, for which the fractional filter provides an analytical solution. The filter is introduced in the following lemma.
The proof is contained in appendix (ref). As can be seen from lemma (ref), the fractional filter provides a solution for $\operatorname{E}_\theta(x_{t+1}| \mathcal{F}_t)$ that only depends on $\theta$ and $y_{1},...,y_t$. By plugging it into (ref), one has the closed-form expression for the prediction error
Based on (ref) the objective function of the CSS estimator for $\theta_0$ is set up
Note that estimating the parameters of the fractional UC model via the CSS estimator (ref) in combination with the fractional filter deviates from the methodological state space literature: There, expectation and variance of $x_{t+1}$ conditional on $\mathcal{F}_t$ are typically obtained from the Kalman filter recursions DurKoo2012. The resulting prediction error and its conditional variance then enter the Gaussian (quasi-)likelihood function that is maximized to estimate $\theta_0$. However, the Kalman filter becomes computationally infeasible when the dimension of the state vector is high, as for fractionally integrated processes. Thus, a computationally simpler filter is required. The fractional filter, as defined in lemma (ref), is a modification of the Kalman filter: Its solution for $\operatorname{E}_\theta(x_{t+1}| \mathcal{F}_t)$ is identical to the Kalman filter DurKoo2012, but it avoids the Kalman recursions for the conditional variance of $x_{t+1}$. While the conditional variance is necessary for (quasi-)maximum likelihood estimation, the CSS estimator only requires a closed-form expression for the prediction error, for which the fractional filter is sufficient. The objective function of the CSS estimator in (ref) is of course proportional to the exponent in the conditional Gaussian (quasi-)likelihood function. However, CSS estimation is computationally much simpler due to the fractional filter. Together, the fractional filter and the CSS estimator provide a computationally feasible alternative to the Kalman filter and the (quasi-)maximum likelihood estimator, particularly for UC models with richer long-run dynamics.
While the asymptotic theory of the CSS estimator is well established for ARFIMA models, see HuaRob2011 and Nie2015, it has not yet been derived for structural UC models. To fill this gap, theorems (ref) and (ref) summarize the asymptotic estimation theory for the CSS estimator for fractional UC models. In addition, the finite sample properties are addressed by a Monte Carlo study in appendix (ref). For consistency and asymptotic normality of the CSS estimator, the moment assumptions on the shocks $\eta_t$, $u_t$ need to be strengthened.
The proof of theorem (ref) is given in appendix (ref) and is carried out as follows: First, the model in (ref) is shown to be identified. Next, $v_t(\theta)$, as given in (ref), is shown to be integrated of order $d_0 - d$, and thus is stationary for $d_0 - d < 1/2$ and nonstationary for $d_0 - d > 1/2$. As the asymptotic behavior of the objective function changes around the point $d_0 - d = 1/2$, the objective function does not uniformly converge in probability on $\varTheta$. Adopting the results of HuaRob2011 and Nie2015, who show for ARFIMA models encompassing the reduced form of (ref) that the probability of the CSS estimator to stay in the region of the parameter space where $v_t(\theta)$ is nonstationary is asymptotically zero, the relevant region of $\varTheta$ asymptotically reduces to the region where $d_0 - d < 1/2$ holds. Within the relevant region of the parameter space this paper then proves weak convergence of the objective function by showing the objective function to satisfy a uniform weak law of large numbers. This yields consistency of the CSS estimator, see Woo1994.
The proof of theorem (ref) is again contained in appendix (ref). Since the CSS estimator is consistent, the asymptotic distribution theory is inferred from a Taylor expansion of the score function about $\theta_0$. A central limit theorem is shown to hold for the score function at $\theta_0$, together with a uniform weak law of large numbers for the Hessian matrix. The latter allows to evaluate the Hessian matrix in the Taylor expansion of the score function at $\theta_0$. Thus, the asymptotic distribution of the CSS estimator, as given in theorem (ref), can be inferred from solving the Taylor expansion for $\sqrt{n}(\hat \theta - \theta_0)$. As usual in the state space literature, no analytical solution to the asymptotic variance of the CSS estimator can be provided. The parameters of the reduced form depend non-trivially on $\theta$, so that the partial derivatives of the reduced form cannot be analytically derived. However, from theorem (ref) it follows that an estimate for the parameter covariance matrix can be obtained from the negative inverse of the Hessian matrix computed in the numerical optimization.
Estimation of the latent component $x_t$ in (ref) is considered next. In line with the methodological literature on state space models, $x_t$ is estimated by plugging the CSS estimates $\hat \theta$ into the projection
where $y_{n:1} = (y_n,...,y_1)'$, $\varSigma_{\eta_{t:1}y_{n:1}} = \operatorname{Cov}_\theta(\eta_{t:1}, y_{n:1})$, and $\varSigma_{y_{n:1}} = \operatorname{Var}_\theta(y_{n:1})$. The superscript in $\varSigma_{\eta_{t:1}y_{n:1}}^{(i, \cdot)}$ denotes the $i$-th row of the matrix, and
while the entries of $\varSigma_{y_{n:1}}$ follow from lemma (ref) by setting $t=n$. For $\theta = \theta_0$, $x_{t|n}(\theta_0)$ is the minimum variance linear unbiased estimator given $y_1,...,y_n$, see DurKoo2012. Due to theorem (ref), this property holds asymptotically for $x_{t|n}(\hat \theta)$. Note that (ref) is identical to the Kalman smoother, see DurKoo2012. However, (ref) is computationally much simpler, as it avoids the computationally intensive Kalman recursions for the conditional variance. In line with lemma (ref), (ref) is the fractional filter for $x_t$ given $\mathcal{F}_n$.
Finally, estimation of the seasonal components $\alpha_{i,0}$, $i=1,...,7$, and $\mu_0$ is considered. Theoretically, all parameters of the model in (ref) could be estimated jointly by the CSS estimator. But as TscWebWe2013b explain, including deterministic terms in the optimization can lead to poor results in finite samples for fractionally integrated processes, particularly when $d_0$ is close to unity, as the deterministic terms suffer from poor identification. They provide simulation evidence and a line of reasoning explaining why the following two-step estimator is more robust: In the first step, the integration order $d_0$ is estimated using the exact local Whittle estimator of Shi2010, which allows for unknown deterministic terms and yields $\hat d_{EW}$. Based on $\hat d_{EW}$, the deterministic terms $\mu_0, \alpha_{0,1},...,\alpha_{0,7}$ in
are estimated by ordinary least squares. In the second step, the objective function of the CSS estimator in (ref) is minimized for the adjusted $\log \tilde Y_t - \hat \mu - \sum_{i=1}^7 \hat\alpha_i s_{i, t}$.
As an alternative to (ref), one could also eliminate the seasonal components by averaging over seven neighboring observations, as $\sum_{i=1}^7 \alpha_{i, 0} s_{i, t} = 0.$ The intercept in
could then be estimated by ordinary least squares. $q$ determines whether averages are calculated solely based on past data ($q = 6$), based on centered data around $t$ ($q=3$), or based on future data ($q=0$). While the second approach does not require to estimate $\alpha_{1,0},...,\alpha_{7,0}$, averaging over seven days smooths out potential kinks in the contact rate which is problematic. Furthermore, averaging may pollute the estimates of $x_t$ and induce spurious long memory. Finally, the choice of $q$ is not trivial: While for forecasting purposes $q=6$ is adequate, choosing $q=0$ is likely to account best for the delay in reporting of case numbers, and obviously $q=3$ may be a good compromise between the two options. In the applications $\mu_0, \alpha_{0,1},...,\alpha_{0,7}$ will be estimated via (ref).
In this section, estimation results for the time-varying contact rate $\beta_t$ are presented for Canada, Germany, Italy, and the United States. The underlying data on confirmed, recovered, and deceased cases stems from the JHU CSSE. As in LeeLiaSe2021 and LiuMooSc2021, $t=1$ is set once the number of cumulative cases reaches $100$. Prior smoothing as suggested by LeeLiaSe2021 and LiuMooSc2021, who use one-sided three-day rolling averages to smooth the data, is avoided, as this likely pollutes the kinks in the contact rate that occur due to containment measures.\footnote{LiuMooSc2021 argue that one-sided three-day rolling averages smooth out noise generated by the timing of the reporting. However, the opposite should be the case, as a one-sided smoothing shifts case numbers from past to present, while a delay in reporting shifts case numbers from present to the future. To fix the latter, a forward-looking filter is required, not a backward-looking one, see the discussion at the end of section (ref).}
Instead of smoothing out seasonality, the data is adjusted for weekly seasonal patterns as described at the end of section (ref) using (ref). The bandwidth for the exact local Whittle estimator in (ref) is set to $m=\lfloor n^{0.65} \rfloor$, which is justified by the Monte Carlo study in appendix (ref). Based on the seasonally adjusted data, the parameters $\theta_0$ are estimated via the CSS estimator (ref), where $100$ combinations of starting values for $\theta_0$ are drawn from uniform distributions with appropriate support (in particular $d \in [0.5, 2]$) to avoid convergence to a local optimum. However, due to the parsimonious parametrization of the model, all combinations of starting values converged to virtually identical optima, implying that the procedure is robust to the choice of starting values. Plugging the CSS estimates into (ref), together with $\hat \mu$ in (ref), yields the log contact rate estimate $\log \hat \beta_t$.
The average infected period is required for $\mathcal{R}_t$ and is estimated by solving (ref) for $\gamma$ and taking the average
This reflects that the definition of recovered varies over the countries under study, particularly as non-hospitalized persons are typically assumed to have recovered $h$ days after they tested positive, and $h$ varies over the countries under study. The choice of $h$ proportionally affects the number of currently infected $I_t$, and thus $\beta_t$ is inversely proportional to $h$ by (ref). Consequently, only the dynamics of $\beta_t$ should be compared over the different countries, not the absolute numbers. In contrast, the reproduction rate $\mathcal{R}_t = \beta_t / \gamma$ accounts for the different $h$ when $\gamma$ is estimated via (ref). If instead $\gamma = 1/18$ is fixed as in LeeLiaSe2021, the dependence on $h$ is not resolved and countries with a higher $h$ will exhibit a smaller reproduction rate by construction. This is precisely the reason for the implausible estimates for $\mathcal{R}_t$ in LeeLiaSe2021, and is solved by accounting for different $h$ via (ref).
Results are reported for Canada, Germany, and Italy in subsection (ref). They are selected as they are all members of the G7 and have implemented containment measures of different strength, duration, and at different points in time, thus making a comparison interesting. For the selected countries, there exist reliable data on confirmed, recovered, and deceased cases provided by the JHU CSSE. As will be shown, the latter is not the case for the US, where data on recovered subjects suffers heavily from under-reporting, yielding a severe downward-bias for the estimated contact rate and the resulting reproduction rate $\mathcal{R}_{t}$ as reported by LeeLiaSe2021. The problem is fixed by an assumption on the average duration of an infection, and results for the US are presented in subsection (ref). As will become apparent there, the fractional filter is quite robust to under-reporting of recovered cases. To monitor the pandemic in real-time, subsection (ref) examines the precision of the fractional filter at the end of the sample.
For Canada, figure (ref) sketches the estimated log contact rate and the resulting reproduction rate $\hat{\mathcal{R}}_{t} = \hat \beta_t/\hat\gamma$ in the first row. The average duration of an infection is estimated to be $1/\hat \gamma = 18.29$ days. The second row of figure (ref) displays the estimated prediction error $v_t(\hat \theta)$ and its estimated autocorrelation function.
Based on the top-left panel of figure (ref), several turning points of the contact rate can be identified using a simple algorithm that defines a minimum (maximum) whenever the contact rate $\beta_t$ at $t$ is smaller (greater) than all $\beta_{t+1},...,\beta_{t+10}$, the contact rates of the next ten days. These periods correspond to several policy regimes characterized by the strengthening and easing of containment measures. While a small selection of policy measures is presented below, a detailed overview is given by McCSmiAn2020.
The estimated contact rate and the resulting reproduction rate are well in line with the chronology of policy interventions. In particular, the fractional filter allows to identify turning points of the contact rate that are not visible from the raw data that is plotted in gray color in figure (ref).
The two graphs at the bottom of figure (ref) illustrate how well the Canadian data fits the model assumptions. Assumptions (ref) and (ref) assume the measurement error $u_t$ and the log contact rate shock $\eta_t$ to be homoscedastic white noise processes. Since $\hat \theta$ is consistent, see theorem (ref), by (ref) the prediction error $v_t(\hat \theta)$ becomes a white noise process as $n \to \infty$ if assumptions (ref) and (ref) hold. While some outliers exist, about $95\%$ of the prediction errors lie within two standard deviations, as the bottom-left panel illustrates. The bottom-right panel shows that there is not much autocorrelation left in the prediction error. Given the parsimonious parametrization of the fractional UC model, this is surprising. The two panels at the bottom of figure (ref) thus substantiate that the dynamics of the log contact rate are well captured by a fractionally integrated process.
The estimated integration order is $\hat d = 1.2166$, which implies that a unit shock on the contact rate growth $\Delta \log \beta_t$ retains $21.66\%$ of its impact in $t+1$, $13.17\%$ in $t+2$, and $9.73\%$ in $t+3$. After one week, the impact is still $5.10\%$, after two weeks $2.98\%$, and after three weeks $2.17\%$, which is due to the strong persistence of fractionally integrated processes, see assumption (ref) for the formula for $\pi_i(d-1)$. The slow decay may very well describe the persistent impact of past information shocks on today's social behavior.
For Germany, figure (ref) plots the empirical results. The average infected period is estimated to be $1/\hat \gamma = 21.27$ days, which is slightly greater compared to Canada and likely results from different algorithms to estimate the number of recovered individuals. The integration order estimate is of similar size as for Canada ($\hat d = 1.2693$, see table (ref)).
From the top-left panel of figure (ref) the following contact regimes can be identified:
Similar to Canada, the two panels at the bottom of figure (ref) indicate that the model assumptions are largely satisfied, despite little remaining autocorrelation in the prediction error.
For Italy, a slight adjustment of the data is required, as $\Delta C_t = 0$ on June 19 and thus $\log \tilde Y_t$ is not defined, see (ref). To adjust the single observation, averages from the neighboring observations are used (i.e.\ $\Delta C_t = 1/3 (\Delta C_{t-1} + \Delta C_{t+1})$), while the adjusted cases in $t-1$ and $t+1$ will equal $2/3$ of the reported cases. Thus, the cumulated number of reported cases is unaffected by the adjustment.
The empirical results are displayed in figure (ref). An estimate for the average infected period is $1/\hat \gamma = 35.92$ days, which is significantly higher than for Canada and Germany. The discussion below (ref) gives an explanation for the high variation in $\gamma$ over the different countries. The estimated integration order is $\hat d =1.4304$ (see table (ref)), which is somewhat greater than the estimates for Germany and Canada, implying that a shock $\eta_t$ in Italy will yield a more persistent effect on the contact rate. This may be explained by the severity of the pandemic in Italy in spring 2020, which likely had a long-lasting impact on social behavior. Based on the top-left panel of figure (ref) the following regimes and turning points are visible:
The two panels at the bottom of figure (ref) are similar to Canada and Germany, and indicate that the prediction errors are rather homoscedastic, although outliers exist, and little autocorrelation is left.
The US is treated separately, since data on recovered cases reported by the JHU CSSE seem heavily downward-biased. To see this, consider the difference between lagged cumulative confirmed, cumulative recovered, and cumulative deceased cases for different lags $h$
For $h=0$, (ref) measures the number of currently infected subjects. For small $h$, (ref) should be positive, as it takes some time for the infected subjects to either recover or die. As $h$ increases, (ref) should turn negative, as an increasing number of subjects contained in the cumulative cases $C_{t-h}$ and subjects infected between $t-h$ and $t$ (and thus contained in $C_t - C_{t-h}$) either recover or die. The turning point, denoted by $\bar h$, should be close to the average infected period $1/\gamma$, as long as new confirmed cases between $t-h$ and $t$, i.e.\ $C_t - C_{t-h}$, do not explode. If they do, then $\bar h$ should be smaller than $1/\gamma$, as outflows from $C_t - C_{t-h}$ disproportionally increase $R_t$ and $D_t$.
Figure (ref) plots (ref) in case numbers for lags $h=15,18,...,42,45$. As can be seen, even after $45$ days the difference between lagged cumulative confirmed, cumulative recovered, and cumulative deceased cases is predominantly positive. This is at odds with the average infected periods for Canada, Germany, and Italy as found in subsection (ref), and indicates that data on $R_t$, $D_t$ may suffer from under-reporting. As stated by the JHU CSSE, data on recovered cases are based on local media reports as well as on state reporting when available. US state-level recovered cases stem from the COVID Tracking Project (www.covidtracking.com). CTR2020 recently pointed out that several states and territories, including California and Florida, do not report data on recovered cases, which may explain the downward-bias. In addition, definitions of recovered cases differ considerably across states, and the majority of the states consider a case as recovered if a certain number of days (generally between 10 and 30) after a positive test result or symptom onset have passed and the patient has not died.
Since no reliable data on recovered cases is available for the US, an approximation is required. In the following, it will be assumed that individuals either recover or die $\bar{h} = 21$ days after they tested positive, $C_{t-21} = D_t + R_t$. The assumption is justified as follows: First, it is similar to the average infected period estimated for Germany and more conservative than the estimate for Canada. And second, it is centered in the range of definitions for recovered individuals by the federal states. In addition, estimates for $\bar{h} = 18$ and $\bar h = 24$ days are presented, which gives a reasonable interval for the contact rate.
Under the assumption of $\bar h = 21$, the estimated integration order equals $\hat d = 1.2499$ (see table (ref)) and is very similar to the results in subsection (ref). Estimates for contact and reproduction rate are visualized in figure (ref) and again allow to decompose the chronology of the pandemic into different regimes:
As can be seen from figure (ref), estimates for the contact rate are rather robust to the choice of $\bar h$. They are slightly greater for $\bar h = 18$, as the number of currently infected $I_t$ is smaller by construction and thus additional contacts are required to explain new confirmed cases, while they are slightly smaller for $\bar h = 24$ exactly for the opposite reason. The estimated reproduction rate is virtually identical among the three scenarios, as it is normalized by the average infected period $\hat{\mathcal{R}}_{ t} = \hat \beta_t / \hat \gamma$. For $\bar h=21$, the two panels at the bottom of figure (ref) indicate that the model assumptions are largely satisfied, despite some weak correlation in the prediction errors. The plots are very similar for $\bar h=18$ and $\bar h=24$ and thus not shown. The reproduction rate $\hat{\mathcal{R}}_{t}$ is greater than unity during the whole sample, which contradicts the results of LeeLiaSe2021 who rely on the downward-biased data on recovered subjects from the JHU CSSE while fixing the average infected period to be $1/\gamma = 18$.
This subsection investigates the end-of-sample properties of the fractional filter for real-time estimation of the contact rate. Reliable contact rate estimates at the current frontier of the data would allow to real-time monitor the state of the pandemic and can serve as a surveillance measure for future outbreaks. Based on reliable real-time estimates for the contact rate, policy rules can be implemented to prevent an exponential growth of case numbers. Acting early reduces economic and social costs of containment measures, and consequently a well-designed policy rule will be beneficial, given that the fractional filter yields a reliable estimate for the current level of the contact rate. Drawing inference on the latter is the focus of this subsection.
In detail, real-time monitoring is simulated by truncating the sample at a certain point $t$, $r \leq t \leq n$, where $r$ is the minimum sample size for the CSS estimator to produce reasonable estimates. The parameters $\theta_0, \mu_0, \alpha_{1,0},...,\alpha_{7, 0}$ of (ref) are then estimated as described in section (ref) using the information available at time $t$, $\mathcal{F}_t$, and the resulting parameter estimates are denoted as $\hat \theta^{(t)}$, $\hat \mu^{(t)}$, etc. To take into account reporting lags, and to be robust against outliers at the end of the sample, a little backward-smoothing is allowed by reporting the smoothed estimate for the log contact rate at period $t-3$ given the information available at period $t$. From (ref), the smoothed estimates are
As (ref) only depends on information available at $t$, it mimics the situation of a policy maker at $t$ and can be used to draw inference on the monitoring properties of the fractional filter at time $t$. Based on $\hat \beta_{t-3|t}$, policy rules to prevent an exponential spread of the virus can be designed. Such rules could, for instance, define a threshold for $\hat{\mathcal{R}}_{t-3}$ at which additional containment measures are implemented. As the threshold should naturally depend on the number of currently infected, current hospital capacities, and other parameters, the precise design of such a policy rule is left to the experts, and only a primitive policy rule will be introduced later for illustrative purposes.
The reliability of real-time estimates for the contact rate is assessed by the following experiment: First, an estimation sample that consists of information available until May 31 is defined, for which $\theta_0, \mu_0, \alpha_{1,0},...,\alpha_{7, 0}$ are estimated. It consists of at least 80 observations, which is considered as a reasonable sample size for the estimation sample. Based on these estimates, $\log \hat \beta_{r-3|r}$ is obtained as described above. In a second step, information available on June 1 is added to the sample and parameter estimates are updated using the $\hat \theta^{(r)}$ from the estimation sample as starting values for the CSS estimator, which gives $\hat \theta^{(r+1)}$. As before, the estimate for $\log \hat \beta_{r-2|r+1}$ is stored. The procedure repeats for all $t$, $r < t \leq n$, where in every step $t$ the CSS estimator is initialized by $\hat \theta^{(t-1)}$. The resulting real-time estimates for the contact rate are then compared to those of subsections (ref) and (ref) to draw inference on their reliability. In addition, a primitive policy rule is introduced. It assumes governments to take action as soon as $\hat{\mathcal{R}}_{t-3} = \hat{\beta}_{t-3|t} / \hat \gamma > 1.2$. The latter is motivated by the observation that preventing an exponential propagation (i.e.\ $\mathcal{R}_{ t-3} > 1$) is desirable, and a margin of $0.2$ is included to be robust against outliers. Finally, the real-time contact rate estimates are compared to a rolling seven-day average
which includes three forward-looking observations and should smooth out the seasonality.
The real-time experiment considered in this paper deviates from LeeLiaSe2021, who suggest to monitor the current state of the pandemic by fitting a linear time trend with structural breaks to $\log \tilde Y_t$. LeeLiaSe2021 evaluate the monitoring properties of their contact rate estimate ex-post, using all information available in their sample. Consequently, their estimates at point $t$ depend on information that was not available to policy makers at period $t$ whenever $t < n$. As structural breaks are not well identified at the end of the sample, recent changes in the contact rate cannot be expected to be found by the estimator of LeeLiaSe2021.
The results of the real-time experiment are visualized in figure (ref). The four panels on the left side sketch the resulting real-time estimates for the contact rate for the four countries under study, together with the (full sample) results of subsections (ref) and (ref). As can be seen, the real-time estimates almost perfectly overlap with estimates using the full sample information. This implies that the real-time estimates are well suited for monitoring the current state of the contact rate. Minor deviations are visible for Canada at the end of July, and for Germany at the middle of June, while no greater deviations are visible for Italy and the US. The real-time estimates exceed the threshold $\hat{\mathcal{R}}_{t-3} = \hat{\beta}_{t-3|t} / \hat \gamma > 1.2$ on the same day as estimates based on the full sample for Italy (August 13) and for the US (May 31), where May 31 is the first observation for the real-time estimates. For Canada, the reproduction rate based on the full sample exceeds the threshold on July 21, one day after the real-time estimate. For Germany, the real-time estimates exceed the threshold for the reproduction rate on June 19, two days before those based on the full sample. These findings again substantiate the usefulness of the fractional filter as a surveillance measure to monitor the current state of the pandemic. From the results in figure (ref), it follows that the primitive policy rule would have required governments to take action during the summer where case numbers were comparably low across the four countries under study. Such an early intervention would have likely reduced the economic and social costs compared to the lockdown measures as implemented at the end of year 2020.
The four panels on the right side of figure (ref) display the deviations of $\log \hat \beta_{t-3|t}$ and $\log \hat \beta_{t-3|t}^{benchmark}$ from the log contact rate estimates based on the full sample information $\log \hat \beta_{t-3|n}$. Thus, they shed light on whether the fractional filter improves estimates for the contact rate compared to a rolling seven-day average that uses three forward-looking observations. For Italy and the US, the advantages of the fractional filter directly become apparent, as the benchmark exhibits greater deviations. For Canada and Germany, the fractional filter performs comparably well when large outliers occur, e.g.\ around July 20 for Canada and around June 20 for Germany.
To extract a time-varying signal for the COVID-19 contact rate from daily data on confirmed, recovered, and deceased cases, this paper introduces a novel unobserved components model. It models the log contact rate as a fractionally integrated process of unknown integration order. A computationally simple modification of the Kalman filter is introduced and is termed the fractional filter. It provides a closed-form expression for the prediction error that allows to estimate the model parameters by a conditional-sum-of-squares (CSS) estimator. The asymptotic theory for the CSS estimator is provided. For the countries under study, estimation results are well in line with the chronology of the pandemic. They allow to draw inference on the impact of policy measures such as contact restrictions. The new filtering method bears great potential as a monitoring device for the current state of the pandemic, as it yields reliable contact rate estimates at the current frontier of the data.
As vaccines become more and more available, future research can generalize the model to include the number of vaccinated. For instance, this can be done by decomposing $1 = S_t + I_t + R_t + D_t + V_t$, where $V_t$ is the fraction of vaccinated. The states $R_t$ and $V_t$ should be non-overlapping as long as vaccines are not rolled out to recovered subjects. While vaccine recommendations vary over the different countries, some assign a lower priority to recovered subjects, so that $R_t$ and $V_t$ are non-overlapping at the early stage of the vaccine roll-out. Furthermore, mutations of the Coronavirus can be taken into account e.g.\ by allowing for a smooth transition between a contact rate with a low probability of virus transmission and one with a high probability.
For applications beyond COVID-19 related data, the fractional filter offers a robust, flexible, and data-driven way for signal extraction of data of unknown persistence. It requires no prior assumptions on the integration order of a process, and thus provides a solution to model specification in the unobserved components literature. Due to its computational advantages compared to the classic Kalman filter, it allows to estimate unobserved components models with richer dynamics.
The author thanks Nicolas Apfel, Uwe Hassler, Timon Hellwagner, Roland Jucknewitz, Alina Prechtl, Veronika P\"uschel, Lars Schlereth, Rolf Tschernig, Enzo Weber, and the participants of the Department Seminar at the University of Regensburg for very helpful comments.