Recursive VAR Forecasts and Prediction Errors

Why Recursive Forecasts Need an Uncertainty Calculation

A fitted vector autoregression (VAR) produces forecasts for several related time series at once. At a one-step horizon, the forecast follows directly from the model. At longer horizons, future values on the right-hand side are unknown, so forecasts must be built recursively: use the first forecast to form the second, the first two forecasts to form the third, and so on. The same recursion determines how future shocks accumulate into prediction errors.

The distinction matters in practice. A sequence of forecast values describes the model’s conditional mean path, not the range of plausible future outcomes. To quantify that range, we need the covariance matrix of the forecast error at each horizon. The formulas below condition on the fitted coefficients and use the model’s innovation covariance; they do not include uncertainty from estimating the coefficients.

Model and Recursive Forecasts

Let \(y_t\) be a \(k\)-dimensional vector, and consider a fitted VAR(\(p\)):

Here \(c\) is a \(k\)-vector, each \(A_j\) is a \(k\times k\) coefficient matrix, and \(u_t\) is the innovation vector. For the forecast-error formulas, assume that for all future dates \(s,r\geq1\), \(E(u_{T+s}\mid\mathcal{F}_T)=0\) and \(E(u_{T+s}u_{T+r}'\mid\mathcal{F}_T)=\mathbf{1}\{s=r\}\Sigma_u\). The information available through date \(T\) is denoted by \(\mathcal{F}_T\). Future innovations are unpredictable given this information.

In practice, the coefficients and innovation covariance are estimated. In the equations below, hats on the coefficients indicate fitted values; the forecasts are conditional on those fitted values. If the model also contains deterministic terms that are known at future dates, include their forecast values in the recursion along with the intercept.

For a fitted VAR(\(p\)), the recursion is

Thus each forecast uses the actual observed history available at \(T\), then replaces unavailable future values with forecasts. This is a conditional-mean forecast under the stated zero-mean innovation assumption.

Forecast Errors and Their Covariance

Define the \(h\)-step forecast error as the realized future vector minus its forecast:

The coefficient matrices that describe how innovations enter forecasts are defined recursively. Set \(\Psi_0=I_k\), the \(k\times k\) identity matrix, and for \(m\geq1\) define

These are the moving-average coefficients of the VAR. For a VAR(1), for example, \(\Psi_m=A_1^m\). For a known-parameter VAR whose forecast recursion uses the same coefficients and intercept as the data-generating process, the expansion below is exact. For a fitted model, write \(\hat\Psi_i\) for the moving-average coefficients computed from the fitted coefficients; the displayed formula is a plug-in approximation that treats those coefficients as the working model’s true coefficients, not an exact expansion of the actual forecast error when fitted and true parameters differ.

This expression includes the innovation at the forecast target date, multiplied by \(\Psi_0=I_k\), and earlier future innovations multiplied by the appropriate propagation matrices. Only innovations after \(T\) appear because observed values through \(T\) are already known when the forecast is made.

For estimated models, the practical calculation substitutes the fitted \(\hat\Psi_i\) and an estimate \(\hat\Sigma_u\). The covariance generally grows with the horizon because each additional step introduces another future innovation. It need not grow in every individual variance or covariance entry, but it is nondecreasing in the positive-semidefinite matrix ordering: \(V_{h+1}-V_h=\Psi_h\Sigma_u\Psi_h'\), which is positive semidefinite.

Why the Covariance Formula Holds

For a known-parameter VAR with the same parameters in the data-generating equation and forecast recursion, the first-horizon error is \(e_{T+1\mid T}=u_{T+1}\), so its coefficient on the innovation is \(\Psi_0=I_k\). At later horizons, subtracting the forecast recursion from the VAR equation cancels the intercept and known terms, leaving the current innovation plus coefficient matrices multiplying errors at earlier future dates. Substituting those earlier error expansions yields the coefficients \(\Psi_m\) defined above and proves the exact expansion. For a fitted model, this derivation applies to the working model that treats the fitted parameters as true; the resulting expansion and covariance are plug-in approximations, not exact statements about the actual forecast errors.

To obtain the covariance, substitute the error expansion into its variance. Terms involving innovations at different dates have zero covariance by the serial-uncorrelatedness assumption. The remaining term for innovation \(u_{T+h-i}\) is \(\Psi_i\Sigma_u\Psi_i'\). Summing over \(i=0,\ldots,h-1\) gives \(V_h\). In particular, for the first two horizons,

This result concerns forecast uncertainty under the specified innovation process. It does not require a stable VAR to calculate a finite-horizon forecast, although stability is important for the long-horizon behavior of the process. If innovations are correlated across time, the omitted cross-covariance terms must also be included.

Prediction Intervals and Joint Prediction Regions

If future innovations are jointly Gaussian conditional on the information at \(T\), the forecast error is Gaussian with mean zero and covariance \(V_h\). A marginal prediction interval for component \(r\) uses the corresponding diagonal entry \((V_h)_{rr}\). At the 95% level, it is the forecast plus or minus \(1.96\sqrt{(V_h)_{rr}}\). This is an interval for a future realized value, in the units of that variable.

A joint region for all \(k\) forecast components requires care about the rank of \(V_h\). If \(V_h\) is positive definite, the Gaussian squared distance \((y-\hat y)'V_h^{-1}(y-\hat y)\) has a chi-squared distribution with \(k\) degrees of freedom. The 95% joint prediction region is therefore the ellipsoid defined by a cutoff of \(\chi^2_{k,0.95}\). If \(V_h\) has rank \(r<k\), the inverse and the \(k\)-degree-of-freedom result do not apply. A Gaussian forecast error then lies on the \(r\)-dimensional support of \(V_h\); use its Moore–Penrose pseudoinverse \(V_h^+\), restrict the region to that support, and use a \(\chi^2_{r,0.95}\) cutoff.

Interpretation and Common Pitfalls

The formulas give forecast uncertainty conditional on the fitted model and its estimated coefficients. In an applied analysis, replacing \(\Sigma_u\) with \(\hat\Sigma_u\) and using fitted coefficients makes the calculated interval or region a plug-in approximation. It reflects future innovation uncertainty, but not uncertainty in the estimated VAR parameters. If the innovation distribution is not Gaussian, the covariance calculation can still describe second-moment uncertainty under the assumptions above, but Gaussian prediction cutoffs need not give exact coverage.

• Using actual future values during recursion. The forecast at horizon \(h\) cannot use \(y_{T+1},\ldots,y_{T+h-1}\), which are not yet observed. Use the forecasts previously produced for those dates.
• Using only the current innovation covariance at every horizon. At horizon \(h>1\), earlier future innovations also affect the outcome. Include each term \(\Psi_i\Sigma_u\Psi_i'\) for \(i=0,\ldots,h-1\).
• Calling a marginal interval a joint region. Separate 95% intervals for each component do not in general provide 95% simultaneous coverage for the whole vector. Use a joint region when the claim concerns all components together.
• Inverting a singular covariance matrix. A Gaussian assumption alone does not guarantee that \(V_h\) is nonsingular. Check its rank; for rank \(r<k\), use the pseudoinverse on the covariance support and a \(\chi^2_r\) cutoff.
• Treating a plug-in interval as if it included every source of uncertainty. The covariance formula here accounts for future innovations conditional on fitted coefficients; it does not automatically account for coefficient-estimation uncertainty.