**R on Rob J Hyndman**, and kindly contributed to R-bloggers)

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)
(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 98.11
```

I’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) + 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).

**leave a comment**for the author, please follow the link and comment on their blog:

**R on Rob J Hyndman**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...