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.
51,280 characters · 22 sections · 34 citation commands
Bayesian Shrinkage in High-Dimensional VAR Models: A Comparative Study
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 {
\affil[1]{Department of Statistics, UCLA} \affil[2]{Forecasting, Data Science, Airbnb}
\affil[3]{Department of Biostatistics, UCLA Fielding School of Public Health } } \fi
\if10 {
} \fi
Vector autoregressive (VAR) models remain a cornerstone of multivariate time‑series analysis Sims1980, chan2020large, KOOP2013185. A $d$‑dimensional VAR($p$) posits that an observed series $\{\mathbf{y}_t\}_{t=1}^{T}$ satisfies
where $\boldsymbol{\varepsilon}_t\!\sim\!\mathcal{N}(\mathbf{0},\Sigma_\varepsilon)$ is white noise. VARs are ubiquitous in macroeconomics and finance; see stock2002macroeconomic,CrumpEtAl2021_FRBNY,Carriero2022,ZhouChan2023 for recent surveys and methodological refinements.
A VAR becomes non‑stationary when the roots of $(I_d-\mathbf{A}_1z-\cdots-\mathbf{A}_p z^{p})$ lie on or inside the unit circle. Standard remedies include differencing, error‑correction representations, or Bayesian priors that down‑weight explosive parameter draws.
Because a VAR($p$) contains $d^{2}p$ coefficients, dimensionality rises quickly with either $d$ or $p$ Banbura2010,koop2013forecasting,korobilis2019adaptive. Shrinkage—penalising or down‑weighting coefficients toward simpler structures—therefore plays a critical role in stabilising estimation and forecasts stock2002macroeconomic,chan2021minnesota, especially when $d^{2}p \gg T$ bai2022macroeconomic.
\noindentBayesian shrinkage. Local–global priors such as the horseshoe Carvalho2010 or Bayesian lasso Park2008 shrink most coefficients strongly while allowing a few to escape toward their likelihood values. Variants include the spike‑and‑slab george1997variable, hierarchical horseshoes makalic2016simple,pruser2021horseshoe, order‑invariant priors ChanKoopYu2024, dynamic sparsification for hybrid TVP‑VARs Chan2023, and structured factor‑augmented horseshoes ZhouChan2023. Semi‑global or block‑specific shrinkage extends the idea by letting related groups share a common scale parameter GRUBER2025,PrueserBlagov2022. Computational advances now make these priors feasible in very large systems: fast variational Bayes Bernardi2024, approximate Bayes for huge multi‑country VARs Huber2022, and conjugate subspace shrinkage HuberKoop2023. Pandemic‑robust “outlier‑aware’’ priors improve real‑time performance during extreme events CascaldiGarcia2022. Other recent Bayesian shrinkage methods in time series have been explored HuberKoop2023,kowal2019dynamic,KATZ20241556,forecast7030032.
\noindentFrequentist and machine‑learning regularisation. Ridge regression Doan1984 remains popular, while equation‑by‑equation Lasso or elastic‑net penalties Tibshirani1996,zou2005elasticnet,SanchezGarcia2022 yield sparse high‑dimensional VARs with frequentist guarantees under weak dependence Masini2022. Envelope and reduced‑rank ideas provide additional dimension‑reduction tools SamadiHerath2023,CubaddaHecq2022, and non‑linear factor compression has recently been explored Klieber2024. Structured penalties—group Lasso, hierarchical lag shrinkage, SCAD, MCP—continue to improve point forecasts and interpretability NicholsonMattesonBien2017,Basu2019,Nicholson2020bigvar,SongBickel2019,ChenChen2021. New inferential results ensure valid impulse‑response analysis even after sparse estimation Krampe2023.
\noindentForecasting and empirical performance. Time‑varying parameter BVARs with horseshoe or normal‑gamma shrinkage adapt well to structural change BittoFruhwirthSchnatter2019,Feldkircher2024. Factor‑adjusted network VARs (FNETS) separate common and idiosyncratic dynamics to enhance both forecasts and network interpretation Barigozzi2024. Large‑scale empirical comparisons consistently find that local–global or semi‑global priors outperform purely global shrinkage and many frequentist alternatives in density and point forecasts HuberFeldkircher2019,Aprigliano2020,GefangKoopPoon2023.
In this paper, we compare five shrinkage approaches for high-dimensional VAR estimation: three Bayesian priors (horseshoe, lasso, and normal) and two frequentist estimators (ridge and nonparametric shrinkage). We evaluate each method’s handling of over-parameterization in both low- and high-dimensional settings, using root mean squared error, interval coverage, interval length, and RMSE to assess parameter recovery and forecast performance.
We find that local-global priors—particularly the horseshoe—strike a strong balance between parsimony and flexibility, delivering accurate estimation and consistent coverage even in heavily overfitted scenarios. While ridge regression frequently provides competitive point forecasts, it underestimates uncertainty when the parameter space grows large. Meanwhile, nonparametric shrinkage, though computationally efficient, suffers from undercoverage in complex models.
Application to Canadian macroeconomic time series illustrates how each method behaves with various lag choices and a relatively small sample. These empirical findings corroborate the simulation evidence: local-global priors, especially the horseshoe, remain robust even when the chosen lag order exceeds what is strictly necessary.
The remainder of the paper is structured as follows. Section 2 defines the VAR($p$) framework and outlines both Bayesian and frequentist shrinkage estimators. Section 3 describes our simulation designs, metrics, and results. Section 4 provides the real-data application to Canadian macroeconomic variables. Finally, we offer concluding remarks.
We study a $d$-dimensional VAR($p$) of the form
where each $\mathbf{A}_i$ is a $d \times d$ coefficient matrix, and $\boldsymbol{\varepsilon}_t$ is a white-noise process following \[ \boldsymbol{\varepsilon}_t \;\sim\; \mathcal{N}\!\bigl(\mathbf{0}, \,\Sigma_\varepsilon\bigr), \quad \mathrm{Cov}\bigl(\boldsymbol{\varepsilon}_t, \boldsymbol{\varepsilon}_s\bigr) \;=\; \mathbf{0} \text{ for } t \neq s. \] We treat each $\mathbf{y}_t$ as a $d\times1$ column vector.
Vectorizing the Coefficients.\; Let \[ \mathbf{B} \;=\; \bigl[\mathbf{A}_1\;\mathbf{A}_2\;\cdots\;\mathbf{A}_p\bigr] \;\in\;\mathbb{R}^{\,d\times(d\,p)}, \] i.e., the horizontal concatenation of the $p$ coefficient matrices with elements $\beta_j$, $j=1, \dots, d \times dp$. Then define the lagged-regressor vector \[ \mathbf{X}_t \;=\;
\;\in\;\mathbb{R}^{\,d\,p\times 1}, \] so that \[ \mathbf{B}\,\mathbf{X}_t \;\in\;\mathbb{R}^{\,d\times1}. \] The VAR($p$) model in (ref) can thus be written as \[ \mathbf{y}_t \;=\; \mathbf{B}\,\mathbf{X}_t \;+\; \boldsymbol{\varepsilon}_t, \quad \boldsymbol{\varepsilon}_t \sim \mathcal{N}\!\bigl(\mathbf{0}, \,\Sigma_\varepsilon\bigr). \]
In a fully Bayesian treatment, we specify priors for both the coefficient matrix and the error covariance. We gather the coefficients into a matrix \( \mathbf{B}\in\mathbb{R}^{d \times (d\,p)}, \) so that \[ \mathbf{y}_t \;\sim\; \mathcal{N}\!\bigl(\mathbf{B} \mathbf{X}_t^\prime ,\;\Sigma_\varepsilon\bigr), \] where $\mathbf{X}_t$ is the row vector of the lagged responses at time $t$. To ensure $\Sigma_\varepsilon$ is positive-definite, we use a Cholesky-factor parameterization: \[ \Sigma_\varepsilon \;=\; \mathbf{L}\,\mathbf{L}^\top, \quad \mathbf{L} = \text{diag}(\sigma) \,\mathbf{L}_\Omega, \] where $\mathbf{L}_\Omega$ is the Cholesky factor of a correlation matrix with an LKJ prior Lewandowski2009, and each component of $\sigma$ follows a half-Cauchy prior. This flexible structure permits correlation among the $d$ error components.
We then place shrinkage priors on each coefficient in $\mathbf{B}$. Below, we detail three such priors---normal (ridge), horseshoe, and Bayesian lasso---all of which can be combined with the same LKJ-based prior for $\Sigma_\varepsilon$:
\paragraph{Normal Prior (Bayesian Ridge).} A normal (Gaussian) prior imposes a global $\ell_2$ penalty. For each coefficient $\beta_j$ we set priors, \[ \beta_j \;\sim\; \mathcal{N}\!\bigl(0,\,1\bigr), \quad j=1,\ldots,d^2 p, \] so that most coefficients are moderately shrunk towards zero. This parallels the frequentist ridge penalty, and one can include an additional scale factor if stronger or weaker global shrinkage is desired depending on the data set.
\paragraph{Horseshoe Prior.} The horseshoe prior Carvalho2010, makalic2016simple introduces more adaptive shrinkage via a local-global hierarchy. Each $\beta_j$ is modeled as \(\beta_j = B_{\mathrm{raw},j}\,\lambda_j\,\tau\), where \(B_{\mathrm{raw},j} \sim \mathcal{N}(0,1)\), \(\lambda_j \sim \mathrm{C}^+(0,1)\) (local scale), and \(\tau \sim \mathrm{C}^+(0,1)\) (global scale). Small coefficients are heavily shrunk by small local scales, while the heavy-tailed Cauchy priors allow some large signals to remain.
\paragraph{Bayesian Lasso Prior.} Finally, the Bayesian lasso Park2008 imposes a Laplace (double-exponential) prior, \[ \beta_j | \eta \;\sim\; \mathrm{Laplace}(0,\eta), \] which corresponds to an $\ell_1$ penalty in a frequentist setting. As with the horseshoe, this encourages coefficients to be near zero, possibly leading to sparsity in $\boldsymbol{\beta}$.
In all three cases, the covariance matrix $\Sigma_\varepsilon$ is handled by the same LKJ-based prior, thus capturing potential correlations in the innovation process. The posterior distribution factors as \[ \pi\bigl(\bf{B}, \Sigma_\varepsilon \mid \mathbf{y_t}\bigr) \;\propto\; \ell\bigl(\mathbf{y_t} \mid \bf{B}, \Sigma_\varepsilon\bigr)\,\pi(\bf{B})\,\pi(\Sigma_\varepsilon), \] where $\ell$ is the Gaussian likelihood induced by the VAR model, and $\pi(\bf{B})$, $\pi(\Sigma_\varepsilon)$ encode the chosen shrinkage and covariance priors, respectively. By jointly estimating $\mathbf{B}$ and $\Sigma_\varepsilon$, this framework avoids the assumption of uncorrelated errors and allows us to examine how different shrinkage priors influence coefficient estimation in a fully multivariate setting.
\paragraph{How the priors influence estimation.} The three priors differ in the shape and hierarchy of their scale parameters, and these choices translate directly into the amount and selectivity of shrinkage:
These qualitative differences predict the empirical patterns we later observe: ridge delivers the narrowest—but sometimes under‑covering—intervals, the horseshoe attains the best balance of low RMSE and near‑nominal coverage, and the lasso sits in between. Section (ref) returns to this point in light of the simulation and data results.
\paragraph{Ridge Regression.} Classical ridge regression for a VAR($p$) solves
where $\bf{B}$ is just the vectorized collection of $\{\mathbf{A}_i\}$. We set the regularization parameter $\lambda=.1$ in our analysis.
\paragraph{Nonparametric Shrinkage (NS).} We use a James--Stein-like shrinkage approach for VAR coefficients Giannone2015,DelNegro2015, implemented in \verb|R| via \verb|VARshrink| with method="ns". Instead of explicitly solving (ref), the NS method estimates the necessary sample covariances of ${\mathbf{y}_{t-i},\mathbf{y}_t}$ and then applies a closed-form shrinkage rule to these covariance estimates, thereby producing a shrunk solution for $\bf{B}$. In this paper, we rely on the default choice for the shrinkage parameter, which \verb|VARshrink| selects via a moment-based (empirical Bayes) formula akin to Stein’s unbiased risk estimate.
We examine three Monte--Carlo scenarios that differ in dimension ($d$) and in whether the fitted lag order $p$ coincides with the true lag order $p^{\star}$. Table (ref) records the fixed design choices; the recipe that follows is applied independently in each of the $N_{\text{rep}}=50$ replications.
\paragraph{Stationarity margin.} Following the textbook discussion in Lutkepohl2005, we scale the simulated companion matrix so that its spectral radius is comfortably inside the unit circle. Specifically, we divide the initial draw $A_{1}$ by $1.1\,\rho_{\max}$, where $\rho_{\max}$ is the largest eigenvalue (in modulus). This simple rescaling keeps all roots of the characteristic polynomial strictly below one—thereby ensuring covariance–stationarity—while preserving empirically plausible coefficient magnitudes for macro‑VAR applications.
\paragraph{Design Matrix Setup.} To estimate a VAR($p$) in a linear regression framework, we arrange lagged responses into a design matrix \(\mathbf{X}\in\mathbb{R}^{(T_{\mathrm{train}} - p)\times(d\,p)}\). For \(t = p+1,\dots,T_{\mathrm{train}}\), the \(t\)-th row of \(\mathbf{X}\) (denoted \(\mathbf{X}_t^\top\)) is formed by horizontally concatenating the transposes of the \(p\) lagged column vectors \[ \mathbf{y}_{t-1}, \mathbf{y}_{t-2}, \ldots, \mathbf{y}_{t-p} \quad (\text{each } d\times1). \] Hence, each row of \(\mathbf{X}\) is a \(1\times(d\,p)\) vector. Likewise, the \(t\)-th row of the response matrix \(\mathbf{Y}\in\mathbb{R}^{(T_{\mathrm{train}} - p)\times d}\) is simply the transpose \(\mathbf{y}_t^\top\in\mathbb{R}^{1\times d}\). Once \(\mathbf{X}\) and \(\mathbf{Y}\) are formed, any penalized or Bayesian regression method can be applied directly, and the estimated coefficient matrix is then reshaped to match \(\mathbf{B}\in\mathbb{R}^{d\times(d\,p)}\) as defined in Section (ref).
\paragraph{Frequentist fits and block-bootstrapped standard errors.} We estimate the VAR coefficients in Ridge (glmnet with \(\alpha=0\)) and NS (VARshrink with method="ns") by penalized least squares, and then obtain empirical standard errors via a block bootstrap to better respect local time dependence. Specifically we
This procedure retains within-block autocorrelations (up to 4 lags) while randomly mixing which blocks are selected, preserving important time-series structure better than naive row-wise (i.i.d.) resampling. As a result, the resulting intervals yield more realistic coverage for dependent data.
\paragraph{Bayesian fits.} We fit the three Bayesian models by calling Stan with 4 parallel Markov chains, each run for 2000 total iterations (the first 500 of which are warm-up). We fix seed=123 for reproducibility and use $\{\texttt{adapt\_delta}=0.9, \;\texttt{max\_treedepth}=12\}.$ In Stan’s Hamiltonian Monte Carlo (HMC) framework, adapt_delta is the target acceptance probability, and increasing it to 0.9 aims for smaller step sizes and more conservative sampling. The max_treedepth parameter caps the depth of the binary tree in each iteration’s leapfrog integrator, preventing extremely long trajectories.
For each chain, we obtain posterior draws of the coefficient vector \(\boldsymbol{\beta}\). We summarize each coefficient by its posterior mean and 95% central credible interval (2.5% and 97.5% quantiles).
We evaluate each method along two dimensions: parameter estimation and forecast performance.
\paragraph{Parameter estimation.} Let \(\boldsymbol{\beta}_{\mathrm{true}}\) denote the true parameters in $\bf{B}$. Each frequentist method (Ridge or Nonparametric Shrinkage) estimates \(\boldsymbol{\beta}\) by minimizing a penalized least squares criterion, whereas each Bayesian method (Normal/Ridge, Lasso, Horseshoe) uses the posterior mean from MCMC samples as \(\widehat{\boldsymbol{\beta}}\). We then compute the root mean squared error (RMSE), \[ \mathrm{RMSE} \;=\; \sqrt{\frac{1}{d^2 p} \sum_{j=1}^{d^2 p} \bigl(\widehat{\beta}_j - \beta_{j,\mathrm{true}}\bigr)^2}, \] to measure how closely \(\widehat{\boldsymbol{\beta}}\) matches \(\boldsymbol{\beta}_{\mathrm{true}}\).
Next, we construct 95% intervals for each coefficient by applying a block bootstrap to estimate standard errors and forming approximate normal intervals of the form \(\widehat{\beta}_j \pm z_{0.975}\,\mathrm{SE}_j\) for the frequentist approaches, whereas for the Bayesian methods we use the 2.5% and 97.5% posterior quantiles from the MCMC samples. We record the empirical coverage (the percentage of intervals that contain the true value) and the average interval length to assess how well each approach quantifies uncertainty.
\paragraph{Forecasting performance.} To assess predictive accuracy, we reserve the final 20 observations as a test set. Each method then produces sequential one-step-ahead forecasts by estimating \(\mathbf{y}_{t+1}\) at time \(t\) based on all data up to \(\mathbf{y}_{t}\), avoiding the accumulation of multi-step errors. We compute the average root mean squared forecast error (RMSE) \[ \mathrm{Forecast\,RMSE} \;=\; \sqrt{ \frac{1}{20\,d} \sum_{t=T_{\mathrm{train}}+1}^{T_{\mathrm{train}}+20} \left\lVert\mathbf{y}_t - \hat{\mathbf{y}}_t\right\rVert_2^2 }, \] where \(\hat{\mathbf{y}}_t\) is the forecast at time \(t\). After 50 replications per scenario, we summarize the average forecasting RMSE and coverage to compare each method’s predictive capabilities.
The results from the three simulation studies are shown in tables (ref)--(ref) and figures (ref)--(ref). The percentage of replications in which each method has the lowest forecast RMSE or parameter RMSE is shown in table (ref).
\paragraph{Forecasting} The top block of Table (ref) has each method’s mean forecast RMSE. Horseshoe has the smallest value (0.211), followed by ns and Ridge (0.213), Lasso (0.214), and Normal (0.215). Table (ref) has the proportion of replications in which each method has the best forecast: Horseshoe leads with 60%, ns has 20%, Normal 10%, Ridge 6%, and Lasso 4%.
\paragraph{Parameter Estimation} Horseshoe has the lowest overall parameter RMSE (0.0434) and exceeds average nominal coverage (97.2%), with intervals about 8% shorter than those of the next-best method. Lasso (0.0803) and Normal (0.0838) have higher RMSEs but maintain coverage near 94--95%. Both ns (0.0693) and Ridge (0.0730) occupy a middle tier; ns has coverage of 85.7% but yields narrower intervals (mean length 0.204). Horseshoe has the best parameter RMSE in all replications (100%), giving strong shrinkage without sacrificing coverage.
\paragraph{Forecasting} All methods have similar forecasting accuracy (middle block of Table (ref)). Horseshoe has the smallest mean forecast RMSE (0.325), followed by Lasso (0.326) and Normal, ns, and Ridge (0.327). Horseshoe is the top forecaster in 48% of replications, Lasso and ns each in 20%, Ridge in 10%, and Normal in 2% (Table (ref)).
Horseshoe again has the lowest parameter RMSE (0.0536). Lasso, Normal, ns, and Ridge cluster between 0.0568 and 0.0598. Coverage remains high (94--95%) for Horseshoe, Lasso, Normal, and ns, but dips to 84% for Ridge. Table (ref) has results for zero coefficients, where Horseshoe has an RMSE of 0.0357 and 99.0% coverage. The ns method handles zero parameters well but sometimes underperforms on nonzeros. Horseshoe has a nonzero RMSE of 0.0596, higher than Lasso and Normal (0.0571--0.0576), while maintaining overall coverage of 92.9%. Horseshoe has the best parameter RMSE in 90% of replications, followed by Ridge in 10% (Table (ref)).
\paragraph*{Forecasting} In the bottom block of Table (ref), Horseshoe has the smallest mean forecast RMSE (0.342), followed by ns and Ridge (0.365--0.366), Lasso (0.404), and Normal (0.418). Horseshoe is the top forecaster in 100% of replications (Table (ref)), indicating a strong ability to handle overfitting.
\paragraph*{Parameter Estimation} Horseshoe has the lowest parameter RMSE (0.0394), with better-than-nominal average coverage (97.5%) and intervals about 8% shorter than those of the next-best method. Lasso (0.104) and Normal (0.117) have higher RMSEs but maintain nominal coverage (95--96%) through wider intervals (0.432--0.464). Both ns and Ridge have moderate RMSEs (0.0619--0.0635) but show lower coverage (88.2% and 84.5%) and narrower intervals (0.182--0.197). Horseshoe remains the top performer in parameter RMSE for 100% of the replications.
Overall, Horseshoe consistently has excellent forecast accuracy and parameter recovery, including the lowest RMSE, high coverage, and moderate interval lengths. Ridge occasionally has strong forecasts but frequently undercovers in higher dimensions. Lasso and Normal have intermediate performance for both forecasting and parameter estimation, ensuring reasonable coverage by using somewhat larger intervals. The ns approach is computationally efficient and sometimes has precise point estimates, but coverage can be volatile due to overly narrow intervals. These findings reinforce the advantages of local-global shrinkage (Horseshoe) in moderate- and high-dimensional VAR contexts, especially when the lag order is inflated.
We illustrate our methods on the Canada dataset from the R package vars JSSv027i04, which provides quarterly macroeconomic observations on four Canadian variables spanning $T=84$ quarters (1980Q1--2000Q4): employment (e, in log-index form), productivity (prod, in log-index form measuring labor productivity), real wages (rw, in log-index form), and the unemployment rate (\emph{U}, in percent). Economic considerations suggest these variables are jointly dependent, making a Vector Auto Regression (VAR)-based approach suitable.
\paragraph{Differencing and Stationarity.} To reduce nonstationarity, we difference each series once \[ \Delta \mathbf{y}_{t} \;=\; \mathbf{y}_{t} \;-\; \mathbf{y}_{t-1} \quad (t=2,\dots,T). \] We then estimate the VAR on these differenced observations. To obtain forecasts on the original scale, we invert the differencing by recursively summing the predicted differences \[ \widehat{\mathbf{y}}_{T+1} \;=\; \mathbf{y}_{T} \;+\; \widehat{\Delta \mathbf{y}}_{T+1}, \quad \widehat{\mathbf{y}}_{T+2} \;=\; \widehat{\mathbf{y}}_{T+1} \;+\; \widehat{\Delta \mathbf{y}}_{T+2}, \quad \text{etc.} \] This ensures the final forecasts reflect the original scale of the data, while the differencing step helps achieve stationarity and limit spurious trends when fitting the VAR.
\paragraph{Lag–order experiment ($p=1,\dots,12$).} For each shrinkage method we fit twelve separate VAR($p$) models ($p=1,\dots,12$). Every model is trained on the first $T-4$ quarterly differences and produces one‑step‑ahead level forecasts for 2000Q1–Q4. For a given lag order $p$ we compute \[ \mathrm{RMSE}_{p} \;=\; \sqrt{\frac{1}{4}\sum_{h=1}^{4}\bigl(y_{T+h}-\hat y^{(p)}_{T+h}\bigr)^{2}}, \qquad \mathrm{MAPE}_{p} \;=\; \frac{100}{4}\sum_{h=1}^{4} \Bigl|\frac{y_{T+h}-\hat y^{(p)}_{T+h}}{y_{T+h}}\Bigr|. \] Hence each method yields a set of twelve numbers $\{\mathrm{RMSE}_{p}\}_{p=1}^{12}$ and $\{\mathrm{MAPE}_{p}\}_{p=1}^{12}$. We summarize these sets by their means \[ \overline{\mathrm{RMSE}} = \frac{1}{12}\sum_{p=1}^{12}\mathrm{RMSE}_{p}, \quad \overline{\mathrm{MAPE}} = \frac{1}{12}\sum_{p=1}^{12}\mathrm{MAPE}_{p}, \] and by their sample standard deviations \[ \mathrm{SD\,RMSE} = \sqrt{\frac{1}{11}\sum_{p=1}^{12} \bigl(\mathrm{RMSE}_{p}-\overline{\mathrm{RMSE}}\bigr)^{2}}, \quad \mathrm{SD\,MAPE} = \sqrt{\frac{1}{11}\sum_{p=1}^{12} \bigl(\mathrm{MAPE}_{p}-\overline{\mathrm{MAPE}}\bigr)^{2}}. \] These four statistics: $\overline{\mathrm{RMSE}}$, $\overline{\mathrm{MAPE}}$, $\mathrm{SD\,RMSE}$ and $\mathrm{SD\,MAPE}$ provide, respectively, the average forecast accuracy and the across‑lag variability reported in Table (ref).
Overall, Horseshoe achieves the most consistent and accurate forecasts, attaining the lowest mean $ \overline{\mathrm{RMSE}}$ (0.51). Ridge ($ \overline{\mathrm{RMSE}}$ = 0.56) and NS ($ \overline{\mathrm{RMSE}}$ = 0.60) provide somewhat intermediate performance, while Lasso (RMSE = 0.60) and Normal ($ \overline{\mathrm{RMSE}}$ = 0.63) tend to yield slightly larger prediction errors. In terms of variability, Horseshoe’s standard deviation of RMSE (0.22) is comparable to that of Lasso and Ridge, whereas Normal exhibits the largest SD (0.26).
For $ \overline{\mathrm{MAPE}}$, Horseshoe again stands out with an average of 0.71%, followed by NS (0.97%), Ridge (1.06%), Lasso (1.24%), and Normal (1.37%). This pattern suggests that Horseshoe effectively suppresses many small coefficients without overshrinking the larger signals, yielding robust relative‐error forecasts even as \(p\) grows. By contrast, Normal’s wider prior and Lasso’s strong \(\ell_1\) shrinkage can lead to higher MAPE in certain lag settings (see Figure (ref)).
Shrinkage patterns for a few representative lag orders (\(p=3,6,9,12\)) appear in Figure (ref). Horseshoe and Lasso consistently display heavy shrinkage toward zero for small coefficients, Normal has moderate Gaussian-like shrinkage, NS exhibits a broader coefficient spread, and Ridge pulls estimates closer to zero but never to an exact zero. Across the different $p$ values, Horseshoe’s local-global prior structure appears to adapt more flexibly, resulting in better overall forecasting metrics.
To illustrate performance at a higher lag, we examine the VAR(11) specification, which delivers the lowest average forecasting errors among the tested orders. Table (ref) shows each method’s RMSE and MAPE on the final four quarters of the holdout set. Notably, Horseshoe again achieves the best performance on both measures (RMSE = 0.51, MAPE = 0.60%). The next closest method is Ridge (RMSE = 0.61, MAPE = 1.66%) and NS (RMSE = 0.65, MAPE = 1.26%), while Lasso and Normal both exhibit slightly higher errors (RMSE = 0.70--0.78, MAPE = 1.66--1.81%). This example highlights Horseshoe’s ability to preserve large coefficients and aggressively shrink small ones, maintaining strong predictive accuracy even at high lag orders.
The 1-step-ahead forecasts of the observed data across the four holdout quarters are shown in figure (ref). Horseshoe, ns, and Ridge all track the actual values fairly closely. Lasso and Normal lag behind somewhat. Overall, these VAR($11$) results echo our broader simulation findings: Horseshoe’s adaptive local-global prior can maintain strong performance at high lag orders, and Ridge remains reasonably robust as well, whereas ns, Lasso, and Normal can become less accurate or more variable depending on the specific error metric.
For unemployment, the Horseshoe prior cuts the one‑step‑ahead RMSE from 0.63 to 0.51 (–19%), reducing the average error at a 6 % jobless rate from roughly 0.45 pp to 0.36 pp. System‑wide, it lowers the mean RMSE from 0.63 to 0.51 and roughly halves MAPE from 1.37 % to 0.71 % (Table (ref)), delivering noticeably sharper short‑term forecasts for all four Canadian macro‑series.
Relative to the strongest non‑Bayesian competitor, Ridge, Horseshoe still trims about 0.05 RMSE points and cuts MAPE by nearly one‑third. Forecast variability across lag orders remains comparable (SD RMSE 0.22 for Horseshoe vs 0.24 for Ridge; Table (ref)), underscoring that the accuracy gains are achieved without sacrificing stability.
Our simulation results lead to several important takeaways about shrinkage estimation in VAR models under varying dimensionality and lag orders. First, the Horseshoe prior stands out for consistently achieving the lowest parameter RMSE and near-nominal coverage, particularly in the most challenging high-dimensional or overfit scenarios. This local-global prior structure successfully suppresses small coefficients while preserving truly large effects, thereby producing stable estimates and competitive forecasts across the board. By contrast, Lasso and Normal priors often deliver mid-range forecast accuracy and parameter estimation, but they maintain coverage near or above the 95% target, albeit with wider intervals in some cases.
Ridge regression remains effective for forecasting in low- to moderate-dimensional scenarios (where the ratio of parameters to observations is not excessively large), frequently ranking second or third in terms of forecast RMSE. However, it underestimates parameter uncertainty in high-dimensional settings e.g., when the dimension-to-sample-size ratio is particularly large, leading to undercoverage. Similarly, Nonparametric Shrinkage (ns) provides very short intervals and can achieve strong point forecasts, but it exhibits markedly low coverage in the same high-dimensional regimes, suggesting that its narrower intervals are overly optimistic about uncertainty in heavily over-parameterized models.
In the Canadian macro-economic data application, similar patterns emerge: Horseshoe and Ridge each exhibit strong one-step-ahead forecast accuracy, while Lasso, Normal, and ns occasionally lag behind, particularly when the model order is large. Overall, these findings reinforce the benefits of using local-global shrinkage to adapt to large model spaces, especially for practitioners seeking reliable inference and coverage. Frequentist options like Ridge can still perform competitively in lower-dimensional or less overfit settings but risk severe undercoverage when the parameter space grows.
Taken together, these results underscore that when parameter interpretation and interval validity are paramount, Horseshoe or other local-global Bayesian priors are well-suited to handle high-dimensional or inflated-lag VAR models. If short-term predictive performance alone is the principal goal, Ridge can remain attractive, provided one is willing to accept somewhat lower coverage in complex settings. The Lasso and Normal priors offer middle-ground alternatives, balancing coverage and moderate forecasting performance without fully matching Horseshoe’s combination of shrinkage strength and coverage reliability.
All R scripts and Stan model files used in this study are publicly available at \\ \href{https://github.com/harrisonekatz/BayesVAR-SimStudy}{https://github.com/harrisonekatz/BayesVAR-SimStudy}. In particular, the main simulation script var_three_sim_script.R (which orchestrates data generation, frequentist and Bayesian estimation, and result collation) may be found in the repository’s R/ directory. The repository also includes each of the Stan model files (var_normal.stan, \texttt{var_lasso.stan}, \texttt{var_horseshoe.stan}), along with examples illustrating their usage. All results and figures in this manuscript can be reproduced by running the scripts found in that repository.
The authors declare that there are no conflicts of interest.