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?
 A: 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.
