How does R calculate prediction intervals in the forecast package?

Question

I have a large dataset with different factors that I want to forecast to the future. These forecasts I will then later on use as inputs for a Monte Carlo simulation. My idea would be to use arima forecasting on the different variables. Subsequently, I would use the resulting prediction interval as inputs for the Monte Carlo simulation.

Using R, (I think) I get what I want by using the following.

First, I set some parameters FC_years <- 4 FC_boundaries <- 68 # This would be 1 sd

Next, the forecasting is done by: fit <- auto.arima(POP) forecast <- data.frame(forecast(fit,FC_years,level=FC_boundaries))

What is now very important for me in order to use these results for a MC simulation, is to know how R calculates the prediction interval (I have been searching for quite a while now, but I can't get a clear answer), and how this interval is distributed (I assume it is normal, but again, I can't get a clear answer).

Just as an example, the following dataset could be used:

1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 68235 72498 76700 80326 83195 85447 87276 89004 90858 92894 94995 97015 98742 100031 100830 101219 101344 101418 101597 101932 102384 102911

Can anybody give me any hints?

Answer 1

Unless you are using the forecast function's bootstrap = TRUE option the forecast package's ARIMA intervals are calculated by passing an ARIMA object to predict(). The intervals assume that residuals are normally distributed . If you look at the R documentation for predict.Arima() you will see that it uses KalmanForecast() to produce them.

See the following link for a good description of how to calculate a prediction interval for an ARIMA model manually: https://onlinecourses.science.psu.edu/stat510/node/66

The code is available on Rob Hyndman's github page here: https://github.com/robjhyndman/forecast/blob/master/R/arima.R

An exception to the resulting prediction intervals being normal around the point estimate is if you have set lambda (are using a Box-Cox transformation) in which case they will not be normal after the back-transformation.

Answer 2

https://www.otexts.org/fpp/2/7 Goes into some further detail.

95% prediction interval is $y^t \pm 1.96 \sigma$

and

80% prediction interval is $y^t \pm 1.28 \sigma$

Where $y^t$ is a forecast value and $\sigma$ "is an estimate of the standard deviation of the forecast distribution".