Monday, November 25, 2013

The ARIMAX model muddle

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 y_1,\dots,y_n. First, we will define an ARMA(p,q) model with no covariates:
    \[ y_t = \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} - \theta_1 z_{t-1} - \dots - \theta_q z_{t-q} + z_t, \]
where z_t 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:
    \[ y_t = \beta x_t + \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} - \theta_1 z_{t-1} - \dots - \theta_q z_{t-q} + z_t \]
where x_t is a covari­ate at time t and \beta 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 \beta is not the effect on y_t when the x_t 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 \beta 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
    \[ \phi(B)y_t = \beta x_t + \theta(B)z_t \qquad\text{or}\qquad y_t = \frac{\beta}{\phi(B)}x_t + \frac{\theta(B)}{\phi(B)}z_t, \]
where \phi(B)=1-\phi_1B -\cdots - \phi_pB^p and \theta(B)=1-\theta_1B-\cdots-\theta_qB^q. 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.
    \begin{align*} y_t &= \beta x_t + n_t\\ n_t &= \phi_1 n_{t-1} + \cdots + \phi_p n_{t-p} - \theta_1 z_{t-1} - \dots - \theta_q z_{t-q} + z_t \end{align*}
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
    \[ y_t = \beta x_t + \frac{\theta(B)}{\phi(B)}z_t. \]

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:
    \[ y_t = \frac{\beta(B)}{v(B)} x_t + \frac{\theta(B)}{\phi(B)}z_t. \]
This allows for lagged effects of covari­ates (via the \beta(B) oper­a­tor) and for decay­ing effects of covari­ates (via the v(B) 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 \phi(B) with \nabla^d\phi(B) where \nabla=(1-B) denotes the dif­fer­enc­ing oper­a­tor. Notice that this is equiv­a­lent to dif­fer­enc­ing both y_t and x_t 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)
Posted by at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)