key: cord-280683-5572l6bo
authors: Liu, Laura; Moon, Hyungsik Roger; Schorfheide, Frank
title: Panel forecasts of country-level Covid-19 infections()
date: 2020-10-16
journal: J Econom
DOI: 10.1016/j.jeconom.2020.08.010
sha: 
doc_id: 280683
cord_uid: 5572l6bo

We use a dynamic panel data model to generate density forecasts for daily active Covid-19 infections for a panel of countries/regions. Our specification that assumes the growth rate of active infections can be represented by autoregressive fluctuations around a downward sloping deterministic trend function with a break. Our fully Bayesian approach allows us to flexibly estimate the cross-sectional distribution of slopes and then implicitly use this distribution as prior to construct Bayes forecasts for the individual time series. We find some evidence that information from locations with an early outbreak can sharpen forecast accuracy for late locations. There is generally a lot of uncertainty about the evolution of active infection, due to parameter and shock uncertainty, in particular before and around the peak of the infection path. Over a one-week horizon, the empirical coverage frequency of our interval forecasts is close to the nominal credible level. Weekly forecasts from our model are published at https://laurayuliu.com/covid19-panel-forecast/.

This paper contributes to the rapidly growing literature on generating forecasts related to the current Covid-19 pandemic. We are adapting forecasting techniques for panel data that we have recently developed for economic applications such as the prediction of bank profits, charge-off rates, and the growth (in terms of employment) of young firms; see Liu (2020) , Liu, Moon, and Schorfheide (2020) , and Liu, Moon, and Schorfheide (2019) . We focus on the prediction of the smoothed daily number of active Covid-19 infections for a crosssection of approximately one hundred countries/regions, henceforth locations. The data are obtained from the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University. While we are currently focusing on country-level aggregates, our model could be easily modified to accommodate, say, state-or county-level data.

In economics, researchers distinguish, broadly speaking, between reduced-form and structural models. A reduced-form model summarizes spatial and temporal correlation structures among economic variables and can be used for predictive purposes assuming that the behavior of economic agents and policy makers over the prediction period is similar to the behavior during the estimation period. A structural model, on the other hand, attempts to identify causal relationships or parameters that characterize policy-invariant preferences of economic agents and production technologies. Structural economic models can be used to assess the effects of counterfactual policies during the estimation period or over the out-of-sample forecasting horizon.

The panel data model developed in this paper to generate forecasts of Covid-19 infections is a reduced-form model. It processes cross-sectional and time-series information about past infection levels and maps them into predictions of future infections. While the model specification is motivated by the time-path of infections generated by the workhorse compartmental model in the epidemiology literature, the so-called susceptible-infected-recovered (SIR) model, it is not designed to answer quantitative policy questions, e.g., about the impact of social-distancing measures on the path of future infection rates.

Building on a long tradition of econometric modeling dating back to Haavelmo (1944) , our model is probabilistic. The growth rates of the infections are decomposed into a deter-about model parameters and uncertainty about future shocks. We model the growth rate of active infections as autoregressive fluctuations around a deterministic trend function that is piecewise linear. The coefficients of this deterministic trend function are allowed to be heterogeneous across locations. The goal is not curve fitting -our model is distinctly less flexible in samples than some other models -but rather out-of-sample forecasts, which is why we prefer to project growth rates based on autoregressive fluctuations around a parsimonious linear time trend with a single break.

A key feature of the Covid-19 pandemic is that the outbreaks did not take place simultaneously in all locations. Thus, we can potentially learn from the speed of the spread of the disease and subsequent containment in country A, to make forecasts of what is likely to happen in country B, while simultaneously allowing for some heterogeneity across locations. In a panel data setting, one captures cross-sectional heterogeneity in the data with unit-specific parameters. The more precisely these heterogeneous coefficients are estimated, the more accurate are the forecasts. A natural way of disciplining the model is to assume that the heterogeneous coefficients are "drawn" from a common probability distribution. If this distribution has a large variance, then there is a lot of country-level heterogeneity in the evolution of Covid-19 infections. If instead, the distribution has a small variance, then the path of infections will be very similar across samples, and we can learn a lot from, say, China, that is relevant for predicting the path of the disease in South Korea or Germany. Formally, the cross-sectional distribution of coefficients can be used as a so-called a priori distribution (prior) when making inference about country-specific coefficients. Using Bayesian inference, we combine the prior distribution with the unit-specific likelihood functions to compute a posteriori (posterior) distributions. This posterior distribution can then be used to generate density forecasts of future infections. Unfortunately, the cross-sectional distribution of heterogeneous coefficients is unknown. The key insight in the literature on Bayesian estimation of panel data models is that this distribution, which is called random effects (RE) distribution in the panel data model literature, can be extracted through simultaneous estimation from the cross-sectional dimension of the panel data set. There are several ways of implementing this basic idea.

reflect parameter uncertainty as well as uncertainty about shocks that capture deviations from the deterministic component of our forecasting model.

Our empirical analysis makes the following contributions. First, we present estimates of the RE distribution as well as the distribution of location-specific coefficient estimates.

Second, we document how density forecasts from our model have evolved over time, focusing on the forecasts for three countries in which the level of infections peaked at different points in time: South Korea, Germany, and the U.S. Due to the exponential transformation from growth rates to levels, density forecasts can feature substantial tail risk by assigning nontrivial probability to very high infection levels which materialized in the U.S. but not in Germany and South Korea.

Third, we evaluate one-week and four-week ahead density forecasts based on the continuous ranked probability score and interval forecasts based on cross-sectional coverage frequency and average length. In addition to forecasts from our panel data model, we also consider forecasts based on location-level time series estimates of our trend-break model and a simple SIR model.

Once we decompose the set of locations into those that experienced the Covid-19 outbreak early (prior to 2020-03-28) and those that experience the outbreak later on, then we find some evidence that for the late group the panel density forecasts are more accurate than the time-series forecasts. However, because of the substantial heterogeneity in our panel and the poor data quality for some countries, the empirical evidence in favor of the panel approach is not as tidy as the simulation evidenced provided in the Monte Carlo section of this paper.

Over time, in particular after the infection level has peaked and started to fall, forecast accuracy increases. The timing of the peak appears to be very difficult to forecast. Prior to the middle of May the panel and time-series forecasts from our trend-break model are substantially more accurate than the forecasts from a simple time-varying coefficient SIR model. For subsequent forecast origins, the accuracy across the three forecasting procedures becomes much more similar. Weekly real-time forecasts are published on the companion website https://laurayuliu.com/covid19-panel-forecast/.

In terms of interval forecasts we find that over a one-week horizon the empirical coverage frequency of the trend-break model forecasts is close to the nominal coverage level based on which the forecasts were constructed. Moreover, in April and May, the average interval lengths of the panel model forecasts are slightly smaller than the time-series intervals. At the four-week horizon the coverage frequency is considerably smaller than the nominal level and it deteriorates further for longer horizons. This paper is connected to several strands of the literature. The panel data forecasting approach is closely related to work by Gu and Koenker (2017a,b) and our own work in Liu (2020), Liu, Moon, and Schorfheide (2020) , Liu, Moon, and Schorfheide (2019) . All five papers focus on the estimation of the heterogeneous coefficients in panel data models. The forecasting model for the Covid-19 infections is based on the alternative parametric model considered in Liu (2020) and tailored to the specifics of the Covid-19 pandemic. The approach has several desirable theoretical properties. For instance, Liu, Moon, and Schorfheide (2020), building on Brown and Greenshtein (2009) , show that an empirical Bayes implementation of the forecasting approach based on Tweedie's formula can asymptotically (as the crosssectional dimension tends to infinity) lead to forecasts that are as accurate as the so-called oracle forecasts. Here the oracle forecast is an infeasible benchmark that assumes that the distribution of the heterogeneous coefficients is known to the forecaster. Liu (2020) shows that the density forecast obtained from the full Bayesian analysis converges strongly to the oracle's density forecast as the cross-section gets large.

The piecewise linear conditional mean function for the infection growth rate resembles a spline; see de Boor (1990) for an introduction to spline approximation. Unlike a typical spline approximation in which the knot locations are free parameters and some continuity of smoothness restrictions are imposed, the knot placement in our setting is closely tied to the first component of the spline, and we do not impose continuity. Our model could be generalized by adding additional knots in the deterministic trend component of infection growth rates, but the extension is not pursued in this paper. Other authors have explored alternative forms of nonlinearity which are often tied to the object that is being modeled, e.g., active infections, cumulative infections, new infections, deaths. For instance, Li and Linton (2020) model the logarithm of country-level new infections and new deaths via a quadratic trend, using rolling samples. Ho, Lubik, and Matthes (2020) model the cumulative number of infections using a very flexible nonlinear parametric function.

An important aspect of our modeling framework is that the panel model is specified in event time, i.e., time since the level of infections in a particular location exceeds 100.

The forecasts, however, are generated based on calendar time. This allows us to sharpen forecasts for countries/regions that experienced an outbreak at a late stage (in terms of calendar time), based on information from locations with an early outbreak. This idea J o u r n a l P r e -p r o o f Journal Pre-proof is also utilized by Larson and Sinclair (2020) who use state-level panel data to nowcast unemployment insurance claims during Covid-19.

A growing number of researchers with backgrounds in epidemiology, biostatistics, machine learning, economics, and econometrics are engaged in forecasting aspects of the Covid-19 pandemic. Because this is a rapidly expanding and diverse field, we do not attempt to provide a meaningful survey at this moment. Instead, we simply provide a few pointers.

Forecasts are reported in the abovementioned papers by Li and Linton (2020) and Ho, Lubik, and Matthes (2020) . The paper by Avery, Bossert, Clark, Ellison, and Fisher Ellison (2020) The remainder of this paper is organized as follows. Section 2 provides a brief survey of epidemiological models with a particular emphasis on the SIR model. The specification of our panel data model is presented in Section 3. Section 4 contains a small-scale Monte Carlo study and the empirical analysis is conducted in Section 5. Finally, Section 6 concludes.

There is a long history of modeling epidemics. A recent survey of modeling approaches is provided by Bertozzi, Franco, Mohler, Short, and Sledge (2020) . The authors distinguish three types of macroscopic models: 1 (i) the exponential growth model; (ii) self-exciting point processes / branching processes; (iii) compartmental models, most notably the SIR model that divides a population into susceptible (S t ), infected (I t ), and resistant (R t ) individuals. Our subsequent discussion will focus on the exponential growth model and the SIR model. While epidemiological models are often specified in continuous time, we will con-sider a discrete-time specification in this paper because it is more convenient for econometric inference.

The exponential model takes the form I t = I 0 exp(γ 0 t). The number of infected individuals will grow exponentially at the constant rate γ 0 . This is a reasonable assumption to describe the outbreak of a disease, but not the subsequent dynamics because the growth rate will typically fall over time and eventually turn negative as more and more people become resistant to the disease. The SIR model dates back to Kermack and McKendrick (1927) . In its most elementary version it can be written in discrete-time as follows:

where N is the (fixed) size of the population, β is the average number of contacts per person per time, and γ is the rate of recovery or mortality. The model could be made stochastic by assuming that β and γ vary over time, e.g.,

In response to the recent Covid-19 pandemic, several introductory treatments of SIR models have been written for economists, e.g., Avery, Bossert, Clark, Ellison, and Fisher Ellison (2020) and Stock (2020) . Moreover, there is a growing literature that combines compartmental models with economic components. In these models, economic agents account for the possibility of contracting a disease when making their decisions about market participation.

This creates a link between infection rates and economic activity through the frequency of interactions. Examples of this work in macroeconomics include Eichenbaum, Rebelo, and Trabandt (2020) , Glover, Heathcote, Krueger, and Rios-Rull (2020) , and Krueger, Uhlig, and Xie (2020). The advantage of models that link health status to economic activity is that they can be used to assess the economic impact of, say, social distancing measures.

We now simulate the constant-coefficient SIR model in ( model. It is a monotonically decreasing function of time that we approximate by fitting a piecewise linear least-squares regression line with a break point at t * which is the point in time when the infections peak and the growth rate transitions from being positive to being negative. Under the second parameterization the transmission rate β = 0.06 is much lower and the recovery rate is slightly faster. This leads to an almost bell-curve shaped path of infections. While the resulting growth rate of the infections is not exactly a linear function of time t, the break at t * is much less pronounced. While the piecewise-linear regression functions do not fit perfectly, they capture the general time-dependence of the growth-rate path implied by the SIR model. In particular, they allow for a potentially much slower change in the growth rate of infections after the peak.

We use these simulations as a motivation for the subsequent specification of our empirical model. 2 This model assumes that the growth rate of infections is a decreasing piecewiselinear function of time with a break when the growth rates cross zero and the infections peak. This deterministic component is augmented by a stochastic component that follows a first-order autoregressive, AR(1), process. We refer to the model as trend-break model. We will revisit a stochastic version of the SIR model that comprises (1) and (2) in Section 5.4

where we compare its forecasts to the proposed trend-break model.

We now describe our empirical model in more detail. We begin with the specification of a regression model for the growth rate of infections in Section 3.1. Our model features location-specific regression coefficients and heteroskedasticity. The prior distribution for the Bayesian analysis is summarized in Section 3.2. Posterior inference is implemented through a Gibbs sampler that is outlined in Section 3.3. Further computational details are provided through replication files on the companion webpage. The algorithm to obtain simulated infection paths from the posterior predictive distribution is outlined in Section 3.4.

We specify a panel data model for infection growth rates y it = Δ ln I it , i = 1, . . . , N and t = 1, . . . , T . We assume that

where

] captures the size of the break in the regression coefficients at t = t * i . The deterministic part of y it corresponds to the piecewise-linear regression functions fitted to the infection growth paths simulated from the SIR in Figure 1 .

The serially-correlated process u it generates stochastic deviations from the deterministic path γ i x t of the infection growth rate. The u it shocks may capture time variation in the (β, γ) parameters of the SIR model or, alternatively, model misspecification. In Section 2 the break point t * i was given by the peak of the infection path. Abstracting from a potential discontinuity at the kink, we define t * i as

which implies that E[y it |t = t * i ] = 0. Because of the AR(1) process u it , t * i is not the peak of the observed sample path, nor is it an unbiased or consistent estimate of the period in which the infections peak. For δ i = 0, the model reduces to

Note that the break date t * i is identified in this model even if δ i = 0, because we assume the break occurs when the deterministic component of the growth rate falls below zero.

To construct a likelihood function we define the quasi-difference operator Δ ρ = 1 − ρL such that Δ ρ u it = it . Thus, we can rewrite (3) as follows

Now let λ i = [γ i , δ i ] and n λ be the dimension of λ. The parameters of the panel data model are (ρ, λ 1:N , σ 2 1:N ). Here, we use the notation Z 1:L to denote the sequence z 1 , . . . , z L . Using this notation, we denote the panel observations by Y 1:N,1:T . We will subsequently condition J o u r n a l P r e -p r o o f Journal Pre-proof on Y 1:N,0 to initialize conditional likelihood function. Finally, from the growth-rates y it we can easily recover the level of active infections as

To conduct Bayesian inference, we need to specify a prior distribution for ( ρ, λ 1:N , σ 2 1:N ). We do so conditional on a vector of hyperparameters ξ that do not enter the likelihood function.

Our prior distribution has the following factorization:

where ∝ denotes proportionality and f (•) is an indicator function that we will use to impose the following sign restrictions on the elements of λ i :

The restriction γ 1i < 0 ensures that the growth rates are falling over time. After the break point the rate of decline decreases (δ 1i > 0), but stays negative (γ 1i + δ 1i < 0). In addition we assume that the decrease in the rate of decline is associated with a downward shift, i.e., δ 0i < 0, of the intercept as shown in the SIR simulation.

Because of the presence of the indicator function f (•) the right-hand side of (8) is not a properly normalized density. In view of the indicator function f (•) we define the RE distribution of λ i given ξ as

In turn, the marginal prior distribution of the hyperparameters is given by

Building on Liu (2020), we use the following densities p(•) in (8) for ρ, λ i , and σ 2 i :

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

Thus, the vector of hyperparameters is ξ = (μ, Σ, a, b). We decompose p(ξ) = p(μ, Σ)p(a, b).

The density p(μ, Σ) is constructed as follows:

The degrees of freedom for the Inverse Wishart distribution is set to

The matrix W 0 is constructed to align the scale of the variance of μ i with the cross-sectional variance of the data, adjusting for the average magnitudes of the regressors that multiply the λ i elements.

To obtain the density p(a, b), we follow Llera and Beckmann (2016) 

The parameters (α a , β a , γ a , α b , β b ) need to be chosen by the researcher. We use α a = 1, β a = γ a = α b = β b = 0.01, which specifies relatively uninformative priors for hyperparameters a and b.

Posterior inference is based on an application of Bayes Theorem. Let p(Y 1:N,1:T |λ 1:N , σ 2 1:N , ρ) denote the likelihood function (for notational convenience we dropped Y 1:N,0 from the conditioning set). Then the posterior density is proportional to p(ρ, λ 1:N , σ 2 1:N , ξ|Y 1:N,0:T ) ∝ p(Y 1:N,1:T |λ 1:N , σ 2 1:N , ρ)p(ρ)p λ 1:N , σ 2 1:N , ξ ,

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

where the prior was given in (8). To generate draws from the posterior distribution we use a Gibbs sampler that iterates over the conditional posterior distributions

The Gibbs sampler generates a sequence of draws ρ s , λ s 1:N , (σ 2 1:N ) s , ξ s , s = 1, . . . , N sim , from the posterior distribution. The implementation of the Gibbs sampler closely follows Liu (2020) .

For the Gibbs sampler to be efficient, it is desirable to have a model specification in which it is possible to directly sample from the conditional posterior distributions in ( 15). Unfortunately, the exact likelihood function leads to a non-standard conditional posterior distribution for λ 1:N |(Y 1:N,0:T , ρ, σ 2 1:N , ξ) because γ i enters the indicator function in (3) through the definition of t * i . Thus, rather than using the exact likelihood function, we will use a limited-information likelihood function of the form

The densities p l (Y i,1:T |λ i , σ 2 i ) are constructed as follows. Let Δ be some positive number, e.g., three or five time periods. Given a sample (Y i,1:T , ln I i,1:T ) we define

On the other hand, if t i,max < T , then it is likely that t * = t i,max . Thus, we distinguish two cases:

Because δ i does not enter the likelihood function, its posterior is p(δ

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

Now δ i does enter the likelihood function and its prior gets updated in view of the data.

Bayesian forecasts reflect parameter and shock uncertainty. We simulate trajectories of infection growth rates from the posterior predictive distribution using the Algorithm 1. The simulated growth rates can be converted into simulated trajectories for active infections using (7).

Algorithm 1 (Simulating from the Posterior Predictive Distribution)

2. Based on the simulated paths I s 1:N,T +1:T +H , s = 1, . . . , N sim , compute point, interval, and density forecasts for each period t = T + 1, . . . , T + H.

We now conduct a small Monte Carlo experiment that compares the forecasts derived from the panel data model to time-series forecasts generated for each location separately. The experiment shows that in our environment forecasts for locations that experience an outbreak at a later point in time are more accurate than forecasts for locations that have an early outbreak because the early outbreaks facilitate learning about the RE distribution that benefits the forecasts for the remaining locations. The data generating process (DGP) is described in Section 4.1 and the results are summarized in Section 4.2.

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

Initial infection level: I i0 = 101 for all i.

The data generating process (DGP) is given by the trend-break model (3) for the growth rates of infections. For the simulation experiment we assume that the innovations it are homoskedastic, i.e., σ 2 i = σ 2 for all i. The DGP matches certain aspects of the empirical application in Section 5, but it is more stylized in other dimensions. The time period t is a day. The number of locations in our simulation is N = 150. We split the locations into two groups: N 1 = 75 locations experience an early outbreak, starting at t = 1, and N 2 = 75 locations experience a late outbreak, starting at t = 56. We refer to these groups as "early" and "late." For the early group calendar time and event time are identical. For the late group, the event time is calendar time minus T δ = 56 (8 weeks).

The parameters of the DGP are summarized in Table 1 . The persistence of the growth rates is set to ρ = 0.8. The dispersion of the parameters λ i is controlled by a vector of means, λ and a covariance matrix Σ. Both are calibrated to match some stylized facts about the cross-sectional distribution of the country-level data used in Section 5. We then draw the λ i s independently from the N (λ, Σ) distribution. The innovation variance σ 2 corresponds to a high-density value of the estimated density σ 2 i ∼ IG(a, b). We assume that the outbreak starts in each geographical location i with I i0 = 101. There is also heterogeneity in the timing of the peak, which is illustrated in Figure 2 . The figure shows the percentage of locations that have peaked in or prior to period t. By construction, infections in the early-group locations tend to peak sooner than in the late-group locations. However, the peak dates in each group are quite dispersed: only 20% of early locations have peaked after 60 days. It takes more than 100 days for the remaining early locations to peak.

Forecasting Models. We report results for two forecasting models: (i) the panel data Forecast Evaluation. Because of the exponential transformation in (7) from growth rates to levels, there is a large degree of cross-sectional heterogeneity among the levels of infection. Locations with larger numbers of infections tend to be associated with larger forecast errors. If we simply average forecast errors or forecast interval lengths across locations, the results will be driven by a few locations with a high level of infections. Therefore, we are standardizing all level-forecast evaluation statistics by the level of infections at the forecast origin, I iT , i.e., we are reporting results for the forecast of I iT +h /I iT .

We will report measures of density and interval forecasting performance below. We do not consider point forecasts because we strongly believe that due to the highly uncertain path of infections during a pandemic it is essential for forecasters to report forecasts that convey the degree of uncertainty in the predictive distribution. The density forecast performance is evaluated based on continuous ranked probability scores (CRPS). The CRPS measures the L 2 distance between the cumulative density functionF iT +h|T (x) associated with a predictive distribution for location i at forecast origin T and a perfect probability forecast that assigns probability one to the realized x iT +h :

The CRPS is a proper scoring rule, meaning that it is optimal for the forecaster to truthfully reveal her predictive density. Here x iT +h could either be a growth rate y iT +h or a relative level I iT +h /I T .

For interval forecasts we will report the cross-sectional coverage frequency and the average length separately. As discussed in more detail in Askanazi, Diebold, Schorfheide, and Shin We begin with the top left panel of Figure 3 . In most early-group locations, the infections tend to peak between the forecast origin T = 63 and the four-week-ahead forecast target T + h = 91. In the late-group locations, the peak occurs after the forecast target date.

For T = 63 the three important findings emerge. First, the panel forecasts clearly dominate the time-series forecasts. The discrepancy is particularly large for locations in the late group. Second, while for the early group the CRPS based on the panel forecasts seem to be unrelated to the peak date, the accuracy of the time-series forecasts is substantially worse for early-group locations that peak between periods 63 and 91 than it is for locations that peak prior to period 63. Third, the four-week-ahead panel forecasts for the late group are much more accurate than the panel forecasts for the early group.

These findings can be explained as follows. First, in a panel setting, the experience of the early locations allows for relatively precise inference about the RE distribution, which then sharpens the posterior inference for the late locations because the uncertainty about the prior distribution is reduced. Note that the time series dimension for the late group is only 7. Second, due to the structural break in the growth rate at the peak infection level, it is very difficult to predict how quickly the infections will die out after they have peaked. This makes it easier to predict infections for the late group which includes the locations that are still far away from the peak than for the early group in which infection levels are relatively close to the peak.

The top right panel of Figure 3 indicates that after 18 weeks (T = 126) the benefit of the panel approach is a lot smaller, both for the early group and the late group. Because more time series information is available to estimate the location-specific parameters, the benefit from using prior information is significantly diminished. The bottom panel of the figure shows CRPS for levels rather than growth rates. The key message remains the same: early on in the pandemic, the panel approach substantially improves forecasts for locations that experience a delayed outbreak, because there is some learning from the locations in which the outbreak occurred early on.

In Figure 4 we plot the group-specific average CRPS as a function of the forecast origin are similar to the messages from Figure 3 , but now the results span a broad range of forecast origins. First, the panel forecasts are (at least weakly) more accurate than the time-series forecasts. However, the accuracy differential vanishes as the time-series dimension of the estimation sample increases over time. Second, the benefit from using a panel approach is more pronounced for the locations that experience a late outbreak than those that experience an early outbreak.

Interval Forecast Accuracy. Finally, we report results on the interval forecast performance for infection growth rates and levels in Figure 5 . 

We apply the panel forecasting techniques to country/region-level data on active Covid-19

infections. The data set used in the empirical analysis is described in Section 5.1. We discuss the posterior estimates in for the 2020-04-18 forecast origin in Section 5.2. In Section 5. 

The data set is obtained from CSSE at Johns Hopkins University. 4 We define the total number of active infections in location i and period t as the number of confirmed cases minus the number of recovered cases and deaths. Throughout our study we use country-level aggregates. The time period t corresponds to a day and we fit our model to one-sided three-day rolling averages to smooth out noise generate by the timing of the reporting. In a slight abuse of notation, the time subscript t in (3) 

Before discussing the forecasts, we will examine the parameter estimates for one of the early samples, namely 2020-04-18.

Heterogeneous Slope Coefficients. Our Gibbs sampler generates draws from the joint posterior of (ρ, λ 1:N , σ 2 1:N , ξ)|Y 1:N,0:T . We begin with a discussion of the estimates of γ 1i and δ 1i , which affect the speed at which the growth rates are expected to change on a daily basis. γ 1i measures the average daily decline in the growth rate of active infections. For instance, suppose the at the beginning of the outbreak, in event time t = 0, the growth rate ln(I t /I t−1 ) = 0.2, i.e., approximately 20%. A value of γ 1i = −0.02 implies that, on average, the growth rate declines by 0.02, meaning that after 10 days it is expected to reach zero and turn negative subsequently. A positive value of δ 1i = 0.01 implies that after the growth rate becomes negative, its decline is reduced (in absolute value) to γ 1i + δ 1i = −0.01. peaked will take considerably longer than the rise to the peak.

Distribution. An important component of our model is the RE distribution π(λ i |ξ) defined in (9). Prior and posterior uncertainty with respect to the hyperparameters ξ generate uncertainty about the RE distribution. In the remaining panels of Figure 6 we plot draw from the posterior (center column) and prior (right column) distribution of the RE density π(λ i |ξ). Each draw is represented by a hairline. Because the normalization constant C(ξ) of π(λ i |ξ) is difficult to compute due to the truncation of a joint Normal distribution, we show kernel density estimates obtained from draws from π(λ i |ξ). 

We now turn to density forecasts generated from the estimated panel data model. For now, we will focus on the early stage of the pandemic. We use Algorithm 1 to simulate trajectories of infection growth rates which, conditional on observations of the initial levels I iT , we convert into stocks of active infections. For each forecast horizon h we use the values y s iT +h and I s iT +h , s = 1, . . . , N sim to approximate the predictive density. Strictly speaking, we are not reporting complete predictive densities. Instead, we plot medians and construct equal-tail-probability bands that capture the range between the 20-80% and 10-90% quantiles. The wider the bands, the greater the uncertainty. infections. The path of active infections broadly resembles the paths simulated with the SIR model in Section 2. The rise of infections during the outbreak tends to be faster than the subsequent decline, which is a feature that is captured by the break in the conditional mean function of our model for the infection growth rate y it in (3). The difference between the bands depicted in the second and third rows is that the former reflects parameter uncertainty only (we set future shocks equal to zero), whereas the latter reflects parameter and shock uncertainty. In the case of Germany, shock uncertainty increases the width of the bands by approximately 30%.

Due to the exponential transformation that is used to recover the levels, the predictive densities are highly skewed and exhibit a large upside risk. This is particularly evident for the U.S. The growth rate prediction in the first row indicates that there is an approximately 20% probability of a positive infection growth rate throughout April and at least a 10% probability until the middle of June. Converted into levels, temporarily positive growth rates of infections can generate a rise of infections from less than one million in April to more than five million two months later.

In the bottom row of Figure 8 we plot cumulative density function for the date of recovery, which we define as the first date when the infections fall below the initial level I i0 . The density function is calculated by examining each of the future trajectories I s iT +h for h = 1, . . . , 60 generated by Algorithm 1. For South Korea the probability that the infection rate will fall below I i0 over the two month period is close to 80%, whereas for Germany and the U.S. the probability is approximately 50% and 60%, respectively.

In Figure 9 we overlay eight weeks of actual infections onto density forecasts generated 

We now turn to a more systematic evaluation of the forecasts and will assess the accuracy of density and interval forecasts represented by the bands in Figures 8 and 9 . For reasons previously discussed in Section 4.2, we standardize future infections I iT +h by the level of infections I iT at the forecast origin. A closer inspection of the forecasts for more than 100 countries/regions reveals that the long-run forecasting performance is not particularly good. This is not just a feature of our panel trend-break model, but also a feature of other epidemiological models such as the SIR model for which we will report results below. Thus, in this section we will focus on one-week and four-week ahead forecasts and not report results for an eight-week horizon.

Alternative Models. In addition to the panel model forecasts, we consider two alternative forecasts. First, as in Section 4, we generate time-series forecasts based on the trend-break model (3) for each location. Second, we estimate a version of the simple SIR in ( 1) with time-varying parameters β t and γ t . Notice, that by rewriting (1) we can express β t and γ t directly as a function of the observables (here we are omitting i subscripts): 8

This allows us to estimate the AR(1) law of motion in (2) for each country using Bayesian techniques. The AR(1) models are then used to simulate trajectories ( β T +1:T +H , γ T +1:T +H ) from the posterior predictive distribution. 9 For each parameter sequence, we iterate the SIR model (1) forward to obtain a predictive distribution of the active infections.

Density Forecast Accuracy. Figure 10 summarizes the one-week-ahead density forecasting performance for once-a-week forecast origins starting on 2020-04-18 and ending on 2020-07-04. For each location, we compute the probability score CRPS i,T +h|T . The top row shows the cross-sectional median as a function of the forecast origin, whereas the center and the bottom row show the cross-sectional empirical distribution for two forecast origins: 2020-04-18 and 2020-06-06.

The panels in the left column of Figure 10 cover all locations, whereas the panels in the right column distinguish between early-group and late-group locations. The early group comprises locations that experienced more than 100 infections before 2020-03-28. The remaining 8 The following additional variables are obtained from the JHU CSSE dataset: N is the total population of each country. S t is computed as N -I t -recovered cases -deaths. 9 Based on the specification of the SIR model, we let β t , γ t > 0 and 0 ≤ S t , I t , R t ≤ N , for all t. The panels in the right column of Figure 10 distinguish between locations that experienced the Covid-19 outbreak at an early stage and locations that were hit by the pandemic at a later stage. The key result is that for forecast origins dated 2020-05-09 or earlier, the panel forecasts for the late group are more accurate than the time series forecasts from the trend-break model. This result confirms the basic intuition that the panel approach can be advantageous during a slowly spreading pandemic because the experience of the early-group countries can sharpen inference on the RE distribution for the latter countries.

Unfortunately, because the time series approach dominates the panel approach for the early countries, in the aggregate there is no clear advantage to the panel analysis in our data set.

The left panels of Figure 10 The panel data forecasts have a smaller average length than the individual-level forecasts for both groups and in the aggregate. Thus, on balance, in terms of interval forecasting, the panel approach comes out slightly ahead. Finally, the bottom right panel shows that the interval forecasts for the late group are generally wider than for the early group. The additional uncertainty is caused by the difficulty of predicting the change in infection growth rates around the peak. Figure 13 displays results for a four-week horizon. Over this longer horizon, the coverage frequency is generally poor. As for the shorter horizon, the SIR model interval forecasts are substantially worse in terms of coverage frequency and interval length than the panel and time-series forecasts from the trend-break model.

We adopted a panel forecasting model initially developed for applications in economics to forecast active Covid-19 infections. A key feature of our model is that it exploits the experience of countries/regions in which the epidemic occurred early on, to sharpen forecasts and parameter estimates for locations in which the outbreak took place later in time. At the core of our model is a specification that assumes that the growth rate of active infections Coverage Probability

Notes: The nominal coverage probability is 80%. Left column panels: solid is panel, dashed is country-level, and dashed-dotted is SIR. Right column panels: solid is panel, dashed is country-level. Blue lines correspond to early group and orange lines to late group.

trend function with a break. Our specification is inspired by infection dynamics generated from a simple SIR model. According to our model, there is a lot of uncertainty about the evolution of infection rates, due to parameter uncertainty and the realization of future shocks. Moreover, due to the inherent nonlinearities and exponential transformations, predictive densities for the level of infections are highly skewed and exhibit substantial upside risk. Consequently, it is important to report density or interval forecasts, rather than point forecasts.

A natural extension of our model is to allow for additional, data-determined breaks in the deterministic trend function as the pandemic unfolds and countries/regions are adopting 

On the Comparison of Interval Forecasts

Policy Implications of Models of the Spread of Coronavirus: Perspectives and Opportunities for Economists

The Challenges of Modeling and Forecasting the Spread of COVID-19

Nonparametric Empirical Bayes and Compound Decision Approaches to Estimation of a High-dimensional Vector of Normal Means

Splinefunktionen

The Macroeconomics of Epidemics

Estimating and Simulating a SIRD Model of COVID-19 for Many Countries, States, and Cities

Health versus Wealth: On the Distributional Effects of Controlling a Pandemic

Unobserved Heterogeneity in Income Dynamics: An Empirical Bayes Perspective

The Probability Approach in Econometrics

Going Viral: Forecasting the Coronavirus Pandemic Across the

Containing papers of a mathematical and physical character

Macroeconomic Dynamics and Reallocation in an Epidemic

Nowcasting unemployment insurance claims in the time of COVID-19

When Will the Covid-19 Pandemic Peak?

Density Forecasts in Panel Data Models: A Semiparametric Bayesian Perspective

Forecasting with Dynamic Panel Data Models

Estimating an Inverse Gamma distribution

Forecasting the Impact of the First Wave of the COVID-19 Pandemic on Hospital Demand and Deaths for the USA and European Economic Area Countries

Dealing with Data Gaps