Stochastic Differential Equations: Lab 1

This background for these exercises is article of D Higham, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Review 43:525-546 (2001). Higham provides Matlab codes illustrating the basic ideas at http://personal.strath.ac.uk/d.j.higham/algfiles.html, which are also given in the paper.

For random processes in python you should look at the numpy.random module. To set the initial seed (which you should not do in a real simulation, but allows for reproducible testing), see numpy.random.seed.

Brownian processes

A random walk or Brownian process or Wiener process is a way of modelling error introduced by uncertainty into a differential equation. The random variable representing the walk is denoted $W$. A single realization of the walk is written $W(t)$. We will assume that

1. The walk (value of $W(t)$) is initially (at $t=0$) $0$, so $W(0)=0$, to represent "perfect knowledge" there;
2. The walk is on average zero, so $\mathbb{E}[W(t+h) - W(t)] = 0$, where the expectation value is $$ \mathbb{E}[W] = \int_{-\infty}^{\infty} t W(t) \, \text{d}t $$
3. Any step in the walk is independent of any other step, so $W(t_2) - W(t_1)$ is independent of $W(s_2) - W(s_1)$ for any $s_{1,2} \ne t_{1,2}$.

These requirements lead to a definition of a discrete random walk: given the points $\{ t_i \}$ with $i = 0, \dots, N$ separated by a uniform timestep $\delta t$, we have - for a single realization of the walk - the definition $$ \begin{align} \text{d}W_i &= \sqrt{\delta t} {\cal N}(0, 1), \\ W_i &= \left( \sum_{j=0}^{i-1} \text{d}W_j \right), \\ W_0 &= 0 \end{align} $$ Here ${\cal N}(0, 1)$ means a realization of a normally distributed random variable with mean $0$ and standard deviation $1$: programmatically, the output of numpy.random.randn.

When working with discrete Brownian processes, there are two things we can do.

1. We can think about a single realization at different timescales, by averaging over more points. E.g. $$ W_i = \left( \sum_{j=0}^{i_1} \sum_{k=0}^{p} \text{d}W_{(p j + k)} \right) $$ is a Brownian process with timestep $p \, \delta t$.
2. We can think about multiple realizations by computing a new set of steps $\text{d}W$, whilst at the same timestep.

Both viewpoints are important.

Tasks

1. Simulate a single realization of a Brownian process over $[0, 1]$ using a step length $\delta t = 1/N$ for $N = 500, 1000, 2000$. Use a fixed seed of 100. Compare the results.
2. Simulation different realizations of a Brownian process with $\delta t$ of your choice. Again, compare the results.

Evaluate the function $u(W(t)) = \sin^2(t + W(t))$, where $W(t)$ is a Brownian process, on $M$ Brownian paths for $M = 500, 1000, 2000$. Compare the average path for each $M$.

The average path at time $t$ should be given by $$ \begin{equation} \int_{-\infty}^{\infty} \frac{\sin(t+s)^2 \exp(-s^2 / 2t)}{\sqrt{2 \pi t}} \,\text{d}s. \end{equation} $$

Stochastic integrals

We have, in eg finite elements or multistep methods for IVPs, written the solution of differential equations in terms of integrals. We're going to do the same again, so we need to integrate random variables. The integral of a random variable with respect to a Brownian process is written $$ \int_0^t G(s) \, \text{d}W_s, $$ where the notation $\text{d}W_s$ indicates that the step in the Brownian process depends on the (dummy) independent variable $s$.

We'll concentrate on the case $G(s) = W(s)$, so we're trying to integrate the Brownian process itself. If this were a standard, non-random variable, the answer would be $$ \int_0^t W(s) \, \text{d}W_s = \frac{1}{2} \left( W(t)^2 - W(0)^2 \right). $$

When we approximate the quadrature numerically than we would split the interval $[0, T]$ into strips (subintervals), approximate the integral on each subinterval by picking a point inside the interval, evaluating the integrand at that point, and weighting it by the width of the subinterval. In normal integration it doesn't matter which point within the subinterval we choose.

In the stochastic case that is not true. We pick a specific point $\tau_i = a t_i + (1-a) t_{i-1}$ in the interval $[t_{i-1}, t_i]$. The value $a \in [0, 1]$ is a constant that says where within each interval we are evaluating the integrand. We can then approximate the integral by

\begin{equation} \int_0^T W(s) \, dW_s = \sum_{i=1}^N W(\tau_i) \left[ W(t_i) - W(t_{i-1}) \right] = S_N. \end{equation}

Now we can compute (using that the expectation of the products of $W$ terms is the covariance, which is the minimum of the arguments)

\begin{align} \mathbb{E}(S_N) &= \mathbb{E} \left( \sum_{i=1}^N W(\tau_i) \left[ W(t_i) - W(t_{i-1}) \right] \right) \\ &= \sum_{i=1}^N \mathbb{E} \left( W(\tau_i) W(t_i) \right) - \mathbb{E} \left( W(\tau_i) W(t_{i-1}) \right) \\ &= \sum_{i=1}^N (\min\{\tau_i, t_i\} - \min\{\tau_i, t_{i-1}\}) \\ &= \sum_{i=1}^N (\tau_i - t_{i-1}) \\ &= (t - t_0) a. \end{align}

The choice of evaluation point matters.

So there are multiple different stochastic integrals, each (effectively) corresponding to a different choice of $a$. The two standard choices are There are two standard choices of stochastic integral.

1. Ito: choose $a=0$.
2. Stratonovich: choose $a=1/2$.

These lead to $$ \int_0^t G(s) \, \text{d}W_s \simeq_{\text{Ito}} \sum_{j=0}^{N-1} G(s_j, W(s_j)) \left( W(s_{j+1}) - W(s_j) \right) = \sum_{j=0}^{N-1} G(s_j) \text{d}W(s_{j}) $$ for the Ito integral, and $$ \int_0^t G(s) \, \text{d}W_s \simeq_{\text{Stratonovich}} \sum_{j=0}^{N-1} \frac{1}{2} \left( G(s_j, W(s_j)) + G(s_{j+1}, W(s_{j+1})) \right) \left( W(s_{j+1}) - W(s_j) \right) = \sum_{j=0}^{N-1} \frac{1}{2} \left( G(s_j, W(s_j)) + G(s_{j+1}, W(s_{j+1})) \right) \text{d}W(s_{j}). $$ for the Stratonovich integral.

Tasks

Write functions to compute the Itô and Stratonovich integrals of a function $h(t, W(t))$ of a given Brownian process $W(t)$ over the interval $[0, 1]$.

Test the functions on $h = W(t)$ for various $N$. Compare the limiting values of the integrals.

Euler-Maruyama's method

Now we can write down a stochastic differential equation.

The differential form of a stochastic differential equation is $$ \frac{\text{d}X}{\text{d}t} = f(X) + g(X) \frac{\text{d}W}{\text{d}t} $$ and the comparable (and more useful) integral form is $$ \text{d}X = f(X) \, \text{d}t + g(X) \text{d}W. $$ This has formal solution $$ X(t) = X_0 + \int_0^t f(X(s)) \, \text{d}s + \int_0^t g(X(s)) \, \text{d}W_s. $$

We can use our Ito integral above to write down the Euler-Maruyama method

$$ X(t+h) \simeq X(t) + h f(X(t)) + g(X(t)) \left( W(t+h) - W(t) \right) + {\cal{O}}(h^p). $$

Written in discrete, subscript form we have

$$ X_{n+1} = X_n + h f_n + g_n \, \text{d}W_{n} $$

The order of convergence $p$ is an interesting and complex question.

Tasks

Apply the Euler-Maruyama method to the stochastic differential equation

$$ \begin{equation} dX(t) = \lambda X(t) + \mu X(t) dW(t), \qquad X(0) = X_0. \end{equation} $$

Choose any reasonable values of the free parameters $\lambda, \mu, X_0$.

The exact solution to this equation is $X(t) = X(0) \exp \left[ \left( \lambda - \tfrac{1}{2} \mu^2 \right) t + \mu W(t) \right]$. Fix the timetstep and compare your solution to the exact solution.

Vary the timestep of the Brownian path and check how the numerical solution compares to the exact solution.

Convergence

We have two ways of thinking about Brownian paths or processes.

We can fix the path (ie fix $\text{d}W$) and vary the timescale on which we're looking at it: this gives us a single random path, and we can ask how the numerical method converges for this single realization. This is strong convergence.

Alternatively, we can view each path as a single realization of a random process that should average to zero. We can then look at how the method converges as we average over a large number of realizations, also looking at how it converges as we vary the timescale. This is weak convergence.

Formally, denote the true solution as $X(T)$ and the numerical solution for a given step length $h$ as $X^h(T)$. The order of convergence is denoted $p$.

Strong convergence

$$ \mathbb{E} \left| X(T) - X^h(T) \right| \le C h^{p} $$

For Euler-Maruyama, expect $p=1/2$!.

Weak convergence

$$ \left| \mathbb{E} \left( \phi( X(T) ) \right) - \mathbb{E} \left( \phi( X^h(T) ) \right) \right| \le C h^{p} $$

For Euler-Maruyama, expect $p=1$.

Tasks

Investigate the weak and strong convergence of your method, applied to the problem above.