OLS - Non stationary variables but stationary residuals - is this OK or not? I am running an OLS on which the dependent variable (Y) and the independent variables (X1, X2, X3, ...) are non-stationary. But the residuals are found to be stationary. Does this mean my regression is OK? Or is it spurious? I tried searching online but I got a bit confused.
 A: 
...the dependent variable (Y) and the independent variables (X1, X2,
X3, ...) are non-stationary. But the residuals are found to be
stationary. Does this mean my regression is OK?...

This suggests there is a cointegrating relationship among the dependent and independent variable.
In this case, the OLS estimates are (super)-consistent.
For example, consider the simplest cointegration model:
$$
y_t = \beta x_t + \epsilon_t 
$$
where $\epsilon_t \stackrel{i.i.d.}{\sim} (0,1)$ and $u_t \stackrel{i.i.d.}{\sim} (0,1)$ are two independent white noise and
$$
x_t = \sum_{s = 0}^t u_s
$$ is a random walk with $\{ u_t \}$ as innovations.
Both $x_t$ and $y_t$ are not non-stationary but the population error term $y_t - \beta x_t$ is stationary. There is a cointegrating relationship with cointegrating vector $(1, -\beta)$.
The OLS estimate $\hat{\beta}$ satisfies
$$
T ( \hat{\beta} - \beta ) = \frac{\sum_1^T x_t \epsilon_t}{\sum_1^T x_t^2}
\stackrel{d}{\rightarrow} \frac{\int_0^1 W^{(1)}_t dW^{(2)}_t}{\int_0^1 \left( W^{(1)}_t \right)^2 dt},
$$
where $W^{(1)}$ and $W^{(2)}$ are independent standard Brownian motions.
Therefore $\hat{\beta} - \beta = O_p(\frac{1}{T})$, i.e. $\hat{\beta}$ is T-consistent. In contrast, in the stationary case, $\hat{\beta}$ is $\sqrt{T}$-consistent.
Here (T-)consistency holds even when $\epsilon_t$ and $u_t$ are correlated.
This is in contrast to the stationary case where consistency requires $E[\epsilon_t x_t] = 0$.
Unlike the stationary case, where inference is obtained via asymptotic normal distributions, here the asymptotic distribution of $T ( \hat{\beta} - \beta )$ is non-normal. Same goes for the asymptotic distribution of the usual Wald/F-statistics (critical values can be obtained by simulation, if needed).
Notice that sum of squared residuals of this cointegration regression is
$$
\sum_1^T \epsilon_t ^2 - \frac{( \sum_1^T x_t \epsilon_t)^2}{\sum_1^T x_t^2} = O_p(T) - O_p(1) = O_p(T),
$$
which implies the residuals are stationary.
On the other hand, if the regression is spurious, the sum of squared residuals is
$$
\sum_1^T y_t ^2 - \frac{( \sum_1^T x_t y_t)^2}{\sum_1^T x_t^2} = O_p(T^2) - O_p(T^2) = O_p(T^2)
$$
which implies the residuals are non-stationary.
The fact that you found residuals to be stationary
suggests your regression is cointegrated, rather than spurious.
Further Comments

*

*In applying unit root tests to residuals to check for non-stationarity, standard critical values cannot be used. For example, for the Augmented Dickey-Fuller test-statistic computed using residuals, the Engel-Granger critical values should be used.


*As already pointed out in comments, for a cointegration model there exists a corresponding error correction model (ECM). For the above simple example, the corresponding ECM is
$$
\Delta y_t = \Delta x_t + \underbrace{ (y_{t-1} - \beta x_{t-1})}_\text{$\epsilon_t$}.
$$
While the cointegration model describe the long-run equilibrium relationship between $x_t$ and $y_t$, the corresponding ECM describes the short-run relationship between $\Delta y_t$, $\Delta x_t$, and the deviation $\epsilon_t$ from long-run level. Notice that all variables in the ECM are stationary. Which model you estimate depends on the question of interest. To estimate the ECM, regress $\Delta y_t$ on $\Delta x_t$ and $e_t$ where $e_t$ is the residual from the cointegration regression.
