# Estimating prediction interval of ARMA process using R forecast function

the theme is forecasting with ARMA models.

I'm trying to understand how the R forecast function works if applied to an Arima object and, in particular, how the prediction interval is computed.

In the following code I'm fitting a signal called variable to predict the next 4 points (from $$104$$ to $$107$$).

x.fit.c <-Arima(variable, order = c(1, 0, 3),include.mean = FALSE, include.drift = FALSE)
forecast<-forecast(x.fit.c, h = 4, level=0.999)

> forecast
Lo.99.9       Point.Forecast  Hi.99.9
104   -37.84        -23.65          -9.47
105   -52.70        -22.06           8.58
106   -70.39        -20.50           29.40
107   -90.51        -19.42           51.67


Now, how does the forecast function work? For each step $$h$$, there are two output:

• prediction (Point.Forecast, in the code);
• prediction interval (Lo 99.9 and Hi 99.9 in the code).

A generic ARMA model, built on a sample of $$n$$ observations, has the following equation:

$$y_t = \mu + \sum \phi_iy_{t-i}+ \sum \theta_j\epsilon_{t-j}$$

then, how to calculate these outputs? I describe it in detail, hoping to be useful, and then ask the question at the end.

## How to compute prediction

The prediction is obtained by updating $$t$$ with $$t+h$$ and by using the following rules:

• For any $$\epsilon_j$$ with $$j \in [1,n]$$ , use the sample residual for time point $$j$$;
• For any $$\epsilon_j$$ with $$j > n$$ , use 0 as the value of $$\epsilon_j$$;
• For any $$y_j$$ with $$j \in [1,n]$$ , use the observed value of $$y_j$$;
• For any $$y_j$$ with $$j > n$$ , use the forecasted value of $$y_j$$;

If we assume that $$\epsilon \sim N(0,\hat{\sigma}_{residuals})$$ , a step-by-step Montecarlo simulation to obtain the first 4 forecasted points is the following:

  N.repl <- 1000000
set.seed(1234)

sim.cond <-
rnorm(N.repl,mean=0,sd = sqrt(x.fit.c$$sigma2)) + x.fit.c$$coef[1]*x.fit.c$$x[103] + x.fit.c$$coef[2]*x.fit.c$$residuals[103] + x.fit.c$$coef[3]*x.fit.c$$residuals[102] + x.fit.c$$coef[4]*x.fit.c$residuals[101] sim.cond_2 <- rnorm(N.repl,mean=0,sd = sqrt(x.fit.c$$sigma2)) + x.fit.c$$coef[1]*mean(sim.cond) + x.fit.c$$coef[2]*0 + x.fit.c$$coef[3]*x.fit.c$$residuals[103] + x.fit.c$$coef[4]*x.fit.c$residuals[102]

sim.cond_3 <-
rnorm(N.repl,mean=0,sd = sqrt(x.fit.c$$sigma2)) + x.fit.c$$coef[1]*mean(sim.cond_2) +
x.fit.c$$coef[2]*0 + x.fit.c$$coef[3]*0 +
x.fit.c$$coef[4]*x.fit.c$$residuals[103]

sim.cond_4 <-
rnorm(N.repl,mean=0,sd = sqrt(x.fit.c$$sigma2)) + x.fit.c$$coef[1]*mean(sim.cond_3) +
x.fit.c$$coef[2]*0 + x.fit.c$$coef[3]*0 +
x.fit.c$coef[4]*0  Clearly the prediction, at step $$h$$, is the mean of the relatve sim.cond. So far so good: this works pretty good (and hope could be usefull to beginners!). ## Prediction Interval Here there is a problem! The literature states that to calculate the interval prediction we have to writing the ARMA model (p, q) in alternative form, a Moving Average (MA) model of infinite order: $$y_t - \mu = 1+\sum_{i=1}^{\infty} \psi_i\epsilon_{t-i}$$ where $$\psi_0 =1$$ by definition. Then the prediction interval at $$h$$ step is: $$[\hat{y}_{t+h} - q_{\alpha/2} \hat{\sigma}(h) , \hat{y}_{t+h} + q_{\alpha/2} \hat{\sigma}(h)]$$ where is clear that $$\hat{y}_{t+h}$$ is the forecasted value at $$h$$ step, $$q_{\alpha}$$ is the $$\alpha$$-level quantile of the error distribution with variance $$\hat{\sigma}^2(h)$$, and $$\hat{\sigma}^2(h)=\hat{\sigma}^2 \sum_{i=0}^{h-1}\hat{\psi}^2_i$$ Here, $$\hat{\sigma}^2$$ is the estimated error variance and $$\hat{\psi}_i$$ are the estimated coefficients of a moving average (MA) representation of the ARMA(p,q) process. Under normally distributed errors, the interval is: $$\hat{y}_{t+h} \pm q_{\alpha/2} \hat{\sigma}(h) ,$$ where $$q_{\alpha/2}$$ is the $$\alpha$$-level quantile of $$N(0,1)$$ distribution (es, qnorm(0.999)=3.090232). The code which compute this interval for value of $$h>1$$ is: h=2 #from h+1 to h sqrt(x.fit.c$$sigma2)* qnorm(0.999)* sqrt(1+ sum(ARMAtoMA(ar=x.fit.c$$coef[1], ma=c(x.fit.c$$coef[2], x.fit.c$$coef[3], x.fit.c$coef[4]), lag.max = h-1)^2))


(While for $$h=1$$ is simply qnorm(0.999)*sqrt(x.fit.c$sigma2)). ## Results and final question Following the two previous section, we obtain:  Lo.99.9 Point.Forecast Hi.99.9 104 -36.96 -23.65 -10.35 105 -50.80 -22.06 6.68 106 -67.32 -20.50 26.33 107 -86.16 -19.42 47.32  Please, note that these results are a bit different from those obtained using forecast function. Why? Computing the implicit standard deviation from the results of the forecast function, by inverting $$\mu + q_{\alpha/2}\sigma=Hi.99.9$$, I observe that forecast estimates a different variance with respect to the one estimated from regression residuals (x.fit.c$sigma2). The first is a bit grather. Why this difference?

I also observed that reducing the interval (that is, taking $$\alpha$$ lower than 0.999, like 0.9 or 0.75) I get a greater implicit variance value and, vice versa, as $$\alpha$$ tends to 1, the value of variance used by forecasts tends to those estimated from the regression residuals.

I am interested in knowing which of the two practices is correct (the one I dictated or that implemented by forecast) and, if the best practice is the forecast one, how this function calculates $$\hat{\sigma}$$ in $$\hat{\sigma}^2(h)=\hat{\sigma}^2 \sum_{i=0}^{h-1}\hat{\psi}^2_i$$.

@ IrishStat

It is important to understand that what brought me to this question is not a comparison between two methods (# 1 and # 2, to use your taxonomy), but a way to understand well what the forecast function does and if what it does it's correct.

1) Different literature can be found simply by searching on Google: for this I'll avoid to write down links or names of books easily identifiable looking on the web. I note, however, that the image you reported is taken from the statistics course of the Eberly College of Science (https://newonlinecourses.science.psu.edu/stat510/lesson/3/3.3) which can serve as an example.

Instead, what approaches did you see in the forecasting of ARMA models?

2) Yes, is exactly what I wrote. Or I'm missing something, according to you?

3) Is true. In fact I have no pretense of generality with this method, although I can say that it works in many cases, because often the residues - normal or not - are symmetrically distributed around an average value. If they are not, there could be a problem in the model, which does not capture the information well, leaving a substantial part in the residues.

• please post the 103 historical values and the estimated model parameters and summary statistics. Commented Jul 26, 2019 at 13:27
• A couple of reflections .. you point to two solutions ... forecast package is # 1 . 1) Can you provide a citation for what you did in your solution # 2 . I have never seen nor do I understand this approach. 2) The Psi Weights that are used in solution #1 are obtained by expressing the arima model as a pure moving-average model rather than as a mixed model or a pure autoregressive model. 3) Both results can be seriously flawed by a distribution of errors that deviated from normality possibly due to anomalies or an inherently non-normal distribution as both formulas premise symmetry in computing Commented Jul 26, 2019 at 15:18
• 1) You should check the stats::arima function for the intervals. The forecast package Arima is just a wrapper to select suitable hyperparameters. The actual estimation, including the intervals, is done in stats::arima, so you should check that package documentation. 2) The estimation is done using Kalman Filter: in order to find how variance is estimated refer to the 13th chapter of Hamilton TSA. Commented Jul 28, 2019 at 16:26