Order 1.5 strong Taylor scheme for 2D Geometric
Brownian Motions
Contents
• Introdution
• Algorithm
• Implementation
• Summary
• Subroutines
Introdution
This project demonstrates simulation of a 2D geometric Brownian motions with
a constant drift vector and correlation matrix using an order 1.5 strong Taylor
scheme. We consider the geometric Brownian Motion (Black-Scholes Model)
dX1(t) = a1X
1(t)dt + b11X
1(t)dW 1(t) + b12X
1(t)dW 2(t)
dX2(t) = a2X
2(t)dt + b21X
2(t)dW 1(t) + b22X
2(t)dW 2(t)
Here W 1(t), W 2(t) are independent standard Brownian motions. The exact
solution is given by
X1(t) = X1(0) exp((a1 − 1
2
(b211 + b
2
12))t + b11W
1(t) + b12W
2(t)).
X2(t) = X2(0) exp((a2 − 1
2
(b221 + b
2
22))t + b21W
1(t) + b22W
2(t)).
Algorithm
We simulate the geometric Brownian motion using order 1.5 strong Taylor
scheme, on a uniform numerical grid of time of step ∆t, {tn = n∆t, n =
0, 1, ..., N}, with T = N∆t. To simplify notation without loss of generality,
we consider the one-step simulation from time 0 to time t. The Wagner-Platen
1
expansion of X1(t) is given by
X1(t) = X1(0)× (
1 + a1t + b11W
1(t) + b12W
2(t)
+b211I(1, 1) + b11b12(I(2, 1) + I(1, 2)) + b
2
12I(2, 2)
+a21I(0, 0) + a1b11(I(1, 0) + I(0, 1)) + a1b12(I(2, 0) + I(0, 2))
+b311I(1, 1, 1) + b
2
11b12(I(2, 1, 1) + I(1, 2, 1) + I(1, 1, 2))
+b11b
2
12(I(2, 1, 2) + I(1, 2, 2) + I(2, 2, 1)) + b
3
12I(2, 2, 2))
+O(t2)
following (Platen & Bruti-Liberati 2010). Here
I(i, j) =
∫ t
0
W i(s)dW j(s)
I(i, j, k) =
∫ t
0
∫ s
0
W i(u)dW j(u)dW k(s)
are multiple regular or Ito integrals, for i, j, k = 0, 1, 2, where W 0(t) = t for
simplicity. Simulation scheme for X2 is analogous.
Following Lemma 4.2.3 of (Platen & Bruti-Liberati 2010), the relevant multiple
integrals or sums of multiple integrals are given explicitly as
I(0, 0) = 12 t
2
I(1, 1) = 12W
1(t)2 − 12 t
I(2, 2) = 12W
2(t)2 − 12 t
I(1, 0) + I(0, 1) = tW 1(t)
I(2, 0) + I(0, 2) = tW 2(t)
I(2, 1) + I(1, 2) = W 1(t)W 2(t)
I(1, 1, 1) = 16W
1(t)3 − 12 tW 1(t)
I(2, 2, 2) = 16W
2(t)3 − 12 tW 2(t)
I(1, 1, 2) + I(2, 1, 1) + I(1, 2, 1) = 12 (W
1(t)2 − t)W 2(t)
I(2, 2, 1) + I(2, 1, 2) + I(1, 2, 2) = 12 (W
2(t)2 − t)W 1(t)
2
thanks to commutativity of the 2D GBM. For example
I(1, 2) =
∫ t
0
W 1(s)dW 2(s)
= W 1(s)W 2(s)|t0 −
∫ t
0
W 2(s)dW 1(s)
= W 1(t)W 2(t)− I(2, 1),
implying I(1, 2) + I(2, 1) = W 1(t)W 2(t).
It is worth mentioning that without symmetry, not all relevant multiple integrals
can be obtained as simple functions of increments of the Brownian motions. In
this case the multiple integrals may be represented by Stratonovich integrals and
approximated by the trapezoidal rule over [0, t] on a finer grid, as in contrast to
the Ito integrals, which should be approximated using a rectangular rule with
lower accurracy. For example:
I(2, 1, 1) + I(1, 2, 1) =
∫ t
0
W 1(s)W 2(s)dW 1(s)
=
∫ t
0
W 1(s)W 2(s) ◦ dW 1(s)− 12
∫ t
0
W 2(s) ◦ ds
where the second integral with respect to time can also be understood as a
Stratonovich integral. To simulate the Stratonovich integrals on a numerical grid
from tn to tn+1, take tn as 0, tn+1 as t, simulate incremets of W
1 and W 2 over
[0, t] on a finer sub-grid over the time step, and approximate the Stratonovich
integrals over [0,t] using the trapezoidal rule. For the 2D GBM we consider
here, this is not needed. Nevertheless, simulations using this approximations
are included in the following demonstrations as a comparison.
Implementation
function GBMTaylor15()
Simulation with Commutativity
This simulation makes full use of commutativity.
rng(500);
simErrPlot(0,[1;1],1,[10,50,100,200,400],@simfunCom,100);
3
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
1
2
3
4
5
6
7
∆ t = 0.1
sim
sim
exact
exact
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
∆ t = 0.02
sim
sim
exact
exact
4
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
∆ t = 0.01
sim
sim
exact
exact
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.2
0.4
0.6
0.8
1
1.2
1.4
∆ t = 0.005
sim
sim
exact
exact
5
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.8
1
1.2
1.4
1.6
1.8
2
∆ t = 0.0025
sim
sim
exact
exact
log(∆ t)
-6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 -2
lo
g(E
rro
r)
-12
-11
-10
-9
-8
-7
-6
-5
Fitted slope=1.5146
error
fitted
6
Simulation with Trapezoidal Approximation
This simulation makes some use of commutativity and uses Stratonovich integral
and trapezoidal approximation.
rng(500);
simErrPlot(0,[1;1],1,[10,50,100,200,400],@simfunStr,100);
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
1
2
3
4
5
6
7
∆ t = 0.1
sim
sim
exact
exact
7
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
∆ t = 0.02
sim
sim
exact
exact
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
∆ t = 0.01
sim
sim
exact
exact
8
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.2
0.4
0.6
0.8
1
1.2
1.4
∆ t = 0.005
sim
sim
exact
exact
Time
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.8
1
1.2
1.4
1.6
1.8
2
∆ t = 0.0025
sim
sim
exact
exact
9
log(∆ t)
-6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 -2
lo
g(E
rro
r)
-12
-11
-10
-9
-8
-7
-6
-5
Fitted slope=1.5144
error
fitted
end
Summary
Order 1.5 convergence of the strong Taylor scheme was successfully achieved.
The trapezoidal approximations seem to be fairly accurate compared with the
closed-form simulations.
Subroutines
function order = simErrPlot(T0,X0,DT,NT,simfun,NS)
Nd = length(NT);
Dt = DT./NT;
Err = zeros(1,Nd);
Dim = size(X0,1);
for n = 1:Nd
XTS = zeros(Dim,NS);
XTE = zeros(Dim,NS);
for s = 1:NS
N = NT(n);
dt = Dt(n);
XS = zeros(Dim,N+1);
XS(:,1)=X0;
10
XE = XS;
t0 = T0;
for k = 1:N
[XS(:,k+1),XE(:,k+1)] = ...
simfun(t0,XS(:,k),XE(:,k),dt);
t0 = t0+dt;
end
XTS(:,s) = XS(:,N+1);
XTE(:,s) = XE(:,N+1);
end
figure(n); clf
plot((0:N)*dt,XS,’.-’); hold on
plot((0:N)*dt,XE,’-.’); hold off
xlabel(’Time’)
leg=[repmat({’sim’},1,Dim),...
repmat({’exact’},1,Dim)];
legend(leg{:})
title([’\Delta t = ’,num2str(dt)])
Err(n) = mean(sum(abs(XTS-XTE),1),2);
end
LDt = log(Dt);
LErr = log(Err);
fit = polyfit(LDt,LErr,1);
figure(Nd+1); clf
plot(LDt,LErr,’o-’); hold on
fLErr = fit(1)*LDt+fit(2);
plot(LDt,fLErr,’-.’); hold off
xlabel(’log(\Delta t)’)
ylabel(’log(Error)’)
legend(’error’,’fitted’)
title([’Fitted slope=’,num2str(fit(1))])
order = fit;
end
function [xts,xte] = simfunCom(~,x0s,x0e,dt)
ddt=1e-6;
ndt=round(dt/ddt);
ddt=dt/ndt;
ddw = sqrt(ddt)*randn(2,ndt);
dw = sum(ddw,2);
I11=0.5*dw(1)^2-0.5*dt;
I21I12 = dw(1)*dw(2);
I22=0.5*dw(2)^2-0.5*dt;
11
I00=0.5*dt^2;
I10I01=dt*dw(1);
I20I02=dt*dw(2);
I111 = dw(1)^3/6-0.5*dt*dw(1);
I222 = dw(2)^3/6-0.5*dt*dw(2);
I221I212I122=0.5*(dw(2)^2-dt)*dw(1);
I112I211I121=0.5*(dw(1)^2-dt)*dw(2);
xts=zeros(2,1);
xte=zeros(2,1);
a=0.1; b1=0.2; b2=0.3;
xts(1)=x0s(1)*(1+a*dt+b1*dw(1)+b2*dw(2)+ ...
b1^2*I11+b1*b2*I21I12+b2^2*I22+...
a^2*I00+a*b1*I10I01+a*b2*I20I02+...
b1^3*I111+b1^2*b2*I112I211I121+...
b1*b2^2*I221I212I122+...
b2^3*I222);
xte(1)=x0e(1)*exp(a*dt+b1*dw(1)+b2*dw(2)-...
0.5*(b1^2+b2^2)*dt);
a=0.2; b1=-0.4; b2=0.7;
xts(2)=x0s(2)*(1+a*dt+b1*dw(1)+b2*dw(2)+...
b1^2*I11+b1*b2*I21I12+b2^2*I22+...
a^2*I00+a*b1*I10I01+a*b2*I20I02+...
b1^3*I111+b1^2*b2*I112I211I121+...
b1*b2^2*I221I212I122+...
b2^3*I222);
xte(2)=x0e(2)*exp(a*dt+b1*dw(1)+b2*dw(2)-...
0.5*(b1^2+b2^2)*dt);
end
function [xts,xte] = simfunStr(~,x0s,x0e,dt)
ddt=1e-6;
ndt=round(dt/ddt);
ddt=dt/ndt;
ddw = sqrt(ddt)*randn(2,ndt);
cdw=[[0;0],cumsum(ddw,2)];
cdt=ddt*(0:ndt);
dw = cdw(:,end);
12
I11=0.5*dw(1)^2-0.5*dt;
I21=stratonovich(cdw(2,:),cdw(1,:));
I12=stratonovich(cdw(1,:),cdw(2,:));
I21I12 = I21+I12;
I22=0.5*dw(2)^2-0.5*dt;
I00=0.5*dt^2;
I10=stratonovich(cdw(1,:),cdt);
I01=stratonovich(cdt,cdw(1,:));
I10I01=I10+I01;
I20=stratonovich(cdw(2,:),cdt);
I02=stratonovich(cdt,cdw(2,:));
I20I02=I20+I02;
I111 = dw(1)^3/6-0.5*dt*dw(1);
I222 = dw(2)^3/6-0.5*dt*dw(2);
cI22=0.5*cdw(2,:).^2-0.5*cdt;
I221=stratonovich(cI22,cdw(1,:));
cI21I12=cdw(1,:).*cdw(2,:);
I212I122=stratonovich(cI21I12,cdw(2,:))-0.5*I10;
I221I212I122=I221+I212I122;
cI11=0.5*cdw(1,:).^2-0.5*cdt;
I112=stratonovich(cI11,cdw(2,:));
I211I121=stratonovich(cI21I12,cdw(1,:))-0.5*I20;
I112I211I121=I112+I211I121;
xts=zeros(2,1);
xte=zeros(2,1);
a=0.1; b1=0.2; b2=0.3;
xts(1)=x0s(1)*(1+a*dt+b1*dw(1)+b2*dw(2)+ ...
b1^2*I11+b1*b2*I21I12+b2^2*I22+...
a^2*I00+a*b1*I10I01+a*b2*I20I02+...
b1^3*I111+b1^2*b2*I112I211I121+...
b1*b2^2*I221I212I122+...
b2^3*I222);
xte(1)=x0e(1)*exp(a*dt+b1*dw(1)+b2*dw(2)-...
0.5*(b1^2+b2^2)*dt);
13
a=0.2; b1=-0.4; b2=0.7;
xts(2)=x0s(2)*(1+a*dt+b1*dw(1)+b2*dw(2)+...
b1^2*I11+b1*b2*I21I12+b2^2*I22+...
a^2*I00+a*b1*I10I01+a*b2*I20I02+...
b1^3*I111+b1^2*b2*I112I211I121+...
b1*b2^2*I221I212I122+...
b2^3*I222);
xte(2)=x0e(2)*exp(a*dt+b1*dw(1)+b2*dw(2)-...
0.5*(b1^2+b2^2)*dt);
end
function J = stratonovich(X,W)
%STRATONOVICH Summary
dW=diff(W,[],2);
n=size(dW,2);
J=sum((X(:,1:n)+X(:,2:n+1))/2.*dW,2);
end
Copyright: Jin Sun (2017)
Dated: January 26, 2017
14