Bayesian Vector Autoregression

Model

A vector autoregression is the workhorse of empirical macroeconomics: each variable is regressed on $p$ lags of every variable, letting the data speak about dynamic interdependence. This example fits the classic small US monetary VAR — GDP growth, CPI inflation, and the 3-month T-bill (1959–2009) — three ways: the classical frequentist VAR for the identification machinery, a Minnesota-prior Bayesian VAR, and the Villani steady-state BVAR, all on a from-scratch Gibbs engine (bvar.py).

y_t = c + \Pi_1 y_{t-1} + \dots + \Pi_p y_{t-p} + \varepsilon_t, \qquad \varepsilon_t \sim N(0, \Sigma)

$y_t$ — the $m$ -vector of variables (here GDP growth, inflation, T-bill)
$\Pi_l$ — the lag- $l$ coefficient matrix; stationarity ⇔ companion-matrix eigenvalues inside the unit circle
$\Sigma$ — innovation covariance; a Cholesky factor gives the structural shocks

Classical VAR & structural identification

The classical VAR. With the same regressors in every equation, equation-by-equation OLS is efficient. The reduced-form innovations are correlated, so structural shocks are recovered by a Cholesky (recursive) identification — the ordering [GDP, Inflation, T-bill] puts the monetary shock last (the policy rate reacts within-quarter but affects output/inflation only with a lag). This yields Granger-causality tests, orthogonalized impulse responses, and the forecast-error variance decomposition. The IRF even reproduces the famous price puzzle (inflation rising after a tightening) — a known artefact of recursive identification, and a reminder that an IRF is only as credible as its identifying assumption.

Minnesota-prior BVAR

The Minnesota prior (Litterman 1986). An unrestricted VAR overfits (21 coefficients here at $p=2$ ), so the Minnesota prior shrinks each equation toward a parsimonious benchmark, governed by three hyperparameters: $\lambda_1$ (overall tightness — the master knob, $\lambda_1\to 0$ forces the benchmark, $\to\infty$ recovers OLS), $\lambda_2$ (cross-variable tightness — shrink other variables' lags harder, set per-pair to $0.8\,\widehat{\text{corr}}$ ), and $\lambda_3$ (lag decay). An independent Normal–Inverse-Wishart version gives a conjugate 2-block Gibbs sampler (the coefficient block is the SUR/GLS update tying the equations through $\Sigma^{-1}$ ).

Villani steady-state BVAR

The Villani steady-state prior (2009) is the macro headline. A reduced-form VAR's long-run mean is the tangled $(I-\sum_l\Pi_l)^{-1}c$ — you cannot state a belief about it. Villani re-parameterises around the unconditional mean $\psi$ directly, $y_t = \psi + \sum_l \Pi_l(y_{t-l}-\psi) + \varepsilon_t$ , so an informative prior on the steady state can pin the long run — e.g. a 2% inflation target and trend growth — and forecasts revert to that, not to whatever the (crisis-laden) sample average happens to be. A 3-block Gibbs sampler cycles the dynamics, the steady state, and the covariance.

Notebooks

Part 1 — the monetary VAR. The structural story is textbook (Taylor-rule-like policy reaction in the FEVD; output falling ~3 quarters after a tightening). The steady-state prior earns its keep: where the Minnesota BVAR reverts to the sample-implied mean (baking in the high inflation of the 1970s–80s, ~3.5–4%), the Villani BVAR pulls long-run inflation down to a pinned ~2.6%, expressing a modern inflation-targeting view the data alone would not give — and out of sample through the 2008 collapse it has the lowest RMSE (3.90), beating the OLS-VAR (4.20), the Minnesota BVAR (4.15), and even a random walk (4.12). Part 2 — the Treasury yield curve. The same three models on four constant-maturity yields (3m, 5y, 10y, 30y; monthly 1994–2024) — now near-unit-root levels where the natural steady state is a neutral curve. Policy pass-through decays up the curve: a short-rate shock explains 45% of the 3-month's variance but only 11% of the 30-year's, and the IRF shows the front end tracking policy while the long end stays anchored by expectations and the term premium — the textbook curve-flattening picture. The instructive twist: here the pinned neutral curve almost coincides with the sample means, so the steady-state and Minnesota forecasts nearly overlap. That is precisely the lesson — the steady-state prior earns its keep only when the sample average is unrepresentative of the future (as in Part 1's high-inflation sample); on this yield sample, by luck, it isn't. Yields are near-random-walk, so out of sample the random walk wins, but among the models the steady-state BVAR is best. Both applications cross-checked in R (vars, BVAR).

Python Part 1 · Monetary VAR — classical + Minnesota & Villani steady-state BVAR from scratch: IRFs, FEVD, forecasting (Python) Open → R Part 1 · Monetary VAR — R cross-check: vars (classical, Granger, IRFs) + BVAR (Minnesota) (R) Open → Python Part 2 · Yield-curve VAR — the same three models on the Treasury curve, policy pass-through across maturities (Python) Open → R Part 2 · Yield-curve VAR — R cross-check (vars + BVAR) (R) Open →

Downloads

bvar.py From-scratch Minnesota + Villani steady-state BVAR Gibbs samplers (NumPy / scipy) bvar_macro_data.csv US quarterly monetary VAR, 1959Q2–2009Q3 — GDP growth, CPI inflation, 3-month T-bill bvar_yields_data.csv US Treasury yields, monthly 1994–2024 — 3-month, 5-year, 10-year, 30-year

References

Sims, C. A. (1980). Macroeconomics and reality. Econometrica 48(1), 1–48. — the vector autoregression and structural identification
Doan, T., Litterman, R. & Sims, C. (1984). Forecasting and conditional projection using realistic prior distributions. Econometric Reviews 3(1), 1–100. — the Minnesota prior
Litterman, R. B. (1986). Forecasting with Bayesian vector autoregressions — five years of experience. Journal of Business & Economic Statistics 4(1), 25–38. — Minnesota-prior BVAR forecasting
Kadiyala, K. R. & Karlsson, S. (1997). Numerical methods for estimation and inference in Bayesian VAR models. Journal of Applied Econometrics 12(2), 99–132. — the Normal–Inverse-Wishart / Gibbs machinery
Villani, M. (2009). Steady-state priors for vector autoregressions. Journal of Applied Econometrics 24(4), 630–650. — the steady-state (unconditional-mean) prior
Koop, G. & Korobilis, D. (2010). Bayesian multivariate time series methods for empirical macroeconomics. Foundations and Trends in Econometrics 3(4), 267–358. — a survey of BVAR methods

BVAR Module — Source Code

"""From-scratch Bayesian VAR with Minnesota prior and a Gibbs sampler.

Follows Meseguer (2010): mortality is modelled in **first differences of log
rates** (Delta log m), separately by sex, over the 22 abridged age groups. Two
models (both estimated by a hand-written Gibbs sampler, general lag order p):

  M1  "Minnesota"    : reduced-form VAR with intercept, Minnesota prior on the AR
                       coefficients (own first lag centred on 0 for differenced
                       data), independent Normal-inverse-Wishart -> 2-block Gibbs.
  M2  "steady state" : Villani mean-adjusted VAR y_t = psi + sum Pi_l (y_{t-l}-psi)
                       + e, with the same Minnesota prior on Pi and an informative
                       Normal prior on the steady state psi -> 3-block Gibbs.

Minnesota prior standard deviations (Meseguer eq. 16), for coefficient on variable
j at lag l in equation i:
    j = i : lam1 / l**lam3
    j != i: lam1 * lam2(i,j) / l**lam3 * (s_i / s_j)
with s_i the residual std of an AR(p) fit of series i, and lam2(i,j) set to
0.8 * corr(Delta log y_i, Delta log y_j).
"""

import numpy as np
from scipy.stats import invwishart


# --------------------------------------------------------------------------- #
# data / prior construction
# --------------------------------------------------------------------------- #

def build_regressors(Y, p, intercept):
    """Stack a VAR(p) design. Y is (T, m). Returns (Yt, X) with X lag-blocks
    ordered [lag1 (m cols), lag2, ...], optionally prefixed by an intercept."""
    T, m = Y.shape
    rows = T - p
    lags = [Y[p - l:T - l] for l in range(1, p + 1)]     # each (rows, m)
    X = np.hstack(lags)
    if intercept:
        X = np.hstack([np.ones((rows, 1)), X])
    return Y[p:], X


def _series_scale(Y, p):
    """s_i = residual std of an AR(p) regression of each series on its own lags."""
    T, m = Y.shape
    s = np.empty(m)
    for i in range(m):
        yi = Y[p:, i]
        Xi = np.column_stack([np.ones(T - p)] + [Y[p - l:T - l, i] for l in range(1, p + 1)])
        beta, *_ = np.linalg.lstsq(Xi, yi, rcond=None)
        s[i] = (yi - Xi @ beta).std(ddof=len(beta))
    return s


def minnesota_moments(Y, p, lam1=0.2, lam2=0.5, lam3=1.0, own_mean=0.0,
                      intercept=True, intercept_sd=1e3):
    """Prior mean B0 and prior variance V0, each shaped (k, m) with k regressors
    per equation and columns = equations. `lam2` may be a scalar or an (m, m)
    matrix lam2(i,j) (e.g. 0.8 * corr)."""
    T, m = Y.shape
    s = _series_scale(Y, p)
    k = (1 if intercept else 0) + m * p
    B0 = np.zeros((k, m))
    V0 = np.zeros((k, m))
    off = 1 if intercept else 0
    if intercept:
        V0[0, :] = intercept_sd ** 2                     # diffuse intercept
    lam2_mat = lam2 if np.ndim(lam2) == 2 else np.full((m, m), lam2)
    for i in range(m):                                   # equation i
        for l in range(1, p + 1):
            for j in range(m):                           # regressor variable j
                r = off + (l - 1) * m + j
                if j == i and l == 1:
                    B0[r, i] = own_mean
                if j == i:
                    V0[r, i] = (lam1 / l ** lam3) ** 2
                else:
                    V0[r, i] = (lam1 * lam2_mat[i, j] / l ** lam3 * s[i] / s[j]) ** 2
    return B0, V0, s


def corr_lam2(Y, p, weight=0.8):
    """lam2(i,j) = weight * corr(Delta log y_i, Delta log y_j), used for cross shrinkage."""
    C = np.corrcoef(Y[p:].T)
    return weight * C


# --------------------------------------------------------------------------- #
# Gaussian draw from a given precision matrix
# --------------------------------------------------------------------------- #

def _draw_from_precision(rhs, Prec, rng):
    """Sample x ~ N(Prec^{-1} rhs, Prec^{-1}) via Cholesky."""
    L = np.linalg.cholesky(Prec)                         # Prec = L L^T
    mean = np.linalg.solve(L.T, np.linalg.solve(L, rhs))
    z = rng.standard_normal(rhs.shape[0])
    return mean + np.linalg.solve(L.T, z)


# --------------------------------------------------------------------------- #
# Model 1 - Minnesota reduced-form BVAR (2-block Gibbs)
# --------------------------------------------------------------------------- #

def gibbs_minnesota(Y, p=1, lam1=0.2, lam2=None, lam3=1.0, own_mean=0.0,
                    ndraw=3000, burn=1000, thin=1, seed=0, iw_df=None, iw_scale=None):
    """2-block Gibbs for y_t = c + sum_l Pi_l y_{t-l} + e, e~N(0,Sigma)."""
    rng = np.random.default_rng(seed)
    Yt, X = build_regressors(Y, p, intercept=True)
    rows, m = Yt.shape
    k = X.shape[1]
    if lam2 is None:
        lam2 = corr_lam2(Y, p)
    B0, V0, s = minnesota_moments(Y, p, lam1, lam2, lam3, own_mean, intercept=True)
    b0 = B0.reshape(-1, order="F")                       # vec by equation (column-major)
    V0inv = np.diag(1.0 / V0.reshape(-1, order="F"))
    XtX, XtY = X.T @ X, X.T @ Yt

    v0 = iw_df if iw_df else m + 2
    S0 = iw_scale if iw_scale is not None else np.diag(s ** 2)
    B = B0.copy()
    Sig = np.diag(s ** 2)

    kept_B, kept_S = [], []
    for it in range(ndraw + burn):
        Sinv = np.linalg.inv(Sig)
        Prec = V0inv + np.kron(Sinv, XtX)
        rhs = V0inv @ b0 + (XtY @ Sinv).reshape(-1, order="F")
        beta = _draw_from_precision(rhs, Prec, rng)
        B = beta.reshape(k, m, order="F")
        E = Yt - X @ B
        Sig = invwishart.rvs(df=v0 + rows, scale=S0 + E.T @ E, random_state=rng)
        if it >= burn and (it - burn) % thin == 0:
            kept_B.append(B.copy()); kept_S.append(Sig.copy())
    return {"B": np.array(kept_B), "Sigma": np.array(kept_S), "p": p,
            "intercept": True, "kind": "minnesota"}


# --------------------------------------------------------------------------- #
# Model 2 - Villani steady-state (mean-adjusted) BVAR (3-block Gibbs)
# --------------------------------------------------------------------------- #

def gibbs_steady_state(Y, p=1, lam1=0.2, lam2=None, lam3=1.0, own_mean=0.0,
                       psi_prior_mean=None, psi_prior_sd=10.0,
                       ndraw=3000, burn=1000, thin=1, seed=0, iw_df=None, iw_scale=None):
    """3-block Gibbs for y_t = psi + sum_l Pi_l (y_{t-l}-psi) + e.

    psi_prior_mean : (m,) informative steady-state centre (e.g. SSA Alt-2 ultimate
        Delta log m). Default 0 (diffuse). psi_prior_sd : scalar or (m,) prior sd.
    """
    rng = np.random.default_rng(seed)
    T, m = Y.shape
    rows = T - p
    if lam2 is None:
        lam2 = corr_lam2(Y, p)
    B0, V0, s = minnesota_moments(Y, p, lam1, lam2, lam3, own_mean, intercept=False)
    b0 = B0.reshape(-1, order="F")
    V0inv = np.diag(1.0 / V0.reshape(-1, order="F"))
    k = m * p

    psi0 = np.zeros(m) if psi_prior_mean is None else np.asarray(psi_prior_mean, float)
    psd = np.full(m, psi_prior_sd) if np.ndim(psi_prior_sd) == 0 else np.asarray(psi_prior_sd)
    Lam0inv = np.diag(1.0 / psd ** 2)

    v0 = iw_df if iw_df else m + 2
    S0 = iw_scale if iw_scale is not None else np.diag(s ** 2)

    Ylag_full = np.hstack([Y[p - l:T - l] for l in range(1, p + 1)])    # (rows, m*p) original units
    Yt = Y[p:]
    mu = Y.mean(axis=0).copy() if psi_prior_mean is None else psi0.copy()
    Pi = B0.copy()                                        # (k, m)
    Sig = np.diag(s ** 2)

    kept_Pi, kept_mu, kept_S = [], [], []
    for it in range(ndraw + burn):
        Sinv = np.linalg.inv(Sig)
        # --- block 1: Pi | mu, Sigma  (regress demeaned y on demeaned lags) ---
        Z = Y - mu                                       # (T, m) demeaned
        Zt, Zlag = build_regressors(Z, p, intercept=False)
        ZtZ, ZtY = Zlag.T @ Zlag, Zlag.T @ Zt
        Prec = V0inv + np.kron(Sinv, ZtZ)
        rhs = V0inv @ b0 + (ZtY @ Sinv).reshape(-1, order="F")
        Pi = _draw_from_precision(rhs, Prec, rng).reshape(k, m, order="F")

        # --- block 2: mu | Pi, Sigma  (u_t = (I - sum Pi_l) mu + e) ---
        U = Yt - Ylag_full @ Pi                          # (rows, m)
        sumPi = sum(Pi[l * m:(l + 1) * m, :] for l in range(p))   # sum of blocks
        D = np.eye(m) - sumPi.T
        Prec_mu = Lam0inv + rows * (D.T @ Sinv @ D)
        rhs_mu = Lam0inv @ psi0 + D.T @ Sinv @ U.sum(axis=0)
        mu = _draw_from_precision(rhs_mu, Prec_mu, rng)

        # --- block 3: Sigma | Pi, mu ---
        Z = Y - mu
        Zt, Zlag = build_regressors(Z, p, intercept=False)
        E = Zt - Zlag @ Pi
        Sig = invwishart.rvs(df=v0 + rows, scale=S0 + E.T @ E, random_state=rng)

        if it >= burn and (it - burn) % thin == 0:
            kept_Pi.append(Pi.copy()); kept_mu.append(mu.copy()); kept_S.append(Sig.copy())
    return {"Pi": np.array(kept_Pi), "mu": np.array(kept_mu), "Sigma": np.array(kept_S),
            "p": p, "kind": "steady_state"}


# --------------------------------------------------------------------------- #
# forecasting -> draws of future Delta log m
# --------------------------------------------------------------------------- #

def stationary_mask(fit):
    """Boolean mask of posterior draws whose VAR is covariance-stationary
    (all companion-matrix eigenvalues inside the unit circle)."""
    p = fit["p"]
    coefs = fit["B"] if fit["kind"] == "minnesota" else fit["Pi"]
    off = 1 if fit["kind"] == "minnesota" else 0
    m = fit["Sigma"].shape[1]
    mask = np.ones(len(coefs), bool)
    eye = np.eye(m * (p - 1)) if p > 1 else None
    for d, B in enumerate(coefs):
        Pis = [B[off + (l - 1) * m: off + l * m].T for l in range(1, p + 1)]  # Pi_l = block^T
        comp = np.zeros((m * p, m * p))
        comp[:m] = np.hstack(Pis)
        if p > 1:
            comp[m:, :m * (p - 1)] = eye
        mask[d] = np.max(np.abs(np.linalg.eigvals(comp))) < 1.0
    return mask


def filter_stationary(fit):
    """Return a copy of `fit` keeping only covariance-stationary draws."""
    mask = stationary_mask(fit)
    out = dict(fit)
    for key in ("B", "Pi", "mu", "Sigma"):
        if key in fit:
            out[key] = fit[key][mask]
    out["stationary_frac"] = float(mask.mean())
    return out


def forecast(fit, Y, H, seed=1):
    """Simulate H-step-ahead paths of the modelled series (Delta log m) for every
    posterior draw. Returns array (ndraw, H, m)."""
    rng = np.random.default_rng(seed)
    p = fit["p"]
    m = Y.shape[1]
    hist0 = Y[-p:]                                        # last p observations
    if fit["kind"] == "minnesota":
        Bs, Ss = fit["B"], fit["Sigma"]
    else:
        Bs, Ss, mus = fit["Pi"], fit["Sigma"], fit["mu"]
    nd = len(Ss)
    out = np.empty((nd, H, m))
    for d in range(nd):
        Sig = Ss[d]
        chol = np.linalg.cholesky(Sig)
        hist = [hist0[i].copy() for i in range(p)]       # hist[-1] = most recent
        for h in range(H):
            eps = chol @ rng.standard_normal(m)
            if fit["kind"] == "minnesota":
                B = Bs[d]
                x = B[0].copy()                          # intercept
                for l in range(1, p + 1):
                    x += B[1 + (l - 1) * m:1 + l * m].T @ hist[-l]
                y = x + eps
            else:
                Pi, mu = Bs[d], mus[d]
                x = mu.copy()
                for l in range(1, p + 1):
                    x += Pi[(l - 1) * m:l * m].T @ (hist[-l] - mu)
                y = x + eps
            out[d, h] = y
            hist.append(y)
        # keep only last p in hist via slicing not needed (indexing by -l)
    return out