Can I difference after fitting a time series regression model?

Suppose that I have a time series that exhibits a notable trend, and I want to test a hypothesis that a second variable is related to that trend. I fit a linear regression model with that second variable as a predictor, and it accounts for a significant amount of variance. However, there is some part of the trend that it does not account for.

I then want to include AR or MA terms to the model to account for the autocorrelated residuals, but ARIMA models assume a stationary series, and visual inspection of the series and an ADF test indicate that the residual error series is not stationary in the mean (i.e., there is some trend left unaccounted for by the regression model). Is it fine to difference after the regression model has been fit to remove this source of non-stationarity, or should one always try and find a regression model that explains all the trend in the data?

• You can fit an ARIMA model along with external regressors. If you see that the residuals from the linear regression are not stationary, then this may confirm that you should use an ARIMA model (not an ARMA model) with external regressors. You could split the process in two steps and fit an ARMA model to the differenced residuals, but a joint estimate of the parameters may yield better results. – javlacalle Jul 17 '14 at 20:21
• Thanks, but how could an ARMA (after differencing the residuals) lead to different results than an ARIMA model--doesn't the ARIMA model essentially split the process in two steps as well, as the AR/MA terms require an already stationary series? – ATJ Jul 17 '14 at 21:37

It may depend on the implementation and method. For example the stats::arima function in the R software considers three methods: 1) minimisation of the conditional sum of squares (CSS), 2) maximum likelihood (ML) and 3) a combination of the other methods (CSS-ML).
The first method, CSS, splits the process in two steps: First, the coefficients of the regressors are estimated in a regression of the differenced series on the differenced regressors. (The differencing filter of the selected ARIMA model is applied to the observed series and to the regressors.) Then, the CSS is minimized for the series adjusted for the effects of the external regressors, $y_t^a = y_t - X_t \delta$.
In the second method, ML, the vector of parameters that is passed to the optimization algorithm contains the AR and MA coefficients as well as the coefficients of the external regressors. All of them take part in the optimization algorithm and in the objective function. The function that is minimised removes the effects of the external regressors $X_t$ from the observed series, $y_t$, given the vector of coefficients at the current iteration of the optimization algorithm, $\delta^{(i)}$, i.e.: $y_t^a = y_t - X_t \delta^{(i)}$. Then the objective function evaluates the negative of the likelihood function for the adjusted series $y_t^a$ and for the corresponding ARIMA model. Thus, all the parameters are estimated in a single step (ignoring a preliminary step to initialise the coefficients of the regressors).
In general I would say it is better to follow the one-step approach or to somehow refine the values of the two-steps process, for example iterating the process or using those values as starting values for the one-step procedure (this last option is the default in stats::arima).