# Vector Autoregression (VAR)
Vector Autoregression (VAR) with `numpyro_forecast`
This notebook ports the blog post [**Bayesian VAR in NumPyro**](https://juanitorduz.github.io/var_numpyro/) to the [`numpyro_forecast`](https://github.com/juanitorduz/numpyro_forecast) package. A vector autoregression (VAR) models several time series jointly: each series is regressed on the past values of all series, and the shocks are correlated across series. We fit a VAR with two lags to the quarterly growth rates of US real GDP, consumption and investment, sample the posterior with NUTS, forecast 30 quarters ahead, and compute impulse response functions (IRFs), the standard tool to read a VAR.
The package provides the VAR pieces as reusable components. You do not write the lag recursion, the forecast loop or the IRF recursion yourself:
- [`var_step`](https://juanitorduz.github.io/numpyro_forecast/reference/var.var_step.html) turns sampled coefficients into a step for the [`ssoe`](https://juanitorduz.github.io/numpyro_forecast/reference/models.ssoe.html) building block. The block runs the in-sample recursion and the generative forecast.
- [`impulse_response`](https://juanitorduz.github.io/numpyro_forecast/reference/var.impulse_response.html) computes the responses for all posterior draws at once, with optional orthogonalization and cumulation.
- [`companion_matrix`](https://juanitorduz.github.io/numpyro_forecast/reference/var.companion_matrix.html) gives the stability check.
- [`minnesota_prior`](https://juanitorduz.github.io/numpyro_forecast/reference/priors.minnesota_prior.html) returns the moments of the Minnesota shrinkage prior. It lives in a separate module and is independent of the VAR code: the prior is always your own `numpyro.sample` call.
The components are deliberately minimal. If you need a complete Bayesian VAR toolkit (identification schemes, variance decompositions, lag selection), see [Impulso](https://github.com/thomaspinder/impulso) by Thomas Pinder. Its `MinnesotaPrior` parameterization and its batched moving-average recursion inspired the two helpers used here.
# Prepare notebook
``` python
import datetime as dt
import itertools
import warnings
import arviz as az
import jax.numpy as jnp
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import numpyro
import numpyro.distributions as dist
import polars as pl
import xarray as xr
from jax import random
from jaxtyping import Float
from numpyro.infer import MCMC, NUTS
from numpyro_forecast import Horizon, eval_crps, predictions_to_datatree, ssoe, to_datatree
from numpyro_forecast.arrays import pad_future
from numpyro_forecast.priors import minnesota_prior
from numpyro_forecast.typing import Array, ForecastModel
from numpyro_forecast.var import companion_matrix, impulse_response, var_step
az.style.use("arviz-darkgrid")
plt.rcParams["figure.figsize"] = [10, 6]
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.facecolor"] = "white"
warnings.filterwarnings(
"ignore", message="When multiple credible intervals are plotted", category=UserWarning
)
numpyro.set_host_device_count(n=4)
rng_key = random.PRNGKey(seed=42)
%load_ext autoreload
%autoreload 2
%load_ext jaxtyping
%jaxtyping.typechecker beartype.beartype
```
# Read data
We use the `macrodata` dataset shipped with `statsmodels`: quarterly US macroeconomic data from 1959Q1 to 2009Q3, compiled from the Federal Reserve Bank of St.Β Louis (FRED) and released in the public domain. We read the CSV from the `statsmodels` repository and keep three series in billions of chained 2005 dollars: real GDP, real personal consumption and real gross private domestic investment.
The levels trend upward and are not stationary. We take log differences, which turn levels into quarter-on-quarter growth rates, and multiply by 100 to read them in percent. A VAR assumes stationarity, and growth rates are a standard way to get there for macro aggregates.
``` python
def quarter_label(d: dt.date) -> str:
"""Format a date as e.g. ``1959Q2``."""
return f"{d.year}Q{(d.month - 1) // 3 + 1}"
url = (
"https://raw.githubusercontent.com/statsmodels/statsmodels/main/"
"statsmodels/datasets/macrodata/macrodata.csv"
)
macro_df = pl.read_csv(url)
names = ["realgdp", "realcons", "realinv"]
y_pct = macro_df.select(
pl.date(pl.col("year"), (pl.col("quarter") - 1) * 3 + 1, 1).alias("date"),
*[(pl.col(name).log().diff() * 100).alias(name) for name in names],
).drop_nulls()
print(
f"shape: {y_pct.shape}, "
f"from {quarter_label(y_pct['date'][0])} to {quarter_label(y_pct['date'][-1])}"
)
y_pct.head()
```
shape: (202, 4), from 1959Q2 to 2009Q3
shape: (5, 4)
| date | realgdp | realcons | realinv |
|------------|-----------|----------|------------|
| date | f64 | f64 | f64 |
| 1959-04-01 | 2.494213 | 1.528611 | 8.021268 |
| 1959-07-01 | -0.119295 | 1.038598 | -7.213104 |
| 1959-10-01 | 0.349453 | 0.108401 | 3.442511 |
| 1960-01-01 | 2.219018 | 0.953415 | 10.266377 |
| 1960-04-01 | -0.468455 | 1.257243 | -10.669385 |
``` python
stats = ["mean", "std", "min", "max"]
y_pct.select(names).describe().filter(pl.col("statistic").is_in(stats)).with_columns(
pl.col(names).round(3)
)
```
shape: (4, 4)
| statistic | realgdp | realcons | realinv |
|-----------|---------|----------|---------|
| str | f64 | f64 | f64 |
| "mean" | 0.776 | 0.837 | 0.814 |
| "std" | 0.88 | 0.694 | 4.685 |
| "min" | -2.071 | -2.296 | -19.316 |
| "max" | 3.859 | 2.773 | 12.209 |
Investment growth is about five times more volatile than GDP or consumption growth. Keep this in mind for the Minnesota prior section: the three series are not on a common scale.
``` python
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(12, 8), layout="constrained")
for ax, name, color in zip(axes, names, ("C0", "C1", "C2"), strict=True):
ax.plot(y_pct["date"], y_pct[name], color=color, lw=1.2, label=name)
ax.axhline(0.0, color="gray", lw=0.8, ls="--")
ax.set(ylabel="percent")
ax.legend(loc="upper right")
axes[-1].set(xlabel="date")
fig.suptitle("Quarterly growth rates (100 x log difference)", fontsize=16, fontweight="bold");
```
## Stability
A VAR is stable when all eigenvalues of its companion matrix have modulus below one. Stability is what makes the forecast revert to a finite unconditional mean and the impulse responses die out. Because the impulse responses are a nonlinear function of the coefficients, a single explosive posterior draw would dominate their posterior mean at long horizons, so we check the share of stable draws before the IRF section and mask the unstable draws if there are any.
``` python
phi_draws = jnp.asarray(posterior["phi"]) # (4000, 2, 3, 3)
radius = np.abs(np.linalg.eigvals(np.asarray(companion_matrix(phi_draws)))).max(axis=-1)
stable = radius < 1.0
print(
f"stable draws: {stable.mean():.3f}, median spectral radius: {np.median(radius):.3f}, "
f"max: {radius.max():.3f}"
)
```
stable draws: 1.000, median spectral radius: 0.597, max: 0.827
# In-sample fit
The posterior predictive of the `obs` site gives the one-step-ahead predictive distribution for every in-sample quarter. We plot the 50\\ and 94\\ HDI bands per series and score the fit with the continuous ranked probability score (CRPS).
``` python
def hdi_label(prob: float, prefix: str = "") -> str:
r"""Legend label for an HDI band, e.g. ``$94\%$ HDI``."""
percent = f"{prob:.0%}".replace("%", r"\%")
return f"{prefix}${percent}$ HDI"
def stack_draws(group: str, tree_: xr.DataTree) -> Float[np.ndarray, " sample time series"]:
"""Flatten ``(chain, draw)`` of the ``obs`` variable of a tree group into a sample axis."""
da = tree_[group].dataset["obs"]
return da.stack(sample=("chain", "draw")).transpose("sample", "time", "obs_dim").to_numpy()
hdi_probs = (0.5, 0.94)
hdi_alphas = [0.6, 0.3] # 50% band darker, 94% band lighter
dates_num = mdates.date2num(dates)
future_dates_num = mdates.date2num(future_dates)
def plot_fit_and_forecast(
train_draws: Float[np.ndarray, " sample time series"],
future_draws: Float[np.ndarray, " sample future series"] | None,
title: str,
n_last: int | None = None,
) -> None:
"""Facet the in-sample predictive (and optionally the forecast) per series.
Parameters
----------
train_draws
In-sample posterior predictive draws ``(sample, time, series)``.
future_draws
Forecast draws ``(sample, future, series)``, or ``None`` for the in-sample plot only.
title
Figure title.
n_last
Plot only the last ``n_last`` in-sample quarters (``None`` for all).
"""
start = 0 if n_last is None else train_draws.shape[1] - n_last
x_train = dates_num[start:]
observed = np.asarray(data)[start:]
pc = az.plot_lm(
predictions_to_datatree(train_draws[:, start:], x_train, names, observed=observed),
y="obs",
x="t",
plot_dim="time",
ci_kind="hdi",
ci_prob=hdi_probs,
smooth=False,
col_wrap=1,
visuals={
"ci_band": {"color": "C0"},
"observed_scatter": False,
"pe_line": False,
"xlabel": False,
"ylabel": False,
},
aes={"alpha": ["prob"]},
alpha=hdi_alphas,
figure_kwargs={"figsize": (12, 9)},
)
train_bands = pc.viz["ci_band"]["t"].sel(series=names[0])
handles = [train_bands.sel(prob=prob).item() for prob in (0.94, 0.5)]
for handle, prob in zip(handles, (0.94, 0.5), strict=True):
handle.set_label(hdi_label(prob, prefix="in-sample " if future_draws is not None else ""))
if future_draws is not None:
az.plot_lm(
predictions_to_datatree(future_draws, future_dates_num, names),
y="obs",
x="t",
plot_dim="time",
plot_collection=pc,
ci_kind="hdi",
ci_prob=hdi_probs,
smooth=False,
visuals={
"ci_band": {"color": "C1"},
"observed_scatter": False,
"pe_line": False,
"xlabel": False,
"ylabel": False,
},
)
future_bands = pc.viz["ci_band"]["t"].sel(series=names[0])
for prob in (0.94, 0.5):
band = future_bands.sel(prob=prob).item()
band.set_label(hdi_label(prob, prefix="forecast "))
handles.append(band)
for i, name in enumerate(names):
ax = pc.get_target("t", {"series": name})
(obs_line,) = ax.plot(x_train, observed[:, i], color="black", lw=1.2, label="observed")
ax.axhline(0.0, color="gray", lw=0.8, ls="--")
ax.set_title(name, fontsize=11)
locator = mdates.AutoDateLocator()
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(mdates.ConciseDateFormatter(locator))
handles.append(obs_line)
pc.get_target("t", {"series": names[0]}).legend(
handles=handles, loc="center left", bbox_to_anchor=(1, 0.5), fontsize=9
)
fig = pc.viz["figure"].item()
fig.supxlabel("date")
fig.supylabel("growth rate (percent)")
fig.suptitle(title, fontsize=18, fontweight="bold", y=1.02)
train_pp = stack_draws("posterior_predictive", tree)
crps_train = eval_crps(train_pp, data)
print(f"in-sample CRPS: {float(crps_train):.4f}")
```
in-sample CRPS: 0.9674
``` python
plot_fit_and_forecast(
train_pp, None, title=f"One-step-ahead in-sample fit (CRPS: {float(crps_train):.3f})"
)
```
# Forecast
The `predictions` group of the tree holds the 30-quarter forecast paths. Each path draws its own correlated shocks and feeds them back through the lag window, so the uncertainty compounds over the horizon. We show the last 40 in-sample quarters for context.
``` python
forecast_draws = stack_draws("predictions", tree)
plot_fit_and_forecast(train_pp, forecast_draws, title="VAR(2) forecast, 30 quarters", n_last=40)
```
The forecast bands widen over the first few quarters and then settle. For a stable VAR the forecast error covariance \sum\_{s \< h} \Psi_s \Sigma \Psi_s^\top converges to the unconditional covariance of the process, and the forecast mean converges to the unconditional mean (I - \sum_l \Phi_l)^{-1} c. The posterior bands also carry parameter uncertainty, so they are a mixture over draws, but with every draw stable the same picture holds. The table shows the width of the 94\\ HDI per series at a few horizons, and the printout compares the unconditional mean implied by the posterior means with the sample means and the mean forecast at the last horizon.
``` python
def hdi_width(draws: Float[np.ndarray, " sample time series"], prob: float) -> np.ndarray:
"""Width of the HDI of ``draws`` per time step and series."""
da = xr.DataArray(np.asarray(draws), dims=["sample", "time", "series"])
hdi = az.hdi(da, prob=prob, dim="sample") # (time, series, ci_bound)
return (hdi.sel(ci_bound="upper") - hdi.sel(ci_bound="lower")).to_numpy()
width_94 = hdi_width(forecast_draws, 0.94)
horizons = [1, 5, 10, 20, 30]
width_94_h = width_94[[h - 1 for h in horizons]]
pl.DataFrame({"h": horizons} | dict(zip(names, width_94_h.T, strict=True))).with_columns(
pl.col(names).round(3)
)
```
shape: (5, 4)
| h | realgdp | realcons | realinv |
|-----|---------|----------|---------|
| i64 | f64 | f64 | f64 |
| 1 | 2.897 | 2.499 | 14.846 |
| 5 | 3.219 | 2.735 | 16.676 |
| 10 | 3.249 | 2.697 | 16.686 |
| 20 | 3.14 | 2.717 | 15.905 |
| 30 | 3.094 | 2.735 | 15.976 |
``` python
phi_mean = np.asarray(posterior["phi"]).mean(axis=0)
c_mean = np.asarray(posterior["intercept"]).mean(axis=0)
unconditional_mean = np.linalg.solve(np.eye(k) - phi_mean.sum(axis=0), c_mean)
pl.DataFrame(
{
"series": names,
"unconditional mean": unconditional_mean,
"sample mean": np.asarray(data).mean(axis=0),
"mean forecast at h=30": forecast_draws.mean(axis=0)[-1],
}
).with_columns(pl.exclude("series").round(3))
```
shape: (3, 4)
| series | unconditional mean | sample mean | mean forecast at h=30 |
|------------|--------------------|-------------|-----------------------|
| str | f64 | f32 | f32 |
| "realgdp" | 0.789 | 0.772 | 0.79 |
| "realcons" | 0.837 | 0.832 | 0.824 |
| "realinv" | 0.947 | 0.818 | 1.02 |
# Impulse response functions
A forecast tells you where the system goes on average. An impulse response tells you how a shock to one series propagates to all series over time. For a stable VAR the moving-average (Wold) representation
y_t = \mu + \sum\_{h=0}^{\infty} \Psi_h \\ \varepsilon\_{t-h}
exists, and the coefficient matrices follow the recursion
\Psi_0 = I, \qquad \Psi_h = \sum\_{j=1}^{\min(h, p)} \Phi_j \\ \Psi\_{h-j} \quad (h \geq 1).
The entry \Psi_h\[i, j\] is the response of series i, h quarters after a unit shock to the reduced-form residual \varepsilon\_{t, j}, with the other residuals held at zero. `impulse_response(phi, horizon)` runs this recursion for all posterior draws at once (the draws pass through the leading batch axis; no `vmap` is needed) and returns an array of shape `(draws, horizon + 1, series, series)`, indexed as `[draw, h, response, shock]`.
The recursion exists for any coefficients, but the representation and the decay \Psi_h \to 0 need stability, which we checked above. If some draws were unstable we would mask them here; the mask below is the identity when all draws are stable.
Writing out the recursion for our VAR(2) (p = 2) makes it concrete. The sum only ever has one or two terms, because \min(h, p) \leq 2:
\Psi_0 = I, \qquad \Psi_1 = \Phi_1 \\ \Psi_0 = \Phi_1, \qquad \Psi_2 = \Phi_1 \\ \Psi_1 + \Phi_2 \\ \Psi_0 = \Phi_1^2 + \Phi_2, \qquad \Psi_3 = \Phi_1 \\ \Psi_2 + \Phi_2 \\ \Psi_1.
Each \Psi_h only ever combines the two lag matrices \Phi_1, \Phi_2 (sampled once per posterior draw) with the previously computed \Psi\_{h-1}, \Psi\_{h-2}: this is exactly what [impulse_response](../../reference/var.impulse_response.md#numpyro_forecast.var.impulse_response) scans over, and it is also what [companion_matrix](../../reference/var.companion_matrix.md#numpyro_forecast.var.companion_matrix) block-multiplies in one shot when we only need the stability check rather than every intermediate \Psi_h.
``` python
n_irf_steps = 10
irf_labels = [f"{response} response to {shock} shock" for response in names for shock in names]
irf_draws = impulse_response(phi_draws[stable], n_irf_steps) # (draws, 11, 3, 3)
print(f"irf_draws: {irf_draws.shape}")
print("posterior mean responses at h = 0, 1, 2 (rows: response, columns: shock):")
print(np.round(np.asarray(irf_draws.mean(axis=0)[:3]), 3))
```
irf_draws: (4000, 11, 3, 3)
posterior mean responses at h = 0, 1, 2 (rows: response, columns: shock):
[[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]
[[-0.06 0.443 0.01 ]
[-0.057 0.228 0.021]
[-0.498 2.827 0.072]]
[[-0.028 0.367 0.008]
[-0.134 0.314 0.029]
[ 0.055 1.27 -0.009]]]
\Psi_0 is the identity, as expected: at h = 0 every series responds only to its own shock, one for one. \Psi_1 and \Psi_2 show how the shock starts to spread to the other two series through the estimated \Phi_1, \Phi_2 coefficients; the full grid below plots this spread out to h = 10 with posterior uncertainty.
``` python
def plot_irf_grid(
irf: Float[Array, " sample steps series series"],
title: str,
ylabel: str,
overlay: Float[Array, " sample steps series series"] | None = None,
overlay_label: str = "",
legend_loc: str = "upper right",
) -> None:
"""Plot a ``series x series`` grid of impulse responses with HDI bands.
Parameters
----------
irf
Impulse response draws ``(sample, steps, response, shock)``.
title
Figure title.
ylabel
Shared y-axis label.
overlay
Optional second set of draws whose posterior mean is overlaid as a line.
overlay_label
Legend label of the overlaid mean.
legend_loc
Legend location, forwarded to `matplotlib.axes.Axes.legend`.
"""
n_draws, n_steps = irf.shape[:2]
steps = np.arange(n_steps, dtype=float)
pc = az.plot_lm(
predictions_to_datatree(np.asarray(irf).reshape(n_draws, n_steps, -1), steps, irf_labels),
y="obs",
x="t",
plot_dim="time",
ci_kind="hdi",
ci_prob=hdi_probs,
smooth=False,
point_estimate="mean",
col_wrap=3,
visuals={
"ci_band": {"color": "C0"},
"observed_scatter": False,
"pe_line": {"color": "C0", "alpha": 1.0, "width": 1.5},
"xlabel": False,
"ylabel": False,
},
aes={"alpha": ["prob"]},
alpha=hdi_alphas,
figure_kwargs={"figsize": (14, 10)},
)
bands = pc.viz["ci_band"]["t"].sel(series=irf_labels[0])
handles = []
for prob in (0.94, 0.5):
band = bands.sel(prob=prob).item()
band.set_label(hdi_label(prob))
handles.append(band)
mean_line = pc.viz["pe_line"]["t"].sel(series=irf_labels[0]).item()
mean_line.set_label("posterior mean")
handles.append(mean_line)
if overlay is not None:
az.plot_lm(
predictions_to_datatree(
np.asarray(overlay).reshape(overlay.shape[0], n_steps, -1), steps, irf_labels
),
y="obs",
x="t",
plot_dim="time",
plot_collection=pc,
ci_kind="hdi",
ci_prob=hdi_probs,
smooth=False,
point_estimate="mean",
visuals={
"ci_band": False,
"observed_scatter": False,
"pe_line": {"color": "C1", "alpha": 1.0, "width": 1.5},
"xlabel": False,
"ylabel": False,
},
)
overlay_line = pc.viz["pe_line"]["t"].sel(series=irf_labels[0]).item()
overlay_line.set_label(overlay_label)
handles.append(overlay_line)
for label in irf_labels:
ax = pc.get_target("t", {"series": label})
ax.axhline(0.0, color="gray", lw=0.8, ls="--")
ax.set_title(label, fontsize=10)
pc.get_target("t", {"series": irf_labels[0]}).legend(
handles=handles, loc=legend_loc, fontsize=8
)
fig = pc.viz["figure"].item()
fig.supxlabel("quarters after the shock")
fig.supylabel(ylabel)
fig.suptitle(title, fontsize=18, fontweight="bold", y=1.02)
plot_irf_grid(
irf_draws,
title="Impulse responses to a unit reduced-form shock",
ylabel="response (percentage points)",
)
```
The `realinv response to realcons shock` panel dominates the grid: its y-axis runs past 3, while every other off-diagonal panel stays within about \pm 0.6. The printed h = 0, 1, 2 snapshot above puts a number on it: a one-unit consumption-growth shock moves investment growth by 2.827 percentage points at h = 1 and 1.27 at h = 2, six to ten times larger than any other off-diagonal entry at the same horizons, consistent with investment's well-known sensitivity to demand shocks. Every panel decays toward zero by h \approx 8-10, as stability requires; the own-shock (diagonal) panels start at exactly 1 by construction and decay the fastest, some dipping briefly negative (e.g.Β `realgdp response to realgdp shock` around h = 1-2) before settling.
## Orthogonalized and cumulative responses
A unit shock to one reduced-form residual with the others held at zero is not an experiment we can observe when \Sigma is not diagonal: the residuals move together. The standard fix is to rewrite the shocks as \varepsilon_t = L \\ u_t with u_t \sim \text{MultivariateNormal}(0, I) and L the Cholesky factor of \Sigma, and to report the responses to the *orthogonalized* shocks u_t:
\Theta_h = \Psi_h \\ L.
A unit shock to u\_{t, j} is a one-standard-deviation shock. Because L is lower triangular, the first series in the ordering responds only to its own shock in the impact quarter, the second series to the first two shocks, and so on. This is the recursive identification, and the ordering `realgdp, realcons, realinv` is part of the model: a different ordering gives different orthogonalized responses. [impulse_response](../../reference/var.impulse_response.md#numpyro_forecast.var.impulse_response) takes the factor through `scale_tril`, here built from the posterior draws of \sigma and L\_\Omega.
Our series are growth rates in percent, g_t = 100 \\ \Delta \log Y_t, so the running sum \sum\_{s=0}^{h} \Theta_s is the response of the *log level* in percent, approximately the percent change of the level. We request it with `cumulative=True` and extend the horizon to 20 quarters. For a stable VAR the cumulative response converges to the long-run effect (I - \sum_l \Phi_l)^{-1} L.
``` python
# Reassemble the Cholesky factor L = diag(sigma) @ L_Omega per draw, the same construction
# the model uses for the likelihood's scale_tril, so the orthogonalized shocks are one
# posterior-consistent standard deviation of the fitted shock covariance.
sigma_draws = jnp.asarray(posterior["sigma"])[stable]
l_omega_draws = jnp.asarray(posterior["l_omega"])[stable]
scale_tril_draws = sigma_draws[..., :, None] * l_omega_draws # (draws, 3, 3)
# cumulative=True sums Theta_h = Psi_h @ L over h, turning the growth-rate response into a
# level response; horizon 20 is long enough for the stable draws to approach that long-run sum.
irf_level_draws = impulse_response(
phi_draws[stable], 20, scale_tril=scale_tril_draws, cumulative=True
)
print(f"irf_level_draws: {irf_level_draws.shape}")
print("posterior mean cumulative response at h = 20 (rows: response, columns: shock):")
print(np.round(np.asarray(irf_level_draws.mean(axis=0)[-1]), 3))
```
irf_level_draws: (4000, 21, 3, 3)
posterior mean cumulative response at h = 20 (rows: response, columns: shock):
[[1.26 0.612 0.139]
[0.762 0.897 0.174]
[5.187 1.377 2.671]]
``` python
plot_irf_grid(
irf_level_draws,
title="Cumulative responses to a one standard deviation orthogonalized shock",
ylabel="level response (percent)",
legend_loc="lower right",
)
```
# Minnesota prior
A VAR has many coefficients for its sample size: here 18 lag coefficients plus 3 intercepts for 200 quarters, and the count grows with p k^2. Unregularized fits overfit and forecast poorly. Litterman (1986) and Doan, Litterman and Sims (1984) proposed the *Minnesota prior*, which encodes three beliefs:
1. Each series is close to a univariate process: the prior mean of the first own lag is m\_{\text{own}} and every other coefficient is centered at zero. For series in levels m\_{\text{own}} = 1 (a random walk); for differenced or otherwise stationary series, as here, m\_{\text{own}} = 0.
2. Longer lags matter less: the prior standard deviation decays with the lag, d(l) = 1/l (harmonic) or 1/l^2.
3. Other series matter less than the own past: cross-variable coefficients get a tighter prior by a factor \kappa \in \[0, 1\].
Together, for the coefficient of series j at lag l in the equation of series i,
\Phi\_{l, ij} \sim \text{Normal}\left(m\_{l, ij}, \\ \lambda \\ d(l) \\ \kappa^{\[i \neq j\]}\right), \qquad m\_{l, ij} = m\_{\text{own}} \\ \[l = 1\] \\ \[i = j\],
with an overall tightness \lambda. `minnesota_prior(n_lags, n_obs, tightness, cross_shrinkage, decay, own_lag_mean)` returns the `loc` and `scale` arrays in the `(lags, series, series)` layout of [var_step](../../reference/var.var_step.md#numpyro_forecast.var.var_step), and we pass them to `dist.Normal(...).to_event(3)`. Nothing else changes: the prior is an argument of the model factory.
This parameterization follows [Impulso's `MinnesotaPrior`](https://thomaspinder.github.io/Impulso/reference/generated/impulso.priors.MinnesotaPrior.html): the same three knobs (`tightness`, `decay`, `cross_shrinkage`) and a tightness that is fixed rather than estimated. There is no closed-form marginal likelihood for the independent-normal prior, so Impulso treats the tightness as a modeling choice, and so do we. Two differences: Impulso calls the 1/l^2 decay "geometric" (in Doan, Litterman and Sims it is the harmonic decay with exponent two), and it fixes the own-lag mean at one, while we expose `own_lag_mean` because differenced data call for zero.
## Series on different scales
The classic Litterman formulation multiplies the standard deviation of the cross-variable coefficients by \sigma_i / \sigma_j, the ratio of the residual standard deviations of the two series. Impulso omits this factor and asks for pre-scaled data. Our series are not on a common scale: investment growth is about five times more volatile than GDP or consumption growth, so the investment equation carries coefficients about five times larger, and a common tightness would shrink them five times too hard. We apply the classic correction with the sample standard deviations as a proxy for the residual scales. The table shows the ratios \sigma_i / \sigma_j (rows: equation, columns: lagged series).
``` python
sample_sd = y_pct.select(names).std().to_numpy()[0]
scale_ratio = sample_sd[:, None] / sample_sd[None, :]
pl.DataFrame({"equation": names} | dict(zip(names, scale_ratio.T, strict=True))).with_columns(
pl.col(names).round(2)
)
```
shape: (3, 4)
| equation | realgdp | realcons | realinv |
|------------|---------|----------|---------|
| str | f64 | f64 | f64 |
| "realgdp" | 1.0 | 1.27 | 0.19 |
| "realcons" | 0.79 | 1.0 | 0.15 |
| "realinv" | 5.33 | 6.75 | 1.0 |
With \lambda = 0.5 and \kappa = 0.5 the prior standard deviation of a first-lag cross coefficient is 0.25 before scaling. The printout compares, for the investment equation at lag one, the unscaled (Impulso-style) prior standard deviations, the scaled ones we use, and the posterior under the weak prior: the coefficient on lagged consumption growth has a posterior mean near 2.8 with a standard deviation near 0.5, so an unscaled prior with standard deviation 0.25 sits more than ten prior standard deviations away from what the data say and would dominate the posterior, while the scaled prior is compatible with it.
``` python
tightness = 0.5
loc_mn, scale_unscaled = minnesota_prior(
p, k, tightness=tightness, cross_shrinkage=0.5, own_lag_mean=0.0
)
scale_mn = scale_unscaled * jnp.asarray(scale_ratio, dtype=scale_unscaled.dtype)
minnesota = dist.Normal(loc_mn, scale_mn).to_event(3)
phi_weak_mean = np.asarray(posterior["phi"]).mean(axis=0)
phi_weak_sd = np.asarray(posterior["phi"]).std(axis=0)
pl.DataFrame(
{
"lagged series": names,
"unscaled prior sd": np.asarray(scale_unscaled[0, 2]),
"scaled prior sd": np.asarray(scale_mn[0, 2]),
"weak prior posterior mean": phi_weak_mean[0, 2],
"weak prior posterior sd": phi_weak_sd[0, 2],
}
).with_columns(pl.exclude("lagged series").round(3))
```
shape: (3, 5)
| lagged series | unscaled prior sd | scaled prior sd | weak prior posterior mean | weak prior posterior sd |
|----|----|----|----|----|
| str | f32 | f32 | f32 | f32 |
| "realgdp" | 0.25 | 1.331 | -0.498 | 0.605 |
| "realcons" | 0.25 | 1.687 | 2.827 | 0.498 |
| "realinv" | 0.5 | 0.5 | 0.072 | 0.104 |
``` python
var_model_mn = make_var_model(minnesota, y_init)
rng_key, rng_subkey = random.split(rng_key)
mcmc_mn = fit_nuts(rng_subkey, var_model_mn, data, covariates_train)
posterior_mn = mcmc_mn.get_samples()
print(f"divergences: {n_divergences(mcmc_mn)}")
rng_key, rng_subkey = random.split(rng_key)
tree_mn = export(rng_subkey, var_model_mn, posterior_mn)
summary_mn = az.summary(
tree_mn, var_names=["intercept", "sigma", "phi"], ci_kind="hdi", ci_prob=0.94
)
summary_mn
```
divergences: 0
| | mean | sd | hdi94_lb | hdi94_ub | ess_bulk | ess_tail | r_hat | mcse_mean | mcse_sd |
|----|----|----|----|----|----|----|----|----|----|
| intercept\[realgdp\] | 0.25 | 0.092 | 0.079 | 0.43 | 2494 | 2731 | 1.00 | 0.0018 | 0.0013 |
| intercept\[realcons\] | 0.548 | 0.088 | 0.38 | 0.71 | 3226 | 3077 | 1.00 | 0.0015 | 0.0011 |
| intercept\[realinv\] | -1.69 | 0.46 | -2.6 | -0.81 | 3052 | 2944 | 1.00 | 0.0084 | 0.0061 |
| sigma\[realgdp\] | 0.746 | 0.0353 | 0.68 | 0.82 | 3335 | 2824 | 1.00 | 0.00062 | 0.00044 |
| sigma\[realcons\] | 0.658 | 0.0331 | 0.6 | 0.72 | 4268 | 3306 | 1.00 | 0.00051 | 0.00038 |
| sigma\[realinv\] | 3.867 | 0.183 | 3.5 | 4.2 | 3918 | 2927 | 1.00 | 0.0029 | 0.0021 |
| phi\[1, realgdp, realgdp\] | -0.056 | 0.119 | -0.28 | 0.17 | 2048 | 2568 | 1.00 | 0.0026 | 0.0019 |
| phi\[1, realgdp, realcons\] | 0.467 | 0.102 | 0.28 | 0.66 | 2304 | 2662 | 1.00 | 0.0021 | 0.0015 |
| phi\[1, realgdp, realinv\] | 0.0117 | 0.0183 | -0.023 | 0.045 | 2207 | 2590 | 1.00 | 0.00039 | 0.00028 |
| phi\[1, realcons, realgdp\] | 0.001 | 0.104 | -0.2 | 0.19 | 2800 | 2672 | 1.00 | 0.002 | 0.0015 |
| phi\[1, realcons, realcons\] | 0.195 | 0.094 | 0.017 | 0.37 | 3085 | 2960 | 1.00 | 0.0017 | 0.0012 |
| phi\[1, realcons, realinv\] | 0.0155 | 0.0162 | -0.015 | 0.046 | 2764 | 2955 | 1.00 | 0.00031 | 0.00023 |
| phi\[1, realinv, realgdp\] | -0.86 | 0.65 | -2.1 | 0.4 | 2067 | 2401 | 1.00 | 0.014 | 0.011 |
| phi\[1, realinv, realcons\] | 3.27 | 0.56 | 2.2 | 4.3 | 2496 | 2707 | 1.00 | 0.011 | 0.0078 |
| phi\[1, realinv, realinv\] | 0.119 | 0.101 | -0.077 | 0.31 | 2228 | 2297 | 1.00 | 0.0021 | 0.0016 |
| phi\[2, realgdp, realgdp\] | 0.06 | 0.084 | -0.096 | 0.22 | 2956 | 3155 | 1.00 | 0.0015 | 0.0011 |
| phi\[2, realgdp, realcons\] | 0.168 | 0.083 | 0.0083 | 0.32 | 2534 | 2936 | 1.00 | 0.0017 | 0.0011 |
| phi\[2, realgdp, realinv\] | -0.0078 | 0.0134 | -0.033 | 0.018 | 3293 | 2976 | 1.00 | 0.00023 | 0.00017 |
| phi\[2, realcons, realgdp\] | -0.004 | 0.07 | -0.13 | 0.13 | 3726 | 3191 | 1.00 | 0.0011 | 0.00083 |
| phi\[2, realcons, realcons\] | 0.126 | 0.085 | -0.028 | 0.29 | 3347 | 2830 | 1.00 | 0.0015 | 0.001 |
| phi\[2, realcons, realinv\] | 0.0083 | 0.0114 | -0.014 | 0.03 | 4401 | 2756 | 1.00 | 0.00017 | 0.00013 |
| phi\[2, realinv, realgdp\] | 0.19 | 0.44 | -0.64 | 1 | 2776 | 2965 | 1.00 | 0.0084 | 0.0059 |
| phi\[2, realinv, realcons\] | 0.44 | 0.44 | -0.36 | 1.3 | 3003 | 2792 | 1.00 | 0.0081 | 0.0057 |
| phi\[2, realinv, realinv\] | -0.06 | 0.079 | -0.21 | 0.088 | 3063 | 3013 | 1.00 | 0.0014 | 0.001 |
## Shrinkage of the coefficients
The table compares the posterior standard deviation of every coefficient under the two priors. The ratio column is below one where the Minnesota prior tightened the posterior. The effect is largest on the second lag, where the harmonic decay halves the prior standard deviation, and on the GDP and consumption equations. The investment equation at lag one is unchanged within Monte Carlo error: after the scale correction its prior is wide relative to what the data say, so the data decide.
``` python
lag_labels, equation_labels, lagged_labels = zip(
*itertools.product(range(1, p + 1), names, names), strict=True
)
phi_sd = pl.DataFrame(
{
"lag": lag_labels,
"equation": equation_labels,
"lagged series": lagged_labels,
"weak prior": phi_weak_sd.reshape(-1),
"minnesota prior": np.asarray(posterior_mn["phi"]).std(axis=0).reshape(-1),
}
).with_columns((pl.col("minnesota prior") / pl.col("weak prior")).alias("ratio"))
own_mean = phi_sd.filter(pl.col("equation") == pl.col("lagged series"))["ratio"].mean()
cross_mean = phi_sd.filter(pl.col("equation") != pl.col("lagged series"))["ratio"].mean()
print(f"mean ratio on own lags: {own_mean:.2f}, on cross lags: {cross_mean:.2f}")
phi_sd.with_columns(pl.exclude("lag", "equation", "lagged series").round(3))
```
mean ratio on own lags: 0.78, on cross lags: 0.76
shape: (18, 6)
| lag | equation | lagged series | weak prior | minnesota prior | ratio |
|-----|------------|---------------|------------|-----------------|-------|
| i64 | str | str | f32 | f32 | f32 |
| 1 | "realgdp" | "realgdp" | 0.141 | 0.119 | 0.843 |
| 1 | "realgdp" | "realcons" | 0.114 | 0.102 | 0.894 |
| 1 | "realgdp" | "realinv" | 0.023 | 0.018 | 0.804 |
| 1 | "realcons" | "realgdp" | 0.146 | 0.104 | 0.713 |
| 1 | "realcons" | "realcons" | 0.115 | 0.094 | 0.818 |
| β¦ | β¦ | β¦ | β¦ | β¦ | β¦ |
| 2 | "realcons" | "realcons" | 0.123 | 0.085 | 0.69 |
| 2 | "realcons" | "realinv" | 0.022 | 0.011 | 0.525 |
| 2 | "realinv" | "realgdp" | 0.621 | 0.442 | 0.712 |
| 2 | "realinv" | "realcons" | 0.541 | 0.442 | 0.818 |
| 2 | "realinv" | "realinv" | 0.103 | 0.079 | 0.765 |
## Forecast bands
Tighter coefficients mean less parameter uncertainty in the forecast. The table reports the mean width of the 94\\ HDI over the 30 forecast quarters, per series and per prior. The change is small: a few percent for GDP and consumption and none for investment. With 200 quarters for 18 coefficients, the forecast uncertainty comes from the shock covariance, not from the coefficients, and a prior of this tightness cannot move it much.
``` python
forecast_draws_mn = stack_draws("predictions", tree_mn)
pl.DataFrame(
{
"series": names,
"weak prior": width_94.mean(axis=0),
"minnesota prior": hdi_width(forecast_draws_mn, 0.94).mean(axis=0),
}
).with_columns(pl.exclude("series").round(3))
```
shape: (3, 3)
| series | weak prior | minnesota prior |
|------------|------------|-----------------|
| str | f64 | f64 |
| "realgdp" | 3.21 | 3.168 |
| "realcons" | 2.705 | 2.63 |
| "realinv" | 16.489 | 16.536 |
## Impulse responses
The grid overlays the posterior mean responses under the Minnesota prior (orange) on the bands and mean of the weak prior fit (blue). The two means agree closely and the orange line stays inside the 50\\ band of the weak prior fit in every panel. The Minnesota mean is smoother at two and three quarters after the shock, where the prior halves the standard deviation of the second-lag coefficients and irons out the wiggle that the weak prior fit shows there. With 200 quarters for 18 coefficients the data dominate a prior of this tightness; a smaller `tightness` trades this agreement for more shrinkage.
``` python
phi_draws_mn = jnp.asarray(posterior_mn["phi"])
radius_mn = np.abs(np.linalg.eigvals(np.asarray(companion_matrix(phi_draws_mn)))).max(axis=-1)
stable_mn = radius_mn < 1.0
print(f"stable draws (Minnesota prior): {stable_mn.mean():.3f}")
irf_draws_mn = impulse_response(phi_draws_mn[stable_mn], n_irf_steps)
plot_irf_grid(
irf_draws,
title="Impulse responses: weak prior (bands) vs Minnesota prior (orange mean)",
ylabel="response (percentage points)",
overlay=irf_draws_mn,
overlay_label="posterior mean (Minnesota prior)",
)
```
stable draws (Minnesota prior): 1.000
The tightness \lambda is a modeling choice, not an estimate: there is no closed-form marginal likelihood for the independent-normal prior to optimize it, so pick it from the scale of the coefficients you find plausible, or compare forecast scores across a few values with [backtest](../../reference/evaluate.backtest.md#numpyro_forecast.evaluate.backtest).
# References
- Orduz, J. [Bayesian VAR in NumPyro](https://juanitorduz.github.io/var_numpyro/). The source of this notebook.
- LΓΌtkepohl, H. (2005). *New Introduction to Multiple Time Series Analysis*. Springer. Chapters 2 and 5 cover the moving-average representation, impulse responses and 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.
- Doan, T., Litterman, R. B. and Sims, C. A. (1984). Forecasting and conditional projection using realistic prior distributions. *Econometric Reviews*, 3(1), 1-100.
- Pinder, T. [Impulso](https://github.com/thomaspinder/impulso): a Bayesian VAR package for Python ([documentation](https://thomaspinder.github.io/Impulso/), [`MinnesotaPrior` reference](https://thomaspinder.github.io/Impulso/reference/generated/impulso.priors.MinnesotaPrior.html)). The Minnesota prior parameterization and the batched moving-average recursion used here follow its design.
- statsmodels. [macrodata](https://www.statsmodels.org/stable/datasets/generated/macrodata.html): United States macroeconomic data, 1959Q1 to 2009Q3, public domain.
[Source: Vector Autoregression (VAR) with `numpyro_forecast`](_src/var-preview.html#50138624)
array([0, 1, 2, 3])
draw
(draw)
int64
0 1 2 3 4 5 ... 995 996 997 998 999
