Vector Autoregression

Vector Autoregression (VAR)

A vector autoregressive process of order $p$ generalizes the AR(p) model from a single series to $k$ series whose components depend on each other through time. The state vector ${\mathbf{x}}_t=(x_{t,1},\dots ,x_{t,k}{)}^{\top }$ evolves as
${\mathbf{x}}_t=\mathbf{\alpha }+{\mathbf{\Phi }}_1{\mathbf{x}}_{t-1}+\cdots +{\mathbf{\Phi }}_p{\mathbf{x}}_{t-p}+{\mathbf{w}}_t,{\mathbf{w}}_t\sim {\mathrm{WN}}_k(0,\mathbf{Q}),$
where each ${\mathbf{\Phi }}_j$ is a $k\times k$ transition matrix that says how much the lag-$j$ value of every component predicts every other component, and $\mathbf{Q}$ is the contemporaneous covariance of the vector white-noise innovations. The simplest case VAR(1) reads ${\mathbf{x}}_t=\mathbf{\alpha }+\mathbf{\Phi }{\mathbf{x}}_{t-1}+{\mathbf{w}}_t$.

The off-diagonal entries of ${\mathbf{\Phi }}_j$ are what make VAR more than a stack of independent AR fits. The entry $({\mathbf{\Phi }}_j{)}_{i,m}$ measures how much the lag-$j$ value of series $m$ predicts today’s value of series $i$, the cross-influences that a univariate AR cannot capture. Many real systems require this: interest rates and inflation move together; blood markers in a patient affect each other day-to-day; neighboring brain regions co-fluctuate.

VAR is stationary when all roots of $\det (\mathbf{I}-{\mathbf{\Phi }}_{1}z-\cdots -{\mathbf{\Phi }}_p{z}^p)=0$ lie strictly outside the unit circle in the complex plane. The condition reduces to “all eigenvalues of the companion matrix have modulus less than one,” generalizing the univariate root condition. The standard estimator is equation-by-equation OLS: each row of the system is a separate regression of one component on all lagged components, and OLS on each row produces the maximum-likelihood estimate of the corresponding row of ${\mathbf{\Phi }}_1,\dots ,{\mathbf{\Phi }}_p$.

Coefficient Count Grows Quadratically $p$) on $k$ series has $k^{2}p+k$ coefficients in $\{{\mathbf{\Phi }}_j,\mathbf{\alpha }\}$ plus $k(k+1)/2$ in $\mathbf{Q}$. With $k=20$ series and $p=4$ lags that is over $1,600$ parameters before the noise covariance, far too many to estimate without regularization or sparsity assumptions. Classical unregularized VAR is practical only for $k$ up to around $5$ or $6$; modern high-dimensional VAR uses lasso or group lasso to identify which cross-influences are non-zero.

A VAR(

Intuition

A VAR is a coupled set of AR regressions, one per series, in which each series’s regressors include not only its own past but also the past of every other series in the system. The diagonal entries of the transition matrices ${\mathbf{\Phi }}_j$ govern how much each series predicts itself at lag $j$, the same self-feedback that a univariate AR captures. The off-diagonal entries govern how much each series predicts the others, the new information that VAR captures and that motivates fitting it in the first place. Reading the fitted ${\mathbf{\Phi }}_1$ matrix entry-by-entry tells you the system’s interaction structure: who drives whom, on what timescale, and how strongly.

The fitting story is unusually clean for a multivariate model. Each row of the VAR equations is a separate regression problem with the same regressors (lagged copies of all series), so OLS row-by-row is consistent and equivalent to the joint maximum-likelihood estimate. The dependence between equations shows up only in the residual covariance $\mathbf{Q}$, which is just the sample covariance of the residuals across the $k$ equations. The hard work is interpreting the fit. A $k\times k$ transition matrix can be enormous, and the meaningful summaries (impulse responses, variance decompositions, Granger causality tests) require additional steps beyond the raw coefficients.

Properties

(Companion Form Reduces VAR(p) to VAR(1)) Stack the lags into a $kp$-dimensional state.

Define ${\tilde{\mathbf{x}}}_t=({\mathbf{x}}_t^{\top },{\mathbf{x}}_{t-1}^{\top },\dots ,{\mathbf{x}}_{t-p+1}^{\top }{)}^{\top }$, a vector of length $kp$ that holds the most recent $p$ states stacked. The VAR(p) equations rewrite as ${\tilde{\mathbf{x}}}_t=\tilde{\mathbf{\Phi }}{\tilde{\mathbf{x}}}_{t-1}+{\tilde{\mathbf{w}}}_t$ where $\tilde{\mathbf{\Phi }}$ is the companion matrix with ${\mathbf{\Phi }}_1,\dots ,{\mathbf{\Phi }}_p$ in the top block-row and identity blocks below. Stationarity of the VAR(p) is equivalent to the spectral radius of $\tilde{\mathbf{\Phi }}$ being less than one. This is also the form used internally by software (statsmodels.tsa.vector_ar) and by the Kalman-filter machinery in DLMs when treating the VAR as a state-space model.

(OLS Equation-by-Equation) Each row of the system is a separate regression with the same design matrix.

Stacking lagged copies of all $k$ series gives a single design matrix that the same regression is run against $k$ times, once per output series. Provided the design has more rows than columns ($n\ge kp+1$, with substantial cushion for stable estimation), OLS on each row produces unbiased coefficients and the joint MLE for $({\mathbf{\Phi }}_1,\dots ,{\mathbf{\Phi }}_p)$. The residual covariance $\mathbf{Q}$ is then the sample covariance of the $k$-dimensional residual vector across observations.

(Granger Causality) Series $j$ Granger-causes series $i$ when lagged $j$ improves the prediction of $i$ beyond $i$‘s own past.

Test this with a standard $F$-test on the joint significance of $({\mathbf{\Phi }}_1{)}_{i,j},\dots ,({\mathbf{\Phi }}_p{)}_{i,j}$, the entries of the lag matrices that send series $j$ into series $i$. Granger causality is a statement about predictive content, not actual causality: lagged smoke can Granger-cause fire if smoke detectors are slow, even though fire causes smoke. It is useful for screening which cross-series links carry predictive information, and misleading if treated as a causal claim without further structural assumptions.

(Impulse Response) The effect on every series of a unit shock to one component.

For the MA($\infty$) representation ${\mathbf{x}}_t={\sum }_{j\ge 0}{\mathbf{\Psi }}_j{\mathbf{w}}_{t-j}$ derived from the inverted operator, the matrix ${\mathbf{\Psi }}_h$ gives the response of ${\mathbf{x}}_{t+h}$ to a unit innovation at time $t$. Because shocks are typically correlated across components ($\mathbf{Q}$ is non-diagonal), interpreting “a shock to series $j$” requires an orthogonalization, typically by Cholesky factorization of $\mathbf{Q}$, which depends on the order of the variables. Different orderings give different orthogonalized impulse responses, the well-known identification problem in structural VAR analysis.

(Variance Decomposition) Forecast error variance attributable to each shock.

The $h$-step-ahead forecast error variance of component $i$ decomposes as a sum over the $k$ orthogonalized shock channels and the $h$ time steps, and the decomposition tells you what fraction of the unpredictable variation in series $i$ is driven by each shock source. Reporting “shock to series 2 accounts for 30% of the 12-step forecast variance of series 1” is a staple of macroeconomic VAR analysis. Like impulse responses, the decomposition depends on the orthogonalization order.

(Lag Selection by Information Criteria) Choose $p$ by AIC, BIC, or HQ on the multivariate log-likelihood.

Compute the criterion for each candidate order and pick the minimum. BIC tends to select smaller orders than AIC; HQ (Hannan-Quinn) is intermediate. For modest $k$ all three criteria usually agree; for large $k$ they diverge and BIC’s preference for parsimony is usually the right default to avoid overfitting.

Examples

Bivariate VAR(1)

${\mathbf{x}}_t=\mathbf{\Phi }{\mathbf{x}}_{t-1}+{\mathbf{w}}_t$ with $\mathbf{\Phi }=\left(\begin{array}{cc}0.7& 0.1\\ 0.0& 0.5\end{array}\right)$, $\mathbf{Q}=\mathbf{I}$. What does the fitted matrix tell you?

Reading the structure $0.7$ (more persistent), series 2 with $0.5$ (less). The off-diagonal ${\Phi }_{1,2}=0.1$ says past values of series 2 do influence series 1, so there is a small spillover. The off-diagonal ${\Phi }_{2,1}=0$ says series 1 does not influence series 2 at all. Information flows one-way: $x_2\to x_1$. Stationarity check: eigenvalues of $\mathbf{\Phi }$ are $0.7$ and $0.5$, both inside the unit circle, so the process is stationary. A Granger-causality test would reject "series 2 does not Granger-cause series 1" (since ${\Phi }_{1,2}\ne 0$) but not "series 1 does not Granger-cause series 2."

The diagonal entries say each series has its own AR(1)-like persistence: series 1 with parameter

Bone-marrow transplant biomarkers

Three biomarkers ${\mathbf{y}}_t=(\log {\text{WBC}}_t,\log {\text{PLT}}_t,{\text{HCT}}_t{)}^{\top }$ are modeled jointly as a VAR(1).

What the fitted matrix reveals ${\Phi }_{\text{WBC},\text{PLT}}>0$ would mean that platelet recovery yesterday predicts white-blood-cell recovery today, consistent with the underlying bone-marrow function recovering as a single coordinated process. The $\mathbf{Q}$ matrix captures shocks that hit multiple markers on the same day. To handle the heavy missing-data pattern in this dataset (about 40% missing, especially after day 35), the VAR has to be embedded in a state-space model and inferred through the Kalman filter. VAR alone has no built-in mechanism for missing observations.

The diagonal entries say how persistent each marker is from one day to the next (close to one means slow change; close to zero means rapid turnover). The off-diagonal entries reveal cross-marker coupling. For example,

Lag selection on a 3-variable VAR

For three quarterly macro series with $n=200$ observations, BIC is computed for $p=1,2,3,4$.

Procedure $p$, fit the VAR by OLS, compute the residual covariance ${\hat{\mathbf{Q}}}_p$, and evaluate $\mathrm{BIC}(p)=\log \det {\hat{\mathbf{Q}}}_p+\frac{(\log n)k^{2}p}{n}.$ Pick the $p$ that minimizes BIC; macro data typically gives $\hat{p}=2$ or $\hat{p}=3$. The chosen lag order is then used for impulse-response analysis, Granger-causality tests, and forecasting.

For each candidate