# Implementing Ornstein–Uhlenbeck in Matlab

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

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.

• @student441: I'm not sure why you'd want to do that. I also think you mean $\mathrm{d}W_s$, not $ds$. The integrand, $\frac{\sigma e^{\theta s}}{\sqrt{2 \theta}} \mathrm{d}W_s$, is already an integrated solution – you'd be integrating this again. To do it numerically, you can look at #2 above (the for loop case is simpler to adapt). In this case you have Wiener increments with $0$ mean and variance $e^{2 \theta t}-1$ (instead of $0$ mean and variance $t$). The diffusion function only differs by a sign in the exponent. Be careful if your exponential terms get very large or small. May 31, 2015 at 23:22
• I am simulating your part1 and I agree completely on the procedure. However plotting seems to show the process has much higher variance than the plot in the wikipedia article, regardless of the realisation of the process. Would you be able to explain why? Are the parameters given in the wikipedia article possibly wrong? Nov 21, 2017 at 21:36
• @semola: Thanks for pointing this out. I looked into it and it appears that the plot on Wikipedia is not being calculated properly. The effective standard deviation being used there is $\sigma^2$ rather than just $\sigma$. This can be seen by examining the "R-Quelltext" source code section for the original image. The code sigma*rnorm(n=1,mean=0,sd=sqrt(sigma^2*(to-from)/steps) should either be sigma*sqrt((to-from)/steps)*rnorm(n=1,mean=0,sd=1) or rnorm(n=1,mean=0,sd=sqrt(sigma^2*(to-from)/steps) or another equivalent form. Nov 22, 2017 at 13:45
• The plot also uses different time step (twice as long) than my code above, but that doesn't impact the variance of the process. Nov 22, 2017 at 13:46
• @horchler If we assume that $x_0=0$ and $\mu=0$ then $x_t = \sigma \frac{\sqrt{1-e^{-2 \theta t}}}{\sqrt{2 \theta t}} W_t$, correct? So, to simulate the paths, I first simulate a path of $W_t$ and then scale it by $\sigma \frac{\sqrt{1-e^{-2 \theta t}}}{\sqrt{2 \theta t}}$, right? This approach, however, seems to produce different results to the one outlined in your post and I am not quite sure why this is. Could you help explain why this might be the case? Thanks. Jun 7, 2018 at 9:46

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.

• Thanks, but can you elaborate some more? Do I understand correctly that the equation you would implement is: $x_t = x_0 e^{-\theta t} + \mu(1-e^{-\theta t}) + \sigma \sum_{i=1}^{[tn]}f(s_i^*)\Big(W(s_i)-W(s_{i-1})\Big) e^{\theta (s)}\,$ ? Can you please explain a bit more? May 18, 2015 at 15:10
• you are welcome, the solution below by horchler is actually quite extensive. May 19, 2015 at 9:12
• Hey, I know this was a long time ago but could you answer the question @student441 asked? Essentially, are you proposing 𝑥𝑡=𝑥0𝑒−𝜃𝑡+𝜇(1−𝑒−𝜃𝑡)+𝜎∑[𝑡𝑛]𝑖=1𝑓(𝑠∗𝑖)(𝑊(𝑠𝑖)−𝑊(𝑠𝑖−1))𝑒𝜃(𝑠) ? Mar 9, 2020 at 0:04