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?
 A: 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, does 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.
A: 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)

