Stationarity requirements of using regression with ARIMA errors - ROUND 2 Referring to this Post: What are the stationarity requirements of using regression with ARIMA errors for inference? I would like to seek a confirmation of the below practice:

The situation is as below:
1) I am trying to form a Multiple Regression with ARIMA errors
2) Y = B0 + B1X1 + B2X2 + e, where e is modelled with ARIMA.
3) Y, X1, and X2 are all NON-STATIONARY series.
4) Y and X1 are cointergated; but Y and X2 are NOT cointergated.

The Questions:
a)-  Do I merely need to do differencing (or do what so ever) to make X2 become stationary, but no need to do so for Y and X1 ?
So the model will simply become  Y = B0 + B1X1 + B2D(X2,1) + e , where e is modelled with ARIMA?

b)- If this is the case, does it mean that: as long as the variables are cointegrated, ARIMA errors can fix the non-stationary problem? So the ARIMA error is somehow another way of "Error Correction"?
But from my gut feeling, regression between non-stationary series may generate a spurious regression problem, it makes R2 extremely high, which implies less error is composed in such model. With less error, even though an ARIMA-error term is added, it may not solve the problem as error is merely contributing a small part...
So, I am not convinced with this method and it seems making all Y and Xs become stationary is the only way to do regression, no matter with or without ARIMA error terms. What do you guys think?
 A: I think you should estimate an error correction model (ECM) of the form 
$$ \Delta y_t = \alpha (y_{t-1} - \beta x_{t-1}) + \sum_{i=1}^{p}{\Delta y_{t-i}} + \sum_{j=1}^{q}{\Delta x_{1,t-j}} + \sum_{i=k}^{r}{\Delta x_{2,t-k}} + \varepsilon_t $$
where the possibly long sums of lagged $\Delta y$, $\Delta x_1$ and $\Delta x_2$ are used to approximate the ARMA (not ARIMA) structure in the original error term (here $\varepsilon_t$ has no ARMA structure). Also, a constant and/or a time trend may or may not be included in the error correction term. Such a model can be estimated by OLS. 
Alternatively, you could perhaps have shorter sums of lagged variables but then include some lagged errors (MA elements); the model would have to be estimated by ML. I am not sure if there are standard routines that can do that (you would perhaps have to write your own function).
The benefit of an ECM over the model for levels is that the model for levels has some non-linear parameter constraints (easily derived from the ECM representation by adding $y_{t-1}$ to both sides of the equation above) that may be difficult to impose given the off-the-shelf software. The benefit of ECM over the model in first differences is that is has the error correction term which is omitted in the model in first differences.

I would not say that ARMA (not ARIMA) errors fix the problem of non-stationarity. It is cointegration and ECM that are key here. Even with simple i.i.d. errors the ECM is the solution to deal with non-stationarity when the variables are cointegrated; ARMA errors are a side issue.

Regarding $R^2$, if you already know that $y$ and $x_1$ are cointegrated, you need not worry about spurious regression problem. Although the $R^2$ may be high, just keep in mind that this is due to cointegration. From the subject-matter perspective, this $R^2$ then should be of little interest; you could be more interested in the $R^2$ associated with the change in $y$ versus the change in $x_1$ and the error correction term.
