Fit an ARMAX model in R

Question

I would like to fit an ARMAX model in R of the form that is mostly used in literature:

$$y_t = \beta_1 x_t+\cdots+ \beta_{k} x_{t_{k-1}}+ \phi_1 y_{t-1}+\cdots+\phi_p y_{t-p}+ \theta_1 z_{t-1}+\cdots \theta_{q} z_{t-q} + z_t$$ where $x$ is exogenous input and $z$ is white noise. Using the backshift operator $B$ with $B^k(y_t)=y_{t-k}$ the equation becomes

\begin{align}\tag{1}\label{a} \phi(B)y_t = \beta(B)x_t + \theta(B)z_t \end{align}

I know that R's built-in arima function includes exogenous input as follows:

$$ y_t = \beta x_t + n_t$$ $$ \varphi(B) n_t = \theta(B) z_t $$

But this is not want I want. I did some research and found out that there are (at least) three possible functions that fit ARMA models with exogenous variables:

1) stats:::arima (built-in)

2) forecast:::Arima

3) TSA:::arima/arimax (2 identical functions with different names)

As @Rob Hyndman writes in 2010 (https://robjhyndman.com/hyndsight/arimax/), functions 1) and 2) build a model of the form I don't want, whereas function 3) builds a transfer function model of the form

$$ y_t = \frac{\beta(B)}{\nu(B)}x_t + \frac{\theta(B)}{\phi(B)}z_t$$

An ARMAX model is a special case of the above transfer function model with $\nu=\phi$. However, (as far as I know) the TSA:::arimax function does not provide a possibility to respect this constraint. Therefore, it is not possible to fit a "true" ARMAX model.

My questions are:

1) Is what I write correct?

2) Since most of what I have read about this is 6-8 years old, is there a way to fit an ARMAX model in the form of (\ref{a}) in R?

Answer 1

As far as my research on this topic has taken me, I agree that that the arima/Arima functions from the stats and forecast packages do not fit transfer functions as you mention, but instead a linear model with ARMA errors.

I don't see the possibility to tell the TSA::arimax function that the $\nu(B)$ should be equal to $\phi(B)$. But you can force the order of all the individual polynomial in $B$ to be what you want, including $0$. But that does not really give you your ARMAX model.

The last thing I can suggest is to take a look at the MARIMA package. It should somehow be able to fit an ARMAX model, but I am not 100 % sure about the procedure.

Answer 2

As of 2022, Hyndman is using fable::ARIMA, but in his otherwise excellent time series regression guide (https://otexts.com/fpp3/) is still only showing how to do a linear regression with arima errors. However, it appears that as the fable package lets you specify the model via a formula, then you can do the following, which I think is a ARMAX(1,0,1) model (seems to work):

ARIMA(y~lag(y)+x+pdq(p=0,d=0,q=1)