Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Forecasting error in R when passing around arguments in forecast() and ar()

When trying to compose a function from smaller ones using Rob Hyndman's forecast library, like so:

> library('forecast')
> arf <- function(data, ...) forecast(ar(data, order.max=1, method="ols"), ...)

I get an error when trying to plug in some data:

> arf(ts(1:100, start=c(2000,1), frequency=4))
Error in ts(x, frequency = 1, start = 1) : object is not a matrix

However, using the body of arf directly works perfectly:

> forecast(ar(ts(1:100, start=c(2000,1), frequency=4), order.max=1,method="ols"))
        Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2025 Q1            101   101   101   101   101
2025 Q2            102   102   102   102   102
2025 Q3            103   103   103   103   103
2025 Q4            104   104   104   104   104
2026 Q1            105   105   105   105   105
2026 Q2            106   106   106   106   106
2026 Q3            107   107   107   107   107
2026 Q4            108   108   108   108   108
2027 Q1            109   109   109   109   109
2027 Q2            110   110   110   110   110

Why does arf not work as it should?

like image 933
CrystalCuckoo Avatar asked Sep 14 '14 11:09

CrystalCuckoo


People also ask

What is the forecast function in R?

forecast is a generic function for forecasting from time series or time series models. The function invokes particular methods which depend on the class of the first argument. For example, the function forecast. Arima makes forecasts based on the results produced by arima .

What does auto Arima do in R?

The auto. arima() function in R uses a variation of the Hyndman-Khandakar algorithm (Hyndman & Khandakar, 2008), which combines unit root tests, minimisation of the AICc and MLE to obtain an ARIMA model. The arguments to auto. arima() provide for many variations on the algorithm.

Which package is auto Arima in?

In this case, auto. arima from the forecast package in R allows us to implement a model of this type with relative ease.


2 Answers

This is a problem (not really a bug) in forecast.ar(). All forecast.xxx() functions try to store the data used to estimate the time series model as it is required for plots and accuracy calculations. However, ar() does not return the data, so forecast.ar() attempts to find the data in the calling environment, or in the parent environment. When you call forecast(ar(...)) directly, the function manages to find the data, but arf() places the call to ar() one level deeper making it that much harder for forecast to figure out what data was being used.

I can modify the function to make it look harder for the data (i.e., look in the grandparent environment also), but the construction will still fail because predict.ar() (part of the stats package) will then cause a similar error. It would be much better if ar() returned the data, but ar() is part of the stats package and I have no control over it.

There are several possible solutions.

  1. You could replace ar with Arima:

    arf <- function(dat, ...) forecast(Arima(dat, order=c(1,0,0)), ...)
    

    That should return the same model if the data are stationary (although the parameter estimates will be slightly different). It won't return the same answer for the example in your question because the time series is non-stationary.

  2. You could use auto.arima() instead if you were willing to use more general ARIMA models than AR(1).

    arf <- function(dat, ...) forecast(auto.arima(dat, ...)
    
  3. (Based on a suggestion from @agstudy). A workaround is to ensure the data is stored within the ar object:

    arf <- function(dat, ...) 
    {
      object <- ar(dat, order.max=1, method="ols")
      object$x <- dat
      forecast(object,...)
    }
    
like image 176
Rob Hyndman Avatar answered Oct 17 '22 01:10

Rob Hyndman


the problem is a bug in the S3 method predict for ar class. predict.ar try to evaluate the newdata parameter using ar object. Showing the first lines of getS3method('predict','ar')

function (object, newdata, n.ahead = 1L, se.fit = TRUE, ...) 
{
    if (n.ahead < 1L) 
        stop("'n.ahead' must be at least 1")
    if (missing(newdata)) {
        newdata <- eval.parent(parse(text = object$series))
        if (!is.null(nas <- object$call$na.action)) 
            newdata <- eval.parent(call(nas, newdata))
    }

  .....
}

the relevant/bugged line is:

 newdata <- eval.parent(parse(text = object$series))

But object$series don't have the right expression/character. since it is hidden by the new level of wrapper function arf. Here a workaround, is to set the right expression for this term:

arf <- function(dat, ...) 
  {
    object <- ar(dat, order.max=1, method="ols")
    object$series <- as.character(as.expression(as.list(match.call())$dat))
    forecast(object,...)
}
arf( ts(1:100, start=c(2000,1), frequency=4)

Note that ; this solution works also with:

aa <- ts(1:100, start=c(2000,1), frequency=4)
arf(aa)
like image 41
agstudy Avatar answered Oct 17 '22 03:10

agstudy