IEOR E4703: Monte Carlo Simulation Columbia University Generating .

Generating Random Variables and Stochastic Processes7Example 11 (Splitting)Suppose161[ 1,0] (x) 1[0,2] (x).515How do we simulate a value of X using vertical splitting? How would horizontal splitting work?fx (x) 2.3The Acceptance-Rejection AlgorithmLet X be a random variable with density, f (·), and CDF, Fx (·). Suppose it’s hard to simulate a value of Xdirectly using either the inverse transform or composition algorithm. We might then wish to use theacceptance-rejection algorithm. Towards this end let Y be another random variable with density g(·) andsuppose that it is easy to simulate a value of Y . If there exists a constant a such thatf (x) a for all xg(x)then we can simulate a value of X as follows.The Acceptance-Rejection Algorithmgenerate Y with PDF g(·)generate Uf (Y )while U ag(Y)generate Ygenerate Uset X YQuestion: Why must we have a 1?We must now prove that this algorithm does indeed work. We define B to be the event that Y has beenaccepted in the while loop, i.e., U f (Y )/ag(Y ). We need to show that P(X x) Fx (x)Proof: First observeP(X x) P(Y x B)P ((Y x) B).P(B) We can compute P(B) as P(B)while the numerator in (1) satisfiesZP ((Y x) B) f (Y ) P U ag(Y ) 1a P ((Y x) B Y y) g(y) dy f (Y ) Y y g(y) dyP (Y x) U ag(Y ) Z xf (y) g(y) dy (why?)P U ag(y) Z Fx (x)a (1)

Generating Random Variables and Stochastic Processes8Therefore P(X x) Fx (x), as required.Example 12 (Generating a Beta(a, b) Random Variable)Recall that X has a Beta(a, b) distribution if f (x) cxa 1 (1 x)b 1wish to simulate from the Beta(4, 3) so thatf (x) 60x3 (1 x)2for 0 x 1. Suppose now that wefor 0 x 1.We could, for example, integrate f (·) to find F (·), and then try to use the inverse transform approach.However, it is hard to find F 1 (·). Instead, let’s use the acceptance-rejection algorithm:1. First choose g(y): let’s take g(y) 1 for y [0, 1], i.e., Y U (0, 1)2. Then find a. Recall that we must havef (x) a for all x,g(x)which implies60x3 (1 x)2 a for all x [0, 1].So take a 3. It is easy to check that this value works. We then have the following algorithm.Algorithmgenerate Y U (0, 1)generate U U (0, 1)while U 20Y 3 (1 Y )2generate Ygenerate Uset X YEfficiency of the Acceptance-Rejection AlgorithmLet N be the number of loops in the A-R algorithm until acceptance, and as before, let B be the event that Yhas been accepted in a loop, i.e. U f (Y )/ag(Y ). We saw earlier that P(B) 1/a.Questions:1: What is the distribution of N ?2: What is E[N ]?How Do We Choose a?E[N ] a, so clearly we would like a to be as small as possible. Usually, this is just a matter of calculus.Example 13 (Generating a Beta(a, b) Random Variable continued)Recall the Beta(4, 3) example with PDF f (x) 60x3 (1 x)2 for x [0, 1]. We chose g(y) 1 for y [0, 1]so that Y U (0, 1). The constant a had to satisfyf (x) ag(x)for all x [0, 1]and we chose a 3. We can do better by choosingf (x) 2.073.x [0,1] g(x)a max

Generating Random Variables and Stochastic Processes9How Do We Choose g(·)?We would like to choose g(·) to minimize the computational load. This can be achieved by taking g(·) ‘close’ tof (·). Then a will be close to 1 and so fewer iterations will be required in the A-R algorithm. There is a tradeoff,however: if g(·) is ‘close’ to f (·) then it will probably also be hard to simulate from g(·). So we often need tofind a balance between having a ‘nice’ g(·) and a small value of a.Acceptance-Rejection Algorithm for Discrete Random VariablesSo far, we have expressed the A-R algorithm in terms of PDF’s, thereby implicitly assuming that we aregenerating continuous random variables. However, the A-R algorithm also works for discrete random variableswhere we simply replace PDF’s with PMF’s. So suppose we wish to simulate a discrete random variable, X,with PMF, pi P(X xi ). If we do not wish to use the discrete inverse transform method for example, thenwe can use the following version of the A-R algorithm. We assume that we can generate Y with PMF,qi P(Y yi ), and that a satisfies pi /qi a for all i.The Acceptance-Rejection Algorithm for Discrete Random Variablesgenerate Y with PMF qigenerate UpYwhile U aqYgenerate Ygenerate Uset X YGenerally, we would use this A-R algorithm when we can simulate Y efficiently.Exercise 1 (From Simulation by Sheldon M. Ross)Suppose Y Bin(n, p) and that we want to generate X whereP(X r) P(Y r Y k)for some fixed k n. Assume α P(Y k) has been computed.1. Give the inverse transform method for generating X.2. Give another method for generating X.3. For what values of α, small or large, would the algorithm in (2) be inefficient?Example 14 (Generating from a Uniform Distribution over a 2-D Region)Suppose (X, Y ) is uniformly distributed over a 2-dimensional area, A. How would you simulate a sample of(X, Y )? Note first that if X U ( 1, 1), Y U ( 1, 1) and X and Y are independent then (X, Y ) is uniformlydistributed over the regionA : {(x, y) : 1 x 1, 1 y 1}.We can therefore (how?) simulate a sample of (X, Y ) when A is a square. Suppose now that A is a circle ofradius 1 centered at the origin. How do we simulate a sample of (X, Y ) in that case?Remark 1 The A-R algorithm is an important algorithm for generating random variables. Moreover it can beused to generate samples from distributions that are only known up to a constant. It is very inefficient inhigh-dimensions, however, which is why Markov Chain Monte Carlo (MCMC) algorithms are required.

Generating Random Variables and Stochastic Processes310Other Methods for Generating Univariate Random VariablesBesides the inverse transform, composition and acceptance-rejection algorithms, there are a number of otherimportant methods for generating random variables. We begin with the convolution method.3.1The Convolution MethodSuppose X Y1 Y2 . . . Yn where the Yi ’s are IID with CDF Fy (·). Suppose also that it’s easy togenerate the Yi ’s. Then it is straightforward to generate a value of X:1. Generate Y1 , . . . , Yn that have CDF Fy2. Set X Y1 . . . YnWe briefly mentioned this earlier in Example 7 when we described how to generate a Gamma(λ, n) randomvariable. The convolution method is not always the most efficient method. Why?More generally, suppose we want to simulate a value of a random variable, X, and we know thatX g(Y1 , . . . , Yn )for some random variables Yi and some function g(·). Note that the Yi ’s need not necessarily be IID. If we knowhow to generate (Y1 , . . . , Yn ) then we can generate X by generating (Y1 , . . . , Yn ) and settingX g(Y1 , . . . , Yn ). We saw such an application in Example 7.Example 15 (Generating Lognormal Random Variables)Suppose X N(µ, σ 2 ). Then Y : exp(X) has a lognormal distribution, i.e., Y LN(µ, σ 2 ). (Note E[Y ] 6 µand Var(Y ) 6 σ 2 .) How do we generate a log-normal random variable?Example 16 (Generating χ2 Random Variables)Suppose X N(0, 1). Then Y : X 2 has a a chi-square distribution with 1 degree of freedom, i.e., Y χ21 .Question: How would you generate a χ21 random variable?Suppose now that Xi χ21 for i 1, . . . , n. Then Y : X1 . . . Xn has a chi-square distribution with ndegrees-of-freedom, i.e., Y χ2n .Question: How would you generate a χ2n random variable?Example 17 (Generating tn Random Variables)Suppose X N(0, 1) and Y χ2n with X and Y independent. ThenXZ : qYnhas a a t distribution with n degrees of freedom , i.e., Z tn .Question: How would you generate a tn random variable?

Generating Random Variables and Stochastic Processes11Example 18 (Generating Fm,n Random Variables)Suppose X χ2m and Y χ2n with X and Y independent. Then XZ : m Ynhas an F distribution with m and n degrees of freedom, i.e., Z Fm,n .Question: How would you generate a Fm,n random variable?4Generating Normal Random VariablesWhile we typically rely on software packages to generate normal random variables for us, it is nonethelessworthwhile having an understanding of how to do this. We first note that if Z N(0, 1) thenX : µ σZ N(µ, σ 2 )so that we need only concern ourselves with generating N(0, 1) random variables. One possibility for doing thisis to use the inverse transform method. But we would then have to use numerical methods since we cannot findFz 1 (·) : Φ 1 (·) in closed form. Other approaches for generating N(0, 1) random variables include:1. The Box-Muller method2. The Polar method3. Rational approximations.There are many other methods such as the A-R algorithm that could also be used to generate N (0, 1) randomvariables.4.1The Box Muller AlgorithmThe Box-Muller algorithm uses two IID U (0, 1) random variables to produce two IID N(0, 1) random variables.It works as follows:The Box-Muller Algorithm for Generating Two IID N(0, 1) Random Variablesgenerate U1 and U2 IID U (0, 1)setppX 2 log(U1 ) cos(2πU2 ) and Y 2 log(U1 ) sin(2πU2 )We now show that this algorithm does indeed produce two IID N(0, 1) random variables, X and Y .Proof: We need to show that 2 2 1x1y exp f (x, y) exp 222π2πFirst, make a change of variables:R: θ: pX2 Y 2 Y 1tanX

Generating Random Variables and Stochastic Processes12so R and θ areppolar coordinates of (X, Y ). To transform back, note X R cos(θ) and Y R sin(θ). Notealso that R 2 log(U1 ) and θ 2πU2 . Since U1 and U2 are IID, R and θ are independent. Clearly2θ U (0, 2π) so fθ (θ) 1/2π for 0 θ 2π. It is also easy to see that fR (r) re r /2 for r 0, so thatfR,θ (r, θ) 1 r2 /2re, 0 θ 2π, r 0.2πThis impliesP(X x1 , Y y1 ) P(R cos(θ) x1 , R sin(θ) y1 )Z Z1 r2 /2 redr dθA 2π(2)where A {(r, θ) : r cos(θ) x, r sin(θ) y}. We now transform back to (x, y) coordinates withx r cos(θ) and y r sin(θ)and note that dx dy rdr dθ, i.e., the Jacobian of the transformation is r. We then use (2) to obtain Z x1 Z y1(x2 y 2 )1exp dxdyP(X x1 , Y y1 ) 2π 2Z x1Z y111exp( x2 /2) dx exp( y 2 /2) dy 2π 2π as required.4.2The Polar MethodOne disadvantage of the Box-Muller method is that computing sines and cosines is inefficient. We can getaround this problem using the polar method which is described in the algorithm below.The Polar Algorithm for Generating Two IID N(0, 1) Random Variablesset S 2while S 1setgenerate U1 and U2 IID U (0, 1)set V1 2U1 1, V2 2U2 1 and S V12 V22rr 2 log(S) 2 log(S)V1 and Y V2X SSCan you see why this algorithm1 works?4.3Rational ApproximationsLet X N(0, 1) and recall that Φ(x) P(X x) is the CDF of X. If U U (0, 1), then the inverse transformmethod seeks xu Φ 1 (U ). Finding Φ 1 in closed form is not possible but instead, we can use rationalapproximations. These are very accurate and efficient methods for estimating xu .1 SeeSimulation by Sheldon M. Ross for further details.

Generating Random Variables and Stochastic ProcessesExample 19 (Rational Approximations)For 0.5 u 1xu t 13a0 a1 t1 b1 t b2 t 2pwhere a0 , a1 , b1 and b2 are constants, and t 2 log(1 u). The error is bounded in this case by .003. Evenmore accurate approximations are available, and since they are very fast, many packages (including Matlab) usethem for generating normal random variables.5The Multivariate Normal DistributionIf the n-dimensional vector X is multivariate normal with mean vector µ and covariance matrix Σ then we writeX MNn (µ, Σ).The standard multivariate normal has µ 0 and Σ In , the n n identity matrix. The PDF of X is given byf (x) 11(2π)n/2 Σ 1/2 e 2 (x µ)Σ 1 (x µ)where · denotes the determinant, and its characteristic function satisfiesh i 1 φX (s) E eis X eis µ 2 s Σs .(3)(4)Recall again our partition of X into X1 (X1 , . . . , Xk ) and X2 (Xk 1 , . . . , Xn ) . If we extend thisnotation naturally so that µ1Σ11 Σ12µ andΣ .µ2Σ21 Σ22then we obtain the following results regarding the marginal and conditional distributions of X.Marginal DistributionThe marginal distribution of a multivariate normal random vector is itself multivariate normal. In particular,Xi MN(µi , Σii ), for i 1, 2.Conditional DistributionAssuming Σ is positive definite, the conditional distribution of a multivariate normal distribution is also amultivariate normal distribution. In particular,X2 X1 x1 MN(µ2.1 , Σ2.1 ) 1where µ2.1 µ2 Σ21 Σ 111 (x1 µ1 ) and Σ2.1 Σ22 Σ21 Σ11 Σ12 .Linear CombinationsLinear combinations of multivariate normal random vectors remain normally distributed with mean vector andcovariance matrix given byE [AX a] AE [X] aCov(AX a) A Cov(X) A .

Generating Random Variables and Stochastic Processes14Estimation of Multivariate Normal DistributionsThe simplest and most common method of estimating a multivariate normal distribution is to take the samplemean vector and sample covariance matrix as our estimators of µ and Σ, respectively. It is easy to justify thischoice since they are the maximum likelihood estimators. It is also common to take n/(n 1) times the samplecovariance matrix as an estimator of Σ as this estimator is known to be unbiased.Testing Normality and Multivariate NormalityThere are many tests that can be employed for testing normality of random variables and vectors. These includestandard univariate tests and tests based on QQplots, as well omnibus moment tests based on whether theskewness and kurtosis of the data are consistent with a multivariate normal distribution. Section 3.1.4 of MFEshould be consulted for details on these tests.5.1Generating Multivariate Normally Distributed Random VectorsSuppose that we wish to generate X (X1 , . . . , Xn ) where X MNn (0, Σ). Note that it is then easy tohandle the case where E[X] 6 0. Let Z (Z1 , . . . , Zn ) where the Zi ’s are IID N(0, 1) for i 1, . . . , n. If C isan (n m) matrix then it follows thatC Z MN(0, C C).Our problem therefore reduces to finding C such that C C Σ. We can use the Cholesky decomposition of Σto find such a matrix, C.The Cholesky Decomposition of a Symmetric Positive-Definite MatrixA well known fact from linear algebra is that any symmetric positive-definite matrix, M, may be written asM U DUwhere U is an upper triangular matrix and D is a diagonal matrix with positive diagonal elements. Since Σ issymmetric positive-definite, we can therefore write Σ U DU (U D)( DU) ( DU) ( DU). The matrix C DU therefore satisfies C C Σ. It is called the Cholesky Decomposition of Σ.The Cholesky Decomposition in Matlab and RIt is easy to compute the Cholesky decomposition of a symmetric positive-definite matrix in Matlab and R usingthe chol command and so it is also easy to simulate multivariate normal random vectors. As before, let Σ be an(n n) variance-covariance matrix and let C be its Cholesky decomposition. If X MN(0, Σ) then we cangenerate random samples of X in Matlab as follows:Sample Matlab Code Sigma [1.0 0.5 0.5;0.5 2.0 0.3;0.5 0.3 1.5]; C chol(Sigma); Z randn(3,1000000); X C’*Z; cov(X’)ans 4971

Generating Random Variables and Stochastic Processes15We must be very careful in Matlab2 and R to pre-multiply Z by C and not C. We have the following algorithmfor generating multivariate random vectors, X.Generating Correlated Normal Random Variablesgenerate Z MN(0, I)/ Now compute the Cholesky Decomposition /compute C such that C C Σset X C Z6Simulating Poisson ProcessesRecall that a Poisson process, N (t), with intensity λ is a process such thatP (N (t) r) (λt)r e λt.r!For a Poisson process the numbers of arrivals in non-overlapping intervals are independent and the distributionof the number of arrivals in an interval only depends on the length of the interval.The Poisson process is good for modeling many phenomena including the emission of particles from aradioactive source and the arrivals of customers to a queue. The ith inter-arrival time, Xi , is defined to be theinterval between the (i 1)th and ith arrivals of the Poisson process, and it is easy to see that the Xi ’s are IID Exp(λ). In particular, this means we can simulate a Poisson process with intensity λ by simply generating theinter-arrival times, Xi , where Xi Exp(λ). We have the following algorithm for simulating the first T timeunits of a Poisson process:Simulating T Time Units of a Poisson Processset t 0, I 0generate Uset t t log(U )/λwhile t Tset I I 1, S(I) tgenerate Uset t t log(U )/λ6.1The Non-Homogeneous Poisson ProcessA non-homogeneous Poisson process, N (t), is obtained by relaxing the assumption that the intensity, λ, isconstant. Instead we take it to be a deterministic function of time, λ(t). More formally, if λ(t) 0 is theintensity of the process at time t, then we say that N (t) is a non-homogeneous Poisson process with intensityλ(t). Define the function m(t) byZ tm(t) : λ(s) ds.02 Unfortunately,C some languages taketo be the Cholesky Decomposition rather C. You must therefore always be awareof exactly what convention your programming language / package is using.

Generating Random Variables and Stochastic Processes16Then it can be shown that N (t s) N (t) is a Poisson random variable with parameter m(t s) m(t), i.e.,P (N (t s) N (t) r) exp ( mt,s ) (mt,s )rr!where mt,s : m(t s) m(t).Simulating a Non-Homogeneous Poisson ProcessBefore we describe the thinning algorithm for simulating a non-homogeneous Poisson process, we first needthe following3 proposition.Proposition 1 Let N (t) be a Poisson process with constant intensity λ. Suppose that an arrival that occursat time t is counted with probability p(t), independently of what has happened beforehand. Then the process ofcounted arrivals is a non-homogeneous Poisson process with intensity λ(t) λp(t).Suppose now N (t) is a non-homogeneous Poisson process with intensity λ(t) and that there exists a λ such thatλ(t) λ for all t T . Then we can use the following algorithm, based on Proposition 1, to simulate N (t).The Thinning Algorithm for Simulating T Time Units of a NHPPset t 0, I 0generate U1set t t log(U1 )/λwhile t Tgenerate U2if U2 λ(t)/λ thenset I I 1, S(I) tgenerate U1set t t log(U1 )/λQuestions1. Can you give a more efficient version of the algorithm when there exists λ 0 such that min0 t T λ(t) λ?2. Can you think of another algorithm for simulating a non-homogeneous Poisson process that is not based onthinning?6.2Credit Derivatives ModelsMany credit derivatives models use Cox processes to model company defaults. A Cox process, C(t), is similar toa non-homogeneous Poisson process except that the intensity func

Generating Random Variables and Stochastic Processes 4 The Inverse Transform Method for Continuous Random Variables Suppose now that Xis a continuous random variable and we want to generate a value of X. Recall that when Xwas discrete, we could generate a variate by rst generating Uand then setting X x j if F(x j 1) U F(x j). This suggests .