SSVS for VARs — searching restrictions on a vector autoregression

Bayesian variable selection, Part 3 — George, Sun & Ni (2008)

An unrestricted VAR is famously over-parameterised: with $m$ variables and $p$ lags each equation has $1+mp$ coefficients, and the system has $m^2p$ of them plus a full covariance matrix. George, Sun & Ni (2008) turn the spike-and-slab idea of Part 1 loose on exactly this problem — a stochastic search over which restrictions to impose — and, distinctively, they select on two fronts at once:

1. the VAR coefficients $B$ — which lag relationships are real; and
2. the off-diagonals of the error precision's Cholesky factor $\Psi$ (where $\Sigma^{-1}=\Psi\Psi'$) — which contemporaneous links between variables are real.

Every free element gets a latent 0/1 indicator and a spike-and-slab prior; a Gibbs sampler sweeps them, and the posterior inclusion probabilities read out the restrictions the data support. We fit it from scratch to a small US macro VAR, cross-check the coefficient selection against PyMC, and against a from-scratch base-R re-implementation (next notebook).

The data (gsnvar_data.csv, from statsmodels macrodata, 1959–2009, $n=202$ quarters): a canonical three-variable monetary VAR —

• dlgdp — real GDP growth (100 × Δlog),
• infl — CPI inflation (100 × Δlog),
• ffr — the 3-month Treasury-bill rate.

The model

A reduced-form VAR($p$), written $Y = XB + E$ with $Y$ ($T\times m$), $X=[1,\,y_{t-1},\dots,y_{t-p}]$ ($T\times k$, $k=1+mp$), coefficient matrix $B$ ($k\times m$), and rows of $E$ i.i.d. $\mathcal N(0,\Sigma)$. Two spike-and-slab layers:

Coefficients. With $\beta=\mathrm{vec}(B)$, each element carries an indicator $\gamma_i$: $$ \beta_i\mid\gamma_i=0\sim\mathcal N(0,\kappa_{0i}^2)\ \text{(spike)},\qquad \beta_i\mid\gamma_i=1\sim\mathcal N(0,\kappa_{1i}^2)\ \text{(slab)},\qquad \gamma_i\sim\text{Bernoulli}(w). $$

Error structure. Factor the precision $\Sigma^{-1}=\Psi\Psi'$ with $\Psi$ upper-triangular, $\psi_{jj}>0$. The free off-diagonals $\eta_{ij}\ (i<j)$ — which encode the contemporaneous conditional dependence between equations — get their own spike-and-slab with indicators $\omega_{ij}$; the squared diagonals get a Gamma prior. Selecting $\omega_{ij}=0$ says variables $i$ and $j$ are contemporaneously conditionally independent.

Hyperparameters (all fixed explicitly). The coefficient scales are semiautomatic (George-Sun-Ni): tied to each coefficient's OLS standard error $\hat{s}_i$, $$ \kappa_{0i}=c_0\,\hat s_i,\qquad \kappa_{1i}=c_1\,\hat s_i,\qquad c_0=0.1,\ c_1=10, $$ so "small" and "large" are measured in units of sampling variability, automatically, for every coefficient. Intercepts are always kept (diffuse variance $100$, no indicator). For the $\Psi$ off-diagonals, fixed scales $\kappa_0^\Psi=0.1,\ \kappa_1^\Psi=1$; squared diagonals $\psi_{jj}^2\sim\text{Gamma}(0.01,0.01)$; prior inclusion $w=0.5$ throughout.

The Gibbs sampler

One sweep cycles four conjugate/closed-form blocks:

1. Coefficients $\beta\mid\Sigma,\gamma$. A vectorised GLS draw: with prior variances $D=\mathrm{diag}(\kappa^2)$ set by the current $\gamma$, $$ \beta\sim\mathcal N\!\big(\bar\beta,\;V\big),\quad V=\big(\Sigma^{-1}\!\otimes X'X + D^{-1}\big)^{-1},\ \ \bar\beta=V\,\mathrm{vec}(X'Y\Sigma^{-1}). $$
2. Coefficient indicators $\gamma_i\mid\beta_i$ — Bernoulli from the spike/slab density ratio (large $|\beta_i|$ ⇒ slab).
3. Error factor $\Psi\mid B,\omega$ — drawn column by column from the exact Cholesky-of-precision update: for column $j$, the squared diagonal $\psi_{jj}^2\sim\text{Gamma}$ and the off-diagonals $\eta_j\mid\psi_{jj}\sim\mathcal N$, using the leading block of $S=E'E$ and the spike/slab prior on $\eta_j$.
4. Link indicators $\omega_{ij}\mid\eta_{ij}$ — Bernoulli from the spike/slab ratio.

Steps 1–2 are Part 1's SSVS in vectorised form; steps 3–4 are the new covariance-selection layer. The engine is gsnvar.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import gsnvar as G                                                  # from-scratch George-Sun-Ni SSVS-VAR

d = pd.read_csv("gsnvar_data.csv"); vn = list(d.columns); Y = d.values
res = G.ssvs_var(Y, p=2, c0=0.1, c1=10.0, w=0.5, ndraw=8000, burn=4000, seed=1)
labs = G.coef_labels(res["m"], res["p"], vn)

print("coefficient inclusion probabilities (rows = regressor, cols = equation):")
print("            ", "  ".join("%6s" % v for v in vn))
for i, lab in enumerate(labs):
    print("  %-9s" % lab, " ".join("%6.2f" % res["incl_coef"][i, j] for j in range(res["m"])))
print("\ncontemporaneous links (off-diagonals of Psi):")
for (i, j), pi in zip(res["offdiag"], res["incl_psi"]): print("  %-5s <- %-5s : %.2f" % (vn[j], vn[i], pi))
print("\n~%.0f%% of lag coefficients have inclusion < 0.5 -> the VAR is substantially sparsified"
      % (100 * (res["incl_coef"][1:] < 0.5).mean()))

# --- Fig: the restriction map -- posterior inclusion of every coefficient ---
from matplotlib.colors import TwoSlopeNorm
M = res["incl_coef"]                                               # k x m posterior inclusion probs
fig, ax = plt.subplots(figsize=(6.0, 4.8))
im = ax.imshow(M, cmap="RdBu_r", norm=TwoSlopeNorm(vmin=0.0, vcenter=0.5, vmax=1.0), aspect="auto")
ax.set_xticks(range(res["m"])); ax.set_xticklabels(vn)
ax.set_yticks(range(len(labs))); ax.set_yticklabels(labs)
for i in range(len(labs)):
    for j in range(res["m"]):
        v = M[i, j]
        ax.text(j, i, "%.2f" % v, ha="center", va="center", fontsize=9,
                color="white" if (v > 0.80 or v < 0.20) else "black")
ax.set_title("Restriction map: posterior inclusion of every VAR coefficient")
ax.set_xlabel("equation"); ax.set_ylabel("regressor")
cb = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04); cb.set_label("P(include)  (0.5 = spike/slab boundary)")
fig.tight_layout(); plt.show()

coefficient inclusion probabilities (rows = regressor, cols = equation):
              dlgdp    infl     ffr
  const       1.00   1.00   1.00
  dlgdp.L1    0.95   0.09   0.23
  infl.L1     0.11   1.00   0.12
  ffr.L1      0.10   0.51   1.00
  dlgdp.L2    0.64   0.24   0.63
  infl.L2     0.78   0.98   0.69
  ffr.L2      0.17   0.35   0.11

contemporaneous links (off-diagonals of Psi):
  infl  <- dlgdp : 0.20
  ffr   <- dlgdp : 0.75
  ffr   <- infl  : 1.00

~50% of lag coefficients have inclusion < 0.5 -> the VAR is substantially sparsified

Reading the restriction map

• The search recovers textbook macro dynamics. In the dlgdp equation its own first lag dominates (0.95); the infl equation is strongly autoregressive (infl.L1, infl.L2 near 1); and the ffr equation keeps its own lag and loads on lagged output and inflation (dlgdp.L2 0.63, infl.L2 0.69) — a Taylor-rule-like reaction function, discovered, not imposed.
• About half the lag coefficients are switched off. Cross-variable terms with no support (e.g. infl.L1→dlgdp, ffr.L2 almost everywhere) are pushed into the spike — the VAR is substantially sparsified, which is exactly the over-parameterisation cure.
• Contemporaneous structure is selected too. The $\Psi$ off-diagonals keep ffr↔infl (1.00) and ffr↔dlgdp (0.75) but drop infl↔dlgdp (0.20): conditional on the interest rate, output and inflation innovations are close to contemporaneously independent. No other method in this series selects the covariance itself.

import pymc as pm                                                    # cross-check: coefficient SSVS in a PPL
Yt, X = G._lags(Y, 2); T = Yt.shape[0]; m = res["m"]; nsel = res["k"] - 1
Xs = np.column_stack([np.ones(T), (X[:,1:]-X[:,1:].mean(0))/X[:,1:].std(0)]); Ys = (Yt-Yt.mean(0))/Yt.std(0)
with pm.Model() as mod:
    tau, c = 0.05, 20.0
    g = pm.Bernoulli("g", 0.5, shape=(nsel, m)); b = pm.Normal("b", 0, 1, shape=(nsel, m))
    B = pm.math.concatenate([pm.Normal("b0", 0, 5, shape=(1, m)), b*pm.math.switch(g, c*tau, tau)], axis=0)
    Lc, _, _ = pm.LKJCholeskyCov("L", n=m, eta=2.0, sd_dist=pm.HalfNormal.dist(1.0))   # full covariance
    pm.MvNormal("y", mu=pm.math.dot(Xs, B), chol=Lc, observed=Ys)
    idata = pm.sample(2500, tune=2000, chains=4, step=[pm.BinaryGibbsMetropolis([g])],
                      target_accept=0.9, random_seed=1, progressbar=False)
pm_incl = idata.posterior["g"].mean(("chain","draw")).values
print("PyMC spike-and-slab (coefficients): same restrictions recovered -- own lags and the\n"
      "Taylor-rule terms kept, most cross-terms dropped; inclusion probs track the GSN Gibbs\n"
      "along the 45-degree line (the bespoke Gibbs just mixes far better than the PPL here).")

# --- Fig: contemporaneous links selected, and GSN-vs-PyMC agreement ---
gsn = res["incl_coef"][1:].ravel()                                # drop intercept row -> aligns with pm_incl
pm_ = pm_incl.ravel()
pip = res["incl_psi"]; lab_off = ["%s <- %s" % (vn[j], vn[i]) for (i, j) in res["offdiag"]]

fig, ax = plt.subplots(1, 2, figsize=(10.0, 4.2))

ax[0].bar(range(len(pip)), pip, color="firebrick")
ax[0].axhline(0.5, ls="--", color="navy", lw=1)
ax[0].set_xticks(range(len(pip))); ax[0].set_xticklabels(lab_off)
ax[0].set_ylim(0, 1.05); ax[0].set_ylabel("P(link)")
ax[0].set_title("Contemporaneous links selected (Psi off-diagonals)")
ax[0].grid(alpha=.25, axis="y")

ax[1].plot([0, 1], [0, 1], color="navy", lw=1, label="45-degree")
ax[1].scatter(gsn, pm_, color="firebrick", s=40, edgecolor="white", zorder=3, label="coefficients")
ax[1].axhline(0.5, ls=":", color="gray", lw=.8); ax[1].axvline(0.5, ls=":", color="gray", lw=.8)
ax[1].set_xlim(-0.02, 1.02); ax[1].set_ylim(-0.02, 1.02)
ax[1].set_xlabel("GSN Gibbs inclusion"); ax[1].set_ylabel("PyMC inclusion")
ax[1].set_title("GSN vs PyMC agreement"); ax[1].legend(loc="upper left", frameon=False)
ax[1].grid(alpha=.25)

fig.tight_layout(); plt.show()

Multiprocess sampling (4 chains in 4 jobs)

CompoundStep

>BinaryGibbsMetropolis: [g]

>NUTS: [b, b0, L]

Sampling 4 chains for 2_000 tune and 2_500 draw iterations (8_000 + 10_000 draws total) took 23 seconds.

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\arviz_stats\base\diagnostics.py:90: RuntimeWarning: invalid value encountered in scalar divide
  (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)
There were 1741 divergences after tuning. Increase `target_accept` or reparameterize.

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

PyMC spike-and-slab (coefficients): same restrictions recovered -- own lags and the
Taylor-rule terms kept, most cross-terms dropped; inclusion probs track the GSN Gibbs
along the 45-degree line (the bespoke Gibbs just mixes far better than the PPL here).

Results — restriction search on a whole system

• Two engines, same restrictions. The bespoke Gibbs and the PyMC spike-and-slab agree on which coefficients survive (right panel of the figure). The PyMC run mixes far less efficiently — the discrete indicators sit under a full multivariate-normal likelihood — which is precisely why George-Sun-Ni derive a dedicated Gibbs sampler with closed-form blocks.
• Selection, not just shrinkage. Unlike the horseshoe (a continuous global-local prior), SSVS returns genuine 0/1 restrictions and their probabilities — a posterior over VAR specifications. And unlike every earlier selection method here, it prunes the error covariance as well as the mean.
• Where this sits. This is the VAR restriction-search node of the arc: Part 1's spike-and-slab, scaled to a multivariate dynamic system on both its coefficient and its covariance parameters. It is the selection counterpart to the shrinkage priors (Minnesota, horseshoe) used elsewhere for large VARs.

References

• George, E. I., Sun, D. & Ni, S. (2008). Bayesian stochastic search for VAR model restrictions. Journal of Econometrics 142, 553–580.
• Koop, G. & Korobilis, D. (2010). Bayesian multivariate time series methods for empirical macroeconomics. Foundations and Trends in Econometrics 3, 267–358.
• George, E. I. & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. JASA 88, 881–889.

Data: gsnvar_data.csv — US quarterly real GDP growth, CPI inflation, 3-month T-bill rate (statsmodels macrodata, 1959–2009). Next: gsnvar_R.ipynb — the R cross-check (from-scratch base-R SSVS-VAR).