Simulation Ornstein-Uhlenbeck continuous I'd like to simulate a continuous Ornstein-Uhlenbeck-process.
Just recently I  tried to simulate gaussian-white-noise and had the theoretical result by hand. I was told that, in order to simulate, for example, $\int_0^t\epsilon(t')\mathrm{d}t'$, where $\epsilon(t)$ is gaussian-distributed $\mathcal{N}(0,\sigma^2)$ and fullfilling the gaussian-noise condition, i had to do the step that:
$\int_0^t\epsilon(t')\mathrm{d}t'\longrightarrow \Delta t\sum_i \gamma_i$,
where $\gamma_i$ is distributed like $\mathcal{N}(0,\sigma^2/\Delta t)$.
In this case, i found my mistake because i could calculate the process also theoretically and found the discrepancies.
In the case of the Ornstein-Uhlenbeck-process (or possibly others) I have no clue how to compare my simulated results to 'the real ones', especially because my function-depencendence on the stochastic variables becomes  more complex.
So my question is: 
How can i simulated the continuous Ornstein-Uhlenbeck-process in discrete time-steps.
Thanks already!
 A: For the Langevin equation
\begin{align*}
\dot{x} = h(x, t) + g(x,t)\Gamma(t),
\end{align*}
where $\Gamma(t)$ is such that
\begin{align*}
\langle \Gamma(t) \rangle = 0, \qquad\langle \Gamma(t) \Gamma(t') \rangle = 2\delta(t-t'),
\end{align*}
then the infinitesimal properties of the trajectory are given by
$$
\begin{align*}
D^{(1)}(x,t) &:= \lim_{\Delta t \rightarrow 0} \frac{1}{\Delta t}\langle X(t+\delta t)- X(t) \rangle_{X(t) = x(t) } = h(x,t) + \frac{\partial g(x,t)}{\partial x}g(x,t) \\
D^{(2)}(x,t) &:= \frac{1}{2}\lim_{\Delta t \rightarrow 0} \frac{1}{\Delta t} \langle \left[X(t+\Delta t) - X(t) \right]^2 \rangle_{X(t)=x(t)} = g^2(x,t)
\end{align*}
$$
Then we can approximate a discrete sample from the stochastic differential equation by
$$
\begin{align*}
x_{n+1} \approx x_n + D^{(1)}(x_n,t_n)\Delta t + \left( \sqrt{D^{(2)}(x_n,t_n) \Delta t} \right)w_n, \qquad w_n \sim \mathcal{N}(0, 2)
\end{align*}
$$
so that as $\Delta t \rightarrow 0$ the approximated process has the same infinitesimal properties as the theoretical process. Now for an OU process let us say
$$
x(t) = -\gamma x + \sigma \Gamma(t), \qquad \gamma > 0
$$
so $D^{(1)} = -\gamma x$ and $D^{(2)}=\sigma^2$ which you can now simulate as described above. However for this process you can also solve the transition probability exactly to give 
$$
p(x_{n+1}, t_{n+1} | x_n , t_n ) = \sqrt{\frac{\gamma}{2 \pi \sigma^2 (1 - e^{-2 \gamma(t_{n+1}-t_n)} )}} \exp\left( -\frac{\gamma(x_{n+1}-e^{-\gamma(t_{n+1}-t_n)}x_n)^2}{2\sigma^2(1-e^{-2\gamma(t_{n+1}-t_n)}) }\right)
$$
whichis a Gaussian random variable with mean $e^{-\gamma (t_{n+1}-t_n)}x_{n}$ and variance $\sigma^2 (1-e^{-2\gamma(t_{n+1}-t_{n})})/\gamma$
So that is two ways to simulate your trajectories, the first an approximation and the second exact.
