1 Simulating Brownian motion (BM) and …

1 Copyrightc 2013 by Karl Sigman1 Simulating Brownian motion (BM) and geometric Brownianmotion (GBM)For an introduction to how one can construct BM, see the Appendix at the end ofthese stochastic processB={B(t) :t 0}possessing (wp1) continuous sample paths is calledstandard Brownian motion (BM) (0) = both stationary and independent (t) B(s) has a normal distribution with mean 0 and variancet s,0 s < )and3)together can be summarized by: Ift0= 0< t1< t2< < tk, then the incrementrvsB(ti) B(ti 1), i {1,..k}, are independent withB(ti) B(ti 1) N(0,ti ti 1)(normal with mean 0 and varianceti ti 1). In particular,B(ti) B(ti 1) is independent ofB(ti 1) =B(ti 1) B(0).If we only wish to simulateB(t) at one fixed valuet, then we need only generate a unitnormalZ N(0,1) and setB(t) = tZ. But typically, we will want to simulate (say)kvaluesat timest0= 0< t1< t2< < tkto get the entire vector (with correlated coordinates):(B(t1).

2 ,B(tk)).We can easily do so by using the fact thatB(s+t) =B(s) + (B(s+t) B(s)), andB(s)andB(s+t) B(s) are independent we generatekiid unit normalsZ1,Z2,..,Zk, then we can construct the independentincrements viaB(ti) B(ti 1) = ti ti 1Zi, i= 1,.., to simulate the valuesB(t1),..,B(tk), we can use the recursionB(ti+1) =B(ti) + (B(ti+1) B(ti)) =B(ti) + ti+1 tiZi+1, i {0,..k 1}.This yields: Simulating Standard BM at times0 =t0< t1< t2< < tk:Sequentially generate unit normalsZ1,Z2,..,Zk, and recursively defineB(t1) = t1Z1B(t2) =B(t1) + t2 t1Z2= t1Z1+ t2 (tk) =k i=1 ti ti the end, to simulate BM then, we need only generate unit normals, which we alreadylearned how to do using (for example) the polar method or the acceptance rejection BM with driftX(t) = B(t) + twill denote the BM with drift Rand variance term >0. It hascontinuous sample paths and is defined (0) = both stationary and independent (t) X(s) has a normal distribution with mean (t s) and variance 2(t s),0 s < (t) X(s) thus can be constructed (simulated) by generating a standard normal rvZand settingX(t) X(s) = t sZ+ (t s).

3 Again, by the stationary and independentincrements, we can simulate such a BM at times 0 =t0< t1< t2< < tk, by generatingkiid unit normalsZ1,Z2,..,Zkand using the recursionX(ti+1) =X(ti) + (X(ti+1) X(ti)) =X(ti) + ti+1 tiZi+1+ (ti+1 ti). Simulating BM with drift and variance term at times0 =t0< t1< t2< < tk:Sequentially generate unit normalsZ1,Z2,..,Zk, and recursively defineX(t1) = t1Z1+ t1X(t2) =X(t1) + t2 t1Z2+ (t2 t1) = t1Z1+ t1+ t2 t1Z2+ (t2 t1)..X(tk) =k i=1( ti ti 1Zi+ (ti ti 1)). Geometric BMGeometric Brownian motion (GBM) is given byS(t) =S(0)eX(t), t 0,whereX(t) = B(t) + t, t 0,is a (t)has a lognormal distribution for each fixedt >0. In general ifY=eXis lognormal withX N( , 2), then we can easily simulateYviasettingY=e Z+ , withZ N(0,1).Moreover, for any 0 s < tit holds thatS(t) =S(0)S(s)S(0) S(t)S(s)=S(0)eX(s) eX(t) X(s),and since the incrementX(s) is independent of the incrementX(t) X(s), we conclude thatthe consecutive ratiosS(s)S(0)andS(t)S(s)are independent lognormals.

4 We can thus simulate the pair(S(s),S(t)) by generating two iidN(0,1) rvs,Z1,Z2and settingS(s) =S(0)e sZ1+ s, S(t) =S(s)e t sZ2+ (t s)=S(0)e sZ1+ s e t sZ2+ (t s).More generally, for 0 =t0< t1< t2< < tk, defineYi=S(ti)/S(ti 1), i {1,2,..,k}.Then we can write2S(t1) =S(0)Y1S(t2) =S(t1)Y2=S0Y1 (tk) =S(tk 1)Yk=S0Y1 Y2 independent lognormal rvs and can be constructed by generatingkiidN(0,1) rvs,Z1,Z2,..Zkand settingYi=e ti ti 1Zi+ (ti ti 1), i {1,2,..,k}.(1) Simulating paths of GBM is thus an easy consequence of our algorithm for Simulating BMsince for 0 =t0< t1< t2< < tk, the following recursion holdsS(ti+1) =S(ti)eX(ti+1) X(ti), i {0,1,..,k 1}. Simulating Geometric BM (with drift and variance term ) at times0 =t0< t1< t2< <tk:Sequentially generate unit normalsZ1,Z2,..,Zk, and set theYias in (1). Then recursivelydefineS(t1) =S(0)Y1S(t2) =S(t1)Y2=S0Y1 (tk) =S(tk 1)Yk=S0Y1 Y2 Applications in Financial EngineeringHere we letS(0) =S0denote the price per share of a risky asset (stock) initially, andS(t) =S0eX(t)as the price at timet.

5 The classic European call option with expiration dateTand strikepriceKhas payoff at timeTofCT= (S(T) K)+. The famous Black-Scholes-Merton optionpricing theory/formula makes this option s price known explicitly, but other options(derivativesof the stock) are typically impossible to compute exactly; Monte Carlo simulation is thuscommonly used to do estimate the call optionA variation on a European call option (that is cheaper) is to average the price of the stock overthe time interval [0,T] and use that instead ofS(T) yielding payoff(1T T0S(t)dt K)+.In practice one can t compute (or simulate exactly) such an average as1T T0S(t)dtto offer suchan option, and instead one samples the price at a sequence of times such as the beginning of3each day. So let s assume that time is in days, thatTis an integer and thus we will considerthe payoffCT=(1TT n=1S(n) K)+,(2)and our objective is to estimate the expected payoff,E(CT).

6 We thus need to simulate iidcopies ofCTand then theYirvs from the previous Section, they are now iid sinceti ti 1= 1, and thuseachX(ti) X(ti 1)has aN( , 2) the recursionS(n+ 1) =S(n)eX(n+1) X(n), n {0,..T 1},we conclude that we can represent the stock prices by introducingTiid unit normalsZ1,..,ZT,and rewritingS(n+ 1) =S(n)e Zn+1+ .We thus can construct a copy call optionAnother variation on a European call option (that is cheaper) is the introduction of a upperbarrierb >0 which, in order that the holder receive payoff (S(T) K)+at timeT, the stockprice must remain below the barrier at pre-specified times 0< t1< t2< < tk< T. Thusthe payoff isCT= (S(T) K)+I{S(t1)< b,..,S(tk)< b};ifS(ti) bat any one of thekcheck timest1,..,tk, thenCT= (ti+1) =S(ti)eX(ti+1) X(ti), and noting thatX(ti+1) X(ti) N( (ti+1 ti), 2(ti+1 ti)), we can construct sequentially theS(ti) and check ifS(ti) b; if so then we stop and setCT= 0; otherwise we continue to the next check point.

7 If allkcheck points are passed with-out violating the barrier constraint, then we finally constructS(T) =S(tk)eX(T) X(tk)withX(T) X(tk) N( (T tk), 2(T tk)), and setCT= (S(T) K)+. This then gives us ourfirst sample illustrate this below:Algorithm for generating one copy ofCT:Inputkand thektimes 0< t1< t2< < tk< , S=S0, i= GenerateZ. ResetS=Se tZ+ IfS bandi k, then setCT= 0 and Otherwise ifS < bandi k, reseti=i+ 1,t=ti ti 1and go back to(1).4. Otherwise ifi=k+ 1, then resett=T tk, generateZ, resetS=Se tZ+ t,and setCT= (S K)+and Other Applications in FinanceMonte Carlo simulation can also be used to estimate other quantities of interest in finance thatdo not involve derivatives. For example, suppose you invest in two different stocks,S1(t) andS2(t), buyingN1shares of the first andN2of the second. At timet= 0 you payV(0) =N1S1(0) +N2S2(0) for this portfolio, and at timetit will be worthV(t) =N1S1(t) +N2S2(t).

8 Although each stock price on its own has a lognormal distribution, the sum of the two does not;computations involving this sum can be intractable. The stocks are typically correlated you wish to compute the probability that at timet=T, the value of your investmenthas increased by at least 10%. That is, you wish to computeP(V(T) ( )V(0)).This is equivalent to taking expected values of the indicator functionX=I{N1S1(T) +N2S2(T) ( )V(0)}.For illustration, let s suppose the two stocks are independent in which case by lettingZ1andZ2denote two generated iid unit normals, we can constructXviaX=I{N1S1(0)e 1 TZ1+ 1T+N2S2(0)e 2 TZ2+ 2T ( )V(0)}.In general, some correlation exists between the two stocks, and this can be incorporated bygenerating the unit normals with a desired correlation:One can constructcorrelatedunit normalsZ1andZ2with any desired correlationcoefficient 1< <1 by first generatingY1andY2iid unit ones and settingZ1=Y1andZ2= Y1+ 1 , simulatingniid copiesX1,X2.

9 ,XnofXyields our estimatorP(V(T) ( )V(0)) 1nn i= : Taking two independent standard Brownian motions,W1(t), W2(t) we can constructa correlated two-dimensional Brownian motion via definingB1(t) =W1(t), B2(t) = W1(t) + 1 2W2(t). APPENDIX: Construction of Brownian motion from the simple sym-metric random walkRecall the simplesymmetricrandom walk,R0= 0,Rn= 1+ + n=n i=1 i, n 1,where the iare iid withP( = 1) =P( = 1) = thusE( ) = 0 andV ar( ) =E( 2) = view timenin minutes, andRnas the position at timenof a particle, moving on IR,which every minute takes a step, of size 1, equally likley to be forwards or backwards. BecauseE( ) = 0 andV ar( ) = 1, it follows thatE(Rn) = 0 andV ar(Rn) =n, n a large integerk >1, if we instead make the particle take a step every 1/kminutesand make the step size 1/ k, then by timetthe particle will have taken a large number,n=tk,of steps and its position will beBk(t) =1 ktk i=1 i.

10 (3)(By convention iftkis not an integer then we replace it by the largest integer less than orequal to it; denoted by [tk].) This leads to the particle taking many many iid steps, but eachof small magnitude, in any given interval of time. We expect that ask , these small stepsbecome a continuum and the process{Bk(t) :t 0}should converge to a process{B(t) :t 0}with continuous sample paths. We call this process Brownian motion (BM) after the Scottishbotanist Robert properties will be derived that for fixedk, any incrementBk(t) Bk(s) =1 ktk i=sk+1 i,0 s < t,has a distribution that only depends on the length,t s, of the time interval (s,t] because itonly depends on the number,k(t s), of iid imaking up its construction. Thus we deducethat the limiting process (ask ) will possessstationary increments:The distribution ofany incrementB(t) B(s)has a distribution that only depends on the length of the time intervalt s.)