1
$\begingroup$

I am using an ARIMA model to forecast a few periods ahead of my data in R. To do this I begin by estimating the ARIMA model using the forecast::auto.arima function and then compute predictions using the forecast::forecast function. How can I obtain the standard errors of these predictions from the forecast function output?

Here is a simple example to illustrate my problem. I begin by loading the AirPassengers data that comes with R and transforming it into a ts object:

data(AirPassengers)
air.ts <- ts(AirPassengers, freq=12)

Then I estimate an ARIMA model:

air.arima <- auto.arima(air.ts)

I then forecast three months ahead of the data, which ends in December 1960:

air.fc <- forecast(air.arima, h=3)

> air.fc
       Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jan 13       445.6349 430.8903 460.3795 423.0851 468.1847
Feb 13       420.3950 403.0907 437.6993 393.9304 446.8596
Mar 13       449.1983 429.7726 468.6241 419.4892 478.9074

As you can see, the output includes point forecasts as well as 80% and 95% confidence intervals, but not standard errors.

I could obtain these standard errors using the stats::predict function:

> predict(air.arima, n.ahead=3, se.fit=TRUE)
$pred
        Jan      Feb      Mar
13 445.6349 420.3950 449.1983

$se
        Jan      Feb      Mar
13 11.50524 13.50261 15.15799

However, the output of the function predict is not compatible with the forecast::accuracy function which I need to use to compute the Mean Absolute Scaled Errors of my predictions, and this is why I would like to obtain the standard errors from the forecast function instead.

Added in response to comments: Just to clarify, I don't intend to use the standard errors of the predictions as inputs into the accuracy function or any other analysis, I am simply interested in them as a measure of the uncertainty associated with those predictions.

At the moment I need to use both the predict and the forecast function, the former to obtain the standard errors (one of the outputs of my analysis) and the latter as an input into the accuracy function in order to compute the MASE (a second, separate output).

Improve this question
$\endgroup$
8
  • 3
    $\begingroup$ The way I understood your question is: I need to get standard errors -> so I can give them to the function accuracy as an argument -> so I can get MASE from the accuracy function. My question is: why do you want the standard errors by themselves if the function accuracy actually takes the whole forecast object as an argument? $\endgroup$
    – suckrates
    Jul 26, 2018 at 11:30
  • $\begingroup$ +1 to @AlvaroFuentes. You can simply do accuracy(forecast(air.arima, h=3)) to get your MASE. In addition, the MASE does not have anything to do with the standard errors, it only looks at your point predictions. Why do you want SEs? $\endgroup$
    – Stephan Kolassa
    Jul 26, 2018 at 11:40
  • $\begingroup$ Thank you for your comments and apologies if I was unclear. The reason I wanted the SEs is that I am interested not only in the model's forecasting performance (measured through MASE) but also in the forecasts themselves. Since I will also be looking at the forecasted values, it's useful to have a measure of the uncertainty associated with them. I know I can get confidence intervals from the forecast output but SEs would be a more compact measure. $\endgroup$
    – Lino Ferreira
    Jul 26, 2018 at 11:54
  • 1
    $\begingroup$ OK, but in this case, predict.Arima() does what you want, doesn't it? Is your question already answered? (Note that forecast() does out output confidence intervals but prediction intervals; there is a difference.) $\endgroup$
    – Stephan Kolassa
    Jul 26, 2018 at 12:20
  • 1
    $\begingroup$ You can always define your own function wrapper that calls both predict.Arima and accuracy if you want to make your code simpler. I doubt that the function calls take up a lot of time, especially compared to auto.arima. $\endgroup$
    – Stephan Kolassa
    Jul 26, 2018 at 13:41

1 Answer 1

Reset to default
1
$\begingroup$

The prediction intervals are based on a normal distribution, so if you want to find the standard error of the forecast distribution, you can just compute it from the intervals.

If fc is a forecast object, then

fsd <- (fc$upper[,1] - fc$lower[,1]) / (2 * qnorm(.5 + fc$level[1] / 200))

gives you the forecast standard deviations for all horizons. In your example:

library(forecast)
air.arima <- auto.arima(AirPassengers)
air.fc <- forecast(air.arima, h=3)
predict(air.arima, n.ahead=3, se.fit=TRUE)
#> $pred
#>           Jan      Feb      Mar
#> 1961 445.6349 420.3950 449.1983
#> 
#> $se
#>           Jan      Feb      Mar
#> 1961 11.50524 13.50261 15.15799
(air.fc$upper[,1] - air.fc$lower[,1]) / (2 * qnorm(.5 + air.fc$level[1] / 200))
#>           Jan      Feb      Mar
#> 1961 11.50524 13.50261 15.15799

By the way, this is not a good model as the data should be logged first.

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.