Simulating an Ornstein-Uhlenbeck Process

From Qwiki

$dX_t = - \gamma X_t dt + \sqrt{2 D} dW_t$

Simulating a one-dimensional Ornstein-Uhlenbeck process in MATLAB is quick and easy. At first glance, it appears that one needs to calculate the updates in a FOR loop, since there's an $X t$ on the right-hand side of the differential update rule. We can get around this problem with MATLAB's wonderful "filter" function. Given a vector 't' of times, linearly spaced at intervals 'dt', with 'g' and 'D' your damping parameter and diffusion coefficient respectively, then the following lines of code will simulate the Ornstein-Uhlenbeck process:

X_free = cumsum(sqrt(2*D*dt)*randn(size(t)));
X_filt = filter([0 g*dt], [1 -1+g*dt] , X_free);
X_OU = X_free-X_filt;

'X_OU' is a realization of the Ornstein-Uhlenbeck process sampled at the times in 't'. The simulation is exact, because a realization of an Ornstein-Uhlenbeck process is equal to a realization of free diffusion minus a low-pass filtered copy of the same realization. It is easy to see why this trick works. Let $X f r e e$ be free diffusion and $X f i l t$ be a low-pass filtration of $X f r e e$ :. Then

$dX_{free} = \sqrt{2D}dW$
$dX_{filt} = -\gamma\left(X_{filt}-X_{free}\right) dt$