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?