## The ARIMAX model muddle

Published on 4 October 2010
There is often con­fu­sion about how to include covari­ates in ARIMA mod­els, and the pre­sen­ta­tion of the sub­ject in var­i­ous text­books and in R help files has not helped the con­fu­sion. So I thought I’d give my take on the issue. To keep it sim­ple, I will only describe non-​​seasonal ARIMA mod­els although the ideas are eas­ily extended to include sea­sonal terms. I will include only one covari­ate in the mod­els although it is easy to extend the results to mul­ti­ple covari­ates. And, to start with, I will  assume the data are sta­tion­ary, so we only con­sider ARMA models.
Let the time series be denoted by . First, we will define an ARMA model with no covariates:

where is a white noise process (i.e., zero mean and iid).

### ARMAX mod­els

An ARMAX model sim­ply adds in the covari­ate on the right hand side:

where is a covari­ate at time and is its coef­fi­cient. While this looks straight-​​forward, one dis­ad­van­tage is that the covari­ate coef­fi­cient is hard to inter­pret. The value of is not the effect on when the is increased by one (as it is in regres­sion). The pres­ence of lagged val­ues of the response vari­able on the right hand side of the equa­tion mean that can only be inter­preted con­di­tional on the value of pre­vi­ous val­ues of the response vari­able, which is hardly intuitive.
If we write the model using back­shift oper­a­tors, the ARMAX model is given by

where and . Notice how the AR coef­fi­cients get mixed up with both the covari­ates and the error term.

### Regres­sion with ARMA errors

For this rea­son, I pre­fer to use regres­sion mod­els with ARMA errors, defined as follows.

In this case, the regres­sion coef­fi­cient has its usual inter­pre­ta­tion. There is not much to choose between the mod­els in terms of fore­cast­ing abil­ity, but the addi­tional ease of inter­pre­ta­tion in the sec­ond one makes it attractive.
Using back­shift oper­a­tors, this model can be writ­ten as

### Trans­fer func­tion models

Both of these mod­els can be con­sid­ered as spe­cial cases of trans­fer func­tion mod­els, pop­u­lar­ized by Box and Jenkins:

This allows for lagged effects of covari­ates (via the oper­a­tor) and for decay­ing effects of covari­ates (via the operator).
Some­times these are called “dynamic regres­sion mod­els”, although dif­fer­ent books use that term for dif­fer­ent models.
The method for select­ing the orders of a trans­fer func­tion model that is described in Box and Jenk­ins is cum­ber­some and dif­fi­cult, but con­tin­ues to be described in text­books. A much bet­ter pro­ce­dure is given in Pankratz (1991), and repeated in my 1998 fore­cast­ing text­book.

### Non-​​stationary data

For ARIMA errors, we sim­ply replace with where denotes the dif­fer­enc­ing oper­a­tor. Notice that this is equiv­a­lent to dif­fer­enc­ing both and before fit­ting the model with ARMA errors. In fact, it is nec­es­sary to dif­fer­ence all vari­ables first as esti­ma­tion of a model with non-​​stationary errors is not con­sis­tent and can lead to “spu­ri­ous regression”.

### R func­tions

The arima() func­tion in R (and Arima() and auto.arima() from the fore­cast pack­age) fits a regres­sion with ARIMA errors. Note that R reverses the signs of the mov­ing aver­age coef­fi­cients com­pared to the stan­dard para­me­ter­i­za­tion given above.
The arimax() func­tion from the TSA pack­age fits the trans­fer func­tion model (but not the ARIMAX model). This is a new pack­age and I have not yet used it, but it is nice to finally be able to fit trans­fer func­tion mod­els in R. Some­time I plan to write a func­tion to allow auto­mated order selec­tion for trans­fer func­tions as I have done with auto.arima() for regres­sion with ARMA errors (part of the fore­cast pack­age)