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

  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] +

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] +

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 +

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 +

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  
         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.

  • $\begingroup$ please post the 103 historical values and the estimated model parameters and summary statistics. $\endgroup$
    – IrishStat
    Commented Jul 26, 2019 at 13:27
  • $\begingroup$ 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 $\endgroup$
    – IrishStat
    Commented Jul 26, 2019 at 15:18
  • $\begingroup$ 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. $\endgroup$ Commented Jul 28, 2019 at 16:26


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.