There is often confusion about how to include covariates in ARIMA models, and the presentation of the subject in various textbooks and in R help files has not helped the confusion. So I thought I’d give my take on the issue. To keep it simple, I will only describe non-seasonal ARIMA models although the ideas are easily extended to include seasonal terms. I will include only one covariate in the models although it is easy to extend the results to multiple covariates. And, to start with, I will assume the data are stationary, so we only consider ARMA models.
Let the time series be denoted by 
. First, we will define an ARMA
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, \]](http://robjhyndman.com/researchtips/wp-content/ql-cache/quicklatex.com-fc08bbfd6fd4309a25481dadea9317b7_l3.png "Rendered by QuickLaTeX.com")
where 
is a white noise process (i.e., zero mean and iid).
ARMAX models
An ARMAX model simply adds in the covariate 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 \]](http://robjhyndman.com/researchtips/wp-content/ql-cache/quicklatex.com-906204c6ba691ca0db596b607609a608_l3.png "Rendered by QuickLaTeX.com")
where 
is a covariate at time 
and 
is its coefficient. While this looks straight-forward, one disadvantage is that the covariate coefficient is hard to interpret. The value of is not the effect on 
when the is increased by one (as it is in regression). The presence of lagged values of the response variable on the right hand side of the equation mean that can only be interpreted conditional on the value of previous values of the response variable, which is hardly intuitive.
If we write the model using backshift operators, 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, \]](http://robjhyndman.com/researchtips/wp-content/ql-cache/quicklatex.com-fcb536f4ac62c3451f92f0f881bf914c_l3.png "Rendered by QuickLaTeX.com")
where 
and 
. Notice how the AR coefficients get mixed up with both the covariates and the error term.
Regression with ARMA errors
For this reason, I prefer to use regression models with ARMA errors, defined as follows.
In this case, the regression coefficient has its usual interpretation. There is not much to choose between the models in terms of forecasting ability, but the additional ease of interpretation in the second one makes it attractive.
Using backshift operators, this model can be written as
![\[ y_t = \beta x_t + \frac{\theta(B)}{\phi(B)}z_t. \]](http://robjhyndman.com/researchtips/wp-content/ql-cache/quicklatex.com-e6c3d58af72edcaf46f14fca9b2dfe1e_l3.png "Rendered by QuickLaTeX.com")
Transfer function models
Both of these models can be considered as special cases of transfer function models, popularized by Box and Jenkins:
![\[ y_t = \frac{\beta(B)}{v(B)} x_t + \frac{\theta(B)}{\phi(B)}z_t. \]](http://robjhyndman.com/researchtips/wp-content/ql-cache/quicklatex.com-f22708e52c2b57af4470f37f892ff058_l3.png "Rendered by QuickLaTeX.com")
This allows for lagged effects of covariates (via the 
operator) and for decaying effects of covariates (via the 
operator).
Sometimes these are called “dynamic regression models”, although different books use that term for different models.
The method for selecting the orders of a transfer function model that is described in Box and Jenkins is cumbersome and difficult, but continues to be described in textbooks. A much better procedure is given in Pankratz (1991), and repeated in my 1998 forecasting textbook.
Non-stationary data
For ARIMA errors, we simply replace 
with 
where 
denotes the differencing operator. Notice that this is equivalent to differencing both and before fitting the model with ARMA errors. In fact, it is necessary to difference all variables first as estimation of a model with non-stationary errors is not consistent and can lead to “spurious regression”.
R functions
The arima() function in R (and Arima() and auto.arima() from the forecast package) fits a regression with ARIMA errors. Unfortunately there is a bug in the code so that if a non-stationary error is specified, the data are not differenced before estimation of the regression coefficients. Consequently, you should do your own differencing when required. Note also that R reverses the signs of the moving average coefficients compared to the standard parameterization given above.
The arimax() function from the TSA package fits the transfer function model (but not the ARIMAX model). This is a new package and I have not yet used it, but it is nice to finally be able to fit transfer function models in R. Sometime I plan to write a function to allow automated order selection for transfer functions as I have done with auto.arima() for regression with ARMA errors (part of the forecast package)
Related Posts:
- Statistical tests for variable selection
- Exponential smoothing and regressors
- Cyclic and seasonal time series
- Why every statistician should know about cross-validation
- Initializing the Holt-Winters method
Rob J Hyndman