Simulating systems of SDEs

Setup

Assume that a general two dimensional system of SDE is given by

σσμ

σσμ

++=

++=

)(),,()(),,(),,(

)(),,()(),,(),,(

tdztYXtdztYXdttYXdY

tdztYXtdztYXdttYXdX

2121111

ttttttt

2221212

ttttttt

With starting values

)0(

)0(

=

=

YY

XX

where

X

0

0

0

and Y

could possibly be stochastic variables, we will here assume that they

are real numbers.

The Wiener processes z

0

1

and z

are possibly correlated according to

dtdzdzE ρ=][

21

2

where we here assume that ρ is a real constant (which cannot have absolute value

greater than 1, by definition of correlation).

The Euler-Maryama Method

The Euler-Maryama Method is a basic method of simulating SDEs. The idea is

simple. It is in essence the same as the Euler method for ODEs, that is the simulation

starts at the starting points and walks stepwise to the endpoint. The important

difference here is how the noise term dz is interpreted in this situation.

Say that you want to simulate your system of SDEs on an interval [0,T]. First the

interval is discretized into N (positive integer) number of time steps, ie Δt = T/N. We

will only consider equally spaced time steps. The interpretation of the noise term dz is

due to Ito stochastic integral in which the forward increment of the standard Wiener

process is considered. Regarding our simulation the noise increment is simulated by

Δ=Δ

)1,0()( NisZwheretZtz

ji

(For more information about Ito integral go to any book treating stochastic calculus.)

The system is therefore discretisized according to

Δ+Δ+Δ=Δ

Δ+Δ+Δ=Δ

σσμ

σσμ

)(),,()(),,(),,(

)(),,()(),,(),,(

tztYXtztYXttYXY

tztYXtztYXttYXX

2121111

ttttttt

2221212

ttttttt

We will here denote the simulated values of X(t

j

) and Y(t

j

) as x

j

and y

. The iteration

scheme will look like

YYXXtSet

ttt

loopNjfor

jj

1

)1,0()(

tztyxtztyxttyxxx

NisZwheretZtz

ji

−=

===

jjjjjjjjjjjjj

tztyxtztyxttyxyy

jjjjjjjjjjjjj

loopend

1…0

)0(,)0(,0

000

Δ=Δ

Δ+=

+

+

Δ+Δ+Δ+=

Δ+Δ+Δ+=

)(),,()(),,(),,(

)(),,()(),,(),,(

21211111

+

σσμ

σσμ

22212121

Here it is important to recall that when generating random numbers should be

adjusted to incorporate possible correlation. Normally a computer generates

independent variables and in that case an easy way of transforming independent

random variables Δw

1

and Δw

2

into correlated Δz

is

Δ−+Δ=Δ

Δ=Δ

1 wwz

wz

11

ρρ

12

2

2

1

and Δz

2

One simulation only give you one possible path on [0,T] and so doesn’t yield much

information. To be able to extract some useful information about the statistical

properties of the model a large number of simulations are made, we will denote the

number of simulations by M,

j

Example with Matlab code: A stochastic volatility model

Setup of the model

As an example of a two dimensional system of SDE consider the stochastic volatility

model given by

+−=

⋅+−=

)0(

)0(

υκ

dzrdtrrdr

υυ

)2()(

)1()(

υσυυκυ

rr

dzdtd

=

=

rr

υυυ

0

0

In this model r is the short rate and

υ

is the instantaneous variance of the short rate

process,

υ will then of course be the instantaneous correlation. Equation (1) can

therefore be viewed as

υσσκ =+−=

dzrdtrrdr ,)(

rrrr

Where

υσ =

and υ is stochastic according to SDE (2).

r

A little bit about Matlab

Matlab is an interactive environment for numerical computing. The programming in

Matlab is done is what is referred to as m-files, which is textfiles with the extension

.m. There are two types of M-files, script-files and functions. The most important

difference is how they treat variables. Script-files use global variables defined in the

workspace. A function-file has local variables and therefore input and output

parameters need to be specified. There is no need to compile m-files since Matlab

interpret the code (m-files can compiled by using a separate toolbox). Running scriptfiles

are the

same

as writing

the

commands

directly

in

the

command

window.

Executing

a script-file

is done

by

having

the

current

directory

containing

the

file

and

typing

the name

of the

file

(without

the

extension

.m)

in

the

command

window.

The

example

program

containing

the

volatility

model

is a

script-file.

The

file the

stochastic

volatility

model

is

called ‘Simulating_SV_model’.

This file

contains

the

programme

and

sufficient comments

so that

the

code

should be

reasonably

easy

to

understand

for

someone

used

to reading

code

in any

common

language.

A

few

comments

regarding

the

syntax

might

be useful

though

randn(‘state’,101)

Sets the state for the random number generator. This means

that the same random numbers will be generated every time

the simulation is repeated for the same state.

randn(N,M)

Generates a NxM random matrix with independent N(0,1)

variables.

A(X,:)

Given a matrix A the command A(N,:) = X will give the

value X to all elements on row N.

zeros(N,M)

Creates a NxM matrix with all elements equal zero. In

Matlab it is important to preallocate arrays for efficiency.

for m = 1:M

For loop, i.e. makes a loop repeated M times with the value

m will have the value of the current loop. Suitable in

situations where a vector or matrix is indexed within the

loop, see code for example of this.

A note on computational efficiency in Matlab. For loops is not the most

computationally efficient way of doing simulations in MATLAB. It’s actually quite

slow, vectorized code is usually much faster in MATLAB. But the code is often easier

to read with loops rather than vectorized code and that is the reason it has been kept

that way in the example here. Especially those who are used to read code from some

other common language such as C/C++, FORTRAN etc would normally find loops

easier to read. Often it is possible to vectorize code in situations like this simulation.

Matlab code for the stochastic volatility model

% ——————————–

% – Simulating_SV_model –

% ———————————

%

% This program shows the simulated distribution of the stochastic

% volatility model given by:

%

% dr = K(r_bar – r(t))dt + sqrt(v(t)*r(t))*dz1

% r(0) = r0

%

% dv = K_sig(v_bar – v(t))dt + sig_v*sqrt(v(t))*dz2

% v(0) = v0

%

% E[dz1*dz2] = rho*dt

%

% This program is now setup for graphing process r(t) and

% volatility process v(t) respectively and their simulated distributions.

% USER DEFINED: General simulation parameters

T = 1; % End time

N = 200; % Number of time steps

M = 10000; % Number of simulations

randn(‘state’,101) % Setting the seed for the random number generator

% USER DEFINED: Model parameters for short rate processes r(t)

r0 = 0.03; % Starting value

K = 0.6;

r_bar = 0.07; % Long run level

% USER DEFINED: Parameter for volatility process v(t) and correlation coeff

v0 = 0.01; % Starting value

K_sig = 2;

v_bar = 0.01; % Long run level

sig_v = 0.03; % Standard deviation

rho = 0.5; % Instantaneous correlation

% Setting parameters, generate random numbers

dt = T/N;

sqrtdt = sqrt(dt);

dw_setup = randn(N,2*M); % Generates normally distributed RVs

% Transforming independent random varibales to correlated

dz_1 = dw_setup(:,1:M);

dz_2 = rho*dw_setup(:,1:M) + sqrt(1-rho^2)*dw_setup(:,M+1:2*M);

% Simulating the volatility SDE volatility and short rate

% Allocating memory space for SDE v(t) and r(t) paths

% Each column represents one simulated path

sde_v = zeros(N+1,M);

sde_r = zeros(N+1,M);

% Setting initial value for each simulation in first row

sde_v(1,:) = v0;

sde_r(1,:) = r0;

%Generates the M paths for volatility and short rate v(t) and r(t)

for m = 1:M

for n = 1:N

v = sde_v(n,m); % Set new variable for readability in loop

r = sde_r(n,m); % Set new variable for readability in loop`

sde_v(n+1,m) = v + K_sig*(v_bar-v)*dt+sig_v*sqrt(v)*dz_2(n,m)*sqrtdt;

sde_r(n+1,m) = r + K*(r_bar-r)*dt+sqrt(v*r)*dz_1(n,m)*sqrtdt;

end;

end;

Example output from the simulated stochastic volatility model

r(t)

0.2

Sample paths of r(t)

0.15

0.1

0.05

0 0.5 1

0

t

v(t)

Freq

150

Sample distribution of r(t)

100

50

0

0 0.05 0.1 0.15 0.2

r(t)

Freq

0.02

Sample paths of v(t)

0.015

0.01

0.005

0 0.5 1

0

100

50

0

t

150

Sample distribution of v(t)

0.005 0.01 0.015 0.02

v(t )