Classical structural models imposed a priori exclusion restrictions that were “incredible.” Sims proposed the unrestricted VAR as an atheoretical alternative, letting the data speak. All variables are treated symmetrically — there is no arbitrary distinction between endogenous and exogenous.
Tip
MA(\(\infty\)) representation
For a stable VAR(\(p\)), the Wold representation is: \[\mathbf{y}_t = \boldsymbol{\mu} + \sum_{h=0}^{\infty}\boldsymbol{\Phi}_h\,\mathbf{u}_{t-h}, \qquad \boldsymbol{\Phi}_0 = \mathbf{I}_K, \quad \boldsymbol{\Phi}_h = \sum_{j=1}^{p}\mathbf{A}_j\boldsymbol{\Phi}_{h-j}\]
The matrices \(\boldsymbol{\Phi}_h\) are the impulse-response coefficients at horizon \(h\).
Criterion
Formula
Penalty
AIC
\(\ln|\hat\Sigma_u| + \frac{2}{T}K^2 p\)
Lightest
HQ
\(\ln|\hat\Sigma_u| + \frac{2\ln\ln T}{T}K^2 p\)
Medium
BIC / SC
\(\ln|\hat\Sigma_u| + \frac{\ln T}{T}K^2 p\)
Heaviest
Rules of thumb: AIC tends to over-fit; BIC is consistent; HQ intermediate. Always confirm lag adequacy via residual autocorrelation (Portmanteau test).
\(\omega_{jk,h}\): fraction of the \(h\)-step forecast error variance of variable \(j\) explained by shocks to variable \(k\).
Tip
Reading a FEVD table
If \(\omega_{jj,h} \approx 1\) at all horizons, variable \(j\) is nearly exogenous — its variance is dominated by own shocks. If \(\omega_{jk,h}\) grows with \(h\) for \(k \ne j\), the inter-variable linkage matters more at longer horizons (typical in monetary policy settings).
Variable \(x\)Granger-causes\(y\) if past values of \(x\) improve forecasts of \(y\) beyond what past \(y\) alone provides.
Test in VAR:\(H_0\): \(A^{xy}_{1} = A^{xy}_{2} = \cdots = A^{xy}_{p} = \mathbf{0}\) — the block of coefficients of \(x\) in the \(y\)-equation are jointly zero.
Implemented via Wald test: causality(var_fit, cause = "x") in R; vargranger in Stata.
Warning
Granger causality is a predictive concept, not structural causality. Rejection means \(x\) contains incremental predictive information about \(y\) — it does not imply \(x\) causes \(y\) in a counterfactual sense.
library(vars)data(Canada)# 1. Lag selectionvs <-VARselect(Canada, lag.max =8, type ="const")vs$selection # AIC, HQ, SC, FPE selections
AIC(n) HQ(n) SC(n) FPE(n)
3 2 1 3
# 2. Estimate VAR(2) — chosen by AICvar_fit <-VAR(Canada, p =2, type ="const")# 3. Coefficient table for the employment (e) equationround(coef(var_fit)$e, 4)
Solid line: point estimate of the IRF (the OIRF or generalised IRF)
Shaded band / dashed lines: bootstrap 95% confidence band
Horizontal dashed red line at zero: if the band crosses zero, the response is not statistically significant at 5%
Horizon: measured in the model’s time unit (quarters here)
Magnitude: in units of the response variable per one standard deviation shock
The Cholesky ordering matters: here, productivity (prod) is ordered first, so its shock is contemporaneously uncorrelated with shocks to the other three variables by construction.
Identification problem: OLS gives \(\hat\Sigma_u = \mathbf{A}_0^{-1}\mathbf{A}_0^{-\prime}\). This gives \(\frac{K(K+1)}{2}\) equations but \(K^2\) unknowns in \(\mathbf{A}_0\). We need \(\frac{K(K-1)}{2}\)identifying restrictions.
Blanchard & Quah (1989) identify supply vs demand shocks by imposing that demand shocks have no long-run effect on output (but can affect unemployment temporarily).
Restriction: the (1,2) element of the cumulated IRF matrix ((1) = _{h=0}^{}_h) is constrained to zero — the long-run multiplier of the demand shock on output is zero.
In R: BQ(var_fit) implements Blanchard-Quah decomposition automatically for a bivariate VAR.
Key reference:Blanchard, O.J. & Quah, D. (1989). “The Dynamic Effects of Aggregate Demand and Supply Disturbances.” American Economic Review, 79(4), 655–673.
Sign Restrictions (Uhlig 2005)
Instead of zero restrictions, Uhlig proposes restrictions on the sign of IRFs at certain horizons. Example for a contractionary monetary policy shock:
Interest rate rises (response > 0)
Money supply falls (response < 0)
Output falls (response < 0)
No restriction on the price level
Implemented via rejection sampling over the rotation matrix (): draw random orthogonal (), accept if all sign conditions are met, else reject. In R: uhlig.reject() from the VARsignR package.
Advantage: avoids controversial zero restrictions; relies on theory-consistent qualitative features only.
Cholesky Ordering Matters
The Cholesky decomposition imposes an asymmetric contemporaneous structure: the first variable is affected only by its own shock at (t); the last variable can be contemporaneously affected by all others.
For the Canada data, a defensible ordering based on information availability and sluggishness of adjustment:
prod → e → rw → U
Productivity is the most exogenous (determined by technology), employment responds next, then real wages adjust, and unemployment is the most endogenous sluggish outcome. Always motivate your ordering from theory.
Problem: many macro series move together but VAR(\(p\)) has a curse of dimensionality. Running a VAR on all \(N\) series is infeasible for large \(N\).
Solution (Stock & Watson 2005; Bernanke, Boivin & Eliasz 2005): extract a small number of common factors \(\mathbf{F}_t\) from a large panel \(\mathbf{X}_t\), then build a VAR on \((\mathbf{F}_t, \mathbf{Y}_t)\):
Factors estimated by PCA in Step 1 (two-step estimator)
Two-step procedure:
Factor extraction: Apply PCA to (standardised) \(\mathbf{X}_t\) → first \(r\) principal components \(\hat\mathbf{F}_t\)
FAVAR estimation: Run VAR on \((\hat\mathbf{F}_t, \mathbf{Y}_t)\) by OLS equation-by-equation
Note
Choosing \(r\) (number of factors)
Use the information criteria of Bai & Ng (2002): ICp1, ICp2, BIC3 — all implemented in the FactoMineR and factoextra R packages. A scree plot of eigenvalues provides a graphical guide: look for the “elbow” where the slope flattens.
Tip
One-step (Bayesian) alternative
Bernanke, Boivin & Eliasz (2005) propose a one-step Gibbs sampler that jointly estimates factors and VAR parameters, yielding proper posterior uncertainty over both. Implemented in the BFARpack R package and in MATLAB.
A VAR(\(p\)) with \(K\) variables has \(K^2 p + K\) free parameters.
\(K\) (variables)
\(p = 2\) lags
\(p = 4\) lags
4
36
68
10
210
410
20
820
1620
100
20 100
40 100
With \(T = 200\) quarterly obs, a VAR(4) with \(K = 10\) already uses 410 parameters — more than 2 parameters per observation.
Bayesian solution: Shrink coefficients toward a prior belief (e.g. each variable follows a random walk). The posterior combines prior + likelihood, effectively borrowing strength across equations.
\(\delta_i = 1\) for variables in levels (random walk prior), \(= 0\) for stationary variables
\(\lambda\) (overall tightness): small \(\lambda\) → strong shrinkage toward prior; large \(\lambda\) → dominated by data
\(l^2\) in denominator: higher lags are shrunk more strongly (coefficient decay)
\(\sigma_i/\sigma_j\): rescales for different units across variables
Modern variants (Giannone, Lenza & Primiceri 2015) treat \(\lambda\) as a hyperparameter and integrate it out — the data choose the shrinkage intensity.
Code — BVAR
MCMC Settings for BVAR
n_draw: total Gibbs sampler draws (recommend ≥ 10 000 for publication)
n_burn: warm-up draws discarded (usually 50% of n_draw)
n_thin: keep every \(k\)-th draw to reduce autocorrelation
Convergence diagnostics: Gelman-Rubin \(\hat{R}\) across multiple chains (values < 1.05 indicate convergence); trace plots for visual inspection
Parallel chains: set n_chains = N_CORES and use future::plan(multisession) for efficient computation
import numpy as npimport pandas as pdfrom statsmodels.datasets import macrodata# Manual Minnesota prior (Ridge-regularised VAR) as a lightweight approximationfrom numpy.linalg import solvedata = macrodata.load_pandas().dataY_b = pd.DataFrame({"dlgdp": np.log(data["realgdp"]).diff(),"dlcons": np.log(data["realcons"]).diff(),"dlinv": np.log(data["realinv"]).diff(),}).dropna().to_numpy()T_, K_ = Y_b.shape; p_ =2# Build regressor matrixZ = np.column_stack([np.ones(T_-p_)] + [Y_b[p_-j:T_-j] for j inrange(1, p_+1)])Y_dep = Y_b[p_:, :]# Tikhonov (L2) shrinkage — approximates Minnesota prior lam =0.2# overall tightnessI_ = np.eye(Z.shape[1])I_[0, 0] =0# do not shrink the constantB_bayes = solve(Z.T @ Z + lam * T_ * I_, Z.T @ Y_dep)print("Bayesian (Minnesota-like) coefficient matrix:\n", np.round(B_bayes, 4))
* Native Bayesian VAR is not available in Stata 19 without add-ons.* Use R (BVAR package) or Python (PyMC / bvar package via pip) instead.* For a quick frequentist shrinkage alternative in Stata:* ridge regression on the VAR stacked system gives a similar effect.
\(q_{t-d}\): threshold variable (observed at delay \(d\); can be one of the system variables or an external indicator)
\(\gamma\): threshold value — estimated by grid search over the sorted values of \(q_{t-d}\)
Regime-specific coefficient matrices \(\mathbf{A}_j^{(r)}\) allow fully different dynamics in each regime
Trim parameter (e.g. 10%): excludes extreme quantiles of \(q\) to ensure a minimum number of observations per regime
Economic motivation: monetary transmission may differ in recessions vs expansions; financial stress can trigger non-linear adjustment; credit constraints bind asymmetrically.
Testing linearity:Hansen (1996) sup-Wald test (bootstrap \(p\)-value required, since the threshold \(\gamma\) is unidentified under \(H_0\)).
Note
Regime assignment
The threshold \(\hat\gamma\) is a super-consistent estimator (Chan 1993) — it converges at rate \(T\), faster than the \(\sqrt{T}\) rate for the slope coefficients. This means in large samples, regime membership is determined nearly without error, and the subsequent inference on \(\mathbf{A}_j^{(r)}\) can proceed as if the regimes were known.
Warning
Persistent threshold variable
If \(q_{t-d}\) is highly persistent (e.g. a near-unit-root variable), the effective sample in each regime can be very small and the threshold estimate unreliable. Check stationarity of \(q\) and the balance of observations across regimes before reporting.
\(w_{ij}\): trade-share weights (bilateral import shares, typically averaged over 3 years)
Each country-specific VARX* estimated separately by OLS; global solution via link matrix
Global solution: stack all country VARXs → solve simultaneously for \(\mathbf{x}_t = (\mathbf{x}_{0t}',\ldots,\mathbf{x}_{Nt}')'\) — a global VAR with \(\sum_i K_i\) variables.
Tip
Key applications
Shock transmission in global supply chains; contagion of financial crises; commodity price pass-through; multi-country monetary policy spillovers. A companion dataset with 33 countries and 7 macro variables for 1979Q1–2019Q4 is available at mohaddes.com.
library(BGVAR)# eerDatasmall: 3 regions (US, Euro Area, Rest of World)# Variables: y (output), Dp (inflation), r (interest rate)data("eerDatasmall")# Trade-weighted cross-sectional weights (2000–2012 average)# W.trade0012 ships with BGVARgvar_fit <-bgvar(Data = eerDatasmall,W = W.trade0012,draws =500,burnin =500,thin =2,plag =1, # one lag per country VARX*prior ="MN", # Minnesota prioreigen =1.05, # stability cutoffverbose =FALSE)# IRF: US output shock → Euro Area outputus_shock <-list(var ="y", country ="US")irf_gvar <-IRF(gvar_fit, shock = us_shock, n.ahead =20,ident ="chol", verbose =FALSE)plot(irf_gvar, resp.var ="y",main ="US output shock → Global output (BGVAR)")
Note
Practical guidance
Estimating a GVAR requires:
A balanced panel of country-level macro variables (monthly or quarterly)
A matrix of bilateral trade (or financial) weights — must be row-stochastic
Country-specific lag selection; typically 1–2 lags per country VARX*
Global stability check: eigenvalues of the global companion matrix should lie inside the unit circle
The BGVAR package ships with the eerDatasmall dataset (3 regions, quarterly 1995Q1–2019Q4) for replication of the main US/EA linkage results. For a full 33-country dataset see Mohaddes & Raissi (2024).
\(\boldsymbol{\alpha}_i\): unit-specific fixed effects (allow for unobserved heterogeneity)
Slope matrices \(\mathbf{A}_j\) are homogeneous across units (common dynamics assumption)
\(\mathbf{u}_{it}\): idiosyncratic errors; \(\Sigma_u\) may differ across units
Estimation challenge: FE estimation is inconsistent for short \(T\) (Nickell bias). Solution: first-difference to remove \(\boldsymbol{\alpha}_i\), then use lagged levels as GMM instruments (Arellano-Bond / Blundell-Bond).
Panel IRFs: compute as for a standard VAR using the pooled coefficient estimates. Bootstrap uncertainty by block-resampling panel units.
Dahlberg, Mörk, Rattsø & Ågren (2008) — Swedish local government finance:
Variable
Description
expenditures
Per-capita total expenditure (SEK)
revenues
Per-capita total revenues (SEK)
grants
Per-capita central government grants (SEK)
290 Swedish municipalities, 1979–1987. Classic dataset for testing flypaper effect: do grants raise public spending one-for-one, or is there fiscal substitution?
library(panelvar)# Dahlberg et al. (2008) — Swedish local government financedata("ex1_dahlberg_data")# Panel VAR(1) via Arellano-Bond first-difference GMMpvar_fit <-pvars( ex1_dahlberg_data[, c("expenditures", "revenues", "grants")],lags =1,transformation ="fd", # first-difference to remove FEsexo =TRUE,panel.id = ex1_dahlberg_data$id,panel.T = ex1_dahlberg_data$year)summary(pvar_fit)# Generalised Impulse Response Functions (Pesaran & Shin 1998)set.seed(14159)girf_fit <-girf(pvar_fit, n.ahead =8, ma.approx.n =20)plot(girf_fit, main ="Panel GIRF — Swedish municipalities")
[panelvar::pvars] not available — please install panelvar from CRAN
Code
* Panel VAR is not a built-in Stata 19 command.* The pvar user-written package (Abrigo & Love, 2016) must be installed:* ssc install pvar* Once installed:* webuse dahlberg, clear (Stata's own Dahlberg data)* xtset id year* pvar expenditures revenues grants, lags(1) gmmopts(twostep)* pvarirf, mc(200) step(8) impulse responses* pvarfevd, mc(200) step(8) variance decomposition* For native Stata, first-difference GMM is available via:xtabond2 expenditures L.expenditures L.revenues L.grants, /// gmm(L.expenditures, lag(1 3)) iv(L.revenues L.grants) ///tworobust small
where \(\text{FCI}_t\) is a financial conditions index. The distribution of \(h\)-step-ahead growth is skewed left during financial stress — the downside tail thickens.
Important
Why quantile responses differ from mean responses
OLS-VAR only characterises the conditional mean response. QVAR reveals:
Asymmetric shocks: a positive productivity shock raises the median of employment, but the 5th percentile may barely move (downside protection)
Crisis regimes hidden in the mean: the mean IRF for GDP growth in a financial-crisis period averages across normal-expansion paths and rare-crash paths, masking the left-tail risk
Tail spillovers: a shock can be neutral at the median (\(\beta^{(0.5)} \approx 0\)) yet strongly affect the 5th percentile (\(\beta^{(0.05)} \ll 0\)) — detected only by QVAR
Key idea: instead of iterating a VAR forward to compute the \(h\)-step-ahead response, regress the actual future value \(y_{t+h}\) directly on current and lagged information:
Ramey (2016, Handbook of Macroeconomics) shows that LP-IRFs and VAR-IRFs agree asymptotically when the VAR is correctly specified and correctly lag-augmented. Differences in finite samples arise from the LP’s higher variance (no cross-horizon efficiency gains) vs the VAR’s bias when the true model is not finite-order. Nakamura & Steinsson (2018) advocate LP-IV for externally identified shocks.
Code — Local Projections
HAC Lag Selection in LP
At horizon \(h\), the LP residual \(\varepsilon_{t+h}\) has an MA(\(h-1\)) structure by construction. Use at least \(h\) Newey-West lags for valid SE. In R lpirfs, this is handled automatically. In Stata: newey ..., lag(h). In Python: cov_type="HAC", cov_kwds={"maxlags": h}. Using too few HAC lags produces under-sized confidence bands — a common error in applied work.
Need structural shocks? → SVAR (Cholesky, BQ, sign)
Many indicators, small model? → FAVAR
Large \(K\), small \(T\)? → BVAR or Sparse VAR
Non-linear / regime change? → TVAR
International spillovers? → GVAR
Firms, regions, or countries? → PVAR
Tail risk or GaR? → QVAR
Robust IRFs, external shocks? → LP / LP-IV
High-dimensional, sparsity suspected? → Sparse VAR
Common pitfalls:
Problem
Symptom
Fix
Too many lags
Overfit, noisy IRF
Use BIC
Wrong ordering
Cholesky artefacts
Use BQ or sign restrictions
Ignoring non-stationarity
Spurious IRFs
Difference or VECM
No HAC in LP
Under-sized CI
newey, lag(h)
BVAR with flat prior
Same as OLS
Tune \(\lambda\) by ML or hierarchical
GVAR stability
Explosive IRFs
Check companion eigenvalues
Tip
Golden rule
Always compare VAR-IRF and LP-IRF. If they agree, the parametric VAR structure is likely correct. If they diverge at long horizons, the LP is more robust but less efficient — prefer it for policy inference.
Further Reading
Textbooks
Lütkepohl (2005) — New Introduction to Multiple Time Series Analysis. Springer. The standard reference for VAR, SVAR, VECM theory.
Hamilton (1994) — Time Series Analysis. Princeton. Chapters 10–11: VAR; Chapter 20: state-space models.
Nakamura & Steinsson (2018). “Identification in Macroeconomics.” JEP 32(3), 59–86. doi:10.1257/jep.32.3.59
Nicholson, Matteson & Bien (2017). “VARX-L: Structured Regularization for Large Vector Autoregressions with Exogenous Variables.” International Journal of Forecasting 33(3), 627–651. doi:10.1016/j.ijforecast.2017.01.003
Journal of Applied Econometrics · Journal of Econometrics · Review of Economics and Statistics · Journal of Monetary Economics · Economic Modelling
Exercises
Lag sensitivity. Using the Canada dataset, estimate VAR(1), VAR(2), VAR(4). Compare the IRF of prod → e at horizons 1, 4, 8 and 20 across the three models. At which horizon does lag order matter most?
Ordering sensitivity. In the Canada SVAR, swap the Cholesky ordering to U → rw → e → prod. How do the structural IRFs change? Interpret the result in terms of what contemporaneous exogeneity you are imposing.
Blanchard-Quah. Apply BQ(var_fit) to the Canada data. Identify the permanent and transitory shocks. Which variable absorbs the permanent shock? Does this match economic intuition?
LP vs VAR. For the Canada data, plot the LP-IRF and the VAR-IRF (both prod → e) on the same graph with their 95% CI. At which horizons do the bands diverge? What does divergence at long horizons suggest about the VAR model?
BVAR hyperparameter. Re-estimate the BVAR with \(\lambda \in \{0.01, 0.1, 0.2, 1.0\}\) (overall tightness). Plot the posterior mean IRFs for prod → e. How much does the choice of \(\lambda\) affect inference at short vs long horizons?
TVAR regime asymmetry. In the TVAR(Canada), extract regime-1 and regime-2 coefficient matrices. Compare the impulse response of prod → U in each regime. Is the productivity shock more recessionary in the high-unemployment regime?
GaR replication. Using US macro data (statsmodels.macrodata), regress 4-quarter-ahead GDP growth on a financial conditions proxy (e.g. interest rate spread) at \(\tau \in \{0.05, 0.25, 0.5, 0.75, 0.95\}\). Plot \(\hat\beta(\tau)\) and interpret the quantile heterogeneity.
Sparse VAR selection. Using BigVAR with struct = "SparseLag" on Canada, identify which cross-variable links survive penalisation at the optimal \(\lambda\). Do the surviving links match those identified by Granger causality tests in the full VAR?
Thank You
Athanassios Stavrakoudis Applied Informatics and Computational Economics Lab Department of Economics University of Ioannina, Greece
