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?
Thanks in advance!
 A: Answering the question in your comment:
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).
Under some circumstances splitting the vector of parameters and design 
a two-steps procedure may yield very similar results to a one-step procedure where the vector of parameters is not split.
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). 
