The nnetar function in the forecast package for R fits a neural network model to a time series with lagged values of the time series as inputs (and possibly some other exogenous inputs). So it is a nonlinear autogressive model, and it is not possible to analytically derive prediction intervals. Therefore we use simulation.
Suppose we fit a NNETAR model to the famous Canadian lynx data:
library(forecast)
set.seed(2015)
(fit <- nnetar(lynx, lambda=0.5))Series: lynx
Model: NNAR(8,4)
Call: nnetar(y = lynx, lambda = 0.5)
Average of 20 networks, each of which is
a 8-4-1 network with 41 weights
options were - linear output units
sigma^2 estimated as 95.77I’ve used a Box-Cox transformation with \lambda=0.5 to ensure the residuals will be roughly homoscedastic.
The model can be written as y_t = f(\boldsymbol{y}_{t-1}) + \varepsilon_t where \boldsymbol{y}_{t-1} = (y_{t-1},y_{t-2},\dots,y_{t-8})' is a vector containing lagged values of the series, and f is a neural network with 4 hidden nodes in a single layer.
The error series \{\varepsilon_t\} is assumed to be homoscedastic (and possibly also normally distributed).
We can simulate future sample paths of this model iteratively, by randomly generating a value for \varepsilon_t, either from a normal distribution, or by resampling from the historical values. So if \varepsilon^*_{T+1} is a random draw from the distribution of errors at time T+1, then y^*_{T+1} = f(\boldsymbol{y}_{T}) + \varepsilon^*_{T+1} is one possible draw from the forecast distribution for y_{T+1}. Setting \boldsymbol{y}_{T+1}^* = (y^*_{T+1}, y_{T}, \dots, y_{T-6})', we can then repeat the process to get y^*_{T+2} = f(\boldsymbol{y}^*_{T+1}) + \varepsilon^*_{T+2}. In this way, we can iteratively simulate a future sample path. By repeatedly simulating sample paths, we build up knowledge of the distribution for all future values based on the fitted neural network. Here is a simulation of 9 possible future sample paths for the lynx data. Each sample path covers the next 20 years after the observed data.
sim <- ts(matrix(0, nrow=20, ncol=9), start=end(lynx)[1]+1)
for(i in seq(9))
sim[,i] <- simulate(fit, nsim=20)
library(ggplot2)
autoplot(lynx) + forecast::autolayer(sim)
If we do this a few hundred or thousand times, we can get a very good picture of the forecast distributions. This is how the forecast.nnetar function produces prediction intervals:
fcast <- forecast(fit, PI=TRUE, h=20)
autoplot(fcast)
Because it is a little slow, PI=FALSE is the default, so prediction intervals are not computed unless requested. The npaths argument in forecast.nnetar controls how many simulations are done (default 1000). By default, the errors are drawn from a normal distribution. The bootstrap argument allows the errors to be “bootstrapped” (i.e., randomly drawn from the historical errors).