Arimax Prediction – Using Forecast Package for Arimax Prediction

arimaforecastingintervention-analysisrtime series

The arimax function in the TSA package is to my knowledge the only R package that will fit a transfer function for intervention models. It lacks a predict function though which is sometimes needed.

Is the following a work-around for this issue, leveraging the excellent forecast package? Will the predictive intervals be correct? In my example, the std errors are "close" for the components.

  1. Use the forecast package arima function to determine the pre-intervention noise series and add any outlier adjustment.
  2. Fit the same model in arimax but add the transfer function
  3. Take the fitted values for the transfer function (coefficients from arimax) and add them as xreg in arima.
  4. Forecast with arima
library(TSA)
library(forecast)
data(airmiles)
air.m1<-arimax(log(airmiles),order=c(0,0,1),
              xtransf=data.frame(I911=1*(seq(airmiles)==69)),
              transfer=list(c(1,0))
              )

air.m1

Output:

Coefficients:
  ma1  intercept  I911-AR1  I911-MA0
0.5197    17.5172    0.5521   -0.4937
s.e.  0.0798     0.0165    0.2273    0.1103

sigma^2 estimated as 0.01223:  log likelihood=88.33
AIC=-168.65   AICc=-168.09   BIC=-155.02

This is the filter, extended out 5 more periods that the data

tf<-filter(1*(seq(1:(length(airmiles)+5))==69),filter=0.5521330,method='recursive',side=1)*(-0.4936508)
forecast.arima<-Arima(log(airmiles),order=c(0,0,1),xreg=tf[1:(length(tf)-5)])
forecast.arima

Output:

Coefficients:
         ma1  intercept  tf[1:(length(tf) - 5)]
      0.5197    17.5173                  1.0000
s.e.  0.0792     0.0159                  0.2183

sigma^2 estimated as 0.01223:  log likelihood=88.33
AIC=-168.65   AICc=-168.28   BIC=-157.74

Then to Predict

predict(forecast.arima,n.ahead = 5, newxreg=tf[114:length(tf)])

Best Answer

Prediction intervals are based on residual variance as estimated by the maximum likelihood optimization.

I like forecast function from the forecast package:

fc <- forecast(forecast.arima,h = 5, xreg=tf[114:length(tf)])

which gives me the following predictions: 

  Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jun 2005       17.61393 17.47222 17.75565 17.39720 17.83067
Jul 2005       17.51725 17.35753 17.67697 17.27299 17.76152
Aug 2005       17.51725 17.35753 17.67697 17.27299 17.76152
Sep 2005       17.51725 17.35753 17.67697 17.27299 17.76152
Oct 2005       17.51725 17.35753 17.67697 17.27299 17.76152

the way prediction interval works is by following formula:

$y_t \pm z\sigma$ where $z$ is a multiplier which takes values such as 1.96 for 95% prediction interval and 1.28 for 80% prediction interval $\sigma$ is the standard deviation of the residual also forecast distribution which is also the square root of sigma^2 in the maximum likelihood estimate. As you show sigma^2 ( 0.01223) is identical in both your models, prediction intervals will also be same and will match.

If you want to check it,

Upper Limits:

> fc$mean[1]+sqrt(forecast.arima$sigma2)*1.96
[1] 17.83068

Lower Limits:

> fc$mean[1]-sqrt(forecast.arima$sigma2)*1.96
[1] 17.39719

which matches the prediction interval provided by the forecast function. To answer your question, yes leveraging forecast package will work in this case and the prediction intervals will be correct.

Related Question