MonteCarlo Intro in Mtlab

8/10/2019 MonteCarlo Intro in Mtlab

1/26

An Introduction to

Monte Carlo Simulations

Notes by

Mattias Jonsson

Version of October 4, 2006

Math 472, Fall 2006

Numerical Methods with Financial Applications


2/26

2

Contents

1. Introduction 4

2. Computing areas, volumes and integrals 42.1. Examples 42.1.1. Example A: Area of a rectangle 42.1.2. Example B: Area of a disc 42.1.3. Example C: Area of a superellipse 42.1.4. Example D: Volume of ann-dimensional ball 42.1.5. Example E: Computing a triple integral 42.2. Methods 42.3. First method: exact computation 52.4. Second method: the lattice method 52.4.1. Example A 52.4.2. Examples B and C 52.4.3. Example D 62.4.4. Example E 62.5. Third method: the Monte Carlo method 72.5.1. Examples B and C 72.5.2. Example D 82.5.3. Example E 93. A crash course in probability theory 93.1. Outcomes, events and probabilities 93.2. Random variables 103.3. Discrete random variables 103.4. Continuous random variables 10

3.5. Expectation and variance 113.6. Some important distributions 113.6.1. Binomial random variables 123.6.2. Poisson random variables 123.6.3. Normal random variables 123.7. Random vectors 133.8. Independence 133.9. Conditional expectation and variance 134. Computing expected values 144.1. General strategy 144.2. Sampling 144.2.1. Sampling discrete random variables 15

4.2.2. Sampling continuous random variables 164.3. The law of large numbers 164.4. Presenting the results 174.5. Examples 184.5.1. Flipping a coin 184.5.2. Rolling dice 194.5.3. Normal random variables 19


3/26

3

4.5.4. Computing integrals 204.6. Variance reduction techniques 20

4.6.1. Antithetic variables 214.6.2. Moment matching 215. Examples from Financial and Actuarial Mathematics 225.1. Option prices in the Black-Scholes model 225.2. The collective loss model 235.3. The binomial tree model 245.3.1. Calibration 255.3.2. Pricing a barrier option 25


4/26

4

1. Introduction

These notes are meant to serve as an introduction to the Monte Carlo

method applied to some problems in Financial and Actuarial Mathematics.Roughly speaking, Monte Carlo simulation is a method for numerically

computing integrals and expected values. It is a probabilistic method atheart, based on simulating random outcomes or scenarios. Such a methodmay bring up the analogy with gambling at the casino. At any rate themethod is named after the casino in Monte Carlo, Monaco.

2. Computing areas, volumes and integrals

2.1. Examples. In order to introduce the Monte Carlo method, we startby giving a few examples of problems where the method applies. First welist the problems, then we discuss methods of solving them.

2.1.1. Example A: Area of a rectangle. We start by a seemingly trivial prob-lem: compute the area of the rectangle

= [0.2, 0.8] [0.1, 0.6]sitting inside the square [0, 1]2.

2.1.2. Example B: Area of a disc. The next example seems essentially aseasy: compute the area of the unit disc

={(x, y) | x2 + y2


5/26

5

2.3. First method: exact computation. Of course, if we can computean area, a volume or an integral exactly, then there is usually no need to

employ a numerical method. Consider the examples above. In Example A, the area equals 0.3. In Example B, the area is . In Example C, it is possible but not so easy to compute the exact

area. In Example D, it is also possible to find the exact answer. For exam-

ple, the answer is 1, and 4/3, 2/2 forn = 1, 2, 3, 4, respectively. In Example E, however, it is (as far a I know) not possible to compute

the exact answer as a closed formula.

2.4. Second method: the lattice method. The lattice method essen-tially exploits the definition of the (multiple, Riemann) integral. Let us

illustrate it on the five examples above.

2.4.1. Example A. Fix an integer J 1 and set N = J2. Divide the unitsquare [0, 1]2 intoNequally large subsquares Rk,k = 1, . . . , N . Each squarehas area 1/N = 1/J2. Let pk be the center of Rk. These centers have

coordinates (2i12J ,

2j12J ), where 1 i, j J. Now we use the approximation

Area() AN=pk

Area(Rk) = 1

N#{k |pk }

Exercise: Check that A16 = 0.25 but thatAN= 0.3 whenever Jis divisibleby 10.

2.4.2. Examples B and C. The lattice method can be used to compute thearea of quite general domains in the plane. We proceed as follows:

First find a rectangle Rthat you know contains the region . Pick a (large) integer J 1 and set N=J2. Divide the rectangle R into Nequally large subrectangles Rk, 1

k N. Let pk be the center ofRk. The centers form a lattice, explaining

the name of the method! Approximate the exact areaA= A() of by

AN= pk

Area(Rk) =Area(R)

N #{k | pk }.

The method is easy to implement as long as it is easy to decidewhether a given point belongs to or not. This is the case inExamples B and C.

It follows by the definition ofA = Area() that AN converges toAas N (i.e.J ).


6/26

6

One can check that ifis a nice domain (more precisely, if it hassufficiently smooth boundary), then the convergence rate is equal to

O( 1J

) =O( 1N

).

This convergence is not very fast. Basically, to get one more signifi-cant digit, we have to increase Nby a factor 100.

Here is a simple matlab function that uses the lattice method to computean approximation to the area of the unit disk.

function area=da_lattice(J)

v=(2*[1:J]-1)/J-1;

[x,y]=ndgrid(v);

area=4*sum(sum(x.^2+y.^2


7/26

7

Approximate the exact value of the integral I by

IN= pkVol(Rk) f(pk), =Vol(R)

N pk f(pk),wheref is the function f(x1, x2, x3) =e

x1+x2+x3 .

It should be clear how to modify the algorithm in order to compute moregeneral integrals

f(x) dx,

f(x1, x2) dx1dx2,

f(x1, x2, x3) dx1dx2dx3, . . .

The implementation is easy as long as we can evaluate the function f, anddecide whether a given point belongs to the domain .

It follows by the definition of the integral that INconverges toIasN

(i.e. J

). As before, for nice domains and nice functions f,

the convergence rate is equal to

O(1

J) =O(

1

N1/n),

wheren is the dimension.

2.5. Third method: the Monte Carlo method. The idea of the MonteCarlo method is quite simple. Instead of usingallpoints on a lattice, wepick them totally randomly. We illustrate the method on ExamplesB Eabove.

2.5.1. Examples B and C. Let us use the Monte Carlo method to compute

the are of a domain

in the plane. First find a rectangleR = RxRy that you know contains the region.

Pick a (large) integer N(not necessarily of the form N=J2. Pick 2N independent samples x1, . . . , xN and y1, . . . , yN from the

uniform distribution on Rx and Ry, respectively. Setpk = (xk, yk).Thenp1, . . . , pNare independent samples of points from the uniformdistribution ofR.

Approximate the exact area A= A() of exactly as before, thatis, by AN, where

AN= pkArea(Rk) =Area(R)

N

#{k | pk

}.

As before, the method is easy to implement as long as it is easy todecide whether a given point belongs to or not.

Notice that AN is a random variable. Different runs will give rise todifferent approximations!

We hope that AN converges (in some sense) as N . More onthis later.


8/26

8

The following simple matlab function computes an approximation to thearea of the unit disk using the Monte Carlo method.

function area=da_mc(N)

x=2*rand(1,N)-1;

y=2*rand(1,N)-1;

area=4*sum(x.^2+y.^2


9/26

9

2.5.3. Example E. Having worked through the previous examples, it shouldnow be straightforward to modify the lattice method for the multiple integral

in Example E to instead use the Monte Carlo method. We approximate the(multiple) integral

I=

f

by

IN=Vol(R)

N

pk

f(pk),

wherepkare points taken from a uniform distribution on a cube Rcontainingthe domain .

The following matlab function computes the volume of the unit ball indimension dim using the Monte Carlo method.

function int3d=int3d_mc(N)

x=2*rand(3,N)-1;

int3d=(2^3/N)*sum((sum(x.^4)


10/26

10

event in F, since it corresponds to knowing that the we rolled a 1, andnothing else.

A probability measure is a function P : F [0, 1], satisfying certainnatural conditions. It should be though of as a rule that assigns a likelihoodto an event. An event happens almost surely if it has probability one. Thetriple (, F, P) is called aprobability space.

3.2. Random variables. A random variable is, formally, a function X : R which is F-measurable, that is, for any x R, the set { | X() x} is in F.

The intuition behind this is that the information encoded in Fabout theexperiment allows us to tell us whether or not Xis less thanx. By propertiesof-algebras this means that we are able to tell many other things aboutX, for instance whether Xbelongs to a given interval.

For many purposes, what we need to know about a random variable isencoded in its distribution functionFX :R [0, 1], defined by

FX(x) =P{X x}.Notice that there is no explicit reference to orFhere.

A natural example is to think of Xas the price of an asset tomorrow.Given information today, we typically want to know the probability of theasset price tomorrow being, say, less that $100.

3.3. Discrete random variables. A random variable X is discrete if itcan take only finitely or countably many values x1, x2,.. . . The simplestnontrivial example is that when X is the number of heads after flipping

a fair coin once. In this case X takes only the values 0 and 1. A morecomplicated example is that of a Poisson random variable X, which cantake any nonnegative integer as its value.

Consider a discrete random variable Xtaking only the nonnegative inte-gers as values, i.e.X= 0, 1, 2, . . . . Then the distribution ofX is convenientlydescribed by theprobability mass function(pmf)fX. This function is simplydefined by

fX(n) =P{X=n}.

We can easily obtain the cdf (cumulative distribution function) from thepdf:

FX(x) = nx fX(n)3.4. Continuous random variables. Together with the discrete ones, thecontinuousrandom variables are the most important ones. While a discreterandom variable can only take (at most) countably many values, the prob-ability of a continuous random variable taking a specific value is alwayszero.They can be defined as random variables X admitting a probability


11/26

11

density function (pdf)fX. The (pdf) is related to the cumulative distribu-tion function by

FX(x) = x

fX(u) du.

Thus a random variable Xcan be though of as a random variable for whichthe cdf is (more or less) differentiable. By contrast, the cdf of a discreterandom variable is a step function.

3.5. Expectation and variance. Intuitively, it is natural to talk aboutthe typical or average value of a random variable X, as well as thesize of the typical fluctuation ofX around this value. The correspond-ing mathematical notions are those ofexpected value (or expectation) andvariance.

IfXis discrete, taking values X= 0, 1, 2, . . . then the expected value is

given by

E[X] =

n=0

nfX(n).

More generally, ifg : N R is a function, then g(X) is again a discreterandom variable (not necessarily integer-valued) with expectation

E[g(X)] =

n=0

g(n)fX(n).

When instead Xis continuous, the expected value is given by

E[X] = xfX(x) dx.More generally, ifg : RR is a function, then g(X) is a (not necessarilycontinuous) random variable with expectation

E[g(X)] =

g(x)fX(x) dx.

For an arbitrary random variable X, the varianceofXis defined by

Var[X] =E[(X E[X])2] =E[X2] E[X]2.In general, there is no guarantee that the integrals or sums defining the

exposed value or variance are convergent. There are random variables forwhich the expected value or variance are not defined. However, these randomvariables will not play an important role in these notes.

3.6. Some important distributions. There are many important familiesof distributions that are commonly used throughout applied mathematics.We shall make no attempt to list them all, but content ourselves with re-calling a few selected distributions.


12/26

12

3.6.1. Binomial random variables. A binomial random variable with pa-rameters (m, p) can be thought of as the number of heads after tossing m

independent biased coins, where the probability of heads isp for each singlethrow. IfXis binomial, then we write

X Bin(m, p).ClearlyXthen is discrete, takes values 0, 1, . . . , mand the pmf ofXis givenby

fX(k) =

m

k

pk(1 p)mk.

Sometimes a binomial random variable with parameter (1, p) is called aBernoullirandom variable. A binomial random variable can be seen as thesum ofm i.i.d. Bernoulli random variables.

A binomial random variable X Bin(m, p) has expected value and vari-ance given by

E[X] =mp and Var[X] =mp(1 p).Binomial random variables are often used to model the number of events

during a fixed time period.

3.6.2. Poisson random variables. Another distribution commonly used tomodel the number of events is the Poisson distribution. We write

X

Poisson()

for a Poisson random variable with parameter . Then X takes valuesX= 0, 1, 2, . . . and has pmf

fX(k) =e

k

k!.

A Poisson random variable X Poisson() has expected value E[X] =and Var[X] =.

It is possible to view the Poisson random variable as the limit of binomialrandom variables, but we shall not dwell on this here.

3.6.3. Normal random variables. A very important example is a normal (orGaussian) random variable, having density function

fX(x) = 1

22exp

(x m)

2

22

,

for some constantsm R and > 0. We write X N(m,2).IfX N(m,2), then E[X] =m and Var[X] =2.


13/26

13

3.7. Random vectors. Let X1, . . . , X n be random variables on the sameprobability space (, F, P). We may think of this as ann-dimensionalran-

dom vector X= (X1, . . . , X n).The joint distribution function of X is the function FX : Rn [0, 1]

defined byFX(x1, . . . , xn) =P{X1 x1, . . . , X n xn}.

We say that X is continuous if there exists a function fX 0, the jointdensity function, such that

FX(x1, . . . , xn) =

x1

. . .

xn

fX(u1, . . . , un) du1 . . . dun.

A normal random vector is specified by its with mean m = (m1, . . . , mn)and its positive definite variance-covariance matrix R= (Rij)1i,jn. It isa continuous random vector Xwith joint density function

fX(x) = 1

(2)n/2

detRexp

1

2(x m)TR1(x m)

.

3.8. Independence. There is an intuitive (but imprecise) notion of tworeal-world events happening independently of each other. In probabilitytheory, this is formulated as follows: given a probability space (, F, P),two events A, B F are independent ifP(A B) = P(A)P(B). For ourapplications, it is more important to talk about independent random vari-ables. We say that the random variables X1, . . . , X n are independent if, forany events Ai F, we have

P{X1 A1, . . . , X n An}= P{X1 A1} . . . P {Xn An}.This means that the joint distribution function FX of (X1, . . . , X n) is givenin terms of the individual distribution functions as

FX(x1, . . . , xn) =FX1(x1) . . . F Xn(xn).

IfX is continuous with density function fX, then this is in turn equivalentto

fX(x1, . . . , xn) =fX1(x1) . . . f Xn(xn).

3.9. Conditional expectation and variance. Let (, F, P) be a proba-bility space. IfB F is an event with nonzero probability, then the we candefine a new probability measure on (, F) by

P(A|B) =

P(A

B)

P(B) .

This is called the probability ofA conditioned onB , or theconditional prob-ability ofA givenB. IfA and B are independent, then P(A|B) =P(A).

The expectation of a random variable Xwith respect to this measure isgiven by

E[X|B] =E[X 1A]

P(B) ,


14/26

14

where 1A is the (F-measurable) random variable whose value at is 1 if Aand 0 if / A. The interpretation of this is the expected value ofXgiven that the event B occurred.IfX and Y are F-measurable random variables, then we want to defineE[X|Y], the conditional expectation ofX given Y, to be a function ofY.Its value at Y = y should mean the expected value ofX given that Y istaking value y. Since the probability of the latter event could be zero, onehas to be careful in defining E[X|Y]. One way is to define E[X|Y] to bethe functiong (Y) ofYthat minimizes the value ofE[(X g(Y))2].

We may also define the conditional varianceofXgiven Y :

Var[X|Y] =E[(X E[X|Y])2|Y] =E[X2|Y] E[X|Y]2.An important property is the law of iterated expectations, asserting that

E[X] =E[E[X|Y]].

A consequence of this formula is

Var[X] =E[Var[X|Y]] + Var[E[X|Y]].

4. Computing expected values

The Monte Carlo method, as far as we will use it, is a method for comput-ing (numerically) the expected value E[X] of a random variable X. Noticethat both the probability of an event and the variance of a random variablecan be viewed as the expectation of some other random variable, so theMonte Carlo method will allows us to estimate probabilities and variances,too.

We shall also see that the computations of area, volumes and multipleintegrals in Section 2.5 can be cast as computations of expected values.

4.1. General strategy. In principle, the Monte Carlo method is very sim-ple. To compute the expected value E[X] ofXone can do as follows:

(i) Pick a large number N(ii) Generate (somehow) N independent samples 1, 2, . . . ,N ofX.

(iii) Compute EN= 1

N

Nn=1 n.

(iv) We hope thatEN E[X] ifN is large enough.We shall discuss the convergence in (iv) (and the choice ofN) later on. Fornow we shall briefly discuss how to sample the random variables (ii) andhow to present and interpret the results from a Monte Carlo simulation.

4.2. Sampling. It can be tricky to sample for an arbitrary distribution. Inmatlab, there are two basic commands for generating samples of randomvariables:

The command rand samples from U(0, 1), the uniform distributionon the interval [0, 1].

The command randn samples from N(0, 1), the standard normaldistribution.


15/26

15

By adding arguments to these functions, matlab will produce matrices ofrandom numbers. For instance, rand(2,3) will produce 2x3 matrix whose

elements are independent samples fromU(0, 1).It should be noted the numbers produced by matlab or any similar soft-ware are not true random numbers but merely pseudo-randomnumbers. Wewill not discuss what this means here. The algorithms used by matlab seemfairly robust and tend to produce good result if used properly.

One more remark about matlab is in order. Every time you start a matlabsession, the seed used by the random number generator is reset. Thismeans that you may get the same result when running the same programseveral times. Occasionally, this may be desirable, but typically it is not agood thing. A standard way around this problem is to reset the seed usingthe computers clock. For instance, one can use the command

rand(state,sum(100*clock))This command should then be invoked at the start of the matlab session,

butnotevery time a particular program is run. Resetting the seed too oftenwill compromise the samples.

4.2.1. Sampling discrete random variables. One can use the function randin matlab to sample discrete random variables.

For example, to generate N independent samples of a Bernoulli randomvariable, i.e. a random variable X with pmffX(1) = p, fX(0) = 1 p forsome 0< p


16/26

16

4.2.2. Sampling continuous random variables. In principle one can sampleany continuous random variable X(with values in R) as follows:

X=F1X (Y), that is, Y =FX(X),

whereY is uniformly distributed on (0, 1) (hence can be sampled usingrandand FX :R [0, 1] is the cdf for X.

For example, when X is exponentially distributed with expected valueone, then

FX(x) = 1 exwhich implies

F1X (y) = log(1 y).In general, however, it can be difficult to apply this procedure in a com-

putationally effi

cient way, unless there is an explicit formula for the inversecdfF1X .A special (and important) case that matlab knows how to handle well is

that of a standard normal distribution. The command

x=randn(1,N)

generates a vector ofN independent samples from the standard normaldistribution function.

Sometimes one can do simple transformations to get new distributions.For instance,

x=mu+sigma*randn(1,N)

generates samples from the normal distribution N(,2) with mean and variance 2, and

x=exp(mu+sigma*randn(1,N))

generates samples from a log-normal distribution.

4.3. The law of large numbers. The reason why Monte Carlo works isthe Law of Large Numbers.

Theorem(Kolmogorovs Strong Law of Large Numbers). IfX1, X2, . . . , X n, . . .are i.i.d. random variables with (well-defined) expected value E[Xn] = ,

then

1

NNn=1 Xn a.s. asN .This version of the law of large numbers does not tell us how fast the

convergence is. In general, we have a rule of thumb: In Monte Carlo Simu-lations, the error isO( 1

N) when usingN samples.

Note: this is quite slow: to get 1 more decimal, need 100 times moresimulations!

The reason for this rule of thumb is


17/26

17

Lemma. If X1, X2, . . . , X n, . . . are i.i.d. random variables with expectedvalueE[Xn] = and varianceVar[Xn] =

2, then

E 1

N

Nn=1

Xn

= and Var

1N

Nn=1

Xn

=2

N

The formula for the variance implies that the random variable N1N

n=1 Xnhas standard deviation/

N.

While the strong law of large numbers is rather difficult to prove, thetwo formulas above are elementary. A more precise (and non-elementary!)version of the rule of thumb is given by

Theorem (Central Limit Theorem). If X1, X2, . . . , X n, . . . are i.i.d. ran-dom variables with mean and variance2: i.e. E[Xn] = andVar[Xn] =2, then N

1 Xn N

N N(0, 1) asN .

In other words, after subtracting the mean and dividing by the standarddeviation, the sum converges to a standard normal/Gaussian.

In particular we get

1

N

N1

Xn error

N(0,2

N)

so the error has standard deviation N

.

4.4. Presenting the results. The Monte Carlo method differs from mostother methods by its probabilistic nature. In particular, the approximationwe get will vary from simulation to simulation. For this reason it is importantto present the results of a Monte Carlo simulation in an instructive way.

It is a good habit to always include the following elements in a presenta-tion of the results of a Monte Carlo simulation:

N: the number of paths xN: the Monte Carlo estimate N: the standard error A convergence diagram

Let us explain the items above. We have a random variable X with(unknown) expected value E[X] = and variance Var[X] = 2 and weare interested in computing . To this end, we use independent samplesx1, . . . , xN fromXand form the average

xN= 1

N(x1+ x2+ + xN).

This is the Monte Carlo estimate ofreferred to above.


18/26

18

The standard error Nis supposed to measure the error in the Monte Carlosimulation. Of course, we cannot know the error exactly, or else there would

be no point in doing the Monte Carlo simulation in the first place. What wewould like to know is the approximate size of the error. Now, if we think ofxNas a random variable rather than a number, then by the Central LimitTheorem,xN is approximately distributed asN(,

2/

N). Unfortunately,we do not know ! However, we can use the following estimate for :

N=

1N 1

Ni=1

x2i Nx2N

The standard error is then defined as

N= N

N.

Finally, aconvergence diagramis a plot of xn againstn for 1 n N. Itis a good visual tool for determining whether a given Monte Carlo simulationis close to having converged. See the examples below for more details onhow to do this.

4.5. Examples. Let us illustrate the Monte Carlo method for approximat-ing expected values in a few simple cases. More interesting examples will begiven later on.

4.5.1. Flipping a coin. Let X be a Bernoulli random variable, i.e. X = 0or X= 1 with probability fX(0) = 1 p and fX(1) = p. ThenX can bethough of as the number of heads after a single toss of an unfair (ifp = 0.5)coin. We wish to use Monte Carlo simulations to approximate the expectedvalue ofXand compare it with the exact value E[X] =p.

We saw above how to sample Xusing matlab. The following matlab codegeneratesNsamples, computes a Monte Carlo estimate approximate valueofE[X], computes the standard error, and produces a convergence diagram.In the code, we have set p= 0.4 and N= 100000 but this can of course bemodified.

clear;

N=100000;p=0.4;

x=(rand(1,N)


19/26

19

figure(1);plot([Nmin:N],hx(Nmin:N));

title(MC convergence diagram: expectation of Bin(1,0.4));

xlabel(No paths);ylabel(estimate);

4.5.2. Rolling dice. Next we want to use Monte Carlo simulations to com-pute the probability of rolling a full house in one roll in the game ofYahtzee. This means rolling five (fair, independent, six-sided) dice and get-ting a pattern of type jjjkk, where 1 j, k 6 and j=k.

Recall that the probability of an even A is the same as the expected valueof the random variable 1A taking value 1 if A occurs and zero otherwise.We can use the following matlab code to compute a Monte Carlo estimateof the probability of a full house.

clear;

N=100000;

x=sort(ceil(6*rand(5,N)));

y=((x(1,:)==x(2,:)==x(3,:))&(x(4,:)==x(5,:))&(x(3,:)~=x(4,:)))|...

((x(3,:)==x(4,:)==x(5,:))&(x(1,:)==x(2,:))&(x(2,:)~=x(3,:)));

hy=cumsum(y)./[1:N];

eps=sqrt((sum(y.^2)-N*hy(N)^2)/(N-1))/sqrt(N);

fprintf(No paths: %i\n,N);

fprintf(MC estimate: %f\n,hy(N));

fprintf(standard error: %f\n,eps);

Nmin=ceil(N/100);figure(1);plot([Nmin:N],hy(Nmin:N));

title(MC convergence diagram: probability of full house);

xlabel(No paths);

ylabel(estimate);

4.5.3. Normal random variables. LetXi,i = 1, 2, 3 be independent randomnormal variables with meaniand variance

2i, respectively. Let us compute

the expected value of the difference between the maximum and the minimumof the Xi. In other words, we want to compute E[Y], where

Y= max{X1, X2, X3} min{X1, X2, X3}For this we use the normal random number generator discussed above. A

possible partial matlab code is as follows:

clear;

N=100000;

mu=[0.1 0.3 -0.1];

sigma=[0.1 0.2 0.5];


20/26

20

x=randn(3,N);

z=repmat(mu,1,N)+repmat(sigma,1,N).*x;

y=max(z)-min(z);hy=cumsum(y)./[1:N];

eps=sqrt((sum(y.^2)-N*hy(N)^2)/(N-1))/sqrt(N);

fprintf(No paths: %i\n,N);

fprintf(MC estimate: %f\n,hy(N));

fprintf(standard error: %f\n,eps);

(We have not included the code for plotting the convergence diagram.)

4.5.4. Computing integrals. Our first examples of Monte Carlo simulationsconcerned the computation of multiple integrals and, as special cases, areasand volumes of regions in the plane or space. Let us show how these com-putations fit into the framework of computing expected values of randomvariables.

Let us stick to dimension three for definiteness. Suppose we want tocompute the triple integral

I=

f(x1, x2, x3) dx1dx2dx3,

where is a (bounded) region in the plane and f is a (sufficiently nice)function defined on .

Pick a large boxRof the form

R= [a1, b1]

[a2, b2]

[a3, b3] ={ai

xi

bi, i= 1, 2, 3}

such that is contained in R. Define a trivariate function f on R by

f(x1, x2, x3) =

f(x1, x2, x3) if (x1, x2, x3) 0 if (x1, x2, x3) /

Then we have

I=E[f(X1, X2, X3)],

whereX1,X2,X3are i.i.d. random variables andXi U(ai, bi) is uniformlydistributed on [ai, bi] for i = 1, 2, 3.

This point of view leads to the Monte Carlo simulation algorithm that wediscussed in connection with Example E above.

4.6. Variance reduction techniques. The Monte Carlo method is, gen-erally speaking, easy to implement and does not suffer from the curse of di-mensionality. It does however, have the drawback that it is relatively slow,with rate of convergence (in the probabilistic sense) equal to O(1/

N).

Here we will briefly discuss a couple of methods for speeding up the conver-gence. They will typically still have a convergence rate ofO(1/

N) but the

constant hidden in the O()-terminology may be smaller. More precisely,


21/26

21

the variance of the Monte Carlo estimator is small, explaining the namevariance reduction methods for these techniques.

In these notes we shall cover two variance reduction methods: antitheticvariables and moment matching. There are many more variance reductiontechniques that we have to leave out due to limited space. Some of the moreimportant ones that we have left out are stratified sampling, importancesampling, control variates and low-discrepancy sequences (or quasi-MonteCarlo methods).

4.6.1. Antithetic variables. The method of antithetic variables is based onthe symmetry of the distributions U(0, 1) and N(0, 1):

ifY U(0, 1) then 1 Y U(0, 1). ifY N(0, 1) thenY N(0, 1).

This can be exploited as follows. Suppose we need to generate 2Nsamples

of a random variable Y, and that this is done through 2N samples of auniform random variable Y U(0, 1), say X=f(Y) for some function f.

Normally we would use 2N independent samples

y1, . . . , y2N

ofYand covert these to samples xn= f(yn) ofX.With antithetic variables, we would instead only use Nindependent vari-

ables

y1, . . . , yN

and generate the 2Nsamples ofY as

f(y1), f(1

y1), f(y2), f(1

y2), . . . , f (yN), f(1

yN).

If the underlying random variable Ywas normal rather than uniform,sayY N(0, 1), then we would generate 2Nsamples ofXas

f(y1), f(y1), f(y2), f(y2), . . . , f (yN), f(yN).The method of antithetic variables typically reduces the variance of the

Monte Carlo estimator, but it is hard to quantify how large the reductionwill be. Generally speaking, antithetic variables are good to use unless thereis a specific reason not to.

4.6.2. Moment matching. Moment matching is somewhat similar in spiritto antithetic variables. Suppose we want to sample a random variable X=

f(Y), where Y N((0, 1). Let y1, . . . , yNbe iid samples ofy and think ofthe yns as random variables. Then the random variables= 1N

Nn=1 yn

2 = 1N1N

n=1(yn )2

have expected values

E[] = 0 and E[2] = 1.


22/26

22

Of course, there is no reason why for a particular sample, we would have= 0 and 2 = 1.

The idea is now to correct the numbersy1, . . . , yNto match the first twomoments:

zn= 1

(yn ), n= 1, . . . , N

Now usey1, . . . , yNas the generated samples ofN(0, 1), i.e. set xn = f(yn).

5. Examples from Financial and Actuarial Mathematics

To end these notes we give outline how to apply the Monte Carlo techniqueon three different problems in Financial and Actuarial mathematics.

5.1. Option prices in the Black-Scholes model. Let us consider a (sim-plified) model for the price of a call option on a stock. We have a so called

one-period model with two time points: today (t = 0) and one time unit(e.g. one year) from now (t= 1). The stock priceS0 today is known. Thestock price att = 1 is a random variable. We assume that it is a log-normalrandom variable. This means that

S1=S0eY,

whereY N(r 2/2,2) is a normal random variable with meanr 2/2and variance. Herer is the interest rate and is thevolatilityof the stock.

We are interested in computing the price of a call option expiring at timet = 1 and with strike price K. Such an option gives the holder (i.e. theowner) the right, but not the obligation, to buy one stock at time t = 1 forthe predetermined price K. At time t = 1, the value of this option will be

C1 = max{0, S1 K}=

S1 K ifS1> K0 otherwise.

One can then argue that the market price of the option at time t = 0 (today)should be

C0= er[C1] =er max{0, S1 K}.

To implement this model we could do as follows:

Fix a large number N. Generate N independent samples z1, . . . , zNof an standard normal

random variable Z N(0, 1). Setyn = zj+ (r

2/2). Theny1, . . . , yNare independent samples

of a random variable Y N(r 2/2,). Set s1,n =s0e

yn . Then s1,1, . . . , s1,Nare independent samples fromthe distribution ofS1.

Set c1,n = max{0, s1,n K}. These are the option values at timet= 1.

Set cn = er 1

n

nk=1 c1,k. These are the successive Monte Carlo

estimates of the option price at time t= 0.


23/26

23

ReportN, the final MC estimate cN, the standard errorNand plota convergence diagram.

The reader should note that there is an exact formula, the celebratedBlack-Scholes formulafor the option price C0. In practice, one would there-fore not use Monte Carlo simulations to compute the option prices. However,it is well known that the Black-Scholes model does not accurately reflect allthe features of typical stock prices. In more advanced (and realistic) models,it may not be possible to produce an exact formula for the option price andthen numerical methods are crucial. Monte Carlo simulations is one suchmethod that is easy to implement.

5.2. The collective loss model. Next we study an example from ActuarialMathematics. Pretend that you are an insurance company insuring, say,cars. There are a lot of policy holders (car owners) and each time a policy

holder has an accident, loses his/her car to theft or fire etc, he/she makes aclaim and you have to pay some money. In compensation for this, you receivemoney in the forms of insurance premia. How big premia is it reasonable tocollect?

A general way to look at this problem is to say that the total lossS, i.e.the total amount of money that you (the insurer) must pay out during, say, ayear, is a random variable. Definitely, you would like to set the premium sothat the total premium collected is at least the expected valueE[S], or elseyou be in the red on average. In practice, premia are higher than that. Letus for simplicity say that the total premium collected is twicethe expectedvalue E[S]. Let us ask what the probability is that the collected premiaexceeds the total (oraggregate) loss.

To attack this problem, we must know the distribution of the aggregateloss S(during a fixed time period, say a year). A common model for S isthe so called collective risk model, in which

S= X1+ X2+ + XN,

whereNis the total number of accidents (the frequency) andXj is what youhave to pay to the insured in accident j. Here both N and Xj are randomvariables. We assume that:

the random variables X1, X2,. . . are i.i.d. ., so they have the samedistribution as a random variable Xcalled the severity;

N and X1,X2,. . . are all independent.

One can then show that

E[S] =E[N]E[X],

and we are interested in computing the probability

P{S >2E[S]}.

To do this we must know the distributions ofN and S above, but even ifwe know them it can be difficult to calculate the probability above exactly.This is where Monte Carlo simulations come into the picture.


24/26

24

Let us consider the case when frequency Nis a Poisson random variableand severity X 0 is a Pareto random variable. This means thatN takesvalues 0, 1, 2, . . . and has a pmf of the form

fN(k) =ek/k!,

for some > 0 and that Xis a continuous random variable with pdf

fX(x) =

(x + ),

for some constants > 0 and > 1. Then we have

E[S] =E[N]E[X] =

1 .We can generate samples of the Poisson distribution using the matlab

command poissrnd in the statistics toolbox. It can also be done by hand,but that is a little tricky.

As for the Pareto distribution, it can be sampled by inverting the cdf.

y= FX(x) = 1

(x + ) x=

1

(1 y)1/ 1

.

To obtain a Monte Carlo estimate for the probability P{S >2E[S]} weproceed as follows:

Pick a large number M(we do not call it N this time...) Generate M independent samples n1, n2, . . . , nMfrom the Poisson

distribution. Generate n1+ + nm independent samples

xmj, 1 m M, 1 j nmfrom the Pareto distribution.

Setxm= xm1+ + xmnm for 1 m M. Setzm= 1 ifxm> 2 E[S] =

1 andzm= 0 otherwise. Set zm=

1

m

mj=1 for 1 m M.

Report M, the final Monte Carlo estimate zM, the standard errorMand draw a convergence diagram.

5.3. The binomial tree model. Finally we consider another model forthe evolution of stock prices: the binomial tree model. In this model, timestarts att = 0 and ends at some timet = Tin the future. In between, thereare Mtime steps, so time is discrete,

t= 0,t, 2t , . . . , M t= T .

We consider a market with a single stock. The price of the stock at timemtis denoted by Sm. It is a random variable. Given its value Sm at timemt, its value Sm+1 is equal to

Sm+1=

Sm u with probabilitypu

Sm d with probabilitypd


25/26

25

Here u and d are constants with d 0, pu+pd= 1.This implies that given the value S0 of the stock at time 0, Sm can take

the m + 1 possible valuesS0d

m, S0udm1, . . . , S 0um.

It is possible to transform Sm to a binomial randdom variable.

Now, we wish to do two things with this model:

Calibrate the parameters u, d, pu and pd. Compute the price of a so called Barrier option in this model.

5.3.1. Calibration. We wish to use the binomial tree to model the behaviorof a real stock. How should we pick the parametersu,d,puand pd. Basicallywe need four conditions to nail down these four unknown. One condition is

(1) pu+pd= 1,so we need three more.

We can get two conditions by specifying the conditional expected valueand variance ofSm+1 giveen Sm. For financial reasons it is natural to de-mand

E[Sm+1|Sm] =Sm(1 + rt)

Var[Sm+1|Sm] =S2m

2t,

where r is the interest rate and is the volatilityof the stock (assumedconstant). This leads to the equations

puu +pdd= 1 + rt(2)

pu(u (1 + rt))2 +pd(1 + rt d)2 =2t(3)Finally, it is common to put in a condition of symmetry, namely

(4) ud= 1

It is in general not possible to find an analytic expression for the valuesofu, d, pu, pd determined by the four conditions (1)-(4). However, we canreduce the four equations to one, and solve the latter using e.g. Newtonsmethod!

5.3.2. Pricing a barrier option. We wish to use Monte Carlo simulations toprice an up-and-out barrier call option in the binomial tree model. Suchan option works in the same way as a (vanilla) call option, but it becomesworthless if the stock price exceeds a certain barrier B prior to expiry.

Today is t= 0 and the option expires at time t= T. At the latter time,the option is worth

CM=

0 ifSM K0 ifSm B for some m, 0< m MSM K otherwise


26/26

26

Here SM and Sm, 1 m Mare of course unknown at time t= 0.One can then argue that at time 0, the option is worth

C0= (1 + rt)ME[max{SM K, 0} 1{max1mMSm

MonteCarlo Intro in Mtlab

Documents

Transcript of MonteCarlo Intro in Mtlab