Felisa J. V azquez-Abad Department of Computer Science and - PowerPoint PPT Presentation

Felisa J. V´ azquez-Abad and Daniel Dufresne 0 Accelerated Simulation fo r Pricing Asian Options Accelerated Simulation fo r Pricing Asian Options Felisa J. V´ azquez-Abad Department of Computer Science and Operations Research University of Montreal, Montr´ eal, Qu´ ebec H3C 3J7 e-mail: vazquez@iro.umontreal.ca and Daniel Dufresne Department of Economics University of Melbourne e-mail: dufresne@myriad.its.unimelb.edu.au Invited paper to the Winter Simulation Conference , Washington, D.C., December 1998.

Felisa J. V´ azquez-Abad and Daniel Dufresne 1 Introduction Two assets: a risk-free asset yielding interest r , and a risky asset, with price Geometric Brownian Motion S t = S 0 e µt + σW t . S t • European Call Option: Op- K tion to buy one unit of stock at strike price K at the exercise time T . Only if S T > K the holder of the option will exercise it. Profit: ( S T − K ) + . • European Average Call Op- tion: Instead of the terminal S 0 value S T , the average is compared T t=0 to K . P , discounted asset prices { e − rt S t } form Asset Pricing: Under the risk-neutral measure a martingale. The price π is the expectation w.r.t. P of the discounted gain. Discrete model:   √  ( r + σ 2 h = T   S i +1 = S i exp 2 ) h + σ hZ i  , Z i i.i.d. ∼ N (0 , 1) , N      1 N �  e − rT   π = i =1 S i − K If K >> S 0 out of the money ⇒ rare event estimation. E   N +

Felisa J. V´ azquez-Abad and Daniel Dufresne 2 Simulation Methods I Z i ∼ N (0 , 1) , i = 1 , . . . , N, are independent standard normal variables. √ 1 ≤ i ≤ N, where u σ = u − σ 2 X u ⇒ X u i ∼ N ( u σ h, σ 2 h ) , = u σ h + σ hZ i , i 2 S u = S u i − 1 exp( X u i ) , 1 ≤ i ≤ N, are the asset prices, i 1 N � i =1 S u A u = is the arithmetic average of asset prices, i N Na¨ ıve Estimation: Take u = r above, and use D 0 = Y r 1 = e − rT ( A r − K ) + . Confidence Intervals: Approximate Confidence Interval: apply the Central Limit Theo- rem to M replications of the simulation. Use a confidence level α = 0 . 05. � � V a r [ D 0 ] π ± 1 . 96 ˆ M The precision of the estimation can be reduced if we use other estimators with smaller variance than the na¨ ıve estimator. We shall discuss the methods: • Control Variable • Change of Measure • Hybrid Estimation

Felisa J. V´ azquez-Abad and Daniel Dufresne 3 Simulation Methods II √ X u Z i ∼ N (0 , 1) , i = 1 , . . . , N, i = u σ h + σ hZ i , S u = S u i − 1 exp( X u i ) , 1 ≤ i ≤ N, i 1 N � i =1 S u A u = is the arithmetic average of asset prices, i N 1   N N � i =1 S u G u = is the geometric average of asset prices.   i The Control Variable Method: Take u = r , and let Y r 2 be some other estimator Y r 2 whose expected value is known. Use: 2 − Y r D 1 = Y r E Y r 1 + α ( 2 ) , which implies that E D 1 = π Cov ( Y r 1 , Y r r Y r It can be shown that α = 2 ) / V a 2 minimizes the variance. Use as a control variable the option price when the geometric mean is used in the average, that is: Y r 2 = e − rT ( G r − K ) + . Mean and variance of Y r 2 are known. Variance reduction:   Cov ( Y r 1 , Y r 2 )    1 −  ≤ V a r [ D 1 ] = V a r [ D 0 ] V a r [ D 1 ]   � r ( Y r r ( Y r V a 1 ) V a 2 )

Felisa J. V´ azquez-Abad and Daniel Dufresne 4 Simulation Methods III Z i ∼ N (0 , 1) , i = 1 , . . . , N, B i = Z 1 + · · · + Z i , 1 ≤ i ≤ N, B 0 = 0 √ A u = 1 N � X u S u i = S u i − 1 exp( X u i =1 S u = u σ h + σ hZ i , i ) , 1 ≤ i ≤ N, i . i N The Change of Measure Approach: Girsanov’s Theorem for one-dimensional Brow- nian Motion implies for the discrete model: E e − Nv 2 2 − vB N f ( Z 1 + v, . . . , Z N + v ) . ∀ v ∈ I R, E f ( Z 1 , . . . , Z N ) = √ Let v = ( u − r ) h/σ , then for f ( Z 1 , . . . , Z N ) = A r = 1 f ( Z 1 + v, . . . , Z N + v ) = 1 N N � i − 1 e X r � i − 1 e X u i =1 S r i , i =1 S u i = A u N N √ √     2    ( u − r ) − ( u − r )  − N  h h  Call L u = e − Nv 2   2 − vB N = exp   B N  , then it follows that:    2 σ σ   D 2 ( u ) = L u Y u 1 = L u ( A u − K ) + is unbiased for π , ∀ u ∈ I R. Now we simulate at drift u σ instead of r σ . The optimal value of u is that which minimizes V a r D 2 ( u ). We shall write: D 2 = D 2 ( u ∗ ) Likelihood Ratio Estimators, rare events: Importance Sampling.

Felisa J. V´ azquez-Abad and Daniel Dufresne 5 Simulation Methods IV The Hybrid Estimators: Adding a control variable to D 2 , we simulate at u and thus we can calculate G u , yielding another unbiased estimator: D 3 = D 3 ( u ∗ ) D 3 ( u ) = L u Y u E Y u 2 − Y u 1 + α ( 2 ) , where u ∗ minimizes the variance. Finally, consider changing the measure of the controlled estimator D 1 : D 4 ( u ) = L u Y u E Y r 2 − L u Y u D 4 = D 4 ( u ∗ ) 1 + α ( 2 ) , where u ∗ minimizes the variance. Summarizing: • D 0 : Na¨ ıve Estimator • D 1 : Controlled Estimator • D 2 : Na¨ ıve Estimator under Importance Sampling • D 3 : Controlled Likelihood Ratio Estimator • D 4 : Likelihood Ratio Estimator under Importance Sampling

Felisa J. V´ azquez-Abad and Daniel Dufresne 6 Efficient Importance Sampling The change of measure approach yields estimators that can improve efficiency. When rare events are involved, Importance Sampling will simulate more often the rare situations and L u weighs appropriately the estimate to yield unbiasedness. Problem: For a likelihood estimator D ( u ), the optimal value of u ∗ is problem-dependent.  ⇒ Find u ∗ : ∂   Minimize: V a r [ D ( u )] , V a r [ D ( u )] u = u ∗ = 0  ∂u Solutions: • The first solution and most commonly used is to perform pilot simulations. We used Functional Estimation to save computational time. Use the same random numbers to evaluate in parallel the Likelihood Ratio Estimators for different values of u and each value of K . Z i ∼ N (0 , 1) , B i = Z 1 + · · · + Z i , 1 ≤ i ≤ N, i = 1 , . . . , N, B 0 = 0 √ � � − ( u − r )2 T √ A u = S 0 − ( u − r ) hBN N i =1 e i ( u + σ � 2 ) h + σ hB i , 2 σ 2 σ L u = e N • New Approach: optimize the parameter u at the same time D ( u ) is simulated.

Felisa J. V´ azquez-Abad and Daniel Dufresne 7 Functional Estimation Functional Estimation: We performed functional estimation for 10 values of u ∈ [0 , 1], with r = 0 . 05 , σ 2 = 0 . 2 , S 0 = 50 , T = 1 . 0 and M = 10 , 000. It took 20-28 seconds for each K , both for D 2 and D 4 . 70 1 60 0.8 50 0.6 40 30 0.4 20 0.2 10 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 Estimated Variance of D 2 and D 4 as a function of u . � � Notice the difference in scale of V a r ( D 2 ) compared with V a r ( D 4 ). Conjecture: The variance is locally convex and has a unique minimum, which depends on K/S 0 . The optimal values of u for D 2 and D 4 seem to be the same (or very close).

Felisa J. V´ azquez-Abad and Daniel Dufresne 8 Optimal values of u ∗ Optimal values of u ∗ found by pilot simulations with Functional Estimation ( ≈ 24 secs CPU). Optimal Values for Change of Measure: u ∗ Method K= 30 K= 45 K=50 K= 55 K = 75 D 2 0.25 0.40 0.50 0.60 0.80 D 3 0.07 0.07 0.07 0.07 0.07 D 4 0.25 0.40 0.50 0.60 0.80 Summarizing: The behaviour of the 20 variances of the estimators: • r D 2 ( u ) and r D 4 ( u ) are similar, V a V a 15 but D 4 ( u ) has variances several orders 10 of magnitud smaller. • For D 2 and D 4 the optimal u ∗ seems to 5 be the same and it depends on K/S 0 . • D 3 does not reduce the variance consid- -0.05 0.05 0.1 0.15 0.2 erably, but u ∗ seems to be independent Estimated Variance of D 3 of K/S 0 .

Felisa J. V´ azquez-Abad and Daniel Dufresne 9 Comparison of the Methods D 3 = L u Y u α ∗ ( E Y r 2 − Y r Y u D 1 = Y r α ∗ ( E Y r 2 − Y r 1 + ˆ 2 ) , 1 + ˆ 2 ) , 1 = ( A u − K ) + D 4 = L u Y u α ∗ ( E Y r 2 − L u Y u 2 ) , Y u D 0 = Y r D 2 = L u Y u 1 , 1 , 1 + ˆ 2 = ( G u − K ) + α ∗ . D 1 , D 3 , D 4 : Control Variable methods, estimated parameter ˆ D 2 , D 3 and D 4 : Likelihood Ratio estimators, estimated optimal value u ∗ . r = 0 . 05 , σ 2 = 0 . 2 , S 0 = 50 , T = 1 . 0 and M = 10 , 000 Estimators Method K= 30 K= 45 K=50 K= 55 K = 75 D 0 20.46 ± 0.26 8.45 ± 0.216 5.80 ± 0.189 3.83 ± 0.160 0.630 ± 0.068 D 2 20.34 ± 0.137 8.32 ± 0.115 5.66 ± 0.096 3.74 ± 0.075 0.583 ± 0.020 D 1 20.31 ± 0.016 8.28 ± 0.013 5.64 ± 0.012 3.72 ± 0.011 0.585 ± 0.010 D 3 20.31 ± 0.014 8.28 ± 0.011 5.64 ± 0.010 3.71 ± 0.010 0.583 ± 0.009 D 4 20.31 ± 0.014 8.27 ± 0.009 5.62 ± 0.008 3.70 ± 0.006 0.573 ± 0.003 Variance Method K= 30 K= 45 K=50 K= 55 K = 75 D 0 176.09 121.70 92.58 66.28 12.04 D 2 49.04 34.59 23.76 14.95 1.07 D 1 0.64 0.42 0.36 0.33 0.25 D 3 0.48 0.28 0.25 0.24 0.23 D 4 0.53 0.207 0.150 0.095 0.028 α ∗ ˆ 0.998 1.05 1.07 1.10 1.20