SMM313 Numerical Methods: Applications

1. Random Number Generation

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

1/37

Ioannis Kyriakou 1. Random Number Generation 1 / 37

Topics

• General introduction to Monte Carlo simulation

• The inverse transform sampling method
• Sampling univariate normal random variates
• Simulating a probability distribution
• Sampling multivariate normal random variates
• Relevant Matlab applications
• Students’ self-studying. Glasserman (2004): Chapter 2 (Sections 2.2, 2.3)

2/37

Ioannis Kyriakou 1. Random Number Generation 2 / 37

Notation

• N (µ, σ 2 ): normal distribution with mean µ and variance σ 2

• CDF: Cumulative Distribution Function
d
• =: ‘equality in distribution’
• ∼: [random variable] ‘distributed as’

3/37

Ioannis Kyriakou 1. Random Number Generation 3 / 37

General considerations

• What is Monte Carlo (MC) simulation?

◦ A computational algorithm that relies on repeated random sampling to obtain
numerical results.
◦ The idea is to use randomness to solve problems that might be deterministic in
principle.
◦ It is often the method of choice when it is difficult or impossible to use other
approaches.
• At the core of MC simulation is a sequence of apparently random numbers used to
drive simulation.
◦ Example. Simulation of Brownian motion and stock price process (see later).
• MC simulation is like statistical sampling:
◦ We sample from a population (e.g., stock prices, option payoffs) according to a given
model (e.g., Brownian motion, geometric Brownian motion, etc.)
◦ Given the samples, we estimate the various quantities of interest, such as the
population mean, variance, etc.

4/37

Ioannis Kyriakou 1. Random Number Generation 4 / 37

Uniform distribution

Assume available sequence U (1) , U (2) , U (3) . . . of independent random variables from
Unif[0, 1], each satisfying

   0, u<0
(j)
P U ≤ u = u, 0≤u≤1. (1)

1, u>1


5/37

Ioannis Kyriakou 1. Random Number Generation 5 / 37

Uniform random numbers

• The Matlab command

rand
generates a uniform random variate from [0, 1].
• We can repeat the simulation M times and check if the simulator is good enough:
◦ Compare the theoretical mean (E (U) = 1/2 = 0.5) with the mean of the samples
◦ Compare the theoretical variance (Var (U) = 1/12 = 0.0833) with the variance of the
samples

6/37

Ioannis Kyriakou 1. Random Number Generation 6 / 37

Inverse transform sampling method I

• Most simulations entail sampling random variables from distributions other than
the uniform.
• A typical simulation uses methods for transforming samples from the uniform
distribution to samples from other distributions.
• Suppose we want to sample from a CDF F :
◦ We want to generate a random variable X with the property that

P(X ≤ x) = F (x)
for all x
• Inverse transform method sets

X = F −1 (U), U ∼ Unif[0, 1], (2)

where F −1 is the inverse of F .

7/37

Ioannis Kyriakou 1. Random Number Generation 7 / 37

Inverse transform sampling method II

• Assume a continuous and strictly increasing CDF F .

• Figure below illustrates an example of mapping u = F (x) : values 0 ≤ u < F (0)
(e.g., u2 ) are mapped to x < 0 (e.g., x2 ), whereas values F (0) ≤ u < 1 (e.g., u1 )
are mapped to x ≥ 0 (e.g., x1 ).

8/37

Ioannis Kyriakou 1. Random Number Generation 8 / 37

Inverse transform sampling method III

Assume X ∼ F .
To verify that the inverse transform (2) generates samples from F , check that:

P(X ≤ x) = P F −1 (U) ≤ x = P(U ≤ F (x)) = F (x).

• First equality follows from (2).

• Second equality follows from the fact that the events F −1 (u) ≤ x

and
{u ≤ F (x)} coincide for all u and x.
• Third equality follows from (1) as 0 ≤ F (x) ≤ 1.

9/37

Ioannis Kyriakou 1. Random Number Generation 9 / 37

Sampling using the inverse transform method: example I

• Exponential distribution. This is the distribution of the time T between consecutive

jumps of a Poisson process with rate 1/θ.
• Exponential distribution has mean θ and its CDF is

F (τ ) = 1 − e −τ /θ , τ ≥ 0. (3)

• Use U = F (T ) and invert (3) to get

U = 1 − e −T /θ ⇒ T = −θ ln(1 − U).

• This can also be implemented as

T = −θ ln(U),
d
as U = 1 − U ∼ Unif[0, 1].
• Note. Recall that the Matlab command rand generates a uniform random variate
from [0, 1].

10/37

Ioannis Kyriakou 1. Random Number Generation 10 / 37

Sampling using the inverse transform method: example II

• Arcsine law. The random time t = X at which the standard Brownian motion
(BM) Wt attains its maximum over the time interval [0, 1] has CDF

2 √
F (x) = arcsin( x), 0 ≤ x ≤ 1. (4)
π
• Use U = F (X ) and invert (4) to get
√
 
2 Uπ
U = arcsin( X ) ⇒ X = sin2 .
π 2

• Using the identity 2 sin2 (t) = 1 − cos(2t) for 0 ≤ t ≤ π/2 further yields
 
2 Uπ 1 1
X = sin = − cos(Uπ).
2 2 2

11/37

Ioannis Kyriakou 1. Random Number Generation 11 / 37

Sampling using the inverse transform method: example III

• Rayleigh distribution. If we condition a standard Brownian motion Wt starting at

the origin to be at b at time 1, i.e., conditional on W0 = 0 and W1 = b, then its
maximum Y over the time interval [0, 1] has a Rayleigh distribution with CDF

F (y ) = 1 − e −2y (y −b) , y ≥ b.

• Equation u = F (y ), where u ∈ (0, 1), gives

e −2y (y −b) = 1 − u ⇒ −2y (y − b) = ln(1 − u) ⇒

1
y 2 − by + ln(1 − u) = 0,
2
which solves at
1 1
q
y = b± b 2 − 2 ln(1 − u).
2 2

12/37

Ioannis Kyriakou 1. Random Number Generation 12 / 37

Sampling using the inverse transform method: example III cont’d I

• We require that y ≥ b, as the maximum of the Brownian path must be at least as

large as the terminal value b.
• For u ∈ (0, 1) ⇒ (1 − u) ∈ (0, 1), hence

ln(1 − u) < 0 ⇒ −2 ln(1 − u) > 0 ⇒

1s 2 1
b − 2 ln(1 − u) ≥ b
2 | {z } 2
≥b 2

and
1 1
q
y = b+ b 2 − 2 ln(1 − u) ≥ b (5)
2 |2 {z }
≥ 21 b

which is acceptable.
• We reject
1 1
q
y = b− b 2 − 2 ln(1 − u) < b.
2 2

13/37

Ioannis Kyriakou 1. Random Number Generation 13 / 37

Sampling using the inverse transform method: example III cont’d II

• Applying U = F (Y ) and from (5) we get

1 1
q
−1
Y =F (U) ⇒ Y = b + b 2 − 2 ln(1 − U)
2 2
or, alternatively,

1 1p 2
Y = b+ b − 2 ln U,
2 2
d
as U = 1 − U ∼ Unif[0, 1].
• Exercise. Generate samples from the Rayleigh distribution with parameter b = 2,
compute its mean and variance, and compare against the true mean
E (Y ) = 2.21068 and variance Var (Y ) = 0.03424.

14/37

Ioannis Kyriakou 1. Random Number Generation 14 / 37

Sampling using the inverse transform method: example IV

• Even if the inverse F −1 is not known explicitly, the inverse transform method is still
applicable through numerical evaluation of F −1 .
• Computing x = F −1 (u) is equivalent to finding the root x of the equation

F (x) = u ⇒ F (x) − u = 0.

• This is possible using the Matlab function fsolve:

◦ Fix generated value u from Unif[0, 1]
◦ Implementing fsolve(@(x)F (x) − u, x0 ) returns the corresponding value x from the
distribution F
• Exercise. Check with exponential distribution with mean θ = 3.
• Note. Problem can be dealt with even if F is not known explicitly (see Ballotta and
Kyriakou, 2014) – not examinable.

15/37

Ioannis Kyriakou 1. Random Number Generation 15 / 37

Sampling using the inverse transform method: example V

• Discrete distributions. Consider the discrete random variable N taking possible

values c0 < · · · < cn
◦ cn may be infinite, e.g., see Poisson distribution
• Probability mass function: pj := P (N = cj ) , j = 0, 1, 2, . . . , n
• CDF: Fi := P (N ≤ ci ) = ij=0 pj , i = 0, 1, 2, . . . , n
P

• CDF is non-decreasing step function: for ci ≤ c ′ < ci+1 ,

P{i:c ≤c ′ <ci+1 }
P (N ≤ c ′ ) = j=0i pj = Fi

16/37

Ioannis Kyriakou 1. Random Number Generation 16 / 37

Sampling using the inverse transform method: example V cont’d I

To sample from F of the discrete distribution:

• Generate a uniform U ∼ Unif[0, 1]
• Find K ∈ {1, . . . , n} such that FK −1 < U ≤ FK
• Set N = cK (cK corresponds to FK )

17/37

Ioannis Kyriakou 1. Random Number Generation 17 / 37

Sampling using the inverse transform method: example V cont’d II

Example of discrete distribution. Random variable N ∼ Poisson(θ) :

c0 = 0 < c1 = 1 < c2 = 2 < · · ·

e −θ θj
pj = P(N = j) = , j = 0, 1, 2, . . .
j!
X i
Fi = P(N ≤ i) = pj , i = 0, 1, 2, . . .
j=0
⌊k⌋
X
P(N ≤ k) = pj , for any k > 0
j=0

18/37

Ioannis Kyriakou 1. Random Number Generation 18 / 37

Sampling using the inverse transform method: example V cont’d III

• In the figure above (obtained with parameter θ = 3), if, for example, u = 0.35, then

N = 2 as F1 < u < F2

(see slide 17).

19/37

Ioannis Kyriakou 1. Random Number Generation 19 / 37

Sampling using the inverse transform method: example V cont’d IV

• In cases where the CDF has flat sections (see figure above), i.e., there are multiple
values along the step matching to a CDF value, we choose the smallest.
• In the figure above, if, for example, u = 0.2, then

N = inf {k : Fk ≥ u} = 1 as u = F1 .

• From slide 17, F0 < u ≤ F1 ⇒ N = 1.

20/37

Ioannis Kyriakou 1. Random Number Generation 20 / 37

Sampling using the inverse transform method: example V cont’d V

• Exercise. Sample from Poisson(θ = 3) in Matlab.

• Alternatively, in Matlab, random numbers from the Poisson distribution with mean
parameter θ can be generated using

poissrnd(theta).

• The extended command

poissrnd(theta,m,n)
generates an m-by-n array.

21/37

Ioannis Kyriakou 1. Random Number Generation 21 / 37

More about the inverse transform method, etc.

• The inverse transform method maps the input U to the output X ∼ F .

• This can be useful in the implementation of antithetic sampling and stratified
sampling (see variance reduction techniques later).
• The inverse transform method requires just one uniform random variable per
sample output, so other methods requiring more than one uniform for each sample
generated implicitly result in higher-dimensional representations (e.g.,
Acceptance-Rejection techniques –not examinable– see Glasserman, 2004).
• Alternative methods should be considered only when the use of the inverse
transform method is restricted.
• Note. In Matlab,
random(pd)
returns a random number from the distribution specified by the probability
distribution object pd. For a list of the acceptable values for pd, as well as the
parameters for that distribution, refer to Matlab’s documentation (https:
//[Link]/help/stats/[Link]).

22/37

Ioannis Kyriakou 1. Random Number Generation 22 / 37

Generating univariate normal: linear transformation property

• If Z ∼ N (0, 1) (standard normal distribution), then the random variable

X = µ + σZ

is normal with mean µ and variance σ 2 .

• It therefore suffices to consider methods for sampling from N (0, 1).
• Note. In Matlab, command randn generates a random variate from N (0, 1).

23/37

Ioannis Kyriakou 1. Random Number Generation 23 / 37

Building simulated distributions
M
1. Generate a collection of M random numbers X (j)

j=1
.
2. Fix the number of classes (bins). You may use the rule:
√
nbins = floor( M)

(floor rounds to nearest integer less than or equal to that element).

3. The Matlab command

histogram(x, nbins, ‘Normalization’, ------)

specifies the normalization scheme of the histogram values.

• The normalization scheme affects the scaling of the histogram along the vertical
axis. For example:
◦ ‘pdf’: probability density function estimate. The area, Ai , of each bar is the relative
number
P of observations in the bin, fi /M, and the sum of the bar areas is
all i fi /M = 1
◦ ‘cdf’: cumulative distribution function estimate. The height, Hi , of each bar is the
cumulative
P relative number of observationsPin each bin and all previous bins,
k≤i fk /M. The height of the last bar is all k fk /M = 1
24/37

Ioannis Kyriakou 1. Random Number Generation 24 / 37

Multivariate normal distribution I

• A d-dimensional normal distribution is characterized by a mean d-column vector µ

and a d × d covariance matrix Σ.
• To qualify as a covariance matrix, Σ must be:
◦ Symmetric (Σ and Σ⊤ must be equal), and
◦ Positive semidefinite (all eigenvalues of Σ must be nonnegative – as a symmetric
matrix, Σ automatically has real eigenvalues)
• The eigenvalues of Σ are precisely the solutions λ to the equation

det (Σ − λId ) = 0,
where ‘det’ is the determinant of the matrix formed by Σ − λId and Id the d × d
identity matrix.
• Note. In Matlab, command eig(Σ) returns the column vector of the eigenvalues of
Σ.

25/37

Ioannis Kyriakou 1. Random Number Generation 25 / 37

Multivariate normal distribution II

• If the d-random vector X has multivariate normal distribution, then the i th

component Xi has normal distribution with mean µi and variance σi2 = Σii ,
i = 1, . . . , d.
• Components Xi and Xj have covariance

Cov (Xi , Xj ) := E ((Xi − µi ) (Xj − µj )) = Σij .

• Correlation between the components Xi and Xj is given by

ρij := Σij / (σi σj ) ,

hence Σij = ρij σi σj and Σii = ρii σi σi = σi2 .

• Covariance matrix Σ :
σ12
 
ρ12 σ1 σ2 ··· ρ1d σ1 σd
 ρ21 σ2 σ1 σ22 ρ2d σ2 σd 
Σ= .
 
.. .. ..
 . . . 
ρd1 σd σ1 ρd2 σd σ2 · · · σd2

26/37

Ioannis Kyriakou 1. Random Number Generation 26 / 37

Multivariate normal distribution III

• If Σ is positive definite (all eigenvalues of Σ are positive), then the normal

distribution has a density (see Glasserman, 2004, p. 64).
• If the d × d symmetric matrix Σ is positive semidefinite but not positive definite,
then:
◦ the rank of Σ is less than d (Σ is rank-deficient),
◦ Σ fails to be invertible,
◦ there is no normal density with covariance matrix Σ.

27/37

Ioannis Kyriakou 1. Random Number Generation 27 / 37

Multivariate normal distribution IV

• Let the d-dimensional random vector Z := (Z1 , . . . , Zd ) be standard normal with

mean d-vector µ = 0 and covariance matrix Id .
• The components Z1 , . . . , Zd of Z are independent.
• We define the d-dimensional normal distribution characterized by some mean
d-vector µ and d × d covariance matrix Σ as the distribution of the random vector

X = µ + AZ

for any d × d matrix A that satisfies

AA⊤ = Σ.

• The resulting distribution is independent of which such A is chosen.

28/37

Ioannis Kyriakou 1. Random Number Generation 28 / 37

Generating multivariate normal

To sample X from the d-dimensional normal distribution with mean vector µ and
covariance matrix Σ:
• Generate independent standard normal variables Z1 , . . . , Zd and assemble them into
the d-random vector Z with mean vector equal to 0 and covariance matrix given by
Id
• Find a matrix A for which AA⊤ = Σ
• Apply the linear transformation X = µ + AZ

29/37

Ioannis Kyriakou 1. Random Number Generation 29 / 37

Generating multivariate normal using Cholesky factorization I

• With A being lower triangular, the representation

AA⊤ = Σ

is called a Cholesky factorization of Σ.

• For a positive definite d × d covariance matrix, Cholesky factorization yields
  
A11 0 ··· 0 A11 A21 · · · Ad1
 A21 A22 0   0 A22 Ad2 
Σ= . ..  .
  
. ..   .. ..
 .. .. .  . . . 
Ad1 Ad2 · · · Add 0 0 · · · Add
• If Σ is positive definite, it has a Cholesky factorization and the matrix A is unique.
• Note. In Matlab, command A =chol(Σ,‘lower’) produces the lower triangular
matrix A satisfying AA⊤ = Σ.
Matrix Σ must be positive definite, otherwise Matlab displays an error message.

30/37

Ioannis Kyriakou 1. Random Number Generation 30 / 37

Generating multivariate normal using Cholesky factorization II

Example. Consider the 2 × 2 covariance matrix

σ12
 
σ1 σ2 ρ
Σ= .
σ1 σ2 ρ σ22

• For σ1 , σ2 > 0, the Cholesky factor is

 
σ1 0
A= p .
σ2 ρ σ2 1 − ρ2

• Thus, we can sample X = (X1 , X2 ) by setting

X1 = µ1 + A11 Z1 = µ1 + σ1 Z1
p
X2 = µ2 + A21 Z1 + A22 Z2 = µ2 + σ2 ρZ1 + σ2 1 − ρ2 Z2

with Z1 , Z2 independent standard normals.

31/37

Ioannis Kyriakou 1. Random Number Generation 31 / 37

Generating multivariate normal using Cholesky factorization III

• In general, a lower triangular matrix A is particularly convenient because it reduces

the calculation of X = µ + AZ to:

X1 = µ1 + A11 Z1 ,
X2 = µ2 + A21 Z1 + A22 Z2 ,
..
.
Xd = µd + Ad1 Z1 + Ad2 Z2 + . . . + Add Zd .

• A full multiplication of the vector Z by the matrix A would require approximately

twice as many multiplications and additions.

32/37

Ioannis Kyriakou 1. Random Number Generation 32 / 37

Generating multivariate normal using Cholesky factorization IV

Exercise. Implement Cholesky factorization in Matlab.

33/37

Ioannis Kyriakou 1. Random Number Generation 33 / 37

Generating multivariate normal using Cholesky factorization V

• If Σ is merely positive semidefinite, then it is rank-deficient.

• It follows that any matrix A satisfying AA⊤ = Σ must also be rank-deficient.
• If A is lower triangular and rank-deficient, then at least one element of the diagonal
of A must be zero resulting in det(A) = 0.
• Having at least one element of the diagonal of A equal to zero results in failure of
the algorithm that attempts the Cholesky factorization (see Glasserman, 2004, p.
73–74).

34/37

Ioannis Kyriakou 1. Random Number Generation 34 / 37

Eigenvector factorization I
• Alternatively to Cholesky factorization, the equation AA⊤ = Σ can be solved by
diagonalizing Σ.
• As a symmetric d × d matrix, Σ has real eigenvalues λ1 , . . . , λd .
• Additionally, by positive definiteness or semidefiniteness of Σ, the λi are
nonnegative.
• It follows that
Σ = V ΛV ⊤ ,
where matrix V is orthogonal VV ⊤ = Id and Λ is the diagonal matrix with

diagonal entries λ1 , . . . , λd .
• Hence, if we choose
 √ 
λ1 √0 ··· 0
 0 λ2 0 
A = V Λ1/2 = V  , (6)
 
.. .. ..
. .
√.
 
0 0 ··· λd

then
AA⊤ = V ΛV ⊤ = Σ.
35/37

Ioannis Kyriakou 1. Random Number Generation 35 / 37

Eigenvector factorization II

• Note. In Matlab, using [V , Λ] =eig(Σ) produces the matrices V and Λ. We obtain

Λ1/2 using sqrt(Λ).
• Unlike the Cholesky factor, the matrix A in (6) has no particular structure that can
provide a computational advantage in evaluating AZ (see slide 30), nor is this
matrix faster to compute than the Cholesky factor.
• However, generating the normal random vector X using (6) is sometimes useful in
designing variance reduction techniques (the particular application is beyond the
scope of this syllabus; refer to Glasserman, 2004, p. 88–91).
• Exercise. Implement eigenvector factorization in Matlab.

36/37

Ioannis Kyriakou 1. Random Number Generation 36 / 37

References

Ballotta, L., I. Kyriakou. 2014. Monte Carlo simulation of the CGMY process and
option pricing. Journal of Futures Markets 34(12) 1095–1121.
Glasserman, P. 2004. Monte Carlo Methods in Financial Engineering. Springer,
New York.

37/37

Ioannis Kyriakou 1. Random Number Generation 37 / 37