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!
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.