Simulation and Monte Carlo
With applications in finance and MCMC
J. S. Dagpunar
School of Mathematics
University of Edinburgh, UK
b 21 2006 il /S i ffi
Simulation and Monte Carlo
b 21 2006 il /S ii ffi
b 21 2006 il /S iii ffi
Simulation and Monte Carlo
With applications in finance and MCMC
J. S. Dagpunar
School of Mathematics
University of Edinburgh, UK
b 21 2006 il /S i ffi
Copyright © 2007 John Wiley & Sons Ltd, The Atrium, Southern Gate, Chichester,
West Sussex PO19 8SQ, England
Telephone +44 1243 779777
Email (for orders and customer service enquiries): [email protected]
Visit our Home Page on www.wiley.com
All Rights Reserved. No part of this publication may be reproduced, stored in a retrieval system or transmitted
in any form or by any means, electronic, mechanical, photocopying, recording, scanning or otherwise, except
under the terms of the Copyright, Designs and Patents Act 1988 or under the terms of a licence issued by the
Copyright Licensing Agency Ltd, 90 Tottenham Court Road, London W1T 4LP, UK, without the permission
in writing of the Publisher. Requests to the Publisher should be addressed to the Permissions Department,
John Wiley & Sons Ltd, The Atrium, Southern Gate, Chichester, West Sussex PO19 8SQ, England, or
emailed to [email protected], or faxed to +44 1243 770620.
Designations used by companies to distinguish their products are often claimed as trademarks. All brand names
and product names used in this book are trade names, service marks, trademarks or registered trademarks of
their respective owners. The Publisher is not associated with any product or vendor mentioned in this book.
1.2 Monte Carlo is integral estimation 4
1.3 An example 5
1.4 A simulation using Maple 7
1.5 Problems 13
2 Uniform random numbers 17
2.1 Linear congruential generators 18
2.1.1 Mixed linear congruential generators 18
2.1.2 Multiplicative linear congruential generators 22
2.2 Theoretical tests for random numbers 25
2.2.1 Problems of increasing dimension 26
2.3 Shuffled generator 28
2.4 Empirical tests 29
2.4.1 Frequency test 29
2.4.2 Serial test 30
2.4.3 Other empirical tests 30
2.5 Combinations of generators 31
2.6 The seed(s) in a random number generator 32
2.7 Problems 32
3 General methods for generating random variates 37
3.1 Inversion of the cumulative distribution function 37
3.2 Envelope rejection 40
3.3 Ratio of uniforms method 44
3.4 Adaptive rejection sampling 48
3.5 Problems 52
4 Generation of variates from standard distributions 59
4.1 Standard normal distribution 59
4.1.1 Box–Müller method 59
4.1.2 An improved envelope rejection method 61
4.2 Lognormal distribution 62
viii Contents
6.3.5 Discrete hedging 116
6.4 Asian options 118
6.4.1 Naive simulation 118
6.4.2 Importance and stratified version 119
6.5 Basket options 123
6.6 Stochastic volatility 126
6.7 Problems 130
7 Discrete event simulation 135
7.1 Poisson process 136
7.2 Time-dependent Poisson process 140
7.3 Poisson processes in the plane 141
7.4 Markov chains 142
7.4.1 Discrete-time Markov chains 142
7.4.2 Continuous-time Markov chains 143
Contents ix
7.5 Regenerative analysis 144
7.6 Simulating a G/G/1 queueing system using the three-phase method 146
7.7 Simulating a hospital ward 149
7.8 Problems 151
8 Markov chain Monte Carlo 157
8.1 Bayesian statistics 157
8.2 Markov chains and the Metropolis–Hastings (MH) algorithm 159
8.3 Reliability inference using an independence sampler 163
8.4 Single component Metropolis–Hastings and Gibbs sampling 165
8.4.1 Estimating multiple failure rates 167
8.4.2 Capture–recapture 171
8.4.3 Minimal repair 172
8.5 Other aspects of Gibbs sampling 176
8.5.1 Slice sampling 176
8.5.2 Completions 178
all mathematics students are conversant with it from year 1. I believe this is true of many
other mathematics departments. The disadvantage of slow numerical processing in Maple
is neutralized by the wide range of probabilistic, statistical, plotting, and list processing
functions available. A large number of specially written Maple procedures are available
on the website accompanying this book (www.wiley.com/go/dagpunar_simulation). They
are also listed in the Appendices.
1
The content of the book falls broadly into two halves, with Chapters 1 to 5 mostly
covering the theory and probabilistic aspects, while Chapters 6 to 8 cover three application
areas. Chapter 1 gives a brief overview of the breadth of simulation. All problems at the
end of this chapter involve the writing of Maple procedures, and full solutions are given
in Appendix 1. Chapter 2 concerns the generation and assessment of pseudo-random
numbers. Chapter 3 discusses three main approaches to the sampling (generation) of
random variates from distributions. These are: inversion of the distribution function, the
envelope rejection method, and the ratio of uniforms method. It is recognized that many
other methods are available, but these three seem to be the most frequently used, and
they have the advantage of leading to easily programmed algorithms. Readers interested
in the many other methods are directed to the excellent book by Devroye (1986) or an
earlier book of mine (Dagpunar, 1988a). Two short Maple procedures in Appendix 3
allow readers to quickly ascertain the efficiency of rejection type algorithms. Chapter 4
deals with the generation of variates from standard distributions. The emphasis is on
short, easily implemented algorithms. Where such an algorithm appears to be faster
than the corresponding one in the Maple statistics package, I have given a listing in
Appendix 4. Taken together, I hope that Chapters 3 and 4 enable readers to understand
how the generators available in various packages work and how to write algorithms for
distributions that either do not appear in such packages or appear to be slow in execution.
Chapter 5 introduces variance reduction methods. Without these, many simulations are
incapable of giving precise estimates within a reasonable amount of processing time.
Again, the emphasis is on an empirical approach and readers can use the procedures in
1
of simulation.
Chapter 8 deals with the other burgeoning area of simulation, namely Markov chain
Monte Carlo and its use in Bayesian statistics. Here, I have been influenced by the
works of Robert and Casella (2004) and Gilks et al. (1996). I have also included several
examples from the reliability area since the repair and maintenance of systems is another
area that interests me. Maple has been quite adequate for the examples discussed in this
chapter. For larger hierarchical systems a purpose-built package such as BUGS is the
answer.
There are problems at the end of each chapter and solutions are given to selected
ones. A few harder problems have been designated accordingly. In the text and problems,
numerical answers are frequently given to more significant figures than the data would
warrant. This is done so that independent calculations may be compared with the ones
appearing here.
I am indebted to Professor Alastair Gillespie, head of the School of Mathematics,
Edinburgh University, for granting me sabbatical leave for the first semester of the
2005–2006 session. I should also like to acknowledge the several cohorts of simulation
students that provided an incentive to write this book. Finally, my thanks to Angie for
her encouragement and support, and for her forbearance when I was not there.
Glossary
beta
1
A random variable that is beta distributed with p.d.f. f
x
=
Standard Brownian motion.
c.d.f. Cumulative distribution function.
C
The transpose of a matrix C.
Cov
f
X Y
The covariance of X and Y where f
x y
is the joint p.d.f./p.m.f.
of X and Y (the subscript is often dropped).
e.s.e. Estimated standard error.
Exp
A r.v. that has the p.d.f. f
x
= e
−x
x≥ 0, where >0.
E
.
gamma
A gamma distributed r.v. with the p.d.f. f
x
=
x
−1
e
−x
/
x≥ 0, where >0>0.
gig
A r.v. distributed as a generalized inverse Gaussian distribution with
the p.d.f. f
0 1 2, where 0 <p<1.
N
2
A normal r.v. with expectation and variance
2
, or the density
itself.
N
A vector r.v. X distributed as multivariate normal with mean and
covariance matrix .
p.d.f. Probability density function.
p.m.f. Probability mass function.
Poisson
A r.v. distributed as a Poisson with the p.m.f. f
x
=
x
e
−x
0 1
A continuous r.v. that is uniformly distributed in the interval (0, 1).
Var
f
X
The variance of a random variable X that has the p.d.f. or p.m.f f
(the subscript is often dropped).
v.r.r. Variance reduction ratio.
Weibull
A Weibull distributed random variable with the p.d.f. f
x
=
x
−1
e
−
x
2
n
An r.v. distributed as chi-squared with n degrees of freedom, n =
1 2. Therefore,
2
n
= gamma
n/2 1/2
.
1
P
Equals 1 if P is true, otherwise equals 0.
∼ ‘Is distributed as’. For example, X ∼ Poisson
indicates that X
has a Poisson distribution.
= y In Maple or in pseudo-code this means ‘becomes equal to’. The
value of the expression to the right of = is assigned to the variable
or parameter to the left of =.
1
Introduction to simulation and
Monte Carlo
A simulation is an experiment, usually conducted on a computer, involving the use of
random numbers. A random number stream is a sequence of statistically independent
random variables uniformly distributed in the interval [0,1). Examples of situations where
simulation has proved useful include:
I
=
0
x
−1
e
−x
dx (1.1)
for a specific value of >0. Consider a random variable X having the probability density
function (p.d.f.) f on support 0 where
f
x
= e
−x
Then from the definition of the expectation of a function of a random variable
Equation (1.1) leads to
I
= E
f
X
−1
and assuming the
X
i
are independent, the variance
of
I
is given by
Var
f
I
=
Var
f
X
−1
n
Thus, the standard deviation of the sampling distribution of the statistic
I
0
x
09
e
−x
dx
Evaluating a definite integral 3
Firstly, we need to know how to generate values from the probability density function
f
x
= e
−x
. It will be seen in Chapter 3 that this is done by setting
X
i
=−ln R
i
(1.3)
where
R
i
i= 0 1
is a random number stream with R
i
I
19
when n =4. This is unknown
as
f
X
09
is unknown. However, the sample standard deviation of X
09
is s where
s =
1
4 −1
⎛
⎝
4
i=1
x
09
i
6
. We learn from this that an uncritical design
and analysis of a simulation will often lead to a vast consumption of computer time.
Is it possible to design the experiment in a more efficient way? The answer is ‘yes’.
Rewrite the integral as
I
19
=
0
xe
−x
1
x
01
dx (1.5)
There is a convenient method of generating variates
x
from the probability density
function
g
x
= xe
where X has the density of Equation (1.6). Given the same four random numbers, two
random values (variates) can be generated from Equation (1.6). They are x
1
=−ln R
1
R
2
=
49206 and x
2
=−ln R
3
R
4
= 72928. Therefore,
I
19
=
1
2
49206
−01
+72928
−01
=08363 (1.7)
4 Introduction to simulation and Monte Carlo
This is a great improvement on the previous estimate. A theoretical analysis shows that
g
<<
f
X
09
.
Now try to estimate I
199
using both methods with the same random numbers as before.
It is found that when averaging 1/X
001
sampled from g
I
199
= 098226, which is very
close to the true value, I
199
=099581. This is not the case when averaging X
099
sampled
from f.
This simple change in the details of the sampling experiment is an example of a
variance reduction technique. In this case it is known as importance sampling, which is
explored in Chapter 5.
1.2 Monte Carlo is integral estimation
How is it that the Monte Carlo method evolved from its rather specialized use in integral
evaluation to its current status as a modelling aid is used to understand the behaviour of
1
w
i
n
i= 1m
. However, it should not be too difficult to
generate a realization of the waiting time W
1
of the first customer. Given this, the known
structure of the system is then used to generate a realization of W
2
given a knowledge
of the state of the system at all times up to the departure of customer 1. Note that it is
often much easier to generate a realization of W
2
(given W
1
and the state of the system
up to the departure of customer 1) than it is to write down the conditional distribution
of W
2
given W
1
. This is because the value assumed by W
2
can be obtained by breaking
down the evolution of the system between the departures of customers 1 and 2 into easy
stages. In this way it is possible to obtain realizations of values of W
. Here is another
example. Suppose we wish to estimate the expectation of the time average of queue length
in
0T
. Monte Carlo is used to estimate E
f
1/T
T
0
Q
t
dt
, where
Q
t
t≥ 0
is a stochastic process giving the queue length and f is the probability density for the
paths
. It can be seen that
P
X<x
=
−
ft
1 −1
x>t
dt
where f is the probability density of the length of the shortest path and 1
x>t
= 1if
x>t; else 1
x>t
= 0. Since the probability can be expressed as an integral and since
realizations of 1
x>t
can be simulated by breaking down the calculation into easier parts
using the structure of the network, the probability can be estimated using Monte Carlo.
In fact, if ‘minimum’ is replaced by ‘maximum’ we have a familiar problem in project
planning. This is the determination of the probability that the duration of a project does
not exceed some specified value, when the individual activity durations are random and
perhaps statistically dependent. Note that in all these examples the integration is over
many variables and would be impossible by conventional numerical methods, even when
of day n and let H
n
K
be the number of new hires during day n. To simplify notation
the dependence upon K will be dropped for the moment.
The problem is to find the optimal value of K. It is known that each new hire generates
a fixed revenue of £c
f
per skip and a variable revenue of £c
v
per skip per day or part-day.
The K skips have to be bought at the outset and have to be maintained irrespective of
how many are on hire on a particular day. This cost is equivalent to £c
0
per skip per day.
Firstly, the stochastic process
X
n
n= 0 1
is considered. Assuming skips are
returned at the beginning of a day before hiring out,
Y
n
= Poisson
depends on X
n
and Y
n+1
only, and since Y
n+1
is independent of X
n−1
X
n−2
,
it follows that
X
n
n= 0 1
is a discrete-state, discrete-time, homogeneous Markov
chain. The probability transition matrix is P =
p
ij
ij= 0K
where p
ij
=
P
X
ij
=
min
ij
r=0
i
r
1 −p
r
p
i−r
j−r
e
−
j −r
!
(1.8)
For the case j = K,
p
iK
H = j
as n →. Let
denote the row vector
0
K
. is the stationary distribution of the chain
X
n
n= 0 1
. Suppose we wish to maximize the long-run average profit per day.
Then a K is found that maximizes Z
K
where
Z
K
= c
f
E
H
n
K
= lim
n→
E
X
n
K
−E
binomial
X
n−1
K
1 −p
= E
X
K
pc
f
+c
v
−c
0
K
The first difference of Z
K
is
D
K
= Z
K +1
−Z
K
−E
X
K
is decreasing in K.
In that case Z
K
will have a unique maximum.
A simulation using Maple 7
To determine the optimal K we require E
X
K
for successive integers K.IfK is
large this will involve considerable computation. Firstly, the probability transition matrix
would have to be calculated using Equations (1.8) and (1.9). Then it is necessary to invert
a
K +1
pn
x=X
0
K
For i =1n
y=Poisson
r=binomial
x 1 −p
x=min
r +yK
Output x
Next i
Now E
X
K
=
K
. There is a minor inconvenience
in that unless X
0
K
is selected randomly from , the stationary distribution, then
X
1
K
X
n
K
will not be precisely from , and therefore the estimate will be
slightly biased. However, this may be rectified by a burn-in period of b observations,
say, followed by a further n −b observations. We then estimate E
X
K
K +1
−E
X
K
with some
precision. We can reduce the sampling variation in our estimate of this by inducing positive
correlation between our estimates of E
X
K +1
and E
X
K
. Therefore, we might
consider making Y
n
‘skipexample.mws’. It explores the problem considered in Section 1.3. It can also be
downloaded from the website accompanying the book.
It will now be shown how the algorithm for simulation of skip hires can be
programmed as a Maple procedure. Before considering the procedure we will start with
8 Introduction to simulation and Monte Carlo
a fresh worksheet by typing ‘restart’ and also load the statistics package by
typing ‘with(stats)’. These two lines of input plus the Maple generated response,
‘[anova, , transform]’ form an execution group. In the downloadable file this is
delineated by an elongated left bracket, but this is not displayed here.
> restart;
with (stats);[anova, describe, fit, importdata, random, statevalf,
statplots, transform]
A Maple procedure is simply a function constructed by the programmer. In this case
the name of the procedure is skip and the arguments of the function are pK x0 and
n, the number of days to be simulated. Subsequent calling of the procedure skip with
selected values for these parameters will create a sequence; hire1hiren where
hirei =i x. In Maple terminology hirei is itself a list of two items: i the day number
and x the total number of skips on hire at the end of that day. Note how a list is enclosed
by square brackets, while a sequence is not.
Each Maple procedure starts with the Maple prompt ‘>’. The procedure is written
within an execution group. Each line of code is terminated by a semicolon. However,
anything appearing after the ‘#’ symbol is not executed. This allows programmer
comments to be added. Use the ‘shift’ and ‘return’ keys to obtain a fresh line within
the procedure. The procedure terminates with a semicolon and successful entry of the
procedure results in the code being ‘echoed’ in blue type. The structure of this echoed
code is highlighted by appropriate automatic indentation.
> skip: =proc (lambda, p, K, x0, n) local x, i, y, r, hire;
randomize (5691443); # An arbitrary integer sets the ‘seed’
for the U(0, 1) random number generator.
x:=x0;#xisthetotal number on hire, initially set to x0
end proc
The ‘echoed’ code is now examined. Local variables are those whose values cannot be
transmitted to and from the procedure. The ‘randomize’ statement uses an arbitrarily
chosen integer to set a seed for a U0 1 random number generator within Maple. More
will be said about such generators in Chapter 2. Initially, the total number of skips
x0 on hire in the current day (day 0) is assigned to the variable x. This is followed
by a ‘for i to…do…end do’ loop. The statements within this loop are executed for
i = 1n. Maple contains a ‘stats’ package and ‘random’ is a subpackage of this.
Note the use of the function ‘poisson’ within this subpackage. The generated Poisson
variate is assigned to the variable y. Following this is an example of a conditional
‘if…then…else…endif’ statement. If the total number of skips on hire on the previous
day is zero, then the number remaining on hire today r must be zero; otherwise r equals
a random binomial variate with parameters x and 1 −
p. Following this, the value of
the variable x is updated. Then the list i x is assigned to the variable hirei. The last
executable statement forms a sequence hire1hiren. Maple procedures use the
convention that the result of the last statement that is executed is assigned to the name of
the procedure. So when skip is subsequently called it will output a random realization of
this sequence.
Some results will now be obtained when = 5 per day, p = 02K = 30x0 = 0,
and n = 100 days. Calling ‘skip’ and assigning the results to a variable res gives the
sequence res.
> res :=skip (5, 0.2, 30, 0, 100);
res := [1, 6], [2, 6], [3, 13], [4, 17], [5, 16], [6, 19], [7, 25], [8, 25], [9, 24], [10, 27],
[11, 29], [12, 30], [13, 30], [14, 28], [15, 25], [16, 30], [17, 26], [18, 30], [19, 26],
[20, 23], [21, 25], [22, 25], [23, 28], [24, 30], [25, 24], [26, 24], [27, 28], [28, 24],
[29, 21], [30, 22], [31, 21], [32, 20], [33, 22], [34, 25], [35, 23], [36, 29], [37, 29],
[38, 29], [39, 26], [40, 25], [41, 27], [42, 27], [43, 30], [44, 30], [45, 27], [46, 30],
[47, 25], [48, 23], [49, 21], [50, 16], [51, 20], [52, 20], [53, 18], [54, 22], [55, 23],
[56, 26], [57, 25], [58, 29], [59, 27], [60, 26], [61, 30], [62, 27], [63, 29], [64, 27],
Copyright: Tài liệu đại học ©