# How is Arima from forecast/stats package with external regressors (dynamic regression) evaluated?

I use the Arima function from the forecast package in R. I also took a look at this short introduction to the topic (author of the forecast package): https://otexts.org/fpp2/dynamic.html.

However, I am not sure, how are coefficients for external regressors (xreg specified) + coefficients for AR and MA terms calculated. I have 2 theories:

1) First, a normal regression with external regressors is fit. Next, we fit ARMA errors to the OLS errors we got from the regression model, and find AR and MA coefficients from there. This one contradicts what is written in the source I mentioned:

"When we estimate the parameters from the model, we need to minimise the sum of squared $$\epsilon$$ values. If we minimise the sum of squared $$\eta$$ values instead (which is what would happen if we estimated the regression model ignoring the autocorrelations in the errors), then several problems arise"

2) Regressor coefficients and AR + MA coefficients are fit at the same time together. However, OLS cannot be used, because clearly the error is not i.i.d. Gaussian. Thus, is GLS used?

I noticed that forecast package references stats::arima internally, where this estimation happens. However, I cannot figure from code how all coefficients (regressor + AR + MA) are estimated. Can anyone give a hint? I at least would like to know which method is used there: 1 or 2, and if 2, what is the name, is it GLS?

## 1 Answer

Neither 1 nor 2, though 2 is closer to the truth. Regressor coefficients and AR + MA coefficients are fit simultaneously by maximum likelihood estimation. Indeed, OLS cannot be used there for the reason you give (though Gaussianity is not a necessary prerequisite for OLS estimation to yield consistent and minimum-variance unbiased estimators; it is only needed to make OLS coincide with maximum likelihood estimation). GLS does not work either because some of the regressors, namely, lagged errors are unobservable. Therefore, there is no way to form the design matrix $$X$$ that is an essential element in both OLS and GLS estimation. Thus we are left with maximum likelihood.

How is the maximum likelihood implemented? Here is an answer from the help file of the arima function:

the exact likelihood is computed via a state-space representation of the ARIMA process, and the innovations and their variance found by a Kalman filter,

And here is the Gartner et al. (1980) paper describing the algorithm used in R.

• Wow, that was fast. Thank you. Btw, could you maybe give me some reference to ML estimation in this type of model (tutorial will do too)? – SWIM S. Oct 30 '18 at 18:30
• That is the hard part. Estimation of ARMA-type models is nasty. There are many ways to choose from, none of which is easy (in my opinion). For example, you can write down the likelihood function and submit it to a generic optimizer; that will not be computationally efficient, as far as I know. You can use Kalman filter which is (in my opinion) a complicated machine but it does the job. EM algorithm is another possibility, probably slower than Kalman filter but better than direct optimisation of the likelihood. Hamilton's "Time Series Analysis" (1994) has a good overview of some common methods. – Richard Hardy Oct 30 '18 at 18:36
• I see. Do you know if the stats package use the Kalman filter, or likelihood formula? – SWIM S. Oct 30 '18 at 18:40
• It seems to be using Kalman filter: The exact likelihood is computed via a state-space representation of the ARIMA process, and the innovations and their variance found by a Kalman filter, see the help file for the arima function. And here is the Gartner et al. (1980) paper describing the algorithm used in R. – Richard Hardy Oct 30 '18 at 18:43