I have a system of two inputs and one output that I'd like to model using the following Box-Jenkins transfer function ("dynamic regression") structure:
$$y_t=\frac {\beta_1(B)}{\nu_1(B)}\omega_1Bx_{1,t}+\frac{\beta_2(B)}{\nu_2(B)}\omega_2Bx_{2,t}+\frac {\theta(B)}{\phi(B)}z_t+c$$
(using notation to conform with this Rob Hyndman post: http://robjhyndman.com/hyndsight/arimax/ where $B$ is the backshift operator, numerator polynomials like $\theta(B)$ stand for $(1+\theta_1 B+\theta_2 B^2+\cdots +\theta_nB^n)$, and denominator polynomials like $\phi(B)$ stand for $(1-\phi_1 B-\phi_2 B^2-\cdots -\phi_nB^n)$ )
My problem is this:
My $y$ is a stationary variable. However, both of my $x$s are not. (or to be precise, unlike my $y$, my $x$s are only very weakly stationary and have significant autocorrelation.) In discussions of pre-whitening that I've read, it's recommended to pre-whiten the $x$ variable (and then apply that same filter to $y$) before looking at the correlation structure. However, those examples only cover cases with only one $x$. On the other hand, this section of Hyndman's book: https://www.otexts.org/fpp/9/1 recommends repeatedly applying a difference to $y$ and all $x$s until all of the variables are stationary.
My question is:
Should I difference $y$ and all $x$s once before looking at the cross-correlations, even though this will make $y$ over-differenced? Or, should I pre-whiten $x_1$ and $x_2$, then apply each filter to the other variables, and then look at the cross-correlations? Or is another method recommended?