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.
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.
Identification — what is a shock?
2 Structural VARrecursive, sign, instrument, long-run — four schemes on one posterior 8 Oil-market SVARand 8b, which stress-tests the assumption the whole result rests on 10 Via heteroskedasticityimposes nothing — and so can test what the others assume 13 Fiscal multiplieridentification from institutional timing and a calibrated elasticityScale & shrinkage — how many variables?
3 Large BVAR1,620 coefficients on 240 observations — shrink harder as it grows 7 Hierarchical priorstop guessing the tightness; estimate it by marginal likelihood 9 Factor augmentationthe other route to big data — compress 169 series instead of shrinking 12 Horseshoea prior that learns which coefficients matter — and costs forecast accuracyTwo 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.
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).
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).
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 prior —
tightened 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.
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.
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.
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.
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.
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).
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).
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).
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).