forecast package v4.0

Published on 3 December 2012

A few days ago I released version 4.0 of the forecast package for R. There were quite a few changes and new features, so I thought it deserved a new version number. I keep a list of changes in the Changelog for the package, but I doubt that many people look at it. So for the record, here are the most important changes to the forecast package made since v3.0 was released.

Exponential smoothing

The ets() function handles exponential smoothing in the ETS framework, but for teaching purposes it is sometimes useful to discuss a more traditional implementation that does not involve optimizing the initial values. The HoltWinters() function in the stats package attempts to do this, although the initialization used for the Holt and Holt-Winters methods is non-standard, making it unsuitable for teaching.

Consequently, I’ve now modified the ses(), holt() and hw() functions to allow both traditional heuristic initialization and ETS-style optimization of initial states. These functions also now allow damped trend and multiplicative trend.

Multiple seasonality and BATS and TBATS models

Time series with multiple seasonal patterns (e.g., hourly data that contains a daily pattern, weekly pattern and an annual pattern) now have their own model class msts. This is useful when plotting the data, and in using the dshw() double seasonal Holt-Winters function for example.

The BATS and TBATS models (fitted using bats() and tbats()) also handle multiple seasonalities as well as complex seasonalities (such as non-integer seasonality, non-nested seasonality and large period seasonality). The models are in the exponential-smoothing framework and were proposed by De Livera, Hyndman and Snyder (JASA 2011).

BATS and TBATS models allow a form of time series decomposition, and so the fitted objects can now plotted in the same way as other decomposition objects (e.g., stl objects).

These models are quite computationally intensive, and so the functions will use parallel processing when it is available.

auto.arima()

This is probably the most popular function in the package, and I am always looking at ways to speed it up and improve the results.

Default model selection has been switched from AIC to AICc (the bias-corrected form of AIC) which may affect model selection for very short time series.
The maximum orders are now restricted to less than ¹⁄₃ of the length of the data.
Parallel processing is used where possible when stepwise=FALSE.

Accuracy calculations

The accuracy function now works with both time series forecasts and cross-sectional forecasts.

The MASE for seasonal time series is now defined using seasonal naïve forecasts for the in-sample scaling. The MASE for cross-sectional forecasts is defined using the mean forecasts for in-sample scaling. This is consistent with my new book.

Neural network AR models

A new experimental function nnetar has been introduced to fit a neural network with lagged values of the time series as inputs. It is a feed-forward network with one hidden layer, specified using the notation NNAR( $p,k$ ) to indicate there are $p$ lagged inputs and $k$ nodes in the hidden layer. For example, a NNAR(9,5) model is a neural network with the last nine observations $(y_{t-1},y_{t-2},\dots,y_{t-9}$ ) used as inputs to forecast the output $y_t$ , and with five neurons in the hidden layer. A NNAR( $p,0$ ) model is equivalent to an ARIMA( $p,0,0$ ) model but without the restrictions on the parameters to ensure stationarity.

With seasonal data, it is useful to also add the last observed values from the same season as inputs. For example, an NNAR(3,1,2) $_{12}$ model has inputs $y_{t-1}$ , $y_{t-2}$ , $y_{t-3}$ and $y_{t-12}$ , and two neurons in the hidden layer. More generally, an NNAR( $p,P,k$ ) $_m$ model has inputs $(y_{t-1},y_{t-2},\dots,y_{t-p},y_{t-m},y_{t-2m},y_{t-Pm})$ and $k$ neurons in the hidden layer. A NNAR( $p,P,0$ ) $_m$ model is equivalent to an ARIMA( $p,0,0$ )( $P,0,0$ ) $_m$ model but without the restrictions on the parameters to ensure stationarity.

The nnetar() function fits an NNAR( $p,P,k$ ) $_m$ model. If the values of $p$ and $P$ are not specified, they are automatically selected. For non-seasonal time series, the default is the optimal number of lags (according to the AIC) for a linear AR( $p$ ) model. For seasonal time series, the default values are $P=1$ and $p$ is chosen from the optimal linear model fitted to the seasonally adjusted data. If $k$ is not specified, it is set to $k=(p+P+1)/2$ (rounded to the nearest integer).

By default, 25 networks with random starting values are trained and their predictions averaged.

Like many of the functions in the forecast package, the default arguments usually give reasonably good forecasts without the user needing to know much. For example:

fit <- nnetar(lynx)
plot(forecast(fit,h=20))

Because random starting values are used, the forecasts will be slightly different each time the function is called.

It is possibly to compute prediction intervals from NNAR models via simulation, but I have not yet implemented this. The function will only give point forecasts for now.

New default colors for forecast plots

The old yellow and orange colors have been replaced by shades of blue-grey which are a little easier on the eye. The darker the color, the higher the density. Here is an example with 80% and 95% prediction intervals.

library(fpp)
fit <- Arima(euretail, order=c(0,1,3), seasonal=c(0,1,1))
plot(forecast(fit,h=12))

Future development

Earlier this year I started using the package for an undergraduate forecasting class which gave it a good work-out, and I’m also using it as the basis of my new book with George Athanasopoulos. As a result, a lot of new functionality has been added, and lots of bugs have been discovered and fixed.

The process is continuing and I have lots of ideas I want to implement, and no doubt more bugs will be discovered. The development of the forecast package is now maintained on github where you can see what is in development and report any new problems. Please let me know if you find anything you think is an error, or if you have a suggestion.

Tags: computing, forecasting, R, statistics

22 Comments

fernando

It would be nice to have an automatic cross-validation function, given the fitted model and the validation parameters (training-set size, max Horizon, etc). Anyway, great job in the new version!
- http://www.facebook.com/zach.a.mayer Zachary Mayer
  
  There’s some code here you can try out: http://robjhyndman.com/researchtips/tscvexample/
  
  You can also take a look at the R package I am developing to do this: https://github.com/zachmayer/cv.ts
Samik

Hi Rob, thanks for the notification. Indeed, the ‘What’s new’ section rarely gets read. Anways, I was wondering if you have an updated version of the original paper which detailed the auto.arima() algorithm, that includes the improvements and changes?
- http://robjhyndman.com Rob J Hyndman
  
  The description at http://otexts.com/fpp/8/ is up-to-date.
http://www.logisticadescomplicada.com/ Leandro Callegari Coelho

Rob, congrats for the development. I like this package a lot. By the way, whenever a new version is available, the first thing I look for is the changelog/what’s new section!
Shea Parkes

Thanks for all your hard work Rob! I appreciate the subtler default colors, but you’ll need to go back and update your online book to stop mentioning yellow/orange now: http://otexts.com/fpp/1/4/
Jerome

Thank you Rob for the function dshw. When I get the point of forecast, how could I get the forecast interval (low 0.80 and upper 0.80)? Thank you
- http://robjhyndman.com Rob J Hyndman
  
  I have not implemented prediction intervals for dshw(). Use bats() to get PI for an equivalent model.
BMiner

Hey Rob, You mention prediction intervals with ANN via simulation. Do you have any sources to reference?
- http://robjhyndman.com Rob J Hyndman
  
  No, but a neural network is essentially a nonlinear autoregression and simulation is the most common way to get PI for other nonlinear autoregression models.
  - Bminer
    
    Do you think bootstrapping would be a viable method (using the boot package and tsboot)? I am thinking of bootstrapping the series and then using each of the n models created to predict a new case — and thus getting a prediction distribution.
    - http://robjhyndman.com Rob J Hyndman
      
      No. The bootstrapped series won’t allow prediction conditional on the last few values of the observed series.
Joseph Johnson

Dr. Hyndman,

Does the ets function in forecast package v4.0 constrain the parameters to be;
0<α<1, 0<β<α, 0<γ<1-α ?
Would you recommend this constraint for all models?
Does the package go even further to constraint the certain models as mentioned in Chapter 10 of your exponential smoothing book?

Thank you for the package. It has been very helpful to me!
- http://robjhyndman.com Rob J Hyndman
  
  Try reading the help file. It is all explained. The constraints you mention are implemented using bounds=“usual”. The constraints in ch10 of our book are implemented using bounds=“admissible”. The default is bounds=“both” meaning both sets of constraints are applied.
- Joseph Johnson
  
  Thank you.
Francesco Gianferrari Pini

Hi! How can I extract the time series for the seasonal component(s) on a TBATS forecast? I can see them in the plot!
- http://robjhyndman.com Rob J Hyndman
  
  Look at plot.tbats() and you will see how they are constructed.
Jerome

Bonjour Rob! I am running a bats function and I get this error:

Error in checkForRemoteErrors(val) : 6 nodes produced errors; first error: argument is of length zero
I have no clue to solve it. Thank you in advance for your tips.
Jerome
- http://robjhyndman.com Rob J Hyndman
  
  Try using debug() to track where the error is occurring.
SimonMN

Hi Prof Hyndman,

Is there any way to add exogenous inputs to the nnetar function in the same way you can add exogenous variables to the Arima or auto.arima functions?
- http://robjhyndman.com/ Rob J Hyndman
  
  No, but it isn’t hard to do it using the avNNet function from the caret package. Basically nnetar() is forecasting using an average of 20 neural nets where the inputs are lagged values of the time series. You would just need to add the xreg variables as inputs as well.
Christian Miller

I don’t have anything substantive to add, but I want to give a huge thanks to you Rob. This has been a helpful and useful package in learning about forecasting and implementing ideas. And you continue to add functionality faster than I can keep up: it’s a great problem to have.

Greatly appreciated (and the book too!)

Hyndsight

A blog by Rob J Hyndman