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?
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 .
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.
In this case, auto. arima from the forecast package in R allows us to implement a model of this type with relative ease.
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.
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.
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, ...)
(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,...)
}
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)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With