Estimate a regression model with ARMA errors?

Mathematica provides a good estimation method for pure time series process and regression model. But how to estimate a regression model with ARIMA errors? The model is

y = x*b + u,
where u ~ ARIMA process.

I think it is big shortage for Mathematica to estimate some regression model. I believe it can be improved a lot if it allows the users to specify some exogenous variables in its existing distribution functions.

Thanks.