VAR Models: Vector Autoregression

Before Sims (1980), macroeconomists built giant simultaneous-equations models and used “theory” to justify their identification restrictions. Sims pointed out that the justifications mostly did not hold up: the restrictions were untestable, the models fit poorly, and the policy conclusions depended on assumptions nobody could defend. His remedy was the unrestricted vector autoregression — regress every variable on the lags of every variable, and impose no “incredible restrictions” at all. The proposal was simple enough to force the entire field to rethink how it built models.

But what this article wants you to take away is not just how to derive a VAR. There is a thread running through all of it: “atheoretical” only moves the theory from the estimation stage to the interpretation stage. Estimating a VAR genuinely requires no economics; the moment you ask it anything interesting — how shocks propagate, who drives whom, what happens if policy changes — identification assumptions come back through the side door, and every one of them is a choice you make. See that clearly and you see what the entire spectrum from Cholesky orderings to sign restrictions is actually arguing about. The references are Hamilton (1994, Ch. 11), Lütkepohl (2005), and Tsay (2010, Ch. 4).

The Basic VAR

A $K$-variable VAR of order $p$:

where $\mathbf{y}_t = (y_{1t}, y_{2t}, \ldots, y_{Kt})'$ is a $K$-vector, $\mathbf{c}$ an intercept vector, the $\mathbf{A}_i$ are $K \times K$ coefficient matrices, and $\mathbf{u}_t = (u_{1t}, \ldots, u_{Kt})'$ is the innovation vector, satisfying:

Note that $\boldsymbol{\Sigma}$ is not required to be diagonal — contemporaneous innovations may be correlated, and that “may be correlated” is where half the story starts.

Stationarity requires all roots of $\det(\mathbf{I}_K - \mathbf{A}_1 z - \cdots - \mathbf{A}_p z^p) = 0$ to lie outside the unit circle. The companion form makes this transparent. Define the $Kp \times 1$ state vector:

Then:

where

Any VAR($p$) collapses into one big VAR(1), which is stationary if and only if all eigenvalues of $\mathbf{F}$ have modulus less than 1.

A stationary VAR can be inverted into its VMA($\infty$) representation:

with $\boldsymbol{\mu} = (\mathbf{I} - \mathbf{A}_1 - \cdots - \mathbf{A}_p)^{-1} \mathbf{c}$ and the $\boldsymbol{\Psi}_j$ obeying the recursion $\boldsymbol{\Psi}_j = \mathbf{A}_1 \boldsymbol{\Psi}_{j-1} + \cdots + \mathbf{A}_p \boldsymbol{\Psi}_{j-p}$, with $\boldsymbol{\Psi}_0 = \mathbf{I}_K$ and $\boldsymbol{\Psi}_j = \mathbf{0}$ for $j < 0$. These moving-average matrices record how a shock from $j$ periods ago propagates to today — hold on to them, because impulse responses are built entirely out of them.

Estimation: Suspiciously Cheap

Estimation is the least painful part of a VAR. Lay out the $K$ equations and run OLS on each; moreover, equation-by-equation OLS is equivalent to system GLS — every equation has the same regressors (lagged $\mathbf{y}$ plus an intercept), so cross-equation error correlation costs no efficiency. For equation $i$, regress $y_{it}$ on a constant and $\mathbf{y}_{t-1}, \ldots, \mathbf{y}_{t-p}$, collecting $\hat{\mathbf{A}}_j$ and residuals $\hat{\mathbf{u}}_t$; the covariance matrix is the sample covariance of the residuals with a degrees-of-freedom correction ($Kp + 1$ regressors per equation):

Under standard regularity conditions the estimator is consistent and asymptotically normal (see the results around Hamilton 1994, Theorem 11.4). The stacked coefficient vector $\hat{\boldsymbol{\beta}} = \text{vec}([\hat{\mathbf{c}}, \hat{\mathbf{A}}_1, \ldots, \hat{\mathbf{A}}_p])$ has asymptotic covariance:

where $\boldsymbol{\Gamma} = E[\mathbf{x}_t \mathbf{x}_t']$ is the second-moment matrix of the regressors ($\mathbf{x}_t$ stacking the intercept and the lags of $\mathbf{y}$), built from the autocovariances of $\mathbf{y}$ at the various lags. In practice the software does all of this — statsmodels.tsa.api.VAR in Python, the vars package in R.

The cheapness collects its fee in parameters. Each equation has $Kp + 1$ of them, the system $K^2 p + K$ in total, plus $K(K+1)/2$ in the covariance matrix. With $K = 5$ and $p = 4$ you are already holding a hundred-plus coefficients, and degrees of freedom get tight immediately. So lag selection is not a formality: AIC, BIC, and HQ all work, with BIC leaning toward smaller $p$ than AIC; Sims (1980) advised generosity with lags to avoid dynamic misspecification, and Lütkepohl (2005, Ch. 4) surveys the asymptotics. The more radical cure for the parameter curse is Bayesian shrinkage, treated below.

Granger Causality: The Easiest Test in the VAR, the Most Misread Word

Variable $y_{1t}$ Granger-causes $y_{2t}$ if, beyond the own lags of $y_{2t}$, the lags of $y_{1t}$ add incremental predictive value for $y_{2t}$. Inside a VAR this is a plain coefficient test. Write out the equation for $y_{2t}$:

The null is that the lags of $y_{1t}$ do not enter: $H_0: \beta_{2,1,1} = \cdots = \beta_{2,1,p} = 0$, tested by Wald or F. Equivalently, compare the likelihoods of the unrestricted and restricted models:

Reject when $LR$ is large. Run it for every ordered pair and you have the full Granger causality matrix. A separate question is instantaneous causality: after controlling for the dynamics, is the contemporaneous correlation between $u_{1t}$ and $u_{2t}$ zero — look at the off-diagonal elements of $\hat{\boldsymbol{\Sigma}}$; under the null the sample covariance is asymptotically normal with variance of order $1/T$.

Now the word itself. “Granger-cause” is not “cause” in any philosophical sense; its entire content is “carries incremental predictive information”. A variable can Granger-cause another without being its cause (both driven by a third variable), and a variable can be a cause yet fail to Granger-cause (the effect is instantaneous and gets absorbed into contemporaneous correlation). This distinction is not pedantry — it decides what a test result can be used to claim. The standard empirical questions: does money supply Granger-cause inflation (in macro data, usually yes); does the yield curve Granger-cause recessions (the predictive power of the term spread is a stable finding in the literature); oil prices for stock returns (contested); volume for volatility (a common positive finding). Every one of these is a statement about prediction, not mechanism.

The test itself bites in a few places. Omit a common cause and two variables will appear to Granger-cause each other. With I(1) variables that are not cointegrated, the test’s distribution is non-standard and prone to over-rejection — back to the stationarity foundations. Relationships also drift across samples: what held in 1960-1990 may fail in 1990-2020. And the standard test is linear; nonlinear predictability slips right past it.

Impulse Responses: Where Identification Enters

The IRF answers: hit the system with a shock today, and how do all the variables move over the next $h$ periods? The answer lives in the VMA matrices. Shock the innovations at time 0 by $\boldsymbol{\delta}$, holding all other periods at zero:

The $h$-period response to the $i$-th innovation is the $i$-th column of $\boldsymbol{\Psi}_h$. So far, no economics needed.

The problem arrives immediately: $\boldsymbol{\Sigma}$ has off-diagonal elements, so innovations are contemporaneously correlated across equations. The thought experiment “shock only $u_{1t}$ and hold the rest at zero” does not exist in the data — they typically move together. To pry the shocks apart, first make them orthogonal. The standard device is the Cholesky decomposition:

with $\mathbf{P}$ lower triangular and positive on the diagonal. Define the orthogonalized shocks $\mathbf{e}_t = \mathbf{P}^{-1} \mathbf{u}_t$, mutually uncorrelated with unit variance:

The IRF to the $j$-th orthogonal shock is the $j$-th column of $\boldsymbol{\Psi}_h \mathbf{P}$.

Notice what just happened. Cholesky is not a neutral piece of linear algebra — it imposes a recursive structure: shocks to variable 1 affect all variables contemporaneously, shocks to variable 2 affect variables 2 through $K$ but not 1, and so on. That is a causal ordering: variable 1 is the most contemporaneously “exogenous”, variable $K$ the most endogenous. This is the thread from the introduction cashing out — the “atheoretical” VAR needs theory again the instant you plot your first orthogonalized IRF, and the theory hides in the ordering. Convention orders by speed of reaction: slow variables (real GDP, capital stock) first, fast ones (interest rates, monetary aggregates) last. But convention is not truth: put monetary policy before output and after output, and you get different answers. So the first thing to do with an IRF is re-run it under a few orderings. If the results are qualitatively stable, the identification has some spine; if they flip, you are reporting the ordering, not the data.

Forecast Error Variance Decomposition

The FEVD asks for the long-run ledger: of variable $i$‘s $h$-step forecast error variance, how much belongs to each orthogonal shock? The $h$-step forecast error is:

The MSE matrix is:

Substituting $\boldsymbol{\Sigma} = \mathbf{P}\mathbf{P}'$ and rewriting in orthogonal shocks:

The share of shock $j$ in variable $i$‘s MSE at horizon $h$:

which sums to 1 over $j$. Interpretation: $\omega_{1,2}(24) = 0.6$ means 60% of variable 1’s 24-step forecast error variance is attributable to shocks in variable 2 — the technical foundation under claims like “variable 2 drives variable 1 in the long run”. Naturally, it inherits every identification dependency of the IRF: change the ordering and the decomposition changes faces.

Structural VARs: Identification on the Table

Cholesky is only the entry point of the identification spectrum. Formalize the problem: a structural VAR posits economically meaningful structural shocks $\boldsymbol{\varepsilon}_t$, with structural form:

where $E[\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_t'] = \mathbf{I}_K$ (orthogonal, unit variance), $\mathbf{A}_0$ captures contemporaneous interactions, and $\mathbf{B}$ maps structural shocks to reduced-form errors. Multiply through by $\mathbf{A}_0^{-1}$ to recover the reduced form:

with $\mathbf{A}_i = \mathbf{A}_0^{-1} \mathbf{A}_i^*$, $\mathbf{u}_t = \mathbf{A}_0^{-1} \mathbf{B} \boldsymbol{\varepsilon}_t$, and hence $\boldsymbol{\Sigma} = \mathbf{A}_0^{-1} \mathbf{B} \mathbf{B}' (\mathbf{A}_0^{-1})'$. The identification arithmetic: the mapping to be pinned down has $K^2$ free elements, while $\boldsymbol{\Sigma}$ supplies only $K(K+1)/2$ identifiable quantities, leaving a gap of $K^2 - K(K+1)/2 = K(K-1)/2$ restrictions — restrictions the data cannot provide and only theory can. Where you go to get them separates the three main schools.

Recursive identification (Sims, 1980): $\mathbf{A}_0$ lower triangular with ones on the diagonal, $\mathbf{B}$ diagonal. This is exactly Cholesky — the triangular structure supplies precisely $K(K-1)/2$ restrictions, and the ordering problem of the previous section is the choice of that triangle.

Long-run restrictions (Blanchard & Quah, 1989): constrain nothing contemporaneously; constrain the long run instead. The identifying assumption is that demand shocks have no long-run effect on output, while supply shocks can have permanent effects — the restrictions land on the long-run multiplier $\sum_{j=0}^\infty \boldsymbol{\Psi}_j \mathbf{A}_0^{-1} \mathbf{B}$. Suited to settings where theory speaks clearly about what is transitory and what is permanent.

Sign restrictions (Uhlig, 2005): the weakest tier. Require only that certain shocks move certain variables in certain directions over certain horizons — a contractionary monetary shock should not raise prices, say. The cost is retreating from point identification to set identification: what is consistent with the restrictions is a region of shocks, from which Bayesian methods sample. In current empirical macro this is standard equipment.

None of the three schools rides free: the recursive ordering, the long-run theory, the sign priors — each is a choice you make on the data’s behalf. What the VAR did was move those choices from the shadows of the simultaneous-equations era into the open. That, arguably, is the real legacy of the Sims revolution.

Unit Roots in Levels: Go to the VECM

Everything above assumed stationarity. If the levels $\mathbf{y}_t$ are non-stationary but some linear combination $\boldsymbol{\beta}' \mathbf{y}_t$ is stationary, then $\mathbf{y}_t$ is cointegrated, a VAR in levels is misspecified, and the finite-sample standard errors are distorted. Differencing everything rescues stationarity but throws away the long-run information; the right framework is the vector error correction model:

with $\boldsymbol{\beta}' \mathbf{y}_{t-1}$ the cointegrating relation and $\boldsymbol{\alpha}$ the adjustment speeds. Cointegration earns two full articles of its own, immediately after this one.

Bayesian VARs: The Cure for the Parameter Curse

Back to the mine laid in the estimation section: $K^2 p + K$ parameters. With $K = 10$ and $p = 4$ that is 410, against macro samples of $T \approx 200$ — unrestricted estimation simply cannot stand. The Bayesian VAR of Litterman (1986) and Sims & Zha (1998) shrinks the coefficients toward sensible defaults with a prior. The workhorse is the Minnesota prior, whose default belief is “each variable is approximately its own random walk”: the own first lag gets a prior mean of 1, everything else a prior mean of 0; prior variances shrink as the lag order grows (distant lags are held closer to zero a priori), and cross-variable lags are additionally scaled by the ratio of the two variables’ residual variances, so the prior’s tightness is comparable across variables of different scales. The posterior mean then becomes a weighted average of the OLS estimate and the prior mean, with weights set by the prior’s precision relative to the data. The payoff is concrete: BVARs forecast well out of sample, especially when $T$ is not rich relative to the parameter count — they sidestep the unrestricted VAR’s overfitting.

Forecasting and Putting It to Work

The $h$-step VAR forecast is computed recursively:

For $h \leq p$ the right-hand side mixes observed values and already-computed forecasts; beyond $h > p$ it is forecasts all the way down. The forecast error:

with MSE matrix $\text{MSE}(h) = \sum_{j=0}^{h-1} \boldsymbol{\Psi}_j \boldsymbol{\Sigma} \boldsymbol{\Psi}_j'$, and Gaussian forecast intervals:

The code side is entirely off the shelf:

import pandas as pd
from statsmodels.tsa.api import VAR


model = VAR(y)
lag_order = model.select_order(maxlags=8).aic
print("Selected lag order:", lag_order)

results = model.fit(lag_order)
print(results.summary())

# Granger causality
print(results.test_causality('gdp', caused=['inflation', 'interest_rate'], kind='f'))

# IRF and FEVD
irf = results.irf(periods=20)
irf.plot(orth=True)

fevd = results.fevd(periods=20)
fevd.plot()

Where it gets used. Macro forecasting: BVARs are the standard for nowcasting and short-horizon forecasts of GDP, inflation, and unemployment, and central-bank research departments use them for policy analysis. Monetary policy analysis: estimate a VAR on output, prices, interest rates, and money; identify a contractionary shock (via sign restrictions, say) and trace the response of output and prices — the monetary policy shock is the cornerstone of modern empirical monetary economics. Fiscal policy works the same way, treating government spending and taxes as shocks; multiplier estimates (how much does output respond to another dollar of spending?) come out of fiscal SVARs. On the finance side: VARs on yield-curve factors (level, slope, curvature) for term-structure dynamics, VARs on Fama-French factors for risk premia, VARs on returns and macro variables for return predictability. Mixed-frequency VARs (MF-VARs) use weekly and monthly data to update estimates of the current quarter in real time — that is the nowcasting machine.

Wherever you use it, run this article’s thread once more before signing off on a conclusion: you used no theory at the estimation stage — which theory did you use at the interpretation stage? Who chose the lag order, did you try other orderings, are there unit roots in the levels, and does the word “cause” mean prediction or mechanism? The honesty of the VAR is that each of those questions has a definite place where the answer can be checked — provided you actually go check.

References

• Blanchard, O. J., & Quah, D. (1989). The dynamic effects of aggregate demand and aggregate supply disturbances. American Economic Review, 79(4), 655-673.
• Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press. Ch. 11.
• Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer.
• Sims, C. A. (1980). Macroeconomics and reality. Econometrica, 48(1), 1-48.
• Sims, C. A., & Zha, T. (1998). Bayesian methods for dynamic multivariate models. International Economic Review, 39(4), 949-968.
• Stock, J. H., & Watson, M. W. (2001). Vector autoregressions. Journal of Economic Perspectives, 15(4), 101-115.
• Tsay, R. S. (2010). Analysis of Financial Time Series (3rd ed.). Wiley. Ch. 4.
• Uhlig, H. (2005). What are the effects of monetary policy on output? Results from an agnostic identification procedure. Journal of Monetary Economics, 52(2), 381-419.
• Litterman, R. B. (1986). Forecasting with Bayesian vector autoregressions. Five years of experience. Journal of Business & Economic Statistics, 4(1), 25-38.