The conjugate VAR: fast Bayesian estimation

When to reach for ConjugateVAR instead of the NUTS VAR

Impulso estimates a reduced-form VAR two ways. Both return the same FittedVAR, so identification, impulse responses, FEVDs, and forecasts are byte-for-byte identical downstream. What differs is the mode of inference:

• The NUTS VAR (VAR) places independent-Normal priors on the coefficients and samples the full posterior with Hamiltonian Monte Carlo. Maximally flexible — it admits per-equation shrinkage, stochastic volatility, sign restrictions, and external instruments — but every coefficient is a sampled latent, so large systems are slow.

• The conjugate VAR (ConjugateVAR) places a Normal-Inverse-Wishart prior, which is conjugate to the VAR likelihood. The coefficient/covariance posterior is then available in closed form: we draw \((\beta, \Sigma)\) analytically and reserve Monte Carlo for a single low-dimensional hyperparameter — the Minnesota tightness \(\lambda\) — which the data selects by marginal likelihood (Giannone et al. [2015]).

This notebook fits both on the same series, shows they reach the same structural conclusions, times them, and ends with a rule for choosing between them. We use an environmental system — a German climate–energy VAR — rather than the usual macro data, to show the machinery is domain-agnostic.

Scope

This is the estimator-first tour. For the conjugate VAR wearing a deterministic volatility break — the COVID application it was built for — see Estimating a VAR after March 2020. Why a conjugate estimator is a sibling of VAR rather than a mode of it is recorded in ADR 0004.

import time

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from impulso import Cholesky, ConjugateVAR, MinnesotaPrior, NIWPrior, VAR, VARData
from impulso.samplers import NUTSSampler

Data: a German climate–energy system

Weather is where the environment meets the economy: temperature drives heating and cooling demand, sunshine and wind set renewable supply, and rainfall feeds hydro and agriculture. We assemble a four-variable monthly system for Berlin (1980–2024) from the ERA5 reanalysis (Hersbach and others [2020]), served by the Open-Meteo historical archive.

Variable	ERA5 series (monthly mean of daily)	Unit	Economic reading
temperature	temperature_2m_mean	°C	heating / cooling demand
radiation	shortwave_radiation_sum	MJ/m²	solar-PV potential
wind	wind_speed_10m_mean	km/h	wind-power potential
precipitation	precipitation_sum	mm/day	hydro / runoff

The committed CSV is produced once by scripts/fetch_berlin_climate.py and read offline here — no network call at render time. That script targets Open-Meteo’s free archive endpoint, so anyone can reproduce the file without credentials.

Fig. 7 Raw monthly ERA5 series for Berlin, 1980–2024. Temperature and radiation are dominated by the seasonal cycle.

The raw series are overwhelmingly seasonal — a VAR fit on them would spend its coefficients re-learning the calendar. We model anomalies instead: each observation minus its month-of-year climatological mean, standardised to unit variance. The result is stationary, comparable across variables, and lets impulse responses read in standard deviations.

climatology = raw.groupby(raw.index.month).transform("mean")
anomalies = raw - climatology
anomalies = (anomalies - anomalies.mean()) / anomalies.std()
anomalies.describe().round(2)

	temperature	radiation	wind	precipitation
count	540.00	540.00	540.00	540.00
mean	0.00	0.00	-0.00	0.00
std	1.00	1.00	1.00	1.00
min	-4.30	-4.51	-2.74	-2.10
25%	-0.58	-0.45	-0.69	-0.67
50%	0.00	-0.01	-0.05	-0.13
75%	0.63	0.43	0.56	0.54
max	2.79	3.72	4.02	6.04

Fig. 8 Standardised monthly anomalies — the seasonal cycle removed. This is what the VAR sees.

Fitting the conjugate VAR

We use twelve lags — enough to capture up to a year of dynamic feedback in monthly data — giving \(4 \times 12 = 48\) coefficients per equation. The prior is the conjugate Minnesota prior NIWPrior; select=True asks the estimator to choose the overall tightness \(\lambda\) by maximising the marginal likelihood and then sample its posterior, rather than fixing it by hand (Giannone et al. [2015]; the Minnesota shrinkage idea goes back to Doan et al. [1984] and Litterman [1986]).

LAGS = 12
data = VARData.from_df(anomalies, endog=list(anomalies.columns))

conjugate_prior = NIWPrior(select=True, decay=2.0, cross_shrinkage=1.0)
start = time.perf_counter()
fitted_conjugate = ConjugateVAR(lags=LAGS, prior=conjugate_prior, draws=2000, tune=1000, seed=0).fit(data)
conjugate_seconds = time.perf_counter() - start

lambda_hat = float(fitted_conjugate.idata.posterior["lambda_"].median())
print(f"data-selected Minnesota tightness  lambda = {lambda_hat:.3f}")
print(f"conjugate fit wall-clock                  = {conjugate_seconds:.2f} s")

data-selected Minnesota tightness  lambda = 0.271
conjugate fit wall-clock                  = 3.82 s

The estimator reports a posterior for \(\lambda\) (not a fixed value): the data speak to how much shrinkage the system needs. Everything downstream — coefficients, covariance, the base Cholesky factor — was drawn in closed form conditional on those hyperparameter draws.

The same model by NUTS

To make this a clean inference-mode comparison, we fit the NUTS VAR at the same tightness the conjugate estimator just selected (MinnesotaPrior(tightness=lambda_hat)). Now the only differences are the prior family (independent-Normal vs conjugate NIW) and the sampler — not the amount of shrinkage.

if ci:
    sampler = NUTSSampler(draws=10, tune=50, chains=1, cores=1, random_seed=0)
else:
    sampler = NUTSSampler(draws=1000, tune=1500, chains=2, cores=1, random_seed=0)

start = time.perf_counter()
fitted_nuts = VAR(lags=LAGS, prior=MinnesotaPrior(tightness=lambda_hat)).fit(data, sampler=sampler)
nuts_seconds = time.perf_counter() - start
print(f"NUTS fit wall-clock = {nuts_seconds:.2f} s")
az.summary(fitted_nuts.idata, var_names=["intercept"], kind="diagnostics")

NUTS fit wall-clock = 69.96 s

	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
intercept[0]	0.001	0.001	3916.0	1763.0	1.01
intercept[1]	0.001	0.001	3430.0	1530.0	1.00
intercept[2]	0.001	0.001	3190.0	1759.0	1.00
intercept[3]	0.001	0.001	3421.0	1689.0	1.00

Speed

Both estimators fit the identical 4-variable, 12-lag system. The conjugate path spends Monte Carlo only on a single scalar; NUTS explores a ~200-dimensional coefficient posterior. At full render the gap is an order of magnitude or more.

if ci:
    print("CI smoke mode: NUTS is shrunk to a few draws — these timings are NOT representative.")
else:
    print(
        f"conjugate: {conjugate_seconds:.2f} s | "
        f"NUTS: {nuts_seconds:.2f} s | "
        f"speed-up: {nuts_seconds / conjugate_seconds:.0f}x"
    )

conjugate: 3.82 s | NUTS: 69.96 s | speed-up: 18x

Downstream parity: identical structural machinery

Because both estimators return a FittedVAR, identification is the same call on each. We apply a Cholesky scheme with the ordering radiation → temperature → wind → precipitation (solar forcing is the most exogenous; rainfall the most responsive). The ordering encodes real assumptions — see Monetary Policy Analysis for how much it can matter — but here we hold it fixed and vary only the estimator.

ordering = ["radiation", "temperature", "wind", "precipitation"]
irf_conjugate = fitted_conjugate.set_identification_strategy(Cholesky(ordering=ordering)).impulse_response(horizon=24)
irf_nuts = fitted_nuts.set_identification_strategy(Cholesky(ordering=ordering)).impulse_response(horizon=24)

Fig. 9 Conjugate-VAR impulse responses (Cholesky). Column shock → row response, over 24 months.

Now overlay the two estimators on the same axes. If the conjugate VAR is a legitimate estimator and not a shortcut, its responses should track the NUTS responses in shape and sign, with band widths of the same order.

def irf_band(irf_result, shock, response, prob=0.9):
    """Return (horizons, median, hdi_low, hdi_high) for one shock→response pair."""
    draws = irf_result.idata.posterior_predictive["irf"].sel(shock=shock, response=response)
    median = draws.median(dim=("chain", "draw")).values
    hdi = az.hdi(draws, hdi_prob=prob)["irf"]
    return np.arange(median.shape[0]), median, hdi.sel(hdi="lower").values, hdi.sel(hdi="higher").values

Fig. 10 Conjugate vs NUTS impulse responses at the same tightness. Medians (lines) and 90% bands (shaded).

The two estimators tell the same structural story: a positive radiation (sunshine) shock warms temperature; a warmth shock is followed by calmer winds. The medians track closely and the bands overlap. They are not identical — the conjugate NIW prior imposes a symmetric Kronecker structure across equations while the NUTS prior is independent-Normal — and that is exactly the point: the inference mode is a modelling choice, not a source of contradiction.

if not ci:
    correlations = [
        np.corrcoef(irf_band(irf_conjugate, shock, response)[1], irf_band(irf_nuts, shock, response)[1])[0, 1]
        for shock in ordering
        for response in ordering
    ]
    print(f"median-IRF shape correlation across all 16 shock/response pairs: {np.nanmean(correlations):.2f}")

median-IRF shape correlation across all 16 shock/response pairs: 0.98

When to reach for which

Both estimators share the entire post-fitting pipeline — identification, IRFs, FEVDs, forecasts — so the choice is purely about the estimation path.

Reach for the conjugate VAR when…	Reach for the NUTS VAR when…
speed matters — hyperparameter selection, model comparison, or many refits	you need per-equation or asymmetric cross-variable shrinkage
you want the tightness \(\lambda\) chosen by the data (hierarchical)	you need stochastic volatility, sign restrictions, or external instruments
the conjugate NIW (symmetric, Kronecker) prior suits the problem	you need arbitrary or non-conjugate priors
the system is large and full MCMC over every coefficient is costly	you want full HMC convergence diagnostics on all coefficients

The conjugate VAR trades flexibility for closed-form speed and a data-driven prior. When your problem fits inside that trade — as macro and climate systems with symmetric Minnesota shrinkage usually do — it is the sharper tool. When you need volatility that moves or priors that bend per equation, the NUTS VAR is there, and everything you build on top is the same.

We currently have some availability for consulting on how Bayesian modelling, vector autoregressions, and impulso can be integrated into your team's macroeconomic, financial, and environmental forecasting work. If this sounds relevant, book an introductory call. These calls are for consulting inquiries only. For technical usage questions and free community support, please use GitHub Discussions and the documentation.

References

The works cited above are collected on the project bibliography page.