Suppose that we have a regression model with ARIMA errors. We need to determine the appropriate level of differencing, so that the errors from the regression are actually stationary. In Hyndman's
auto.arima function (https://github.com/robjhyndman/forecast/blob/master/R/newarima2.R), the appropriate level of differencing is determined using the kpss test on the residuals from the usual regression
# Finally find residuals from regression in order # to estimate appropriate level of differencing j <- !is.na(x) & !is.na(rowSums(xregg)) xx[j] <- residuals(lm(x ~ xregg)) ... d <- do.call("ndiffs", c(list(dx, test = test, max.d = max.d), test.args))
However, the coefficients of the regression, and, of course also the errors, are actually not consistent because they do not take the coefficients of ARMA into consideration as mentioned in Hyndman's book (https://otexts.com/fpp2/estimation.html):
When we estimate the parameters from the model, we need to minimise the sum of squared $\epsilon_t$ values. If we minimise the sum of squared $\eta_t$ values instead (which is what would happen if we estimated the regression model ignoring the autocorrelations in the errors), then several problems arise.
If the errors are not consistent with the underlying model, then how can we apply any test to the simple regression errors? These errors (as above) are actually non-existent in the final estimation, which is done with combined regression + Kalman Filter at the same time.
Is it legal to perform such a test as above? What is the theoretical justification, given that it is not actually allowed to first regress on covariates to find coefficients, and then use ML estimation for ARMA errors.
PS. Any suggestion is good, but a proof/equation would be nice too.