7
$\begingroup$

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

Plot

$\endgroup$
5
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you so much Rob. I really appreciate your prompt reply. $\endgroup$ – Stat Sep 18 '12 at 1:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.