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")