Intereting Posts

Induction: $\sum_{k=0}^n \binom nk k^2 = n(1+n)2^{n-2}$
In how many ways can you select one of the two but not both?
Does every inner product space has an orthogonal basis?
How could we define the factorial of a matrix?
Proof of analogue of the Cauchy-Schwarz inequality for random variables
Inverting formal power series wrt. composition
$fg\in L^1$ for every $g\in L^1$ prove $f\in L^{\infty}$
$\int_{-\infty}^{+\infty} e^{-x^2} dx$ with complex analysis
If $\displaystyle \lim _{x\to +\infty}y(x)\in \mathbb R$, then $\lim _{x\to +\infty}y'(x)=0$
Prove that the only homomorphism between two cyclic groups with distinct prime orders is the trivial one
Proof – Square Matrix has maximal rank if and only if it is invertible
Polynomials over finite fields
The adjoint of finite rank operator is finite rank
When to learn category theory?
Intersection and union of simply connected domain

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…

- Strong Markov property of Bessel processes
- Product rule with stochastic differentials
- Brownian motion - Hölder continuity
- fractional Brownian motion is not a semimartingale. How to apply Ergodic theorem in the proof of this theorem?
- To confirm the Novikov's condition
- Expected value of average of Brownian motion

Can someone please help me out?

I am new to stochastic calculus, so please help me understand step by step.

- How to compute elliptic integrals in MATLAB
- Robust Control VS Optimal Control
- Fast Matlab Code for hypergeometric function $_2F_1$
- fractional Brownian motion is not a semimartingale. How to apply Ergodic theorem in the proof of this theorem?
- Brownian motion introduction
- Expectation of Stopping Time w.r.t a Brownian Motion
- Hitting times for Brownian motions
- Convergence of exponential Brownian martingale to zero almost surely
- Brownian motion martingale
- Brownian Motion Conditional Expectation Question

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). 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.),43525–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.

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.

- How to solve a bivariate quadratic (not necessarily Pell-type) equation?
- Proofs that involve Tricks
- Why if $f'$ is unbounded, then $f$ isn't uniformly continuous?
- How to prove $\int_0^1\tan^{-1}\left\frac{dx}{x}=\frac{\pi}{8}\ln\frac{\pi^2}{8}?$
- 4 Element abelian subgroup of S5.
- Can the same subset be both open and closed?
- Prove that $f(1999)=1999$
- Word Question Involving the Definition of the Derivative
- How to Visualize points on a high dimensional (>3) Manifold?
- Which of the following is not true?
- Monte Carlo double integral over a non-rectangular region (Matlab)
- Prove that ${p_1},\cdots ,\sqrt{p_n}):\mathbb{Q}]=r^n$
- How to prove that $\lim_{x\rightarrow \infty}\dfrac{x^2}{e^x}=0$?
- How to get started on this statistics problem?
- Number of birational classes of dimension d, geometric genus 0 varieties?