stat.CO 2012-11-27 0

epsilon-Strong simulation of the Brownian path

We present an iterative sampling method which delivers upper and lower bounding processes for the Brownian path. We develop such processes with particular emphasis on being able to unbiasedly simulate them on a personal computer. The dominating proce…

Authors: Alex, ros Beskos, Stefano Peluchetti

epsilon-Strong simulation of the Brownian path
Bernoul li 18 (4), 2012, 1223–12 4 8 DOI: 10.315 0/11-BEJ 383 ε -Strong sim ulation of the Bro wnian path ALEXANDR OS BESKOS 1 , STEF ANO PELUCHETTI 2 and GARETH R OBE R TS 3 1 Dep artment of St atistic al Scienc e, UCL, Gower St r e et, L ondon, WC1E 6BT, UK. E-mail: alex@stats.ucl.ac.uk 2 HSBC Bank, 8 Canada Squar e, L ondon, E14 5HQ, UK. E-mai l: stefano.p eluchett i@hsb cib.c om 3 Dep artment of Statistics, University of Warwick, Coventry, CV4 7AL, U K. E-mail: gar eth.o.r ob erts@warwick.ac.uk W e present an iterativ e sampling method whic h deliv ers upper and low er b ounding p rocesses for th e Brow nian path. W e develop s uch pro cesses with particular emphasis o n b eing able to unbiasedly simulate th em on a p ersonal computer. The dominating pro cesses converge almost surely in the supremum and L 1 norms. In particular, the rate of conv erge in L 1 is of the order O ( K − 1 / 2 ), K denoting the compu t ing cost. The a.s. enfolding of th e Brownian path can b e exploited in Monte Carlo applications inv olving Brownian paths whence our algorithm (termed the ε -strong al gorithm) can deliver un biase d Mon te Carlo estimators ov er path expectations, o vercoming discretisation errors characterising stand ard approac hes. W e will show analytical results from applications of the ε -strong algorithm for estimating exp ectations arising in option pricing. W e will also illustrate that individual steps of th e al gorithm can b e of separate interest, giving new simulation metho ds for interesting Brow nian distributions. Keywor ds: Brownia n bridge; inters ection la yer; iterative algorithm; option p ricing; pathwise conv ergence; unbiased sampling 1. In tro duction Brownian motion (BM) is a n ob ject of pa r amoun t significance in sto c hastic mo delling. Starting from its original mathematical f or m ulation by [ 2 ], its pr operties ar e still under meticulous in vestigation by contempor ary resea rc hers. Relev ant to the purp oses of this pap er, considerable w or k has fo cused on v arious constr uctions and repre sen tations of BM paths. Leaving aside the simple finite-dimensional Ga ussian s tr ucture of BM, resea r c hers hav e o ften b een interested on mor e complex functiona ls. Hitting times, extremes, lo cal times, reflections and other characteristics of BM hav e b een investigated (for a genera l exp osition see [ 18 ]). F or simulation purp oses, man y of the relev ant distributions are easy to sample from on a computer [ 10 ]. Several c onditione d constructions o f B M are also known rela ting B M with the Bessel pro c e ss, the Rayleigh dis tribution and other sto c hastic ob jects (see , e.g., [ 3 ]). This is an electronic reprint of the original article published by the I SI/BS in Bernoul li , 2012, V ol. 18, No. 4, 1223–1248 . This reprint differs from the original in pagination and typographic detail. 1350-7265 c  2012 ISI/BS 2 A. Beskos, S. Peluchetti and G. R ob erts This pap er presents a co n tribution o f our own at simulation metho ds for Brownian dynamics. W e develop an itera tiv e sampling a lgorithm, the ε -str ong algo rithm, which simulates upper and low er paths env eloping a .s . the Brownian path. T o meet this ob- jective, we co lle c t a num be r of characterisations a nd combine them in a way that they can deliver s imple sampling methods implementable on a p ersonal computer. W e will show that after O ( K )-computational e ffort, the dominating pro cess hav e L 1 -distance of O ( K − 1 / 2 ). This a.s. enfolding of the Brownian path ca n b e e x ploited in Monte Carlo ap- plications inv olving Brownian motion in tegr als, minima, maxima or hitting times; in such scenaria, the ε -strong alg orithm ca n deliver unbiased Monte Carlo estimators ov er Brow- nian exp ectations, ov erco ming discretizatio n e rrors characterising standard appro ac hes (for the latter a ppr oac hes, see, for instanc e , the ex position in [ 12 ] in the context of ap- plications in finance). W e will show applications of the alg orithm and exp erimen tally compare the r equired computing r esources against t ypica l alternatives employ ed in the literature inv olving Euler approximation. Our exa mples will inv olve a collection of double-barrie r option pricing problems in a Black and Sch ole s framework arising in fina nce. Also, we will demonstrate that individual steps of the alg orithm can b e o f s e parate interest, g iving new simulation metho ds for interesting Brownian distr ibutions. The ε -strong algo r ithm delivers a pa ir of dominating pro cesses, denoted by X ↓ ( n ) = { X ↓ u ( n ); u ∈ [0 , 1 ] } and X ↑ ( n ) = { X ↑ u ( n ); u ∈ [0 , 1 ] } , that ca n b e simulated on a p ersonal computer without any discretisa tion erro r, with the prop ert y: X ↓ u ( n ) ≤ X ↓ u ( n + 1) ≤ X u ≤ X ↑ u ( n + 1) ≤ X ↑ u ( n ) (1.1) for all instances u ∈ [0 , 1] ; here, X is the Brownian path. The tw o dominating pro cesses will co n verge in the limit: w . p . 1 , lim n →∞ sup u ∈ [0 , 1] | X ↑ u ( n ) − X ↓ u ( n ) | → 0 . (1.2) The algor ithm builds o n the no tion of the interse ction layer , a collective informa tion, containing the s tarting and ending p oin ts of a Br o wnian pa th together with information ab out its extrema. A num b er of op erations (bisection, refinemen t, see main text) can be a pplied on this infor mation, explicitly on a c o mputer, allowing the s ampler to iterate itself to get closer to X . W e sho uld no te here that the metho ds desc r ibed in this pap er will b e relev ant a lso for nonlinear Sto c has tic Differential Equations (SDEs ). Recent developmen ts in the simula- tion of SDEs under the framework of the so-calle d ‘Exa ct Algorithm’ (see [ 4 – 7 , 9 , 1 3 ]) build up on the result that, conditionally o n a collection of randomly sampled po in ts, the path o f the SDE is made of indep enden t Brownian paths. Onc e this colle c tion of p o in ts is sampled, the metho dology of this pap er can then be applied separ ately on each of the constituent Brownian sub-paths. The str ucture o f the pap er is a s following. In Sectio n 2 , we pr esen t the notion of the interse ction layer which will b e critical for our metho ds. In Section 3 , we pres e n t the individual s teps forming the ε -strong alg orithm; they will require origina l simulation tech- niques for some Br o wnian distributions. Once we identify in Section 4 the ζ - function , a n ε -St r ong simulation of the Br ownian p ath 3 alternating mono tone ser ies at the co r e of Brownian dynamics, we exploit its structur e in Section 5 to a nalytically develop these new s a mpling metho ds. In Section 6 , we apply the ε -str ong alg orithm to unbiasedly estimate some path ex pectations aris ing when pricing options in finance. W e will contrast the computationa l cost of the algo rithm with Euler approximation a lter nativ es to g et a b etter understanding of its pra ctical co mp etitiveness. In Sectio n 7 , we sketc h some other p otential applications of the ε -s trong algo r ithm. W e finish with some discussion a nd conclusions in Section 8 . 2. In tersection la y er and op erations W e will, in general, write paths as X = { X u ; u ∈ [ s, t ] } for s < t . A Brownian bridge on [ s, t ] is a B r o wnian motion conditio ne d to start a t X s and end at X t , for some presp ecified X s , X t ; its finite-dimensio nal dynamics a re easily deriv able following this interpretation (see, for instance, [ 18 ]). Instrumental in our considera tions is the notion of (what we call) the interse ction layer . Consider a Brownian br idge X on [ s, t ] . L e t m s,t , M s,t be the extrema of X : m s,t = inf { X u ; u ∈ [ s , t ] } , M s,t = sup { X u ; u ∈ [ s , t ] } . The ε -strong algor ithm will requir e some information on b oth m s,t and M s,t . W e will ident ify interv als: U s,t = [ U ↓ s,t , U ↑ s,t ] , L s,t = [ L ↓ s,t , L ↑ s,t ] , such that: M s,t ∈ U s,t , m s,t ∈ L s,t . W e will write simply m , M , U ↑ , U ↓ , L ↑ , L ↓ ignoring the s, t -subscripted versions when the time interv al under consideratio n is clea rly implied by the context. The intersection lay er idea refers to the collective information I s,t = { X s , X t , L s,t , U s,t } , (2.1) that is the s tarting a nd ending p oints o f the br idge together with interv als that contain its maximum and minimum . Figur e 1 (a) pre s en ts a gra phical illustratio n of the intersection lay er: the extrema of an under lying Brownian bridg e lie in the s haded recta ngles. W e will lo ok now a t tw o simple op erations on the informatio n I s,t which nonetheless will b e the building blo cks of the complete ε - strong a lgorithm des c ribed in the next section. 2.1. Refining the information I s,t During the iter ations at the ex ecution of the ε -strong algor ithm, for each piece of infor - mation I s,t we will need to control the width of the lay ers L s,t , U s,t relatively to the size t − s of the time interv al to ensur e convergence of the b ounding paths env eloping the 4 A. Beskos, S. Peluchetti and G. R ob erts Figure 1. T op panel (a): the intersecti on lay er information I s,t for a Brownian path . The underlying tra jectory starts at X s and finishes at X t with its ex trema found in th e shaded areas. Bottom panel (b): the bisection of I s,t into I s,t ∗ and I t ∗ ,t . The algorithm simulates X t ∗ and then decides that the ext rema for each of the interv als [ s, t ∗ ] and [ t ∗ , t ] are in th e shaded areas, that is, U s,t ∗ = [ X t ∗ , U ↓ ], U t ∗ ,t = [ U ↓ , U ↑ ], L s,t ∗ = [ L ↓ , L ↑ ] and L t ∗ ,t = [ L ↑ , X t ∗ ]. The algorithm output s the up gra ded information I s,t ∗ = { X s , X t ∗ , L s,t ∗ , U s,t ∗ } and I t ∗ ,t = { X t ∗ , X t , L t ∗ ,t , U t ∗ ,t } . ε -St r ong simulation of the Br ownian p ath 5 T able 1. The pro cedure for b ise cting the information I s,t . It retu rns the intersection laye rs I s,t ∗ and I t ∗ ,t with refin ed information ab out th e u n derlying p ath (compared to I s,t ) Bise ct( I s,t ): 1. Set t ∗ = ( t + s ) / 2. Simulate X t ∗ giv en I s,t . Set U ↓ = U ↓ ∨ X t ∗ , L ↑ = L ↑ ∧ X t ∗ . 2a. D ecide if U s,t ∗ = [ X s ∨ X t ∗ , U ↓ ] or [ U ↓ , U ↑ ]. 2b. D ecide if U t ∗ ,t = [ X t ∗ ∨ X t , U ↓ ] or [ U ↓ , U ↑ ]. 2c. Decide if L s,t ∗ = [ L ↓ , L ↑ ] or [ L ↑ , X s ∧ X t ∗ ]. 2d. D ecide if L t ∗ ,t = [ L ↓ , L ↑ ] or [ L ↑ , X t ∗ ∧ X t ]. 3. Retu rn I s,t ∗ ∨ I t ∗ ,t . underlying Br o wnian pa th. Thus, the refinement of the information I s,t corres p onds to a pr ocedure that up dates I s,t by ha lving the allowed width for the minimum m or the maximum M of the path, thereby corr espondingly up dating the layers L s,t or U s,t . More analy tically , refinement of I s,t corres p onds to deciding whether the minimum m on [ s, t ], alrea dy known to b e in [ L ↓ , L ↑ ], lies in [ L ↓ , ( L ↓ + L ↑ ) / 2] o r [( L ↓ + L ↑ ) / 2 , L ↑ ], that is, whether L s,t is equal to [ L ↓ , ( L ↓ + L ↑ ) / 2] o r [( L ↓ + L ↑ ) / 2 , L ↑ ]; the appa ren t a na logue of such a co nsideration applies for the ma xim um M . The ana lytical metho d of sampling the relev ant bina ry ra ndom v ariables for carry ing out this pro cedure will b e desc ribed in Section 5 . 2.2. Bisecting the information I s,t This is a more inv olved oper a tion on I s,t , and inv olves bisecting I s,t int o the more analytical infor mation I s,t ∗ ∨ I t ∗ ,t for some intermediate time instance t ∗ ∈ ( t, s ) . In particular, we will b e selecting t ∗ = ( t + s ) / 2 within the ε -strong algorithm. The metho d beg ins by sampling the middle p oin t X t ∗ conditionally o n I s,t , and then appr opriately sampling the lay ers for the tw o pieces of informa tio n I s,t ∗ , I t ∗ ,t . The pra cticalities of implemen ting the sec o nd par t of the metho d will dep end o n whether X t ∗ falls within a lay er of I s,t or no t, thus we present the bisection op eration in more detail in T able 1 . Note that if X t ∗ > U ↓ the t wo upp er lay ers (for I s,t ∗ and I t ∗ ,t ) w ill b e dir ectly set to [ X t ∗ , U ↑ ], and we will have to simulate extr a r andomness a bout the underly ing pa th only to determine the lo wer layers. Co rrespo ndingly , if X t ∗ < L ↑ the tw o low er lay ers will immediately b e set to [ L ↓ , X t ∗ ]. In the scena rio when L ↓ < X t ∗ < U ↓ , we will hav e to simu late extra ra ndomness to determine a ll four layers. W e describ e in Section 5 the algorithms for sampling X t ∗ and determining the lay ers. Figure 1 shows a graphica l illustration o f the bisection pro cedure. 3. ε -Strong sim u lat ion of Bro wnian path W e introduce an iter a tiv e simulation algor ithm with input a Brownian bridg e X o n the domain [0 , 1] a nd output, after n iterations , upp e r and lower dominating pr ocesses X ↓ ( n ) = { X ↓ u ( n ); u ∈ [0 , 1 ] } a nd X ↑ ( n ) = { X ↑ u ( n ); u ∈ [0 , 1 ] } satisfying the monoto nicit y 6 A. Beskos, S. Peluchetti and G. R ob erts T able 2. The ε -strong algorithm. It iterativ ely unv eils extra information ab out the un d erlying path. It outputs the collection of intersection la yers P = W 2 n j =1 I ( j − 1)2 − n ,j 2 − n ε -str ong( X 0 , X 1 , n ): 1. Initialize U 0 , 1 , L 0 , 1 , set I 0 , 1 = { X 0 , X 1 , U 0 , 1 , L 0 , 1 } . Set P = {I 0 , 1 } and i = 1. 2. F or each of the 2 i − 1 inters ection lay ers in P , say I s,t , d o th e follo wing: i. Bisect the information I s,t into I s,t ∗ , I t ∗ ,t , where t ∗ = ( t + s ) / 2. ii. Refine I s,t ∗ , I t ∗ ,t until the width of their la yers is not greater than p ( t − s ) / 2. 3. Collect th e up dated information, P = W 2 i j =1 I ( j − 1)2 − i ,j 2 − i . 4. If i < n set i = i + 1 and return to Step 2; otherwise return P . and limiting r equiremen ts ( 1.1 ) a nd ( 1.2 ) r espectively . Note that X here is a contin uous time Br o wnian br idge path, thus a n infinite-dimensional rando m v ariable. How ever, the bo unding pro cesses will b e piece-wise constant, thus inherently finite-dimensional. O ne will b e able to r e alise complete sa mple paths o f X ↓ ( n ) o r X ↑ ( n ) o n a computer without retreating to any sort of discr etization or approximation err ors (apart from those due to finite co mputin g accura cy). 3.1. ε -Strong algorithm Given so me initial intersection lay er info r mation I 0 , 1 , the algorithm will na turally set X ↑ u (0) = U ↑ 0 , 1 and X ↓ u (0) = L ↓ 0 , 1 for all instances u ∈ [0 , 1 ]. It will then iter a tiv ely bise ct the a cquired int er section lay ers, as descr ibed in Section 2.2 , to o bta in mor e infor mation ab out the underlying sample pa th on finer time interv als. T o ens ur e co n vergence of the discrepancy X ↑ ( n ) − X ↓ ( n ) the algor ithm will sometimes r efine the infor mation on s ome int ers ection layers, as descr ibed in Sectio n 2.1 , to reduce the uncertaint y for the extr e ma. W e give the ps e udo co de ab out the algor ithm in T a ble 2 . Utilising the information the ε - strong algorithm returns, w e define the dominating pro cesses as fo llows: X ↑ u ( n ) = 2 n X i =1 U ↑ ( i − 1)2 − n ,i 2 − n · I u ∈ (( i − 1)2 − n ,i 2 − n ] , (3.1) X ↓ u ( n ) = 2 n X i =1 L ↓ ( i − 1)2 − n ,i 2 − n · I u ∈ (( i − 1)2 − n ,i 2 − n ] . The squar e-ro ot rate at Step 2 .ii of the alg orithm in T able 2 is to guar a n tee conv erge nc e of the dominating paths with minimal computing cost: it pr ovides the correct distribu- tion of e ffort b etw een time-interv al a nd extr ema-in terv al bisec tio ns. T o understand this, note that the ra nge of a Brownian motion (or a Brownian bridge) o n [0 , 2 − n ] sca les as O (2 − n/ 2 ); s ee, for instance, [ 18 ]. Thus, had we used the actual Brownian minima and maxima to define do minating pr o cesses for the Br o wnian path in the way of ( 3.1 ) the rate of co n vergence would hav e b een O (2 − n/ 2 ); we cannot exceed such a r ate, but we ε -St r ong simulation of the Br ownian p ath 7 can pres erv e it if our extrema are not further than O (2 − n/ 2 ) from the actual ones. This int uitive statement will b e made rig orous in the sequel, when an explicit r esult on the rate o f convergence of the do minating pr ocess e s in L 1 -norm is g iv en. Figure 2 shows succ e s siv e steps of the ε -strong a lgorithm as implement ed on a com- puter. F or each n , the horizontal black lines show the interv al wher e the maxima and the minima a re lo cated: this infor mation is av ailable for all 2 n sub-interv als bisecting the initial time interv al [0 , 1] . The dashed bla c k line corr esponds to the linear interpola tion of successively unv eiled p ositions of the underlying Brownian path. The last graph (f ) corres p onds to n = 12; in this cas e, we hav e zo omed on a par ticular subinterv al of [0 , 1 ] to b e a ble to visualise the difference betw een the b ounding pa ths a nd the underly ing Brownian o ne . 3.2. Con v ergence prop erties Almost sure co n vergence o f the dominating paths follows dir ectly from the contin uity of the Br o wnian path X . The analytica l pro of is given in the following prop osition. Prop osition 3.1. Consider the c ontinuous-time pr o c esses X ↑ ( n ) , X ↓ ( n ) define d in ( 3.1 ). Then, the c onver genc e in supr emum n orm in ( 1.2 ) wil l hold in the limit n → ∞ . Pro of. F or a B ro wnian bridg e X o n [0 , 1] , we consider : D n := sup 1 ≤ i ≤ 2 n ( M ( i − 1)2 − n ,i 2 − n − m ( i − 1)2 − n ,i 2 − n ) . Uniform contin uity implies that, with probability 1 : lim n →∞ D n = 0 . Now, we hav e that: sup u ∈ [0 , 1] | X ↑ u ( n ) − X ↓ u ( n ) | ≤ D n + 2 · 2 − n/ 2 → 0 , where we hav e used the fact that Step 2.ii of the ε -strong algo rithm guar a n tees that U ↑ ( i − 1)2 − n ,i 2 − n ≤ M ( i − 1)2 − n ,i 2 − n + 2 − n/ 2 , L ↓ ( i − 1)2 − n ,i 2 − n ≥ m ( i − 1)2 − n ,i 2 − n − 2 − n/ 2 .  A mor e inv olved result ca n give the r ate of conv er g ence o f the dominating pro cesses and will b e of practical significa nce for the efficiency of Monte Car lo metho ds based on the ε - s trong algo rithm. Prop osition 3.2. Consider the L 1 -distanc e: | X ↑ ( n ) − X ↓ ( n ) | 1 = Z 1 0 | X ↑ u ( n ) − X ↓ u ( n ) | d u. 8 A. Beskos, S. Peluchetti and G. R ob erts Figure 2. The ε -strong algorithm as applied on a p ersonal computer. F or each step n , the horizon tal black lines show the allo we d interv al for the minima and the maxima: this information is separately av ailable for all 2 n time sub -in terv als partitioning [0 , 1] . Note that the last graph corresponds to n = 12, with the subplot in its frame correspond ing to a zo oming on the p ositio n of the paths on the time interv al [0 . 424 , 0 . 434]. ε -St r ong simulation of the Br ownian p ath 9 Then: 2 n/ 2 × E[ | X ↑ ( n ) − X ↓ ( n ) | 1 ] = O (1) . Pro of. W e pro ceed as follows: | X ↑ ( n ) − X ↓ ( n ) | 1 = 2 n X i =1 ( U ↑ ( i − 1)2 − n ,i 2 − n − L ↓ ( i − 1)2 − n ,i 2 − n ) · 2 − n (3.2) ≤ 2 n X i =1 ( M ( i − 1)2 − n ,i 2 − n − m ( i − 1)2 − n ,i 2 − n + 2 · 2 − n/ 2 ) · 2 − n , the inequality being a direct consequence of Step 2.ii of the ε -s trong algorithm in T able 2 . Consider now the path from X ( i − 1)2 − n to X i 2 − n . Let Z b e a Brownian bridge from Z 0 = 0 to Z 2 − n = 0; we denote by M z and m z its maximum and minim um, resp ectiv ely . Conditionally on X ( i − 1)2 − n and X i 2 − n , a known pro perty of the Br o wnian bridge implies (see, e.g., [ 14 ]) that: X t +( i − 1)2 − n = Z t +  1 − t 2 − n  X ( i − 1)2 − n + t 2 − n X i 2 − n , t ∈ [0 , 2 − n ] , in the sense tha t the pro cesses on the tw o s ides of the ab o ve e q uation hav e the same distribution. It is now clear that: M ( i − 1)2 − n ,i 2 − n − m ( i − 1)2 − n ,i 2 − n ≤ | X i 2 − n − X ( i − 1)2 − n | + ( M z − m z ) . So, taking exp ectations at ( 3.2 ), we get: E[ | X ↑ ( n ) − X ↓ ( n ) | 1 ] ≤ E [ M z − m z ] + 2 · 2 − n/ 2 + 2 n X i =1 E | X i 2 − n − X ( i − 1)2 − n | 2 − n The finite-dimensional distributions o f the initial Brownian br idge from X 0 to X 1 imply that: X i 2 − n − X ( i − 1)2 − n ∼ N (( X 1 − X 0 )2 − n , 2 − n (1 − 2 − n )) , which gives directly that: E | X i 2 − n − X ( i − 1)2 − n | = O (2 − n/ 2 ) . It remains to show that E[ M z − m z ] = O (2 − n/ 2 ) to co mplete the pro of. Now, self- similarity of Brownian motion implies that: Z u = 2 − n/ 2 ˜ Z u/ 2 − n , 10 A. Beskos, S. Peluchetti and G. R ob erts where ˜ Z is a Brownian bridg e from ˜ Z 0 = 0 to ˜ Z 1 = 0. Let ˜ M z , ˜ m z be the maximum and minim um of ˜ Z . Due to the self-similarity , we have M z − m z = 2 − n/ 2 ( ˜ M z − ˜ m z ) . Since ˜ M z − ˜ m z in a random v ariable of finite exp ectation (see, e.g., [ 14 ]), we obta in directly that E[ M z − m z ] = O (2 − n/ 2 ) which completes the pro of.  4. The ζ -function W e hav e yet to present the s ampling metho ds employed when refining or bisecting an int ers ection lay er during the exec ution of the ε -strong algo rithm, thu s constituting the building blo cks of o ur algorithm. All probabilities in volved in these metho ds can b e expressed in terms of a hitting probability of the Brownian path. W e denote by W ( l,x,y ) the proba bilit y law o f a Br o wnian bridge from X 0 = x to X l = y . Let ζ ( L, U ; l , x, y ) , with L < U , b e the pro babilit y that the B r o wnian bridge escap es the interv al [ L, U ] . That is: ζ ( L , U ; l , x, y ) = W ( l,x,y ) [ m 0 ,l < L or M 0 ,l > U ] . W e als o define: γ ( L, U ; l, x, y ) = 1 − ζ ( L, U ; l , x, y ) . (4.1) These pr obabilities ca n be calcula ted analy tically in terms of an infinite series . The res ult is ba sed on a partition of Brownian paths w.r.t. to a trace they leav e on tw o b ounding lines and can b e attributed back to [ 11 ]; for more recent r eferences see [ 1 , 9 , 17 ]. W e define for j ≥ 1, ¯ σ j ( x, y , δ, ξ ) = exp  − 2 l [ δ j + ξ − x ][ δ j + ξ − y ]  , (4.2) ¯ τ j ( x, y , δ ) = exp  − 2 j l [ δ 2 j + δ ( x − y )]  . Then, Theo rem 3 of [ 17 ] yie lds ζ ( L , U ; l , x, y ) =      ∞ X j =1 ( σ j − τ j ) , L < x, y < U, 1 , otherwise , (4.3) where σ j = ¯ σ j ( x, y , U − L , L ) + ¯ σ j ( − x, − y , U − L, − U ) , (4.4) τ j = ¯ τ j ( x, y , U − L ) + ¯ τ j ( − x, − y , U − L ) . ε -Str ong simulation of the Br ownian p ath 11 The infinite series in ( 4.3 ) exhibits a mo notonicit y prop erty which will b e e x ploited by our simulation alg orithms. W e consider the seq ue nce { S n } , with S n = S n ( L, U ; l , x, y ), defined a s: S 2 n − 1 = n − 1 X j =1 ( σ j − τ j ) + σ n , S 2 n = S 2 n − 1 − τ n , (4.5) when L < x, y < U , otherwise S n ≡ 1. Then: 0 < S 2 n ≤ S 2 n +2 ≤ ζ ≤ S 2 n +1 ≤ S 2 n − 1 (4.6) for all n ≥ 1 ; fo r a pro of see [ 9 ] or [ 5 ]. 4.1. ζ -Derived even ts W e can combine ζ - probabilities to calcula te other co nditional probabilities ar ising in the context of the ε - strong algorithm. W e beg in with the following definition: β ( L ↓ , L ↑ , U ↓ , U ↑ ; l , x, y ) := W ( l,x,y ) [ L ↓ < m 0 ,l < L ↑ , U ↓ < M 0 ,l < U ↑ ] . Now, we hav e the set e qualit y: { L ↓ < m 0 ,l < L ↑ , U ↓ < M 0 ,l < U ↑ } (4.7) = { L ↓ < m 0 ,l , M 0 ,l < U ↑ } − { L ↑ < m 0 ,l , M 0 ,l < U ↑ } ∪ { L ↓ < m 0 ,l , M 0 ,l < U ↓ } . Thu s, taking probabilities and reca lling the definition of γ in ( 4.3 ), we find that: β ( L ↓ , L ↑ , U ↓ , U ↑ ; l , x, y ) = γ ( L ↓ , U ↑ ; l , x, y ) − γ ( L ↑ , U ↑ ; l , x, y ) (4.8) − γ ( L ↓ , U ↓ ; l , x, y ) + γ ( L ↑ , U ↓ ; l , x, y ) . Before the next even t, we enrich the notation for the Brownian bridge measur e. W e define (for 0 < q < l ): W ( l,x,y ) ( q,w ) [ · ] = W ( l,x,y ) [ · | X q = w ] . W e set r = l − q . Co nsider now the conditional probability: ρ ( L ↓ , L ↑ , U ↓ , U ↑ ; q , r , x, w , y ) = W ( l,x,y ) ( q,w ) [ L ↓ < m 0 ,l < L ↑ , U ↓ < M 0 ,l < U ↑ ] . Using ag ain the set equality ( 4.7 ), a nd taking proba bilities under W ( l,x,y ) ( q,w ) , we obtain: ρ ( L ↓ , L ↑ , U ↓ , U ↑ ; q , r , x, w , y ) = γ 1 γ 2 − γ 3 γ 4 − γ 5 γ 6 + γ 7 γ 8 , (4.9) 12 A. Beskos, S. Peluchetti and G. R ob erts where we hav e defined: γ 1 = γ ( L ↓ , U ↑ ; q , x, w ) , γ 2 = γ ( L ↓ , U ↑ ; r , w , y ) , γ 3 = γ ( L ↑ , U ↑ ; q , x, w ) , γ 4 = γ ( L ↑ , U ↑ ; r , w , y ) , γ 5 = γ ( L ↓ , U ↓ ; q , x, w ) , γ 6 = γ ( L ↓ , U ↓ ; r , w , y ) , γ 7 = γ ( L ↑ , U ↓ ; q , x, w ) , γ 8 = γ ( L ↑ , U ↓ ; r , w , y ) . Note that the pr oduct terms ar ise due to the indep endency of the B ro wnian bridges on [0 , q ] and [ q , l ]. W e will b e using these expressio ns for β ( · ; · ) a nd ρ ( · ; · ) in the s equel. 4.2. Sim ulation of ζ -deriv ed ev en ts W e will need to b e a ble to decide whether even ts of probability ζ ha ve o ccurred or no t. In a simulation context, this co r respo nds to deter mining the v alue of the bina ry v a r iable I R<ζ for R ∼ Un [0 , 1]. With ( 4.6 ) in mind, we define: J = inf { n ≥ 1 : n o dd, S n < R o r n even, S n > R } . Due to the alternating monotonicity pr operty ( 4.6 ) o f S n : I R<ζ = I J is even . Thu s, we need a.s. finite n umber o f J computations to ev aluate I R<ζ . Note that S n conv erges to its limit ex ponentially fast, so J will b e of small exp ectation; o ne ca n ea sily verify that all its moments ar e finite. Such an approach was also followed in [ 5 ]. In a more gener al context, we will also b e r equired to decide if ev ents of pro babil- it y β ( · ; · ) o r ρ ( · ; · ) hav e taken place or not; w e will in fact be co ns idering even more complex e v ents r elated with the ζ -function. In the most encompass ing scena rio, when executing o ur sampling metho ds, we will b e re q uired to compare a given rea l num ber R with Z ( ζ 1 , ζ 2 , . . . , ζ m ) for some given function Z , with the different ζ i ’s corres ponding to different choices of the ar g umen ts l , x, y , L, U for ζ ( · ; · ). Using the monotonicity pro perty ( 4.6 ), we will b e able to develop cor respo nding a lternating sequences S Z n such that: S Z 2 n ≤ S Z 2 n +2 ≤ Z ( ζ 1 , ζ 2 , . . . , ζ m ) ≤ S Z 2 n +1 ≤ S Z 2 n − 1 ; (4.10) lim n →∞ S Z n = Z ( ζ 1 , ζ 2 , . . . , ζ m ) , and pro ceed as a bov e. Analytica lly , we will determine the v alue of the compa r ison bina ry indicator I R R (whenc e r eturn 1 ) . ε -St r ong simulation of the Br ownian p ath 13 5. Distributions and their sim ulation W e will now desc r ibe a nalytically all simulation algo rithms employ ed at the development of the ε -strong a lgorithm prese nted in T a ble 2 . In particular , one has to develop sampling metho ds to carry out the r efinemen t and bise ction (see Section 2 ) of the intersection lay er I s,t . T o simplify the presentation, when co nditio ning on X s , X t ∗ or X t we will make the corres p ondence: x = X s , w = X t ∗ , y = X t , l = t − s, q = t ∗ − s, r = t − t ∗ . 5.1. Bisection of I s,t : Sampling the middle poin t X t ∗ Bisection of I s,t = { X s , X t , L s,t , U s,t } , with L s,t = [ L ↓ , L ↑ ], U s,t = [ U ↓ , U ↑ ], b egins by sampling a p oint of the Brownian bridge conditionally on the collected infor mation ab out its minimum and maximum; this is Step 1 of T able 1 . Such a conditiona l distributio n is analytically tra ctable via Bayes’ theorem. Prop osition 5. 1. The distribution W [ X t ∗ | I s,t ] , with t ∗ ∈ [ s, t ] , has pr ob ability density: f ( w ) ∝ ρ ( L ↓ , L ↑ , U ↓ , U ↑ ; q , r , x, w , y ) × π ( w ) wher e ρ ( · ; · ) is define d in ( 4.9 ) and π ( w ) = exp  − 1 2  w −  r l x + q l y  2 .  q r l  . Pro of. The function π ( w ) corres ponds to the prio r (unnormalised) density for the middle po in t X ∗ t | X s , X t which is easily found to b e norma lly distributed with mean a nd v a riance as implied by the expres sion for π ( w ) . So, following the definition of ρ ( · ; · ) in ( 4.9 ), the stated res ult is an application of B ayes’ theore m.  W e will develop a metho d fo r sa mpling from f ( w ). It is eas y to co nstruct an alterna ting series b ounding f ( w ). Let: ζ i = 1 − γ i , 1 ≤ i ≤ 8 , for the eight γ -functions app earing at the definitio n of ρ in ( 4.9 ). Let { S i,n } n ≥ 1 be the alternating ser ies ( 4.6 ) for ζ i , for each 1 ≤ i ≤ 8; that is: 0 < S i, 2 n ≤ S i, 2 n +2 ≤ ζ i ≤ S i, 2 n +1 ≤ S i, 2 n − 1 , (5.1) with lim n →∞ S i,n = ζ i . Co nsider the s equence { S Z n } defined a s follows: S Z n = (1 − S 1 ,n +1 − S 2 ,n +1 + S 1 ,n S 2 ,n ) − (1 − S 3 ,n − S 4 ,n + S 3 ,n +1 S 4 ,n +1 ) (5.2) − (1 − S 5 ,n − S 5 ,n + S 5 ,n +1 S 6 ,n +1 ) + (1 − S 7 ,n +1 − S 8 ,n +1 + S 7 ,n S 8 ,n ) . 14 A. Beskos, S. Peluchetti and G. R ob erts Due to ( 5.1 ), one can easily verify that { S Z n } is an a lternating sequence for ρ ( · ; · ) , in the sense that: S Z 2 n ≤ S Z 2 n +2 ≤ ρ ( L ↓ , L ↑ , U ↓ , U ↑ , q , r , x, w , y ) ≤ S Z 2 n +1 ≤ S Z 2 n − 1 (5.3) with lim n →∞ S Z n = ρ ( L ↓ , L ↑ , U ↓ , U ↑ ; q , r , x, w , y ) . W e explo it this structure to build a rejection sa mpler to dr aw from the density f ( w ) in P ropos ition 5 .1 . W e will use prop osals from: f 2 n +1 ( w ) = S Z 2 n +1 ( w ) × π ( w ) , where we hav e emphasized the dep endence of S Z 2 n +1 on the argument w . Note that the domain o f b oth f ( w ), f 2 n +1 ( w ) is [ L ↓ , U ↑ ]. Now, we will illustra te that S Z 2 n +1 ( w ) has a c oncrete structur e that we will exploit for our sa mpler. Consider the fir st of the four terms for ming up S Z 2 n +1 from ( 5.2 ): 1 − S 1 , 2 n +2 − S 2 , 2 n +2 + S 1 , 2 n +1 S 2 , 2 n +1 . (5.4) F o llo wing the analytical definition of the alterna ting sequences in equations ( 4.2 ), ( 4.4 ), ( 4.5 ), b oth S 1 ,n and S 2 ,n , ca n b e express e d a s a s um of 2 n terms each having the exp o- nent ial structur e ± e x p { a + bw } I L ↓ E · e E , E ∼ E xp(1) , (6.1) with E b eing indep enden t of X , is an unbiased estimator o f E[ F ( X )] . The ε -strong algorithm co uld b e utilised her e to unbiasedly obtain the v alue of the binary v ariable I F ( X ) >E in finite co mputations. W e can easily find the seco nd moment o f the unbiased estimator in ( 6.1 ): E[I F ( X ) >E · e 2 E ] = E [ e F ( X ) ] − 1 . (6.2) W e descr ibe for a moment in mor e detail the identification of I F ( X ) >U via the ε - strong algorithm. Utilising the lo wer and upp er con vergent pro cesses X ↓ ( n ), X ↑ ( n ) in ( 3.1 ) one could in many c ases analytica lly identif y qua n tities F ↓ n , F ↑ n (realisable with finite computations) such that: F ↓ n ≤ F ↓ n +1 ≤ F ( X ) ≤ F ↑ n +1 ≤ F ↑ n ; F ↑ n − F ↓ n → 0 . Given enoug h iterations, there will b e a greemen t; for the a .s. finite random instance: κ = inf { n ≥ 0: I F ↓ n >E = I F ↑ n >E } (6.3) we will hav e I F ( X ) >E = I F ↓ κ >E . (6.4) Thu s, combining ( 6.1 ) w ith ( 6 .4 ), we hav e developed an un bias e d es timator o f a path exp ectation, inv olving finite c o mputations. Certainly , the num er ical efficiency o f such a n estimation will r ely heavily on the sto c has tic pro perties of κ and the cost o f gener ating F ↓ n , F ↑ n , a nd of cour se the v a riance of the estimator . The particular deriv ation of the ab o ve un bias ed estimator o f the path ex pectation is by no means restrictive; o ne can genera te unbiased es timators using distr ibutions o ther than the exp onen tial. Consider the following scenar io. W e can genera te some pr eliminary bo unds F ↓ n 0 , F ↑ n 0 up to some fixed or random (dep e nding on X ) instance n 0 . Now, o ne can eas ily check (by considering the conditiona l exp ectation w.r.t. R | X ) that: I F ( X ) >R F ↑ n 0 + I F ( X ) 0, we get that the second mo men t of the estimator ( 6.5 ) will be: E[I F ( X ) >R ( F ↑ n 0 ) 2 + I F ( X ) 0 (volatilit y) and a Brownian motion { W t } . Als o , T ab o ve is the maturity time, K S the strike pr ice and L S , U S the lower and upp er bar r iers resp ectiv ely; suprema and infima ar e considered ov er the time p erio d [0 , T ] . Note tha t E[ F b ( S )] co rresp onds to the price of the Asian option, see, for exa mple, [ 21 ]. The pro cess S t is an 1–1 transfor mation of a Br o wnian motio n with drift. In particular, we ca n rewrite the functionals ( 6.6 )–( 6 .8 ) a s follows: F a ( X ) = e − r T (e σ sup X t − K S ) + I L< inf X t < sup X t R in ( 6.5 ) is dec ided. P ropo sition 3 .2 will be of r elev ance in this context. Rec a ll that κ in ( 6.3 ) denotes the num b er of steps to decide ab out I F ( X ) >R . T able 7. Similar results as for T able 5 , bu t n o w for the case of E[ F c ( S )] in ( 6.7 ) Euler app ro ximation δ Time (secs) 95% Conf. I n t. 1 / 10 0.5 [797 , 807] × 10 − 4 1 / 20 0.9 [822 , 832] × 10 − 4 1 / 40 1.7 [833 , 844] × 10 − 4 1 / 80 3.3 [835 , 846] × 10 − 4 1 / 160 6.6 [846 , 858] × 10 − 4 ε -strong n 0 Time (secs) 95% Conf. Int. 2 178.0 [842 , 854] × 10 − 4 ε -St r ong simulation of the Br ownian p ath 23 The co st of κ iterations of the ε -str ong alg orithm (when its full machinery is req uir ed) is prop ortional to K = 2 κ . In the context of ( 6.5 ), we find: P[ K > 2 n ] = E [P[ κ > n | X ]] = E  F ↑ n − F ↓ n F ↑ n 0 − F ↓ n 0  . Prop osition 3.2 states that | X ↑ ( n ) − X ↓ ( n ) | L 1 = O (2 − n/ 2 ). The same r ate of co n vergence will many times also b e true for E[ F ↑ n − F ↓ n ]: this will be the case for instance when F ( X ) = f ( R 1 0 g ( X s ) d s ) under g eneral assumptions o n f , g (e.g., if | f ( y ) − f ( x ) | ≤ M ( x, y ) | y − x | , for a po lynomial M , and the same for g ; a pr o of is not essential her e ). F or such a rate (and since the user-sp ecified F ↑ n 0 − F ↓ n 0 should b e ea sily c o n trolled), we will get: P[ K > 2 n ] = O (2 − n/ 2 ) giving the infinite exp ectation E [ K ] = ∞ . In a g iven application thoug h, o ne could fix a big eno ugh maximum num ber n max , stop the bisections if that num b er has b een r eac hed and rep ort, say , ( F ↑ n 0 + F ↓ n 0 ) / 2 as the realiza tio n of the estimator if that happ ens, witho ut practical effect on the r e s ults. T o explain this, note tha t we know, from ( 6.5 ), that the actual unbiased v alue is either F ↑ n 0 or F ↓ n 0 , s o we know pr e cisely that the absolute bias from the single realiza tio n when n max was r eac hed canno t be greater than ( F ↑ n 0 − F ↓ n 0 ) / 2. In total, when averaging over a num b er of realiza tio ns we can hav e a precise arithmetic b ound on the absolute v a lue of the bia s of the rep orted average; if n max is ‘big eno ugh’ so that we reach n max only in a sma ll prop ortion of re alizations the (analy tically known) bias could be o f such a magnitude that the rep orted results will b e pr ecisely the sa me as when implementing the regular algor ithm without n max for a r easonably selected deg ree of acc ur acy . F or example, in the case of the estimation of E [ F c ( S )] in T able 6 , we hav e in fact used n max = 10 a nd found tha t the intro duce d bia s was sma ller than 3 × 10 − 5 so avoiding it would not make any difference or whatso ever at the res ults re p orted r igh t now in T able 6 . Note that such an issue did no t aris e in the cas es of E [ F a ( S )] and E[ F c ( S )] when a reduced version of the ε -s tr ong alg orithm was applied. 7. F urther directions for applications W e sketch here some o ther p otential a pplications of the ε -strong algo rithm. In the case of bar rier options, so metimes one needs to ev aluate exp ectations inv olving a Brownian hitting time (see, e.g., [ 19 ]). Given a nonconstant b oundary H : [0 , ∞ ) → R , such that S 0 < H 0 , co nsider: τ H = inf { t ≥ 0: S t ≥ H t } , with S = { S t } be ing the g eometric BM in ( 6.9 ). The price of a related deriv ative will b e E[ F ( S )] where now: F ( S ) = ψ ( S T ) · I τ H

Original Paper

Loading high-quality paper...

Comments & Academic Discussion

Loading comments...

Leave a Comment