Transfer Function Models :: Arimax in TSA

Question

I am trying to replicate the model from forecasting with dynamic regression models and I can't match the output in R using the arimax function from the TSA library. I am able to get pretty close to the result with SAS but I want to use R and hoping someone knows how to code the arimax function to achieve this. I have found the function to have issues (normally rooted from arima and optim) converging properly, but in this case a model is returned, but the parameters are way off.

The data is the first 140 observations from the sale sand lead series from Box and Jenkins in the astsa library.

Here is the snippet from the book showing their results (again, I could get close with SAS) and the code used with R (and the results). One thing I note is in the help file for arimax() there is recommendation to "mean-delete" transfer function covariates. I am not sure what this means and not sure if that is part of the issue.

From the book:

and here is the R-code:

library(TSA)
library(Hmisc)
library(astsa)


sales_140<-window(sales,end=140)
lead_140<-window(lead,end=140)

mod<-arimax(window(sales_140,start=4),order=c(0,1,1),
       xtransf = window(Lag(lead_140,3),start=4),transfer = list(c(1,0)),
       xreg=data.frame(seq(1:137)),method="ML")
mod


#Series: window(sales_140, start = 4) 
#ARIMA(0,1,1)                    

#Coefficients:
#         ma1  seq.1.137.  T1-AR1  T1-MA0
#      0.5974      0.3322  0.0613  2.8910
#s.e.  0.0593      0.1111  0.0275  0.1541

#sigma^2 estimated as 0.6503:  log likelihood=-163.94
#AIC=335.87   AICc=336.34   BIC=350.44

Here is the SAS code:

proc arima data=BL;
identify var=sales(1) crosscorr=lead(1);
estimate q=1 input=( 3 $ ( 0 ) / ( 1) lead) method=ml;
forecast out = out1 lead = 0;

run; 

and the estimates:

Answer 1

ARIMAX models can be a bit difficult to implement/interpret in R. In this case there are a few things that tripped you up. Here they are in no particular order:

"mean-delete" is another way of saying "remove the mean". In this case, it refers to the covariate lead_140. So, start with

lead_140_Z <- lead_140 - mean(lead_140).

The order of the ARIMAX model you are trying to fit is (0,1,1), which is the same as ARMAX(0,1) on the first-differenced data. So, rather than work with the differencing inside of the model, just do so beforehand:

sales_140_D <- diff(sales_140)
lead_140_D <- diff(lead_140_Z)

In this case, the order of the transfer function is actually (1,3), but the first, second and third MA parameters (MA0, MA1 and MA2) are fixed at 0 (ie, only B^3 appears in the numerator). To address this, you need to use the fixed argument in ARIMAX() and specify NA for those params to estimate and 0 for those to omit.

You don't need anything for xreg as the covariate occurs in the transfer.

mod <- arimax(sales_140_D,
   order=c(0,0,1),
   include.mean=TRUE,
   fixed=c(NA,NA,NA,0,0,0,NA),
   xtransf=lead_140_D,
   transfer=list(c(1,3)),
   method="ML")

mod

# Coefficients:
          ma1  intercept  T1-AR1  T1-MA0  T1-MA1  T1-MA2  T1-MA3
      -0.5791     0.0286  0.7255       0       0       0  4.7092
s.e.   0.0756     0.0090  0.0040       0       0       0  0.0551

The results aren't exact, but they are quite close.