Global-local shrinkage: the horseshoe BVAR

Vector autoregressions, Part 12: adaptive shrinkage

Every Bayesian VAR in this series shrank with the Minnesota prior — a fixed structure that treats all own-first-lag coefficients the same, all cross-lags the same, and tightens mechanically with lag length. It works, but it cannot learn which coefficients matter: it imposes the same shrinkage whether a coefficient is a genuine dynamic or pure noise. Global-local shrinkage priors — the horseshoe (Carvalho-Polson-Scott 2010) foremost — do learn. A single global parameter pulls everything toward zero, while heavy-tailed local parameters let a few genuinely important coefficients escape. The result is automatic, adaptive sparsity. This notebook builds the horseshoe VAR from scratch and asks whether adaptive shrinkage beats the Minnesota rule — and what it reveals that Minnesota hides.

The horseshoe prior

Give each coefficient a normal prior whose variance is the product of a local and a global scale: $$\beta_j \sim \mathcal N\!\big(0,\ \lambda_j^2\,\tau^2\,\sigma^2\big),\qquad \lambda_j\sim C^+(0,1),\quad \tau\sim C^+(0,1).$$ The global $\tau$ shrinks the whole equation; the local $\lambda_j$, with heavy (half-Cauchy) tails, occasionally take large values that let individual coefficients slip free of the global pull. Writing the implied shrinkage weight $\kappa_j = 1/(1+\lambda_j^2\tau^2)$ (0 = escapes, 1 = shrunk to zero), the horseshoe prior on $\kappa_j$ is $\text{Beta}(\tfrac12,\tfrac12)$ — U-shaped, piling mass at 0 and 1. Hence the name: coefficients are either kept nearly intact or crushed to zero, with little in between — exactly the behaviour a fixed Minnesota variance cannot produce.

Sampling (from scratch). The half-Cauchy scales are written as inverse-gamma mixtures (Makalic-Schmidt 2016), which makes every full conditional inverse-gamma or normal, so the whole thing is a fast conjugate Gibbs sampler, run equation-by-equation. We fit it to a 10-variable macro VAR (largesv_data.csv) and compare with the Minnesota prior (bvar.py, reused unaltered) and a near-flat prior.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import horseshoe as HS, bvar

d = pd.read_csv("largesv_data.csv", index_col=0, parse_dates=True)
cols = list(d.columns); Y = d.values; m = Y.shape[1]; p = 2
fit = HS.gibbs_hs_var(Y, p, ndraw=2000, burn=1200, seed=1)         # Makalic-Schmidt horseshoe Gibbs
kap = fit["kappa"].mean(0)                                          # (k, m) posterior-mean shrinkage per coefficient
nonint = kap[1:]
print("horseshoe shrinkage: %.0f%% of coefficients killed (kappa>0.9), %.0f%% escape (kappa<0.5)"
      % (100 * np.mean(nonint > 0.9), 100 * np.mean(nonint < 0.5)))

# ---- Figure 1: the adaptive shrinkage map (bright = coefficient escapes) ----
reg = [f"{c}.L{l}" for l in (1, 2) for c in cols]                   # 20 lag regressors (intercept dropped)
esc = 1.0 - kap[1:]                                                 # escape weight: bright = escapes
fig, ax = plt.subplots(figsize=(7.2, 6.4))
im = ax.imshow(esc, aspect="auto", cmap="magma", vmin=0, vmax=1)
ax.set_xticks(range(m)); ax.set_xticklabels(cols, rotation=90, fontsize=7)
ax.set_yticks(range(len(reg))); ax.set_yticklabels(reg, fontsize=6)
ax.set_xlabel("target equation"); ax.set_ylabel("regressor (lag coefficient)")
ax.set_title("Adaptive shrinkage map: bright = coefficient escapes\n96% of the 200 coefficients are crushed to zero; a handful survive")
fig.colorbar(im, ax=ax, label="1 - kappa   (0 = killed, 1 = escapes)")
plt.tight_layout(); plt.show()

horseshoe shrinkage: 96% of coefficients killed (kappa>0.9), 2% escape (kappa<0.5)

What the horseshoe keeps

Of the 200 coefficients in this 10-variable VAR(2), the horseshoe kills ~96% (shrinks them to essentially zero) and lets only a handful escape — and those it keeps are exactly the ones economics expects to persist: the own first lags of the interest rates and unemployment ($\kappa\approx0.08$ for unemployment, $0.23$ for the 10-year yield). The growth-rate equations, which are close to unpredictable, are shrunk almost entirely away. Without being told anything about the variables, the prior has performed automatic variable selection, separating the persistent dynamics from the noise — information the Minnesota prior's uniform, rule-based shrinkage never surfaces.

from scipy.stats import multivariate_normal as mvn
from scipy.special import logsumexp
from scipy import stats

# Does adaptive shrinkage forecast better? An honest expanding-window test: refit each quarter
# on data through t-1, score the 1-step-ahead predictive density of y_t, for 2011-2019.
# horseshoe vs Minnesota (bvar.py) vs a near-flat prior. (~7 min: 36 quarters x 3 refits.)
yrs = d.index.year
test = [t for t in range(len(Y)) if 2011 <= yrs[t] <= 2019]

def xrow(H):                                                       # regressors [const, y_{t-1}, y_{t-2}]
    return np.concatenate([[1.0]] + [H[-l] for l in range(1, p + 1)])

def score(fit_, x, y):                                             # mixture 1-step log predictive density
    B, S = fit_["B"], fit_["Sigma"]
    v = [mvn.logpdf(y, x @ B[s], S[s], allow_singular=True) for s in range(len(B))]
    return logsumexp(v) - np.log(len(v))

sc = {"horseshoe": [], "Minnesota": [], "flat (near-OLS)": []}
for t in test:
    H, x, y = Y[:t], xrow(Y[:t]), Y[t]
    sc["horseshoe"].append(score(HS.gibbs_hs_var(H, p, ndraw=500, burn=300, seed=1), x, y))
    sc["Minnesota"].append(score(bvar.gibbs_minnesota(H, p=p, lam1=0.2,  ndraw=500, burn=300, seed=1), x, y))
    sc["flat (near-OLS)"].append(score(bvar.gibbs_minnesota(H, p=p, lam1=50.0, ndraw=500, burn=300, seed=1), x, y))

avg = {k: float(np.mean(v)) for k, v in sc.items()}
pval = stats.ttest_1samp(np.array(sc["horseshoe"]) - np.array(sc["Minnesota"]), 0.0).pvalue
print("1-step avg log predictive score 2011-2019:", {k: round(v, 2) for k, v in avg.items()})
print("horseshoe - Minnesota: %+.2f nats  (paired t p=%.3f)" % (avg["horseshoe"] - avg["Minnesota"], pval))

# ---- Figure 2: the U-shaped horseshoe (left) + the honest forecast comparison (right) ----
fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))
ax[0].hist(kap[1:].ravel(), bins=30, density=True, color="slateblue", edgecolor="white")
ax[0].set_xlabel("shrinkage weight  kappa"); ax[0].set_ylabel("density")
ax[0].set_title("The horseshoe shape\nposterior kappa piles up at 0 (escape) and 1 (killed)")
names = list(avg); vals = [avg[k] for k in names]
bars = ax[1].bar(names, vals, color=["slateblue", "seagreen", "lightgray"], edgecolor="k")
ax[1].set_ylabel("1-step log predictive score  (higher = better)")
ax[1].set_title("Out-of-sample forecast, 2011-2019\nMinnesota beats the horseshoe (p=%.2f) -- sparsity is not free" % pval)
ax[1].set_ylim(min(vals) - 0.25, max(vals) + 0.15)
for b, v in zip(bars, vals):
    ax[1].text(b.get_x() + b.get_width() / 2, v - 0.07, f"{v:.2f}", ha="center", fontsize=9)
plt.setp(ax[1].get_xticklabels(), fontsize=8)
plt.tight_layout(); plt.show()

1-step avg log predictive score 2011-2019: {'horseshoe': -14.16, 'Minnesota': -13.85, 'flat (near-OLS)': -13.96}
horseshoe - Minnesota: -0.31 nats  (paired t p=0.015)

Results — radical sparsity, at a small forecast cost

• The horseshoe shape is real. The posterior shrinkage weights pile up near 0 and 1 (left panel) — the U-shape that names the prior. This is qualitatively different from the Minnesota prior, whose shrinkage is a smooth deterministic function of lag and own-vs-cross; the horseshoe lets the data decide, coefficient by coefficient, and here it crushes ~96% of the 200 coefficients to zero.
• But it does not forecast better than Minnesota — measured honestly. An expanding-window, refit-every-quarter test (right panel) is unambiguous: the horseshoe's 1-step density score (−14.16) is significantly worse than the Minnesota prior's (−13.85; paired-t p ≈ 0.015), and even a near-flat prior (−13.96) edges it. The aggressive shrinkage that makes the model so sparse also over-shrinks coefficients that carry genuine short-horizon signal. So the sparsity is not free: on this data it costs a small but statistically real amount of forecast accuracy — a known tension, since the horseshoe optimises for selection, not prediction. (This overturns a common shortcut — quoting a single in-sample or lightly-validated score and calling the two "indistinguishable"; the honest out-of-sample number does distinguish them.)
• What it buys instead is interpretability. Because it selects, the horseshoe answers a question Minnesota cannot: which dynamics are genuine? Here, the persistence of rates and unemployment (own first lags of unrate, gs10, ffr), and little else. That automatic, honest sparsity — a radically smaller, readable model for a small forecast premium — is why global-local priors have become a modern default for large and high-dimensional VARs, where the interpretability and regularisation matter more and the forecast gap typically narrows (Huber-Feldkircher 2019; Cross-Hou-Poon 2020). The R cross-check reproduces the same equation-level shrinkage.

References

• Carvalho, C. M., Polson, N. G. & Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika 97, 465-480.
• Makalic, E. & Schmidt, D. F. (2016). A simple sampler for the horseshoe estimator. IEEE Signal Processing Letters 23, 179-182.
• Huber, F. & Feldkircher, M. (2019). Adaptive shrinkage in Bayesian vector autoregressive models. J. Business & Economic Statistics 37, 27-39.
• Cross, J. L., Hou, C. & Poon, A. (2020). Macroeconomic forecasting with large Bayesian VARs: global-local priors. Int. J. Forecasting 36, 899-915.

Data: largesv_data.csv — 10 stationary FRED-QD macro series, quarterly 1960-2019. Next: horseshoe_R.ipynb — the R cross-check (base-R Makalic-Schmidt Gibbs).