Analysis, Geometry, and Modeling in Finance 15

Appendix B Monte-Carlo Methods and Hopf Algebra B.1 Introduction Since the introduction of the Black-Scholes paradigm, several alternative mod- els which allow to better capture the risk of exotic options have emerged : local volatility models, stochastic volatility models, jump-diffusion models, mixed stochastic volatility-jump diffusion models, etc. With the growing complica- tion of Exotic and Hybrid options that can involve many underlyings (equity assets, foreign currencies, interest rates), the Black-Scholes PDE, which suf- fers from the curse of dimensionality (the dimension should be strictly less than four in practice), cannot be solved by finite difference methods. We must rely on Monte-Carlo methods. This appendix is organized as follows: In the first section, we review basic features in Monte-Carlo simulation from a modern (algebraic) point of view: generation of random numbers and dis- cretization of SDEs. A precise mathematical formulation involves the use of the Taylor-Stratonovich expansion (TSE) that we define carefully. Details can be found in the classical references [20], [28]. In the second section, we show that the TSE can be framed in the setting of Hopf algebras.1 In particular, TSEs define group-like elements and solutions of SDEs can be written as exponentials of primitive elements (i.e., elements of the universal Lie algebra associated to the Hopf algebra). This section is quite technical and can be skipped by the reader. This Hopf algebra structure allows to prove easily the Yamato theorem that we explained in the last section. As an application, we classify local volatility models that can be written as a functional of a Brownian motion and therefore can be simulated exactly. The use of Yamato’s theorem allows us to reproduce and extend the results found by P. Carr and D. Madan in [69]. 1The author would like to thank his colleague Mr. C. Denuelle for fruitful collaborations on this subject and for help in the writing of this appendix. 353 © 2009 by Taylor & Francis Group, LLC 354 Analysis, Geometry, and Modeling in Finance B.1.1 Monte Carlo and Quasi Monte Carlo Let Xt be an n-dimensional Itoˆ process. According to (2.31), the pricing of a derivative product requires the evaluation of E[f(XT )] (B.1) Through a discretization scheme discussed in the next paragraph, we assume that XT can be approached by a function g of a multidimensional Gaussian variable G ∼ N(0,Σ). To figure out an approximation of (B.1), one simulates M independent random variables (Gm ∼ N(0,Σ))1≤m≤M and computes the associated values (XmT )1≤m≤M . Under reasonable hypothesis on the payoff f and the function g, the strong version of the law of large numbers [26] entails that the empirical average 1M ∑M m=1 f(X m T ) provides a good proxy for E[f(XT )]. Better yet the central limit theorem [26] states that the error εM = 1 M M∑ m=1 f(XmT )− E[f(XT )] is asymptotically normally distributed, with mean 0 and covariance tΣΣ M . Thus there exists a constant C such that the L2 Monte Carlo error is ||εM ||L2 ≤ C√ M Even though the O( 1√ M ) bound is steady, variance reduction techniques such as control variates can considerably reduce the constant C. We refer to [20] for further discussions on these methods and the way they speed up the con- vergence of MC estimates. As liquidity on the market expands, option pricing requires increasing pre- cision, and the O( 1√ M ) MC bound is simply not good enough. To gain a factor 10 in accuracy the number of simulations must be increased by a fac- tor 100. To elude this pitfall, Quasi Monte-Carlo (QMC) methods have been developed [20]. They consist in low-discrepancy sequences (Van Der Corput, Faure, Sobol...) filling up (0, 1)d uniformly. Contrary to what is often as- sumed these sequences are not at all random. QMC trajectories are much too uniform to be random. Standard transformations convert the (0, 1)d- valued sequences into d-dimensional Gaussian sequences, and lead to well chosen solutions (XmT )1≤m≤M . The error coming from the QMC estimator 1 M ∑M m=1 f(X m T ) can be as good as O( (lnM)d M ), which is nearly O( 1 M ). B.1.2 Discretization schemes For the sake of simplicity, we assume that we are trying to simulate on the time interval [0, T ] a stochastic process driven by the one dimensional SDE : dXt = µ(Xt)dt+ σ(Xt)dWt © 2009 by Taylor & Francis Group, LLC Monte-Carlo Methods and Hopf Algebra 355 Let (tn)0≤n≤N be a subdivision of [0, T ] (∆nt = tn+1 − tn)0≤n≤N−1 and ∆W = (∆nW = Wtn+1 −Wtn)0≤n≤N−1 From the basic properties of the Brownian motion it is clear that ∆W ∼ N(0,∆t) where ∆t is the diagonal matrix whose entries are the (∆nt). MC simulation or QMC can therefore be employed to draw M paths for ∆W . As stated in the previous paragraph, discretization schemes will permit the transformation of these paths into paths of the underlying process. The SDE above yields Xtn+1 = Xtn + ∫ tn+1 tn µ(Xt)dt+ ∫ tn+1 tn σ(Xt)dWt Hence a very intuitive discretization is provided by the well-known Euler scheme ([28],[20]) which generates the following path X conditionally on ∆W : Xt0 = X0 Xtn+1 = Xtn + µ(Xtn)∆ nt+ σ(Xtn)∆ nW However the Euler scheme does have an inconvenient: the two approximations it consists in are not of the same order∫ tn+1 tn µ(Xt)dt = µ(Xtn)∆ nt+ o(∆nt) whereas ∫ tn+1 tn σ(Xt)dWt = σ(Xtn)∆ nW + o( √ ∆nt) It would therefore make sense to develop :∫ tn+1 tn σ(Xt)dWt = ∫ tn+1 tn σ ( Xtn + ∫ t tn µ(Xs)ds+ ∫ t tn σ(Xs)dWs ) dWt = σ(Xtn)∆ nW + σ′(Xtn)σ(Xtn) ∫ tn+1 tn ( ∫ t tn dWs)dWt + o(∆nt) = σ(Xtn)∆ nW + σ′(Xtn)σ(Xtn) 2 ( (∆nW )2 −∆nt)+ o(∆nt) This refinement is at the origin of the one-dimensional Milstein scheme: X̂t0 = X0 X̂tn+1 = X̂tn + µ(X̂tn)∆ nt+ σ(X̂tn)∆ nW + σ′(X̂tn)σ(X̂tn) 2 ( (∆nW )2 −∆nt) © 2009 by Taylor & Francis Group, LLC 356 Analysis, Geometry, and Modeling in Finance To study the performance of discretization schemes the literature distinguishes between strong and weak order of convergence ([28], [20]): DEFINITION B.1 Strong/Weak weak order of convergence Not- ing h = max n=0,..,N−1 (∆nt), a scheme X¯ is said to be of strong order β > 0 if there exists a constant c such that E[‖ X¯(tN )−X(T ) ‖] ≤ chβ for some vector norm ‖ · ‖. It is said to be of weak order β > 0 if there exists a constant c such that | E[f(X¯(tN )]− E[f(X(T ))] |≤ chβ for all f polynomially bounded in C2β+2. Under Lipschitz-type conditions on µ and σ the Euler scheme can be shown to be of strong order 12 and of weak order 1, whereas the Milstein scheme is of strong and weak order 1. Writing a multi-dimensional generalization of the Milstein scheme would in- volve adding antisymmetric Wiener iterated integrals, called Le´vy areas Aij = ∫ 1 0 W it dW j t − ∫ 1 0 W jt dW i t (B.2) The standard approximation for this term requires several additional random numbers [28]. There are however approaches to avoid the drawing of many extra random numbers by using the relation of this integral to the Le´vy area formula [124]: E[eiλA ij |W i1 + ıW j1 = z] = λ sinhλ e− |z|2 2 (λ cothλ−1) From this characteristic function, Le´vy areas and Brownian motion W it can be simulated exactly. However, this can become quite costly in computational expense. B.1.3 Taylor-Stratonovich expansion Using the tools described in the previous sections, we explain the Taylor- Stratonovich expansion (TSE) equivalent to a Taylor expansion in the de- terministic case. Note that a similar expansion, called Taylor-Itoˆ expansion, that uses the Itoˆ calculus exists. In order to be simulated, the discretization scheme eventually needs to be set in Itoˆ form. However, we prefer to present TSE as its algebraic structure is simpler. © 2009 by Taylor & Francis Group, LLC Monte-Carlo Methods and Hopf Algebra 357 Let Xt be an n-dimensional Itoˆ process following the Stratonovich SDE dXt = V0dt+ m∑ i=1 Vi�dW it , Xt=0 = X0 ∈ Rn By introducing the notation dW 0t ≡ dt, we have dXt = m∑ i=0 Vi�dW it (B.3) Without loss of generality, we assume that we have a time-homogeneous SDE. A time-inhomogeneous SDE can be written in this normal form (B.3) by including an additional state Xn+1t dXn+1t = dt The Stratonovich calculus entails that for f in C1(Rm,R) f(XT ) = m∑ i=0 ∫ T 0 Vif(Xt)�dW it By iterating this equation, we obtain f(XT ) = f(X0) + m∑ i=0 ( Vif(X0) ∫ T 0 �dW it + m∑ j=0 ∫ T 0 (∫ t 0 ViVjf(Xs)�dW js ) �dW it  With a repeated application of the Stratonovich formula, the iterated integral defined as ∫ 0≤t1<... r. This definition takes its sense from the fact that the exponential of a primitive element in Gr is a group-like element in Gr: PROPOSITION B.1 Gr = exp(Gr) © 2009 by Taylor & Francis Group, LLC Monte-Carlo Methods and Hopf Algebra 363 PROOF We do the proof in one way: If L ∈ Gr then ∆(exp(L)) = ∆( ∞∑ k=0 Lk k! ) = ∞∑ k=0 1 k! ∆(L)k = ∞∑ k=0 1 k! (L ⊗ ε+ ε⊗ L)k = ∞∑ k=0 1 k! k∑ j=0 k! j!(k − j)! (L j ⊗ Lk−j) = exp(L)⊗ exp(L) B.2.2 Chen series Now that this setting has been introduced we can almost reach our goal and see how it appears in TSE. The idea is to use a morphism Γr : εi1 ...εik .ε ∈ Hr 7−→ Vi1 ...Vik(X0) For a semi-martingale ω we define its Chen series in Hr as the pre-image of the Taylor-Stratonovich expansion: X0,1(ω) = ∑ (i1,...,ik)∈Ar εi1 ...εik ∫ 0≤t1<...