FIM 548, Lecture #12

Quasi Monte Carlo

Chapter 5 in Glasserman and Chapter 15 of Owen

Lecture and Codes by Andrew Papanicolaou

Slides originally of Y. Kitapbayev

Feb 26, 2023

1 / 20
Integration using MC
I So far we have focused on computing
Z 1
θ= ψ(x)dx .
0

I Monte Carlo integration can be especially useful for

estimating high-dimensional integrals.
I Suppose we want to estimate
Z 1Z 1
θ= ψ(x, y )dxdy .
0 0

I Monte Carlo is easy because we just need to sample uniformly

on [0, 1) × [0, 1),
N
1 X
θ̂ = ψ(Xi , Yi ) .
N
i=1

2 / 20
Quasi Monte Carlo (QMC)

I We already saw that √

standard MC methods have a
convergence of O(1/ N).
I To accelerate it, one can apply so-called quasi Monte Carlo
method, or low-discrepancy method.
I These methods seek to increase accuracy by generating points
that are too evenly distributed to be random.
I Due to their deterministic nature, statistical methods do not
apply to QMC methods.
I In this lecture we focus on d-dimensional integrals,

θ = Eψ(U) ,

where U ∼ uniform[0, 1)d . Using inverse CDFs we can apply

to a great many other d-dimensional distributions.

3 / 20
Integration Lattices

I We are already familiar with basic QMC.

I The congruential method for pseudo random numbers.

Xi+1 = mod(aXi , m)
Ui+1 = Xi+1 /m

with m = 231 and a = 16807. Generates a scalar sequence of

pseudo-randoms on (0, 1).
I A similar idea, but perhaps more uniformly separated are
rank-1 lattices,
Xi = mod(iV1 , 1) ,
where V1 ∈ Rd . For
√ example, the Fibonacci lattice in d = 2
1+ 5
takes V1 = (1, 2 ).

4 / 20
Lattices

congruential grid rank 1

1 1 1

0.8 0.8 0.8

0.6 0.6 0.6

0.4 0.4 0.4

0.2 0.2 0.2

0 0 0
0 0.5 1 0 0.5 1 0 0.5 1

Hammersley Latin Hypercube Stratified Random

1 1 1

0.8 0.8 0.8

0.6 0.6 0.6

0.4 0.4 0.4

0.2 0.2 0.2

0 0 0
0 0.5 1 0 0.5 1 0 0.5 1 5 / 20
A Simple QMC Example

I Once we have our QMC points in [0, 1)d , we only need to

plug them into our random number generator in place of our
pseudo-random uniforms.
I For example, we can use an integration lattice in d = 2 to
estimate a function of a jointly-normal pair.
I For example,
N q 
1 X
E cos(kX k) ≈ cos Xi2 + Yi2 ,
N
i=1

where (Xi , Yi ) are given by Box-Muller with points from an

integration lattice.

6 / 20
 1 2
1
cos(kxk)e − 2 kxk dx
R
Example: θ = 2π R2

%% S i m p l e QMC Example
n = 2ˆ8;

p h i = (1+ s q r t ( 5 ) ) / 2 ; % F i b o n a c c i l a t t i c e
v1 = [ 1 / ( n+1) ; p h i ] ;
U qmc = mod ( v1 ∗ [ 1 : n ] , 1 ) ;

X = s q r t (−2∗ l o g ( U qmc ( 1 , : ) ) ) . ∗ c o s ( 2 ∗ p i ∗U qmc ( 2 , : ) ) ;

Y = s q r t (−2∗ l o g ( U qmc ( 1 , : ) ) ) . ∗ s i n ( 2 . ∗ p i ∗U qmc ( 2 , : ) ) ;

t h e t a q m c = mean ( c o s ( s q r t (X.ˆ2+Y . ˆ 2 ) ) ) ;

t h e t a m c = mean ( c o s ( s q r t ( sum ( r a n d n ( 2 , n ) . ˆ 2 , 1 ) ) ) ) ;

7 / 20
Discrepancy

I Let us assume that d is the dimension of a problem, and we

want to estimate the integral of ψ over [0, 1)d .
I Now we aim to fill the cube uniformly using some
deterministic rule.
I There are a couple discrepancy definitions.
I We define discrepancy for point set {x1 , . . . , xn } as
 
 #{xi ∈ [0, a)} 
D(x1 , . . . , x2 ) = sup  − [0, a) ,
 
a∈[0,1)d n

where [0, a) = [0, a1 ) × · · · × [0, an ), where #{xi ∈ [0, a)}


denotes the number of xi contained in [0, a), and [0, a)
 
denotes the measure or volume of the interval.

8 / 20
Low Discrepancy Sets
I For Ui ∼ uniform[0, 1)d it is known that
√
2nD(U1 , U2 , . . . , Un )
lim sup √ =1.
n→∞ log log n
I It is known that deterministic  low-discrepancy
 sets can have
log(n)d
D(u1 , u2 , . . . , un ) that is O n .
I Integration error when using a low-discrepancy set is
quantified with the Koksma-Hlawka inequality,

|θ̂ − θ| ≤ Vhk (ψ)D(u1 , u2 , . . . , un ) ,

where Vhk (ψ) is the Hardy-Krause total variation.

I Examples of sequences with low discrepancy include Halton
sequence, Sobol sequence, Faure sequence, and Niederreiter
sequence.
9 / 20
Digital Nets With Low-Discrepancy
Sobol Sobol Scramble
1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Halton Halton Scramble

1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

Figure: Digital nets don’t have the periodicity of lattices. Integration

lattices are good for integrating smooth functions with smooth periodic
features.
10 / 20
Randomized QMC

I Randomized digital nets offer a way to take advantage of

low-discrepancy sets and avoid some of the pitfalls in basic
QMC.
I The Matousek-Owen scramble is commonly used with the
Sobol set.
I Reverse-Radix scrambling is used for the Halton set.

11 / 20
Example: Scrambled Sobol Set to Generate Correlated
Brownian Motion
%% QMC Sim Brownian Motion
T = 3/12;
dt = 1/365;
N = r o u n d (T/ d t ) ;
Rho = −.7;
dW = z e r o s (N, 2 ) ;

S b l = s o b o l s e t ( 2 , ’ S k i p ’ , 1 e3 , ’ Leap ’ , 1 e2 ) ;
p = s c r a m b l e ( Sbl , ’ MatousekAffineOwen ’ ) ;
p = n e t ( p , N) ’ ;

U1 = p ( 1 : N) ;
U2 = p ( (N+1) : end ) ;
dW( : , 1 ) = s q r t (−2∗ l o g ( U1 ) ) . ∗ c o s ( 2 ∗ p i ∗U2 ) ∗ s q r t ( d t ) ;
dW( : , 2 ) = s q r t (−2∗ l o g ( U1 ) ) . ∗ s i n ( 2 ∗ p i ∗U2 ) ∗ s q r t ( d t ) ;

dW( : , 2 ) = Rho∗dW( : , 1 )+s q r t (1−Rho ˆ 2 ) ∗dW( : , 2 ) ;

12 / 20
Example: Heston Model Call Price
Heston Call Implied Volatility (T=3 months)
0.23
MC
QMC
RQMC
0.22 Quadrature

0.21

0.2

0.19

0.18

0.17

0.16
-0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25

Figure: Using QMC and RQMC to estimate Heston call price, and then
comparing implied vols.

13 / 20
Copulas
I Copulas are the CDFs of marginally uniform[0, 1) random
variables.
I For X ∈ Rd , with X ` denoting `th element.
I Let F` be the (marginal) CDF of X ` .
I We have X ` = F`−1 (U` ) where U` is a marginally uniform[0, 1)
random variable, and is the `th element of U ∈ [0, 1)d .
I The CDF of U is copula C .
I Gaussian copula:
s Z −1 Z Φ−1 (ud )
(2π)−d Φ (u1 ) 1 ∗ −1
C (u1 , . . . , ud ) = ··· e − 2 x ρ x dx ,
|ρ| −∞ −∞

where ρ ∈ Rd×d is a correlation matrix and |ρ| = det(ρ).

I Archimedean copula:
C (u1 , . . . , ud ) = ϕ−1 (ϕ(u1 ) + · · · + ϕ(ud )) ,
where ϕ is the copula the generator.
14 / 20
Example: Large Portfolio of Co-Dependent Insurance
Claims
I Let us consider an insurance example.
I Individual loss distributions have standard normal distribution
(i.e., marginally standardnormal).
I U ` = Φ(X ` ) are jointly dependent through a Clayton copula,
with
1
ϕ(t) = (t −ν − 1) ,
ν
where ν > 0 is the parameter.
I A claim is made if X ` ≥ 3.
I We let the portfolio have ` = 1, 2, . . . , 100. (i.e., 100 policies).
I We want to gain an understanding of the magnitude of the
total claims given that at least one claim is made.
I For example, think about fire insurance for homeowners in
California.
15 / 20
Marshall-Olkin Algorithm

I We can sample from the copula using the Marshall-Olkin

Algorithm
1. V ∼ LS −1ϕ
2. Ũ ∼ uniform[0, 1)d
3. Return U where U = ϕ(− log(Ũ)/V ).
R∞
I LS ϕ (v ) = 0 ϕ(t)e −vt dt is the Laplace-Stieltjes transform of
ϕ, which for Clayton copula has known inverse,

LS −1
ϕ (v ) ∝ v
1/ν−1 −v
e ,

i.e., V is gamma distributed with shape 1/ν and scale 1.

16 / 20
Example: Large Portfolio of Co-Dependent Insurance
Claims
I We can run this 100-dimensional example with Monte Carlo,
QMC and RQMC, and we’ll see generally comparable results.
I We can also run the example with a normal approximation to
the data, which requires some shrinkage of the covariance
matrice’s eigenvalues.
I The reduced error from QMC/RQMC does not pop out at you
in this example, or in most other examples.
I To see error reduction we will need to run many trials, save
the MC/QMC/RQMC estimators, and then compare statistics
of the aggregate distributions to see which method has the
least variation.
I We estimate θ = probability of at least 1 claim, M = the
expectation number of claims given at least 1, and the
expectation of the peaks-over-threshold distribution (POT)
given at least 1 claim.
17 / 20
Example: Large Portfolio of Co-Dependent Insurance
Claims
theta mc theta qmc
40 40

20 20

0 0
0.03 0.035 0.04 0.045 0.05 0.035 0.04 0.045

M mc M qmc
40 40

20 20

0 0
2 2.5 3 3.5 4 4.5 2.8 3 3.2 3.4 3.6 3.8

pot mc pot qmc

40 40

20 20

0 0
8 10 12 14 16 9 10 11 12

Figure: 500 trials of MC and RQMC estimations for the insurance

example using Clayton copula.

18 / 20
Example: Large Portfolio of Co-Dependent Insurance
Claims

Based on the 500 trials from on the previous slide, we have the
following statistics:

E [θmc ] = 0.0409, sd(θmc ) = 0.0032

E [θqmc ] = 0.0409, sd(θqmc ) = 0.0017
E [Mmc ] = 3.3259, sd(Mmc ) = 0.2896
E [Mqmc ] = 3.3202, sd(Mqmc ) = 0.1458
E [potmc ] = 10.9189, sd(potmc ) = 0.9520
E [potqmc ] = 10.9020, sd(potqmc ) = 0.4778

As we can see from this table, error in RQMC estimation is about

half of the standard deviation in MC estimation.

19 / 20
Summary

I QMC is good for some problems in multiple dimensions.

I Implementation is more complicated than standard MC.
I Issues such as selection of a lattice/net, parameters for it, and
the absence of statistical analysis of error (i.e., no CLT) make
it less straightforward.
I In 1 and 2 dimension there are quadrature rules that are
better (e.g., Gauss-Lobatto quadrature).
I The Koksma-Hlawka bound and possible O(log(n)d /n) error
is appealing, but is only theoretical and in reality it requires
some now-how to see reduced error.

20 / 20