Loading [MathJax]/jax/output/HTML-CSS/jax.js

Techniques

In this section, we outline simple ways to simulate random variables and random procecces. At the beginning of each subsection, we indicate to which section of the main theoretical text it corresponds.

1. Bernoulli and binomial random variables

We can simulate a Bernoulli random variable Z, which takes values 1 and 0 with probabilities p and 1p, as follows: Simulate UUni[0,1]. If Up, then Z:=1; else Z:=0. A binomial random variable XBin(N,p) can be simulated as the sum of N i.i.d. Bernoulli random variables: X=Z1++ZN.

2. Discrete random variables

Consider a random variable X which takes values 0,1,2, with probabilities p0,p1,p2, Split [0,1] into [0,p0], [p0,p0+p1], [p0+p1,p0+p1+p2],  Simulate UUni[0,1] and consider these intervals one by one: If 0Up0, then let X:=0; if p0Up0+p1, then let X:=1; etc.

3. Joint discrete distributions

Similarly we can simulate jointly distributed discrete random variables. As an example, simulate (X,Y) with joint distribution P(X=Y=0)=0.5, P(X=0,Y=1)=0.3, P(X=1,Y=0)=0.2. Simulate UUni[0,1]. Next, {0U0.5 X:=0,Y:=0;0.5U0.8  X:=0,Y:=1;0.8U1  X:=1,Y:=0.

4. Continuous random variables: the inverse function method

Consider a random variable X with cumulative distribution function F(x):=P(Xx), with density p(x)=F(x). Let us simulate it. Construct the inverse function: F1(u):=x, which solves F(x)=u. Note that F is increasing, so F1(u)x if and only if uF(x). Take a uniform [0,1] random variable U. Then we can take X=F1(U), because it has the correct cumulative distribution function: P(F1(U)x)=P(UF(x))=F(x), because for this uniform random variable U we have: P(Ua)=a,  0a1.

Example 1.

Try to simulate XExp(λ). It has density p(x)=λeλx and cumulative distribution function F(x)=x0p(y)dy=x0λeλydy=(eλy)|y=xy=0=1eλx,  x0. (And F(x)=0 for x0.) Now, let us find the inverse function. To this end, let us solve F(x)=u: 1eλx=u  eλx=1u  x=λ1ln(1u). Now, let X=λ1ln(1U)

5. Rejection method

Sometimes it is too difficult to simulate X using the inverse function method. For example, let XN(0,1), then F(x)=12πxey2/2dy, which is difficult to calculate and even more difficult to invert. You can do this, using tables, but there is a better way. Let f(x) be the density of X. In this case, it is f(x)=12πex2/2. It is bounded from above by c=(2π)1/2. Moreover, 10X10 with very high probability. Then we can simulate independent XUni[10,10] and YUni[0,c]. The point (X,Y) is then uniformly distributed in the rectangle [10,10]×[0,c]. We accept X if this point lies below the density curve: that is, if f(X)Y. Otherwise, we reject this point and repeat this again, until we get an acceptable point (X,Y). This method works equally well for multivariate distributions.

Example 2.

Simulate (X,Y) with density p(x,y)=x+y for x,y[0,1]. Then p2. Simulate independent XUni[0,1], YUni[0,1], ZUni[0,2], and accept (X,Y) if p(X,Y)=X+YZ.

6. Value at risk

Assume we have i.i.d. claims X1,,Xn, and we would like to find value-at-risk corresponding to the level of confidence 95% by simulation. Simulate these claims, rank them from top to bottom, and find the .05n-th ranked value. This is the value-at-risk for the current simulation. Repeat this a lot of times and average the resulting values-at-risk.

7. Monte Carlo method

Assume you have an integral I=10f(x)dx, which you cannot calculate directly. One way to approximate this integral is to split [0,1] into small subintervals and use Simpson’s rule or other similar rules. Another way is to generate N i.i.d. uniform random variables X1,,XN on [0,1] and let IN:=1NNk=1f(Xk). By Law of Large Numbers, because Ef(Xk)=10f(x)dx, and f(X1),,f(XN) are i.i.d., we have: INI. More generally, if we have I:=f(x)p(x)dx, where p is some probability density, let us generate N i.i.d. random variables X1,,XN with density p, and let IN:=1NNk=1f(Xk)I,  N, because Ef(Xk)=f(x)p(x)dx=I. This convergence is with speed N1/2, as follows from the Central Limit Theorem.

Example 3.

Let I=ex4dx. Generate X1,,XNN(0,1), and represent I as I=2πf(x)p(x)dx,f(x)=ex4+x2/2,p(x)=12πex2/2. Then we can approximate IIN:=1NNk=1f(Xk).

This method works equally well for multiple integrals. This is a particular advantage of Monte Carlo, because the usual methods become slow when the dimension grows curse of dimensionality

Example 4.

Consider the following integral on the 10-dimensional block [0,1]10 in R10: I=101010ex21+x22++x210dx1dx2dx10. To split it into subblocks [0,0.1]10 and similar, you need 1010 blocks. But it is just as easy to use Monte Carlo for this case as for dimension one. Generate i.i.d. 10-dimensional vectors Ui=(Ui1,,Ui10), with each component i.i.d. uniform on [0,1]. Then let IN:=1NNi=1exp(U2i1+U2i2++U2i10)I.

8. Discrete-time Markov chains

Consider the example with McDonalds-Subway. Let 0 and 1 correspond to McDonalds and Subway. Assume that the initial distribution is x(0)=[0.70.3]. Let U0Uni[0,1]. If U00.7, we let X0:=0; otherwise X0:=1. Recall the transition matrix: P=[0.80.20.10.9] Simulate U1Uni[0,1] independently of U0. If X0=0, let X1:=0 if U10.8, X1:=1 otherwise. If X0=1, let X1:=0 if U10.1, and X1:=1 otherwise. If this Markov chain has more than two (say four) states, then we should consider four cases for U1. In the same way, we simulate X2,X3,

Assume the Markov chain is irreducible. Then the occupation time shares converge to the stationary distribution π. One can empirically calculate the stationary distribtuion: Fix a large number N of steps. For each state i, find the number Mi of times n=1,,N when Xn=i. Then pi:=Mi/N is an approximation for the stationary distribution π.

9. Google’s PageRank

Take an oriented graph (edges are arrows), where states are Web pages and arrows are hyperlinks. Let us rank pages. A page is ranked highly if there are a lot of links pointing to it. However, the importance of each link is different. A link from an important Web page is also important. But if this important Web page provides lots of links, each of them is not so important. (A recommendation letter from Bill Gates carries more weight than the one from elementary school teacher. However, if Bill Gates is generous and writes thousands of letters every year, maybe his letter does not carry much weight after all.)

Let r(i) be the rank of the ith Web page. It is called PageRank in honor of Larry Page (who founded Google with Sergey Brin as Ph.D. students at Stanford in 1998). Let |i| be the number of links from page i. Then r(i)=jir(j)|j|. This is a mathematical formulation of the idea above. If we normalize r(i): divide it by R=r(i), the sum over all existing Web pages, then we have a probability distribution r(i) on the set of all Web pages. However, it is very hard to solve this system of equation explicitly. Brin and Page invented a Markov chain, which they called a random crawler: if at a certain moment it is at Web page i, it chooses randomly (uniformly) among the links from this page and moves along this link. This is, in effect, the random walk on a (oriented) graph. The transition matrix A of this Markov chain is given by Aij={1|i|, ij;0, otherwise. But we can rewrite the equation for r above as r(i)=jir(j)Aji=jr(j)Aji    r=rA. Therefore, r is a stationary distribution for A. If we run this random crawler long enough, it converges to r. (In fact, we need to make some adjustments to make sure this process indeed converges.)

10. Random walk

Consider a random walk with p=0.75, starting from S0=1. We simulate i.i.d. random variables Z1,Z2, with distribution P(Zk=1)=.75, P(Zk=1)=.25. We do this by taking i.i.d. U1,U2,Uni[0,1] and letting Zk=1 if Uk.75, Zk=1 if U>.75. Then let S0:=1, S1=S0+Z1, S2=S1+Z2,

11. Poisson process

Simulate i.i.d T1,T2,Exp(λ) and let τ1:=T1,τ2:=τ1+T2,τ3:=τ2+T3, Then we find N(t) by comparing t and τk until t<τk. If k is the first such that t<τk, then N(t)=k1. For example, if τ3t, but t<τ4, then N(t)=3. For a compound Poisson process X(t)=N(t)k=1Zk, we aslo need to simulate i.i.d. steps Z1,Z2, and let X(t)=Z1++Zk1 (and X(t)=0 in case k=0).

12. Continuous-time Markov chains

The difference from the discrete-time case is that we cannot simulate the value X(t) for every time t. Instead, we do exactly as for the Poisson process: We simulate the jump times, and the values of this process at these jump times. Essentially, we switch from continuous to discrete time. As an example, consider a continuous-time Markov chain X=(X(t),t0) with generator A=[312011112] The corresponding discrete-time Markov chain Y=(Yn)n0 has transition matrix P=[01/32/30010.50.50] Simulate Y0=X(0) as shown above. Assume X(0)=Y0=1. Then the Markov chain will spend time T1Exp(3) in this state. After exiting this state, it will go to 2 with probability 1/3 and to 3 with probability 2/3. One can simulate Y1 as shown in the subsection devoted to discrete-time Markov chains: Let U1Uni[0,1]; if U12/3 then Y1:=3; otherwise Y1:=2. Assume Y1=3. Then the Markov chain will spend time T2Exp(2) in this state. After exiting this state, it will go to 1 and 2 with equal probabilities .5. Let τ1:=T1 and τ2:=τ1+T2. In this way we can simulate jump times τk and jump positions Yk. For every fixed t0, find the k such that τkt<τk+1 and assign X(t):=Yk.

Similarly to discrete-time Markov chains, we can estimate the stationary distribution for irreducible continuous-time Markov chains by calculating the share of time spent in each state over a long time interval. To this end, we fix a large number N of jumps. Take each state i, and find all k=1,,N such that Yk=i. Sum Tk+1 for all such k, and divide the sum over N. The result will be pi, and [p1p2pN] is an approximation to the stationary distribution.

13. Brownian motion

It can be approximated by a discrete-time random walk: Take t0 to be the step size. Then W(kt0)t0(Z1++Zk), where Z1,Z2, are i.i.d. taking values ±1 with equal probability. This follows from the Central Limit Theorem: EZk=0 and VarZk=1, therefore Z1++ZkkN(0,1). And Z1++ZkkN(0,1)=N(0,k). But W(kt0)N(0,t0k)=N(0,k)t0. Thus, the simulation method is as follows: take a small t0, say t0=.01, and let W((k+1)t0):=W(kt0)+t0Zk+1, with initial value W(0):=0.

14. Stochastic integrals

We approximate them by taking a small step t0, as in the previous subsection: 20W(s)dW(s)2t01k=0W(kt0)(W((k+1)t0)W(kt0)). Having simulated the Brownian motion as in the previous subsection, we can then calculate the sum in the right-hand side.

15. Stochastic differential equations

Similarly to the Brownian motion, we take a small time step t0: dX(t)=3X(t)dt+X(t)dW(t), X(0)=1. We approximate this equation as follows: X((k+1)t0)X(kt0)=3X(kt0)t0+X(kt0)Zk+1, k=0,1,2, where Zk+1:=W((k+1)t0)W(kt0) can be approximated by t0Zk+1, with Zk+1=±1 with equal probabilities. Start from X(0)=1, simulate Z1 and calculate X(t0), then simulate Z2 and calculate X(2t0), etc.

Problems

Problem 1.

Simulate the random variables (X,Y) 10000 times, if they take values (1,0), (0,1), (1,0), (0,1) with equal probability 0.25. Calculate the empirical covariance Cov(X,Y), and compare with the theoretical value.

Problem 2.

Simulate the die flip 10000 times, and find the empirical mean and variance.

Problem 3.

Simulate a geometric random variable 10000 times with probability of success p=20%, by simulating a sequence of Bernoulli trials and then taking the number of the first successful trial. Compute empirical mean and variance, and compare with theoretical values.

Problem 4.

Simulate 1000 times a Poisson random variable Poi(1.5) by splitting [0,1] according to this distribution. Find its empirical mean, and compare with theoretical value. Simulate also Bin(100,0.15) by 100 Bernoulli trials at a time, repeat 1000 times. Compare empirical distributions with each other: find the difference |P(Poi(1.5)=k)P(Bin(100,0.15)=k)|,,k=0,1,2,3,4. For each distribution, find its empirical mean, and compare with theoretical value.

Problem 5.

Simulate 10000 times the random variable X with density p(x)=2π(1+x2), x0, using the inverse function method. Calculate an empirical probability P(X2), compare with theoretical value.

Problem 6.

Simulate 1000 times the Gamma variable Γ(3,2.6), using that it is the sum of three exponential random variables, which can be generated by the inverse function method. Compare the empirical and theoretical mean and variance.

Problem 7.

Simulate 1000 times the Gamma variable Γ(3,2.6) using rejection method. Cut its tail at some large enough point. Compare the empirical and theoretical mean and variance.

Problem 8.

Simulate the following jointly distributed random variables 10000 times using the rejection method: p(x,y)=y, xyx+1,0x1. Calculate empirical covariance, and compare with the theoretical result.

Problem 9.

Using Monte Carlo method with 10000 tries, simulate the double integral 1010(x+y)dxdy.

Problem 10.

Simulate 100 i.i.d. claims of size ZExp(1.5). Repeat this 1000 times. Find the value-at-risk for confidence level 95%. Compare with theoretical value.

Problem 11.

Simulate 10000 steps of a Markov chain with transition matrix P=[0.50.30.20.40.600.10.90] starting from 1. Find the stationary distribution by calculating the share of time spent in each state for this simulation. Compare with theoretical result.

Problem 12.

Simulate 10000 steps of a Markov chain with transition matrix P=[0.50.30.10.10.30.60.100.20.200.600.20.80] starting from 1. Find empirical probability of hitting 2 before 3. Compare with theoretical result.

Problem 13.

Simulate 1000 times the random walk with p=0.3 and q=0.7, starting from 2. Find the empirical probability that it hits 1 before 4. Compare with theoretical value.

Problem 14.

Simulate 10 steps of a random walk with p=0.6 and q=0.4, starting from 1. Repeat this simulation 10000 times. Find the empirical probability P(S10=3). Compare with theoretical value.

Problem 15.

Simulate a stock price 1000 times given by geometric symmetric random walk, starting from P0=1. Find the empirical fair value of the European option call at maturity T=30 with strike K=2.

Problem 16.

Simulate 1000 times the first 50 jumps of a Poisson process N with intensity λ=2. Calculate the empirical expectation E[N(1)N(2)], and compare it with the true value.

Problem 17.

Simulate 1000 times the first 50 jumps of a compound Poisson process X with increments N(2.5,4), and intensity λ=0.5. Use this to find empirical value VarX(4) and compare this with the true value.

Problem 18.

An insurance company receives claims during the year according to the Poisson process N(t) with intensity λ=2. Each claim is distributed as Exp(2.3). We measure time in days (365 days in a year). From 100 simulations, find the value-at-risk for confidence level 90%, and compare with theoretical value.

Problem 19.

Until time T=100, simulate a continuous-time Markov chain with generator A=[312253505] starting from state 1. Find the stationary distribution by calculating the share of time spent in each state. Compare with the theoretical value.

Problem 20.

Simulate 1000 times, starting from state 2, a continuous-time Markov chain with generator A=[3322] Find the empirical distribution of this chain at time t=2. Compare with the theoretical value.

Problem 21.

Simulate the M/M/1 queue with arrival intensity λ=2 and service intensity μ=4, until time horizon T=100. Find the empirical stationary distribution by calculating the share of time spent in each state. Compare this with the theoretical value.

Problem 22.

Simulate a Brownian motion with step size .01 for time horizon T=10. Repeat this 1000 times. Calculate the empirical probability P(W(10)>2), and compare with true value.

Problem 23.

Simulate a geometric Brownian motion X with drift g=.3 and variance σ2=2, starting from X(0)=0.3, until for time horizon T=1, with step size .01. Repeat this 1000 times. Calculate the empirical value of EX2(1), and compare with the true value.

Problem 24.

Find an empirical fair price of a barrier option which pays 1 if the price of a stock at some time until maturity T=2.5 exceeds K=2, if the initial price is S0=1.4, and the volatility is σ=0.4. Use step size 0.01 and 1000 simulations.

Problem 25.

Simulate the following It^o integral 1000 times with step size 0.01: 20eW(s)dW(s). Calculate the empirical variance and compare with the true result.

Problem 26.

Simulate the following Ornstein-Uhlenbweck process: dX(t)=2X(t)dt+3dW(t), X(0)=2, with step 0.01, for time horizon T=1. Repeat this 1000 times. Find empirical mean EX(1), and compare with the true value.