Dealing with different definitions of the Ornstein-Uhlenbeck process I've run up against a wall in reconciling two different definitions of the Ornstein-Uhlenbeck process, and would appreciate some help.
On the one hand, as discussed here, we can define an Ornstein-Uhlenbeck process as a Gaussian process with a kernel function of the form $k(x, y) \propto \exp(-\left\lVert x - y\right\rVert / \theta)$, which is clearly stationary.
On the other hand, we have the definition of the Ornstein-Uhlenbeck process as the solution to the stochastic differential equation $du(t) = \theta(\mu - u(t))+\sigma dW(t)$, which is given by $u(t)= u(0) \exp⁡(-\theta t)+\mu(1-\exp⁡(-\theta t) )+\sigma \exp⁡(-\theta t) \int_0^t \exp⁡(\theta\tau)dW(\tau)$. Using the Itô isometry, as shown here and here on SE, this process can be shown to have the covariance $k(x, y) = \frac{\sigma^2}{2\theta} \exp(-\theta (x+y)) (1 - \exp(2 \theta (x - y)))$.
Not only is this not equal to the previous definition, this is not even stationary. What gives? Are these two distinct definitions of what "Ornstein-Uhlenbeck process" means, or can these two processes be shown to be equal, possibly modulo certain assumptions?
As a bonus question, just because I've never seen anyone prove it - how do you prove that the solution to the above SDE is a Gaussian process (assuming it is)?
 A: In order to have a stationary solution of the Stochastic Differential
Equation (SDE), you have to start from a random initial value $u(0)$
at the fixed time $t=0$. This value must be drawn from the stationary
distribution, here Gaussian.
An alternative is to consider that the process $u(t)$ "has started in the remote 
past". Although $u(t)$ is not
differentiable, it is appealing to write the SDE as 
$$
   \frac{\text{d} u(t)}{\text{dt}} = - \theta u(t) + \eta(t),
$$
where $\eta(t)$ is a Gaussian white noise with variance $\sigma^2$,
leading to the
following representation of $u(t)$
$$
    u(t) = \int_{-\infty}^t e^{-\theta (t -s)} \, \eta(s) \,\text{d}s
$$
which can be made more rigorous by replacing $\eta(s)\text{d}s$ by
$\sigma \text{d}W(s)$. Using this representation, it is easy to
derive the covariance of $u(t)$. For $r \geq 0$
$$
   \text{E}[u(t) u(t+r)]  = \sigma^2 \, \text{E}
   \left[\int_{-\infty}^t\int_{-\infty}^{t +r}
    e^{-\theta (t -s)} \, e^{-\theta (t + r -s')} \,
    \text{d}W(s) \text{d}W(s') \right],
$$
and we can exchange the expectation and the integral(s), the expectation
$\text{E}[\text{d}W(s) \text{d}W(s')]$ being $0$ for $s \neq s'$
and $\text{d}s$ when $s= s'$. Then by simple integration, we find 
$$
   \text{E}[u(t) u(t+r)]  = \frac{\sigma^2}{ \, 2 \theta} \, e^{-\theta r}
   \quad (r \geq 0).
$$
Most of the literature dedicated to Itô's calculus assumes that $W(0)
= 0$; This can be quite misleading for statistical problems where
there is usually no special time "$t=0$" corresponding to a
deterministic value of the process. Some adaptations are needed to cope with continuous time linear filters.
A: I initially believed that the mean parameter had a part to play as, in this link Rasmussen & Williams, 2006, Appendix B, equation B.27, the OU SDE is written as:
$$ dX_t = - a_0 X(t)dt + b_0 dW_t $$
whereas the traditional equation is:
$$ dX_t = \theta \mu dt - \theta X_tdt + \sigma dW_t $$
But this is not the case, the covariance is not dependant on the parameter $\mu$ as Billy pointed out.
We'd expect the moments not to depend on location in a stationary process. The OU isn't stationary even when the parameter $\mu$ is 0, evident in Wikipedia's solution, given below. This in fact, tends to the RW solution:
$$Cov(X_t, X_s) = \frac{\sigma^2}{2\theta} \left( e^{-\theta|t-s|} - e^{-\theta(t+s)} \right) \rightarrow \frac{\sigma^2}{2\theta} e^{-\theta|t-s|} $$
... as t and/or s increase. So, as the process gets farther and farther from its initial condition, it tends to a stationary process.
I can only speculate that they did not include the full non-stationary solution and provided the limiting case or perhaps the method of obtaining the covariance from the power spectrum assumes stationarity perhaps, I do not know.

As for the proof of the the OU process being a GP, it's possible to solve the Fokker-Planck equation - you end up with the normal density (see the Wikipedia link). Hence, any at any finite set of points $\bf t$, $X$ will have a multivariate normal distribution, which is the definition of a GP.
