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.
117,653 characters · 21 sections · 97 citation commands
Prediction intervals for economic fixed-event forecasts
Economic forecasts are often published in the `fixed-event' format. Consider tagesschau2022, a German news website that maintains a list of recent GDP forecasts made by various institutions. For example, in April 2022, the federal government predicted 2.2% GDP growth for 2022, followed by 2.5% in 2023. In July 2022, the European Commission predicted 1.4% growth for 2022, followed by 1.3% in 2023. These forecasts are called `fixed-event' because the quantity of interest -- the GDP growth rate over the year 2022 or 2023 -- remains fixed, whereas the date at which the forecast is made moves forward in time. This forecasting format is routinely used by institutions that make economic forecasts, and by news media which comment on these forecasts.
Economic fixed-event forecasts are typically published without a quantitative measure of uncertainty. Such a measure would be particularly important in the current setup: By construction, forecast uncertainty varies heavily depending on when the forecast is made and which year it refers to. For example, in July 2022, the EU commission's forecast for 2022 should be less uncertain than its forecast for 2023. But just how much less uncertain depends on the time series properties of GDP, and is hard to grasp intuitively.
Motivated by this situation, the present paper constructs measures for the uncertainty in fixed-event point forecasts. Together with the point forecasts themselves, we can then construct forecast distributions for GDP growth and other economic variables. The principle of forecast post-processing -- using point forecasts and past forecast errors to construct forecast distributions -- is popular in meteorology (see e.g. GneitingRaftery2005, GneitingEtAl2005, RaspLerch2018 and Vannitsem2021), economics (e.g. Knueppel2014, KruegerNolte2016 and Clark2020) and other fields. State-of-the-art point forecasts are often publicly available, so that using them as a basis for forecast distributions is more practical than generating forecast distributions from scratch. Furthermore, assessing forecast uncertainty based on past errors does not require knowledge about how the point forecasts were generated. This is an important advantage in practice, where the forecasting process may be judgmental, subject to institutional idiosyncracies, or simply unknown to the public.
The vast majority of the postprocessing literature considers a `fixed-horizon' forecasting setup where the time between the forecast and the realization remains constant. Examples of the fixed-horizon case include daily forecasts of temperature 12 hours ahead, or quarterly forecasts of the inflation rate between the current and next quarter. In economics, Clements2018 constructs a measure of fixed-event forecast uncertainty that is based on fixed-horizon forecast errors, and thus requires that an appropriate database of fixed-horizon forecasts is available. This is the case in the US Survey of Professional Forecasters analyzed by Clements2018, but not in other situations including the German GDP example mentioned earlier. Existing approaches are thus not applicable to the fixed-event case. Instead, the latter requires different tools which we develop in this paper.
The main idea behind our proposed approach is simple: We model quantiles of the forecast error distribution as a function of the forecast horizon. The latter is defined as the time (measured in weeks) between the forecast and the end of the target year. For example, forecasts made on July 1, 2022 for 2022 and 2023 correspond to horizons $26$ and $78$, respectively. To estimate regression models, we use a dataset of past forecast errors at different horizons. Pooling forecast errors across horizons allows us to estimate uncertainty at any given horizon by considering uncertainty at neighboring horizons. This approach is helpful in the fixed-event case, where only a small number of past errors is typically available for a given horizon. For example, the German data set we consider covers forecast-observation pairs ranging from 1991 to 2022, and includes precise information (daily time stamps) on the forecast horizon $h$. Specifically, the data contains $525$ unique values for $h$, ranging from $h = 0$ to $h = 104$. Note that $h$ need not be an integer; for example, forecasts made on the seven days of the target year's final week correspond to horizons $h \in \{0, 1/7, \ldots, 6/7\}$. Given that the data covers $n = 1\,307$ observations in total, the average number of forecast errors corresponding to each of the $525$ different horizons is about $2.5$. Using only forecast errors that correspond exactly to some horizon of interest ($h = 5$ weeks, say) is hence not a promising strategy, and considering forecast errors from neighboring horizons seems advisable. This aspect is not relevant when postprocessing fixed-horizon forecasts, which is typically based on a time series of past forecast errors for the exact horizon of interest.
The statistical methods we consider incorporate constraints that are motivated by the fixed-event forecasting problem. In particular, it is plausible to assume that forecast uncertainty increases monotonically across horizons, and levels off at some horizon. We consider three main approaches that implement this idea: A Gaussian heteroscedastic model, a decomposition approach that imposes a symmetry assumption, and a flexible approach that is (almost) nonparametric. The first two approaches are parsimonious in light of typically short samples of past forecast errors. Despite their simplicity, they perform well in a cross-validated analysis of German and US data. In particular, the prediction intervals attain coverage close to its nominal level, and they clearly outperform benchmark predictions by a survey of professional economists for the US data.
To illustrate the quantitative implications of our results, consider a hypothetical forecaster who, as of mid-September, issues a point prediction of 2.1% GDP growth for the current year, and 1.7% for the next year. Our results for the German data then suggest that a plausible 80% prediction interval for the current year would be on the order of $2.1 \pm 0.35 = [1.75,2.45]$, compared to $1.7 \pm 2 = [-.3, 3.7]$ for the next year, i.e., the prediction interval for the next year is more than five times as wide. We thus argue that common statements along the lines of `we expect 2.1% GDP growth this year, and 1.7% next year' are not helpful: They implicitly put the current- and next-year forecasts on an equal footing, and over-emphasize the next-year point forecast which is surrounded by considerable uncertainty.
The rest of this paper is structured as follows. Section (ref) introduces notation for the fixed-event setup. In order to illustrate the statistical implications of this setup, Section (ref) considers fixed-event forecasting in an autoregressive time series model. Section (ref) introduces the proposed regression methods for forecast postprocessing. Section (ref) describes methodology for forecast evaluation and model selection. Section (ref) studies the performance of the proposed methods in a simulation experiment, and Section (ref) considers empirical fixed-event forecasts from Germany and the US. Section (ref) concludes with a discussion. The online supplement contains derivations and further analyses. Replication materials are available at \url{https://github.com/FK83/gdp_intervals}.
We consider forecasting $Y_t$, the real GDP growth rate in year $t$. The latter is defined as $100 \times \frac{\text{GDP}_t-\text{GDP}_{t-1}}{\text{GDP}_{t-1}},$ where $\text{GDP}_t$ is the level of real GDP in year $t$, which in turn is taken to be the average of the year's four quarterly levels. For many countries, point forecasts of $Y_t$ are available from various public and private forecasting institutions. We denote the time between the origin date (on which the forecast is issued) and the end of year $t$ as the `forecast horizon', denoted by the symbol $h$ and measured in weeks. Let $X_{t,h}$ denote the point forecast of $Y_t$ at horizon $h$. For example, for year $t = 2020$, a forecast issued on December 17, 2020 corresponds to horizon $h = 2$. We focus on forecast horizons $h \in [0, 104]$, i.e., up to two years ahead, which covers most practical economic forecasts. Even a forecast at horizon $h = 0$ need not be perfect in practice since a precise estimate of the outcome may be available only after the end of year $t$. For simplicity, we treat each year as having $52$ weeks.
The forecast error of $X_{t,h}$ is given by $e_{t,h} = Y_t-X_{t,h}$. Let $\mathcal{F}_{t,h}$ denote a suitable information set that is available $h$ weeks before the end of year $t$. In order to obtain a forecast distribution for $Y_t$ given $\mathcal{F}_{t,h}$, or $Y_t|\mathcal{F}_{t,h}$ in short, we model the distribution of $e_{t,h}|\mathcal{F}_{t,h}$. That is, we take the point forecast $X_{t,h}$ as given, and focus on modeling the error of this point forecast. Since $$\mathbb{P}(Y_t \le y|\mathcal{F}_{t,h}) = \mathbb{P}(e_{t,h} \le (y-X_{t,h})|\mathcal{F}_{t,h}),$$ and $X_{t,h}$ is known given $\mathcal{F}_{t,h}$, the distribution $e_{t,h}|\mathcal{F}_{t,h}$ can be used to construct the desired distribution for $Y_t|\mathcal{F}_{t,h}$.
Our approach of modeling the error of a given point forecast, rather than constructing a forecast distribution from scratch, is called `postprocessing'.\footnote{In meteorology, the term `ensemble postprocessing' is more common since postprocessing is typically applied to a collection (`ensemble') of point forecasts stemming from a numerical weather prediction model GneitingRaftery2005. That said, the broader principle of modeling forecast errors readily transfers to the case of a single point forecast that we consider here.} However, all postprocessing studies we are aware of consider fixed-horizon forecasting. Our fixed-event setup requires different statistical tools in that only a small number of observations is typically available for a given forecast horizon. We thus focus on modeling the properties of $e_{t,h}$ as a continuous function of $h$, thus interpolating across forecast horizons. By contrast, postprocessing methods for fixed-horizon forecasts typically treat each horizon separately, using a time series $(e_{t,h})_{t=1}^n$ of past forecast errors at a single horizon $h$.
In this section, we consider fixed-event forecasting in an autoregressive Gaussian time series model from the econometric literature. The stylized facts illustrated here will later motivate our more general empirical methodology (described in Section (ref)).
Our model is a modified version of the one in PattonTimmermann2011. It assumes that fixed-event forecasts are based upon noisy high-frequency observations, but make correct use of these observations. That is, forecasts are equal to the true conditional mean of the predictand, given their information base. While high-frequency data on GDP is not literally available in practice (where GDP is measured at a quarterly frequency only), there are various efforts at measuring economic activity based on economic variables that are available at monthly, weekly or daily frequency AruobaEtAl2009,BraveEtAl2019,LewisEtAl2021,EraslanGoetz2021. Furthermore, publication lags and ex-post revisions in macroeconomic data Croushore2001 imply that realizations become available with a delay. Taken together, the model's assumption of a noisy proxy observed at high frequency hence provides a plausible (albeit abstract) representation of practical GDP forecasting.
We consider hypothetical weekly observations, whereas PattonTimmermann2011 use hypothetical monthly observations. This increase in granularity is motivated by our empirical data setup, where many forecasts are made within the month, so that a monthly frequency may be too coarse. Section A.1 of the online supplement provides a detailed comparison of our model to the one of PattonTimmermann2011. We approximate $Y_t$, the GDP growth rate from year $t-1$ to year $t$, as the weighted sum of $103$ weekly logarithmic growth rates ranging from the beginning of year $t-1$ to the end of year $t$. As detailed in the online supplement, this setup arises from the definitorial convention ('annual-average') that we use to compute $Y_t$. We denote the weekly logarithmic growth rate by $Y_w^*$, with the understanding that each year $t$ corresponds to a distinct set of 52 index values $w$ (see below for an example). Here and henceforth, we use the `star' superscript notation for weekly random variables. We further assume that $Y_w^*$ follows a first-order autoregression, so that
where $w_\text{last}(t)$ denotes the index of the last week of year $t$, and $\gamma_j = 1-\frac{|52-j|}{52}$. Note that the coefficients $\gamma_j$ form a triangle when plotted against $j$. For example, suppose that the sample starts in year $t = 1$ (containing weeks $w = 1, 2, \ldots, 52$), so that week $w = 104 = w_\text{last}(2)$ is the last week of year $2$. According to Equation ((ref)), the approximate annual growth rate $Y_2$ is a weighted sum of the weekly observations $(Y_w^*)_{w=2}^{104}$. The greatest weight of $1$ is associated with $Y^*_{53}$, and the smallest weight of $1/52$ is associated with $Y^*_{104}$ and $Y^*_{2}$. Furthermore, $\sum_{j=1}^{103} \gamma_j = 52,$ in line with the fact that $Y_w^*$ is measured at a weekly frequency whereas $Y_t$ is measured annually. See Section A.2 of the online supplement for a concise derivation of the approximation at ((ref)), and PattonTimmermann2011 for numerical evidence on the high precision of the approximation.
As noted earlier, we assume that forecasters observe a noisy version of $Y_w^*$. Specifically, let $\tilde Y_w^* = Y_w^* + \eta_w^*,$ where $\eta_w^* \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0,\sigma^2_\eta)$ is an independent and identically distributed (IID) Gaussian measurement error. An $h$-week ahead mean forecast of $Y_t$ is then based on the information set $\mathcal{F}_{t,h}$ generated by the sequence of noisy weekly observations $(\tilde Y_w^*)_{w=1}^{w_\text{last}(t)-h}.$ Due to the presence of measurement error, the optimal forecast of $Y_t$ given $\mathcal{F}_{t,h}$ has no simple closed-form expression. However, the optimal forecast can be computed analytically by means of the Kalman filter. The latter uses the model's linear state space representation, which we describe in Section A.3 of the online supplement.
An important implication of the present model is that forecasts of $Y_t$ at horizon $h = 0$ are not perfect, which is in line with reality. For example, the GDP for 2023 is not yet published on December 31, 2023, so that the $0$-week ahead forecast of $Y_t$ is formed on the basis of preliminary estimates.
The model allows to study the root mean squared forecast error (RMSFE) of fixed-event point forecasts as a function of the forecast horizon. This relationship represents the amount of predictability (or lack thereof) at various points in time, and will be relevant for specifying appropriate postprocessing approaches below. Figure (ref) presents the relationship. We have set the three model parameters $\rho = 0.3$ (persistence of weekly GDP growth), $\sigma^2_\varepsilon = 0.09$ (variance of white noise error in weekly GDP growth) and $\sigma^2_\eta = 0.003$ (variance of measurement error) that yield a plausible approximation of the empirical RMSFE values observed in the German GDP data that we consider below. The figure shows that the model-implied RMSFE increases as the forecast horizon increases. This type of monotonicity is a well-known property of forecasts that are optimal given some information set Patton2012,Krueger2021, which is satisfied here. The shape of the function further indicates a rather steep and roughly linear increase in RMSFEs between $h \approx 20$ and $h \approx 70$, whereas RMSFEs are almost constant for $h \le 10$ and $h \ge 80$. In Section A.4 of the online supplement, we explore the RMSFE curves implied by other values of the model's parameters.
Our empirical analysis below is based on a sample of past forecast errors $(e_{t_i, h_i})_{i=1}^n$, where $t_i$ and $h_i$ denote the target year and horizon of the $i$th forecast error. We next discuss the correlation structure of such a sample in the context of the model considered above. For simplicity, we assume that $\sigma^2_\eta = 0$, such that weekly data $Y_w^*$ are observed without error.
Consider two forecast error observations $e_{t_1, h_1}$ and $e_{t_2, h_2}$. Without loss of generality, let $t_2 \ge t_1$. Furthermore, we focus on the case $1 \le h_2 \le 104$ (forecast errors ranging from one week to two years). The two forecast errors are independent if one of the following sufficient conditions holds:
see Section A.5 of the online supplement for details. Intuitively, both conditions rule out any overlap between the two time intervals $[t_1-h_1, t_1]$ and $[t_2-h_2, t_2]$ that range from the forecast date to the target date of the two forecast errors. If neither of the conditions holds, the two forecast errors will typically be dependent.
More broadly, the dependence structure in a sample of forecast errors $(e_{t_i, h_i})_{i=1}^n$ thus reveals a cluster-type pattern. If $h \le 52$, then the forecast errors belonging to a given year (e.g., all errors $e_{t_i, h_i}$ s.t. $t_i = 2019$) form one cluster, with errors being dependent within the same cluster but independent across different clusters. If $h_i \ge 53$ is possible, as is the case in our empirical analysis, then the dependence structure is more complicated, with possibly nonzero correlation across neighboring years. If we allow for the case $\sigma^2_{\eta} > 0$, the dependence structure becomes even more complicated since any pair of forecast errors may be jointly affected by updated `back-casts' (i.e., corrected assessments of past data, as produced by the Kalman filter) that are relevant in the presence of measurement error.
HansenLee2019 present asymptotic results on samples with cluster dependence, and on estimators based upon such samples. Similar results can be derived in the present setup, and allow to derive asymptotic statements as the number of clusters tends to infinity. However, in view of the rather short data samples used in our empirical analysis (about $32$ years for the German data, and about $42$ years for the US data), it is unclear whether such results provide a good description of our setup. In Section (ref), we hence conduct a simulation study that closely mimics our empirical setup.
The model discussed in Section (ref) is useful to illustrate the main stylized facts of the fixed-event forecasting problem. PattonTimmermann2011 estimate the model's parameters and use the model for predicting the distribution of forecast errors as a function of $h$. While conceptually appealing, estimation is challenging in practice, with parameter estimates differing markedly across estimation methods PattonTimmermann2011. Addressing these challenges seems unnecessary in the current setup, where interest lies on forecasting (as opposed to interpreting the model's structural parameters). In the following, we therefore focus on empirical models that are considerably simpler to implement and are less restrictive in terms of functional form assumptions. In the simulation experiments from Section (ref), we demonstrate that these empirical methods perform well even when the true data-generating process is given by the model from Section (ref).
To simplify model building, we assume that the distribution of $e_{t,h}|\mathcal{F}_{t,h}$ is a function of $h$ alone, i.e. that $\mathbb{P}(e_{t,h} \le z|\mathcal{F}_{t-h}) = \int_{-\infty}^z dF_h(z),$ where $F_h$ is a distribution that depends on $h$ but is constant over time. This assumption is satisfied, for example, for the steady state version of the autoregressive model from Section (ref).\footnote{Here `steady state' means that a sequence of initial observations has been removed, so that the prior mean and variance for the initial state vector become irrelevant.} From a pragmatic perspective, the assumption of a time-invariant conditional distribution function is motivated by the paucity of data in our empirical analysis, which implies that estimating elaborate conditional distributions does not seem promising. However, we retain the dependence on $h$, which we expect to be a major determinant of forecast error distributions, with larger values of $h$ corresponding to more variable forecast errors.
In Sections (ref) to (ref), we describe three modeling approaches, all of which incorporate constraints that are motivated by the fixed-event forecasting setup. On the other hand, the stringency of the imposed constraints differs across approaches.
Our first, most restrictive approach assumes that $e_{t,h} \sim \mathcal{N}(\mu, \sigma_h^2),$ with
note that the second factor in ((ref)) is the cumulative distribution function (CDF) of a logistic random variable with location $\theta_2$ and scale $\theta_3 > 0$, evaluated at $h$. We hence model the standard deviation of $e_{t,h}$ as a constant ($\theta_1 > 0$) times a function that is monotonically increasing in $h$. As noted above, the monotonicity of $\sigma_h$ is an implication of optimal forecasting. In practice, we expect monotonicity to hold also under moderate forms of sub-optimality. The logistic functional form of $\sigma_h$ is motivated by the structural model displayed in Figure (ref). The specification at ((ref)) implies that $\sigma_h \rightarrow \theta_1 \equiv \sigma_\infty$ as $h \rightarrow \infty$, and $\sigma_h \rightarrow 0$ as $h \rightarrow -\infty$. Both implications seem economically plausible, noting that negative horizons ($h < 0$) correspond to back-casts that, for sufficiently small values $h << 0$, are no longer affected by data revisions or missing input data. The parameter $\theta_2$ denotes the horizon $h^{'}$ that satisfies $\sigma_{h^{'}} = 0.5~\sigma_\infty$, i.e., the horizon at which half of the uncertainty is resolved, whereas $\theta_3$ is a shape parameter that provides further flexibility. Figure (ref) illustrates the functional form of $\sigma_h$ for different parameter choices. The mean parameter $\mu$ is a parsimonious way to allow for nonzero forecast errors resulting from biased point forecasts. While more sophisticated specifications are possible (with the bias depending, for example, on the forecast horizon $h$), the empirical results in Section (ref) (in particular, Table S2 in the online supplement) indicate that the empirical importance of the bias is quite limited in our setup.
We estimate the Gaussian model by minimizing the continuous ranked probability score MathesonWinkler1976, as implemented in the scoringRules software package JordanEtAl2019 for R R. The CRPS is a loss function for distributions, and has become a popular alternative to the log likelihood function (also called logarithmic score) in recent years. The CRPS is a strictly proper scoring rule, that is, a forecaster has an incentive to state what they think is the true forecast distribution. A general advantage of the CRPS over the logarithmic score is that it does not require a density, and can easily handle discrete or empirical distributions JordanEtAl2019,KruegerEtAl2021. When used as a criterion for estimation, the CRPS yields consistent parameter estimates under standard conditions including correct specification GneitingRaftery2007. We use the CRPS for parameter estimation since it is related to the interval score that we use for forecast evaluation (see Section (ref) for details), thus broadly aligning the criteria used for model estimation versus evaluation.
Our second model is based on a decomposition of $e_{t,h}$ into its sign and its (absolute) magnitude. Similar decompositions have been considered for modeling financial returns ChristoffersenDiebold2006,AnatolyevGospodinov2010, based on the motivation that the sign of returns is far less predictable than their magnitude. Noting that a similar motivation applies to forecast errors, we adopt a similar decomposition method here. To describe the method, note that $$e_{t,h} = \underbrace{(2~\mathbf{1}(e_{t,h} > 0)-1)}_{=S_{t,h}}~|e_{t,h}|,$$ where $\mathbf{1}(A)$ is the indicator function of the event $A$. We assume that the sign $S_{t,h}$ takes values of $\pm 1$ with equal probability, and is independent of the magnitude $|e_{t,h}|$. We further assume that the CDF of $|e_{t,h}|$, which we denote by $G_h$, is constant across time $t$ but (possibly) different for each forecast horizon $h$. Under these assumptions, we obtain
Finally, we assume that $|e_{t,h}|~{\succsim}_{\text{FSD}}~ |e_{t,h-v}|$ for any $v > 0$, where the notation ${\succsim}_{\text{FSD}}$ indicates first-order stochastic dominance (FSD). This means that the distribution of absolute forecast errors becomes stochastically greater as the horizon increases. This assumption seems plausible in the present context, and can be motivated as a more stringent version of increasing mean squared forecast errors across horizons.\footnote{Specifically, the assumption implies that $e_{t,h}~{\succsim_\text{CX}}~e_{t,h-v}$ for $v > 0$, where $\succsim_\text{CX}$ denotes convex order ShakedShanthiku2007. The latter relation implies that $\mathbb{E}(e_{t,h}^2) \ge \mathbb{E}(e_{t,h-v}^2),$ i.e., increasing mean squared forecast errors.} In order to estimate $G_h$ under the present assumption, we use isotonic distributional regression (IDR) as recently studied by HenziEtAl2021 and implemented in the R package isodistrreg HenziEtAl2022. Briefly, IDR estimates a conditional distribution function, subject to the constraint that the outcome increases (in the sense of FSD) as the predictor vector increases (with various notions of `increases' being covered, including partial orders on the covariate space). Our empirical setup is a fairly simple special case, in that we use a single continuous predictor ($h$). See HenziEtAl2021 for an instructive example of IDR in this case. IDR has two main advantages over unconstrained nonparametric estimators of the distribution of $|e_{t,h}|$ given $h$. First, imposing the assumption serves to regularize the estimator and reduce estimation noise. Second, IDR is free of tuning parameters, whereas bandwidth parameters are often crucial for the performance of conventional nonparametric estimators. Note that for a given horizon $h$, Equation ((ref)) implies that a forecast distribution for $e_{t,h}$ obtained via the decomposition approach is necessarily symmetric around zero.
Our most flexible model estimates $F_h(e_{t,h})$ nonparametrically, subject to the constraint that $e_{t,h}~{\succsim}_{\text{ICX}}~ e_{t,h-v}$ for any $v > 0$, where ICX denotes increasing convex order. For two random variables $V$ and $W$, $V~ {\succsim}_{\text{ICX}}~W$ holds if $\mathbb{E}(\phi(V)) \ge \mathbb{E}(\phi(W))$ for all increasing convex functions $\phi$. Imposing the ICX constraint is hence less restrictive than imposing either FSD (which requires that the inequality holds for all increasing functions $\phi$) or convex order (which requires that the inequality holds for all convex functions $\phi$). To estimate the distribution of $e_{t,h}$ under the ICX constraint, we use the recent approach of Henzi2022. Similar to IDR, a major advantage of this approach is that it does not rely on tuning parameters. Section 1 of Henzi2022 provides further discussion of the ICX constraint, whereas Sections 5 and 6 present simulation and empirical examples. In particular, Equation (2) of Henzi2022 illustrates a data-generating process that satisfies ICX but violates FSD. Note that for a given horizon $h$, a forecast distribution for $e_{t,h}$ obtained via the flexible approach need not be symmetric around zero (or around any other value). This is because the approach does not assume a particular shape of the distribution $F_h$ for any given $h$. Instead, it assumes that two distributions $F_{h_1}$ and $F_{h_2}$ are ordered according to ICX.
The following result establishes that, in terms of the respective assumptions, the Gaussian modeling approach with $\mu = 0$ is a special case of the decomposition approach which, in turn, is a special case of the flexible approach.
In its steady-state version (see Footnote (ref)), the autoregressive model from Section (ref) satisfies the assumptions of the proposition's part (a). More specifically, the forecast at horizon $h$ is optimal given a certain information set (consisting of all noisy observations up to the forecast date). Together with the model's Gaussian setup, this implies that forecast errors $e_{t,h}$ are Gaussian with mean zero and variance $\sigma^2_h$, where the latter is increasing in $h$ Patton2012.
In addition to the three postprocessing methods introduced above, we consider a simple combination method: For a given quantile level (10% or 90%), we use the arithmetic mean of the three methods' forecasts. From a practical perspective, a main appeal of combinations is that their good performance is somewhat predictable, based on a large body of empirical evidence and theoretical findings on the properties of combinations GneitingRanjan2013,Lichtendahl2013,WangEtAl2022. By contrast, predicting the relative performance of several forecasting methods is typically more difficult in practice.
We focus on prediction intervals, and evaluate forecast accuracy with the interval score GneitingRaftery2007, BracherEtAl2021. In particular, we consider quantiles at levels $\alpha \in \{0.1, 0.9\}$, which together form the central 80% prediction interval. Such prediction intervals have been found to be a useful format for communicating probabilistic information to both expert and non-expert users Raftery2016. Specifically, let $l$ and $u$ denote the lower and upper bound of the prediction interval, and $y$ denote the realizing outcome. The interval score is then given by
with smaller scores being preferable. The score thus rewards short prediction intervals (with $u-l$ small) that nevertheless cover the realizing outcome $y$. The weight of ten on the penalty terms for not covering the outcome (as represented by the two indicator functions) ensures that the score is proper. That is, a forecaster minimizes their expected score by stating what they think is the true prediction interval.\footnote{The weight of ten for the penalty terms in Equation ((ref)) reflects the target interval coverage of $80 \%$. More generally, for the central $\kappa \%$ prediction interval, with $0 < \kappa < 100$, the weight is given by $200/(100-\kappa)$.} The interval score relates to the CRPS scoring rule mentioned earlier: While the interval score is proportional to the sum of two quantile scores Koenker1978,Gneiting2011 at levels $0.1$ and $0.9$, the CRPS is the unweighted integral over quantile scores at all levels $(0,1)$ GneitingRanjan2011.
The interval score is a summary measure of forecast performance, reflecting both sharpness (i.e., the informativeness of the forecast) and calibration (i.e., the consistency between the forecast and the outcomes). For prediction intervals, sharpness is represented by the length of the intervals, and calibration is represented by the intervals' coverage rate. In our simulation and empirical analyses, we thus report these two measures in addition to the interval score. For forecast case $i = 1, \ldots, n$, denote the lower and upper bound of the prediction interval by $l_i$ and $u_i$, and the associated realization by $y_i$. The average length of the prediction intervals is then given by $\text{AL} = n^{-1} \sum_{i=1}^n (u_i-l_i)$. The coverage rate is given by $\text{CR} = n^{-1} \sum_{i=1}^n \mathbf{1}(l_i \le y_i \le u_i)$, i.e., the share of observations for which the prediction interval covers the corresponding realization. Assuming that $\text{CR} < 1$, the average interval score for an evaluation sample of size $n$ can be written as
where ${\text{ASF}}_{\text{NC}} = \frac{1}{n(1-\text{CR})} \sum_{i=1}^n \text{sf}_i,$ with
i.e. ${\text{ASF}}_{\text{NC}}$ denotes the average shortfall (distance between observation and nearest end of prediction interval) in case of non-coverage.
A model's forecast performance should generally be assessed on data points that were not used for model fitting. In time series contexts, it is common to use rolling or expanding samples of data for estimation.\footnote{In a rolling sample, a one-step forecast for period $t$ is based on a training sample ranging from period $t-R$ to $t-1$, where $R$ is the length of the rolling sample. In an expanding sample, the traning sample ranges from periods $1$ to $t-1$, thus expanding over time.} Such model validation strategies are attractive if the time series of interest is sufficiently long. For short time series, these strategies are less attractive: Long (rolling or expanding) estimation samples mean that only few observations are left for model evaluation, yielding low power in model comparisons. Short estimation samples, on the other hand, may yield erratic model fits and unstable results. We seek to avoid these drawbacks, but nevertheless achieve a clear separation between training and test data. We thus use a cross-validation approach where predictions referring to year $t$ (i.e., referring to forecast errors $e_{t,h_j}$ for various horizons $h_j$) are based upon an estimation sample that comprises all years other than $t$. BergmeirEtAl2018 argue in favor of similar cross-validation strategies in the context of autoregressive time series models.
We next conduct a simulation study in order to assess the finite-sample properties of our methodology, including the models from Section (ref) and the evaluation techniques from Section (ref). To this end, we simulate data from the model presented in Section (ref), as well as optimal mean forecasts $X_{t,h}$ (see Section A.4 of the online supplement) and associated forecast errors $e_{t,h}$ for various horizons $h$. We then construct a random sample $\{e_{t_i, h_i}\}_{i=1}^n$ of forecast errors, where $t_i$ is drawn uniformly from $\{1, 2, \ldots, T_{\text{max}}\}$, and $h_i$ is drawn uniformly from $\{0, 1, \ldots, 104\},$ corresponding to a maximal forecast horizon of two years. This setup aims to mimic the heterogeneous and overlapping nature of the sample of forecast errors as described in Section (ref). We set $(n, T_\text{max})$ to either $(500, 30),$ $(1000, 30)$ or $(500, 60)$. Compared to the first sample type, the second type contains more observations from the same time period, whereas the third type contains the same number of observations from a longer time period. We further use the parameter values $\rho = 0.3, \sigma^2_\varepsilon = 0.09$ and $\sigma^2_\eta = 0.003$ that were also used in Figure (ref). Finally, we use a burn-in period of $30$ years in order to remove the impact of the prior parameters used to initialize the Kalman filter.
In addition to the error postprocessing methods described in Section (ref), we consider the true forecast error distribution that is implied by the structural model from which the data is simulated. Importantly, this distribution is not realistically available in practice, as it requires knowledge of the model's functional form and its true parameter values. The true distribution hence defines an unattainable gold standard that allows to set the other methods' performance in perspective.
Table (ref) displays the simulation results. The findings on the interval score (rightmost column of the table) can be summarized as follows. First, and as expected, the true model outperforms the three methods that are estimated based on samples of forecast errors. Second, the prediction methods' performance tends to improve as the sample becomes more informative (either by covering a longer time span or by covering more forecast errors from the same time span). The added value of training data is especially large for the most flexible method. Note that the expected performance of the true model is the same across all sample types, and any observed differences are within the range of Monte Carlo error.\footnote{Using two-sample $t$-tests, performance differences of the true model across any pair of sample types are not statistically significant at conventional levels.} Third, the Gaussian, decomposition and combination methods outperform the flexible method for all three sample types. The good performance of the decomposition method can be explained by its use of constraints (in particular, symmetry and a zero mean of forecast errors) that are satisfied by the true model, so that they serve to regularize the estimator. While the assumptions of the Gaussian method (in particular, the specification of the standard deviation $\sigma_h$) are not exactly satisfied by the true model, the degree of misspecification seems minor, so that the Gaussian method also benefits from regularization.
The results in Table (ref) further indicate that the Gaussian, decomposition and combination methods mostly reach good coverage rates of 78% or more, as well as prediction intervals of similar length. The flexible method's coverage rates are slightly worse, and its prediction intervals are somewhat shorter.
In the present simulation setup, the methods we consider thus yield plausible results, despite the challenges posed by small sample sizes and overlapping data.
We present results for two empirical data sets consisting of fixed-event point forecasts and associated realizations.
First, we consider a data set covering forecasts of German real GDP growth, as made by ten institutions: The Bundesbank, European Central Bank (ECB), European Comission (EC), International Monetary Fund (IMF), Organisation for Economic Co-operation and Development (OECD), as well as three German research institutes (DIW Berlin, ifo Munich, and IWH Halle) and two committees that play a prominent role in the policy debate (the German council of economic experts, and the Gemeinschaftsdiagnose, a joint forecast made by several research institutes). We obtain all forecasts and the corresponding first-release outcome data by the IWH Halle's forecasting dashboard dashboard,HeinischEtAl2023, a recent initiative that makes a rich archive of German economic forecasts openly available. The data indicate the origin date (day) and target date (year) of each forecast, so that precise information on the forecast horizon is available. We express the forecast horizon in weeks, noting that all methods we consider can accommodate non-integer forecast horizons (such as $h = 25/7$ weeks representing $25$ days). Various studies DoepkeFritsche2006,KoehlerDoepke2022 consider the properties of such fixed-event point forecasts for Germany. Foltas2022 analyze whether the distribution of forecast errors can be predicted by means of regressor variables. However, their analysis seeks to test (a particular notion of) forecast efficiency, as opposed to constructing forecast distributions. Furthermore, they treat each forecast horizon separately, while interpolating across horizons is the key methodological feature of our approach.
Importantly, we treat the forecasting institutions as exchangeable, that is, we do not attempt to model the forecast error as a function of the institution that made the forecast. This choice is motivated by empirical results from the forecast combination literature which suggest that treating economic forecasters as exchangeable often performs well in terms of the bias-variance trade-off GenreEtAl2013,ClaeskensEtAl2016. We expect similar findings to hold in the present case, especially in view of the rather small sample size. In the cases where two or more institutions make a forecast on the same day (11.2% of observations), we keep all forecasts. This means that we seek to predict the error made by a randomly drawn forecasting institution, as opposed to the forecast error made by an ensemble. This approach seems most practically relevant in the present context, in that most days feature at most one forecast.
Second, we consider US data from the Survey of Professional Forecasters Croushore2019 covering GDP and inflation forecasts from 1981 to 2021. The SPF's point forecasts are widely considered a hard-to-beat benchmark for even sophisticated statistical forecasting models FaustWright2013. Here we use fixed-event point forecasts that are available on a quarterly basis, referring to the present and next year. Since the forecasts are made roughly in the middle of a quarter (with each quarter corresponding to $52/4 = 13$ weeks), the forecast horizons satisfy $h \in \{6.5, 19.5, \ldots, 97.5\}$. A specific pair of two horizons is available each quarter. For example, in the first quarter of a year, the current-year forecast corresponds to $h = 45.5$ weeks, and the next-year forecast corresponds to $h = 97.5$ weeks. Furthermore, we consider the SPF forecast distributions that cover participants' subjective probabilities for various ranges of the outcome variable. We use the average survey response for both point and probabilistic forecasts from the SPF. We drop data from 1985:Q1, 1986:Q1 and 1990:Q1 due to possible technical errors in the associated survey rounds SPF_docu. At each forecast date, we approximate the probabilistic forecasts by a continuous distribution, following EngelbergEtAl2009 and using the implementation of KruegerPavlova2022. We use first-release data provided by the Federal Reserve Bank of Philadelphia to compute the actual GDP growth and inflation outcomes.
Table (ref) provides summary information about the data sets. Given the lack of easily available probabilistic forecasts for German GDP, the first data set is particularly interesting from a policy perspective. The US SPF data provide a useful testbed to assess the performance of our statistical models since they cover a longer time period with more variation in macroeconomic outcomes. Furthermore, the SPF's probabilistic forecasts provide a natural benchmark for our postprocessing methods.
Table (ref) presents our main empirical results, which are pooled across all forecast horizons. For the German data, all forecasting methods attain coverage rates close to the nominal level of $80 \%$. The Gaussian method produces the shortest prediction intervals, while the flexible method produces the widest ones. The interval score is highest (i.e., worst) for the flexible method. For the US data (GDP and inflation), the Gaussian, decomposition, and combination methods again attain coverage rates close to 80%, whereas the flexible method's rater are lower. On the other hand, the flexible method produces shorter prediction intervals, so that its overall performance (in terms of the interval score) is similar to the other three methods.
For the Gaussian method, we also considered a restricted variant with mean $\mu = 0$, i.e., assuming unbiased mean forecasts. Its performance (reported in Table S2 in the online supplement) is very similar to that of the unrestricted model.
In order to investigate whether the methods' performance varies across forecast horizons $h$, Figure S2 in the online supplement plots their coverage rate against $h$. As shown by the figure, the Gaussian, decomposition and combination methods attain good coverage rates (close to 80%) across all horizons. The flexible method's coverage rates tend to be slightly too low for the longest horizons.
We next compare the error-based forecasting methods to the SPF's histogram-type forecasts that are available for the two US data sets. As shown in Table (ref), the histograms' coverage rate exceeds its nominal level of 80%, which by itself is desirable. However, this coverage rate comes at the expense of very wide prediction intervals. This applies especially to inflation, where the average length of the histogram intervals exceeds the average length of the other methods' intervals by about 75%. Furthermore, the histograms attain markedly higher interval scores than the other methods.
In order to investigate whether these performance differences are statistically significant, we conduct DieboldMariano1995 type tests, using the combination method as a natural representative of the postprocessing methods. As noted in Section (ref), the dependence structure of forecast errors is potentially complex when pooling the results across forecast horizons. Similar complexity is to be expected for DieboldMariano1995 type tests, which are based on functions of forecasts and realizations (specifically, interval score differences). In order to reduce complexity, and arrive at the standard time series setup considered by DieboldMariano1995, we compare the combination to the SPF histograms separately for each forecast horizon. This choice comes at the cost of a smaller evaluation sample of a single forecast/observation pair per year for each horizon. Due to some variation in data availability across horizons, the size of the evaluation samples then ranges from $36$ to $42$. The test statistics are based on autocorrelation-consistent standard errors as implemented in the function NeweyWest of the R package sandwich Zeileis2004,ZeileisEtAl2020. Figure (ref) summarizes the results. For GDP, the combination significantly outperforms the SPF at the three shortest horizons, using a 5% significance level and a conservative (Bonferroni) $p$-value adjustment for multiple testing. For inflation, we similarly observe significant outperformance at four of the five shortest horizons. For the other horizons, the performance difference between the combination and the SPF is insignificant. In Section B.2 of the online supplement, we provide additional discussion and results on DieboldMariano1995 type testing in the current setup, including three additional implementation variants for computing the test statistic. As described in the supplement, these variants differ in their treatment of potential autocorrelation in the time series of interval score differences. While some of the variants yield considerably higher $p$-values for GDP at the four shortest horizons, the results for inflation are similar to the ones in Figure (ref).
Figure (ref) further shows the length of the combination and SPF methods' prediction intervals separately for each horizon.\footnote{Each dot or triangle in the figure represents a single forecast case. For the combination method (represented by triangles), the intervals are quite similar in length across forecast cases; this is because of highly overlapping training samples resulting from the employed cross-validation procedure. The variation across forecast cases is hence dwarfed by the amount of variation in the SPF, and is often invisible in the figure.} The SPF's intervals are clearly wider than the combination's, especially at short forecast horizons. Findings from the literature on combining forecast distributions GneitingRanjan2013,Lichtendahl2013 suggest that the width of the SPF's prediction intervals may partly be driven by our use of linear combination, in that we consider the mean of the SPF participants' probability forecasts.\footnote{Note, however, that we fit a continuous distribution to these mean probability forecasts in order to compute quantiles. This extra fitting step is not covered by the theoretical literature on combining forecast distributions.} It is thus natural to ask whether nonlinear combination methods are more successful in the context of the SPF. To this end, we consider two quantile-based combination approaches: First, the simple average of the individual participants' quantile forecasts, as introduced in Section (ref) in the context of combining postprocessing methods. Second, the median of the individual quantile forecasts BracherEtAl2021. Both approaches require us to first fit a parametric distribution to each SPF participant's forecast distribution, for which we again follow KruegerPavlova2022. As shown in Table S2 in the online supplement, the quantile-based combination methods indeed yield shorter prediction intervals. However, this comes at the cost of a lower coverage rate, so that they attain similar interval scores as the linear SPF combination, and thus perform worse than the postprocessing methods.
The findings discussed in the previous two paragraphs are closely in line with Clements2010,Clements2014 who documents that the SPF's fixed-event forecast distributions are implausibly wide at short horizons, perhaps reflecting incoherent updating behavior on the part of forecasters. On the whole, our results thus indicate that even fairly simple models based on past forecast errors are more accurate than the forecasters' own assessment of uncertainty.
Figure (ref) illustrates the Gaussian and combination methods for making prediction intervals, for German GDP (top panel), US GDP (middle) and US inflation (bottom). The figure refers to the cross-validation run that excludes the 2020 data from model fitting. Thus, the figure contains both test-sample observations (for 2020, represented by black triangles) and training-sample observations (other years, represented by grey dots). Partly by construction, the clear majority of historical forecast errors is within the methods' prediction intervals.\footnote{For the training-sample observations, close-to-nominal coverage is encouraged by the criteria used for fitting the models. Hence the models' good coverage properties for these observations arises `partly by construction'. By contrast, note that the coverage rates reported in Table (ref) are computed from test-sample observations exclusively, and hence do not arise by construction.} For the Gaussian method (solid line), the intervals are symmetric around zero. For the combination method (dashed line), the intervals are somewhat asymmetric, especially for German GDP. The figure also shows that the 2020 forecast errors for German and US GDP are very large by the standards of the data set. This applies in particular to the negative forecast errors at horizon $41$ or larger (corresponding to forecasts made before mid March 2020), and can be attributed to the effects of the Covid-19 pandemic. By contrast, some of the 2020 forecast errors at horizons 20-30 are large and positive, indicating that forecasts made around mid-2020 were too pessimistic. It seems unsurprising that the methods' prediction intervals -- which are designed to attain an 80% coverage level -- do not capture forecast errors in an excessively turbulent setup like 2020 GDP. For US inflation (bottom panel), 2020 forecast errors are mostly covered by the methods' prediction intervals.
To illustrate the quantitative interpretation of our results, consider the decomposition method as applied to German GDP, and a hypothetical forecast made in mid-September. The 80% prediction intervals for the current-year forecast (corresponding to horizon $15$) are $0.68$ percentage points wide, and a point forecast of $m \%$ GDP growth translates into a prediction interval of $[(m-0.34)\%, (m+0.34)\%]$. The prediction intervals for the next-year forecast (corresponding to horizon $67$) are more than five times as wide, and are given by $[(m-2.01)\%, (m+2.01)\%]$. In forecast-related press statement or media reports, point forecasts for the current and next are often mentioned next to each other. Our numerical example highlights that this practice is possibly misleading, in that forecast uncertainty may differ substantially across the two horizons.
This paper argues that economic fixed-event forecasts should be accompanied by a numerical measure of uncertainty, and proposes methods for computing such a measure. We conclude by discussing relations to the pertinent literature, as well as practical implications.
Transforming fixed-event forecasts. A number of studies consider transforming fixed-event forecasts into fixed-horizon forecasts, on the grounds that the latter are more convenient and more flexible from a statistical perspective. KnueppelVladu2016 consider transforming fixed-event point forecasts, whereas GanicsEtAl2020 seek to transform the SPF's fixed-event forecast distributions. ClarkEtAl2022 use entropic tilting to incorporate fixed-event information from the SPF (both point forecasts and distributions) into a statistical model that can generate flexible types of forecasts. In contrast to these studies, we estimate the uncertainty of fixed-event point forecasts as a function of their horizon $h$.
Pooling versus not pooling of forecast error data. Our methods are based on pooling past forecast errors across horizons $h$, and fitting the distribution of forecast errors as a function of $h$. This approach is motivated by a small sample of forecast errors for each specific horizon $h$, especially for the German data set we consider. In more data-rich situations, using entirely separate statistical models for each forecast horizon may be preferable. DieboldGoebel2022 follow this alternative approach for constructing fixed-event predictions of arctic sea ice extent. Furthermore, we assume that the distribution of forecast errors depends on $h$ alone, which could be relaxed in more data-rich situations. For example, allowing for heterogeneity across target years (by means of random effects type models, say) is conceivable.
Communication of forecast uncertainty. Many economists agree that forecast uncertainty should be measured and communicated. According to ReifschneiderTulip2019, for example, the fact that `[..] prediction errors -- even on occasion quite large ones -- are a normal part of the process [..]' should be communicated to the public in order to enhance the credibility of future forecasts. Perhaps most prominently, many central banks publish forecast distributions (`fan charts') pioneered by the Bank of England in the 1990s GalbraithNorden2012. That said, most economic point forecasts are not accompanied by a numerical measure of uncertainty. Instead, media reports and publications by forecasting institutions often contain a verbal disclaimer that mentions forecast uncertainty. For example, the European Commission's Summer 2022 forecast lists various `risks to the outlook' EuropeanCommission2022, representing possible economic sources of forecast error. In our view, verbal statements are not a satisfactory alternative to a numerical measure in the present context. Among other problems, words such as `likely' often mean different things to different people, and are hard to falsify ex-post DhamiMandel2022. Communicating forecast uncertainty in a practical yet precise way thus remains an important challenge in economic policy. Prediction intervals are an attractive format for doing so since they are reasonably simple (consisting of two numbers only), can be represented in graphical form, and are easy to evaluate ex-post Raftery2016.
\setcounter{figure}{0} \setcounter{table}{0}