Cross correlations I have a question about cross correlations and multivariate time series.  I've read several articles and posts on SE about how to properly prewhiten data for cross correlation of two time series, and I'm seeking some clarification.   This is a summary of what I have read:
Approach 1: make data stationary through first differencing for x and y, then compute cross correlations.
Approach 2: run arima on x and arima on y, save residuals for each, run cross correlations on residuals.  
Approach 3: run arima on y, save residuals.  Run cross correlation on stationary x and residuals of y.
Approach 4: run arima on x, save residuals.  Running cross correlation between residuals of x and stationary y.
For obvious reasons, each approach yields very different results.  
I'm looking for the best approach for selecting potential exogenous variables for a forecast model.  
Any thoughts are appreciated.
 A: The Problem
To understand which approach to take let's look at a toy example.  Let's assume the true unobserved DGP for $x_t$ and $y_t$ are 
$$
y_t = \alpha y_{t-1} + \mu x_t + \phi_t
$$
$$
x_t = \beta x_{t-1} + \gamma y_t + \epsilon_t
$$
If we calculated a cross-validation of $x_t$ and $y_t$ we would get misleading results in the case.  For example, the first element of the cross variance is $$\frac{\sum_{t=1}^T(x_t-\bar{x})\sum_{t=1}^T(y_{t}-\bar{y})}{\sigma_x\sigma_y}$$ 
Substituting in the DGPs, the numerator is 
$$\sum_{t=1}^T(\beta (x_{t-1} - \bar{x}) + \gamma (y_t - \bar{y}) + \epsilon_t)\sum_{t=1}^T(\alpha(y_{t-1} - \bar{y}) + \mu(x_t-\bar{x}) + \phi_t)$$ 
This is looking at the variance between $x_t$ and $y_t$ (i.e. no lags).  Now let's assume that $\gamma =0$ and $\mu =0$,
$$\sum_{t=1}^T\beta (x_{t-1} - \bar{x})\sum_{t=1}^T\alpha(y_{t-1} - \bar{y}) \neq 0 $$.
This is clearly problematic, there is no relationship between $x_t$ and $y_t$, however, the variance is not zero (though it may collapse to zero asymptotically)!  The auto-correlation incorrectly implies a relationship, when there is not one.
The fix
To fix, we must identify a filter that will get rid of this serial correlation.  Let's assume that if we run an AR(1) model on $x_t$ the resulting errors are not serial correlated (i.e. white noise).  Note, whatever model you use to filter (AR, ARMA, ARIMA) the resulting errors must be white noise for this to work!
The AR(1) is $x_t = bx_{t-1} + u_t$ and the filter is $\text{f}(x) = x_t - \hat{b}x_{t-1}$
You then apply this filter to both $x_t$ and $y_t$.
How does it work?
Let's again look at our example above when $\gamma =0$ and $\mu =0$.  In this case,  $\hat{b}=\beta$, as we are estimating the true DGP,
$$\sum_{t=1}^T\ \bigg(f(x_t) - \frac{\sum_{t=1}^Tf(x_t)}{T}\bigg)\sum_{t=1}^T\ \bigg(f(y_t) - \frac{\sum_{t=1}^Tf(y_t)}{T}\bigg) = 0*\sum_{t=1}^T\ \bigg(f(x_t) - \frac{\sum_{t=1}^Tf(x_t)}{T}\bigg) =  0 $$
Good, that make sense!
Comparing to your approaches
Approach 1:
A first difference filter is close to, but not the same as $f(x)$. It is missing the $b$ coefficient (unless it follows a unit root).  Moreover, this approach does not guarantee that the serial correlation problem is solved.
Approach 2:
May work, but I would not want to filter out the auto correlation in y.  We're trying to explain x with y and don't want to filter our important information.
Approach 3:
This one does not make sense.  If you want to explain x with y, you'll still have serial correlation in x, which is bad.
Approach 4:
This is you winner!  Make sure you that the errors in you ARIMA are white noise.
Unit roots
Several times in your post your refer to making your processes stationary.  If your data follows a stochastic trend (i.e. unit root), I would difference until the data is stationary.  After, you can do approach 4.
