bvarForecast#

Purpose#

Generate posterior predictive forecasts with credible bands from a fitted Bayesian VAR model.

Format#

fc = bvarForecast(result, h)#

fc = bvarForecast(result, h, xreg=X_future)

fc = bvarForecast(result, h, level=0.90)

Parameters:

result (struct) – an instance of a bvarResult structure returned by bvarFit().
h (scalar) – forecast horizon (number of steps ahead).
xreg (hxK matrix) – Optional keyword, future values of exogenous regressors. Required if the model was fit with xreg.
level (scalar) – Optional keyword, credible level for prediction bands. Default = 0.68.
quiet (scalar) – Optional keyword, set to 1 to suppress printed output. Default = 0.

Returns:

fc (struct) –

An instance of a forecastResult structure containing:

fc.forecasts	hxm matrix, point forecasts. hx1 for univariate (ARIMA), hxm for multivariate (VAR/BVAR).
fc.se	hxm matrix, forecast standard errors.
fc.lower	hxm matrix, lower confidence or credible bound.
fc.upper	hxm matrix, upper confidence or credible bound.
fc.level	Scalar, confidence or credible level used (e.g., 0.95).
fc.h	Scalar, forecast horizon.
fc.m	Scalar, number of variables (1 for ARIMA).
fc.var_names	Mx1 string array, variable names. Empty for univariate models.

Examples#

Default BVAR Forecast (68% Credible Bands)#

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

// Estimate and forecast
result = bvarFit(data, quiet=1);

struct forecastResult fc;
fc = bvarForecast(result, 12);

Forecast with 90% Credible Bands#

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));
result = bvarFit(data, quiet=1);

fc = bvarForecast(result, 12, level=0.90);

Optimal Hyperparameters to Forecast Pipeline#

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

// Optimize hyperparameters
ho = bvarHyperopt(data);

// Estimate with optimal lambdas
result = bvarFit(data, ho.ctl, quiet=1);

// Forecast
fc = bvarForecast(result, 24);

Compare Forecasts Across Lag Orders#

new;
library timeseries;

data = loadd(getGAUSSHome("pkgs/timeseries/examples/macro.dat"));

ctl = bvarControlCreate();

ctl.p = 2;
r2 = bvarFit(data, ctl, quiet=1);

ctl.p = 4;
r4 = bvarFit(data, ctl, quiet=1);

fc2 = bvarForecast(r2, 12, quiet=1);
fc4 = bvarForecast(r4, 12, quiet=1);

print "GDP forecast (p=2):" fc2.forecasts[., 1];
print "GDP forecast (p=4):" fc4.forecasts[., 1];

Remarks#

Posterior predictive distribution: For each posterior draw \((B^{(i)}, \Sigma^{(i)})\), the function simulates an h-step forecast path by iterating the VAR forward with innovations drawn from \(N(0, \Sigma^{(i)})\). The reported fc.forecasts is the median of the resulting predictive distribution. The fc.lower and fc.upper bands are quantiles at (1-level)/2 and (1+level)/2.

68% vs 95% bands: The default credible level of 0.68 corresponds to approximately \(\pm 1\) posterior standard deviation and is the convention in macroeconomic BVAR applications. For publication or comparison with frequentist intervals, use level=0.90 or level=0.95.

Conjugate vs Gibbs: For models fit with prior="minnesota" (conjugate), exact posterior draws are used. For prior="flat" (Gibbs), the retained MCMC draws are used. In both cases, the number of forecast draws equals result.n_draws.

Model#

For each posterior draw \((B^{(s)}, \Sigma^{(s)})\), the h-step forecast path is generated by iterating the VAR forward:

\[\hat{y}_{T+h}^{(s)} = \hat{B}_1^{(s)} \hat{y}_{T+h-1}^{(s)} + \cdots + \hat{B}_p^{(s)} \hat{y}_{T+h-p}^{(s)} + \hat{u}^{(s)} + \varepsilon_{T+h}^{(s)}\]

where \(\varepsilon_{T+h}^{(s)} \sim N(0, \Sigma^{(s)})\). The reported point forecast is the median across all draws; the credible bands are posterior quantiles.

This integrates over both parameter uncertainty (different \(B, \Sigma\) draws) and innovation uncertainty (random \(\varepsilon\)), giving a proper Bayesian predictive density. Algorithm ———

For each of the n_draws posterior draws \((B^{(s)}, \Sigma^{(s)})\):
1. Draw \(\varepsilon_{T+1}^{(s)}, \ldots, \varepsilon_{T+h}^{(s)} \sim N(0, \Sigma^{(s)})\).
2. Iterate the VAR forward using past data and previously generated forecasts.
Compute the median and quantiles of \(\{\hat{y}_{T+h}^{(s)}\}_{s=1}^{n\_draws}\) at each horizon.

Complexity: \(O(n\_draws \cdot h \cdot m^2 p)\). Sub-second for typical configurations. Troubleshooting —————

Forecast bands are too wide: Tighten the prior (reduce lambda1 in bvarFit()) or add sum-of-coefficients priors (lambda6 > 0). Wide bands often reflect parameter uncertainty in an over-parameterized model.

Forecasts revert to zero immediately: Check that ar is set correctly. With ar = 0 (white noise prior), the BVAR pulls forecasts toward zero. For levels data, use ar = 1.

Forecasts explode: The posterior mean may be non-stationary. Check result.is_stationary. Add regularization via sum-of-coefficients (lambda6) or single-unit-root (lambda7) priors in bvarFit(). Verification ————

BVAR forecasts verified against R BVAR::predict() and ECB BEAR BEARmain() forecast output. Point forecasts agree within Monte Carlo noise (different RNG streams). Prediction interval widths match R within 5% relative error.

See crossval/02_bvar_crossval.R for R comparison and the Verification and Cross-Validation page. References ———-

Banbura, M., D. Giannone, and L. Reichlin (2010). “Large Bayesian vector auto regressions.” Journal of Applied Econometrics, 25(1), 71-92.
Kadiyala, K.R. and S. Karlsson (1997). “Numerical methods for estimation and inference in Bayesian VAR-models.” Journal of Applied Econometrics, 12(2), 99-132.

Library#

timeseries

Source#

forecast.src