Numerical simulation of SDEs#
This is a reproduction of certain scripts found in Higham, An algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations.[Higham, 2001] This paper is an accessible introduction to SDEs, which is centered around ten scripts. Below are reproductions of these scripts (excluding two on linear stability) and some supplementary notes.
Why stochastic differential equations#
We are often interested in modelling a system whose state takes values in a continuous range, and over a continuous time domain. Whereas ordinary differential equations (ODEs) describe variables which change according to a deterministic rule, SDEs describe variables whose change is governed partly by a deterministic component and partly by a stochastic component. SDEs are therefore an appropriate model for systems whose dynamics involve some true randomness, or some fine grained complexity which we cannot afford to or do not wish to model.
The Wiener process#
In order to define the stochastic component of the transition rule of a stochastic system, we must define an appropriate noise model. The Wiener process is a stochastic process that is commonly used for this purpose.
Definition 146 (Wiener process)
A standard Wiener process over [0, T] is a random variable
with probability 1.For
the random variable has distribution .For
then are independent.
We can imagine the Wiener process as the path followed by a particle that experiences infinitely many, infinitesimal kicks. The size of these kicks,
Sampling from a Wiener process#
Before using the Wiener process to define an SDE, let’s look at the process itself. How can we draw a sample
where
Show code cell source
# Set random seed
np.random.seed(0)
# Integration parameters
T = 1
N = 500
dt = T / N
# Times to sample at
t = np.linspace(0, T, N)[:, None]
# Sample dW's and compute cumulative sum
dW = dt ** 0.5 * np.random.normal(size=(N - 1, 1))
W = np.concatenate([[0], np.cumsum(dW)])
So even though we can’t represent the entirety of the path, we can sample it to arbitrary precision.
Function of a Wiener process#
Suppose we are interested in a stochastic process which is a function of a Wiener process, say
In this case we can sample S = 1000
samples of
Show code cell source
# Set random seed
np.random.seed(0)
# Time to simulate for, discretisation level and number of paths
T = 1
N = 500
S = 1000
dt = T / N
# Times to sample at
t = np.linspace(0, T, N)
# Sample dW's and compute cumulative sum
dW = dt ** 0.5 * np.random.normal(size=(N - 1, S))
W = np.concatenate(
[np.zeros(shape=(1, S)), np.cumsum(dW, axis=0)],
axis=0,
)
# Compute values
X = np.exp(t[:, None] + 0.5 * W)
X_mean = np.mean(X, axis=1)
X_stdev = np.var(X, axis=1) ** 0.5
Evaluating a stochastic integral#
The next thing that we look at is the evaluation of a stochastic integral. Let’s consider the integral
where
where
does not depend on the choice of
where we have defined
is equal to the following values under the Ito and Stratonovich definitions:
Proof: Evaluating under the Ito and Stratonovich integrals
If we use the Ito integral, we have
The distribution of
By contrast, if we use the Stratonovich integral, we have
where
and observing that the distribution above has the form of a product of normal distributions over
Since
Ito integral approximation (exact): 1.419 (1.415)
Stratonovich approximation (exact): 1.917 (1.915)
The difference between the Ito and Stratonovich integrals does not vanish as
Euler-Maruyama method#
The Euler-Maruyama method is the analoge of the Euler method for deterministic integrals, applied to the stochastic case.
Definition 147 (Euler-Maruyama method)
Given a scalar SDE with drift and diffusion functions
the Euler-Maruyama method approximates
where
Let’s look at an implementation of the Euler-Maruyama (EM) method.
The euler_maruyama
below function takes the drift and diffusion functions dW
in advance.
def euler_maruyama(seed, X0, T, N, f, g):
# Set random seed
np.random.seed(seed)
# Set discretisation, initial value and times at which to evaluate
dt = T / N
X = [X0]
t = np.linspace(0, T, N + 1)
# Sample Wiener process dW's
dW = dt ** 0.5 * np.random.normal(size=(N,))
for i in range(N):
# Apply Euler-Maruyama at each point in time
dX = f(X[-1], t[i]) * dt + g(X[-1], t[i]) * dW[i]
# Store the new X
X.append(X[-1] + dX)
# Compute W to return it at the end
W = np.concatenate([[0], np.cumsum(dW)])
return t, X, W
Below is the definition of f
and g
that we will be integrating, namely
known as the Black-Scholes model. This is implemented as a closure, i.e. f_g_black_scholes
takes in the appropriate lambda
and mu
and returns the corresponding f
and g
.
def f_g_black_scholes(lamda, mu):
def f(X, t):
return lamda * X
def g(X, t, grad=False):
return mu if grad else mu * X
return f, g
We choose these drift and diffusion terms because the associated SDE has a closed form solution, with which we can compare our numerical solution. The analytic solution to the Black-Scholes model is
which we implement in exact_black_scholes
below. Unlike ODEs, whose solution is a unique function, the solution of an SDE depends on the random noise sample
def exact_black_scholes(X0, t, W, lamda, mu):
return X0 * np.exp((lamda - 0.5 * mu ** 2) * t + mu * W)
We are can now run the EM method and compare it aginst the exact solution below.
We can change the accuracy of the solution by adjusting
In the case
def f_g_sine(omega, mu):
def f(X, t):
return omega * np.cos(omega * t)
def g(X, t):
return mu * X
return f, g
Strong and weak convergence#
Since the choice of the number of bins
Definition 148 (Strong convergence)
A method for approximating a stochastic process
for any fixed
Strong convergence refers to the rate of convergence of the approximation
Definition 149 (Weak convergence)
A method for approximating a stochastic process
for any fixed
The paper states without proof that, under conditions on
Milstein’s higher order method#
Just as higher order methods for ODEs exist for obtaining refined estimates of the solution, so do methods for SDEs, such as Milstein’s higher order method.
Definition 150 (Milstein’s method)
Given a scalar SDE with drift and diffusion functions
the Milstein method approximates
where
Milstein’s method achieves a strong convergence rate of
Let’s first use Milstein’s method to get a single solution of the Black-Scholes model, as we did for EM.
We can also look at the strong and weak rates of convergence for Milstein’s method.
Milstein’s method achieves a strong convergence rate of
Stochastic chain rule#
Suppose we want to evaluate a function
This way, we could solve the following ODE directly for
For an autonomous SDE however the chain rule takes a different form, which under the Ito definition is as follows.
Theorem 126 (Ito’s result for one dimension)
Let
where
If
For a more formal definition and proof of Ito’s result see [Oksendal, 1992] (Theorem 4.1.8 and pages 44-48). Below is a short sketch proof, which highlights why the additional term appears in the formula.
Informal argument: Ito’s result for one dimension
Writing the infinitesimal difference in
where the
In the above expression, we can neglect the first two terms which are of order
Now consider the contribution to the Ito integral of the last term, i.e. the sum
where
With this corrected rule, we can directly integrate the stochastic process for
References#
D.J. Higham. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review, 43(3):525–546, 2001.
Bernt Oksendal. Stochastic Differential Equations (3rd Ed.): An Introduction with Applications. Springer-Verlag, Berlin, Heidelberg, 1992. ISBN 3387533354.