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.
86,680 characters · 20 sections · 90 citation commands
Solving the Forecast Combination Puzzle
Since their inception (Stone1961 and Bates1969), forecast combination methods have garnered a dedicated following due to their flexibility and accuracy (Timmermann2006; Aastveit2019). Such methods also align with the zeitgeist of modern econometric thought, in that they are designed to accommodate the fact that not all data sources are created equal, and that the models we are working with in economics are at best an approximation to reality.
Generally, forecast combinations are constructed by producing forecasts for individual, or constituent models, and then combining them via some combination function or weighting scheme. Point forecast combinations are typically constructed by taking a weighted average of point forecasts produced by the constituent models (Bates1969; Stock2004; Timmermann2006; Smith2009; Claeskens2016). In the case of distributional forecast combinations, two commonly used approaches are the linear opinion pool ( Stone1961; Hall2007; Geweke2011; Opschoor2017; Martin2021), and the beta-transformed linear opinion pool (Ranjan2010; Gneiting2013; Satopaeae2014; Baran2018). The weighted average, linear pool and beta-transformed linear pool are all combination functions that map a set of forecasts, produced using constituent models, to a single forecast combination.
Even though the last fifty years has seen these methods rise to prominence among empirical forecasters (Makridakis2018; Thorey2018; Wang2018; Makridakis2020; Taylor2020), key issues regarding the use and abuse of the methods remain. One of the most interesting issues is the so-called `forecast combination puzzle', which is a stylized fact that states that predictions produced using complicated combinations of different forecasts, e.g., via optimizing the weights in the combination using some criterion function that encapsulates some aspect of forecast accuracy, do not generally outperform simpler procedures (such as equal-weighted combinations). {For example, see Stock2004, Smith2009, Makridakis2018, Makridakis2020 for empirical evidence of this phenomenon.}
Explanations for the puzzle range from the increased sampling variability of complex weighting schemes (Stock2004; Claeskens2016), to the similar performance of equally-weighted and optimally-weighted combinations (Elliott2011), and to bias in the averages loss functions (Chan2018). However, all of the above explanations are specific to linear combinations of point forecasts evaluated according to mean squared forecast error. To date, there is no single universally accepted answer that sufficiently explains this puzzle across both point and distributional forecasts, or across different performance measures. See graefe2014combining for a historical summary of the puzzle.
Herein, we analyze the forecast combination puzzle in general terms, and through the lens of tests of superior forecast accuracy (White2000 ; Hansen2005). We demonstrate that the standard two-step approach to producing forecast combinations results in tests that have no power against a large class of alternatives, including fixed-, random-, and drifting-weighting schemes. This result is a consequence of the following unintuitive feature of forecast combinations: when produced in the standard manner, the {usual test statistic employed to gauge differences in forecasts is such that estimated combination weights imparts no sampling variability into the statistic at first-order}. This finding extends and renders rigorous previous results documented in a variety of contexts, including by Stock2004, and Smith2009. In particular, Smith2009 argue that \textquotedblleft the forecast combination puzzle rests on a gain ... that has no practical significance\textquotedblright , and our results make rigorous this statement by showing that, in general settings, tests of superior forecast accuracy cannot distinguish between estimated and a large class of forecast combinations, including the equally-weighted combination.
More generally, we show that in order for the difference between two (sets of) forecast combinations to be meaningful, the combination weights must be surprisingly disparate. Critically, the distance between two combination weights needed to yield a test with non-trivial power depends entirely on the chosen loss and the variability in the constituent model forecasts. Therefore, when comparing two forecast combinations obtained using the same constituent models but different combination weights, the variability in the constituent forecasting models can easily swamp large differences in the forecast distribution due to differences in the combination weights. Consequently, a statistically significant difference between two forecast combinations is unlikely to eventuate in practice unless we have: 1) large sample sizes, much larger than typically used in empirical applications; and 2) constituent model forecasts that have low variability, which ultimately requires that the estimators of the underlying parameters of the constituent models also have low variability. Consequently, in empirical applications that require the use of high-variance constituent model forecasts, it is unlikely that we will be able to detect differences amongst different forecast combination methods.
We demonstrate that the poor behavior of forecast accuracy tests in this setting results from the test statistic having a non-standard asymptotic distribution under the null of no inferior predictive accuracy. {That is, w} hile the critical values for such tests are {typically }based on the standard normal distribution, we show that in the case of forecast combination methods, an appropriately scaled version of the test statistic converges in distribution to a generalized chi-squared distribution. Furthermore, we show that this result persists across a large class of forecast combination methods, including those with time-varying weights.
Lastly, we demonstrate that the forecast combination puzzle can be circumvented in cases where it is feasible to produce forecast combinations in a single step. That is, except in trivial cases, forecast combinations produced in a single step do not exhibit the forecast combination puzzle, and will always yield superior forecast accuracy over standard, or equally-weighted, forecast combination schemes. In this way, we build on an extensive investigation of one- and two-step forecast combinations in Zischke2022, which also provides support for the one-step approach. We reiterate that it is the two-step approach to producing forecast combinations that is the standard approach adopted in the literature on (frequentist) forecast combinations and, hence, the reason why the combination puzzle has been so empirically prevalent. Revisiting the S&P500 returns density prediction example considered in Geweke2011, we show that: 1) the forecast combination puzzle is evident in their original example for certain classes of volatility models; {and} 2) the puzzle is entirely resolved if we use forecast combinations that are built in a single step.
Before moving on, we note here the following notations used throughout the remainder of the paper. For a probability measure $P$, a random variable $X$ and a random sequence $X_n$, we write $\mathbb{E}_P[X]$ to denote the expectation of $X$ under $P$, $\operatornamewithlimits{plim\,}_n X_n$ to denote the probability limit of $X_n$ as $n \to \infty$ (if it exists), and $ X_n \Rightarrow X$ if $X_n$ converges in distribution to $X$. For some positive sequence $R_n$, the notations $X_n = o_p(R_n)$ and $X_n = O_p(R_n)$ have their usual definitions, see VanDerVaart1998 for a textbook treatment. We say that $X_n \asymp R_n$ if there exists constants $c$ and $C$ such that $c R_n \leq X_n \leq C R_n$ for all $n$ large enough (with probability one). The gradient and hessian of a functional $f$ of $x$ is written $\nabla_{x} f(x)$ and $\nabla^2_{xx} f(x)$, respectively, where for $ x$ on the boundary of the domain of $f$ these symbols denote the left or right derivatives, whichever of those exist at $x$.
We first motivate our analysis by reconsidering the results and discussion in Smith2009 and demonstrating that the forecast combination puzzle extends far beyond their initial analysis.
Our goal is to produce a point forecast for the random variable $Y_{t}$ at time $t=T+1$ using observed data $\{y_t:1\le t\le T\}$ and models $f_{j}$, $ j=1,2$. Denote the point predictions from these models by $\tilde{y}_{jt}$, and define the linear combination of point predictions by $\tilde{y} _{t}^{\eta}=\eta \tilde{y}_{1t}+(1-\eta)\tilde{y}_{2t}$ where $0\le\eta\le 1$ , and assume that the parameters underlying the models $f_1$ and $f_2$ have been estimated in a first stage. The accuracy of the forecasts is measured according to mean squared forecast error (MSFE).
Under the assumptions in Smith2009, e.g., $\tilde{y}_{jt}$ is unbiased with constant variance, $\sigma^2_j$, and covariance with $\tilde{y} _{it}$, $\sigma_{ji}, j,i = 1,2$, the optimal forecast combination weight $ \eta$ is obtained by solving $\min_{\eta\in[0,1]}\mathbb{E}\left\{y_t-\tilde{ y}_t^\eta\right\}^2$ and is given by $\eta^\star=[\sigma^2_2-\sigma_{12}]/( \sigma^2_1+\sigma^2_2-2\sigma_{12}). $
Following notational conventions (see, e.g., West1996): based on a sample of size $T+1$, we split this sample into $R$ in-sample observations used for model fitting, and $P$ out-of-sample observations used for forecast evaluation, so that $R+P=T+1$. Under this evaluation regime, Smith2009 use the out-of-sample MSFE,
to compare two different forecast combinations: 1) the equally-weighted combination, with $\eta=1/2$; 2) the sample estimate, $\tilde\eta$, of the optimal combination, $\eta^\star$. The key finding of Smith2009 is that, under their assumptions, when $\eta^\star=1/2$, i.e. when the optimal approach (`in population') is to actually equally weight the two forecasts,
That is, when the optimal combination weight is $\eta^\star=1/2$ the additional finite-sample noise introduced via the estimation of $\eta$ induces an additional penalty due to estimation error which, in turn, results in the fixed-weight scheme displaying superior performance. This finding leads Smith2009 to conclude that “The parameter estimation effect [of the weights] is not large, nevertheless it explains the forecast combination puzzle.”
While the analysis of Smith2009 is insightful, their findings are not immediately generalizable to related situations. For instance, the analysis is based on linear combinations of point forecasts, with unbiased constituent forecasts that do not depend on unknown parameters.\footnote{ It is perhaps more accurate to say that the impact of having to estimate unknown parameters in the constituent models does not feature in their analysis.} In addition, and most importantly, the optimal MSFE weight is assumed to be $\eta^\star=1/2$, which coincides with the default (equally-weighted) combination that typically underpins the forecast combination puzzle. Finally, this analysis does not immediately extend to other loss functions.
However, there is now mounting evidence to suggest that the forecast combination puzzle extends beyond this stylized setup. To this end, we explore scenarios in which we allow for a range of values for $\eta^{\star}$ , and a range of values for the fixed weight $\eta$ that differ from $\eta^{\star}$, for the case of both the point forecast combination of Smith2009 and a distributional forecast combination, and where these forecast combinations are estimated using different loss functions.
Following Section 3.1 of Smith2009, consider that the true DGP is from the AR(2) family,
and that our forecasts are based on a linear pool $f^{(t)}$ of two normal constituent distributional forecasts $f^{(t)}_1$ and $f^{(t)}_2$:
where $N\{x;\mu,\Sigma\}$ denotes the normal pdf evaluated at $x$ with mean $ \mu$ and variance $\Sigma$, $\gamma_1$ and $\gamma_2$ are the parameters of the constituent models, and $\eta$ is the weight assigned to the first model.
We will estimate the parameters of two forecast combinations. First, we consider the distributional forecast combination given by the linear pool above, and estimate the parameters ($\eta$, $\gamma_1$, $\gamma_2$) by minimizing the log loss (or equivalently by maximizing the log likelihood). Second, we consider the point forecast combination given by the expectation of the distributional forecast combination, and estimate the parameters by minimizing the MSFE. The second combination is identical to the combination considered by Smith2009 and discussed in the previous section. Parameters are estimated in the standard two-step fashion, so that $ \gamma_1$ and $\gamma_2$ are chosen to minimize the selected loss of the first and second constituent models, respectively, and $\eta$ is estimated by minimizing the loss of the combination given the aforementioned estimates for $\gamma_1$ and $\gamma_2$. Section (ref) discusses such optimal forecast combinations in more detail.
For a given (point or distributional) combination, the parameters $ \phi_1,\phi_2$ and $\sigma^2_\epsilon$ can be chosen so that a desired value of $\eta^{\star}$ is achieved (we refer to Appendix (ref) for details). We then test the null of no inferior forecast accuracy across a variety of benchmark forecasts constructed using a range of fixed weights $\eta \in \{0.25, 0.5, 0.75\}$ (which includes the equally-weighted benchmark $\eta = 0.5$) against the alternative optimally-weighted combination.
Figure (ref) presents the rejection frequency ($y$-axis) of the test across each fixed weight (rows) in the case of the MSFE (for the point forecast combination, left-hand column) and the log loss (for the distributional forecast combination, right-hand column). We present these results for DGP parameter values corresponding to a variety of pseudo-true weights $\eta^{\star} \in \{0, 0.25, 0.5, 0.75, 1\}$ (colors) and across a grid of values for the sample size $T+1$ ($x$-axis), comprising the $R = (T+1)/2$ in-sample observations followed by $P = (T+1)/2$ out-of-sample observations. For the sake of brevity, we leave a more detailed discussion of the specific implementation details for this exercise to Appendix (ref).
Analyzing each figure we see that, across both loss functions, the probability of rejection under the null hypothesis of no inferior forecast accuracy of the benchmark (that is, where the $\eta ^{\star }$ value of the color equals the $\eta $ value of the row) is (virtually) zero, and certainly less than 0.05, the nominal size of the test. Analyzing the rejection rates {in all panels, which are given by the cases where the $\eta $ values of the color are different from the $\eta ^{\star }$ value of the row,} suggests that the forecast combination puzzle persists even when the distance $\lvert \eta ^{\star }-\eta \rvert $ between $\eta ^{\star }$ and the incorrect fixed weight $\eta $ is large. For example, a test of no inferior predictive accuracy of the equally-weighted combination (middle row) against the optimally-weighted combination where the truth is $\eta ^{\star }=0.25$ (dark blue) has a rejection frequency smaller than 50% across all sample sizes less than 1000, whether we are using the MSFE (left-hand panel) or the log loss (right-hand panel). Hence, the power of such tests may be quite low in practice even if there are meaningful differences between forecasts. Later, in Section (ref), we will conduct a similar exercise to illustrate how this issue of low power can be resolved by producing optimally-weighted forecast combinations in a one-step fashion. The production of one-step combinations is described in the next section.
Now using more formal notation, let $Y_1, \cdots, Y_{T},\dots$ denote a sequence of random variables generated from the probability triple $(\mathsf{ Y}, \mathscr{A}, G)$, where $\mathsf{Y}\subseteq\mathbb{R}^d$, $d\ge1$. Since the true measure $G$ is unknown in general, we postulate a class of models $\mathcal{Q}$ on $\mathsf{Y}$, which we identify by their distribution functions $Q\in\mathcal{Q}$.
Given an observed sample, $y_1,\dots,y_T$, our goal is to predict some feature of $Y_{T+h}$ at a forecast horizon of $h\ge1$. Since the true model is unknown, there is no single direction of truth with which to predict features of interest. The literature has generally maintained that in such cases entertaining multiple models is a valid approach that can produce reliable predictions, with the most common approach being to produce forecast combinations, as we have highlighted.
Denote by $F_{\gamma_1}:\mathsf{Y}\times\Gamma_1\rightarrow C[0,1]$ a probability measure on $(\mathsf{Y},\mathscr{A})$ indexed by the parameter $ \gamma_1\in\Gamma_1\subseteq\mathbb{R}^{d_{\gamma_1}}$ and itself lying in the family $\mathcal{F}(\Gamma_1):=\{F_{\gamma_1}:\gamma_1\in\Gamma_1\}$. For any $1\le n\le T$, let $\Omega_{n}$ denote the information available to the forecaster at time $n$, and denote the predictive measure based on time-$ n$ information as $F_{\gamma_1}^{(n)}:=F_{\gamma_1}(\cdot:\Omega_n)$.
In the majority of forecasting settings, the practitioner entertains a collection of $K<T$ possible statistical models that can each describe, with varying accuracy, the movements of the stochastic process $\{Y_t:t\le T\}$. We consider that each model is specified using a (semi-) parametric family $ \mathcal{F}(\Gamma_j)$, which depends on $\Gamma_j\subseteq\mathbb{R} ^{d_{\gamma_j}}$ unknown parameters for each $j=1,\dots,K$. For
denoting the collection of all $K$ constituent models, we can combine these models to produce a forecast combination. To this end, and following Gneiting2013, we consider the combination function
Common choices for the combination family $\mathcal{C}:=\{C_\eta:\eta\in \mathcal{E}\}$ include the linear pool (see, e.g., Geweke2011), and the Beta-transformed linear pool (Gneiting2013).
Given the family of combination functions, $\mathcal{C}$, and the member model family $\mathcal{M}$, the class of probability measures used for prediction is the composition of the two: define $\gamma:=(\gamma_1^{\prime },\dots,\gamma_K^{\prime })^{\prime }$, $\theta:=(\eta^{\prime },\gamma^{\prime })^{\prime }$, $\Theta:=\mathcal{E}\times\Gamma_1\times \dots\times\Gamma_K$ and consider $Q_\theta:\mathsf{Y}\times\mathcal{C}\times \mathcal{M}\rightarrow C[0,1]$ defined by $Q_\theta=C_\eta\circ M_\gamma$. We denote by $\mathcal{Q}$ the class $\{Q_\theta:\theta\in\Theta\}$.
Generally, the parameters of $Q_\theta$ are unknown and must be estimated. Throughout, we consider that the forecaster wishes to obtain `optimal' forecasts in the spirit of Gneiting2007, Gneiting2011, and Martin2021. Herein, we take optimal to mean that the distributions we choose to take out-of-sample are produced by targeting a loss function that measures precisely the features of $Y_{T+h}$ that are of interest. Following Gneiting2011, we consider a decision-theoretic framework for such forecasts. Let $\mathsf{Y}\subseteq\mathbb{R}^d$ denote the observation domain, and for some $k\ge1$ let $\mathsf{A}\subseteq\mathbb{R} ^k$ denote the action space. We consider two mutually exclusive cases: the forecaster is interested in measuring predictive accuracy for a functional of the distribution of $Y_{T+h}$, or the entire distribution.
Recall that $\mathcal{Q}$ is a class of distributions on $\mathsf{Y}$, and consider a functional $U:\mathcal{Q}\mapsto\mathsf{A},\quad Q\mapsto U[Q]\subseteq \mathsf{A}$ that maps a distribution $Q\in\mathcal{Q}$ to a subset $U[Q]$ of the action space. A scoring function is a measurable map $S:\mathsf{A}\times \mathsf{Y}\rightarrow[0,\infty)$. We orient the scoring functions so that a lower score implies a more accurate forecast. The scoring function $S(\cdot,\cdot)$ is $\mathcal{Q}$ -consistent for a functional $U[\cdot]$ if
where $Y$ is a random variable with distribution $Q$, and here and throughout we assume that any stated expectation exists and is finite. We say that $S(\cdot,\cdot)$ is $\mathcal{Q}$-strictly consistent for $T[\cdot]$ if
If $U[\cdot]$ admits a strictly consistent scoring function, then it is called elicitable. Here and throughout, we consider that the functional we are interested in is elicitable. For several functionals such as mean, quantiles and expectiles, this is the case, but there are many functionals, such as variance or expected shortfall, that are not elicitable, {at least on their own. Some of these functionals are elicitable jointly with others, which is the case for the pairs (mean, variance) and (value at risk, expected shortfall), for example ( Fissler2016)}.
Assume the forecaster is interested in the case where $U[\cdot]$ is a particular functional of the distribution of the random variable $Y_{T+h}$, and we have available information $\Omega_T$. Then, we are interested in generating predictions at time $T+h$, $h\ge1$, for the functional $ U[F^{(T)}_{Y_{T+h}}]$, where we recall that $F_{Y_{T+h}}^{(T)}$ signifies the distribution of the random variable $Y_{T+h}$ conditional on information $\Omega_T$. To this end, we follow, among others, Patton2020 and assume that any particular combination model $Q_\theta\in\mathcal{Q}$ admits a model for $U[F_{Y_{T+h}}^{(T)}]$ of the form $m(Z_T;\theta)=U[Q^{(T)}_ \theta]$, where $Z_T\in\Omega_T$ denotes observable variables in the conditioning set and $m:\Omega_T\times\Theta\rightarrow\mathsf{A}$ is known up to the unknown $\theta$. We can then produce point forecasts for $ U[F_{Y_{T+h}}^{(T)}]$ using $m(Z_T;\theta)$, and by replacing the unknown $ \theta$ with
The above approach allows one to produce distributions that generate reliable `point forecasts' for a given functional {, but there is no reason for these distributions to be accurate in any other respect. In cases where we want distributional forecasts that are accurate as complete representations of the uncertainty surrounding an unobserved random variable, we can estimate our parameters according a proper scoring rule.}
A negatively-orientated proper scoring rule is a function $S: \mathcal{Q}\times\mathsf{Y}\mapsto\mathbb{R}$ such that
for all $Q,G\in\mathcal{Q}$. A strictly proper scoring rule is a proper scoring rule that is minimized by $G$ alone.\footnote{ Scoring rules have a deep connection to decision theory, and we do not review this literature here (see, e.g., Pesaran2002 and Granger2006 for a discussion in the context of economics and financial forecasting).} While scoring rules can be used to measure the accuracy of the predictive distribution or densities, they can also be used to measure the accuracy of certain features of the distribution: e.g., quantiles, or predictive intervals (see, Gneiting2007). Throughout the remainder, when we refer to scoring rules, it is meant that the action space is either the class of densities/distributions and not some functional of it such as quantiles or intervals; we keep this distinction as accuracy for such latter quantities can be readily measured using (consistent) scoring functions.
Similar to the case of scoring functions, producing density forecasts requires estimating $\theta$ in the combination model $Q_\theta\in\mathcal{Q} $. In this case, we can define an estimator of this parameter by minimizing the expected scoring rule over our observed sample
Regardless of whether one is producing `point' forecasts of a functional or distributional forecasts, producing the `optimal predictive' is categorically the same. Therefore, to simplify the presentation, we jointly treat both approaches. Given a collection of realized observations $\{y_t: 1\le t\le n\}$, $n\le T$, we search for the most accurate combination predictive distribution $Q_\theta\in\mathcal{Q}$ under the given loss function $L:\mathcal{Q}\times\mathsf{Y}\rightarrow \mathbb{R}$.
The optimal predictive can generally be produced in two possible ways. In cases where it is feasible to estimate the parameters jointly, the combination predictive can be produced by estimating the unknown model parameters,
and the predictive distribution at time $T+h$ can be taken as $ Q^{(T)}_{\widehat\theta_n}$. However, due to the dimensionality of $\theta$, it is often difficult to estimate $\theta$ jointly. Generally, forecast combinations are estimated in two steps: first, the constituent model forecasts are produced, and then the combination weights are estimated conditional on the constituent model forecasts; see, e.g., Hall2007, Geweke2011, Gneiting2013. Throughout we assume that such two-step forecast combinations are carried out in the spirit proposed in Gneiting2007. Namely, for each constituent model we estimate $\gamma_j$, $j = 1, 2, \cdots, K, $ via the estimator
Collecting the $\tilde\gamma_{jn}$ into $\tilde\gamma_n:=(\tilde \gamma_{1n}^{\prime },\dots,\tilde\gamma_{Kn}^{\prime })^{\prime }$, we can then estimate the combination parameters $\eta$ by maximizing $ Q_\theta\mapsto L_n(\theta)$, conditional on $\gamma=\tilde\gamma_n$, which yields
Once $\tilde\theta_n$ have been calculated the predictive distribution $ Q^{(T)}_{\tilde\theta_n}$ can be used to produce predictions at time $T+h$.
Throughout the remainder, we refer to the predictive distributions $ Q^{(T)}_{\widehat\theta_n}$ and $Q^{(T)}_{\tilde\theta_n}$ as the one-step and two-step predictive combinations, respectively.
Under standard regularity conditions, the extremum estimators $\widehat{ \theta }_{n}$ and $\tilde{\theta}_{n}$ will converge to well-defined probability limits, which we denote by $\theta ^{0}:=(\eta ^{0^{\prime }},\gamma ^{0^{\prime }})^{\prime }$, and $\theta ^{\star }:=(\eta ^{\star ^{\prime }},\gamma ^{\star ^{\prime }})^{\prime }$, respectively;\footnote{ For more precise definitions of $\theta ^{0}$ and $\theta ^{\star }$, we refer the interested reader to Appendix (ref).} i.e., under our regularity conditions we will have
Further, it is well-known, see Newey1994 for details, that one-step and two-step estimators do not coincide in general, so that throughout we can assume that $\theta ^{0}\neq \theta ^{\star }.$ If there is no asymptotic difference between the one- and two-step estimators, there will be no asymptotic difference between the resulting forecast combinations based on these estimators. In most empirical settings, however, significant differences between one- and two-step estimators exist, and so we restrict our attention to this case in following analysis.
For a general discussion on two-step estimation see Newey1994; for a more modern treatment see Frazier2017; and for a particular approach to measuring the impact of two-step estimation in the context of distributional forecast combinations see Zischke2022.
Following West1996, White2000, and many others, we measure out-of-sample predictive accuracy using loss differences of forecasts over a given out-of-sample period. For simplicity, and for consistency with the earlier exposition of the analysis of Smith2009, we let $h=1$ denote the horizon over which we will make predictions, but note that our results can also accommodate $h\ge1$ at the cost of additional notation. For ease of exposition, we re-introduce the notation introduced in Section (ref), and assume the sample consists of $T+1$ total observations, which we partition into $R$ in-sample periods, and $P$ out-of-sample periods, across which we evaluate the predictions, where $ R+P=T+1$.
In the approach of West1996, the information used to estimate $\theta$ is increased by one unit for each prediction; i.e., at time $t=R+1$, $R$ observations are used to estimate $\theta$, at time $t=R+2$, $R+1$ observations are used, and so forth. This formulation is useful as it allows one to update the parameter estimates as new information becomes available. However, the resulting out-of-sample loss difference at time $T+1$ is then a complex combination of all previous estimators, and also has variability due to the $P$ out-of-sample observations themselves. Consequently, disentangling variability due to parameter uncertainty, from the innate variability of the average loss difference becomes difficult, and obtaining clear intuition regarding the contribution to each of these pieces to the behavior of the out-of-sample average loss difference becomes difficult.
Furthermore, it is critical for us to understand the precise impact of parameter uncertainty on forecast accuracy since this effect, while not large, “explains the puzzle” according to Smith2009. Therefore, we consider a simple framework that cleanly dissects the two types of sampling variability, loss and parameter estimation. This is accomplished by estimating the unknown parameters once using $R$ observations, with the resulting estimators then held fixed over the $P$ out-of-sample periods; \footnote{ More specifically, instead of the first prediction being based on an estimator $\theta_{R}$, obtained using observations $1,\dots,R$, and the next based on $\theta_{R+1}$, and so forth, we only consider estimators based on $R$ observations (with $R\rightarrow\infty$ as $T\rightarrow\infty$ ).} the $P$ out-of-sample periods are then used for evaluation only; further, we will also maintain that the in-sample period, $R$, and the out-of-sample period, $P$, are in rough proportion.\footnote{ We refer to Section 8 of clark2013advances for a discussion on the benefits and disadvantages of various splitting schemes for forecast evaluation.}
Following White2000 and Hansen2005, we measure the accuracy of forecasts by testing the null hypothesis of no inferior forecast performance using the average loss difference over the $P$ out-of-sample periods. Consider that we wish to test the accuracy of a benchmark forecast distribution $Q_{\vartheta ^{b}}^{(T)}$, indexed by unknown parameters $ \vartheta ^{b}$, against an alternative distribution $Q_{\theta ^{a}}^{(T)}$ , indexed by unknown parameters $\theta ^{a}$. The null hypothesis of no inferior forecast accuracy of the benchmark ($Q_{\vartheta ^{b}}^{(T)}$) over the alternative ($Q_{\theta ^{a}}^{(T)}$) is
To test the null in (ref), we approximate the above expectation, and the unknown $\vartheta^b,\theta^a$, using their sample counterparts. For any $t> R$, and any sequence of consistent estimators $\vartheta_R,\theta_R \in\Theta$ of $\vartheta^b, \theta^a$, respectively, define the average loss difference statistic
where we recall that $\ell_{t}(\theta_{R})= L[Q^{(t-1)}_{\theta_{R}},y_{t}]$ . The null hypothesis in (ref) can then be tested using the standardized statistic
where $\hat{\Omega}_R$ is a consistent estimator of the asymptotic variance of $\Delta_P(\vartheta^b,\theta^a)$.
Our results in Section (ref) demonstrate that even when the optimal combination weight (in population), $\eta ^{\star }$, is very different from a fixed, hypothetical combination weight, the standard approach to testing for differences in forecasting accuracy does not result in meaningful rejection rates; i.e., even though the hypothetical combination weight is inferior to the optimal combination weight, the resulting testing procedure does not reliably detect differences. This finding holds across two particular loss functions and a host of different values of $\eta ^{\star }$. In this section, we demonstrate that this {phenomenon} is present in any class of forecast {combinations produced} in the standard manner, i.e., in two steps. Thus we give, for the first time, a truly generic explanation for the puzzle that is agnostic to the chosen loss, and which is valid under standard regularity conditions.
To state this result, consider the setting where we are given a known, i.e., hypothesised combination weight $\eta _{T}^{\delta }$, which we can always represent as
where $\eta ^{\star }$ {is as defined in ((ref)), and }$\eta ^{\star } $ and $\delta _{T}$ are individually unknown. Consider that the (possibly random) sequence $\{\delta _{T}:T\geq 1\}$ is described by one of the following cases: throughout, let $\delta \in \mathcal{E}\subset \mathbb{R} ^{d_{\eta }}$ be {some non-zero} and bounded vector,
The above class of sequences will allow us to evaluate the behavior of the standard testing framework for accessing accuracy of different forecast combination methods across a wide range of hypothesized combination weights $ \eta _{T}^{\delta }$. In particular, the case $\xi =0$ yields fixed alternatives, while $\xi =1/2$ yields the class of {canonical Pitman} sequences. As we shall see, the behavior of the test depends crucially {on the category within which} $\delta _{T}$ resides.
To this end, consider that our benchmark forecast is $Q_{\vartheta ^{b}}^{(T)}$ {with} $\vartheta ^{b}=(\eta _{T}^{\delta },\gamma ^{\star })$, and we wish to test the null hypothesis that this benchmark forecast has no inferior forecast performance relative to the alternative forecast $Q_{\theta ^{a}}^{(T)}$ {with} $\theta ^{a}=(\eta ^{\star },\gamma ^{\star })$. The null hypothesis in (ref) then becomes
{The null hypothesis} $\text{H}_{0}$ {is then tested using the statistic }in (ref), {where the} infeasible $\theta ^{a}$ is replaced {by the feasible estimator,} $\tilde{\theta}_{R}=(\tilde{\eta}_{R}^{\prime },\tilde{ \gamma}_{R}^{\prime })^{\prime }$,{\ and} the infeasible $\vartheta ^{b}$ by its feasible counterpart $\vartheta _{R}^{\delta }=(\eta _{R}^{\delta }, \tilde{\gamma}_{R}^{\prime })^{\prime }$, {where }$\eta _{R}^{\delta }$ { varies according to the neighbourhood of }$\eta ^{\star }$ {defined \ in (ref). }The rejection region for the test is then defined by the corresponding rejection region
where $\Phi ^{-1}(\alpha )$ denotes the $\alpha $-quantile of the standard normal distribution. The {following theorem describes} the behavior of the resulting test.
The above result implies that if two forecasts have combination weights that are at least $O(T^{-1/4+\varepsilon })$ apart, for $\varepsilon >0$, then the standard testing approach can distinguish between the forecasts. Surprisingly, however, if the combination weights are at a distance of $ CT^{-1/4}$, the test can be arbitrarily over- or under-sized depending on the magnitude of the sequence $\{\delta _{T}:T\geq 1\}$. More surprisingly, if two sets of combination weights are within a distance of $ O(T^{-1/4-\varepsilon })$ from one another, e.g., a parametric neighbourhood of width $O(T^{-1/2})$, then the test has no power to detect differences between the combination forecasts. Note that this third result {also encompasses} the case considered in the illustrative example in Section (ref), in which $\eta _{T}^{\delta }=\eta ^{\star }$ -- i.e. the benchmark fixed weight (denoted by $\eta $ therein) {coincided with} $\eta ^{\star }$-- and the test displayed zero empirical size (up to Monte Carlo error).
As an example of {the phenomenon in Theorem (ref)(iii)}, consider that we have two competing combination forecasts defined by different combination weight schemes such that $\tilde{\eta}_{1}$ and $\tilde{\eta}_{2}$ have distinct asymptotic distributions, but $R^{1/2}(\tilde{\eta}_{1}-\tilde{\eta}_{2})=O_{p}(1)$, then the usual test of the null hypothesis of no inferior predictive accuracy will only detect differences between the two forecasts on very rare occasions, and will detect no statistically significant differences between the forecasts with probability converging to one.\footnote{ We recall that under the maintained assumption on $R,P,T$, we have that $ R\asymp P\asymp T$, so that $O(R^{1/2}/T^{1/2})=O(1)$.}
Practically speaking, Theorem (ref) demonstrates that even if the benchmark forecast, e.g., the equally weighted forecast, is far away from the optimally weighted combination forecast, then the standard testing approach is unlikely to reject the inadequacy of this benchmark (with probability converging to one). In particular, Theorem (ref)(i) demonstrates that the standard test will asymptotically reject the null only when $\sqrt{T}(\eta^\delta_T-\tilde\eta_R)$ diverges faster than $O(T^{1/4})$, i.e., when $\xi<1/4$. Therefore, for all intents and purposes, the standard approach to testing for differences in combination forecasts cannot be trusted to deliver reliable conclusions in the majority of empirical situations where it is applied.
Theorem (ref) demonstrates that when comparing between forecast combinations, tests of forecast accuracy behave in non-standard ways. However, it is important to understand the mechanism causing this behavior. Recall that standard forecast combinations are produced in two steps: first, we estimate the unknown model parameters $\gamma$ via $\tilde\gamma_R$, then the combination weights are estimated via $\tilde\eta_R$. In this section, we show that the two-step nature by which the forecast combination $ Q^{(T)}_\theta$ is produced results in an average loss difference whose limiting distribution is non-standard under a large class of combination weights.
To state the asymptotic distribution of the usual loss-difference test statistic $\Delta_P(\vartheta^b,\theta^a)$, we require the following definitions. Let $\mathcal{L}(\theta):=\operatornamewithlimits{plim\,} _{P\rightarrow\infty} L_P(\theta)/P$, $\mathcal{M}_{\eta\eta}= \nabla_{\eta\eta}\mathcal{L}(\eta^\star,\gamma^\star)$, $\mathcal{M} _{\gamma\eta}= \nabla_{\gamma\eta}\mathcal{L}(\eta^\star,\gamma^\star)$, and let $V_{P,R}:=-\mathcal{M}_{\eta\eta}^{-1/2}\{\nabla_\eta L_P(\eta^\star,\gamma^\star)/P+\mathcal{M}_{\eta\gamma}(\tilde\gamma_R- \gamma^\star)\}$. We note that the above exist under Assumptions (ref)-(ref) in Appendix (ref), and that $V_{P,R}=O_p(1/ \sqrt{P})$.
The expansion in Lemma (ref) clarifies the mechanism behind the behavior exhibited in Theorem (ref). Recall that the behavior of the standard test is driven by the behavior of $\sqrt{P}\cdot \Delta _{P}(\vartheta _{R}^{\delta },\tilde{\theta}_{R})$, see equation (ref). However, Lemma (ref) makes clear that if the sequence $ \{\delta_T:T\ge1\}$ goes to zero fast enough, the limit distribution of $ \sqrt{P}\cdot\Delta _{P}(\vartheta _{R}^{\delta },\tilde{\theta}_{R})$ is degenerate. That is, since $V_{P,R}=O(1/\sqrt{P})$ and $(\tilde\eta_R-\eta^ \star)=O_p(1/\sqrt{P})$ under Assumptions (ref)-(ref), scaling $\Delta _{P}(\vartheta _{R}^{\delta },\tilde{\theta}_{R})$ by $\sqrt{ P}$ results in a degenerate test statistic unless
When $\delta_T=\delta/T^\xi$, $\xi\in(1/4,\infty]$, we have that $ \operatornamewithlimits{plim\,}_P\sqrt{P}\|(\eta^\delta_T-\eta^\star)\|^2=0$ , which yields the result in part (iii) of Theorem (ref). If instead we have $\delta _{T}= T^{-1/4}$, the behavior of $\sqrt{P}\Delta _{P}(\vartheta _{R}^{\delta },\tilde{\theta}_{R})$ is driven by the magnitude of
which yields the second result in Theorem (ref). The above term also drives the power of the statistic $D_P(\vartheta_{R},\tilde\theta_R)$ in the case where $\delta_T\in[0,1/4)$ (i.e., part (i) of Theorem (ref)).\footnote{ We note, however, that Lemma (ref) is not valid, as stated, in the case where $\delta_T\in[0,1/4)$ since the remainder term in (ref) is no longer negligible, since $\sqrt{P}\|\eta^\delta_T-\eta^\star\|^2$ will diverge. Nonetheless, a similar argument to that used to prove Lemma (ref) can be used to deal with this case.}
To obtain the limit distribution of the out-of-sample average loss difference under the null hypothesis in (ref), based on the benchmark forecast combination $Q^{(T)}_{\vartheta^{(b)}}$ with $ \vartheta^{(b)}=(\eta^\delta_T,\gamma^\star)$, and with $\eta^\delta_T$ as in (ref), we require a few additional definitions. Let $ X_{P}:=\nabla_\eta L_P(\eta^\star,\gamma^\star)/P$, $Z_{R,\gamma}:=(\tilde \gamma_P-\gamma^\star)$, so that we can write $V_{P,R}=-\mathcal{M} _{\eta\eta}^{-1/2}(X_P+\mathcal{M}_{\eta\gamma}Z_{R,\gamma})$. Under Assumptions (ref)-(ref) in Appendix (ref), we have that $\sqrt{P} X_{P}\Rightarrow X\sim N(0,\Sigma_X)$, where $\Sigma_X:=\lim_P \text{Var}\{\nabla_\eta L_P(\eta^\star,\gamma^\star)/\sqrt{P}\}$ and $\sqrt{R } Z_{R,\gamma}\Rightarrow Z_\gamma\sim N(0,\Sigma_\gamma)$, and where
Recall that $c:=\lim_T R/P$, with $0<c<\infty$, and by Assumption (ref), for some matrix $Q_V$, we have that $\mathcal{V}:=(\mathcal{V} _1^{\prime },\mathcal{V}_2^{\prime })^{\prime }=\mathcal{L}_{\infty}(c\sqrt{R }V_{R,R},\sqrt{P}V_{P,R})=N(0,Q_V)$, for some matrix $Q_V$, is the limit (joint) law of the terms $c\sqrt{R}V_{R,R}$ and $\sqrt{P}V_{P,R}$.
When $\delta_T\asymp \delta/T^\xi$ and $\xi\in(0,1/2)$, the result of Corollary (ref)(i) demonstrates that under the null hypothesis in (ref) with benchmark forecast combination $Q_{\vartheta ^{b}}^{(T)}$, $\vartheta _{}^{(b)}=(\eta _{T}^{\delta },{\gamma}^{\star})$, and alternative forecast combination $ Q_{\theta ^{a}}^{(T)}$, with $\theta ^{a}=\theta^\star$, the test statistic $ P\cdot \Delta _{P}(\vartheta _{R}^{\delta },\tilde{\theta}_{R})$, where $ \vartheta _{R}^{\delta }=(\eta _{T}^{\delta },\tilde{\gamma}_{R})$, diverges. Hence, if $\eta^\star$ and $\eta^\delta_T$ are sufficiently different, we can accurately learn differences between competing forecast combination methods. In contrast, when $\delta_T\asymp \delta/T^\xi$ , with $\xi\in[1/2,\infty]$, Corollary (ref)(ii)-(iii) demonstrates that the asymptotic distribution of $P\cdot \Delta _{P}(\vartheta _{R}^{\delta },\tilde{\theta}_{R})$ is non-standard.
In the regime where $\delta_T\asymp \delta/T^\xi$, with $\xi\in[1/2,\infty]$ , the test statistic $P\cdot \Delta _{P}(\vartheta _{R}^{\delta },\tilde{ \theta}_{R})$ {converges} in distribution to a random variable with two components. {The first component is a generalized} chi-squared random variable, which does not admit a closed-form formula for its density or distribution function.\footnote{ Since the quadratic form $\Vert X+\mathcal{M}_{\eta \gamma }Z_\gamma\Vert _{ \mathcal{M}_{\eta \eta }^{-1}}^2$ cannot be re-written as a quadratic form with an idempotent weighting matrix, the distribution is not chi-squared.} {The second component is itself possibly comprised of two components:} so long as $\xi\ge1/2$, the second term depends on the difference of two mean-zero but correlated normal {random variables}, which captures the behavior of (a scaled version of) the out-of-sample loss difference due to differences in the combination weights, i.e., the $\eta$-components, and a centering term that captures the difference between $\eta^\delta_T$ and $\eta^\star$.\footnote{ This term results from taking a {second-order} Taylor expansion of the loss difference, and grouping terms appropriately.} In the regime where $\xi=1/2$, an additional term is present that captures the fact that $\sqrt{P}(\eta^\delta_T-\eta_\star)=(\sqrt{P/T})\sqrt{T}(\eta^\delta_T-\eta^\star)=\sqrt{P/(P+R)}\delta$, which, since $c=\lim_T R/P$ converges to $(1/(1+c))^{1/2}\delta$.
However, it is important to note that since the first {component} of the asymptotic distribution is a generalized chi-squared random variable, even if the second {component were} not present it would be infeasible to obtain closed-form quantiles for the null distribution of the test statistic. In addition, it is important to realise that the distribution of this test statistic depends on the loss used in the analysis, the choice of constituent models, and the specific combination function chosen. As such, there is no{ hope for a single set of generally applicable critical values}, and any simulated critical values will need to be application-specific {as a consequence.}
While it is not feasible to obtain closed-form quantiles for the null distribution of the test statistic, since the second component is always negative, it is possible to deduce a conservative test that uses just the critical values of the generalized chi-distribution. That is, since
we can use the quantiles of the generalized chi-distribution to deduce a conservative, but feasible, test of the the null that a benchmark forecast $ Q^{(T)}_{\vartheta^b}$, with $\vartheta^{(b)}=(\eta^\delta_T,\gamma^\star)$, is not inferior to an alternative forecast $Q^{(T)}_{\theta^a}$, with $ \theta^a=\theta^\star$. Such a test will be asymptotically conservative in general, but will not be too conservative so long as the differences between $\sqrt{R} V_{R,R}$ and $\sqrt{P}V_{P,R}$ are small.
To implement such a test, for a {given} benchmark forecast $ Q_{\vartheta^{(b)}}^{(T)}$ based on known combination weight $\eta^\delta_T$ , we must construct a feasible critical value from the generalized chi-distribution $\frac{1}{2}\|X+c^{-1}\mathcal{M}_{\eta\gamma}Z_\gamma\|_{ \mathcal{M}_{\eta\eta}^{-1}}^2$ via simulation. In particular, we start out by first simulating $h=1,\dots,H$ realisations for the random variables $ X^{(h)},Z^{(h)}$ from normal distributions where $\Sigma_X,\Sigma_Z$ are replaced with consistent estimators based on $\eta^\delta_R,\tilde\gamma_R$, and then, for each $h=1,\dots,H$, we form the statistic
where $\widehat{\mathcal{M}}_{\eta\eta}$ and $\widehat{\mathcal{M}} _{\gamma\eta}$ are the usual sample estimate counterparts of the matrices $ \mathcal{M}_{\eta\eta}$ and $\mathcal{M}_{\gamma\eta}$ and calculated at $ \eta^\delta_R,\tilde\gamma_R$. Sorting $\{\Delta^{(h)}:h=1,\dots,H\}$, we can obtain an $\alpha$-level critical value by choosing the $ \lceil(1-\alpha)H \rceil$-th smallest value, and rejecting the null when the observed value of the statistic $P\cdot \Delta_{P}(\vartheta^\delta_R,\tilde\theta_R)$ exceeds this value.
Denoting the above simulated critical value by $\widehat{\mathrm{cv}} _{\lceil(1-\alpha)H\rceil}$, we can define a test of the null hypothesis of no inferior predictive accuracy of the benchmark model $Q^{(T)}_{\eta^ \delta_T,\gamma^\star}$ against the alternative model $Q^{(T)}_{\theta^\star} $ that is estimated in two-steps, via the corresponding rejection region
Unlike the usual test, based on $W_{P}(\alpha)$, the test based on the rejection region ${W}^{2s}_{P}(\alpha )$ explicitly accounts for the two-step nature by which $P\cdot\Delta_{P}(\vartheta _{R}^{\delta },\tilde{ \theta}_{R})$ is constructed. As we have already seen, failure to account for the two-step nature of the forecast combinations does not deliver a test with appropriate size control. Conversely, since a test based on ${W} ^{2s}_{P}(\alpha )$ accounts for the two-step nature of the estimation, the test should deliver appropriate size control.
We note that, while it is feasible to use the above simulation method to construct an appropriate critical value for the testing differences between competing forecast combinations, we do not necessarily advocate for this approach in general. In particular, while the above procedure would deliver an appropriately sized test, it is unclear if the above test is the most powerful approach, and additional research is necessary to determine this. In addition, as we elaborate on in the following Section (ref), in certain cases there is a simpler alternative to {adopting} the above simulated critical value that allows us to entirely avoid the forecast combination puzzle.
{Meanwhile, in this section}, we return to the example setting discussed in Smith2009 and demonstrate {numerically} that the simulated critical {value} suggested in Section (ref) delivers a testing procedure with approximately correct size under the null hypothesis of no inferior forecast accuracy, and also has higher power than the standard approach under the alternative.
To demonstrate {the empirical size of the test}, we return to the example in Section (ref) and consider values of $\phi_1,\phi_2$ and $\sigma^2$ for the AR(2) family such that, under MSFE (respectively, log-loss), the optimal value of the combination weight, $\eta^\star$, obtained under MSFE (respectively, log-loss), is (approximately) equal to the benchmark combination weight of 1/2, i.e., the equally-weighted combination benchmark.\footnote{ Under log score, a DGP that {ensures that} the optimal two-step forecast combination weight is (approximately) equal to 1/2 can be obtained by setting $ (\phi_1,\phi_2,\sigma^2)=(.4000,-.4421,1)^\top$ in the AR(2) model. For MSFE, taking $(\phi_1,\phi_2,\sigma^2)=(.4000,-.4070,1)^\top$ delivers a DGP such that the optimal (two-step) forecast combination weight is (approximately) 1/2. The values of $(\phi_1,\phi_2,\sigma^2)$ in both cases were obtained numerically by maximizing the corresponding loss function using a sample size of {10} million observations generated from the DGP.} Under this DGP, we generate 1000 replications across three different samples sizes $T=1000,2000$ and $5000$. Across each of these datasets, we use the first $R=T/2$ observations for training, and the remaining $P=T-R$ observations for testing.
The first three rows of Table (ref) compare the size of the standard test of no inferior forecast accuracy, based on $W_P(\alpha)$ in (ref), against the testing approach that caters for the two-step nature of the combination forecast construction, via the use of ${W}^{2s}_{P}(\alpha )$ in (ref). In the two-step approach, we use $B=10000$ draws to simulate the critical value in all cases. As a comparator, we also present the empirical rejection frequency of the standard {`t-test' of} the null hypothesis that $\eta=1/2$, {based on} the correct two-step standard error for $\tilde\eta_R$. This additional {comparator serves as a benchmark of sorts, enabling us to {diagnose}, in some sense, if the rejection rates for the tests are due primarily to differences in the predictive ability of the combinations, or are due to differences in the combination weights themselves.}
The {results} in the table demonstrate that, as already highlighted in Section (ref), the standard approach has zero size in (nearly) all cases when the benchmark combination weight, 1/2 in this case, is close to the optimal combination weight. In contrast, even with an asymptotically conservative critical value, accounting for the two-step nature by which the forecasts were produced, via ${W}^{2s}_{P}(\alpha )$ in (ref), results in a test that has sizes {that} are much closer to the nominal level. Depending on the chosen loss, {for the finite samples used, }the test based on ${W}^{2s}_{P}(\alpha )$ can be slightly {over- or under-sized}, but delivers results that are much closer to the nominal level than {does} the standard approach.
The second set of three rows in Table (ref), compares the {empirical} power of the {alternative} testing approaches, under both loss functions, {under} a DGP that is a relatively small distance away from the DGP that delivers {an optimal combination weight that is} equal to the benchmark combination weight of 1/2.\footnote{ Under log-loss, the DGP generating the observed data is fixed at $ (\phi_1,\phi_2,\sigma^2)=(.40,-.50,1)^\top$, while under MSFE we take the DGP generating the data to be $(\phi_1,\phi_2,\sigma^2)=(.40,-.45,1)^\top$. That is, in both cases, the DGP deviates from the version that delivered equivalence between the equal-weighted combination and the optimally-weighted combination by moving the second root of the AR(2) process, $\phi_2$, by -.05.}Across both loss functions, our testing approach has vastly higher power than the standard approach, {with the standard testing approach based on the log-loss having zero power under this minor deviation from the null hypothesis, and with its power under the MSFE being only around one percent. Moreover, the magnitude of the power of the correctly sized test is broadly similar (for any given sample size) under the two losses. Interestingly enough, the comparator `t-test' has quite high power under log loss, but very different, and much lower power under MSFE.}
Consequently, the results in Table (ref) demonstrate empirically that the if one wishes to conduct a test of forecast accuracy based on approaches that use forecast combinations, accounting for the two-step nature by which these forecasts were produced will be critical in producing tests with good power and reliable size.
As noted, the use of a two-step forecast combination, i.e., $ Q^{(T)}_{\tilde\theta_R}$, is much more common in practice than a one-step forecast combination, i.e., $Q^{(T)}_{\widehat\theta_R}$, where $ \widehat\theta_R$ is defined in (ref). However, the analysis in Section (ref) demonstrates that the forecast combination puzzle is entirely due to this (two-step) estimation approach. This then begs the question of whether it may be possible to avoid the puzzle altogether by changing the way forecast combinations are produced.
In this section we compare the accuracy of the one-step forecast combination, $Q^{(T)}_{\widehat\theta_R}$ against the standard two-step forecast combination, $Q^{(T)}_{\tilde\theta_R}$. Since the two-step approach is the standard approach in the literature, we test that the two-step (benchmark) approach is not inferior to the one-step (alternative) approach. The following result demonstrates that if the one-step forecast combination approach is computationally feasible, it will always yield superior forecast performance.
Theorem (ref) yields the following immediate corollary on the forecast accuracy of a benchmark equally-weighted combination against the optimally estimated one-step combination: let $\theta_R^{ew}=(K^{-1}\iota^{ \prime },\tilde\gamma_R^{\prime })^{\prime }$, with $\iota$ a $K$ -dimensional vector of ones.
In Figure (ref), we repeat the simulation of Section (ref), but for different benchmark-alternative pairs (rows), with a view to illustrating the implications of Theorem (ref) and Corollary (ref), and thereby illustrating the benefits of the one-step combination.\footnote{{We omit results for $\eta^{\star} = 0$ and $ \eta^{\star} = 1$ to exclude cases where $\theta^0 = \theta^{\star}$, whereby the one- and two-step parameter estimators converge to the same values in the limit, violating Assumption (ref).}} As before, discussion of implementation details is left to Appendix (ref).
In the first row, we test the benchmark two-step combination against the alternative one-step combination, and find that the rejection frequency ($y$ -axis) quickly converges to one (in favour of the one-step combination) as the sample size ($x$-axis) increases, for combinations optimising both the MSFE (left-hand column) and the log loss (right-hand column) and for all values of $\eta^{\star}$ (colours). This reflects the result of Theorem (ref), and supports preferring the one-step combination over its two-step counterpart when optimising the combination to maximise forecast performance.
In the second row, the hypothesis of no inferior predictive accuracy of the equally-weighted two-step benchmark is tested against the alternative one-step combination. Here, we again find that the rejection frequency ($y$ -axis) rapidly converges to one as the sample size ($x$-axis) increases, for DGPs of all limiting two-step weights $\eta^{\star}$, including $ \eta^{\star} = 0.5$, where $\theta^{ew}_R$ is the best-performing two-step combination. This is due to the fact that the one-step combination $\hat{ \theta}_T$ converges to a higher-performing combination in the limit than is possible for any two-step combination, including $\theta^{ew}_R$ in the case where the optimal two-step weight $\eta^{\star} = 0.5$; this follows since even if $\eta^{\star} = 0.5$ under the two-step approach, the optimal one-step weight $\eta^0 \neq 0.5$.
Displayed in the bottom row are rejection frequencies for the test of no inferior predictive accuracy of the equally-weighted two-step benchmark against the optimally-weighted two-step alternative. Rejection frequencies for the same test were also displayed in the middle row of Figure (ref), and here we see again that for the two-step combination this test is undersized and has low power, even when the equally-weighted vector $ \theta^{ew}_R$ is far from optimal (dark blue and green), leading to the forecast combination puzzle. Comparing the bottom and middle rows of Figure (ref), we find that the power of the test increases dramatically to resolve the puzzle when estimating parameters in one step (middle row) rather than two (bottom row). The increase in power from one-step estimation is seen across all sample sizes ($x$-axis), all pseudo-true weights (colours) and both losses (columns).
In this section, we give an empirical example which demonstrates that one-step combinations resolve the forecast combination puzzle. {Specifically}, we follow Section 3 of Geweke2011 and consider a linear pool comprising the Gaussian exponential GARCH$(1,1)$ {model} (“EGARCH”) and the GARCH$(1,1)$ model with i.i.d. Student $t$ {errors} (“$t$-GARCH”); this pool is then used to produce one-step-ahead distributional forecasts of daily {logarithmic} S&P500 returns. All parameters (including combination weights) are estimated using returns for the 3783 trading days from years 1990 to 2004, inclusive ({the} “training set”). We then evaluate and compare the log-score-based forecasting performance of different combinations using returns for the out-of-sample period comprising the 3772 trading days from years 2005 to 2019, also inclusive ({the} “test set”).
Three different ways of estimating {the forecast combinations} are compared: equally-weighted two-step estimation, optimally-weighted two-step estimation, and one-step estimation. In {the} two-step combinations, the EGARCH parameters are first chosen to maximize the {average log score of the EGARCH one-step-ahead predictive distribution over the training set}; likewise for the $t$-GARCH parameter estimates. To produce the equally-weighted two-step combination, we set the weights to $0.5$. The combination {weights for the} optimally-weighted two-step combination are estimated by maximizing the log score of the combination across the training set, with the EGARCH and $t$-GARCH parameters fixed at their first-step values. {The one-step combination parameters are jointly} chosen to maximize the training-set average log score of the one-step-ahead predictive distribution of the combination density, in a single optimization program.
Table (ref) contains the training-set average log scores of the three combinations (first column). Consistent with the forecast combination puzzle, even across the training-set the equally-weighted two-step combination outperforms the optimally-weighted two-step combination, since it has a higher average log score (second column). As suspected, both two-step combinations are beaten by the one-step approach over the training-set.
In Table (ref) we display the $p$ values (right column) for three tests of the null hypothesis that a benchmark combination (left column) is not inferior to an alternative combination (middle column). The test proceeds according to Section (ref) using the loss differences pertaining to the log scores of the out-of-sample test set. The asymptotic variance of the loss difference is estimated using the method described in Section 4.1 of Okui2010, with the quadratic spectrum kernel and $S = \sqrt{T}$. In the first row we fail to reject the null hypothesis that the equally-weighted two-step combination is not inferior to the optimally-weighted two step combination, reflecting the forecast combination puzzle. The tests displayed in the second and third rows unequivocally reject the null that the benchmark - either {the} equally- or {the optimally-weighted two-step combination} - is not inferior to the one-step {combination, with both null hypotheses} rejected at the 1% level. This perfectly reflects the theoretical {results} in (ref) and Corollary (ref) in Section (ref): we can avoid the puzzle and obtain a higher performing forecast combination by estimating all parameters in a single step, rather than the standard two steps.
In this paper, we investigate the forecast combination puzzle through the lens of hypothesis testing approaches aimed at discriminating between the relative performances of equally-weighted and optimally-weighted forecast combinations. Forecast combination parameters are optimized according to a scoring function or scoring rule for point forecasting and distributional forecasting, respectively, and we thereby demonstrate that the forecast combination puzzle is a phenomenon that extends far beyond point forecasts optimized according to the MSFE.
Our theoretical analysis demonstrates that such hypothesis tests have no local power, and lack size control, when applied to hypothesis tests aiming to distinguish between the performance of optimally- and equally-weighted two-step forecast combinations -- a result that is unusual for tests of this nature. {This perverse behavior is entirely due to the fact that the test statistic used to conduct the test does not account for the two-step nature by which forecast combinations are produced. Consequently, as we show in the paper, under the null hypothesis the limiting distribution of the (appropriately scaled) test statistic is not asymptotically normal, but instead converges in distribution to a generalized chi-distribution.}
The source of this problematic property is that all uncertainty in the performance of optimally-weighted two-step forecast combinations derives entirely, at first order, from sampling variability in the parameter estimates of the constituent models, with no contribution from the estimation of the weights. It is the relative performance of the different approaches to obtaining the weights, however, that is the subject of the hypothesis test -- a subject that does not contribute (asymptotically, at least) to the sampling variability of the performance measure, resulting in a test with no local power.
An extension to the Monte Carlo exercise of Smith2009 illustrates how this lack of local power can permeate throughout the parameter space. By producing the rejection frequencies of a variety of hypothesis tests of no inferior forecast accuracy of a fixed-weight benchmark against the optimally-weighted two-step alternative, we illustrate that such hypothesis tests can require large sample sizes to reject in favor of the optimally-weighted combination, even when the (unknown) best-performing weights are very different from the vector of equal weights. This finding was seen repeatedly, for point and distributional forecast combinations optimized according to different scores, across a diverse range of DGP parameter values, and across several fixed-weighted (and not just equally-weighted) benchmark combinations. We also revisit a two-model distributional forecast combination of S&P500 returns in Geweke2011 to obtain empirical evidence that this phenomenon occurs in practice.
It is shown that, under mild assumptions, optimizing combination parameters in one step will always eventually reject the null hypothesis that it does not have a higher forecast accuracy than an equally-weighted or optimally-weighted two-step benchmark combination. Repeating the Monte Carlo exercise for the one-step alternative under the optimally- and equally-weighted two-step benchmarks reveals that the low power that is characteristic of the forecast combination puzzle in the case of two-step alternatives is absent when the alternative optimally-weighted combination is estimated in one-step. In addition, we verify the superiority of one-step combinations in the S&P500 returns example considered in Geweke2011 by showing that {the one-step density combination delivers superior predictive accuracy relative to two-step benchmarks.}
In this way, we argue that the root-cause behind the lack of evidence for the performance of optimally-weighted two-step combinations against their equally-weighted counterparts is an artefact of the way in which such combinations are generally produced. Consequently, the “puzzle” is evidenced by the low power observed under standard hypothesis testing approaches used to obtain evidence in favor of, or against, optimally-weighted combinations produced in the standard two-step manner. Furthermore, we demonstrate that if it feasible to produce optimal combinations in a single step, the forecasting puzzle can be completely avoided. Hence, if the problem at hand is such that forecast combinations can be produced in a one-step fashion, the practitioner will (always) reap appreciable gains, in terms of forecast accuracy, by undertaking such a strategy.
Alternatively, if a one-step approach is infeasible, or if a two-step approach is simply preferable, we have demonstrates how the usual testing framework must be altered to accommodate the two-step nature by which the combinations were produced, and to ensure that the resulting behavior does not corrupt our testing results. In particular, we have show that under a broad range of DGPs, when using a two-step forecast combination approach, both the test statistic and the critical value employed must be altered in order to deliver a test that has correct size, and meaningful power.
Before concluding, we note that there are many interesting cases where a forecast combination procedure may seem, at the outset, not to be produced in a two-step fashion, but which upon closer inspection reveals that such forecast combinations are actually produced in (at least) a two-step fashion. As an illustrative example, consider the context of volatility forecasting using the HAR model (corsi2009simple; Corsi2012). By viewing HAR model forecasts as the combination of lagged moving average models for realised variance, Clements2021 document the existence of a forecasting combination puzzle in HAR models and show that such models do not generally perform better than a simple weighted average of the constituent forecasts.
Interestingly, the results of Clements2021 seem to document the existence of a forecast combination puzzle without a `first-stage' estimation step being required to produce the forecasts; at face value it then seems that such an example lies outside the scope of our general results. However, recall that the “observed value” of realised volatility is not a genuine realisation of “observed data”, but a nonparametric estimator of integrated-variance, computed using inter-daily returns. That is, the very construction of the realised variance series constitutes a first-stage estimation step, and forecasts produced via HAR models can therefore be viewed as two-stage forecast combination methods: the first stage estimates the realised variance series, and the second the combination scheme. Hence, the lack of power HAR models exhibit to distinguish between equally and optimally weighted forecast combinations is also explained by our theoretical results. We leave a full study on such types of first-stage estimation steps for future research, but remark that, the heavy use of high-frequency returns, and realised variance in particular, in financial forecasting applications would seem to imply the existence of undiscovered combination puzzles.