8
$\begingroup$

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)])
$\endgroup$
  • 1
    $\begingroup$ +1, very clever I see nothing wrong with the approach, s.e looks good, I'm not sure which s.e. You are referring? $\endgroup$ – forecaster Aug 31 '15 at 22:22
  • 1
    $\begingroup$ Ah forecaster! I missed you! :) I meant the SE of the MA and intercept between the two fits. The residual variance is identical. I think that the residual variance is all that matters for the prediction error of the forecast from an ARIMA model (although I am rusty). $\endgroup$ – B_Miner Aug 31 '15 at 22:44
  • $\begingroup$ Now I get it, I would be surprised if the s.e. is not close, being close is good, sigma^2 is used for prediction intervals, there is no difference between two models so prediction intervals will be same for arimax or Arima I believe. $\endgroup$ – forecaster Aug 31 '15 at 23:27
  • $\begingroup$ The way you define tf <- filter(...)... I am lost. Do you have any tips to understand it? What would happen if I had: I911-AR1: 0.55, I911-AR2: 0.66, I911-MA0: 0.49, I911-MA1: 0.39 ? $\endgroup$ – Konstantinos May 24 '17 at 1:52
  • 1
    $\begingroup$ Net stumblers, for the filtering, I found quite an insight here: onlinecourses.science.psu.edu/stat510/node/75 $\endgroup$ – Konstantinos May 24 '17 at 16:51
8
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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