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.
44,479 characters · 16 sections · 56 citation commands
Minnesota BART
\onehalfspacing
\newtheorem*{remark}{Remark}
\newtheorem{assumption}{Assumption}
\newtheorem{lemma}{Lemma} \newtheorem{example}{Example} \newtheorem{theorem}{Theorem} \newtheorem{proposition}{Proposition} \newtheorem{definition}{Definition} \newtheorem{corollary}{Corollary} \newenvironment{assumptionprime}[1] { \addtocounter{assumption}{-1}
}
Keywords: Bayesian non-parametrics, non-linear vector autoregressions, shrinkage prior, forecasting.
\onehalfspacing
This paper introduces a multivariate Bayesian additive regression tree (BART) model for macroeconomic forecasting that incorporates structured priors to allow for variable selection and account for time dependency. We extend the standard BART framework to allow for a sparse vector autoregressions (VARs) estimation by introducing a sparsity-inducing Dirichlet hyperprior on the regression tree’s splitting probabilities. This allows for automatic variable selection, reducing overfitting and improving computational efficiency in large-scale models.
Additionally, we propose a structured shrinkage prior that decreases the probability of splitting on higher-order lags, aligning with the Minnesota prior’s principles. This addresses a fundamental limitation of existing BART-based VAR models, which fail to incorporate economic constraints on lag selection and ignore temporal dependencies. We also analyze how different levels of our shrinkage parameter affects the splitting probabilities, results demonstrate that higher values of the parameter lead to a more gradual decay in posterior inclusion probabilities, preserving the influence of lags and cross-lags for a longer range. This choice also can affect the forecasting performance of the model. By integrating these two enhancements, our approach preserves the flexibility of BART while imposing meaningful structure, leading to improved interpretability and forecasting accuracy.
We evaluate our model through an empirical application to U.S. macroeconomic forecasting. Compared to standard BART and Bayesian VAR (BVAR) models, our approach improves predictive accuracy, particularly in capturing time-dependent relationships and higher-order moments of the predictive distribution. In particular, we show that our method enhances density forecasts for key macroeconomic variables such as the Federal Funds Rate and inflation. These findings highlight the potential of structured nonparametric methods for macroeconomic forecasting.
Vector autoregression (VAR) models have been widely used for forecasting and structural analysis of macroeconomic variables (doan1984forecasting; litterman1986forecasting; banbura2010large; koop2013forecasting; carriero2019large; kastner2020sparse). However, as the number of time series included in the model increases, the number of parameters grows quadratically, leading to concerns about overparameterization and in-sample overfitting. To address these challenges, the Bayesian literature on VARs has developed various shrinkage prior specifications.
Despite these advancements, most Bayesian VAR models still assume a linear relationship between endogenous variables and their lags. While macroeconomic relationships are often stable over time, allowing a linear approximation to fit the data reasonably well, this assumption can break down during shocks that alter the economy’s dynamics (huber2023nowcasting). Failing to account for such events can result in poor out-of-sample performance and misinterpretation of impulse response functions.
More recently, nonparametric approaches such as Bayesian additive regression trees (BART) have gained attention as a flexible alternative (chipman2010bart). BART uses regression trees as weak learners, allowing for complex relationships to be modeled without strong parametric assumptions. However, existing frameworks (huber2022inference; clark2023tail) do not accommodate high-dimensional settings or account for time dependency in the construction of regression tree priors. Our paper directly addresses these gaps.
The remainder of the paper is structured as follows. Section 2 introduces the multivariate BART model and provides a necessary introduction to the BART framework. Section 3 develops the prior construction. Section 4 details the prior setup and posterior computation. Sections 5 and 6 present our empirical results: Section 5 discusses the dataset and provides in- and out-of-sample model evidence, while Section 6 focuses on macroeconomic forecasting results. The final section summarizes and concludes.
Define ${\bm{y}}_t = \left(y_{1t}, \dots, y_{nt}\right)'$ as a vector of endogenous variables of dimension $n\times 1$. Define ${\bm{x}}_t = \left({\bm{y}}_{t-1}', \dots, {\bm{y}}_{t-p}'\right)'$ a $k(= np)$ dimensional vector of covariates, where $p$ is the number of lags. We also define ${\bm{G}}({\bm{x}}_t)=\left(g_1({\bm{x}}_t), \ldots, g_n({\bm{x}}_t)\right)'$ as a n-dimensional vector of non-parametric functions :
Each function in the vector ${\bm{G}}(\cdot)$, where ${\bm{G}}: \mathbb{R}^{k}\rightarrow \mathbb{R}^{n}$, hence $g\left({\bm{x}}_t\right): \mathbb{R}^{k} \rightarrow \mathbb{R}$ will be approximated by a Bayesian Additive Regression Trees (BART) model, which is discussed in detail in Section (ref).
Research on large Bayesian vector autoregressions (VARs) shows evidence that stochastic volatility specifications are well supported by the data (carriero2016common; koop2013forecasting; chan2020large) and achieving precise density forecasts. Recent work by chan2023comparing demonstrates that the factor stochastic volatility and Cholesky stochastic volatility specifications outperform the standard common stochastic volatility model.
Therefore, following huber2022inference, the conditional covariance structure is specified as factor stochastic volatility (FSV) (pitt1999time; aguilar2000bayesian; chib2006analysis; lopes2007factor; kastner2020sparse). More precisely, the error term is decomposed as:
where ${\bm{f}}_t = \left( f_{1,t}, \dots, f_{r,t} \right)$ is a $r\times1$ vector of latent factors and ${\bm{\Lambda}}$ is the associated $n\times r$ factor loading matrix. This factor specification is not identified. The latent factors and the idiosyncratic errors are assumed to be independent and jointly Gaussian:
where ${\bm{\Omega}}_t = diag\left( e^{h_{1,t}}, \dots, e^{h_{n,t}} \right)$ and ${\bm{H}}_t = diag\left( e^{h_{n+1,t}}, \dots, e^{h_{n+r,t}} \right)$ are diagonal matrices. The evolution of the log-variance process for $i = 1, \dots, n+r$ is defined as:
for $t =2, \dots, T$. For $t =1$, we assume a stationary distribution $h_{i,1} \sim \mathcal{N}\left( \mu_i, \frac{\sigma_{i}^{2}}{1-\phi_i^2}\right)$. The number of factors are defined using an upper bound as in aguilar2000bayesian. However, in most practical applications, a small number of factors is sufficient to capture the dynamics of the covariance structure, as seen in bolfarine2024decoupling; fruhwirth2024sparse.
To address this, we impose a shrinkage prior on the columns of the ${\bm{\Lambda}}$ matrix, which pushes irrelevant factors toward zero. For this purpose, we adopt the horseshoe prior proposed by carvalho2010horseshoe. Conditional on the latent factors, this non-parametric VAR becomes $n$ unrelated regressions, estimated equation-by-equation.
Rewriting the model in terms of full-data matrices we have:
where ${\bm{X}} = \left({\bm{x}}_1, \dots, {\bm{x}}_T\right)'$ and ${\bm{Y}} = \left({\bm{y}}_1, \dots {\bm{y}}_T \right)'$, is a $T \times k$ and $T \times n$ matrix respectively.
Our choice estimating ${\bm{G}}:\mathbb{R}^{k}\rightarrow \mathbb{R}^{n}$ is using a sequence of nonparametric decision tree ensemble such that:
where $m$ is the number of trees in the ensemble. The general idea is to aggregate individuals “weak learners" into a unified “strong learner". Each regression tree $g_{j,m}(\cdot)$ defines a piecewise constant function based on the arrangement of split rules $\mathcal{T}_{jm}$, associated with a $b_{jm}$-dimensional vector $\mathcal{M}_{jm} = \left(\mu_{m,1}^{(j)}, \dots , \mu_{m,b_{jm}}^{(j)}\right)'$ of terminal nodes coefficients where $b_{jm}$ is the number of leaves per tree $m$ in equation $j$.
There is established evidence in the literature that regression trees have demostrated strong empirical performance in a wide variety of contexts, including supervised learning (chen2016xgboost; ke2017lightgbm), casual inference (hahn2020bayesian), density regression (orlandi2021density). For more discussion see grinsztajn2022tree and hill2020bayesian.
The model prior structure follows chipman2010bart. Conditioned on the model hyperparameters, BART priors are :
The prior distribution for the trees $p_{\mathcal{T}}$ consists of two components. First a prior on the shape of the tree $\mathcal{T}$ and second a prior on the splitting rules $[x_{q} \leq C_{q}]$ for each branch node of the tree. The prior on the tree structure includes a probability that a node of depth $d$ is “terminal" or does not split. Starting from a root, the node will split with probability probability $\gamma(1 + d)^{-\beta}$. If that is not the case, the root is a terminal node. This iterates until all nodes at certain depth are terminal. We follow the convention introduced by chipman2010bart by taking $\gamma = 0.95$ and $\beta = 0.2$.
If a node is split, the split rule, defined by a variable and a cutpoint, is sampled as follows. A variable index $q \in 1, \dots, k$ is sampled according to $q \sim Categorical\left({\bm{s}}\right)$, such that ${\bm{s}}$ is probability vector where $s_{q} = \frac{1}{k}$. Subsequently, a cutpoint $C_q$ is sampled uniformly on the observed range of values of variable $x_{q,t}$ at the current node. Finally, for each terminal node $b_{jm}$ in the tree, we draw a mean parameter $\mu_{m,b_{jm}}^{(j)} \sim \mathcal{N}\left(0, \sigma_{\mu}^{2}/ M\right)$.
Recent advancements in the non-parametric Vector Autoregression (VAR) literature, particularly the integration of Bayesian Additive Regression Trees (BART) to extend the traditional VAR framework, provide compelling evidence of improved forecasting performance, as demonstrated by huber2022inference. However, these recent developments do not adequately address scenarios in which the true data-generating process (DGP) is sparse. Additionally, the original BART prior used in huber2022inference does not account for the temporal dependency structure present in macroeconomic data. Our contributions to the literature involve addressing these shortcomings within the vector autoregression framework.
A well-documented finding in macroeconomic forecasting—from univariate models for GDP, inflation, and interest rates to multivariate formulations is that simple random-walk or near-unit-root forecasts often perform reasonably well, particularly over short-to medium-term horizons, as initially shown by nelson1982trends. This persistence aligns with the near-unit-root behavior frequently observed in economic time series.
Therefore, constructing a prior that embeds a random-walk assumption—where, in the absence of data, the best initial estimate is that today’s value equals yesterday’s value plus some small drift or shock—provides a natural way to incorporate this empirical regularity into the model.
The Minnesota prior, introduced by litterman1980bayesian and doan1984forecasting, is a shrinkage prior specifically designed to mitigate the issue of overparameterization commonly encountered in large linear VAR models. The design for this prior is based on empirical evidence on macroeconomic time-series behavior, which suggests that variables often exhibit persistent dynamics that can be effectively captured through structured shrinkage. It integrates several reasonable assumptions, including cross-variable shrinkage, where coefficients for lags of different variables are reduced more significantly than those for the variable's own lags. Additionally, it reflects the belief that higher-order lags contribute less to forecasting accuracy. For further details on the Minnesota prior, see kadiyala1997numerical, karlsson2013forecasting, and chan2020large.
While the traditional BART framework has shown promise in time series analysis, It employs a uniform prior on splitting variables when sampling split rules—an assumption that is often unrealistic in multivariate forecasting settings. To address this limitation, we leverage the DART prior framework, presented by linero2018high, to incorporate time dependency into our trees, drawing on insights from prior specification in the linear Vector Autoregression (VAR) literature.
By integrating the principles of the Minnesota prior into this framework, our approach introduces structured shrinkage that respects the temporal dependence inherent in macroeconomic data. This allows the estimation of large dynamic systems within a multivariate BART framework while preserving interpretability and forecasting accuracy.
For equation $n$, the prior for the splits probability is defined:
The scale parameters of the Dirichlet distribution are defined are defined as follows:
$$ \phi_{in} =
$$
where $\rho = \frac{\sigma_{i}^{2}}{\sigma_{j}^{2}}$ represents the variance ratio of an AR($p$) process for each variable. This formulation induces a smooth shrinkage effect on the prior probabilities of splits. The choice of $ \lambda_j$, where $j =1,2$, determines the rate at which these probabilities decay. For our analysis, we set the hyperparameters to fixed values, specifically \( \lambda_1 = 1 \) and \( \lambda_2 = 0.5 \). While \( \lambda_1 \) and \( \lambda_2 \) can be estimated directly from the data, we choose to fix them to ensure analytical tractability and simplify the estimation process. A detailed discussion on prior elicitation is provided in Section (ref).
However, this configuration does not lead to a sparse solution. It is possible to elicit a prior that explicitly favors a more parsimonious model with a smaller number of covariates. Conditional on the tree topology, and for a fixed $\lambda$, substituting $l^2$ with $k$, we obtain:
$$ \left(s_{1n}, \ldots, s_{kn}\right) \sim \operatorname{Dirichlet}\left(\frac{\lambda}{k}, \ldots, \frac{\lambda}{k}\right) $$
The parameter \( \lambda \) governs the degree of sparsity introduced in the tree function. As demonstrated by linero2018high, when both the number of predictors (\( k \)) and the number of branches (\( B \)) in the ensemble are large, the prior distribution of the number of relevant predictors (\( Q - 1 \)) can be approximated by a Poisson distribution with parameter \( \theta_B \), where \( \theta_B = \lambda \sum_{i=0}^{B-1} (\lambda + i)^{-1}\). The level of sparsity can be controlled by setting \( \lambda \) to a predetermined value. Under a fully Bayesian parameter selection framework, we follow the suggested approach of assigning \( \lambda \) a hyperprior, specifically \( \lambda/(\lambda + k)\sim \text{Beta}(0.5,1) \).
Figure (ref) illustrates how the value of $\lambda$ corresponds to different levels of sparsity. This approach may introduce some rigidity, but it aligns more closely with our objective of transparently examining the trade-off between sparsity across variables and to a variable's own lags in tree models. In Section (ref), we discuss in detail how $\lambda$ affect an out-of-sample exercise.
The model is estimated using a combination of traditional Bayesian inference techniques commonly employed in the VAR and BART literature. Tree sampling updates are performed using a Metropolis-Hastings (MH) algorithm, as proposed by chipman2010bart, while most of the remaining steps leverage closed-form Gibbs Update. The conditional posteriors for the factor loadings and factors follow well-known Gaussian distributions. For the stochastic volatility components, present in both the factors and idiosyncratic innovations, we employ the efficient sampler outlined in kastner2014ancillarity. In cases where a homoskedasticity assumption is applied, the model uses the traditional inverse-gamma prior.
As mentioned previously, conditioning on the covariance structure we can estimate the VAR equation-by-equation. The model can be rewritten the as a system of $n$ independent equations. Let ${\bm{Y}}_{\bullet j}$ denote the $j\text{-th}$ column of the matrix $Y$ as defined in (ref). Therefore we can write our dynamic system as :
where ${\bm{\Lambda}}_{\bullet j}'$ is the $j\text{-th}$ column of the factor loading matrix. This formulation reveals that our model is a generalized additive model, where the forest component approximates the relationship of ${\bm{y}}_t$ with its lags, while a shared linear component across all equations captures the relationships among the variables.
The two main departures of the sampling strategy proposed by chipman2010bart is regarding the partial residuals definition, that needs to take into account the factors and loadings structure and an additional step to update the vector of split probabilities ${\bm{s}}$. To sample the trees using Bayesian backfitting as in the likelihood function depends of $\left(\mathcal{T}_{jm}, \mathcal{M}_{jm}\right)$ through the partial residuals that should be defined for our case as:
Therefore we can sample tree structure marginalizing over $\mathcal{M}_{jm}$, such that:
can be obtained in a closed form solution up to a constant. Allowing to carry out each draw from $\left(\mathcal{T}_{jm}, \mathcal{M}_{jm} |{\bm{R}}_{jm}, \sigma_{j,t}\right)$ sequentially.
To draw the probability split vector ${\bm{s}}$, we follow linero2018high and leverage the conjugacy between the Dirichlet prior and multinomial sampling, enabling a full-conditional Gibbs update given by:
where $\phi$ represents the shape parameter of the Dirichlet and Minnesota specification, and $m_k$ denotes the number of splitting rules for predictor $k$ in the ensemble.
The factor loadings ${\bm{\Lambda}}$ are drawn from a full conditional distribution that follows a Gaussian distribution in a standard form. For each row of ${\bm{\Lambda}}$, denoted as ${\bm{\Lambda}}_i$, we sample as follows:
\[ \bar{{\bm{W}}}_i = \left({\bm{F}}_i^\prime {\bm{F}}_i +{\bm{W}}_i^{-1}\right)^{-1}, \]
\[ \bar{{\bm{L}}}_i = \bar{{\bm{W}}}_i \left({\bm{F}}_i^\prime \tilde{y}_i\right). \]
Here, ${\bm{F}}_i$ is the $t$-th row ${\bm{f}}_t / e^{h_it / 2}$, and the $t$-th element of $\tilde{y}_i$ is given by $\left(y_{it} - f_i(x_t) \right) / e^{h_it / 2}$. The matrix ${\bm{W}}_i$ is a prior variance-covariance matrix of dimension $(n r) \times (n r)$. Since the number of factors are determined by an upper bound, we sample from a horseshoe prior using the auxiliary sampler proposed in makalic2015simple for each column of ${\bm{\Lambda}}$. The factors are generated on a $t\text{-by-}t $ basis using Gaussian distributions as in aguilar2000bayesian.
We conduct a forecasting exercise using multivariate Bayesian additive regression trees to compare our proposed sparse and Minnesota-BART priors with the baseline BART prior structure from huber2022inference
We use a dataset consisting of 22 U.S. quarterly variables covering the period from 1965Q1 to 2019Q4. The data is sourced from the FRED-QD database at the Federal Reserve Bank of St. Louis, as described in mccracken2016fred. The dataset includes a range of standard macroeconomic and financial variables, such as real GDP, industrial production, inflation rates, labor market indicators, and interest rates. These variables are transformed to achieve stationarity, typically by computing growth rates. We will include 13 lags of the endogenous variables in our model. A detailed description of the variables and their transformations can be found in Appendix (ref).
We begin by presenting predictive evidence of the effectiveness of our proposed priors and their impact on predictive accuracy. To achieve this, we construct a recursive forecasting design, using 1965Q1 to 2004Q4 as the initial training period. We employ an expanding window strategy: after performing $h$-step-ahead forecasts, we incorporate the next observation into the dataset and re-estimate the model, obtaining a new draw from the predictive density with the updated information. This process is repeated iteratively until all available data has been utilized. The one-step-ahead predictive distribution is given by
where ${\bm{y}}^{t}$ represents the historical time series up to time $t$, i.e., ${\bm{y}}^{t} = \left(y_1, \dots, y_{t}\right)$. The parameter ${\mathbf{\boldsymbol{\vartheta}}}$ encapsulates all unknown parameters of the model. This integral is solved using the standard Monte Carlo approach:
$$ p\left({\bm{y}}_{t+1} | {\bm{y}}^{t}\right) \ \approx \, \frac{1}{M} \sum^M_{m=1} p({\bm{y}}_{t+1}|{\bm{y}}^{t},{\mathbf{\boldsymbol{\vartheta}}}^{(m)}), $$
where the one-step-ahead predictive density $p\left({\bm{y}}_{t+1} | {\bm{y}}^{t}, {\mathbf{\boldsymbol{\vartheta}}} \right)$ is Gaussian, conditional on knowing ${\mathbf{\boldsymbol{\vartheta}}}$. Therefore, for each draw ${\mathbf{\boldsymbol{\vartheta}}}^{(m)}$ from the posterior distribution $p\left({\mathbf{\boldsymbol{\vartheta}}} | {\bm{y}}^{t}\right)$, we obtain the predictive density:
$$ {\bm{y}}_{t+1}|{\bm{y}}^{t}, {\mathbf{\boldsymbol{\vartheta}}}^{(m)} \sim \mathcal{N}\left({\bm{F}}^{(m)} \left({\bm{X}}_{t+1}\right), \Sigma_{t+1}^{(m)}\right), $$
where ${\bm{F}}^{(m)}(\cdot)$ is generated by our tree sampling algorithm, and $\Sigma_{t+1}^{(m)}$ is drawn using the covariance structure specified in Eq. (ref). For the stochastic volatility (SV) specification, as outlined in Eq. (ref), the forecasts of $\Omega_{t+1}$ and $H_{t+1}$ are obtained as follows. Given the posterior draws of $h_{it}^{(m)}$, we simulate $h_{it+1}^{(m)}$ from a conditional normal distribution with mean $\mu^{(m)}_{i}+\phi^{(m)}_{i}(h^{(m)}_{i,t} -\mu^{(m)}_{i})$ and variance $\sigma^{2(m)}_{h}$. Higher-order forecasts are computed recursively. In the following section, note that ${\bm{y}}_{t+1}$ corresponds to a one-quarter-ahead prediction, i.e., three months. We use the Gibbs sampler described in the section (ref) to obtain 5,000 posterior, after a burn-in period of 30,000.
We use a standard small-scale Minnesota BVAR with stochastic volatility (SV) as a benchmark to evaluate our model's performance across three key variables: real GDP growth (GDPC1), inflation (CPIAUCSL), and the Federal Funds Rate (FEDFUNDS). To facilitate this comparison, we compute the (relative) root mean square prediction error (RMSPE). A value below one indicates that the model outperforms the benchmark, while a value above one suggests weaker performance. Importantly, we account for heteroskedasticity in each prior specification. A priori, we expect that the combination of our proposed priors and time variation in volatilities will enhance point forecasts. The results are summarized in Figure (ref).
Our findings reveal a consistent pattern: predictive accuracy improves as the forecasting horizon extends. Compared to the linear baseline, incorporating a richer information set and a non-linear approach enhances point forecasts. The SV variants generally improve precision relative to their homoskedastic counterparts. However, the numerical differences are often small, particularly when comparing sparse and non-sparse versions of the prior. Importantly, for inflation forecasting, the introduction of stochastic volatility (SV) leads to a more pronounced improvement in predictive accuracy for the sparse specification compared to its non-SV counterpart. Regarding prior choice, with the exception of real GDP growth, a smoother shrinkage approach in the variable splits tends to outperform its sparse alternative for the variables of interest. This suggests that, for certain variables, a prior that imposes gradual shrinkage—leveraging more information compared to a sparse prior—is preferable.
Although point forecast exercises are important, they provide only a partial assessment of model performance. To obtain a more comprehensive evaluation, it is essential to consider additional metrics that account for the model's ability to predict higher-order moments of the predictive distribution for the variable of interest. Therefore, is necessary to utilize a metric to access the accuracy of density forecasts. As discussed by geweke2010comparing, log predictive density scores (LPDS) are often used to compare different models. In this paper, we use the LPDS as a metric to evaluate and compare the performance of various BART prior specifications. We will make slightly change of notation when presenting this calculations. The first $t_0$ time series observations, ${\bm{y}}^{tr} = (y_1, \ldots, y_{t_0})$, are designated as the “training sample," while the remaining observations, $y_{t_0+1}, \ldots, y_T$, are used for evaluation based on the log predictive density:
In Equation ((ref)), $p({\bm{y}}_t \mid {\bm{y}}^{t-1})$ represents the one-step-ahead predictive density for time $t$. This density is evaluated at the observed value ${\bm{y}}_t$. Note that this framework not only works when evaluating the joint performance for higher order moments of the predictive distribution for the multivariate model, but also the marginal density scores, i.e, we have the LPDS for the three variables of interest as well.
Figure (ref) shows that while our prior does not perform as well in point forecasting for GDP, it excels in predictive distribution forecasting. The Minnesota prior achieves the best performance, with its FSV specification closely following. Notably, this prior alternative not only outperforms the basic BART prior but also surpasses its simpler linear counterpart.
For the FEDFUNDS variable, the results follow a different pattern. As shown in Figure (ref), the more sparse specification underperforms in point forecasting when compared to the linear benchmark. However, it dominates across all horizons in predictive distribution forecasting. The smooth shrinkage prior emerges as the second-best choice.
For the CPI variable, the results favor the baseline prior with an SV correction, except for the first forecasting horizon, as shown in Figure (ref). Even in this case, our smooth shrinkage option remains a close second in forecasting the predictive density distribution.
This results suggest that, for both point forecasts and marginal density, the alternative prior structures proposed in this paper are competitive and, in some cases, even outperform the baseline BART model. Figure (ref) shows that when evaluating the (joint) predictive density, our alternatives, which incorporate varying levels of shrinkage, consistently outperform the basic prior, particularly when heteroskedasticity is accounted for, as initially hypothesized. Notably, the Minnesota specification outperforms the sparse alternative across all time horizons.
Our prior modifications introduce key advantages to multivariate BART analysis. In the baseline BART framework, variable importance is typically assessed by counting the number of times each feature appears in a splitting rule. However, this approach has limitations. Since the traditional prior imposes a uniform split probability across all features, variable importance can only be inferred if there is a trade-off between the number of trees and interpretability. In other words, improving interpretability often comes at the expense of predictive performance chipman2010bart; bleich2014bayesian; linero2018high. This limitation does not apply to our proposed prior structures. To illustrate, Figure (ref) presents the Posterior Inclusion Probability (PIP) for each prior choice. The PIP is defined as: $$ \text{PIP}_j = P(\text{predictor } j \text{ appears in the ensemble} | \text{Data}) $$
We present the PIP for the in-sample results of our model before applying the expanding window, focusing on the Inflation variable. As shown in Figure (ref), the BART prior specification distributes splits relatively evenly across all features in the model. In contrast, the sparse prior specification shrinks the split probability of most variables to near zero while preserving the variable’s own first lag as the most significant predictor. This aligns with the Bayesian VAR literature, which finds that the AR(1) term explains the largest share of variation in the variable of interest. For the Minnesota prior, the expected decay in PIP follows a structured lag hierarchy, preserving the anticipated importance ordering among the split variables.
As previously discussed, the choice of \( \lambda \) is of critical importance, as it plays a central role in determining the expected level of shrinkage in the model. One approach we have presented is to select \( \lambda \) based on a shrinkage target informed by the econometric literature and subject-matter considerations. To further investigate its impact, we analyze different levels of \( \lambda \), specifically considering a grid of values: \( \lambda_1 = \{1, 3, 5, 10, 20\} \) and \( \lambda_2 = \{0.5, 1, 1.5, 2.5, 5, 10\} \). We then examine its effects on the log-predictive density score in comparison to the standard BART prior.
First, we examine how the choice of the smooth shrinkage parameter influences the “hyperbolic" shape of the posterior inclusion probability (PIP) for the inflation variable in the in-sample analysis. As shown in Figure (ref), an increase in \( \lambda \) values results in a slower decay of PIP, reflecting that the shrinkage of both lags and cross-lags does not reach zero as quickly as with our reference values of \( \lambda_1 = 1 \) and \( \lambda_2 = 0.5 \).
Additionally, when comparing different PIPs for CPI’s own lag, we observe the influence of the hyperparameter choice on the rate at which the inclusion probability decays, as illustrated in Figure (ref). To compare different prior configurations we will use the log predictive density score, calculated as how was presented in the sections before. Figure (ref) shows that besides any configuration that you choose for the prior, it is clear that for accurate density forecasts, its necessary to take into account time dependency when constructing the prior.
Overall, our analysis highlights the crucial role of \( \lambda \) in shaping the degree of shrinkage in the model and its impact on both variable selection and predictive performance. The results demonstrate that higher values of \( \lambda \) lead to a more gradual decay in posterior inclusion probabilities, preserving the influence of lags and cross-lags for a longer range. These findings reinforce the importance of carefully selecting the shrinkage parameter, as it directly affects model interpretability and forecasting accuracy.
The classic BART prior for split rules, which samples split variables with uniform probability, implicitly assumes that the mean function is not sparse or time-dependent in its input features, limiting inference on variable importance and hindering effective high-dimensional analysis. This paper introduces a framework that integrates insights from the literature on sparse priors for BART and Bayesian VARs. The proposed model allows for shrinkage in split probabilities, enabling the estimation of large dynamic systems within a multivariate BART framework, unlike existing approaches. Additionally, we demonstrate that incorporating a prior akin to the linear Minnesota prior introduces smooth shrinkage and time dependence information during prior elicitation.
Our approach was illustrated using a large U.S. dataset, where we showed that the proposed priors yield substantial improvements in forecast accuracy, particularly for higher-order moments. In this context, the Minnesota specification appears to be especially effective in extracting forecasting gains by introducing a smoother shrinkage approach compared to the sparse alternative. While point forecast accuracy varies across the variables of interest, nonlinear models often outperform linear specifications, their gains are substantial when they do, whereas losses to linear models are relatively small. Moreover, our prior structure demonstrates predictive gains for key macroeconomic variables, such as the Federal Funds Rate and inflation, highlighting its practical relevance for economic forecasting.
Although our model focuses on a reduced-form specification for forecasting, the framework can also be adapted for structural analysis using techniques such as Generalized Impulse Response Functions (GIRFs) (koop1996impulse) or Local Projection (LP) estimation (jorda2005estimation). Future work could explore alternative sampling methods or algorithmic optimizations to enhance scalability and reduce computational costs. Furthermore, expanding this approach to account for richer structural dynamics, such as incorporating time-varying parameters or state-dependent effects, could significantly enhance its applicability in macroeconomic modeling.