Interpretation of TSA::arimax output model is presented in R

I have a relatively simple problem, but yet taking some time to solve it. I am using the arimax() function from the TSA package. (Note: not arima() from the stats package.) This is the model:

out <- arimax(sub_s_t_series, order=c(2,0,1), xreg=sub_r_t_series, method=c("ML"))

and these are my coefficients:

Call
arimax(x = sub_s_t_series, order = c(2, 0, 1), xreg = sub_r_t_series, method = c("ML"))

Coefficients:
         ar1      ar2      ma1  intercept     xreg
      1.4825  -0.6613  -0.8516  52745.107  -1.0132
s.e.  0.0295   0.0294   0.0064     40.828   0.0012

sigma^2 estimated as 0.08929:  log likelihood = -105.98,  aic = 221.97

All I am trying to do is to interpret the results. According to my understanding and the help given in the TSA package, the above ARIMAX(2,0,1) model is represented as follows:
$$
{rm sub_s_t_series_hat[k]} = {rm intercept} + xregtimes {rm sub_r_t_series[k]} +
frac{a_{t[k]}+ma1*a_{t[k-1]}}{a_{t[k]}-ar1*a_{t[k-1]}-ar2*a_{t[k-2]}} tag{1}
$$
where $a_t$ are the residuals. When I use e_t = fitted(out)-sub_s_t_series_hat to measure the error / residuals myself, e_t matches exactly to the values obtained by out[["residuals"]].

But when I use (1) as follows: e_t_hat = sub_s_t_series_hat – sub_s_t_series,
e_t_hat does not match with out[["residuals"]], in fact the results deviate by a magnitude of almost 4.

My questions is: did an ARIMAX(2,0,1) fit would result in (1) or am I missing something?