Implementing Ornstein–Uhlenbeck in Matlab

Question

I am reading this article on Wikipedia, where three sample paths of different OU-processes are plotted. I would like to do the same to learn how this works, but I face troubles implementing it in Matlab.

I think I have to discretize this equation somehow: $ x_t = x_0 e^{-\theta t} + \mu(1-e^{-\theta t}) + \int_0^t \sigma e^{\theta (s)}\, \mathrm{d}W_s. \, $, but especially the integral equation confuses me a lot.

I also think I will have to use $W_t = W_t-W_0 \sim N(0,t)$ somehow, but don't know how yet...

Can someone please help me out? I am new to stochastic calculus, so please help me understand step by step.

Answer 1

The Wikipedia article you cite provides everything you need to evaluate the analytical solution of the Ornstein–Uhlenbeck process. However, for a beginner, I agree that it may not be very clear.

1. Simulating the Ornstein–Uhlenbeck process

You should first be familiar with how to simulate this process using the Euler–Maruyama method. The stochastic differential equation (SDE)

$$\mathrm{d}x_t = \theta (\mu - x_t)\mathrm{d}t + \sigma \mathrm{d}W_t$$

can be discretized and approximated via

$$ X_{n+1} = X_n + \theta (\mu - X_n)\Delta t + \sigma \Delta W_n$$

where $\Delta W_n$ are independent identically distributed Wiener increments, i.e., normal variates with zero mean and variance $\Delta t$. Thus, $W_{t_{n+1}}-W_{t_n} = \Delta W_n \sim N(0,\Delta t) = \sqrt{\Delta t} \space N(0,1)$. This can be simulated in Matlab very easily using randn to generate standard normal variates:

th = 1;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:2;             % Time vector
x = zeros(1,length(t)); % Allocate output vector, set initial condition
rng(1);                 % Set random seed
for i = 1:length(t)-1
    x(i+1) = x(i)+th*(mu-x(i))*dt+sig*sqrt(dt)*randn;
end
figure;
plot(t,x);

which will result in a plot something like the blue trace in the Wikipedia article (it won't be identical because different random values were used – see also my comments below this answer). Note that the above is not the most efficient code.

2. Solution in terms of integral

The equation in your question is in terms of a stochastic integral

$$x_t = x_0 e^{-\theta t} + \mu (1-e^{-\theta t}) + \sigma e^{-\theta t}\int_0^t e^{\theta s} \mathrm{d}W_s$$

To obtain a numerical solution in Matlab with this, you'll need need to numerically approximate (discretize) the integral term using an SDE integration scheme like Euler–Maruyama described above:

th = 1;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:2;             % Time vector
x0 = 0;                 % Set initial condition
rng(1);                 % Set random seed
W = zeros(1,length(t)); % Allocate integrated W vector
for i = 1:length(t)-1
    W(i+1) = W(i)+sqrt(dt)*exp(th*t(i))*randn; 
end
ex = exp(-th*t);
x = x0*ex+mu*(1-ex)+sig*ex.*W;
figure;
plot(t,x);

or without a for loop using cumsum:

th = 1;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:2;             % Time vector
x0 = 0;                 % Set initial condition
rng(1);                 % Set random seed
ex = exp(-th*t);
x = x0*ex+mu*(1-ex)+sig*ex.*cumsum(exp(th*t).*[0 sqrt(dt)*randn(1,length(t)-1)]);
figure;
plot(t,x);

3. Analytical solution

To compute the full analytical solution of an Ornstein–Uhlenbeck process for a given time series and corresponding Wiener increments, you'll need use a "scaled time-transformed" Wiener process:

$$x_t = x_0 e^{-\theta t} +\mu (1-e^{-\theta t}) + \frac{\sigma e^{-\theta t}}{\sqrt{2 \theta}}W_{e^{2 \theta t}-1}$$

See Doob 1942 for further details and a derivation. The $W_{e^{2 \theta t}-1}$ may appear confusing, but it's just a Wiener process with zero mean and variance $e^{2 \theta t}-1$. To calculate this in Matlab:

th = 1;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:2;             % Time vector
x0 = 0;                 % Set initial condition
rng(1);                 % Set random seed
W = zeros(1,length(t)); % Allocate integrated W vector
for i = 1:length(t)-1
    W(i+1) = W(i)+sqrt(exp(2*th*t(i+1))-exp(2*th*t(i)))*randn;
end
ex = exp(-th*t);
x = x0*ex+mu*(1-ex)+sig*ex.*W/sqrt(2*th);
figure;
plot(t,x);

This can be implemented without a for loop using cumsum and diff:

th = 1;
mu = 1.2;
sig = 0.3;
dt = 1e-2;
t = 0:dt:2;      % Time vector
x0 = 0;          % Set initial condition
rng(1);          % Set random seed
ex = exp(-th*t);
x = x0*ex+mu*(1-ex)+sig*ex.*cumsum([0 sqrt(diff(exp(2*th*t)-1)).*randn(1,length(t)-1)])/sqrt(2*th);
figure;
plot(t,x);

4. Resources

You can also use my own SDETools Matlab toolbox on GitHub for numerically solving SDEs and computing analytical solutions of common stochastic processes. In particular, see the sde_ou function to calculate analytical solutions for the Ornstein–Uhlenbeck process.

I also recommend reading the following excellent article for further details on SDEs and simulating them in Matlab:

Desmond J. Higham, 2001, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev. (Educ. Sect.), 43 525–46. http://dx.doi.org/10.1137/S0036144500378302

The URL to the Matlab files in the paper won't work – they can be found here now. Note, however, that some of the Matlab syntax (particularly related to random number generation and seeding) is a bit outdated as this was written nearly 15 years ago.

Answer 2

We should focus on simulating $\int_0^t e^{\theta (s)}\, \mathrm{d}W_s$. If $f(s)=e^{\theta (s)}$ is continuously differentiable, you could use the fact that $$\sum_{i=1}^{[tn]}f(s_i^*)\Big(W(s_i)-W(s_{i-1})\Big)\to\int_0^tf(s)dW(s)$$ in quadratic mean, for $s_i^* \in [s_{i-1},s_i]$. Note that you should use $$W(s_i)-W(s_{i-1}) \sim N(0,s_i-s_{i-1})$$ and that these increments are independent.