Simulating systems of SDEs
WE WRITE PAPERS FOR STUDENTS
Tell us about your assignment and we will find the best writer for your project.
Get Help Now!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 )
