How to make $h$-step interval forecasts from an ARMA-GARCH model?

Question

I recently wrote various Python functions to fit ARMA models and make forecasts from them. I am now trying to do the same for ARMA-GARCH models.

To make $h$-step forecasts from ARMA models, I used the innovations algorithm, as described in Brockwell's Introduction to Time Series Analysis and Forecasting. I was wondering whether I could still use this algorithm to make $h$-step forecasts from ARMA-GARCH models?

If not, is there any good literature out there which explains how to make $h$-step interval forecasts from ARMA-GARCH models, and how to compute the prediction error?

Answer 1

An ARMA(p,q)-GARCH(r,s) model specifies the conditional distribution of a time series:

\begin{aligned} x_t &= \mu_t + u_t, \\ \mu_t &= \varphi_0 + \varphi_1 x_{t-1} + \dots + \varphi_p x_{t-p} + \theta_1 u_{t-1} + \dots + \theta_q u_{t-q}, \\ u_t &= \sigma_t \varepsilon_t, \\ \sigma_t^2 &= \omega + \alpha_1 u_{t-1}^2 + \dots + \alpha_s u_{t-s}^2 + \beta_1 \sigma_{t-1}^2 + \dots + \beta_r \sigma_{t-r}^2, \\ \varepsilon_t &\sim i.i.d.(0,1). \end{aligned}

(See this answer for an unorthodox, though hopefully enlightening representation of ARMA-GARCH.)

Point forecast
The conditional mean of the distribution is given solely by the ARMA conditional mean equation -- the equation for $\mu_t$. Hence, if the point forecasts are the predicted conditional means (which is a popular choice and is optimal under square loss), the point forecasts from an ARMA-GARCH model will be determined entirely by the estimated ARMA equation. (However, the estimation of the ARMA equation will be affected by the choice of the conditional variance equation and hence will generally differ in presence vs. in absence of GARCH.)

Interval forecast for GARCH with a constant mean (no ARMA)
By the assumption that the innovations in the GARCH process are uncorrelated, the variance of their sum equals the sum of their variances, $$ \text{Var}(u_{t+1}+\dots+u_{t+h})=\text{Var}(u_{t+1})+\dots+\text{Var}(u_{t+h}). $$ A $(1-\alpha)\times 100\%$ level forecast interval is $$ [\ \hat\mu + q_{\alpha/2}(\sigma^2(h)); \ \hat\mu + q_{1-\alpha/2}(\sigma^2(h)) \ ], $$ where $\hat\mu$ is obtained from the conditional mean equation, $\sigma^2(h)=\sum_{i=1}^h \hat\sigma^2_{t+i}$, $\hat\sigma^2_{t+i}$ is obtained iteratively from the conditional variance equation, and $q_{\alpha}$ is the $\alpha$-level quantile of the relevant distribution (Normal if $\varepsilon$ is assumed to be Normal and potentially ugly otherwise; in the latter case, simulation may be needed to obtain the relevant quantiles).

Interval forecast for ARMA-GARCH
To be added.
A simplified answer is available at Hyndman and Athanasopoulos "Forecasting: Principles and Practice", specifically Section 8.8 "Forecasting"; see also Section 3.5 "Prediction intervals" and the thread "Multistep prediction interval for ARMA(p,q) process".

---

All of the above ignores parameter estimation uncertainty, so the actual intervals should be wider. The problem diminishes with the sample size, though, and vanishes asymptotically.

_{Keywords: multi-step, multi-period, multistep, multiperiod, multiple step, multiple period, steps ahead, periods ahead, forecast, predict, forecasting, prediction, point, interval, ARMA, ARIMA, GARCH, volatility, variance.}

Answer 2

I wanted to stress that Interval forecast for ARMA-GARCH seems incorrect. $\sigma^2(h)$ by definition is always growing in time. However, using that formula could lead to a decreasing $\sigma^2(h)$. Furthermore, if your ARIMA is using less lags (k) than the predicted interval (h), that means that the first h-k predicted variances, are not being used to forecast $\sigma^2(h)$, so could lead to decreasing forecasted $\sigma^2(h)$ (since those h-k coefficients are equal to 0))