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

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)
drf.est=coef(f4)
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") ## 1 Answer 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) drf.est <- coef(m) 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.

• Thank you so much Rob. I really appreciate your prompt reply.
– Stat
Sep 18, 2012 at 1:17