Multistep prediction interval for ARMA(p,q) process

Question

How do I find an $h$-step prediction interval (forecast interval) for a zero-mean ARMA(p,q) process $$ x_t = \varphi_1 x_{t-1} + \dots + \varphi_p x_{t-p} + \varepsilon_t + \theta_1\varepsilon_{t-1} + \dots + \theta_q\varepsilon_{t-q} \ ? $$

Answer 1

The $h$-step-ahead point forecast, or the predicted conditional mean, is obtained by predicting one step ahead, $$ \hat x_{t+1} = \varphi_1 x_{t} + \dots + \varphi_p x_{t-p+1} + 0 + \theta_1 e_{t} + \theta_2 e_{t-1} + \dots + \theta_q e_{t-q+1}, $$ taking that forecast and plugging it in place of the true value $x_{t+1}$ and iterating such as $$ \hat x_{t+2} = \varphi_1 \hat x_{t+1} + \dots + \varphi_p x_{t-p+2} + 0 + \theta_1\times 0 + \theta_2 e_{t} + \theta_3 e_{t-1} + \dots + \theta_q e_{t-q+2}, $$ until $\hat x_{t+h}$ is reached. Here $e_t$ stands for the estimate of $\varepsilon_t$.

---

The $h$-step-ahead $(1-\alpha)$-level prediction interval (large sample approximation) is constructed as $$ [ \ \hat x_{t+h} - q_{\alpha/2}(\hat\sigma^2(h)); \ \hat x_{t+h} + q_{1-\alpha/2}(\hat\sigma^2(h)) \ ] $$ where $q_\alpha$ is the $\alpha$-level quantile of the error distribution with variance $\hat\sigma^2(h)$, $$ \hat\sigma^2(h) = \hat\sigma^2\sum_{j=0}^{h-1} \hat\psi_j^2 $$ where $\hat\sigma^2$ is the estimated error variance and $\hat\psi_j$ are the estimated coefficients of a moving-average representation of the ARMA(p,q) process.

Under normally distributed errors, the interval is $$ \hat x_{t+h} \pm z_{\alpha/2}\hat\sigma(h) $$ where $z_\alpha$ is the $\alpha$-level quantile of the N(0,1) distribution.

This derivation ignores parameter estimation uncertainty, so the actual intervals should be wider (recall the qualifier large sample approximation above). The problem diminishes with the sample size, though, and vanishes asymptotically.

---

(Based on Brockwell and Davis "Introduction to time series and forecasting" (3rd ed., 2016), p. 93-94.)

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

Answer 2

Practically it's the easiest to do it with Monte Carlo. I think that's what most people do anyways. You simulate the error term with and collect quantiles of the forecast. In R this can be done with arima.sim.

The quantiles are not so easy to estimate theoretically especially if the series are further combined with other predictions to obtain a composite forecast or when you model a transformed series, but are no brainer in simulation. The quantiles are popularized by BoE in fan charts in inflation report as shown below

Another application of Monte Carlo is for the case when you model correlated time series separately. Suppose, that for some reason you don't want to do a vector model, e.g. for complexity or lack of data. In this case you could still - similar to seemingly unrelated regression (SUR) - get the correlation of residuals of models, then generate correlated innovations and supply them into a command such as arima.sim to obtain correlated predictions