time series - ARMAX in Mathematica

Mathematica has extensive support for time series, but I haven't found anything with exogenous signal. Is it possible to fit a time series with exogenous signal in Mathematica? For example identifying a discrete linear system?

For example, let's suppose there is a heat exchange with inputs $u_1$ (flow) and $u_2$ (input temperature) and the output $y$ as the output temperature. The problem in concern in to fit an ARMA model with exogenous signal $u_1$ and $u_2$.

Given the sampled signals $u_1(k)$, $u_2(k)$ and $y(k)$, $k=1...n$. Find the vectors $a=[a_1,...a_{n_a}]$, $b=[b_1,...b_{n_b}]$ and $c=[c_1,...c_{n_c}]$ such that the process is modeled by:

$$y(k)=\sum_{i=1}^{n_a} a_iy(k-i)+\sum_{i=1}^{n_b} b_iu_1(k-i)+\sum_{i=1}^{n_c} c_iu_2(k-i)+\varepsilon(k-1)$$ where $\varepsilon$ is a noise.

This is a simple and common problem and it is strange that in Mathematica there is ARMA process but not ARMAX process.