The Seven-League Scheme: Deep learning forlarge time step Monte Carlo simulation of SDEs
Shuaiqiang Liu3 , Lech A. Grzelak2 , Kees Oosterlee1
1

Utrecht University, the Netherlands
Rabobank, Utrecht, the Netherlands
3
Centrum Wiskunde & Informatica, the Netherlands
2

March, 2021

C.W. Oosterlee, Utrecht University0

1 / 34

Outline1

Basics of Stochastic Differential Equations

2

The Seven-League Scheme
Stochastic Collocation Monte Carlo sampler
Neural Networks
Algorithm details
Error analysis
The 7L-CDC variant

3

Numerical results

4

Summary and Outlook

C.W. Oosterlee, Utrecht University0

2 / 34

Introduction

Stochastic Differential Equations (SDEs) are widely used to model
processes, in finance and other fields.
Develop an accurate numerical scheme for SDES (in strong
convergence sense) to carry out large time step simulation.
Employ a neural network to obtain this numerical scheme andsolve the SDEs, efficiently.

C.W. Oosterlee, Utrecht University

1 Basics of Stochastic Differential Equations

3 / 34

Stochastic differential equations

Random variable Y (t), on probability space (Ω, Σ, P) with filtration
F(t)t∈[0,T ] , sample space Ω, σ-algebra Σ and prob. measure P.
Consider the generic scalar Itô SDE,
dY (t) = a(t, Y (t), θ)dt + b(t, Y (t), θ)dW (t),0 ≤ t ≤ T,with drift a(t, Y (t), θ), diffusion b(t, Y (t), θ), parameters θ, Wiener
process W (t), and initial value Y0 := Y (t = 0).

C.W. Oosterlee, Utrecht University

1 Basics of Stochastic Differential Equations

4 / 34

The strong convergence
Solution in integral form
Z t+∆tZ t+∆tY (t + ∆t) = Y (t) +a(s, Y (s), θ)ds +b(s, Y (s), θ)dW (s)ttDefinition
Let the exact solution of an SDE at time ti be given by Y (ti ), its
discrete approximation Ỹ (ti ) with time step ∆t converges in the strong
sense, with order βs ∈ R+ , if
E|Y (ti ) − Ỹ (ti )| ≤ K (∆t)βs .

C.W. Oosterlee, Utrecht University

1 Basics of Stochastic Differential Equations

5 / 34

Numerical discretization
Classical numerical schemes
Euler-Maruyama (strong order βs = 0.5):√
Ŷi+1 = Ŷi + a(ti , Ŷi , θ)∆t + b(ti , Ŷi , θ) ∆t X̂i+1
Milstein (strong order βs = 1.0):
Ŷi+1 = Ŷi√+ a(ti , Ŷi , θ)∆t + b(ti , Ŷi , θ) ∆t X̂i+1
1 0
2+
b (ti , Ŷi , θ)b(ti , Ŷi , θ)∆t(X̂i+1− 1).
2

where Ŷi+1 := Ŷ (ti+1 ) is a realization from Ỹ (ti+1 ), and X̂i+1 is drawn
from X ∼ N (0, 1), b0 (·) the first derivative of b(·).

C.W. Oosterlee, Utrecht University

1 Basics of Stochastic Differential Equations

6 / 34

Higher-order numerical approximation

Common ways to improve numerical accuracy.
Include higher-order terms1 : ODE Runge-Kutta schemes plus
high-order random terms of Itô calculus. For example, there are 8
terms for βs = 1.5, and 12 terms for βs = 2.0.
√
1d
Y (t + ∆t) = Y (t) + a∆t + b ∆tX + b0 b∆t(X 2 − 1) + ...
|{z 2}
12 terms

Reduce time step ∆t, e.g., use a finer time grid.

1

Eckhard Platen (1999). An introduction to numerical methods for stochastic differential equations. Acta Numerica.

C.W. Oosterlee, Utrecht University

1 Basics of Stochastic Differential Equations

7 / 34

Data-driven numerical schemes
The Euler-Maruyama scheme,√Ŷi+1 = Ŷi + a(ti , Ŷi , θ)∆t + b(ti , Ŷi , θ) ∆t X̂i+1 ,
{z
} |
{z
}
|
α0

α1

(
α0 = Ŷi + a(ti , Ŷi , θ)∆t,√
α1 = b(ti , Ŷi , θ) ∆t,
The Milstein scheme,1 0α0 = Ŷi + a(ti , Ŷ√i , θ)∆t + 2 b (ti , Ŷi , θ)b(ti , Ŷi , θ),
α1 = b(ti , Ŷi , θ) ∆t,
α2 = 12 b0 (ti , Ŷi , θ)b(ti , Ŷi , θ).
A numerical scheme is a mapping function.
dỸ (ti+1 )|Ỹ (ti ) =

m−1
X



αj X j , with αj = Hj ti , Ỹ (ti ), ∆t, θ .

j=0C.W. Oosterlee, Utrecht University

1 Basics of Stochastic Differential Equations

8 / 34

The Seven-League Scheme“We march through time with the Seven-League scheme(the 7L scheme).”The 7L scheme consists of two key components,
Stochastic Collocation Monte Carlo sampler (SCMC)1 ,efficient sampling from an ’expensive’ distribution;
Neural network, as a powerful function approximator.
1

L. A. Grzelak et al. (2019). The stochastic collocation Monte Carlo sampler. Quant. Finance.

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

9 / 34

Stochastic Collocation Monte Carlo sampler
Two scalar random variables, Y and X , are connected by,dd

FY (Y ) = U = FX (X ),where U ∼ U([0, 1]), FY (ȳ ) := P(Y ≤ ȳ ), FX (x̄) := P(X ≤ x̄).
When FY (ȳ ) and FX (x̄) are strictly monotonic, we haveȳ = FY−1 (FX (x̄)) := g(x̄),
both distributional and also element-wise.
Choosing optimal collocation points, (x̂j ,ŷj ), to approximate,ȳ = g(x̄) ≈ ĝm (x̄) =mXŷj `j (x̄),

j=1where ŷj = FY−1 (FX (x̂j )) (e.g. x̂j Gauss-Hermite quadraturepoints), `j (·) interpolation basis functions.
C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

10 / 34

SCMC algorithm

Exponential convergence to the target distribution.

SCMC algorithm
Calculate CDF FX (xj ) on (x1 , . . . , xm ), e.g., fromGauss-Hermite quadrature, giving m pairs (xj , FX (xj ));1

Invert the target CDF yj = FY−1 (FX (xj )), j = 1, . . . , m, andform m pairs of collocation points (xj , yj );2

Define the interpolation function, y = gm (x), based on thesem pairs (xj , yj );3Obtain sample Ŷ by applying the mapping Ŷ = gm (X̂ ), where
sample X̂ is drawn from X .
4

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

11 / 34

Conditional distribution by SCMC
Given the current state Y (t), the conditional variable Y (t + ∆t) can be
written as,d
Y (t + ∆t)|Y (t) = g(X ) ≈ gm (X ).
Generating a sample path is drawing samples from the probability
distribution conditional on a previous realization,jŶi+1 |Ŷi = gm (X̂i+1).

In the case of Lagrange interpolation, we arrive at the coefficient form,jŶi+1 |Ŷi = gm (X̂i+1
)=

m−1
X

j
α̂i+1,j X̂i+1,

j=0

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

12 / 34

Learn stochastic collocation points

We use an artificial neural network (ANN) to determine theconditional collocation points, yj (ti+1 ), based on the previous
realization,yj (ti+1 )|Ŷi = Hj (Ŷi , ti , ∆t, θ).

When the conditional stochastic collocation points are known, themapping gm (·) is constructed via interpolation.

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

13 / 34

Artificial Neural Network (ANN)
Artificial neural network is used as a function approximator.
Fully connected neural network is a composite function,Ĥ(x́|Θ̃) = h(LA ) (. . . h(2) (h(1) (x́; θ̃1 ); θ̃2 ); . . . θ̃LA ),with x́ input variables, Θ̃ weights and biases, LA hidden layers.
Minimizing the loss function to approximate the target function,
Xargmin L(Θ̃) := argmin
D(Ĥ(x́i |Θ̃), ýi ),Θ̃Θ̃where D(·, ·) measures the distance between predicted Ĥ(x́i |θ)and true ýi values.

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

14 / 34

The offline stage
The ANN learns conditional stochastic collocation points,yj (τ + dτ )|Ŷ (τ ) = Hj Ŷ (τ ), τ, dτ, θ , j = 1, 2, . . . , m,

1

Generate many data samples for different time steps dτ andmodel parameters θ.

2

Use a numerical scheme (e.g. Euler) with tiny ∆τ to find theconditional distribution Y (τ + dτ )|Y (τ ).

3Calculate stochastic collocation points yj (τ + dτ )|Ŷ (τ ) with
SCMC.

4

Train ANN (supervised learning) to obtain Ĥj ≈ Hj .

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

15 / 34

The online stage

Use the trained ANNs to generate the sample paths.
1

Partition [0, T ], 0 < t1 < t2 < . . . < tN ≤ T , with ∆t = ti+1 − ti .

2

Given a sample Ŷi at time ti , compute m collocation points at time
ti+1 using the trained ANN,ŷj (ti+1 )|Ŷi = Ĥj (Ŷi , ti , ti+1 − ti , θ), j = 1, 2, . . . , m.

3Compute the mapping function gm (·) using SCMC.

4

Sample from X to obtain a sample Ŷi+1 |Ŷi = gm (X̂i+1 ).

5

Return to Step 2 by ti+1 →− ti , iterate until terminal time T .

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

16 / 34

45

SC

40

int
po

0.15
Density

35

0.10

30
Y(t)

0.05
0.00

Path No.2

30
20
Y(t)
10

25

Conditional PDF

Generating sample paths with 7L

Path No.2

20
15
10

Path No.1

0

ti

−1

ti

ti

time

+1

5

Path No.10(a) Sample paths by 7L

ti

−1

time

ti

ti

+1

(b) The 2D projection

Figure: Here conditional SC points, represented by , are conditional on aprevious realization, denoted by F.

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

17 / 34

Error analysis
The error from SCMC:vuu m!√π Ψ(2m) (ξˆ1 )p.
E[|g(x) − gm (x)|] ≤ |m | = t
2m(2m)!
The error from ANNs:Z
E [|gm (x) − g̃m (x)|] =
R

Z
≤

m
X

Aj `j (x) fX (x)dx

j=1

max{|A1 |, . . . , |Am |}fX (x)dx

R= max{|A1 |, . . . , |Am |}.
E [|g(x) − g̃m (x)|] ≤ E [|g(x) − gm (x)|] + E [|gm (x) − g̃m (x)|]p|m | + max{|A1 |, . . . , |Am |}.
≤C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

18 / 34

Strong convergence error of 7L

7L strong convergence
Assuming the approximation errors from ANN and SCMC are
negligible, the strong convergence of the 7L scheme is defined, as
E|Ỹ (ti ) − Y (ti )| ≤ (∆τ )  K (∆t)βs ,where ∆τ is used for offline ANN training, and time step ∆t for online
ANN prediction, with ∆τ  ∆t.
For example, the Euler scheme would need a much finer time step, bya factor κ = ∆t/∆τ to achieve a similar accuracy in the strong sense.

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

19 / 34

7L-CDCThe 7L Compression-Decompression (CDC) variant can further
reduce the computational cost.

45
40

30

0.2

25

Y(t)

density

35
0.25

0.15
1

0.1
0.05

20
15
10

0.5

0
25

20

15

10

Y(t)

5

0

0

5
0

t

(a) Paths for marginal SC pointsti − 1timetiti + 1(b) Sample paths by 7L-CDC

Figure: Schematic diagram of the 7L-CDC scheme at time ti . Left: Marginal
SC points. Right: Sample paths generated by 7L-CDC.
C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

20 / 34

7L-CDC computation time
The time ratio between the 7L-CDC and 7L schemes isγ=tI M + tA Ms
tI
Ms
=
+
,
tA M
tA
M

with tA the computational time of the ANN, tI for the interpolation.
Given the fact that the number of sample paths is typically much larger
than the number of SC points M  Ms ,
γ≈

tI
.
tA

When the employed interpolation is computationally cheaper than theANN, γ < 1, so that the 7L-CDC scheme needs fewer computations
than the 7L scheme.

C.W. Oosterlee, Utrecht University

2 The Seven-League Scheme

21 / 34

4. Numerical experiments

With a(t, Y (t)) = µY (t) and b(t, Y (t)) = σY (t), we have Geometric
Brownian Motion (GBM),
0 ≤ t ≤ T,

dY (t) = µY (t)dt + σY (t)dW (t),where drift µ and volatility σ are constant, with initial value Y0 . The
model parameters are θ := {µ, σ}, and the analytic solution is,dC.W. Oosterlee, Utrecht University

1

2

√

Y (t) = Y0 e(µ− 2 σ )(t−t0 )+σ t−t0 X .

(1)

3 Numerical results

22 / 34

ANN settings and off-line training
We created the data set, using Euler, ∆τ = 0.01. There are around 80,000
data samples.
ANN
input

Ĥ1 (·) output
Ĥ2 (·) output
Ĥ3 (·) output
Ĥ4 (·) output
Ĥ5 (·) output

Parameters
drift, µvolatility, σvalue, Y0time, τmax
point, ŷ1
point, ŷ2point, ŷ3point, ŷ4point, ŷ5

Value range(0.0, 0.10][0.05, 0.60][0.10, 15.0](0.0, 1.60](0.0, 25.65)(0.0, 25.98)(0.0, 27.84)(0.0, 54.67)(0.0, 154.35)

ANN hyper-parameters
Hidden layers
Neurons (per layer)
Activation
Initialization
Optimizer
Batch size
Initial Learning rate
C.W. Oosterlee, Utrecht University

Method
LHS
LHS
LHS
Equidistant
SCMC
SCMC
SCMC
SCMC
SCMC

value
4
50
Softplus
Glorot_uniform
Adam
1024
1e-3

3 Numerical results

23 / 34

Performance of the ANNs
The data set is divided into three parts, 80% for training, 10% forvalidation, 10% for testing. After 1500 training epochs, the ANN
converged.
P
Performance on the test dataset. MAE = M1 j |yj − ŷj |.
SC points
R2
MAE

ŷ1
0.999891
0.026

ŷ2
0.999947
0.027

ŷ3
0.999980
0.021

ŷ4
0.999892
0.071

ŷ50.999963
0.066

Predicted vs true collocation points, for example,(a) ŷ2
C.W. Oosterlee, Utrecht University(b) ŷ4
3 Numerical results

24 / 34

Comparing path-wise error
Monte Carlo path-wise error:Monte Carlo paths

Y(t)

4

0.200
0.175
0.150
0.125
0.100
0.075
0.050
0.025
0.000

Milstein strong conv.

7L-CDC-PCHIP strong conv.
7L-CDC-Lagrange strong conv.
7L-CDC-Chebyshev strong conv.

error

5

7L-CDC-Lagrange
7L-CDC-PCHIP
Milstein
Exact-GBM

32
1
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0t(c) Sample paths, ∆t = 0.5

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Δt(d) Strong convergence

Here σ = 0.3, r = 0.1, S0 = 1.0, T = 4.0. We use the barycentric version of
Lagrange interpolation.
C.W. Oosterlee, Utrecht University

3 Numerical results

25 / 34

CPU time

The CPU time (seconds) to reach the same accuracy (CPU: E3-1240,
3.40GHz): simulating 10,000 sample paths until terminal time T = 4.0, basedon 5 × 5 SC points.
Method /Time (Sec.)
7L-CDC Barycentric
7L-CDC Chebyshev
7L-CDC PCHIP
7L scheme
MilsteinCreate C
0.054
0.054
0.054
-

∆t = 1.0
Decom. C
4.93
9.78
11.39
-

Total
4.98
9.83
11.44
12.80
27.01

Create C
0.027
0.027
0.027
-

∆t = 2.0
Decom. C
2.48
4.93
5.73
-

Total
2.51
4.96
5.76
6.39
27.70

⇒ This is much faster on GPUs (recent results)!
We have essentially the same accuracy results with the OU process.

C.W. Oosterlee, Utrecht University

3 Numerical results

26 / 34

Applications in financeWe can carry out large time step Monte Carlo simulations for
Pricing path-dependent options
Computing the sensitivities of path-dependent options.

C.W. Oosterlee, Utrecht University

3 Numerical results

27 / 34

Bermudan option, big time steps

The option holder has the right to exercise the contract at
pre-specified monitoring dates up to maturity time T .
Relative error: rel =VAref (t0 )−VA (t0 );VAref (t0 )Longstaff-Schwartz method; Y0 = 1.0, r = 0.1, strike priceK̃ = 1.1, T = 4.0, number of monitoring dates Nb , M = 100, 000.

σ=0.30Value
MethodAnalytic MC
Milstein MC
7L-CDC

C.W. Oosterlee, Utrecht University∆t = 1.0, Nb =4∆t = 0.5, Nb =80.15213858 (0.00%)
0.13872771 (8.81%)
0.15234901 (0.14%)

0.16161876 (0.00%)
0.15429369 (4.53%)
0.16196264 (0.21%)

3 Numerical results

28 / 34

Pricing Asian options
Pay-off of a fixed strike Asian option:
VA (T ) = max (A(T ) − K̃ , 0),with A(T ) = N1b∈ [0, T ].

PNb

k =1 Ŷ (tk ) discrete average over Nb dates

The option value VA(t) = e−rT EQ
V (A(T )) F(0) .

Y0 = 1.0, K̃ = Y0 , r = 0.1, T = ∆t × Nb , M = 100, 000.

σ=0.40Value
MethodAnalytic MC
Milstein MC
7L-CDC

C.W. Oosterlee, Utrecht University∆t = 1.0, Nb =4∆t = 0.5, Nb =80.28515109 (0.00%)
0.26394277 (7.44%)
0.28482371 (0.11%)

0.25723594 (0.00%)
0.24717425 (3.91%)
0.25647592 (0.30%)

3 Numerical results

29 / 34

Computing Greeks
The Asian option’s vega is defined as


Nb
X
∂V
(T
,
Y
(t
);
σ)
∂Y
(t
)
∂V
i
i
= e−rT EQ 
F(0) .
∂σ∂Y (ti )∂σ
i=1

For M sample paths and N time points,




Nb
M X
m−1
X
X
∂ Ĥj
∂V
1 1
1
≈ e−rT
pj (X̂k ,i+1 ) F(0) .
A(T
)>
K̃
∂σ
MN
∂σ
k =1 i=1

j=0

The derivative of the ANNs Ĥj (·) w.r.t σ can be computed,∂ Ĥj (Ŷi , σ)∂ Ĥj (σ; Ŷi ) ∂ Ĥj (Ŷi ; σ) ∂ Ŷi
=
+
.
∂σ
∂σ
∂σ
∂ Ŷi
C.W. Oosterlee, Utrecht University

3 Numerical results

30 / 34

Asian option’s vega

0.7

7L Scheme
Exact-GBM

10
8

0.5

6

vega

∂Y(t)
∂σ

Exact-vega
Exact-GBM
7L Scheme

0.6

4

0.4

2

0.3

0
−2
0.0

0.5

1.0

1.5

2.0
t

2.5

3.0

3.5

4.0

(e) Path-wise sensitivity of Y (t)

0.2

10.0

12.5

15.0

17.5 20.0 22.5 25.0
Number of paths (x103)

27.5

30.0

(f) Asian option’s vega,∆t = 1.0

The parameters are Y0 = 1.0, r = 0.05, K̃ = Y0 , σ = 0.3, T = 4.0.

C.W. Oosterlee, Utrecht University

3 Numerical results

31 / 34

5. Summary and Outlook

Conclusion:
The 7L-scheme provides a data-driven numerical solver for SDEs.
Numerical error, in the sense of strong convergence, does not
grow significantly when the step increases.
Various applications in computational finance.
Outlook:
Parallel computing on GPUs (is already almost finalized!)
Solving multi-dimensional SDEs.

C.W. Oosterlee, Utrecht University

4 Summary and Outlook

32 / 34

References
1

S. Liu, et al (2020). The Seven-League Scheme: Deep learning for large time
step Monte Carlo simulations of stochastic differential equations,
arXiv:2009.03202.

2

G. Cybenko, et al(1989). Approximations by superpositions of sigmoidal
functions, Mathematics of Control, Signals and Systems.

3

R. Cameron, et al (1947). The Orthogonal Development of Non-Linear
Functionals in Series of Fourier-Hermite Functionals. Annals of Mathematics.

4

Eckhard Platen (1999). An introduction to numerical methods for stochastic
differential equations. Acta Numerica.

5

L. A. Grzelak, et al (2019). The stochastic collocation Monte Carlo sampler.
Quantitative Finance.

6

Y. Bar-Sinai, et al (2019). Learning data-driven discretizations for partial
differential equations. PNAS.

7

Dmitry Yarotsky (2017). Error bounds for approximations with deep ReLU
networks. Neural Networks.

C.W. Oosterlee, Utrecht University

4 Summary and Outlook

33 / 34

Training ANNs
The layer function:!
(`)
zj = ϕ(`)X(`) (`−1)(`)
wij zi+ bj,

iStochastic Gradient Descent (SGD) optimizes the parameters,(w ← w − η(i) ∂L ,∂wb ← b − η(i) ∂L , i = 0, 1, 2, ...,∂bwhere η is a learning rate.

C.W. Oosterlee, Utrecht University

4 Summary and Outlook

34 / 34