Report

Higher Order Finite Difference Scheme for solving 3D Black-Scholes equation based on Generic Factored Approximate Sparse Inverse Preconditioning using Reordering Schemes E-N.G. Grylonakis, C.K. Filelis-Papadopoulos, G.A. Gravvanis Introduction Financial Instruments Derivatives Futures Forwards Swaps Options Options Financial contracts that give the holder the right but not the obligation to buy (call option) or sell (put option) an underlying asset for a fixed price at a specific date. Problem How much money should one pay to buy a specific option contract? Topic of Interest Accurate option pricing for three underlying assets, using the multiasset Black-Scholes equation. Options Pricing Methods Binomial Options Pricing Model Lattice Methods Monte Carlo Methods Black-ScholesMerton Model Black-Scholes (BS) Equation ∂u + Lu = 0 ∂t Lu = 1 2 n p 2 i, j = 1 ij s i s jS i S j ∂u ∂S i ∂S j n + r Si i=1 ∂u ∂S i Time Dependent Convection-Diffusion-Reaction Partial Differential Equation - ru Pricing with the BS equation Single-Asset Option 1D BS PDE Closed-Form Solutions Multi-Asset Option N-D BS PDE Approximate Solutions Pricing Methodology Option Contract ① ② ③ ④ Number of underlying assets Parameters (Strike Price, Expiration date, etc.) Payoff Function (Initial Condition) Boundary Conditions Numerical solution of the corresponding BS Partial Differential Equation Three-Asset Basket Option Payoff Function Max { w[I(T)-K], 0 } I t = w j I j t n j= 1 where wj is the total investment in asset j ( as a percentage) and Ij(t) is the price of j-th asset. Linear Boundary Conditions ¶2 u(0, y,z,t) ¶2 u(S1max , y,z,t) max max = = 0, 0 £ y £ S ,0 £ z £ S 2 3 ¶x 2 ¶x 2 ¶2 u(x,0,z,t) ¶2 u(x,Smax max max 2 ,z,t) = = 0, 0 £ x £ S ,0 £ z £ S 1 3 2 2 ¶y ¶y ¶2 u(x,y,0,t) ¶2 u(x,y,S3max ,t) max max = = 0, 0 £ x £ S ,0 £ y £ S 1 2 ¶z 2 ¶z 2 Commonly used in practical pricing problems, providing stability when used with the Finite Difference Method Spatial Discretization Finite Difference Schemes (4rth order accuracy) u i-2 - 8u i-1 +8u i+1 - u i+2 u' » 12h or -u i-2 +16u i-1 - 30u i +16u i+1 - u i+2 u'' » or 2 12h 1 12 1 12 2 2 0 3 4 3 3 5 4 2 3 1 12 1 12 Ghost Values Treatment computational domain boundary ghost values Richardson’s extrapolation method (4rth order accuracy) Modified Stencils First Derivatives: Second Derivatives: 1 3 [1 1 1 2 - 2 1] 1 6 The imposition of linear boundary conditions forces the second derivatives to vanish on the boundary. The first order derivatives were discretized by a fourth order one-sided approximation: Leftmost boundary: 25 12 4 3 4 3 1 4 1 We denote by the discretized first order derivative x k for coordinate xk . Then, the stencil of the derivative with respect to coordinate k can be formed in a d-dimensional way: d 1 k 1 m k I d k m m 1 Ik m x k x k d 1 The cross-derivative can be approximated by the following expression: d 1 1 k 1 Ik m Im m I d m m k 1 m 1 x k x x x k 2 d 1 1 The coefficient matrix is then formed by the following tensor product: d -1 where: X = [ X1 , X 2 ,..., X d ] X k e x m k k -1 dkm x k e x km m 1 The above schemes reduce the programming effort substantially while providing a compact method to discretize PDE’s in higher dimensions. Numerical Time Integration After the spatial discretization, a system of Ordinary Differential Equations of the following form, occurs: du + Au = 0 dt u(x,y,z,T) = u 0 This system can be solved by the implicit fourth order backward difference scheme (BDF4): 4 n 2 1 n 3 25 n 1 n n 1 I tA u 4 u 3u u u 3 4 12 It can be observed that the BDF4 scheme requires the discrete solution in three previous time steps. These values can be obtained by the Implicit Runge-Kutta method (4rth order accurate): s y n 1 y n t b ik i , i1 s ki = f (x n + t c i , y n + t a ij k j ) , j= 1 The coefficients of the 2-stage method are: 1 3- 2 3 3+ 2 3 1 a11 = , a12 = , a 21 = , a 22 = , 4 12 12 4 1 3- 3 3+ 3 b1 = b2 = , c1 = , c2 = 2 6 6 Solving the Linear System The arising large, sparse,linear system was solved by the Preconditioned BiConjugate Gradient Stabilized (PBiCG-STAB) method, in conjunction with the Modified Generic Factored Approximate Sparse Inverse (MGenFAspI) scheme. MGenFAspI matrix: M=GH The MGenFAspI matrix is computed by solving the following systems: LH lfill droptol =I UG lfill droptol =I The modified approach minimizes the searches for elements and enhances the performance of the method. Approximate Minimum Degree (AMD) Reordering When attempting to solve large sparse linear systems, reordering schemes can be used in order to minimize the fill-in during the factorization process. The AMD algorithm produces a reordering such that the vertices with minimum degree are to be eliminated first. The degree of each vertex is approximated through an upper bound created by the sum of the weights of the neighboring vertices, increasing the performance of the resulting ordering scheme. Implementation Issues In order to compute the three initial solutions with the Runge-Kutta method, the solution of four linear systems at every time step is required. Recalling the vectors, required by the R-K method: s ki = f (x n + t c i , y n + t a ij k j ) , j= 1 The 2-stage method requires the computation of vectors k1 and k2 which can be obtained by the following system: Dt 2 3 2 1 (I A)k1 - Dt ( )Ak2 = DtAu 4 4 6 3 Dt 2 2 1 -Dt ( + )Ak1 + (I A)k2 = DtAu 4 6 4 The above system can be expressed in the following block form: A1 C 1 t Au n n D 1 k 2 t Au B1 k 1 n t Au A 1 B1 k1 t Au n C A 1 t Au n 0 k S 2 1 1 where: S = ( D1 - C1A1-1B1 ) The computation of k1 and k2 is then performed by solving the following linear systems at every time step: ( n) -1 1 1 Sk2 = DtAu - C A (DtAu ) ( n) ( n) A1k1 = DtAu - B1k2 The Schur complement is computed implicitly, since iterative methods do not require the coefficient matrix explicitly, because the product of a matrix by a vector is only needed. Thus: Sx = ( D1 - C A B ) x = ( D1x - C1y), y = A ( B1x ) -1 1 1 1 A1y = ( B1 x ) -1 1 Numerical Results The estimated price of the basket option: Performance (“seconds.hundreds”) of the PBiCG-STAB, based on the MGenFAspI in conjunction with AMD reordering scheme, for various values of N and droptol: Convergence behavior of the PBiCG-STAB, based on the MGenFAspI in conjunction with AMD reordering scheme, for various values of N and droptol: The number of nonzero elements in the G and H factors of the MGenFAspI for various values of N and droptol: Conclusions 1. The Black-Scholes PDE can be used to price options with many underlying assets, without relying solely on Monte Carlo methods. 2. The high order schemes combined with a multi-dimensional PDE result in a large, sparse, linear system, thus, iterative methods are the best choice. 3. Preconditioners and reordering schemes can be used to enhance the performance of the chosen iterative method. 4. The MGenFAspI matrix has been proved to be an effective preconditioner and combined with various iterative methods has achieved better convergence behavior in comparison with other methods. 5. Moreover, the applicability of the MGenFAspI matrix in conjunction with the PBiCG-STAB method has been evaluated, for various model problems, derived from Computational Fluid Dynamics, Computational Structural Analysis and Plasma Physics. References Achdou, Y., and Pironneau, O. 2005. Computational methods for option pricing. SIAM. Amestoy, P., Davis, T.A., and Duff, I.S. 1996. An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications 17(4), 886-905. Chow, E. 2000. A priori sparsity patterns for parallel sparse approximate inverse preconditioners. SIAM J. Sci. Comput. 21, 1804-1822. Duffy, D.J. 2006. Finite difference methods in financial engineering: A partial differential equation approach. John Wiley and Sons. Filelis-Papadopoulos, C.K., and Gravvanis, G.A. A class of generic factored and multilevel recursive approximate inverse techniques for solving general sparse systems (submitted). Grylonakis E-N.G. 2014. On the study and numerical solution of the Black-Scholes equation. Dissertation Thesis, Department of Electrical and Computer Engineering, Democritus University of Thrace. Gustafsson, B. 2008. High order difference methods for time dependent PDE. Springer Series in Computational Mathematics, Vol.38 Haug, E.G. 2007. The complete guide to option pricing formulas. McGraw-Hill. Hull, J. 2009. Options, futures and other derivatives. Prentice-Hall. Jeong, D., Kim, J., and Wee, I-S. 2009. An accurate and efficient numerical method for BlackScholes equations. Commun. Korean Math. Soc. 24(4), 617–628 Leentvaar, C.C.W, and Oosterlee, C.W. 2008. On coordinate transformation and grid stretching for sparse grid pricing of basket options. J. Comp. Appl. Math. 222, 193-209. Persson, J., and von Sydow, L. 2007. Pricing European multi-asset options using a space-time adaptive FD-method. Comput. Vis. Sci. 10, 173–183. Saad, Y., Schultz, M.H. 1986. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Comp. 7, 856-869. Süli, E., and Mayers, D.F. 2003. An introduction to numerical analysis. Cambridge University Press. Trottenberg, U., Osterlee, C.W., and Schuller, A. 2000. Multigrid. Academic Press. Van der Vorst, H.A. 1992. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. and Stat. Comput. 13(2), 631–644. Wilmott, P., Dewynne, J., and Howison, S. 1994. Option pricing: Mathematical models and computation. Oxford Financial Press. Zhang, P.G. 1998. Exotic Options: A guide to second-generation options. World Scientific, 2nd edition. Thank you for your attention!