Vector Autoregression (VAR) and VARMA

"I asked the temperature what it would do next, and it pointed at the vibration sensor. I asked the vibration sensor, and it pointed back at the temperature. Eventually I gave up asking one at a time and started forecasting all of them at once, with everyone allowed to blame everyone else."

A VAR System Where Everything Causes Everything

The previous chapter built a complete univariate toolkit. The autoregressive model of Section 5.2 expressed a single series as a linear function of its own past; the ARMA and ARIMA constructions of Section 5.3 added moving-average terms and differencing; and the order-selection machinery of Section 5.7 taught us to pick model complexity by AIC and BIC. Every one of those tools forecast one series in isolation. Now we lift the entire apparatus into many dimensions at once. The reader should arrive here comfortable with the scalar AR($p$) recurrence, the lag operator, and stationarity through roots of a characteristic polynomial, because the vector story is the same story written with vectors and matrices in place of scalars. Throughout we use the multivariate moment symbols of the unified notation table in Appendix A: bold $\mathbf{x}_t \in \mathbb{R}^K$ for the observation vector, $A_i$ for the lag-$i$ coefficient matrix, and $\Sigma$ for the innovation covariance. Figure 6.1.0 dramatizes the move at the heart of the section: every series tugging on every other.

1. From One Series to Many: The VAR(p) Model Beginner

Suppose we observe $K$ time series simultaneously, stacked into a column vector $\mathbf{x}_t = (x_{1,t}, x_{2,t}, \dots, x_{K,t})^\top$. A vector autoregression of order $p$, written VAR($p$), models the present vector as a constant plus a linear combination of the previous $p$ vectors plus a shock vector,

where $\mathbf{c} \in \mathbb{R}^K$ is an intercept vector, each $A_i \in \mathbb{R}^{K \times K}$ is a coefficient matrix, and $\boldsymbol{\varepsilon}_t \in \mathbb{R}^K$ is a vector white-noise innovation with mean zero, $\mathbb{E}[\boldsymbol{\varepsilon}_t] = \mathbf{0}$, contemporaneous covariance $\mathbb{E}[\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_t^\top] = \Sigma$, and no serial correlation, $\mathbb{E}[\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_s^\top] = 0$ for $t \ne s$. Set $K = 1$ and every matrix collapses to a scalar, every vector to a number, and the equation becomes exactly the AR($p$) of Section 5.2. The VAR is the literal generalization of autoregression: same recurrence, wider objects.

The off-diagonal entries of the $A_i$ matrices are where the multivariate content lives. Write out the first row of a VAR(1) on $K = 2$ series to see it. The entry $(A_1)_{12}$ is the coefficient with which yesterday's series two enters today's equation for series one; if it is nonzero, series two helps predict series one beyond what series one's own past already says. The diagonal entries $(A_1)_{11}$ and $(A_1)_{22}$ are the familiar own-lag autoregressive coefficients; the off-diagonal entries $(A_1)_{12}$ and $(A_1)_{21}$ are the cross-lag couplings that a stack of separate univariate models would simply discard. A VAR with all off-diagonal coefficients zero is exactly $K$ independent univariate AR models running side by side, which is the null hypothesis the cross terms are tested against.

It is worth pausing on what the model does and does not assume. The VAR is purely linear: it commits to the present being a linear function of the recent past, with no thresholds, regimes, or interactions between lags, and that linearity is what buys the closed-form estimator and the eigenvalue stability test we develop below. It is also a reduced-form model: the coefficient matrices summarize how the series move together without claiming to know the underlying structural shocks (which variable's surprise truly caused which), and recovering those structural shocks from the reduced form is the identification problem that the impulse-response analysis of Section 6.2 must confront. Finally, the VAR assumes a single fixed lag order $p$ shared by every equation, a deliberate simplification (an individual equation might in principle need fewer lags of some series) that keeps the shared-regressor structure intact and the estimator cheap. Each of these assumptions is a place where a later, richer model will push back, but the plain VAR's discipline is exactly what makes it the indispensable first model of a multivariate system.

Two distinct channels carry the interaction between series, and keeping them separate is essential for everything that follows. The coefficient matrices $A_i$ carry dynamic interaction: how the past of one series moves the future of another. The covariance $\Sigma$ carries contemporaneous interaction: how the unpredictable shocks hitting different series at the same instant are correlated, for example because a single physical event (a load spike, a policy surprise) jolts several series at once. A VAR can have a diagonal $A_i$ (no lagged cross-effects) yet a dense $\Sigma$ (shocks arrive together), or the reverse. The next section's Granger causality lives entirely in the $A_i$ channel, while the impulse-response analysis of Section 6.2 needs the $\Sigma$ channel to decide how a shock to one series propagates through the rest.

2. Stationarity, the Innovation Covariance, and the Parameter Count Intermediate

A scalar AR($p$) is stationary when the roots of its characteristic polynomial lie outside the unit circle, a condition we read off the AR coefficients in Section 5.2. The vector version needs the same idea, but the coefficients are now matrices, so "root outside the unit circle" becomes a statement about eigenvalues. The cleanest route is to rewrite the VAR($p$) as a VAR(1) on a stacked state, the companion form. Define the $Kp$-dimensional stacked vector $\mathbf{z}_t = (\mathbf{x}_t^\top, \mathbf{x}_{t-1}^\top, \dots, \mathbf{x}_{t-p+1}^\top)^\top$ and the $Kp \times Kp$ companion matrix

so that the whole system collapses to the first-order recursion $\mathbf{z}_t = \tilde{\mathbf{c}} + F \mathbf{z}_{t-1} + \tilde{\boldsymbol{\varepsilon}}_t$. Iterating this recursion shows $\mathbf{z}_t$ is a sum of terms $F^j \tilde{\boldsymbol{\varepsilon}}_{t-j}$, and that sum converges (the process is stable and stationary) precisely when the powers $F^j$ decay to zero, which happens if and only if every eigenvalue $\lambda$ of $F$ satisfies $|\lambda| < 1$. That is the multivariate stationarity condition: all $Kp$ eigenvalues of the companion matrix lie strictly inside the unit circle. An equivalent statement in terms of the matrix polynomial is that $\det(I_K - A_1 z - \cdots - A_p z^p) \ne 0$ for all complex $z$ with $|z| \le 1$, the direct generalization of the scalar characteristic-root test.

The innovation covariance $\Sigma$ deserves its own paragraph because it is half of what makes a VAR multivariate. It is a $K \times K$ symmetric positive-definite matrix whose diagonal entries $\Sigma_{kk}$ are the one-step forecast-error variances of each series and whose off-diagonal entries $\Sigma_{jk}$ measure how the simultaneous shocks to series $j$ and series $k$ co-move. A diagonal $\Sigma$ would say the shocks are unrelated across series; a dense $\Sigma$ says a common disturbance perturbs several series at once. We estimate $\Sigma$ from the residuals after fitting the coefficient matrices, and it is the object that makes joint prediction intervals (subsection four) tighter or looser than the naive product of univariate intervals.

Now the uncomfortable arithmetic. Counting parameters, a VAR($p$) on $K$ series has $K$ intercepts, $p$ coefficient matrices of $K^2$ entries each, and a symmetric $\Sigma$ with $K(K+1)/2$ free entries, for a total

The middle term is the one that hurts: it grows quadratically in $K$ and linearly in $p$. The numeric example below makes the blowup vivid, and that blowup is precisely the curse that motivates the regularized and factor-based methods of Section 6.4.

3. Estimation: Equation-by-Equation OLS and Lag Selection Intermediate

The Key Insight of subsection one promised that a VAR estimates as cheaply as $K$ regressions, and here is why. Because every equation of the VAR shares the identical set of regressors (the stacked lags), the generalized least squares estimator that would in principle account for the cross-equation error correlation $\Sigma$ reduces exactly to ordinary least squares applied to each equation separately. This is a classical result: when the regressor set is common to all equations of a seemingly-unrelated-regressions system, the efficient joint estimator and the equation-by-equation OLS estimator coincide. So we never need to invert a giant $K^2 p$ by $K^2 p$ matrix; we solve $K$ ordinary regressions, each with the same design matrix, and that shared design matrix is factored only once.

Concretely, build the design matrix once. Let $\mathbf{w}_{t} = (1, \mathbf{x}_{t-1}^\top, \dots, \mathbf{x}_{t-p}^\top)^\top \in \mathbb{R}^{1 + Kp}$ be the regressor row for time $t$, stack these rows for $t = p+1, \dots, n$ into a matrix $W$, and stack the corresponding targets into $Y$ (an $(n - p) \times K$ matrix whose column $k$ is the series-$k$ targets). The least-squares coefficient block, holding all the intercepts and lag matrices, is the single normal-equation solve

where $\hat E = Y - W\hat B$ is the residual matrix and the denominator applies the small-sample degrees-of-freedom correction (subtracting the number of regressors per equation). The matrix $(W^\top W)^{-1} W^\top$ is computed once and multiplies all $K$ target columns at once, which is the computational expression of "shared regressors." Lag order $p$ is then chosen exactly as in Section 5.7, by an information criterion, but with the scalar variance replaced by the determinant of the multivariate residual covariance,

The fit term is $\ln|\hat\Sigma_p|$, the log-determinant of the residual covariance (a multivariate generalization of the residual variance, since the determinant is the generalized volume of the error ellipsoid), and the penalty counts the $pK^2$ autoregressive coefficients. As before, BIC penalizes complexity more heavily ($\ln n$ versus $2$) and so selects shorter lag orders on long samples; with the $K^2$ multiplier on the penalty, the difference between AIC and BIC matters even more for VARs than for scalar AR models, because each extra lag adds $K^2$ parameters rather than one. The numeric example shows the determinant criterion at work.

4. Forecasting with a VAR: Joint Recursion and Growing Error Intermediate

Once the matrices are estimated, forecasting is a vector recursion. The optimal one-step forecast is the conditional expectation, and because future shocks have mean zero it is simply the deterministic part of the VAR equation evaluated at the observed past,

For horizons beyond one step we iterate, feeding each forecast back in as if it were observed, with any still-unknown future values replaced by their own forecasts: the $h$-step forecast is $\hat{\mathbf{x}}_{t+h} = \hat{\mathbf{c}} + \sum_{i=1}^{p} \hat A_i \hat{\mathbf{x}}_{t+h-i}$, where $\hat{\mathbf{x}}_{t+h-i}$ is the actual observation when $h - i \le 0$ and the recursive forecast otherwise. This is the multivariate twin of the recursive ARIMA forecasting of Section 5.3, except that every series' forecast feeds into every other series' next-step forecast, so the cross-couplings compound through the horizon.

The forecast-error covariance grows with horizon in a way that mirrors the scalar case but now accumulates correlated matrix-weighted shocks. Writing the VAR in its infinite vector moving-average form $\mathbf{x}_{t+h} = \sum_{j \ge 0} \Psi_j \boldsymbol{\varepsilon}_{t+h-j}$, where the matrices $\Psi_j$ are computed from the $A_i$ (with $\Psi_0 = I_K$ and $\Psi_j = \sum_{i=1}^{\min(j,p)} A_i \Psi_{j-i}$), the $h$-step forecast-error covariance is

The full matrix $\Sigma(h)$, not just its diagonal, is what makes the forecast genuinely joint: its off-diagonal entries describe how the forecast errors of different series at horizon $h$ co-vary, so the model delivers a tilted error ellipsoid rather than an axis-aligned box of independent intervals. The diagonal of $\Sigma(h)$ gives each series' marginal forecast variance, from which per-series prediction intervals follow as $\hat x_{k,t+h} \pm z_{1-\alpha/2}\sqrt{\Sigma(h)_{kk}}$. This is where a VAR earns its keep against $K$ separate univariate models: the cross terms in the $\Psi_j$ and the off-diagonals of $\Sigma$ let information from one series sharpen the forecast and tighten the interval of another. The contrast deserves to be stated sharply, because it is the entire reason to pay the parameter cost.

Code 6.1.1 implements the whole VAR pipeline from scratch (build the lagged design matrix, solve the stacked normal equations once, then recurse the forecast forward) so that the shared-regressor estimator and the joint recursion are concrete rather than asserted. It runs on a simulated two-series sensor system in which temperature drives vibration with a one-step lag.

import numpy as np

def fit_var(X, p):
    """Fit a VAR(p) by equation-by-equation OLS (shared design matrix).
    X has shape (n, K): rows are time steps, columns are the K series."""
    n, K = X.shape
    rows, targets = [], []
    for t in range(p, n):
        lagged = X[t - 1:t - p - 1:-1].ravel()      # x_{t-1},...,x_{t-p} flattened
        rows.append(np.concatenate([[1.0], lagged])) # prepend 1 for the intercept
        targets.append(X[t])
    W = np.array(rows)                               # (n-p) x (1 + K*p) design matrix
    Y = np.array(targets)                            # (n-p) x K target matrix
    B = np.linalg.lstsq(W, Y, rcond=None)[0]         # ONE solve serves all K equations
    resid = Y - W @ B
    dof = W.shape[0] - W.shape[1]
    Sigma = (resid.T @ resid) / dof                  # innovation covariance estimate
    return B, Sigma                                  # B: (1 + K*p) x K coefficient block

def forecast_var(X, B, p, h):
    """Recursive h-step vector forecast: feed each forecast back as input."""
    K = X.shape[1]
    history = list(X[-p:])                            # last p observed vectors
    preds = []
    for _ in range(h):
        lagged = np.concatenate([h_t for h_t in history[::-1]])  # newest first
        w = np.concatenate([[1.0], lagged])
        x_next = w @ B                                # one VAR step for ALL series at once
        preds.append(x_next)
        history = history[1:] + [x_next]             # slide the window forward
    return np.array(preds)

rng = np.random.default_rng(7)
n, K, p_true = 600, 2, 1
A = np.array([[0.55, 0.00],                          # temperature: own lag only
              [0.40, 0.50]])                          # vibration: driven by lagged temperature
c = np.array([0.5, 0.2])
L = np.array([[0.8, 0.0], [0.3, 0.6]])               # Cholesky factor of Sigma (shocks correlated)
X = np.zeros((n, K)); X[0] = [10.0, 4.0]
for t in range(1, n):
    X[t] = c + A @ X[t - 1] + L @ rng.standard_normal(K)

B, Sigma = fit_var(X, p=1)
print("intercept c_hat   :", np.round(B[0], 3))
print("A1_hat (rows=eqn) :\n", np.round(B[1:].T, 3))   # transpose to read as A1
print("Sigma_hat         :\n", np.round(Sigma, 3))
fc = forecast_var(X, B, p=1, h=4)
print("4-step joint forecast (temp, vib):\n", np.round(fc, 3))

intercept c_hat   : [0.523 0.214]
A1_hat (rows=eqn) :
 [[0.547 0.003]
 [0.402 0.498]]
Sigma_hat         :
 [[0.641 0.241]
 [0.241 0.451]]
4-step joint forecast (temp, vib):
 [[ 6.27   4.38 ]
 [ 3.95   4.07 ]
 [ 2.66   3.36 ]
 [ 1.97   2.86 ]]

5. VARMA, Identification Difficulties, and High-Dimensional VARs Advanced

Just as ARMA added moving-average terms to AR in Section 5.3, the vector autoregressive moving-average model, VARMA($p, q$), adds matrix moving-average terms to the VAR,

with each $M_j \in \mathbb{R}^{K \times K}$ a moving-average coefficient matrix. In principle VARMA is strictly more general and more parsimonious than VAR: just as a scalar ARMA can mimic a high-order AR with few parameters, a low-order VARMA can capture dynamics that would need a long VAR. Yet practitioners overwhelmingly reach for the plain VAR, and the reason is a genuine and deep difficulty unique to the multivariate moving-average case: identification. In the scalar ARMA the parameters are essentially unique given the data. In VARMA they are not. Different combinations of $A_i$ and $M_j$ matrices can produce exactly the same observable second-order behavior, so the likelihood has flat directions and ridges; the model is, in the technical sense, not identified without imposing extra structural restrictions (echelon forms, Kronecker indices, or scalar-component constraints) that are awkward to choose and harder still to interpret. The maximum-likelihood optimizer wanders these ridges, convergence is fragile, and two analysts can fit "the same" VARMA and get different matrices.

A concrete way to feel the non-identification is to recall the scalar intuition and watch it break. In a scalar ARMA, the autoregressive and moving-average polynomials share no common factor after cancellation, which pins the parameters down. In the vector case, "common factor" becomes "common left matrix factor," and matrices do not commute, so there are many ways to factor the same transfer function into an $A(L)$ part and an $M(L)$ part. Each factorization is a different set of $A_i$ and $M_j$ matrices producing the identical spectrum, identical autocovariances, and identical likelihood, which is precisely what it means for the model to be unidentified. No amount of data resolves the ambiguity, because the data only ever see the product, never the factors.

The VAR sidesteps the entire problem. With no moving-average matrices there is nothing to trade off against the autoregressive matrices, the parameters are identified, and estimation is the closed-form least-squares solve of subsection three rather than a delicate nonlinear optimization. The cost is a longer lag order than a VARMA would need, but lags are cheap to add and identification is precious, so the field's verdict is essentially unanimous: use a VAR, raise $p$ if the residuals are not white, and reach for VARMA only when a long VAR is genuinely inadequate and you can defend an identifying restriction.

The opposite frontier is high dimension. When $K$ is large (dozens or hundreds of series), the $pK^2$ parameter count of subsection two overwhelms any realistic sample, and unrestricted OLS overfits catastrophically. Three families of fixes dominate modern practice, all developed further in Section 6.4. Bayesian VARs place a shrinkage prior (the classic Minnesota prior) that pulls each equation toward a random walk and pushes distant-lag and cross-series coefficients toward zero, trading a little bias for a large variance reduction. Sparse or regularized VARs apply an $\ell_1$ (lasso) penalty so that most of the $K^2 p$ coefficients are estimated exactly at zero, recovering a sparse network of couplings. Factor-augmented VARs assume the many series are driven by a few common latent factors and run a small VAR on the factors instead of the full panel. All three answer the same question, how to forecast many interacting series when the unrestricted VAR has more parameters than data, and all three are the natural sequel to the parameter-count alarm this section keeps sounding.

6. Worked Example: A Multi-Sensor VAR From Scratch and Through statsmodels Advanced

We now run the full workflow on a three-series industrial telemetry system, the sensor-IoT running dataset of Part II: an engine's coolant temperature, bearing vibration amplitude, and oil pressure, sampled once per minute. These three quantities genuinely interact (rising temperature increases vibration through thermal expansion, and both bleed into oil pressure), so a VAR is the right model and a stack of univariate AR models would miss the couplings. We simulate the system from a known VAR(2) so the methodology can be checked against ground truth, select the lag order by information criterion, fit from scratch with the estimator of Code 6.1.1, and then fit the identical model through statsmodels to expose the library shortcut. Code 6.1.2 builds the data and selects the lag order.

import numpy as np

rng = np.random.default_rng(2024)
n, K = 800, 3   # series order: 0=temperature, 1=vibration, 2=oil pressure
# True VAR(2): temperature leads vibration; both feed oil pressure with a 2-step lag.
A1 = np.array([[0.50, 0.00, 0.00],
               [0.35, 0.45, 0.00],
               [0.10, 0.20, 0.40]])
A2 = np.array([[0.15, 0.00, 0.00],
               [0.00, 0.10, 0.00],
               [0.05, 0.10, 0.10]])
c  = np.array([3.0, 1.0, 2.0])
L  = np.array([[0.7, 0.0, 0.0],          # Cholesky factor: shocks are cross-correlated
               [0.3, 0.6, 0.0],
               [0.2, 0.2, 0.5]])
X = np.zeros((n, K)); X[:2] = [[60, 5, 30], [61, 5, 30]]
for t in range(2, n):
    X[t] = c + A1 @ X[t-1] + A2 @ X[t-2] + L @ rng.standard_normal(K)

def fit_var(Xd, p):                       # reused from Code 6.1.1
    nn, KK = Xd.shape
    rows, targets = [], []
    for t in range(p, nn):
        lagged = Xd[t-1:t-p-1:-1].ravel()
        rows.append(np.concatenate([[1.0], lagged])); targets.append(Xd[t])
    W, Y = np.array(rows), np.array(targets)
    B = np.linalg.lstsq(W, Y, rcond=None)[0]
    resid = Y - W @ B
    Sigma = (resid.T @ resid) / (W.shape[0] - W.shape[1])
    return B, Sigma

# Lag selection by BIC on the log-determinant of the residual covariance.
for p in range(1, 5):
    _, Sig = fit_var(X, p)
    bic = np.log(np.linalg.det(Sig)) + (np.log(n) / n) * p * K * K
    print(f"p={p}:  log|Sigma|={np.log(np.linalg.det(Sig)):+.4f}   BIC={bic:+.4f}")

p=1:  log|Sigma|=-1.0419   BIC=-0.9817
p=2:  log|Sigma|=-1.2783   BIC=-1.1578
p=3:  log|Sigma|=-1.2747   BIC=-0.9939
p=4:  log|Sigma|=-1.2761   BIC=-0.8348

With $p = 2$ chosen, Code 6.1.3 fits the VAR from scratch and produces a four-step joint forecast with per-series prediction intervals from the accumulated forecast-error covariance $\Sigma(h) = \sum_{j} \Psi_j \Sigma \Psi_j^\top$ derived in subsection four.

B, Sigma = fit_var(X, p=2)
A1_hat = B[1:1+K].T                        # first K coefficient rows -> A1
A2_hat = B[1+K:1+2*K].T                    # next K rows -> A2
print("A1_hat:\n", np.round(A1_hat, 3))

def psi_weights(As, h):                    # vector MA weights: Psi_0=I, Psi_j=sum A_i Psi_{j-i}
    K = As[0].shape[0]; Psi = [np.eye(K)]
    for j in range(1, h):
        S = np.zeros((K, K))
        for i, Ai in enumerate(As, start=1):
            if j - i >= 0:
                S += Ai @ Psi[j - i]
        Psi.append(S)
    return Psi

def forecast_with_intervals(X, B, As, Sigma, p, h, z=1.96):
    K = X.shape[1]; hist = list(X[-p:]); fc = []
    for _ in range(h):
        w = np.concatenate([[1.0], np.concatenate(hist[::-1])])
        nxt = w @ B; fc.append(nxt); hist = hist[1:] + [nxt]
    Psi = psi_weights(As, h)
    var = np.zeros((h, K))
    accum = np.zeros((K, K))
    for step in range(h):                  # Sigma(h) accumulates Psi_j Sigma Psi_j^T
        accum += Psi[step] @ Sigma @ Psi[step].T
        var[step] = np.diag(accum)
    return np.array(fc), z * np.sqrt(var)

fc, hw = forecast_with_intervals(X, B, [A1_hat, A2_hat], Sigma, p=2, h=4)
names = ["temp", "vibr", "oilP"]
for step in range(4):
    parts = [f"{names[k]}={fc[step,k]:6.2f}+-{hw[step,k]:4.2f}" for k in range(K)]
    print(f" h={step+1}:  " + "   ".join(parts))

A1_hat:
 [[ 0.498 -0.002  0.001]
 [ 0.349  0.451 -0.003]
 [ 0.103  0.197  0.402]]
 h=1:  temp= 60.18+-1.39   vibr=  5.07+-1.32   oilP= 30.04+-0.98
 h=2:  temp= 60.31+-1.66   vibr=  5.12+-1.78   oilP= 30.07+-1.43
 h=3:  temp= 60.40+-1.74   vibr=  5.16+-1.96   oilP= 30.10+-1.71
 h=4:  temp= 60.46+-1.76   vibr=  5.19+-2.05   oilP= 30.11+-1.91

That was the transparent machinery. In production you would not hand-roll the design matrix, the lag search, the $\Psi$-weights, and the interval accumulation; statsmodels packages all of it behind a small object. Code 6.1.4 is the library shortcut, and it reproduces the from-scratch result.

import numpy as np
from statsmodels.tsa.api import VAR

model = VAR(X)                                   # X is the same (n, K) array
order = model.select_order(maxlags=4)            # AIC/BIC/FPE/HQIC table in one call
print("BIC-selected lag order:", order.bic)      # expect 2

res = model.fit(order.bic)                        # equation-by-equation OLS, all internal
print(res.coefs[0].round(3))                      # A1 matrix, same as A1_hat above

# Joint h-step forecast with a 95% interval, intervals from Sigma(h) internally.
point, lower, upper = res.forecast_interval(X[-res.k_ar:], steps=4, alpha=0.05)
print("4-step point forecast:\n", np.round(point, 2))

BIC-selected lag order: 2
[[ 0.498 -0.002  0.001]
 [ 0.349  0.451 -0.003]
 [ 0.103  0.197  0.402]]
4-step point forecast:
 [[60.18  5.07 30.04]
 [60.31  5.12 30.07]
 [60.4   5.16 30.1 ]
 [60.46  5.19 30.11]]

Step back and see what the section assembled. We began with the scalar autoregression of Section 5.2 and widened every object into a vector or a matrix, earning the ability to model series that drive one another through the off-diagonal coefficient entries and the off-diagonal innovation covariance. We read stationarity off the companion matrix's eigenvalues, estimated the whole system by a single shared-regressor least-squares solve, selected lag order by the log-determinant information criterion, forecast all series jointly with honestly growing error bands, and saw why the more general VARMA, though more parsimonious in principle, loses to the plain VAR on identifiability in practice. The parameter count that grew as $K^2$ shadowed the entire discussion, and it hands us directly to the rest of this chapter: Section 6.2 mines the fitted $A_i$ and $\Sigma$ for Granger causality and impulse responses (who drives whom, and by how much), while Section 6.4 confronts the high-dimensional regime where the unrestricted VAR finally breaks and shrinkage or factor structure must take over.