# Fit an ARMAX model in R

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?

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.
ARIMA(y~lag(y)+x+pdq(p=0,d=0,q=1)