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.
90,055 characters · 16 sections · 69 citation commands
Predictive Density Combination Using a Tree-Based Synthesis Function
\thispagestyle{empty}
\doublespacing
It is commonplace when forecasting macroeconomic variables, such as output growth or inflation, to consider a large number of competing predictive densities. These density forecasts might come from different reduced-form or structural models, and/or be subjective and come from surveys. How to combine these densities is an open question being addressed by a growing literature.\footnote{See, among many others, mitchell_evaluating_2005, wallis_combining_2005, hall_combining_2007, geweke_optimal_2011, kk2012, billio2013time, aastveit2014nowcasting, conflitti_optimal_2015, chernis_nowcasting_2022, knotek_real-time_2022, aastveit_quantifying_2022, capek2023macroeconomic, diebold2022aggregation.} The literature concludes that combined density forecasts tend to be more accurate and more robust than single-model approaches that ignore model uncertainty; for a review see Aastveit_Review. One issue is that traditional forecast combination techniques are often linear and do not exploit information besides the forecasts and the target variable. Contrast this with a policymaker who combines forecasts nonlinearly and uses external information, such as on the current state of the economy or financial conditions, to help determine how much weight to attach to the different forecasts. We propose a novel technique that mimics this practice. Our approach combines density forecasts nonparametrically while allowing the combination weights to be determined by information that may be external to the forecasting models.
Key to our approach is Bayesian predictive synthesis (BPS). It has emerged, as extended into a time-series context by mcalinn2019dynamic, as a general method of density forecast combination with a strong theoretical basis. BPS draws on an earlier Bayesian literature on agent or expert opinion analysis west_modelling_1992, and provides a formal and theoretically justified method for pooling densities. It can be shown to nest many previous approaches mcalinn2019dynamic and has been used successfully in various applications in economics, such as mcalinn2020multivariate, chernis2023BPS, and aastveit_quantifying_2022. In this paper we develop density forecast combination strategies within the BPS framework.
In existing implementations of BPS, the so-called “synthesis function,” which determines the weight attached to each density, needs to be specified parametrically. Common choices, as made in the aforementioned papers, are to assume that the synthesis function takes the form of a dynamic linear regression, with parameters that are allowed to change over time typically as random walk processes. This specification of the synthesis function thus allows the weights on competing density forecasts to evolve over time as linear Gaussian random walks. But such an assumption may or may not be valid. Misspecification occurs if the weights depend on other factors or if they follow a different law of motion than a random walk.
These considerations motivate the present paper. BPS has theoretically rigorous foundations, but the manner in which it has been implemented in practice risks misspecification due to the adoption of particular and untested parametric assumptions. We therefore propose to use flexible nonparametric techniques to specify the synthesis function. Specifically, we use regression trees. In conventional (single-model) forecasting applications, tree-based models of the conditional mean have proven highly successful clark2021tail, huber2022inference, huber2020nowcasting. A small number of other papers have used nonparametric techniques to combine predictive densities jin2022infinite, bassetti2018bayesian, bassetti2023inference. However, unlike our proposed method, these other papers neither use regression trees nor fit explicitly within the formal BPS framework.
While regression trees have become a popular way to estimate nonparametric regressions, here we propose to use them differently. Similarly to coulombe2020macroeconomy, deshpande2020vcbart, and hhkm2023tvpbart, who provide a nonparametric treatment to the parameters rather than the variables in single-equation and VAR models, we model the coefficients in the synthesis function with regression trees (RT). Accordingly we label our version of BPS, BPS-RT. The synthesis function remains linear in the parameters, which, as we will demonstrate, aids in interpretation. Use of regression-tree methods requires the choice of covariates, which we call “weight modifiers.” These weight modifiers help determine the weights attached to the competing density forecasts. Conventional BPS does not make use of weight modifiers, given that the weights are typically assumed to follow random walks. Thus, in popular implementations of BPS any relevant information in the form of additional covariates is neglected.\footnote{Notable recent exceptions are Villani, who following LI2023, let the weights in linear density forecast combinations depend on (potentially time-varying) exogenous variables. As Villani explain, such linear pools are one specific instance of BPS. Letting the combination weights in linear pools change over time according to these “pooling variables,” as in the more general (nonlinear) BPS framework that we consider, can offer more flexibility than assuming that the combination weights follow an assumed autoregressive process; cf. DELNEGRO2016.} But decision makers when combining competing density forecasts may wish to condition their forecasts on such “outside” information. For example, they may wish to let the weights on the different forecasting models vary with the state of the economy or vary as a function of the features of each forecast density. Our tree-based specification for the synthesis function is able to condition on both “global” (that is, information not associated with a particular forecaster) and “local” (that is, information associated with a given forecaster) variables when determining the weights. In our tree-based synthesis function, the weights on each density forecast are dynamically determined via a sequence of decision rules. BPS-RT allows the decision maker to combine predictive densities in a highly flexible way and to distill optimally all relevant information contained in the predictive densities and weight modifiers. The fact that the synthesis function remains conditionally linear in the parameters helps the decision maker interpret the combined density and better understand the role each individual density is playing in the combination. We will show how BPS-RT can be used to understand the role of model incompleteness, agent clustering, and the time-varying importance of the different weight modifiers.
The next section of the paper introduces and motivates BPS in theory and then discusses how it has been implemented in the existing literature. It then proposes our generalization, BPS-RT, and explores its properties. Section 3 demonstrates the utility of BPS-RT by undertaking two forecasting applications. The first application takes the individual forecaster density forecasts from the European Central Bank Survey of Professional Forecasters (ECB SPF) and combines them. The second application forecasts US inflation using a commonly used large set of indicators. The predictive densities that are synthesized are produced by regression models using the different indicators. We find that BPS-RT produces well-calibrated and accurate forecasts. Notably, we find that single-tree models perform best, in contrast to standard recommendations when using regression trees. This suggests that a relatively parsimonious weight scheme with few changes in weights is supported by the data. The superior performance of BPS-RT stems from its better ability to explain periods of volatility, such as the global financial crisis that affected euro area GDP growth and the post-COVID inflation period in the US. Zooming in on the best performing BPS-RT specification in the US inflation application, we show how the combination forecasts from BPS-RT can be interpreted. BPS-RT can be used to understand the role of model incompleteness, agent (forecast) clustering, and the time-varying importance of the different weight modifiers. We find little model set incompleteness during the post-COVID inflation period, suggesting that BPS-RT's success comes from its ability to successfully forecast inflation using the underlying models with changes in the combination weights driven by a time trend. This contrasts with the earlier period of lower inflation, when business cycle indicators are shown to be more important. Section 4 concludes. Online Appendix A provides details on Bayesian inference of BPS-RT and Appendix B provides additional empirical results, as referenced in the main paper.
In Section (ref), we provide some background on BPS, distinguishing between BPS in theory and its use in practice in extant empirical applications. Then in Section (ref) we explain how regression trees can be used to provide a more flexible way of operationalizing BPS.
BPS is a foundational theoretically coherent Bayesian method for combining predictive densities.\footnote{For a general description of BPS, see mcalinn2019dynamic; specific implementation details related to our applications are discussed below and in Appendix A.} The theory of BPS provides a pooled predictive distribution for the variable being forecast (say, GDP growth) given a set of individual density forecasts. Operationally, this pooled predictive distribution is produced using Markov chain Monte Carlo (MCMC) methods involving two steps. In the first step, draws are taken from the individual predictive densities for GDP growth. These draws are then, in effect, treated in a second step as explanatory variables in a time-series model where the dependent variable is the outcomes for GDP growth. This time-series model amounts to the synthesis function. Standard choices for this function are typically based on linearity, either simply a constant linear relationship or a dynamic relationship where the linear coefficients evolve over time according to a random walk. As pointed out by aastveit_quantifying_2022, this means that BPS can be thought of as a multivariate regression relating the target variable (GDP growth) to the forecasts for GDP growth, which are treated as generated regressors. We make use of this generated regressor interpretation below.
More formally, at time $t$ a decision maker $\mathcal{D}$ is confronted with $h$-step-ahead forecast densities for variable $y_{\tau+h}$ produced by $J$ different agents, experts, or models, where $\tau$ ranges from 1 to $t$. At each forecast origin, $\tau$, we label these predictive densities $\{\pi_{j \tau}(y_{\tau+h})\}_{j=1}^J$. These densities, available from time $1$ through $t$, form the information set $\mathcal{H}_{t}$ of $\mathcal{D}$ and can, in principle, be of any distributional form. $\mathcal{D}$ then forms an incomplete joint prior $p(y_{t+h}, \mathcal{H}_t) = p(y_{t+h}) \times \mathbb{E}\left(\prod_j \pi_{jt}(x_{jt+h|t}) \right)$ with $\bm x_{t+h|t} = (x_{1t+h|t}, \dots x_{Jt+h|t})'$ denoting latent agent states (that is, draws from the agent-specific predictive densities). These agents' forecasts target $t+h$ but, under the prior, are made using information through time $t$. The prior is incomplete, in the sense that $\mathcal{D}$ only forms an expectation of the product of the agent densities. Agent opinion analysis theory west_modelling_1992-1, west_modelling_1992, extended to a time-series context by mcalinn2019dynamic, shows that the posterior conditional density for $y_{t+h}$ under this incomplete prior takes the form:
where $\alpha(y_{t+h}|\bm x_{t+h|t},\bm \Psi_{t+h})$ denotes the synthesis function that reflects how $\mathcal{D}$ combines her prior information with the set of expert-based forecasts; $\bm \Psi_{t+h}$ denotes a matrix of parameters and latent states that control the properties of the synthesis function, $\alpha(.)$.
Theory offers no guide as to the specific choice of the synthesis function, $\alpha(y_{t+h}|\bm x_{t+h|t},\bm \Psi_{t+h})$. But a common choice in empirical applications, used, for example, in mcalinn2019dynamic, mcalinn2020multivariate, and aastveit_quantifying_2022, is to assume a dynamic linear regression model treating the draws from the $J$ competing densities as generated regressors, $\bm x_{t+h|t}$. Our synthesis functions will have a dynamic regression form, but we will use a non-centered parameterization SFS_W:
where $c_{t+h}$ is a time-varying intercept assumed to follow a random walk, $\bm \gamma = (\gamma_{1}, \dots, \gamma_{J})'$ are time-invariant weights, and $\bm \beta_{t+h} = (\beta_{1t+h}, \dots, \beta_{Jt+h})'$ denotes the time-varying combination weights. As discussed above, a common choice in the literature is to assume that the weights, $ \beta_{jt+h}$, evolve as a random walk (RW) with innovation covariance matrix $\bm V$, leading to a version of BPS that we label “BPS-RW.” When implementing BPS-RW in our empirical applications below, we make standard choices for the prior and MCMC method. In particular, they are similar to those used in hauzenberger2022fast. The only difference is that we use the hyperparameter-free horseshoe prior instead of the normal-Gamma prior, so as to have a prior that is comparable to the one used with our regression-tree model. Accordingly, we do not provide additional details here on drawing the time-varying weights for BPS-RW; see hauzenberger2022fast for details.
In all of our implementations of BPS, including BPS-RW, we consider two versions: one that assumes stochastic volatility (SV) and another that is homoskedastic. In the SV case, the error variance, $\sigma_{t+h}^2$, changes over time. We assume that the log-volatilities $\varsigma_{t+h} := \log \sigma_{t+h}^2$ evolve according to an AR(1) model with autoregressive coefficient $\rho_\varsigma$, unconditional mean $\mu_\varsigma$, initial value $\varsigma_0$, and error variance $\sigma^2_\varsigma$. The prior choices for these parameters are given in Appendix (ref). Homoskedastic cases are obtained setting $\sigma^2_\varsigma$ to zero. Below, for notational ease, we do not explicitly note those parameters relating to SV in the conditioning arguments.
All of our implementations of BPS also include a time-varying intercept, $c_{t+h}$, which is assumed to follow a random walk. As discussed below, $c_{t+h}$ is included to allow for model incompleteness. Further econometric details are provided in Appendix (ref).
With these notational conventions established, $\bm \Psi_{t+h} = \left(\bm \gamma, \{c_\tau, \bm \beta_{\tau}, \sigma_\tau \}^{t+h}_{\tau=0}, \bm \theta \right)$, where $\bm \theta$ will be method-specific parameters that define the law of motion of latent states or appear in the hierarchical priors (such as $\bm V$ in the case of BPS-RW).
The synthesis function, $\alpha(y_{t+h}|\bm x_{t+h|t}, \bm \Psi_{t+h})$, is quite flexible, given that the weights it attaches to each of the $J$ densities are dynamic and because it allows for time-varying error variances. We can also see that while Eq. ({(ref)}) implies a Gaussian density conditional on $\bm \beta_{t+h}$, $\bm x_{t+h|t}$, and $\sigma_{t+h}^2$, when carrying out predictive inference we marginalize out the unknowns of the model, leading to a predictive density that can be highly non-Gaussian; see Eq. ((ref)) below.
In contrast with other approaches to combining models and density forecasts, such as Bayesian model averaging (BMA), the weights on each density are restricted neither to lie between zero and one nor to sum to unity. In the case of BPS-RW, the degree of change in the weights will depend on the magnitude of the state innovation variances for these parameters: small values imply slow, smooth adjustment of the weights over time, large values allow for bigger sharper changes.
Two additional aspects of this parameterization of the synthesis function are worth noting before we introduce our regression-tree approach, which provides a more flexible nonparametric representation of the synthesis function.
First, as a special case, we define a static version of BPS that assumes time-invariant weights $\bm \beta_\tau = \bm 0_J$ for all $\tau$ but leaves $\bm \gamma$ unrestricted. We label this instance of BPS, which assumes the combination weights to be constant over time, “BPS-CONST.”
Second, the presence of both an intercept and an error in the synthesis function means that these versions of BPS allow for model set “incompleteness” Geweke2010. That is, they allow the “true” (but unknown) model not to be in $\mathcal{D}$'s model space; see, for example, billio2013time and AastveitJBES_2018. A conventional model combination scheme such as BMA sets both intercept and error to zero. The fact that the intercept, $c_{t+h}$, and error variance, $\sigma^2_{t+h}$, are both time varying provides additional flexibility when modeling the degree of model set incompleteness. Note that these specific assumptions are equivalent to embedding a popular benchmark for forecasting (especially of inflation) -- the unobserved components SV (UCSV) model of SW_UCSV - within our set of now $J+1$ density forecasts. This is also related to an alternative treatment of model set incompleteness in BPS that adds a fictitious baseline predictive density to the set of densities being synthesized TallmanWest2023.\footnote{diebold2022aggregation also add a fictitious forecaster in their ECB SPF application that, like ours below, combines forecaster-level density forecasts.} In our case, this baseline predictive density comes from a UCSV model. But importantly, as when estimating a mixture density, the parameters of the UCSV density are estimated simultaneously with the weights in the synthesis function.
To carry out predictive inference, we need to compute the predictive distribution. We do so by simulation. Let $y_{T+h}$ denote a future realization of our target variable at time $T+h$ and let $\mathcal{H}_{T}$ denote the set of agent densities that are available at time $T$ but target $T+h$. The predictive density, in our case, is obtained as follows:
where $\mathcal{I}_{T}$ indicates the information set up to time $T$ and $\bm \Psi_{T+h}$ are the latent states (projected forward to time $T+h$). We can simulate from ((ref)) by simulating from the joint posterior of the agents and states, projecting the states forward to time $T+h$, and then using the synthesis function in ((ref)) to produce a combined forecast draw. By doing so, we integrate out the unknowns of the model and the resulting predictive density can be highly non-Gaussian and feature heavy tails, multi-modalities, and/or skewness.
In this paper, our proposal is to relax the restrictions in BPS-RW by considering more flexible forms of time variation in $\bm \beta_{t+h}$. Specifically, we use techniques from machine learning to model the dynamic evolution of the weights, $\bm \beta_{t+h}$, in a nonparametric manner as a function of additional weight modifiers. This treatment can be contrasted with the alternative of treating the function, $\alpha$, itself nonparametrically. We follow chipman2010bart and use Bayesian additive regression trees (BART) to estimate the regression trees. BART consists of a set of priors for the tree structure and the terminal nodes (the leaf parameters) and a likelihood for data in the terminal nodes.
BPS-RT differs from existing implementations of BPS through both the hierarchical priors used on elements in $\bm \gamma$ and $\bm \beta_{t+h}$ and by incorporating additional covariates into $\mathcal{D}$'s information set. These are stored in a $K_\gamma$ vector $\bm z_j^\gamma$ and a $K_\beta$ vector $\bm z^\beta_{j t+h|t}$, both containing additional “data” known to $\mathcal{D}$ through period $t$.
We postulate a nonlinear relationship between the weights and the weight modifiers through functions $\mu^{\gamma}_j(\bm z_j^\gamma)$ and $\mu^{\beta}_{j}(\bm z_{jt+h|t}^\beta)$ that determine the state transition equation that can be interpreted as a prior. In particular, we assume:
where $\tau^{\gamma}_j$ and $\tau^\beta_j$ denote prior scaling parameters. For convenience, we define $\mu^{\gamma}_j := \mu^{\gamma}_j(\bm z_j^\gamma)$ and $\mu^\beta_{jt+h} := \mu^{\beta}_{j}(\bm z_{jt+h|t}^\beta)$. The best way to illustrate the effect the scaling parameters have on the actual estimates of the weights is to consider a re-parameterization of the synthesis function. Integrating out $\gamma_j$ and $\beta_{jt+h}$ by plugging Eq. ((ref)) into Eq. ((ref)) yields:
with $\nu^\gamma_j, \nu^\beta_{jt+h} \sim \mathcal{N}(0, 1)$ denoting process innovations. The innovations, $\nu^\gamma_j$ and $\nu^\beta_{jt+h}$, and the corresponding scaling terms control the degree of dispersion of the actual weights from those expected under the prior mean. If the scalings are close to zero, the posterior of $\gamma_j$ and $\beta_{jt+h}$ is pulled toward the prior mean and the resulting estimates will be close to $\mu^{\gamma}_j(\bm z_j^\gamma)$ and $\mu^{\beta}_{j}(\bm z_{jt+h|t}^\beta)$ and so strongly depend on $\bm z_j^\gamma$ and $\bm z_{jt+h|t}^\beta$. If this is not the case, the resulting scaling parameters would be larger so that substantial deviations from the prior means are more likely. Another feature of this representation is worth emphasizing. As opposed to a model that directly approximates the synthesis function nonparametrically, the specification in ((ref)) introduces interaction terms of the form $\mu_j^\gamma(\bm z_j^\gamma) \times x_{jt+h|t}$. This specific form might reduce the risk of overfitting by introducing more structure on the space of functions that we approximate.
We approximate the prior mean functions through a sum-of-trees model with $S$ trees:
where $g$ denotes a tree function that is parameterized by so-called tree structures, $\mathcal{T}^{n}_s$, and terminal node parameters, $\bm \phi^n_s$, for $n \in \{\beta, \gamma\}$. The basic idea behind a single tree is that the tree structures describe sequences of disjoint sets. These sets partition the input space (determined by exogenous covariates, $\bm z^\gamma_{j}$ and $\bm z^\beta_{jt+h|t}$, respectively). Each of these sets is associated with a particular terminal node parameter. In our case, the terminal node parameters serve as prior expectations for the $\gamma_j$s and for the $\beta_{jt+h}$s. The input space is associated with vectors of variables, $\bm z^{\gamma}_{j}$ and $\bm z^{\beta}_{jt+h|t}$, which we refer to as weight modifiers. Note that $\bm z^{\beta}_{jt+h|t}$ could include quantities (such as moments from the agent-specific predictive densities) that explicitly target $t+h$ but are available at time $t$.
There are two main justifications for our BPS-RT modeling approach. First, there is no reason to restrict attention, as in BPS-RW, to random walk specifications for the evolution of $\beta_{t+h}$. BPS-RW implies, at a given point in time, a linear relationship between $y_{t+h}$ and $\bm x_{t+h|t}$. This assumption might be warranted in tranquil periods. But, in unusual times, nonlinearities could be present, and exploiting these might lead to more accurate forecasts. Our regression-tree approach allows for flexibility in the way such nonlinearities are modeled and lets the “data speak.” Second, and this holds across all existing instances of BPS not just BPS-RW, an implicit assumption made is that the information set available to $\mathcal{D}$ comprises exclusively the agent-based forecast densities.\footnote{As mentioned in footnote 2, an exception is Villani, who when combining density forecasts using the linear opinion pool also let the weights depend on exogenous variables. Our BPS-RT model generalizes to consider BPS combinations beyond the linear special case and to allow for nonlinearities in how the weight modifiers affect the weights.} But, in principle, additional unmodeled information is available to $\mathcal{D}$ and might help inform evolution of the weights. In our BPS-RT approach, the weight modifiers, $\bm z^{\gamma}_{j}$ and $\bm z^{\beta}_{jt+h|t}$, comprise this extra information.
These weight modifiers might include characteristics of the agents' forecasts not directly reflected in their predictive distributions or other common (to agents) factors, such as general information about the macroeconomic environment. For example, $\bm z^{\gamma}_{j}$ might contain summary metrics of overall past forecast performance (such as the average historical forecast performance) for each agent. Or, as noted above, $\bm z^{\beta}_{jt+h|t}$ might contain more granular and time-varying information, such as time-varying characteristics of the agent-specific predictive densities (say their higher moments and/or time-varying measures of past forecasting performance). We provide specific context and motivate our choice of weight modifiers in the empirical applications in Section (ref) below.
To return to the regression tree, note that it is defined by disjoint sets that are determined by splitting rules of the form $z^{\beta}_{k,jt+h|t} \le d_k$ or $z^{\beta}_{k,jt+h|t} > d_k$, where $z^{\beta}_{k,jt+h|t}$ is the $k\textsuperscript{th}$ weight modifier for the $j\textsuperscript{th}$ agent/model and $d_{k}$ is a threshold parameter associated with the $k$\textsuperscript{th} effect modifier, which is estimated from the data. It is important to note, however, that any splitting rule associated with the $k\textsuperscript{th}$ effect modifier is common across agents and periods (that is, it is specific neither to agent $j$ nor to period $t$). Hence, the thresholds $d_k$ and thus the tree structures do not have $t$ or $j$ subscripts and are common across agents/models and time. Since these splitting rules effectively govern the prior mean, $\mu^{\beta}_{jt+h}$, this structure in a sense captures the notion of a pooling prior and reflects the situation that $\mathcal{D}$ decides on the weights associated to each of the different agents based on using additional factors $\bm z^{\gamma}_{j}$ and $\bm z^{\beta}_{jt+h|t}$ according to a set of common decision/splitting rules. The same structure also holds for the $\gamma_j$s, with the difference that the splitting rules controlling $\mu^{\gamma}_{j}$ pool exclusively over the cross-section and not over time (since the $\gamma_j$s are time invariant).
To see this pooling feature more clearly, consider a BPS-RT model that assumes $\bm \beta_{t+h} = \bm 0_J$ and features only a time-variant part $\bm \gamma$, for which the prior mean $\bm \mu^{\gamma} = (\mu^{\gamma}_1, \dots, \mu^{\gamma}_J)'$ is defined by a single tree ($S=1$) and by a single effect modifier in $z_j^\gamma$ (that is, $z_j^\gamma$ is a scalar with $K_\gamma = 1$). In this case, the prior on $\gamma_j$ can be written as:
If we now compute the difference between $\gamma_j$ and $\gamma_m$ for distinct agents, $j \neq m$, and assume that $z_j^\gamma$ and $z_m^\gamma$ are similar, in the sense that both imply the same decomposition of the input space and are thus located in the same terminal node of the tree, we end up with:
This equation implies that if the tree suggests that the characteristics between agents are so similar that they are grouped together in the same terminal node, the same prior mean applies and the difference between prior means will be zero. The presence of the prior scaling parameters $\tau_j^\gamma$ and $\tau_m^\gamma$ will then allow for data-based testing of whether that restriction should be strongly enforced or not. Since both prior means would coincide, setting both $\tau_j^\gamma$ and $\tau_m^\gamma$ to values close to zero would induce a clustering of $\gamma_j$ and $\gamma_m$ around $g(z_j^\gamma|\mathcal{T}_s^\gamma, \bm \phi_s^\gamma) = g(z_m^\gamma|\mathcal{T}_s^\gamma, \bm \phi_s^\gamma)$. Hence, the choice of the prior specified on the scaling parameter $\tau^\gamma_j$ is crucial in determining the clustering behavior of BPS-RT.
Another feature of our prior is that $\mathcal{D}$ adjusts her weights on the agents' densities depending on the (common) macroeconomic environment as captured by the weight modifiers, which might include, as discussed, indicators of the state of the business cycle, measures of economic uncertainty, or deterministic trends. For example, in turbulent times larger weights on component densities that are far from Gaussian and feature, say, heavy tails might lead to better combined density forecasts. Our approach can control for this, if supported by the data.
Note that we estimate the tree structures and the terminal parameters alongside all other unknown parameters and therefore also specify priors for them. We follow here the recommendations of chipman2010bart and discuss the remaining model and prior specification issues in detail in Appendix (ref). This technical appendix also describes the MCMC methods used to estimate BPS-RT. In summary, these MCMC methods are straightforward. They require a method for predictive simulation from each individual model (to draw from each agent's forecast density) and a method for drawing from the regression-tree model conditional on the individual-agent draws. For BPS-RT the algorithm is taken directly from chipman2010bart.
We now explain how BPS-RT works and allocates combination weights using an illustrative toy example. Assume that, unknown to $\mathcal{D}$, the “true” data for $y_t$ are generated by the following threshold model:
where $\rho_1 = 0.8, \rho_2 = -0.8$ and $\sigma_0 = 1.2$, $y_0 = y_1 = 0$, $c=1/100$, and $\nu_t \sim \mathcal{N}(0, 1)$.
Then, $J=2$ agents predict $y_t$ as follows (these forecasts are one-step-ahead, $h = 1$):
Both agents use forecasting methods with a different type of misspecification. The first agent's forecast is almost correctly specified for the first part of the sample, but the second agent's is substantially misspecified. In the second part of the sample this switches. We would hope that BPS-RT, when combining these two misspecified densities, would put more weight on the first agent when $t \le 200$, then increase the weight on agent 2 when $t > 200$.
Notice that the structure of the data-generating process (DGP) implies that BPS-RW is severely misspecified, since BPS-RW implies that the combination weights on the two agents evolve smoothly over time. Our more flexible choice of synthesis function, ((ref)), conditional on choosing appropriate effect modifiers, as we shall show, is capable of accommodating the break at $t=200$.
We consider three variables as weight modifiers. The first is a simple deterministic time trend, $z^{\beta}_{1,jt+1|t} = t+1$, that is common to both agents. The remaining two effect modifiers are agent-specific and measure historical forecasting performance. To capture historical point forecasting performance, we consider each agent's squared forecast error (SFE) as recursively computed at time $t-1$: $z^{\beta}_{2,jt+1|t} =(y_{t} - \mathbb{E}(x_{jt|t-1}))^2$ for $j=1,2$. Then to measure past density forecasting performance, we consider each agent's continuous ranked probability score (CRPS).\footnote{If $F$ is the c.d.f. of the forecast and $y$ the subsequent realization, then $\text{CRPS}(F,y) = \int (F(x) -\mathbf{1}_{x \geq y} )^2 dx$. }
Our synthesis function is given by Eq. ((ref)). To facilitate illustration of BPS-RT, we make some simplifying assumptions. We set the time-invariant weights $\bm \gamma=\bm 0$ and, for the prior on $\bm \beta_{t+1}$, set the scaling parameters equal to zero so that the weights and prior means coincide, and we focus on the single-tree case ($S=1$). For expositional ease, we drop corresponding sub- and super-scripts when there is no loss in meaning. Under these simplifying assumptions, the synthesis function, similarly to ((ref)), reduces to:
This equation shows that with the scaling parameters set equal to zero, we end up with a BART model that assumes the weights depend nonlinearily on $\bm z_{t+1|t}$.
Figure (ref) depicts in panel (a) the estimated tree and in panel (b) the temporal evolution of the estimated weights. We emphasize that these weights are in-sample estimates, that is, conditional on data through $T=350$.
The tree in panel (a) can be understood as follows. Let us start at the bottom of the tree. We see five terminal nodes. Hence, we observe five groups/clusters that define the prior mean both over time and across agents. Put differently, there are five “breaks” over time and across agents in the prior mean.
How we pool is defined by the splitting rules. These are understood by turning to the top of the tree. At the root (level 0), the SFE is used as a splitting variable. The threshold parameter is $1.8$ and, hence, if the SFE in $t-1$ is larger than or equal to $1.8$, we move down the left branch of the tree. At the first level, the lagged CRPS shows up as the next threshold variable. If the CRPS is smaller than $1.3$, we end up in a terminal node and set the weight associated with an agent that has an SFE greater than or equal to $1.8$ and a CRPS smaller than $1.3$ equal to $\mathbb{E}(\beta_{jt})=0.054$. These conditions are fulfilled $21$ percent of the time. By contrast, if the CRPS is greater than or equal to $1.3$, we drop down to the second level of the tree. In this segment, time shows up as a splitting variable and if $t \ge 201$, we assign a weight of $0.72$. For $t < 201$ we introduce a further splitting rule that splits the sample once more by testing whether $t < 42$. If this is the case a negative weight of $-0.062$ is applied, whereas if $42 \le t < 201$ the weight is $0.15$. If the past SFE is smaller than $1.8$, we end up in the right branch of the tree and assign a weight equal to $0.7$.
Hence, the tree suggests that, first and foremost, $\mathcal{D}$ selects agents according to the past performance of their forecasts, since both SFE and CRPS are identified in the estimated tree. Under our DGP, this implies that weights dynamically update if a given agent issued a poor prediction in the previous period without taking into account the past performance of her forecasts. To understand how these decision rules translate into the actual evolution of model weights, panel (b) shows the weights over time. These indicate that in the first part of the sample, Agent 1 receives substantial weight, while Agent 2 receives relatively little weight. This makes sense, given that the former is only mildly misspecified, whereas the latter features substantial model misspecification. As expected, given the structural break in the DGP, $\mathcal{D}$ now overweights the second agent, whereas the weight on Agent 1 is now much smaller.
This simple exercise illustrates how $\mathcal{D}$ incorporates additional information (time and past forecast errors in this case) to combine models. In general, though, the prior scaling parameters in BPS-RT are greater than zero, and hence, the decision tree gives rise to prior expectations that, in turn, inform the posterior estimates of the weights. Hence, if there is no relationship between the weights and the weight modifiers, the resulting prior variance would be large and the weights would follow a white noise process.
We investigate the performance of BPS-RT in two forecasting exercises. In the first application, we combine predictive densities of GDP growth for the euro area (EA) produced by individual professional forecasters participating in the ECB Survey of Professional Forecasters (SPF). Beyond its intrinsic interest, this data set is a good testing ground for BPS-RT because it has been used before when comparing alternative density forecast combination methods; see diebold2022aggregation, conflitti_optimal_2015, and chernis2023BPS. Second, we forecast US inflation using a set of autoregressive distributed lag (ADL) regression models. This data set and model set has been used by stock2003forecasting and rossi2014evaluating, the latter using a similar ADL strategy to create each of the agent's forecast densities.
These two applications differ not only geographically and in terms of target variables, but also in the number of agents and the nature of the forecast densities the agents provide. The EA GDP growth application features a relatively small number of subjective, most likely judgment-informed, forecasts european_central_bank_results_2019 that are provided in the form of histograms (with $J = 14$). In contrast, the US inflation application uses a large number of model-based predictive densities, which are continuous and produced with distinct ADL regressions (with $J = 56$). Further details on the design of both applications are provided in the subsequent sub-sections (ref) and (ref). Both applications' evaluation samples cover the global financial crisis, the euro area crisis, and the COVID-19 pandemic. Taken together, these two applications enable a comprehensive assessment of BPS-RT.
We experiment with several different specifications of BPS-RT to draw out how density forecast accuracy varies with the characteristics of the specific synthesis function used. In broad strokes, we look at the importance of time variation, in both weights and volatility, the number of trees, and the choice of weight modifiers. Accommodating temporal instabilities RossiJEL is important in macro-modeling, and so is a natural subject of inquiry, while the number of trees is an important aspect of specifying BART models. Being able to specify weight modifiers is an attractive feature of BPS-RT and allows the combination weights to change based on information exogenous to the individual agents but known to the BPS decision maker. Hence this is also a key area of inquiry.
We accordingly investigate the following four specifications of BPS-RT distinguished by their choice of weight modifier(s) and whether that choice introduces cross-sectional (which we label C) or cross-sectional and time variation (which we label TC) in the combination weights seen in ((ref)).
For each of these four versions of the model, we consider models with SV and homoskedastic errors and we allow the BART specification to either have a single tree ($S = 1$), leading to a Bayesian regression-tree specification chipman1998bayesian, or a large number of trees ($S = 250$), leading to BART. In traditional Bayesian implementations using trees for nonlinear regression, such as chipman2010bart, it is generally found that increasing the number of trees, starting at $S=1$, leads to an improvement in forecast performance. But this improvement tends to peter out when the number of trees gets moderately large. The conventional wisdom is that the precise choice of the number of trees is not that important, provided that too small a value is not chosen. This may not be the case in BPS, since the data may prefer to have weights that are reasonably constant over time and change only occasionally. Hence, we choose to focus on single-tree specifications and BART to model the weights in BPS. As we shall see, we find that single-tree methods tend to forecast more accurately. As benchmarks in the forecasting exercises below, we consider both BPS-CONST and BPS-RW (as defined in Section 2.1.2).
\color{black} The ECB has been producing the SPF since $1999$. The ECB SPF is the longest running EA survey of macroeconomic forecasts. Each quarter, the survey elicits from a panel of professional forecasters point and probability forecasts of EA inflation and GDP growth at various horizons.\footnote{For a full description of the EA SPF, see garcia_introduction_2003.} We consider the two-quarter-ahead forecasts of year-on-year EA GDP growth. On average, there are $50$ responses a quarter from a survey panel of over $100$ professional forecasters.
There are a couple of features of the forecaster-level density forecasts from the ECB SPF that we have to address in order to combine them. First, survey respondents provide their probability forecasts over given (fixed) ranges. That is, they produce histogram rather than continuous density forecasts. For example, in the $1999$Q$1$ survey, forecasters were instructed to provide their probability forecasts over $10$ bins. The first bin was GDP growth less than 0 percent, with the bins then increasing in intervals of $50$ basis points, until the tenth bin of higher than 4 percent growth. To accommodate the discretized nature of these probability forecasts, rather than fit a continuous density to the histogram (that may or may not have a good fit), we use the histogram forecast data “as is." We do this by, within our BPS approach, drawing samples for each forecaster directly from the histograms. Details of our algorithm, which involves a Metropolis-Hastings step, are given in Appendix A.2. Our sampling approach changes over time to capture the fact that the bin definitions have been moved over time. In particular, after shocks such as the global financial crisis and COVID-19, the ECB shifted the bins to allow forecasters to say more about the probabilities in what were, prior to the survey change, the extremes of the distribution. We also have to take a stand on the open intervals at the bottom and top of the histogram. We set the end-points for the histograms equal to the outer bin plus or minus (depending on whether we are at the top or bottom of the histogram) two standard deviations of GDP growth, as estimated using the vintage of GDP data available at the time the forecast was made.
Second, forecasters enter and exit the panel. This means that the panel is unbalanced. We follow diebold2022aggregation in constructing the longest consistent panel possible by dropping forecasters who are regular non-responders and then filling in the occasional missing values for the remaining forecasters. Specifically, we drop forecasters who have not responded for five or more consecutive quarters. This results in a panel of $14$ forecasters. Any missing forecast data for these $14$ forecasters are estimated using a Normal distribution based on the unconditional distribution of GDP growth as estimated in real time.\footnote{We differ from diebold2022aggregation in two ways. First, they interpolate missing forecasts based on historical performance. Second, we have a different number of forecasts because we use a different sample and we forecast GDP growth instead of inflation.}
We then take these 14 forecasters' densities and carry out a recursive out-of-sample evaluation of the alternative BPS specifications over the sample $2005$Q$2$ through $2021$Q$1$. To do this, we first estimate the BPS combinations on a set of training samples that comprise a sequence of expanding windows of GDP and density forecast data. The GDP data used in the training sample are that vintage of GDP data available to the forecasters when they made their forecasts. The first training sample uses forecasts from the five-year period targeting GDP outturns from $1999$Q$3$ through $2004$Q$2$. These forecasts are taken from the surveys administered between $1999$Q$1$ and $2003$Q$4$. Given its publication lags and our desire to approximate the information set available at the time the SPF forecasts are publicly available, the GDP outturns required to estimate the BPS synthesis function over this training sample are taken from the $2004$Q$4$ vintage. This estimated synthesis function then uses the $2004$Q$4$ survey to forecast (out-of-sample) $2005$Q$2$. The training sample and vintage of GDP data are then extended by one quarter, and forecasts are produced for $2005$Q$3$. This process is continued until forecasts are produced for $2021$Q$1$. This set of out-of-sample BPS density forecasts is then evaluated against GDP outturns taken from the June $9$, $2021$ vintage.
We follow rossi2014evaluating and construct density forecasts of US inflation using a set of autoregressive distributed lag (ADL) models. Each ADL model considers 1 of $27$ indicators taken from the FRED-QD data set McCrackenNg, which is commonly used when forecasting macroeconomic aggregates such as inflation in the US. The selected indicators capture movements in assets prices, measures of real economic activity, wages and prices, and money. This rich and diverse set of economic indicators allows the ADL density forecasts of US inflation to display significant heterogeneity. Table (ref) in the Appendix provides an overview of the variables used as exogenous predictors and the transformations applied to ensure their stationarity.
We then use each of these ADL models to produce direct forecasts for quarter-on-quarter consumer price (CPIAUCSL) inflation one-quarter-ahead ($h = 1$) and one-year-ahead ($h = 4$). Specifically, for each indicator, $x_{jt}$, for $j=1,...,27$, we estimate the set of ADL models:
where $\pi_t$ is inflation, $\rho_{\pi}$ is the autoregressive coefficient, and $\alpha_{\pi}$ denotes the coefficient related to the $j\textsuperscript{th}$ exogenous indicator.\footnote{For notational ease, we do not use $j$ subscripts to distinguish parameters in Eq. ((ref)).} We supplement these $j=1,\dots,27$ models with a $28\textsuperscript{th}$ model (the AR(1) model) that sets $x_t = 0$ in Eq. ((ref)). We also allow $\sigma_{\pi, t+h}^2$, the error variance, to be both time-varying and constant. Hence, we estimate $28$ models both with and without SV, delivering, in total, a set of $56$ individual models whose density forecasts we then combine using BPS. All 56 models are estimated using standard Bayesian techniques. Details are provided in Appendix (ref).
We first estimate these models on a training sample from $1970$Q$1$ to $1989$Q$4$. We then iterate forward using a rolling estimation window of $80$ quarters to account for possible structural changes in the US economy. The first ten years of forecasts ($1990$Q$1$ to $1999$Q$4$) are used as a training window to estimate the BPS synthesis functions. The combined forecasts are then assessed on the evaluation sample $2000$Q$1$ to $2022$Q$4$. This evaluation period includes distinct economic periods characterized by different inflation dynamics, including the dotcom crash, the global financial crisis, the COVID-19 period, and the post-pandemic inflationary period.
We break the empirical results into three parts presented in the following three sub-sections. First, we evaluate the relative and absolute density forecast accuracy of BPS-RT. Second, we examine why BPS-RT forecasts more accurately than the benchmarks by comparing features of their forecast densities. Third, we demonstrate aspects of interpretability of BPS-RT by examining how BPS-RT can be used to understand the role of model incompleteness, agent clustering, and the time-varying importance of the different effect modifiers.
We evaluate forecast accuracy in several ways. We first evaluate the point (conditional mean) forecasts, extracted from the combined densities, using the root mean squared forecast error (RMSE) loss function. Second, we evaluate the full predictive densities. We emphasize evaluation of the predictive densities rather than the point forecasts. Since the loss functions of forecast users tend not to be quadratic --- as the density forecast literature Aastveit_Review emphasizes --- it is always important to produce and evaluate complete probabilistic forecasts. We measure the relative forecast accuracy of the forecast densities using two popular metrics: CRPS and a tail-weighted CRPS. Both are loss functions that score the density forecast according to the realization that subsequently materializes. CRPS evaluates the “whole" density, while tail-weighted CRPS focuses on accuracy in the tails GneitingRanjan.\footnote{In the empirical appendix we follow GneitingRanjan and break CRPS tails into their left and right tails. See Figures (ref) and (ref) in Appendix (ref).} We also test the absolute calibration of the combined density forecasts using the rossi2019alternative test on the probability integral transforms (PITs); and we assess the temporal stability of forecast performance using the fluctuation test of giacomini2010forecast. The results of both these tests are summarized below, with full results presented in Appendix (ref).
Figures (ref) and (ref) report the relative forecast performance of the different models in the EA GDP growth and US inflation applications, respectively, using the RMSE, CRPS, and CRPS-tails loss functions. Each row in these figures reports the relative (to the BPS-RW benchmark) performance of the four BPS-RT specifications as differentiated by whether they use a single tree or 250 trees and whether they have SV or homoskedastic errors. The four columns in the figures refer to which set of weight modifiers is used.
Looking first at the RMSE panel in Figure (ref) for EA GDP growth, we see little difference between the alternative BPS-RT specifications in terms of their point forecast accuracy. The accuracy of the BPS-RT specifications also tends to be similar to that of BPS-CONST and BPS-RW, with gains/losses in general only around 3 percent. This supports the stylized fact from the forecasting literature that equal-weighted combinations of point forecasts are hard to beat TIMMERMANN2006. Turning to US inflation (Figure (ref)), we do see in the RMSE panel that some of the tree-based methods now improve upon the point forecast accuracy of both benchmarks and in a manner that is statistically significant. Of particular note is the superior performance of the single-tree models, which almost always outperform the more complicated $250$-tree models. We discuss this finding further below.
The CRPS panels in both Figures (ref) and (ref) reveal yet more of a payoff to using BPS-RT, certainly relative to BPS-RW, when we evaluate the whole density. Many of the forecast accuracy gains for BPS-RT are statistically significant. An implication of this finding is that BPS-RW's assumption that the combination weights follow a random walk is not supported by the data. But BPS-CONST, especially when BPS allows for SV, remains competitive for EA GDP growth.
The CRPS and CRPS tail results echo those under RMSE loss in concluding that single-tree structures, $S=1$, are almost always preferred to $S=250$. The fact that a single-tree model produces more accurate forecasts contrasts with the conventional wisdom in the wider BART literature; see chipman2010bart. In our case, however, we model the weights, rather than the observed outcomes, nonparametrically and hence the implied conditional mean relation (see (ref)) introduces more restrictions relative to a standard BART model and hence lessens the risk of overfitting.
While the benefits of allowing for SV are well established in the density forecast literature clark2011, allowing for SV in the BPS combination does not obviously improve the density forecasts from BPS-RT. But recall, and we touch on this again below when showing that these models in fact receive higher combination weights, in the US inflation application half of the components models themselves allow for SV.
We now focus on comparing forecast accuracy across the first four columns of both Figures (ref) and (ref). This comparison reveals that the choice of weight modifier does affect forecast accuracy. It is not always the case that using more weight modifiers delivers more accurate forecasts. The benefit of different modifiers varies by application and by which row (which of the four BPS-RT specifications) is consulted.
Finally, we summarize the results from both the PITs calibration tests and the fluctuation tests. These results are reported in the online appendix for space reasons. The PITs plots (see Figure (ref)) show that the BPS-RT densities are well calibrated and especially so when forecasting EA GDP growth or US inflation one-quarter-ahead. The fluctuation test of giacomini2010forecast reveals that there is temporal variation in the relative performance (under CRPS loss) of the preferred BPS-RT model and BPW-RW. Results (see Figure (ref)) indicate that the superior performance of BPS-RT in the EA application is due to better forecasting performance toward the end of the global financial crisis. For US inflation, the better accuracy of BPS-RT is explained by its more accurate density forecasts in the post-lockdown inflationary period.
In this section we examine how and why BPS-RT forecasts more accurately. We focus on the best performing (most accurate) model in each application and compare its forecast densities to those of the benchmark model, BPS-RW.\footnote{As seen from Figures (ref) and (ref), in the EA GDP growth application, the “best" BPS-RT specification has a single tree and SV and uses average scores as effect modifiers (i.e., RT(C): AVG.-SCORES). For the US inflation application, the “best” BPS-RT specification has a single tree, homoskedastic errors, and the full set of weight modifiers (i.e., RT(TC): ALL).}
Figure (ref) shows a heat map of the difference in probabilities, in intervals of $1.5$ percentage point for EA GDP growth and of $1$ percentage point for US inflation, between BPS-RT and BPS-RW. Green (red) shading indicates that BPS-RT adds (subtracts) probability relative to BPS-RW in that interval. This is the approach pioneered by diebold2022aggregation as a way of visualizing the differences between competing density forecasts.\footnote{For an alternative but complementary visualization, Figure (ref) in Appendix (ref) shows the temporal evolution of the underlying density forecasts from BPS-RT and the benchmark BPS-RW model over the EA and US evaluation samples.}
Panel (a) of Figure (ref) shows that, in general, BPS-RT predictions are less dispersed than BPS-RW with more mass near the subsequent outcomes. Additionally, the BPS-RT density adds probability to low GDP growth outturns prior to the financial crisis and also forecasts higher growth than BPS-RW in both the post-global financial crisis recovery and the rebound from the COVID-19-induced recession.\footnote{As shown in Figure (ref) in the online appendix, in moving the probability mass from the centers to the left tail of the forecast density, BPS-RT captures asymmetries in the forecast densities. While there is some evidence of heightened downside risk asymmetries to GDP growth in the course of the financial crisis, consistent with the growth-at-risk literature Adrian2019, the evidence for negative skew is stronger still during the COVID-19 pandemic.} Panels (b) and (c) of Figure (ref) show the analogous plots for US inflation. Similar to panel (a), BPS-RT places more mass closer to the outturn and produces forecasts that are, in general, less disperse. Moreover, BPS-RT adjusts much more quickly to the increase in inflation post-pandemic, both one-quarter- and one-year-ahead, attributing a higher probability to these outturns than BPS-RW. Consistent with the evidence in rossi2014evaluating that combinations of predictive densities for US inflation appear to be approximately Gaussian, the inflation forecast densities from BPS-RT also tend to be symmetric (see Figure (ref) in the online appendix), although there is clear evidence of a spike in downside risks in 2011, a time when the Fed was engaged in quantitative easing to combat deflation threats.
This section discusses how $\mathcal{D}$ can interpret the combined forecasts from BPS-RT. In so doing we continue to focus on the preferred BPS-RT specification in the US inflation application, not least because this is where we observe greater differences across the competing combination strategies. We first show how to quantify the degree of model set incompleteness, as a way of assessing how well the agents (the $J$ forecasting models) that BPS-RT is combining are actually able to forecast. Second, we assess the relative importance of individual weight modifiers in driving BPS-RT.
To measure model set incompleteness we compute an $R^2$-type measure. This estimates the proportion of the variation in $y_{t+h}$ that is explained by the $J$ agents. This measure is computed, for a specific period in the evaluation sample, as the ratio between the variation in the conditional mean in Eq. ((ref)) explained exclusively by the BPS-RT component --- which is the conditional mean in Eq. ((ref)) without the time-varying intercept $c_{t+h}$ --- and the overall variation of the target variable, $y_{t+h}$. $R^2$ values close to zero signify a high degree of model incompleteness, which means that the agents' forecasts are not informative about the target variable. Instead, the intercept and error term in the BPS synthesis function, Eq. ((ref)), explain a large portion of the total variation. In contrast, $R^2$ values close to one indicate that the agents' forecasts are informative and account for the majority of the variation, implying a complete model space.
Figure (ref) plots this $R^2$-type estimate over the evaluation sample. Given that it is computed recursively, quarter-by-quarter, it experiences some volatility. But Figure (ref) still evidences meaningful temporal variations in the degree of model set incompleteness at both forecast horizons. We see higher model incompleteness for the one-year-ahead forecasts than for the one-quarter-ahead forecasts. This is not surprising, as producing longer-horizon forecasts is obviously more difficult. At both horizons, we see increases in model incompleteness during the period 2004-2008, a time of extreme oil price volatility as well as the global financial crisis and in the disinflation period after the 2015 oil price shock.
Interestingly, there is no clear evidence of an increase in model incompleteness during the post-pandemic rise in inflation, reinforcing the message from Figure (ref) that BPS-RT was better able to anticipate the 2021 rise in US inflation.
We now turn to assessing the relative importance of the individual weight modifiers in driving the density forecasts from BPS-RT. We do so by looking first at the number of tree splits and then by calculating inclusion probabilities for each weight modifier. Inclusion probabilities are calculated as the number of splits associated with the respective weight modifier divided by the total number of splits. For space reasons, we focus our discussion on Figure (ref), which examines the weight modifiers for forecasting US inflation one-quarter-ahead. Analogous results forecasting inflation one-year-ahead are reported in Figure (ref) and summarized below when the conclusions differ markedly from those discussed in greater detail for the one-quarter-ahead forecasts.
We start in panel (a) of Figure (ref) by plotting the evolution of the total number of tree splits over the evaluation sample. This panel indicates whether variability in the combination weights comes from the time-varying ($\beta_{jt+h}$) or constant component ($\gamma_j$) of BPS-RT. Panel (a) reveals that BPS-RT tends to select a relatively small number of tree splits, especially for the time-invariant weights. Typically for $\gamma_j$ we observe that the posterior mean of the number of tree splits lies between $0.52$ (lower quartile over the evaluation sample) and $1.15$ (upper quartile, with a few more exceptions in the upper tail), while the average over the evaluation sample is $1.28$. On the other hand, the posterior mean number of tree splits for the time-varying weights, $\beta_{jt+h}$, ranges from $1.08$ to $1.68$ (indicating the interquartile range) and has an average of $1.59$ over the evaluation sample. To place these numbers in the context of a single-tree split on, for example, $\gamma_j$ indicates that the combination weights tend to cluster around two distinct prior means. With this in mind, we interpret the results in panel (a) as showing that the combination weights often fall into a handful of clusters that are more likely to be determined by time-specific factors. However, the number of splits is modest, so the weights are relatively stable over time. This finding is consistent with the density forecast combination literature that finds that constant weight combinations can forecast well chernis2023BPS.
Panels (b) and (c) of Figure (ref) then show the inclusion probabilities for each of the constant and time-varying weight modifiers. Panel (b) shows the inclusion probabilities for the weight modifiers (CRPS and MSE) used to model the time-invariant combination weights. Neither CRPS nor MSE is obviously more important. Both weight modifiers receive positive and often fairly similar probabilities of inclusion. This implies that BPS-RT does partition models on the basis of their historical forecast accuracy.
Panel (c) of Figure (ref) shows the importance of both time-varying weight modifiers. The first thing to notice is that there is much more sparsity in terms of the weight modifiers BPS-RT selects. In the first half of the evaluation sample, we see that features of the individual density forecasts drive the posterior inclusion probabilities. Specifically, we see that the moments of the marginal densities and CRPS, lagged by the forecast horizon $h$, are selected. But in the second half of the evaluation sample, we see the largest proportion of tree splits attributed to the NFCI during and immediately after recessions. The Michigan survey expectations measure also receives more weight after the financial crisis. This is evidence that nonlinear features of BPS-RT are driven by weight modifiers related to the business cycle. In other words, our BPS-RT model finds that the data support changing the combination weights abruptly with business cycle fluctuations. Finally, the time trend receives a higher weight in the post-COVID period of higher inflation seen in $2021$ and $2022$. This finding indicates that this inflationary episode ––– unprecedented within the sample ––– requires a substantial and rapid adjustment of the combination weights. These required weight dynamics cannot be fully captured by the business cycle weight modifiers. Instead, a time trend (or, more precisely, a time dummy) is ideal for modeling such a regime shift from low to high inflation during this exceptional period.
Finally we summarize the properties of the posterior median estimates of the combination weights that are plotted over the evaluation sample in Appendix (ref). We draw out two conclusions for the combination weights estimated when forecasting US inflation one-quarter-ahead (see Figure (ref)). First, BPS-RT places more weight on those component models with SV, especially toward the end of the evaluation sample. This corresponds to the period when BPS-RT outperforms the BPS-RW benchmark (see Figure (ref)).
Second, among these SV models only a subset receives large, in absolute value, weights. This indicates that there is some pay off, in terms of forecast accuracy, to occasionally placing a significantly higher weight on a small subset of models. Interestingly, some models get large negative weights. This amounts to short-selling those models as a “hedge” against the models with higher weights. A roughly similar pattern is seen for the one-year-ahead forecast combination weights seen in Figure (ref).\footnote{Figure (ref) in the online appendix provides additional perspective on the temporal stability of the combination weights by plotting their sum over the evaluation sample. We see that when forecasting US inflation this sum becomes negative during the global financial crisis, indicating how BPS-RT is re-weighting most agents' densities in the face of temporal instabilities. The sum of the weights also spikes upward during the 2021-22 inflationary episode, again indicating how BPS-RT can quickly adapt to temporal change. }
While this subsection has focused on the US inflation application, we end by returning briefly to the EA GDP growth forecasting application. Figure (ref) in Appendix (ref), shows that the combination weights on most individual forecasters from the ECB SPF are, as anticipated given our earlier results, closer to equal than in the inflation application, where there was greater sparsity in the weights. This said, we do still see higher weights on a couple of experts (forecasters 6 and 14). We take this contrasting evidence across the two applications as empirical proof that BPS-RT is sufficiently flexible to adjust to forecasting scenarios that exhibit different dependence structures between the agents' forecasts.
Within the general BPS framework of mcalinn2019dynamic, this paper develops a method for nonparametric density forecast combination using regression trees: BPS-RT. While a handful of papers use nonparametric techniques to combine densities, ours is the first to use regression trees. In contrast to most applications of regression trees we model the coefficients, in our case the combination weights, instead of the variables using the regression trees. We show how this aids interpretation, since the combination model remains linear in the parameters. Additionally, regression trees use covariates, or weight modifiers, to drive changes in parameters, in contrast to conventional BPS applications where model parameters follow a random walk. Taken together, our approach is flexible but retains interpretability through linearity and the use of weight modifiers. We explain how BPS-RT can be used to understand the role of model incompleteness, agent (forecast) clustering, and the time-varying importance of the different weight modifiers.
We test the performance of BPS-RT in two different applications -- combining model-based US inflation density forecasts and subjective histogram-based forecasts of euro area GDP growth. We find that, across both applications, BPS-RT forecasts well in terms of both relative and absolute accuracy. Interestingly, and in contrast to standard BART applications, we find that using a parsimonious single-tree specification outperforms models with more trees. Inspecting the best-performing specification, we observe that this superior performance is due to less disperse forecast densities and BPS-RT's ability to better accommodate the shocks associated with the global financial crisis (in the GDP application) and COVID-19 (in the inflation application). Our proposed measure of model set incompleteness suggests that BPS-RT is able to capture much of the post-COVID rise in inflation. Triggered by a rise in the relative importance of the time trend in determining tree splits, itself highlighting the unusual nature of this inflationary period, BPS-RT also shifts its combination weights toward component models with SV. This contrasts with the prior period of lower inflation, when the business cycle indicators were found to be more important weight modifiers.
Future lines of research could involve investigating, in other forecasting applications and contexts, the usefulness of different sets of weight modifiers and the implications for weight structure. For instance, this could draw on the ability of BPS-RT, via its choice of weight modifiers, to capture general patterns of cross-sectional dependence between competing agents' probabilistic forecasts. Additional structure could be given to the clustering by, for example, letting the combination weight on a given individual agent's density forecast depend not only on characteristics of her own forecast (such as its mean or variance) but on characteristics of the other agents' forecasts.
{\setstretch{0.9} \addcontentsline{toc}{section}{References} }