← All examples

Multivariate Time Series

System models across many series at once. The centrepiece is a multi-part Bayesian VAR arc, built up module by module from scratch: the reduced-form and structural VAR; large-scale and factor-augmented systems; stochastic-volatility, time-varying-parameter, and Markov-switching dynamics; cointegration (VECM); hierarchical prior selection by marginal likelihood; global-local (horseshoe) adaptive shrinkage; and structural identification across the monetary, oil-market, and fiscal domains — by recursive zeros, sign restrictions, elasticity bounds, heteroskedasticity, and institutional timing. Alongside it sit multivariate GARCH covariance dynamics and a demographic-forecasting pipeline over full age schedules — each cross-checked across Python and R.

What a VAR is, and the three things it is asked for. A vector autoregression is the simplest honest way to model several series that plainly affect one another: regress each variable on several lags of every variable, including itself. Nothing is designated a cause; the system is left symmetric and the data fill in the dynamics. Fitting one is easy — with the same regressors in every equation, ordinary least squares equation by equation is already efficient. Almost everything below is about what to do after that. Three objects are read out of a fitted VAR and recur on every page. An impulse response traces what happens to each variable over the following quarters when one shock hits today — the headline output of every structural application here. A forecast-error variance decomposition (FEVD) asks the same question as an accountant: of the uncertainty about where inflation will be in two years, what fraction is owed to each shock? And a historical decomposition goes further still, attributing an actual episode — the 2003–08 oil run-up, say — to the specific shocks that produced it.

Identification, which is the whole difficulty. A fitted VAR gives residuals that are correlated across equations, and correlated residuals cannot be shocks: if inflation and the interest rate always move together, no combination of them is yet the monetary shock. What the data supply is the covariance matrix Σ of those residuals. What a structural story needs is the impact matrix P saying how each underlying shock hits each variable on impact, and the two are linked only by PP′ = Σ. That equation is short: Σ has m(m+1)/2 distinct entries while P has m², so for three variables the data pin six numbers and nine are needed. Infinitely many structural stories fit the same data equally well, and choosing between them requires a belief brought in from outside. That is identification, it is the reason six of the pages below exist, and the arc works through five genuinely different kinds of belief you might be willing to defend: a timing zero (this variable cannot react within the quarter), a sign (a tightening cannot raise inflation), an outside instrument (a series correlated with this shock and nothing else), a limit (this shock has no permanent effect), a calibrated elasticity known from institutions rather than estimated — and one scheme that imposes nothing at all, letting a change in volatility do the work. None is universally right. The craft is matching the assumption you can defend to the question you are asking, and several pages here deliberately test one scheme against another.

Why a VAR needs shrinking, and what the Minnesota prior does. VARs grow quadratically: m variables and p lags need m(1+mp) coefficients, so the 20-variable system in Part 3 carries 1,620 of them against roughly 240 quarterly observations. Least squares can still fit that, but most of each estimate is noise, and the noise compounds when the model is iterated forward to forecast. The Minnesota prior (Litterman 1986) is the standard discipline: shrink every coefficient toward a deliberately naive benchmark — each variable follows its own past and ignores everything else — and let the data pull away from it only where the evidence is strong. Three hyperparameters govern how hard: overall tightness, how much harder to shrink other variables' lags than a variable's own, and how fast to shrink longer lags. That single knob for overall tightness recurs constantly below, first fixed by hand, then set by a rule as the system grows, then estimated from the data by marginal likelihood, and finally replaced altogether by a prior that learns which individual coefficients matter.

Why identification is the whole problem

This is the real macro system the arc is built on — US GDP growth, inflation and the T-bill — fitted as a VAR(2). The left panel is what the data actually deliver: a cloud of correlated residuals. The three coloured pairs of axes are three different proposals for which combinations of that cloud are the underlying shocks, and each of them reproduces the fitted covariance exactly — to machine precision, 9 × 10⁻¹⁶. No amount of extra data distinguishes them, because they have identical likelihoods, identical fit and identical forecasts.

What the data see: one cloud of residuals residual, gdp_growth residual, inflation every coloured pair of axes reproduces this cloud exactly — the data cannot choose What follows: three different answers quarters after the shock response of inflation same likelihood, same fit, same forecasts — opposite economics

The right panel is why that matters. It shows what each of those three equally-valid choices implies for one economic question — the response of inflation to the first shock. The impact effects are +0.24, +1.52 and −1.55: the three curves disagree about the magnitude by a factor of six and two of them disagree about the sign. Same data, same fit, opposite economics. (The colours match across the two panels.)

So the structural question cannot be settled by fitting better. It is settled by bringing in a belief that picks one set of axes, and the arc below works through six pages of ways to do that — a timing zero, a sign, an outside instrument, a long-run limit, an elasticity calibrated from institutions, and, in Part 10, the one scheme that imposes nothing at all and lets a change in volatility pin the axes instead. Read that page against Part 1 and the point sharpens: the recursive ordering everyone starts with forces one of these lines to be exactly vertical, and when the data are allowed to speak they place it elsewhere.

The arc, and the threads running through it

Thirteen numbered parts, built on one from-scratch engine that each part extends rather than replaces. They are best read as four threads rather than a single line — the pages cross-reference each other along these, not in numerical order.

Two connections are worth making deliberately. Part 10 identifies shocks using exactly the volatility that Part 4 modelled as a nuisance — the same heteroskedasticity, read as signal instead of noise. And Parts 5, 10 and 11 between them span every way a VAR can stop being constant: smooth drift, an imposed break, and regimes the data locate for themselves.

Multivariate Volatility

The cross-asset branch of the volatility arc — modelling the whole time-varying conditional covariance matrix Ht (how assets co-move, and how that rises in crises). Four model families built from scratch and cross-checked: CCC (Bollerslev — the constant-correlation decomposition), DCC (Engle's workhorse — two scalars drive the whole correlation matrix, so it scales; a Bayesian NumPyro version adds an honest posterior), BEKK (Engle–Kroner — covariance modelled directly, guaranteed positive-definite but N² parameters), and factor stochastic volatility (Aguilar–West — the Bayesian capstone, positive-definite and scalable via a low-rank latent factor). Each cross-checked in R (rmgarch, mgarchBEKK, factorstochvol). A closing portfolio horse race shows a dynamic covariance cuts realized portfolio volatility ~28%, harvested in 2008 and 2020.

View example →

Bayesian Vector Autoregression

The multivariate workhorse of empirical macro — each variable on lags of every variable — fit on the classic US monetary system (GDP, inflation, T-bill) three ways: the classical VAR (impulse responses, variance decomposition, and the famous recursive-identification price puzzle), the Minnesota-prior BVAR (Litterman shrinkage), and the Villani steady-state BVAR — whose informative prior on the long-run mean lets you pin the steady state (e.g. a 2% inflation target) rather than revert to a crisis-laden sample average. The steady-state model has the lowest out-of-sample RMSE through 2008, beating OLS, Minnesota, and a random walk. A second application takes the same three models to the Treasury yield curve (3m–30y), where policy pass-through decays up the curve and the steady-state prior barely moves — the flip-side lesson that it only helps when the sample average is unrepresentative of the future. The from-scratch Gibbs engine is my own Meseguer (2010) mortality-forecasting sampler, repurposed for finance and macro; cross-checked in R (vars, BVAR).

View example →

Structural VAR — Shock Identification

Part 2 of the VAR arc: turning the reduced-form BVAR into structural shocks. Part 1's recursive (Cholesky) ordering produced the notorious price puzzle (inflation rising after a tightening); this example cures it with credible identification layered on the same bvar.py posterior. Sign restrictions (Uhlig 2005) impose only the signs economics is sure of and report the whole identified set as a band — dissolving the price puzzle, and decomposing the Treasury curve into level and slope shocks. External-instrument (proxy, Romer–Romer) and long-run Blanchard–Quah restrictions pin the shock two more ways. The lesson — a menu, not a winner: four identifications (recursive, sign, instrument, long-run), each a different defensible belief, all riding one Bayesian posterior. Cross-checked in R (vars::BQ, base-R rotation search).

View example →

Large Bayesian VAR — Shrinkage that Scales

Part 3 of the VAR arc: scale. Parts 1–2 kept the VAR small (3–5 variables) so every coefficient could be read; this example pushes to a 20-variable, 4-lag FRED-QD macro system — 1,620 coefficients against ~240 observations — and reproduces Bańbura, Giannone & Reichlin's (2010) result. Unrestricted OLS breaks down as the system grows (worst forecaster everywhere, 2.5× the benchmark on the funds rate), but the same 20 variables under the Minnesota priortightened as the system grows — collapse the relative RMSE right back into the pack. The lesson: it is a bias–variance problem, not a curse of dimensionality, and shrinkage, not the variable count, is what matters. The BGR in-sample-fit rule for the overall tightness is cross-checked in R against the BVAR package's marginal-likelihood choice — two different mechanisms, the same shrink-with-size answer.

View example →

Bayesian Demographic Forecasting

The applied capstone of the Bayesian VAR example — a full stochastic US population-forecasting pipeline reproducing and extending my SSA work (Meseguer 2010), on the very problem the bvar.py engine was written for. Mortality and fertility are each forecast by a Bayesian VAR over their full age schedules (Minnesota and steady-state priors), benchmarked head-to-head against the classical Lee-Carter method, and validated by expanding-window pseudo-real-time backtests. The forecasts then drive a cohort-component projection of the whole population, with a full-uncertainty backtest (2001–2024) propagating both vital rates' posteriors jointly — showing continued growth and rapid population aging, and why anchoring the long run (the steady-state prior) matters most for fertility. 11 notebooks, from HMD life tables to the population fan chart.

View example →

BVAR with Stochastic Volatility

Part 4 of the VAR arc: time-varying volatility. Every VAR through the large system of Part 3 assumed a constant shock covariance — treating 1974, 1998 and 2008 as equally turbulent. US macro data plainly disagree: the Great Moderation cut output, inflation and rate volatility sharply around 1984, bracketed by the stormy 1970s and 2008. This example adds stochastic volatility (Cogley–Sargent / Primiceri) to the monetary BVAR from scratch in NumPyro — factoring the covariance into a constant contemporaneous matrix and an m-dimensional random-walk log-volatility surface, sampled jointly by NUTS with the coefficients held fixed so the volatility is the whole story. The fitted paths trace the Great Moderation, the Volcker and 1974 peaks, 2008, and the post-ZLB collapse of funds-rate volatility; SV raises the in-sample log-likelihood by ~300 and — the practical payoff — delivers calibrated forecast intervals that widen in storms and narrow in calm, where a fixed-Σ VAR is wrong in both regimes at once. Cross-checked in R with stochvol (AR(1) SV on the OLS residuals) — a different spec and paradigm, the same volatility history.

View example →

Time-Varying-Parameter VAR with Stochastic Volatility

Part 5 of the VAR arc — the capstone. Part 4 let the shock variances move but held the coefficients fixed; Primiceri's (2005) TVP-VAR-SV lets everything drift — the autoregressive coefficients, the contemporaneous relations, and the volatilities all as random walks — to ask whether the monetary transmission mechanism itself changed across the Great Inflation, Volcker, the Great Moderation and 2008. Built from scratch in NumPyro, sampled jointly by NUTS with tight, calibrated Primiceri drift priors (estimated, not fixed) so the data decide how much time-variation is real. Its deliverable — the time-varying impulse response to a policy shock — shows the real effect of a tightening large in the 1970s–80s (peak output −0.72 in 1975), shrinking to −0.15 by the mid-1990s, then deepening again in the 2008 crisis: structural change a constant-coefficient VAR cannot see. The model nests the whole series — freeze the volatilities → Part 3's BVAR, freeze the coefficients → Part 4's BVAR-SV, identify the shocks → Part 2's structural responses. Cross-checked in R against a rolling-window OLS VAR — the same drifting transmission, only noisier.

View example →

Bayesian Vector Error Correction (VECM)

Part 6 of the VAR arc: cointegration — the case in between stationary and unit-root, where several series each wander like a random walk yet specific linear combinations are stationary, pinned by a long-run equilibrium. The level VAR is differenced into error-correction form, and the reduced-rank restriction Π = αβ′ splits the single levels term into cointegrating vectors (β) and adjustment speeds (α), sampled from scratch by a three-block conjugate Gibbs sampler (a new bvecm.py module). Built as a trilogy that walks the full spectrum — each notebook pre-testing the integration order first: the Treasury term structure (genuinely cointegrated, r = 2 — keep the levels, don't over-difference), the Fisher effect (the premise fails — inflation is already stationary, so no cointegration exists), and money demand (all I(1), yet r = 0 — a classic relation that broke down, so forecast in differences). The lesson: the rank is a question you ask the data, and the honest answer dictates the forecasting model. Cross-checked in R with urca.

View example →

Hierarchical Bayesian VAR — Prior Selection & Density Forecasts

Part 7 of the VAR arc: letting the data choose the shrinkage. Every earlier BVAR fixed the Minnesota tightness by hand; Giannone, Lenza & Primiceri (2015) put a prior on the hyperparameters and estimate them by marginal likelihood — the fully Bayesian answer to "how much should a VAR shrink?". A switch to the natural-conjugate dummy-observation prior makes the evidence closed-form, so the overall tightness (λ), the sum-of-coefficients (μ) and dummy-initial-observation (δ) hyperparameters are sampled by Metropolis and the coefficients follow analytically (a new bvarglp.py module). The data want looser shrinkage than the convention (λ ≈ 0.29 vs 0.2), and the same evidence gives Bayes factors for lag length and prior structure (the unit-root / cointegration components supported beyond decisively). The payoff is a proper density-forecast evaluation — log predictive scores and PIT calibration — where the hierarchical model wins (log score −3.85 vs −4.28 fixed, −4.54 flat): estimating the shrinkage with its uncertainty improves the whole predictive distribution. The measured capstone of the shrinkage / SV / time-varying thread. Cross-checked in R against the BVAR package.

View example →

The Oil-Market SVAR — Supply, Demand & Identification

Part 8 of the VAR arc — its first non-monetary structural application, and a study in how much identification drives the answer. For decades the oil price was treated as an exogenous supply shock; but a price is set by supply and demand, so Kilian (2009) disentangles three shocks inside a structural VAR — a supply disruption, an aggregate-demand expansion, and an oil-specific precautionary demand shock — adding the variance and historical decompositions to the series (new oilsvar.py module). The verdict that reshaped the field: supply explains only ~2% of the real oil price, and the 2003–2008 super-spike was a demand event (the China-led boom, not an OPEC squeeze). A companion notebook then stress-tests that conclusion: relax the recursive within-month zeros and the supply share can jump to ~31% — it comes back down only once the short-run supply elasticity is disciplined (Kilian–Murphy bounds, the Baumeister–Hamilton elasticity prior). The honest lesson: Kilian's demand story survives conditional on inelastic supply, and no scheme is assumption-free. Cross-checked in R (vars, rotation search).

View example →

Big-Data VARs — Factor Augmentation & Large-Scale Stochastic Volatility

Part 9 of the VAR arc: letting a VAR see the whole economy, two complementary ways. Real macro needs a system that is big and has moving volatilities — and there are two routes to big data, shrink or compress. 9a keeps all the variables and shrinks: a 10-variable BVAR with stochastic volatility (Carriero–Clark–Marcellino), made feasible by the triangular factorization, whose by-product is a common macro-uncertainty factor — one principal component explains 73% of the ten volatilities' comovement. 9b compresses instead: the Factor-Augmented VAR (Bernanke–Boivin–Eliasz) summarises a 169-series panel by a few principal-component factors, then recovers the impulse response of every underlying series to a monetary shock through its loadings — an economy-wide contraction from a single small VAR. Together they join the large-VAR and stochastic-volatility strands of the series. Cross-checked in R (vars, stochvol).

View example →

SVAR Identification through Heteroskedasticity

Part 10 of the VAR arc: identification without an economic assumption. Every structural VAR so far bought identification with a belief — a recursive zero, a sign restriction, an instrument, an elasticity bound. Rigobon (2003) needs none: if the structural shock variances change over time while the contemporaneous structure stays put, the change in volatility alone pins down the impact matrix — a purely statistical identification via a generalised eigenproblem (simultaneous diagonalisation of two regime covariances). Taking the cleanest cut, the Great Moderation, the three shock variances fell by three different proportions (to 0.11, 0.26, 0.44 of their pre-1984 levels), so the structure is identified with no restriction — and it rejects the recursive ordering: the contemporaneous inflation→GDP effect is 1.19 [0.38, 1.97], which a Cholesky scheme forces to zero. A fourth kind of identification that imposes nothing, so it can test what the others assume — and closes the loop with Part 4: the stochastic volatility modelled there as nuisance is exactly what identifies the shocks. Cross-checked in R.

View example →

Markov-Switching VAR

Part 11 of the VAR arc: discrete regimes, discovered by the data. Part 5 let the dynamics drift smoothly; the Markov-switching VAR (Hamilton; Krolzig) says the important variation is a jump between states — calm vs turbulent — governed by a latent Markov chain the data infer. Built from scratch with a Gibbs sampler whose centrepiece is the forward-filter–backward-sample of the hidden regime path (the Hamilton filter, then a backward draw — the discrete-state Carter–Kohn). Handed no dates, it rediscovers the Great Moderation: the turbulent regime is active 40% of the time before 1984 but only 3% after — exactly where Part 10 had to assume the break. The turbulent state carries 2–4× the shock volatility and spikes in every NBER recession; the calm regime lasts ~26 quarters, the turbulent one ~6. An endogenous break detector delivering a dated, regime-by-regime chronology of postwar macro — completing the trio (smooth drift, imposed regimes, discovered regimes) of how a VAR can be non-constant. Cross-checked in R (MSwM).

View example →

The Horseshoe BVAR — Global-Local Adaptive Shrinkage

Part 12 of the VAR arc: shrinkage that learns. Every VAR so far shrank with the Minnesota prior — a fixed rule that cannot tell a genuine dynamic from noise. The horseshoe (Carvalho–Polson–Scott 2010) is a global-local prior that does: one global scale pulls every coefficient toward zero, while heavy-tailed local scales let a few important ones escape, giving a U-shaped shrinkage profile — coefficients are either kept nearly intact or crushed to zero. Built from scratch with the Makalic–Schmidt (2016) auxiliary-variable Gibbs (every full conditional inverse-gamma or normal) and fit to a 10-variable VAR(2). Told nothing about the variables, it performs automatic variable selection: it kills ~96% of the 200 coefficients and keeps exactly the ones economics expects to persist — the own first lags of the interest rates and unemployment — while shrinking the near-unpredictable growth equations away. But an honest expanding-window test shows the sparsity is not free: the horseshoe forecasts significantly worse than Minnesota out of sample (p ≈ 0.01) — the aggressive shrinkage that buys interpretability over-shrinks genuine short-horizon signal. What it delivers is a radically smaller, readable model that names which dynamics are real, for a small forecast premium — the reason global-local priors are a modern default for large, high-dimensional VARs. Cross-checked in R.

View example →

The Fiscal Multiplier — a Blanchard–Perotti SVAR

Part 13 of the VAR arc: the fiscal multiplier, and a third structural domain after monetary and oil. If the government spends an extra dollar, how much does GDP rise? The obstacle is simultaneity, and Blanchard–Perotti (2002) solve it with an institutional identification rather than a statistical one: fiscal decision lags mean spending is predetermined within the quarter, and the automatic-stabiliser response of taxes is calibrated from the tax code, not estimated. Built on the reused Part 8 SVAR toolkit, it recovers the canonical results: the spending multiplier is ~0.9 (a near one-for-one, clean causal effect), while tax hikes contract output more slowly, building past −1 over two years. The most economics-driven identification in the series — where the elasticity is genuinely known from institutions rather than assumed. Cross-checked in R (vars).

View example →