Multivariate Time Series Forecasting: What Changes When Series Interact

Forecasting several series separately is not wrong, it is incomplete in two specific ways: it drops the coefficients by which one series predicts another, and it drops the correlation between their forecast errors. Both are visible in a two-series example small enough to work by hand.

What a VAR adds

The vector autoregression is the base case and everything else in this area is a variation on it. As statsmodels writes it, a VAR(p) is

Y_t = nu + A_1 Y_{t-1} + ... + A_p Y_{t-p} + u_t     u_t ~ N(0, Sigma_u)

where Y_t is a K-vector holding all the series at time t, each A_i is a K×K matrix, and Sigma_u is the K×K covariance of the one-step errors. Running K univariate autoregressions instead is the same model with two restrictions imposed silently: every A_i is forced to be diagonal, and Sigma_u is treated as diagonal. The first restriction removes cross-series prediction. The second removes cross-series uncertainty. They are different losses and the second is the one people forget.

Two steps ahead, with and without the cross term

Take two series in deviations from their means: y1 is orders and y2 is site sessions, and sessions lead orders. Suppose the true system is a VAR(1) with

A = [ 0.5  0.4 ]        y1_t = 0.5 y1_{t-1} + 0.4 y2_{t-1} + u1_t
    [ 0.0  0.8 ]        y2_t =              0.8 y2_{t-1} + u2_t

state at time t:  y1 = 100,  y2 = 50

Forecast two steps with the full system:

h=1:  y1 = 0.5(100) + 0.4(50) = 50 + 20 = 70
      y2 =              0.8(50) = 40

h=2:  y1 = 0.5(70)  + 0.4(40) = 35 + 16 = 51
      y2 =              0.8(40) = 32

Now forecast y1 alone with the coefficient on its own lag, which is what a diagonal model gives you: 0.5 × 100 = 50, then 0.5 × 50 = 25. The full system says 70 and 51; the univariate one says 50 and 25. At two steps ahead it is out by a factor of two, and the error is not noise — it is the entire contribution of the 0.4 coefficient, compounded.

The deeper point is that the univariate model is misspecified by construction here, not merely less informed. The matrix A is triangular, so its eigenvalues are 0.5 and 0.8, and the marginal process that y1 follows on its own has both of those roots in it. No AR(1) can represent two roots. Fitting one to real data would land on some compromise coefficient between 0.5 and 0.8 that is right for neither horizon.

The covariance term in the error variance

The second loss is in the interval rather than the point. The h-step forecast error covariance of a VAR(1) is Sigma_u at h = 1 and Sigma_u + A Sigma_u A' at h = 2. Suppose the one-step errors have variances 4 and 9 with a correlation of 0.2, so the covariance is 0.2 × 2 × 3 = 1.2:

Sigma_u = [ 4.0  1.2 ]          A = [ 0.5  0.4 ]
          [ 1.2  9.0 ]              [ 0.0  0.8 ]

A Sigma_u        = [ 2.48  4.20 ]
                   [ 0.96  7.20 ]

A Sigma_u A'     = [ 2.92  3.36 ]
                   [ 3.36  5.76 ]

var(y1 error, h=2) = 4.00 + 2.92 = 6.92     sd = 2.63

now set the off-diagonal 1.2 to zero and repeat:
A Sigma_u A'     = [ 2.44  2.88 ]
                   [ 2.88  5.76 ]
var(y1 error, h=2) = 4.00 + 2.44 = 6.44     sd = 2.54

Ignoring a correlation of only 0.2 understates the two-step standard deviation by about 3.5 per cent, and the gap widens with the correlation, with the horizon and with the size of the off-diagonal coefficients. Modest on its own. It stops being modest the moment you add series together, because the variance of a sum carries the covariance at full weight: var(a + b) = var(a) + var(b) + 2cov(a, b). Sum ten series whose errors are pairwise correlated at 0.2 and treating them as independent understates the variance of the total by a large multiple, which is the same arithmetic that makes hierarchical reconciliation worth doing.

The same matrix is behind impulse response analysis, which is what people usually want a VAR for beyond forecasting: the question “if sessions are shocked by one unit today, what happens to orders over the next ten periods?”. Because the one-step errors are correlated, a shock to one series is never observed alone, so the answer is not well defined until you say which series moves first. statsmodels resolves this by orthogonalising through a Cholesky decomposition of Sigma_u, which imposes exactly that ordering: the first series in the vector is assumed to be able to affect the others contemporaneously while they cannot affect it. Reorder the series and you get different impulse responses from identical data. That is not a defect of the software, it is the identification problem made visible, and it is the reason an impulse response plot should always be published with its ordering.

What it costs to estimate

A VAR(p) on K series has K²p coefficients plus intercepts. Five series at four lags is a hundred coefficients; twenty series at four lags is sixteen hundred, from a sample that is however many observations you have. Parameters grow with the square of the number of series and the sample does not, which is the reason unrestricted VARs are rare above a handful of series.

The standard responses are all forms of restriction. Choose p by an information criterion rather than optimism — statsmodels exposes this as order selection over AIC, BIC, HQIC and FPE, and BIC will usually pick a shorter lag than AIC. Impose sparsity, so most cross terms are forced to zero. Impose structure, by replacing free cross terms with a small number of common factors. Or test first: statsmodels provides a Granger causality test, a Wald or F test of whether one series’ lags add explanatory power to another’s equation, which tells you which cross terms are worth keeping. Granger causality is a statement about predictive content and not about causation, and the name has misled people for fifty years.

When it is not worth it

• When the relationship is contemporaneous only. A VAR predicts from lags. If y2 moves with y1 in the same period and never before it, the lagged cross terms are near zero and the relationship lives entirely in Sigma_u, which improves your intervals and not your point forecasts.
• When you know the driver in advance. If the other series is a price you set or a holiday calendar you already have, treat it as a regressor rather than as a jointly modelled series; see forecasting with exogenous variables. Modelling something you control as random discards information.
• When the series are non-stationary and cointegrated. Differencing everything and fitting a VAR discards the long-run relationship between the levels. That is what a vector error correction model exists for, and using a plain VAR on differences here is a specific, named mistake rather than a rough approximation.
• When you have twelve observations per series. No amount of structure rescues a covariance matrix estimated from fewer observations than it has entries.