Transcript Higher Order Finite Difference Scheme for solving 3D

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
dkm
 x k   e x km
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 ,
i1
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!