Simulation of forecasted values in ARIMA (0,1,1)

Question

I have fitted an ARIMA (0,1,1) model with a drift term using the “forecast” package. I want to perform a simulation study to obtain the mean of the forecasted values and the 95% forecasting intervals and compare them with the “Exact” mean and 95 forecasting intervals. The code below is what I have done. The graph shows the simulated and exact results.

Question: Why the upper and lower confidence intervals based on the simulation are different from the exact values? While the mean are close to the exact value. Where is my mistake?

I add that I know there is a “simulate” function provided in this package to do the simulation. But this function produces the NA for me and is not working properly, so I tried to do it myself!

library(forecast)
f4=Arima(WWWusage,order=c(0,1,1),include.drift=TRUE)
mean.exact.f=forecast(f4,h=21)$mean
U.exact.f=forecast(f4,h=21)$upper[,2]
L.exact.f=forecast(f4,h=21)$lower[,2]
ff.sim=function(m,h,N){
sigma.est=m$sigma
ma1.est=coef(f4)[1]
drf.est=coef(f4)[2]
ff=matrix(0,ncol=N,nrow=(h+1))
for(j in 1:N){
              ff[1,j]=(m$x[length(m$x)])
               e=rnorm(h+1,0,sqrt(sigma.est))
                   for(i in 2:(h+1)){
                     ff[i,j]=ff[i-1,j]+drf.est+e[i-1]+ma1.est*e[i]
                                     }
              }
return(list(ff=ff))}
f.sim=ff.sim(f4,21,10000)$ff[-1,]
mean.sim=apply(f.sim,1,mean)
U.sim=apply(f.sim,1,quantile, probs = c(.95))
L.sim=apply(f.sim,1,quantile, probs = c(.05))
plot(1:21,mean.exact.f,type="l",ylim=range(mean.sim,U.sim,L.sim,mean.exact.f,U.exact.f,L.exact.f),xlab="t",col="blue",main="Simulated (green) & Exact (blue), forecasted mean and 95% C.I.")
points(1:21,U.exact.f,type="l",col="blue")
points(1:21,L.exact.f,type="l",col="blue")
points(1:21,mean.sim,type="l",col="green")
points(1:21,U.sim,type="l",col="green")
points(1:21,L.sim,type="l",col="green")

Answer 1

There are a couple of problems with your code:

1. The default prediction intervals are 80% and 95%, but you simulate 90% intervals.
2. The forecasts should use the last residual in the first forecast.
3. Your ma coefficient is applied to $e_t$ rather than $e_{t-1}$.

Here is some corrected code:

library(forecast)
f4 <- Arima(WWWusage, order=c(0,1,1), include.drift=TRUE)
f4.f <- forecast(f4, h=21, level=90)
mean.exact.f <- f4.f$mean
    U.exact.f <- f4.f$upper
L.exact.f <- f4.f$lower

ff.sim <- function(m,h,N)
{
  sigma.est <- m$sigma
      ma1.est <- coef(m)[1]
      drf.est <- coef(m)[2]
      ff <- matrix(0,ncol=N,nrow=(h+1))
      for(j in 1:N)
      {
        ff[1,j] <- m$x[length(m$x)]
        e <- c(m$residuals[length(m$x)],rnorm(h,0,sqrt(sigma.est)))
    for(i in 2:(h+1))
      ff[i,j] <- ff[i-1,j]+drf.est+ma1.est*e[i-1]+e[i]
  }
  return(list(ff=ff[-1,]))
}

f.sim <- ff.sim(f4,21,10000)$ff
mean.sim <-apply(f.sim,1,mean)
U.sim <-apply(f.sim,1,quantile, probs = c(.95))
L.sim <-apply(f.sim,1,quantile, probs = c(.05))
plot(1:21,mean.exact.f,type="l", 
  ylim=range(mean.sim,U.sim,L.sim,mean.exact.f,U.exact.f,L.exact.f), xlab="t",
  col="blue", main="Simulated (green) & Exact (blue), forecasted mean and 95% C.I.")
points(1:21,U.exact.f,type="l",col="blue")
points(1:21,L.exact.f,type="l",col="blue")
points(1:21,mean.sim,type="l",col="green")
points(1:21,U.sim,type="l",col="green")
points(1:21,L.sim,type="l",col="green")

You should get the blue and green lines almost identical now.