SMM313 Numerical Methods: Applications

3. Variance Reduction

Module Leader
Prof Ioannis Kyriakou
[Link]@[Link]
Bayes Business School
City, University of London

1/50

Ioannis Kyriakou 3. Variance Reduction 1 / 50

Topics

• Introduction and implementation of control variate technique: Application to

arithmetic Asian options
• Other variance reduction techniques:
◦ Antithetic sampling
◦ Stratified sampling
• Relevant Matlab applications
• Students’ self-studying. Glasserman (2004): Chapter 4 (Sections 4.1, 4.2,
4.3)

2/50

Ioannis Kyriakou 3. Variance Reduction 2 / 50

Notation

• N µ, σ 2 indicates a normal distribution with mean µ and variance σ 2 .

• (I)CDF stands for (Inverse) Cumulative Distribution Function.

• IID stands for Independent and Identically Distributed (random variables).

3/50

Ioannis Kyriakou 3. Variance Reduction 3 / 50

Variance reduction techniques

• We have seen previously that the standard error of a price estimate depends on the
number of simulations M:
◦ Higher M results in lower standard error
◦ However, increasing M results in higher computational time
◦ We would like to achieve further variance reduction (reduced standard error) without
significant increase of the computational cost.
• In this module, we consider three popular variance reduction techniques:
◦ Control variate sampling
◦ Antithetic variate sampling
◦ Stratified sampling.

4/50

Ioannis Kyriakou 3. Variance Reduction 4 / 50

Traditional (“crude”) MC estimation
M
• Let C (j) j=1 be IID outputs of M simulations. In our context, C (j) represents the


discounted payoff of a derivative security in the j th simulation.

• The objective is to estimate E C (j) = E (C ) = µC (identical distribution: drop

superscript “j”), which in Monte Carlo we “traditionally” estimate using the sample
mean
M
1 X (j)
µ̂C (M) = C .
M
j=1

• µ̂C (M) is an unbiased estimator of µC as

M
1 X  (j)  1
E (µ̂C (M)) = E C = ME (C ) = µC .
M M
j=1

• Also,
M
1 X 
(j)
 1 1
Var (µ̂C (M)) = Var C = 2 M Var(C ) = σC2 , (1)
M2 M M
j=1

where σC2 = Var C (j)

= Var(C ) (identicalP
distribution: drop superscript “j”;
independence: simply interchange Var and ).
5/50

Ioannis Kyriakou 3. Variance Reduction 5 / 50

Control variates I

• Control variates are among

◦ the most effective,
◦ easy-to-implement,
◦ hence, popular techniques for improving the efficiency of Monte Carlo simulation.
• Mechanism: the method exploits the sampling uncertainty in an estimate of a
known quantity to reduce the sampling uncertainty in an estimate of an unknown
but related quantity of interest.
• As we will see: the higher the correlation, the stronger the variance reduction effect.

6/50

Ioannis Kyriakou 3. Variance Reduction 6 / 50

Control variates II

• Suppose that in each simulation j, we calculate some additional output P (j) along
with C (j) .
• Key features of control variates:

M
◦ (Across j) the pairs C (j) , P (j) j=1 are IID (identical distribution: drop superscript
“j”)
◦ (For given j) C (j) and

P (j) are not identically distributed but correlated
◦ Expectation E P (j) = E (P) =: µP is known.
• For some fixed constant b, construct:
 
(j)
Cb := C (j) − b P (j) − µP .

7/50

Ioannis Kyriakou 3. Variance Reduction 7 / 50

Control variate estimator I

• Calculate the sample mean:

M
1 X (j)
µ̂Cb (M) := Cb
M
j=1
M 
1 X  
= C (j) − b P (j) − µP .
M
j=1

• Call µ̂Cb (M) the control variate estimator: we say that the term µ̂P (M) − µP
serves as control in estimating µC .
• The control variate estimator is still unbiased as
M
1 X  (j)  1
E (µ̂Cb (M)) = E Cb = ME (Cb )
M M
j=1

= E (C ) − b (E (P) − µP ) = µC − b (µP − µP )
= µC .

8/50

Ioannis Kyriakou 3. Variance Reduction 8 / 50

Control variate estimator II

• The variance of the control variate estimator is

M
1 X 
(j)
 1
Var (µ̂Cb (M)) = 2
Var C b = 2 M Var (Cb )
M M
j=1
1
Var(C ) + b 2 Var(P) − 2b Cov(C , P)


=
M
1
σC2 + b 2 σP2 − 2bσC σP ρC ,P


=
M
1 2
= σ , (2)
M Cb
where
σC2 b := σC2 + b 2 σP2 − 2bσC σP ρC ,P ,
σC2 := Var(C ),
σP2 := Var(P),
σC σP ρC ,P := Cov(C , P) ̸= 0.

9/50

Ioannis Kyriakou 3. Variance Reduction 9 / 50

Control variate estimator III

• In practice, σC2 b is not known, hence it is estimated by its sample counterpart:

M
1 X  (j) 2
σ̂C2 b (M) := Cb − µ̂Cb (M) .
M −1
j=1

1 2
• Given the estimates µ̂Cb (M) and M σ̂Cb (M), 100α% confidence intervals can also be
constructed (in virtue of the Central Limit Theorem for large M):
 
−1
 α  σ̂Cb (M)
µ̂Cb (M) ± Φ 0.5 + √ ,
2 M

where Φ−1 denotes the inverse standard normal CDF.

10/50

Ioannis Kyriakou 3. Variance Reduction 10 / 50

Optimizing the control variate coefficient I

• The optimal coefficient b = b ∗ that minimizes variance (2) is given by

σC σC σP ρC ,P Cov(C , P)
b∗ = ρC ,P = = .
σP σP2 σP2

• It follows from solving

dσC2 b !
= 2bσP2 − 2σC σP ρC ,P = 0.
db
• We get from (1) and (2) with b = b ∗ that

σ 2 + b ∗2 σP2 − 2b ∗ σC σP ρC ,P


Var µ̂Cb∗ (M)
= C = 1 − ρ2C ,P . (3)
Var (µ̂C (M)) σC2

• From (3), the effectiveness of a control variate is determined by the strength of

correlation between the quantities C and P:
◦ The closer |ρC ,P | is to 1 ( ρC ,P can be positive or negative), the higher the variance
reduction.
11/50

Ioannis Kyriakou 3. Variance Reduction 11 / 50

Optimizing the control variate coefficient II

• In practice, the optimal coefficient b ∗ is unknown and needs to be estimated by its

sample counterpart:
PM
C (j) − µ̂C (M) P (j) − µ̂P (M)



j=1
b̂ ∗ := . (4)
(j) − µ̂ (M) 2
PM

j=1 P P

• Expression (4) is the slope of the least-squares regression line through the points
 (j) (j)
M
C ,P j=1
.
• Note. b̂ ∗ has to be computed first separately based on an independent set of M
simulations to the ones used to calculate the control variate estimate, otherwise an
amount of bias is introduced.

12/50

Ioannis Kyriakou 3. Variance Reduction 12 / 50

Example of use of control variates: case of Asian option

• The arithmetic Asian option with discounted payoff

n
!+
−rT 1X
Cn := e Si − K
n
i=1

does not admit a true closed form pricing formula: no true explicit formula for
E (Cn ) is known.
• This is because the distribution law of the arithmetic average is not known: a sum
of lognormal random variables is not lognormal!
• Hence, a possible way of computing E (Cn ) is by Monte Carlo simulation as
explained previously.

13/50

Ioannis Kyriakou 3. Variance Reduction 13 / 50

Example of control variate: geometric Asian option

• On the contrary, the geometric Asian option with discounted payoff

 !1/n +
Yn
Pn := e −rT  Si − K
i=1

admits a true closed form pricing formula: E (Pn ) is known explicitly.

• As we show in the Appendix (see slide 47), a (Black–Scholes style) pricing formula
follows from the fact that the geometric average of lognormal random variables is
itself lognormal – non-examinable background details of the formula.

14/50

Ioannis Kyriakou 3. Variance Reduction 14 / 50

Modified Black–Scholes formula for geometric Asian option

• An option on the geometric average can thus be evaluated using the modified
Black–Scholes formula:
h i
µP = E (Pn ) = e −rT S0 e (r −q)T̄ Φ (d1 ) − K Φ (d2 ) ,

where
σ̄ 2
     
1 S0
d1 := √ ln + r −q+ T̄ ,
σ̄ T̄ p K 2
d2 := d1 − σ̄ T̄ ,
Φ is the standard normal CDF, and modified parameters T̄ , σ̄, q are applied where
indicated above.
• In addition, arithmetic and geometric Asian payoff samples are highly correlated,
hence the geometric Asian option can serve as a control variate in the simulation of
the arithmetic Asian option.

15/50

Ioannis Kyriakou 3. Variance Reduction 15 / 50

Implementation of control variate technique

• Exercise 1. Modify the Matlab script file from Exercise 3 in the slides of the previous
topic to include a control variate in the shape of a geometric Asian call option.

16/50

Ioannis Kyriakou 3. Variance Reduction 16 / 50

Implementation of control variate technique: Asian option I

Figure: Asian call option. Comparisons of distributions of different estimators: reducing sampling
uncertainty from “crude” estimator to control variate estimator.

17/50

Ioannis Kyriakou 3. Variance Reduction 17 / 50

Implementation of control variate technique: Asian option II

Figure: Arithmetic vs. geometric option payoff samples. Fitted simple linear regression model:
C = 0.0444 + 1.0522P with R 2 = ρ2C ,P = 0.9992 = 99.8% (strong linear relationship).

18/50

Ioannis Kyriakou 3. Variance Reduction 18 / 50

Antithetic variates I

• The method attempts to reduce variance by introducing negative correlation in

pairs of simulations.
M M
• Let U (j) j=1 be uniformly distributed in [0, 1], then 1 − U (j) j=1 are also
 

uniformly distributed in [0, 1].

• The variables U (j) and 1 − U (j) form an antithetic pair in the sense that a large
value of one is accompanied by a small value of the other, and vice-versa.
• Mechanism: antithetic sampling suggests that an unusually large or small output
computed based on the first set of variables may be balanced by the value
computed from the second (“antithetic”) set, resulting in a reduction in variance.

19/50

Ioannis Kyriakou 3. Variance Reduction 19 / 50

Antithetic variates II

• The technique extends to other distributions through the inverse transform method:
random variables F −1 (U) and F −1 (1 − U) have same distribution F but are
antithetic as F −1 is monotone.
• If F is the standard normal distribution, then F −1 (1 − U) and F −1 (U) have the
same magnitudes but opposite signs (by symmetry of the distribution). See figure:

20/50

Ioannis Kyriakou 3. Variance Reduction 20 / 50

Antithetic variates III

• In Matlab, the plot in the previous slide was produced using

u = 0 : 0.01 : 1;
plot(u, norminv (1-u), u, norminv(u));
• The above translates to plotting F −1 (1 − u) and F −1 (u) as functions of u, where
F −1 corresponds to the inverse standard normal CDF:
◦ Either obtain F −1 yourselves as shown in the first lecture
◦ Or, use the Matlab function icdf. Many different F −1 are available – check Matlab’s
documentation.
For example, for the exponential distribution with mean θ = 0.5 use:
u = 0 : 0.01 : 1;
plot(u, icdf(‘exp’, 1-u, 0.5), u, icdf(‘exp’, u, 0.5));

21/50

Ioannis Kyriakou 3. Variance Reduction 21 / 50

Antithetic variates IV

• Therefore, in a simulation driven by IID standard normal random variables

 (j) M
Z j=1
, antithetic variates can be implemented by pairing with the sequence of
M
IID standard normal random variables −Z (j) j=1 .


22/50

Ioannis Kyriakou 3. Variance Reduction 22 / 50

Antithetic variates V

M
• Let C (j) j=1 be IID outputs of M simulations. In our context, C (j) represents the


discounted payoff of a derivative security in the j th simulation.

• The objective is to estimate E C (j) = E (C ) = µC (identical distribution: drop

superscript “j”).
n oM
• Additionally, let C̃ (j) be IID outputs obtained using an implementation of
j=1
antithetic sampling.
n oM
• In other words, we use as inputs to generate C (j) , C̃ (j) the antithetic pairs:
j=1
 (j) (j)

M
◦ Z , −Z in the standard normal case (see slide 29)
j=1


M
F −1 U , F −1 1 − U (j)
 (j)


◦ j=1
in the general distribution case.

23/50

Ioannis Kyriakou 3. Variance Reduction 23 / 50

Antithetic variates VI

Key features of antithetic variates:

n oM
• (Across j) the pairs C (j) , C̃ (j) are IID (identical distribution: drop
j=1
superscript “j”).
• (For given j) C (j) and C̃ (j) are negatively correlated, though they have the same
distribution:  


◦ E C (j) = E C̃ (j) =: µC

 
◦ Var C (j) = Var C̃ (j) =: σC2
 
◦ Cov C (j) , C̃ (j) = Cov(C , C̃ ) < 0.

24/50

Ioannis Kyriakou 3. Variance Reduction 24 / 50

Antithetic variate estimator I
• The antithetic variate estimator is given by the sample mean:
 
M M M
1 X (j) X (j)  1 X C (j) + C̃ (j)
µ̂CAV (M) := C + C̃ =
2M M 2
j=1 j=1 j=1

• The antithetic variate estimator is unbiased as

 
M 
1  X 
E (µ̂CAV (M)) = E C (j) + C̃ (j) 
2M
j=1
M
1 X  (j) 
= E C + C̃ (j)
2M
j=1
M
1 X   (j)   
= E C + E C̃ (j)
2M
j=1
M
1 X
= (E (C ) + E (C̃ ))
2M
j=1
M
1 X 1
= (µC + µC ) = M2µC = µC .
2M 2M 25/50

Ioannis Kyriakou j=1 Reduction

3. Variance 25 / 50
Antithetic variate estimator II

• The variance of the antithetic variate estimator is

 
M 
1 X 
Var (µ̂CAV (M)) = Var  C (j) + C̃ (j) 
4M 2
j=1
M
1 X  
= Var C (j) + C̃ (j)
4M 2
j=1
M
1 X
= Var(C + C̃ )
4M 2
j=1
1
= M(Var(C ) + Var(C̃ ) + 2 Cov(C , C̃ ))
4M 2
1  2 
= σC + σC2 + 2 Cov(C , C̃ )
4M !
1 σC2 + Cov(C , C̃ )
=
M 2
1
= σ 2 . (5)
M CAV 26/50

Ioannis Kyriakou 3. Variance Reduction 26 / 50

Antithetic variate estimator III
• We get from (1) based on 2M simulations and (5) that
 
1 2 + Cov(C , C̃ )
Var (µ̂CAV (M))
1 2
σC 2M σC Cov(C , C̃ )
= M1 AV = 1
=1+ .
Var (µ̂C (2M)) 2
2M σC
2
2M σC
σC2

• Note. To make the comparison efficient, we use 2M samples of C to obtain the

“crude” Monte Carlo estimate, as the antithetic one by construction is based on M
samples of C and M samples of C̃ .
• We say that antithetics reduce variance if

Var (µ̂CAV (M))

< 1⇒
Var (µ̂C (2M))
Cov(C , C̃ )
1+ < 1⇒
σC2
Cov(C , C̃ ) < 0. (6)

• (6) confirms the importance of negative correlation between samples of C and C̃ in

achieving variance reduction.
27/50

Ioannis Kyriakou 3. Variance Reduction 27 / 50

Antithetic variate estimator IV

• In practice, σC2 AV is not known, hence it is estimated by its sample counterpart:

M
!2
1 X C (j) + C̃ (j)
σ̂C2 AV (M) := − µ̂CAV (M) .
M −1 2
j=1

1 2
• Given the estimates µ̂CAV (M) and M σ̂CAV (M), 100α% confidence intervals can also
be constructed (in virtue of the Central Limit Theorem for large M):
 
−1
 α  σ̂CAV (M)
µ̂CAV (M) ± Φ 0.5 + √ ,
2 M

where Φ−1 denotes the inverse standard normal CDF.

28/50

Ioannis Kyriakou 3. Variance Reduction 28 / 50

Antithetic variates: example of application

• Condition (6) requires that negative correlation in the inputs Z (j) and −Z (j) used
respectively in the simulation of the stock prices in the GBM model
√
  
(j) 1 2 (j)
Sn = S0 exp r − σ T + σ TZ ,
2
√ 
  
(j) 1 2 (j)
S̃n = S0 exp r − σ T + σ T −Z ,
2

produces negative correlation between the outputs of paired simulations

 +
(j) (j)
Cn = e −rT Sn − K ,
 +
(j) (j)
C̃n = e −rT S̃n − K .

• This is possible by monotonicity of the mapping from inputs to outputs.

• Therefore, antithetic sampling can be used for path-independent options (see
Exercise 2).
• The method can be extended to path-dependent options (see Exercise 2).
29/50

Ioannis Kyriakou 3. Variance Reduction 29 / 50

Implementation of antithetic sampling technique

Exercise 2. Write a Matlab script file to price options using antithetic variates:
• Assume a plain vanilla call option with K = S0 = 100, T = 1, σ = 0.3 and r = 0.04.
• Assume an arithmetic average discretely monitored Asian call option with
K = S0 = 100, T = 1, n = 12, σ = 0.3 and r = 0.04.

30/50

Ioannis Kyriakou 3. Variance Reduction 30 / 50

Implementation of antithetic sampling technique: plain vanilla option
I

Figure: Plain vanilla call option. Comparisons of distributions of different estimators: reducing
sampling uncertainty from “crude” estimator to antithetic variate estimator.

31/50

Ioannis Kyriakou 3. Variance Reduction 31 / 50

Implementation of antithetic sampling technique: plain vanilla option
II

Figure: Plain vanilla call option. Comparisons of correlations: (Z , −Z̃ ), (S, S̃), (C , C̃ ). Reducing
correlation (absolute terms): ρZ ,−Z = −1.0 > ρS,S̃ = −0.9 > ρC ,C̃ = −0.4.

32/50

Ioannis Kyriakou 3. Variance Reduction 32 / 50

Asian option: antithetic variates vs. control variates I

Figure: Asian call option. Comparisons of distributions of different estimators: reducing sampling
uncertainty from “crude” estimator to antithetic variate estimator and control variate estimator.

33/50

Ioannis Kyriakou 3. Variance Reduction 33 / 50

Asian option: antithetic variates vs. control variates II

Figure: Asian call option. Comparisons of correlations: antithetic (left) vs. control variate
samples (right). Strong linear relationship in case of control variates: ρC ,P = 0.999. Weaker
relationship in case of antithetic variates: ρC ,C̃ = −0.4 < 0.

34/50

Ioannis Kyriakou 3. Variance Reduction 34 / 50

Variance reduction by stratified sampling

• From the Latin stratum ≡ layer (or subset) (strata ≡ layers/subsets).

• Stratified sampling refers broadly to a sampling mechanism that constrains the
fraction of observations drawn from specific subsets (strata) of the sample space.
• The general idea is to split the domain of the MC variable (say, the stock price)
into subintervals, and perform MC simulation for each subinterval separately.
• With stratification:
◦ We obtain a more representative sample
◦ The variance of the payoff in sub-samples is smaller resulting in smaller total variance.

35/50

Ioannis Kyriakou 3. Variance Reduction 35 / 50

Stratifying uniforms I
• Stratified sampling can be used to stratify random variables from the uniform
distribution Unif [0, 1].
• Partition the unit interval [0, 1] into the following N strata:
     
1 1 2 N −1
A1 = 0, , A2 = , , . . . , AN = ,1 .
N N N N
1
• Each of these intervals has probability N under the uniform distribution:
P i−1 i 1


N < U ≤ N = N i = 1, . . . , N.
,
• For each i = 1, . . . , N, let U (i,1) , . . . , U (i,M̃) be M̃ independent Unif [0, 1] random
variables (total number of simulations is M = N × M̃).
• For given i = 1, . . . , N, let V (i,1) , . . . , V (i,M̃) satisfy
 
(i,j) i −1 i i −1
V = + − U (i,j) , j = 1, . . . , M̃.
N N N

• Then, the M̃ random variables V (i,1) , . . . , V (i,M̃) constitute a stratified sample from
the uniform distribution that falls within the i th stratum, i.e.,
V (i,j) ∼ Unif i−1 i
 
N ,N . 36/50

Ioannis Kyriakou 3. Variance Reduction 36 / 50

Stratifying uniforms II

M̃
• Each of i = 1, . . . , N collections V (i,j)

j=1
is uniformly distributed between
(i − 1)/N and i/N:
 (i,j) M̃  M̃  M̃
◦ V j=1
has the distribution of U (i,j) j=1
conditional on U (i,j) j=1
∈ Ai .
• We may generalize into strata of unequal lengths:
◦ If Ai = (ai−1 , ai ], then the distribution of U (i,·) conditional on U (i,·) ∈ Ai is uniform
between ai−1 and ai
◦ We can sample from this conditional distribution by setting
V (i,·) = ai−1 + (ai − ai−1 ) U (i,·) .

37/50

Ioannis Kyriakou 3. Variance Reduction 37 / 50

Stratifying non-uniforms: example

Example. Consider the problem of generating stratified random numbers from the
standard normal distribution:
• Subdivide the space [0, 1] into, say, N = 10 equally spaced strata and generate, say,
M̃ = 1 uniform number forced to be in each stratum: V (1,1) , V (2,1) , . . . , V (10,1)
• Invert the standard normal
CDF to obtain
10 stratified standard normal samples (1
−1 (1,1) −1 (2,1) −1 (10,1)


per stratum): Φ V ,Φ V ,...,Φ V .
Note. In Matlab, Φ−1 (·) is computed using norminv (·, 0, 1).
See figure in the following slide for the case of the illustrative example N = 10, M̃ = 1:
of course, method applies for any N and M̃.

38/50

Ioannis Kyriakou 3. Variance Reduction 38 / 50

Stratifying non-uniforms: illustration

39/50

Ioannis Kyriakou 3. Variance Reduction 39 / 50

Terminal stratification of Brownian motion I

• Much of the variability in the option’s payoff can potentially be eliminated by

stratifying its terminal value.
• As a step in this direction, we detail the stratification of Brownian motion along its
terminal value.
• In the special case of a stock described by geometric Brownian motion with
constant volatility, we can also then stratify the terminal stock price.
• Assume a finite time period [0, T ].
• We wish to generate via stratified sampling the value of the Brownian motion Wt
at time t = T , that is, WT .

40/50

Ioannis Kyriakou 3. Variance Reduction 40 / 50

Terminal stratification of Brownian motion II

• Consider the case of N strata and M̃ simulations allocated to each stratum.

• For given stratum i, generate M̃ independent random variables U (i,1) , . . . , U (i,M̃)
from Unif[0, 1] and set
 
(i,j) i −1 i i −1
V = + − U (i,j) , j = 1, . . . , M̃.
N N N
 
• Then Φ−1 V (i,1) , . . . , Φ−1 V (i,M̃) form a stratified sample from the normal

distribution N (0, 1).

√ √  
• Then T Φ−1 V (i,1) , . . . , T Φ−1 V (i,M̃) form a stratified sample from the

normal distribution N (0, T ), i.e., the distribution of WT .

• Repeat for all strata i = 1, . . . , N.

41/50

Ioannis Kyriakou 3. Variance Reduction 41 / 50

Terminal stratification of Brownian motion: illustration

Example cont’d for Brownian motion. Assume N = 10 strata and produce M̃ = 1

(1,1) (2,1) (10,1)
sample per stratum: WT , WT , . . . , WT

42/50

Ioannis Kyriakou 3. Variance Reduction 42 / 50

Terminal stratification of stock prices and payoff values

• For an underlying stock S driven by the GBM

  
1
ST = S0 exp r − σ 2 T + σWT ,
2

then by driving the simulation of S with the stratified terminal Brownian motions
(i,j) (i,j)
WT , we also stratify the terminal stock prices ST .
• For terminal discounted path-independent payoff C , we subsequently stratify the
(i,j)
terminal payoff values CT .
• This is a consequence of the fact that ST and CT are monotone transformations of
WT .
• Note. If required to generate a discrete path S1 , . . . , Sn , for the purpose of pricing a
path-dependent option, starting from stratified terminal Sn one has to use the
Brownian bridge (see Glasserman, 2004, Section 4.2.3) – non-examinable.

43/50

Ioannis Kyriakou 3. Variance Reduction 43 / 50

Terminal stratification of stock prices and payoff values: illustration

Example cont’d for plain vanilla call. Assume N = 10 strata and produce M̃ = 1 sample
(1,1) (2,1) (10,1) (1,1) (2,1) (10,1)
per stratum: WT , WT , . . . , WT ⇒ ST , ST , . . . , ST ⇒
(1,1) (2,1) (10,1)
CT , CT , . . . , CT

44/50

Ioannis Kyriakou 3. Variance Reduction 44 / 50

Stratification in derivatives pricing: price estimator
• The stratified Monte Carlo price estimate is given by

N X M̃
1 X
µ̂Cstrat (N × M̃) := C (i,j) .
N × M̃ i=1 j=1

• The standard error of the stratified price estimate is computed as

q P
1 N
N i=1 σ̂Cstrat
2 (i, M̃)
σ̂Cstrat (N × M̃) := p , (7)
N × M̃

where σ̂C2 strat (i, M̃) is the sample variance of the M̃ samples of C (i,·) in the i th
stratum.
• (7) represents the average standard error across all strata.
• The 100α% confidence interval is given by:
h  α i
µ̂Cstrat (N × M̃) ± Φ−1 0.5 + σ̂Cstrat (N × M̃) ,
2
where Φ−1 denotes the inverse standard normal CDF. 45/50

Ioannis Kyriakou 3. Variance Reduction 45 / 50

Implementation of stratified sampling

Exercise 3. Write a Matlab script file to price a plain vanilla call option with
K = S0 = 100, T = 1, σ = 0.3 and r = 0.04 using stratified sampling.

46/50

Ioannis Kyriakou 3. Variance Reduction 46 / 50

Appendix. Geometric average distribution I

• In what follows, we show that the geometric average of lognormal random variables
is itself lognormal.
• We have that
n
!1/n
Y
Si = (S1 S2 · · · Sn )1/n
i=1
n n
! !
σ2
 
1X σX
= S0 exp r− ti + Wi .
2 n n
i=1 i=1

47/50

Ioannis Kyriakou 3. Variance Reduction 47 / 50

Appendix. Geometric average distribution II

• Define:
n n
1X σ2 X 1 2
σ̄ 2 = σ − σ̄ 2 .


T̄ = ti , (2i − 1)tn+1−i , q=
n
i=1
n2 T̄ i=1 2

• From Glasserman (2004, p. 100):

n n
!
X X
Wi ∼ N 0, (2i − 1)tn+1−i
i=1 i=1

then:
◦ v
n u n
σ X d
p u σ2 X p
Wi = T̄ t 2 (2i − 1)tn+1−i Z = σ̄ T̄ Z
n n T̄
i=1 i=1

where Z ∼ N (0, 1).

48/50

Ioannis Kyriakou 3. Variance Reduction 48 / 50

Appendix. Geometric average distribution III

◦
n
!1/n n
! n
!
σ2
 
Y 1X σX
Si = S0 exp r− ti + Wi
2n n
i=1 i=1 i=1
σ̄ 2
  p 
d
= S0 exp r −q− T̄ + σ̄ T̄ Z .
2

49/50

Ioannis Kyriakou 3. Variance Reduction 49 / 50

References

Glasserman, P. 2004. Monte Carlo Methods in Financial Engineering. Springer,

New York.

50/50

Ioannis Kyriakou 3. Variance Reduction 50 / 50