Structural VAR — Shock Identification

Python · R · Part 2 of the VAR arc

Overview

Part 2 of the VAR arc: structural identification. The reduced-form BVAR gives correlated innovations εt\varepsilon_t that are not economically interpretable — a structural shock is εt=Put\varepsilon_t = P u_t for some impact matrix PP with PP=ΣPP' = \Sigma. Part 1 chose PP the easy way — a recursive (Cholesky) ordering — and paid the price: the macro impulse responses showed the notorious price puzzle (inflation rising after a monetary tightening). This example fixes it with credible identification, adding three modern schemes on top of the same bvar.py posterior (the engine is left untouched; each scheme is a separate module).

εt=Put,PP=Σ,P=chol(Σ)Q  (Q orthogonal, set-identified)\varepsilon_t = P\,u_t, \quad PP' = \Sigma, \qquad P = \text{chol}(\Sigma)\,Q \ \ (Q\ \text{orthogonal, set-identified})

Sign restrictions

Sign restrictions (Uhlig 2005; Rubio-Ramírez, Waggoner & Zha 2010). The data do not pin the rotation QQ in P=chol(Σ)QP = \text{chol}(\Sigma)\,Q, so the model is set-identified. Instead of a questionable ordering, impose only the signs economics is confident about (a monetary tightening has the rate up, inflation down) and keep every rotation consistent with them — the identified set, reported as an IRF band. On the macro system this dissolves the price puzzle while honestly reporting that the output effect is uncertain (the agnostic band spans zero). The same machine decomposes the Treasury curve into economically-labelled level and slope shocks — sign-identified structural shocks with dynamics and uncertainty, not static PCA loadings.

External instruments & long-run restrictions

External instruments & long-run restrictions. A proxy / external-instrument SVAR (Stock–Watson 2012/18; Mertens–Ravn 2013) brings in an outside variable correlated with the target shock but nothing else — here the Romer–Romer monetary shock pins the impact vector. Long-run restrictions (Blanchard–Quah 1989) impose what a shock does in the limit: the impact matrix is chosen so a demand shock has no permanent effect on output, i.e. ΞP=chol(ΞΣΞ)\Xi P = \text{chol}(\Xi\Sigma\Xi') with the long-run multiplier Ξ=(IlΠl)1\Xi = (I-\sum_l\Pi_l)^{-1}. The supply shock then lifts output to a permanent ~2.5% higher level while the demand shock's output effect dies out — the clean Blanchard–Quah decomposition.

Four identifications, one posterior

A menu, not a winner. Across Part 2 the same reduced form is identified four ways — recursive, sign, external-instrument, long-run — each imposing a different kind of belief: a timing zero, an IRF sign, an outside instrument, a limit. None is universally right; the craft is matching the restriction you can defend to the question you are asking (the external instrument is the current gold standard for a single well-measured shock; sign and long-run restrictions shine when a credible instrument is unavailable). Crucially, every band rides the same Bayesian posterior bvar.py already produces — identification is a layer on top, not a re-estimation. Both schemes are cross-checked in R (vars::BQ, a base-R rotation search, and a by-hand proxy regression).

Notebooks

Downloads

References

Sign-Restriction Module — Source Code

"""
Sign-restricted structural VAR (Uhlig 2005; Rubio-Ramirez, Waggoner & Zha 2010), on top of `bvar.py`.

The reduced-form BVAR posterior (bvar.gibbs_minnesota / gibbs_steady_state) gives draws of the VAR
coefficients and Sigma. Any structural decomposition is eps_t = P u_t with P = chol(Sigma) Q for some
orthogonal Q; the data do not pin Q, so it is SET-identified. Instead of a recursive ordering we impose
only the *signs* we believe in (e.g. a monetary tightening has rate up, inflation down) and keep every
rotation consistent with them -- the identified set, reported as an IRF band.

This module does NOT modify bvar.py: it consumes a bvar fit dict and adds the rotation search.
"""
import numpy as np


def _pi_blocks(coefs_draw, kind, m, p):
    """Extract the list [Pi_1, ..., Pi_p] (each m x m) from one posterior coefficient draw."""
    off = 1 if kind == "minnesota" else 0
    return [coefs_draw[off + (l - 1) * m: off + l * m].T for l in range(1, p + 1)]


def _ma_irf(Pis, m, p, H):
    """Reduced-form moving-average IRFs Theta_h (Theta_0 = I), shape (H+1, m, m), via the companion matrix."""
    F = np.zeros((m * p, m * p)); F[:m] = np.hstack(Pis)
    if p > 1:
        F[m:, :m * (p - 1)] = np.eye(m * (p - 1))
    J = np.zeros((m, m * p)); J[:, :m] = np.eye(m)
    Theta = np.empty((H + 1, m, m)); Fh = np.eye(m * p)
    for h in range(H + 1):
        Theta[h] = J @ Fh @ J.T; Fh = Fh @ F
    return Theta


def _satisfies(g, restr):
    """g is a candidate structural IRF (H+1, m); restr = list of (var_idx, sign, hmax)."""
    for vi, sgn, hmax in restr:
        if not (sgn * g[:hmax + 1, vi] >= 0).all():
            return False
    return True


def sign_irf(fit, restr, H=20, max_tries=200, seed=0):
    """One structural shock identified by sign restrictions.

    fit   : a bvar fit dict (kind 'minnesota' or 'steady_state').
    restr : list of (variable_index, +1/-1, max_horizon) the shock's IRF must obey.
    Returns dict(irf=(n_accepted, H+1, m), accept_frac). For each posterior draw we draw random
    Haar rotations until one column's IRF satisfies the signs (trying both sign flips), then keep it.
    """
    rng = np.random.default_rng(seed)
    coefs = fit["B"] if fit["kind"] == "minnesota" else fit["Pi"]
    Sig = fit["Sigma"]; m = Sig.shape[1]; p = fit["p"]
    kept = []; tries_hit = 0
    for d in range(len(coefs)):
        Theta = _ma_irf(_pi_blocks(coefs[d], fit["kind"], m, p), m, p, H)
        P0 = np.linalg.cholesky(Sig[d])
        found = None
        for _ in range(max_tries):
            q = rng.standard_normal(m); q /= np.linalg.norm(q)     # Haar-uniform impulse direction
            g = Theta @ (P0 @ q)                                   # candidate structural IRF (H+1, m)
            if _satisfies(g, restr):
                found = g; break
            if _satisfies(-g, restr):
                found = -g; break
        if found is not None:
            kept.append(found)
        else:
            tries_hit += 1
    return {"irf": np.array(kept), "accept_frac": len(kept) / len(coefs)}


def bands(res, qs=(16, 50, 84)):
    """Pointwise percentile bands of the accepted IRFs, shape (len(qs), H+1, m)."""
    return np.percentile(res["irf"], qs, axis=0)