Major changes to the forecast package

The forecast package for R has undergone a major upgrade, and I’ve given it version number 3 as a result. Some of these changes were suggestions from the forecasting workshop I ran in Switzerland a couple of months ago, and some have been on the drawing board for a long time. Here are the main changes in version 3, plus a few earlier additions that I thought deserved a mention.

Box-Cox transformations

A Box-Cox transformation can now be incorporated into the forecasting model when calling Arima(), auto.arima(), ets(), arfima(), tslm(), stlf(), rwf(), meanf(), splinef(). For example:

fit <- Arima(lynx, order=c(2,0,0), lambda=0.5)

The model is based on the transformed data, forecasts are calculated, and then the forecasts and prediction intervals are back-transformed. The point forecasts can be interpreted as medians after back-transforming.

If the transformation is done outside the fitting function, the forecasts can still be back-transformed. For example:

fit <- ar(BoxCox(lynx,0.5))

Back-transforming forecasts like this is now available in forecast.Arima(), forecast.ets(), forecast.fracdiff(),, forecast.StructTS, forecast.HoltWinters().

I have also added a function for automatically choosing the Box-Cox parameter using either Guerrero’s (1993) method or the profile log likelihood method. For example:

fit <- Arima(lynx, order=c(2,0,0), lambda=BoxCox.lambda(lynx))

Note that previously there was a lambda argument in the plot.forecast() function. This is no longer available (and so some old code may break). Instead, back-transform the forecasts within the forecast() function.

Improved auto.arima()

The auto.arima() function is widely used for automatically selecting ARIMA models. It works quite well, except that selection of D, the order of seasonal differencing, has always been poor. Up until now, the default has been to use the Canova-Hansen test to select D. Because the CH test has a null hypothesis of deterministic seasonality based on dummy variables, the function will often select D=0. So I’ve now switched to using the OCSB test for selecting D which has a null hypothesis involving a seasonal difference, so it is much more likely to choose D=1 than previously. I’ve done extensive testing of the forecasts obtained under the two methods, and the OCSB test leads to better forecasts. Hence it is now the default. This means that the function may return a different ARIMA model than previously when the data are seasonal.

A separate function for selecting the seasonal order has also been made visible. So you can now call nsdiffs() to find the recommended number of seasonal differences without calling auto.arima(). There is also a ndiffs() function for selecting the number of first differences. Within auto.arima(), nsdiffs() is called first to select D, and then ndiffs() is applied to diff(x,D) if D>0 or to x if D=0.

Double-seasonal Holt-Winters

The new dshw() function implements Taylor’s (2003) double-seasonal Holt-Winters method. This allows for two levels of seasonality. For example, with hourly data, there is often a daily period of 24 and a weekly period of 168. These are modelled separately in the dshw() function.

I am planning some major new functionality to extend this to the various types of complex seasonality discussed in my recent JASA paper. Hopefully that will be ready in the next few weeks — I have a research assistant working on the new code.

Sub-setting time series

Occasionally you want to extract all the Novembers in a monthly time series (or something similar), but this has been fiddly to do up to now. So I’ve included a new function subset.ts. For example:


Acf() and Pacf()

These were actually added in v2.19 but I’ve not mentioned them anywhere so I thought it would be useful to say something here. The acf() function always includes a spike of length 1 at lag 0. This is pointless because the ACF at lag 0 is 1 by definition. It is also annoying because it forces the scale of the y-axis to include 1 which can obscure smaller correlations that might be of interest. The Acf() function works in the same way as the acf() function except that it omits lag 0. The Pacf() function is included for consistency only — it returns the same object and produces the same plot as pacf().

Time series linear models

Another recent addition, but not new in v3, is the tslm() function for handling linear models for time series. This works in the same way as lm() except that the time series characteristics of the data are preserved in the residuals and fitted values. Also, the variables trend and seasonal can be used without needing to be defined. For example:

y <- ts(rnorm(120,0,3) + 20*sin(2*pi*(1:120)/12), frequency=12)
fit1 <- tslm(y ~ trend + season)

trend takes values 1,2,\dots,n where n is the length of y, and season is a matrix of seasonal dummy variables.


The CV() function is another addition from earlier in the year. It implements the cross-validation statistic, AIC, corrected AIC, BIC and adjusted R2 values for a linear model. For example:

y <- ts(rnorm(120,0,3) + 20*sin(2*pi*(1:120)/12), frequency=12)
fit1 <- tslm(y ~ trend + season)
fit2 <- tslm(y ~ season)

CV works with any lm objects including those produced by tslm() and lm().

Forecasting with STL

This functionality was added earlier in the year, but it is so cool I wanted to mention it here. STL is a great method for decomposing time series into trend, seasonal and irregular components. It is robust and handles time series of any frequency.

To forecast with STL, you seasonally adjust the data by subtracting the seasonal component, then forecast the seasonally adjusted data using a non-seasonal ARIMA or ETS model, then re-seasonalize the forecasts by adding back in the most recent values of the seasonal component (effectively using a seasonal naive forecast for the seasonal component). The whole procedure is handled effortlessly as follows

fit <- stl(USAccDeaths,s.window="periodic")

There is also a new function stlf() which does the STL decomposition as well as the forecasting in one step.

plot(stlf(AirPassengers, lambda=BoxCox.lambda(AirPassengers)))

STL decompositions are always additive, but the inclusion of the Box-Cox parameter as shown here allows non-additive decompositions as well.

For data with high seasonal period (such as weekly data, hourly data, etc.), forecasting with STL is often the simplest approach. It also works amazingly well on a wide range of series. If you apply the forecast() function to a time series (rather than a time series model), the forecasts returned will be from stlf() when the seasonal period is 13 or more.

Other changes

A list of all changes to the forecast package is maintained in the ChangeLog.

With so many changes and new functions, I’ve probably introduced new bugs. Please let me know if you find any problems and I’ll endeavour to fix anything ASAP. (Make sure you get v3.02 or later as v3.00 had some bugs.)

Related Posts:

  • giorgio

    Hi, could you provide us a reference to the algorithm you adopted to implement the OCSB test? Thanks

    • Rob J Hyndman

      See the help files. We give the following references in the help for nsdiffs():

      Osborn DR, Chui APL, Smith J, and Birchenhall CR (1988) “Seasonality and the order of integration for consumption”,Oxford Bulletin of Economics and Statistics 50(4):361-377.

      Osborn, D.R. (1990) “Seasonality and the order of integration in consumption”, International Journal of Forecasting,6:327-336.

      • Giorgio

        Thank you. Best, Giorgio

  • Pingback: Research tips - forecast package v4.0()

  • afsdafsdfsdadfsafa

    Hello, many thanks for your effort. I have a question about ocsb and ch tests. How is it possible to get p values and statistics? What is the reason that these two tests give different answers?

    • Rob J Hyndman

      They have different test statistics and different null hypthoses. Why should they give the same answer? You can get the details within the functions forecast:::CHtest and forecast:::OCSBtest.

      • afsdafsdfsdadfsafa

        Thank you for your answer.
        Another question of mine is about test statistics. Is it anyhow possible to get their values and/or p values?

  • Joe Park

    Thank you for making your book open access. Without any background in statistics, I read the book in four days and I am now making forecasts.

    Could you please tell me why auto.arima() is giving a non-seasonal arima for data that seems to clearly have a seasonality of 12 in my monthly data as demonstrated in the attached picture? Also, is there a better place for me to find answers to questions of this sort without using up your time?
    Thank you,

    • Rob J Hyndman

      You have not specified the frequency attribute in the ts object. Please ask questions on

      • Joe Park

        Thank you

  • Pingback: R Forecasting with STL()

  • Pingback: New in forecast 4.0 | Hyndsight()

  • John

    Hi! I am playing with the lambda parameter for the BoxCox transformation using lynx data! I noticed that if I used a data frame I get different lambda value than using the corresponding “ts” object. For example:
    BoxCox.lambda(lynx) = 0.1521849
    BoxCox.lambda( = 1.
    I can’t figure what I’m doing wrong. Can you help me?
    Thanks [John]

    • Rob J Hyndman

      The function requires a ts object or numerical vector, as specified in the help file. If you give it something else, don’t expect the results to make sense.

      • John

        You are right! I considered that if I use different objects the function could respond the same as with “ts” and numerical vectors. Thank you for your immediate response! [John]