A Flexible Galerkin Scheme for Option Pricing in Lévy Models
aa r X i v : . [ q -f i n . C P ] M a r A Flexible Galerkin Scheme for OptionPricing in L´evy Models
Maximilian Gaß and Kathrin GlauSeptember 19, 2018Technische Universit¨at M¨unchen, Center for [email protected], [email protected]
Abstract
One popular approach to option pricing in L´evy models is through solvingthe related partial integro differential equation (PIDE). For the numerical so-lution of such equations powerful Galerkin methods have been put forward e.g.by Hilber et al. (2013). As in practice large classes of models are maintainedsimultaneously, flexibility in the driving L´evy model is crucial for the implemen-tation of these powerful tools. In this article we provide such a flexible finiteelement Galerkin method. To this end we exploit the Fourier representation ofthe infinitesimal generator, i.e. the related symbol, which is explicitly availablefor the most relevant L´evy models. Empirical studies for the Merton, NIG andCGMY model confirm the numerical feasibility of the method.
Keywords
L´evy processes, partial integro differential equations, pseudo-differentialoperators, symbol, option pricing, Galerkin approach, finite element method
Mathematics Subject Classification (2000)
In computational finance, methods to solve partial differential equations come intoplay, when both run time and accuracy matters. In contrast to Monte Carlo for ex-ample, run-time is very appealing and a deterministic and conservative error analysisis established and well understood. And, compared to Fourier methods, the possi-bility to capture path-dependent features like early exercise and barriers is naturallybuilt in. Within these appealing features lies the capacity to attract interest fromacademia and satisfy the needs of the financial industry alike.In academia, a series of publications by Cont and Voltchkova (2005b), Hilber,Reich, Schwab and Winter (2009), Salmi, Toivanen and Sydow (2014), Itkin (2015)and Glau (2016b) and the monograph of Hilber, Reichmann, Schwab and Winter(2013) have opened the theory to include even more sophisticated models of L´evytype, resulting in Partial
Integro
Differential Equations (PIDEs). The theoreticalresults have been validated by sophisticated numerical studies. In this context,Schwab and his working group in particular have taken the lead and unveiled the1otential of PIDE theory for practical purposes in the financial industry. Combiningstate of the art compression techniques with a wavelet basis finite element setup hasresulted in a numerical framework for option pricing in advanced and multivariatejump models and thereby moved academic boundaries.In the financial industry an awareness of the full potential of these tools is yetto be developed. Advocating the advancement of numerical methods one mustacknowledge what practice cherishes most in them. Due to model uncertainty andbehavioural characteristics of different portfolios, financial institutions need to dealwith a number of different pricing models in parallel. Or, in the words of F¨ollmer(2009): ”In any case, the signal towards the practitioners of risk management isclear: do not commit yourself to a single model, remain flexible, vary the modelsin accordance with the problem at hand, always keeping in mind the worst casescenario.” These features need to be reflected in the numerical environment.In this article, we aim at reconciling capacities of state of the art P(I)DE toolswith the flexibility regarding different model choices as required by industry. Desir-able features that such an implementation must offer include(1) a degree of accuracy that reaches levels relevant to practical applications andthat can be measured and controlled by a theoretical error analysis,(2) fast run times,(3) low and feasible implementational and maintenance cost,(4) a flexibility of the toolbox towards different options and models.Two standard methods are available for solving P(I)DEs, that is the finite differenceapproach and the finite element method. More recently, also radial basis methodshave been pushed forward to solve pricing PIDEs. In principle these concepts canbe implemented in such a way that they achieve the desired features 1.–3. and im-plementations for a variety of models and option types have already been developed:Finite difference schemes solving PIDEs for pricing European and barrier optionswith an implementation for Merton and Variance Gamma are provided in Cont andVoltchkova (2005b) and Cont and Voltchkova (2005a). The method has been furtherdeveloped in different directions, we mention one example, Itkin and Carr (2012),who use a special representation of the equation to derive a finite different schemefor jump diffusions with jump intensity of tempered stable type. Wavelet-Galerkinmethods for PIDEs related to a class generalizing tempered stable L´evy processesare derived in Matache et al. (2005) for American options and see e.g. Marazzinaet al. (2012) for a high-dimensional extension. Radial basis for the Merton and Koumodel, American and European options are provided by Chan and Hubbert (2014)and further developed for CGMY models by Brummelhuis and Chan (2014).An implementation that is flexible in the driving model as well as in the optiontype first of all requires a problem formulation covering the collectivity of envisagedmodels and options. In view of feature (1), a unified approach to the error analysisof the resulting schemes is of equal importance. Galerkin methods, accruing fromthe Hilbert space formulation of the Kolmogorov equation, seem to be predestinedto deliver the adequate level of abstraction for this task. It is this abstract level Tranlsated from German. the development of a model-independent approach toset up a FEM solver for option pricing in L´evy models that we call symbol method .We address this goal by expressing the operator in the Fourier space. This meansaccessing the model specific information via the symbol. In contrast to the operator,the symbol is explicitly available for a variety of models and is thus numericallyaccessible. Further advantages will be highlighted in subsequent sections.Section 2 introduces the theoretical framework for our PIDEs of interest andtheir weak formulation. The next section describes the solution scheme, that is theGalerkin approximation in space. We investigate the scheme with regard to thenumerical challenges arising during its implementation. Section 4 introduces thesymbol method itself. All components of the FEM solver are expressed in Fourierspace. The subsequent numerical evaluation of the stiffness matrix entries is sup-ported by an elementary approximation result. Several examples of symbols forwell-known L´evy models confirm the wide applicability of the method and its nu-merical advantages. Two proposals for the implementation of basis functions arepresented. The numerical studies in Section 5 confirm theoretical prescribed ratesof convergence and validate the claim of numerical feasibility.
We first introduce the underlying stochastic processes, the Kolmogorov equation, itsweak formulation as well as the solution spaces of our choice.
Let a stochastic basis (Ω , F , ( F t ) ≤ t ≤ T , P ) be given and let L be an R d -valued L´evyprocess with characteristics ( b, σ, F ; h ), i.e. for fixed t ≥ E e i h ξ,L t i = e − tA ( − iξ ) d s for every ξ ∈ R d , (1)where the symbol of the process is defined as A ( ξ ) := 12 h ξ, σξ i + i h ξ, b i − Z R d (cid:16) e − i h ξ,y i − i h ξ, h ( y ) i (cid:17) F (d y ) . (2)Here, σ is a symmetric, positive semi-definite d × d -matrix, b ∈ R d , and F is aL´evy measure, i.e. a positive Borel measure on R d with F ( { } ) = 0 and R R d ( | x | ∧ F (d x ) < ∞ . Moreover, h is a truncation function i.e. h : R d → R such that R {| x | > } h ( x ) F (d x ) < ∞ with h ( x ) = x in a neighbourhood of 0. The Kolmogorovoperator of a L´evy process L with characteristics ( b, σ, F ; h ) is given by A ϕ ( x ) := − d X j,k =1 σ j,k ∂ ϕ∂x j ∂x k ( x ) − d X j =1 b j ∂ϕ∂x j ( x ) − Z R d (cid:16) ϕ ( x + y ) − ϕ ( x ) − d X j =1 ∂ϕ∂x j ( x ) h j ( y ) (cid:17) F (d y ) (3)for every ϕ ∈ C ∞ ( R d ), where h j denotes the j -th component of the truncationfunction h . Key for the variational formulation of Kolmogorov equation ∂ t u + A u = f (4) u (0) = g (5)is the definition of the bilinear form a ( ϕ, ψ ) := Z R d ( A ϕ )( x ) ψ ( x ) d x for all ϕ, ψ ∈ C ∞ ( R d ) . (6)It is one of the major advantages of variational formulations of evolution equationsthat solution spaces of low regularity, as compared to the space C for example,are incorporated in an elegant way. Departing from the space C ∞ ( R d ) of smoothfunctions with compact support, we can select from a large variety of function spaces V that is characterized by the following assumption.(A1) V and H are Hilbert spaces such that C ∞ ( R d ) lies dense in V and there existsa continuous embedding from V into H .Existence and uniqueness of a variational solution critically hinges on the followingtwo properties of the bilinear form:(A2) Continuity : There exists a constant
C > (cid:12)(cid:12) a ( ϕ, ψ ) (cid:12)(cid:12) ≤ C k ϕ k V k ψ k V for all ϕ, ψ ∈ C ∞ ( R d ) . G˚arding inequality : There exists constants
G > G ′ ≥ a ( ϕ, ϕ ) ≥ G k ϕ k V − G ′ k ϕ k H for all ϕ ∈ C ∞ ( R d ) . We observe that due to (A1) and (A2), the bilinear form a possesses a uniquecontinuous bilinear extension a : V × V that is continuous, i.e. for a constant C > (cid:12)(cid:12) a ( ϕ, ψ ) (cid:12)(cid:12) ≤ C k ϕ k V k ψ k V for all v ∈ V . Also (A3) holds for all v ∈ V .As V is separable, this is also true for H and one can find a continuous embeddingfrom H to the dual space V ∗ of V , i.e. ( V, H, V ∗ ) is a Gelfand triplet. We thendenote by L (cid:0) , T ; H (cid:1) the space of weakly measurable functions u : [0 , T ] → H with R T k u ( t ) k H d t < ∞ and by ∂ t u the derivative of u with respect to time in thedistributional sense. For a detailed definition which relies on the Bochner integralwe refer to Section 24.2 in Wloka (1987). The Sobolev space W (0 , T ; V, H ) := n u ∈ L (cid:0) , T ; V (cid:1) (cid:12)(cid:12)(cid:12) ∂ t u ∈ L (cid:0) , T ; V ∗ (cid:1)o , (7)will play the role of the solution space in the variational formulation of Kolmogorovequation (4), (5). Definition 2.1.
Let f ∈ L (cid:0) , T ; V ∗ (cid:1) and g ∈ H . Then u ∈ W (0 , T ; V, H ) is a variational solution of Kolmogorov equation (4) , if for almost every t ∈ (0 , T ) , h ∂ t u ( t ) , v i H + a ( u ( t ) , v ) = h f ( t ) | v i V ∗ × V for all v ∈ V (8) and u ( t ) converges to g for t ↓ in the norm of H . Remark 2.2.
Assumptions (A1)–(A3) guarantee the existence and uniqueness ofa variational solution u ∈ W (0 , T ; V, H ) of (8) , see for instance Theorem 23.A inZeidler (1990). Definition 6 is based on the L -scalar product and is appropriate for variationalequations in Sobolev spaces. Then, typically H = L . For Kolmogorov equationsfor option prices the initial condition g in (5) plays the role of the (logarithmicallytransformed) payoff function of the option. For a call option with strike K it is ofthe form x ( S e x − K ) + , for a digital up and out option it is given by x x
1. For functions in H α with α < H , also in the more challenging case of solution spaceswith fractional order derivatives we derive the validity of the representation underappropriate assumptions. Lemma 3.1.
Let d = 1 . Let a be defined by (9) . Assume (A1)–(A3) for a , V and H and denote by a : V × V its unique bilinear continuous extension. If H η ( R ) ⊂ V ,we have for every ϕ, ψ ∈ H η ( R ) , a ( ϕ, ψ ) = σ Z R ϕ ′ ( x ) ψ ′ ( x ) e ηx d x − b ( η, σ, F ) Z R ϕ ′ ( x ) ψ ( x ) e ηx d x − Z R Z | y | < Z y Z z ϕ ′ ( x + v ) d v d zF (d y ) (cid:0) ψ ′ ( x ) + 2 ηψ ( x ) (cid:1) e h η,x i d x (17) − Z R Z | y | > (cid:0) ϕ ( x + y ) − ϕ ( x ) (cid:1) F (d y ) ψ ( x ) e h η,x i d x ith b ( η, σ, F ) = b − ση + Z | y | < (cid:0) y − h ( y ) (cid:1) F (d y ) − Z | y | > h ( y ) F (d y ) . The proof is provided in Section A.1.Inspecting the expression for the bilinear form we encounter several numericalchallenges due to the integral part—stemming from the jumps of the process:1. The appealing tridiagonal structure of the stiffness matrix related to the Black-Scholes equation does not extend to the general L´evy setting. Instead, thestiffness matrix is densely populated. Pleasantly, it is still a Toeplitz matrix.2. For some choices of L´evy measures and bases the stiffness matrix entries may bederived in closed form. This is for instance the case for the Merton model andpiecewise linear basis functions. Following Section 10.6.2 in Hilber et al. (2013),the stiffness matrix entries may be derived in semi-closed form expressions fora further group of jump intensities including tempered stable, CGMY andKoBoL processes and the choice of piecewise linear basis functions.An implementation that is flexible in the driving L´evy process therefore hasto rely on numerical approximations of the entries of the stiffness matrix. Theseapproximations inevitably affect the accuracy of the solution to the scheme (13)–(14). The question rises:
How accurate does the integration routine have to bechosen in order to meet a desired accuracy of the solution V ? In order to gain afirst practical insight in the magnitude of the error resulting from an inaccuracy inthe stiffness matrix entries, we next conduct an empirical investigation.
To this end we conduct an accuracy study for the stiffness matrix entries simulatingthe propagation of integration errors in the stiffness matrix up until computed optionprices. We choose the Black-Scholes model for which the precise values of all matrixentries are known and assume a market volatility of σ = 0 . r = 0 .
01 for this study. With a classic FEM solver we price a put option with strike K = 1 and maturities T ∈ [0 ,
3] for current values of the stock S ∈ [ S min , S max ] with S min = 0 .
01 and S max = 10. We set the number of involved FEM hat functions to N = 150, resulting in a mesh with 150 inner grid nodes and mesh fineness h = 0 . M = h · · ·
01 4 1 . . . ...0 . . . . . . . . . 0... . . . 1 4 10 · · · ∈ R N × N , and the stiffness matrix to be given by A = A bs = A (1) + A (2) + rM ∈ R N × N , (18)8here A (1) = 12 (cid:18) r − σ (cid:19) − · · ·
01 0 − − · · · , A (2) = σ h − · · · − − − − · · · − . With these matrices we set up a theta scheme, θ = 0 .
5, and derive Black-Scholesput option prices. We use these matrices to simulate how the resulting pricingsurface is affected by inaccuracies that might occur when these integrals are solvednumerically, instead. To this extent we take the correct stiffness matrix given by (18)and distort each of its entries randomly at different positions D ∈ N after thedecimal point by adding ε Di = 10 − ( D − ε i with random ε i ∈ ( − ,
1) for i ∈ {− ( N − , . . . , − , , , . . . , ( N − } onto the (side) diagonal i of Matrix A . Each individual(side) diagonal of the original stiffness matrix is thus affected evenly, keeping theToeplitz structure of the matrix intact. Since the value of A ij is only determined bythe value of j − i , this distortion mimics the influence that integration inaccuracieswould have.Hence, for D ∈ N we define the distorted stiffness matrix by A D distort = A + ε D ∈ R N × N , (19)with ε D = 10 − ( D − ε ε ε · · · · · · · · · ε N − ε − ε ε . . . . . . . . . ... ε − ε − . . . . . . . . . . . . ...... . . . . . . . . . . . . . . . ...... . . . . . . . . . . . . ε ε ... . . . . . . . . . ε − ε ε ε − ( N − · · · · · · · · · ε − ε − ε ∈ R N × N , with uniformly distributed ε i ∈ ( − , i ∈ {− ( N − , . . . , − , , , . . . , ( N − } ,that are drawn independently from each other. Using these distorted stiffness ma-trices A D distort for different values D ∈ N , we derive again price surfaces of the putoption in the Black-Scholes model and compare the difference between the pricescoming from the distorted stiffness matrix A D distort ∈ R N × N to the prices from theintact stiffness matrix A ∈ R N × N .Figure 1 shows the results of the study and emphasizes the necessity for highnumerical accuracy in the computation of the stiffness matrix. We observe thatthe absolute price differences decrease almost linearly in D . An accuracy of D = 3corresponds to integration results that are exact up to the third digit after thedecimal point. Pricing resulting from stiffness matrices computed with such a lowintegration accuracy are unacceptable. The respective pricing errors observable inthe top left corner of Figure 1 indicate relative errors of several hundred percentpoints. With more precise integration results, the error decreases in D until highly9 D = 3 S -5-10-151050 a b s . p r i ce d i ff . T D = 5 S -0.0100.010.02 a b s . p r i ce d i ff . T D = 7 S -5051015 × -5 a b s . p r i ce d i ff . T D = 10 S × -7 a b s . p r i ce d i ff . Figure 1: Absolute price differences resulting from a distortion of the stiffness matrix A . True and distorted prices describe the market value of a put option in the Black-Scholes model with the parametrization from Section 3.3. We compare the pricesurfaces coming from a theta scheme using the stiffness matrix A given by (18) tothe respective pricing surface when A is replaced by A D distort , the distorted versionof A as defined in (19), for different values of D ∈ N . The influence of the distortiondecreases in D .appealing pricing results are achieved for D = 7 and beyond. The magnitude of thepricing error resulting from a distorted stiffness matrix emphasizes the necessity ofbeing able to derive the stiffness matrix entries as accurately as possible.10 Fourier approach to the Kolmogorov equation
In regard to the high accuracy the approximation of the stiffness matrix entries needsto achieve, we would like to avoid numerical evaluations of the stiffness matrix entrieson the basis of representation (17). Seeking for alternative representations of thestiffness matrix, let us point out that the symbol A of the L´evy process is alwaysavailable. Even more, it is an explicit function of the parameter of the process andthus can be seen as the modelling quantity of the process as the examples 4.6–4.9below show. We therefore take a Fourier perspective on the variational formulationof the Kolmogorov equation. This is especially promising since the Kolmogorovoperator A of a L´evy process is a pseudo differential operator with with symbol A , A ϕ = F − ( A F ( ϕ )) for all ϕ ∈ C ∞ ( R d ) , (20)as elementary manipulations show. Now Parseval’s identity yields a ( ϕ, ψ ) = 1(2 π ) d Z R d F ( A ϕ )( ξ ) F ( ψ )( ξ ) d ξ (21)for all ϕ, ψ ∈ C ∞ ( R d ), respectively, a ( ϕ, ψ ) = 1(2 π ) d Z R d A ( ξ ) F ( ϕ )( ξ ) F ( ψ )( ξ ) d ξ. (22)This well-known identity has already proven highly beneficial for the analysis ofthe variational solutions of the Komogorov equations, compare e.g. Hilber et al.(2013), Glau (2016a) and Glau (2016b). Here we exploit this representation for thenumerical implementation. Lemma 4.1 (Continuous extension of bilinear forms) . Let A be the symbol of a L´evyprocess given by the characteristic triplet ( b, σ, F ) . Denote by A : C ∞ ( R d , C ) → C ∞ ( R d , C ) the pseudodifferential operator associated with symbol A . Furthermore,denote by a : C ∞ × C ∞ → C the bilinear form associated with the operator A . Let η ∈ R d . Ifi) the exponential moment condition Z | x | > e −h η ′ ,x i F (d x ) < ∞ (23) holds for all η ′ ∈ sgn( η )[0 , | η | ] × · · · × sgn( η d )[0 , | η d | ] andii) there exists a constant C > with | A ( z ) | ≤ C (1 + k z k ) α (24) for all z ∈ U − η where U − η = U − η × · · · × U − η d (25) with U − η j = R − i sgn( η j )[0 , | η j | ) , hen a ( · , · ) possesses a unique linear extension a : H α/ η × H α/ η → C which can bewritten as a ( ϕ, ψ ) = 1(2 π ) d Z R d A ( ξ − iη ) b ϕ ( ξ − iη ) b ψ ( ξ − iη ) d ξ (26) for all ϕ, ψ ∈ H α/ η ( R d ) .Proof. The proof can be found in Eberlein and Glau (2011) using Theorem 4.1therein and Parseval’s identity.In order to gain first insight in the convergence analysis, we fix a level N inthe Galerkin scheme and derive conditions for the convergence of the sequence ofweak solutions that we obtain by approximating the stiffness matrix entries. Inthe implementation in Section 5 below we will also approximately compute the righthand side of the equation. We therefore more generally consider sequences of stiffnessmatrices, right hand sides and initial conditions.As usual, we denote for a given bilinear form a : V × V → R the associatedoperator A : V → V ∗ defined by A ( u )( v ) := a ( u, v ) for all u, v ∈ V . Lemma 4.2.
Let V , H and a : V × V → R satisfy (A1)–(A3). Let X :=span { ϕ , . . . , ϕ N } ⊂ V and for each n ∈ N let(An1) f n , f ∈ L (0 , T ; H ) with f n → f in L (cid:0) , T ; X ∗ ) ,(An2) g n , g ∈ H with g n → g in H ,(An3) a n : V × V → R be a bilinear form such that for all j, k ≤ N , (cid:12)(cid:12) ( a n − a )( ϕ j , ϕ k ) (cid:12)(cid:12) → . (27) Then the sequence of unique weak solutions v n ∈ W (0 , T ; X, H ) of ˙ v n + A n v n = f n , v n (0) = g n (28) converges strongly in L (cid:0) , T ; X ) ∩ C (0 , T ; H ) to the unique weak solution v ∈ W (0 , T ; X, H ) of ˙ v + A v = f, v (0) = g. (29)The proof is provided in Section A.2.Profound convergence analysis providing asymptotic convergence rates for thefully discrete scheme can be derived based on the techniques presented by von Pe-tersdorff and Schwab (2003), which is beyond the scope of the present article. The key component of a Galerkin FEM solver is the model dependent stiffnessmatrix. Using expression (16) of section 3.2 above, the entries of that matrix canbe derived. The existence of the L´evy measure F in that expression, however,renders the numerical derivation of the matrix rather cumbersome. Additionally,the empirical accuracy study of section 3.3 emphasize that utmost care must betaken when the stiffness matrix entries are numerically derived. Consequently, inthis section we approach the calculation of that FEM solver components differently.12he Fourier approach indicated by Lemma 4.1 will allow us to access the modelinformation required for the stiffness matrix and all other FEM solver componentsnot via the operator but rather via the associated symbol, instead. In stark contrastto the operator, the symbol of a L´evy model is numerically accessible in a unifiedway for a large set of underlying models and we will present several examples.Let us state the core lemma of this section. Here we concentrate on basis func-tions obeying a simple nodal translation property, which is in particular satisfiedfor the classical piecewise polynomial basis functions. The translation property canalso easily been extended to that one satisfied by biorthogonal wavelet bases thathaven prove to be useful for solving Komlogorov PIDEs for option pricing, compareMatache et al. (2004). Lemma 4.3 (Symbol method for bilinear forms) . Let the assumptions of Lemma 4.1be satisfied with η = 0 . Assume further for N ∈ N a set of functions ϕ , ϕ , . . . , ϕ N ∈ H α/ ( R ) and nodes x , . . . , x N ∈ R , such that for all j = 1 , . . . , Nϕ j ( x ) = ϕ ( x − x j ) , ∀ x ∈ R , holds. Then we have a ( ϕ l , ϕ k ) = 12 π Z R A ( ξ ) e iξ ( x l − x k ) | c ϕ ( ξ ) | d ξ. (30) for all k, l = 1 , . . . , N . If additionally ℜ ( A ( ξ )) = ℜ ( A ( − ξ )) and ℑ ( A ( ξ )) = −ℑ ( A ( − ξ )) , (31) then a ( ϕ l , ϕ k ) = 1 π Z ∞ ℜ (cid:16) A ( ξ ) e iξ ( x l − x k ) (cid:17) | c ϕ ( ξ ) | d ξ (32) for all k, l = 1 , . . . , N .Proof. Elementary properties of the Fourier transform yield c ϕ j ( ξ ) = e iξx j c ϕ ( ξ ) , ∀ ξ ∈ R , (33)where c ϕ ( ξ ) = 2 ξ h (1 − cos( ξh )) , ∀ ξ ∈ R . (34)Since ϕ j ∈ H α/ ( R ), for all j = 1 , . . . , N , the identity (30) follows from Lemma 4.3above or Theorem 4.1 and Remark 5.2 and the lines thereafter in Eberlein and Glau(2011), respectively. The second claim (32) is then elementary. Corollary 4.4 (Symbol method for stiffness matrices) . Let A ∈ S α be a univariatesymbol with associated operator A . Denote by ϕ j ∈ L ( R ) , j ∈ , . . . , N the basisfunctions of a Galerkin scheme associated with an equidistantly spaced grid Ω = { x , . . . , x N } possessing the property ϕ j ( x ) = ϕ ( x − x j ) , ∀ x ∈ R , (35) for some ϕ : R → R . Then, the stiffness matrix A ∈ R N × N of the scheme can becomputed by A kl = 12 π Z R A ( ξ ) e iξ ( x l − x k ) | c ϕ ( ξ ) | d ξ (36) for all k, l = 1 , . . . , N . roof. The proof is an immediate consequence of Lemma 4.3.
Remark 4.5 (On the symbol method for bilinear forms) . From a numerical per-spective, the representations of the stiffness matrix entries provided in Lemma 4.3and Corollary 4.4 are highly promising:1. Instead of the double integrals appearing in (16) , only one dimensional integralneed to be computed.2. The model specific information is expressed via the symbol ξ A ( ξ ) , whichfor a large set of models is available in form of an explicit function of ξ andthe model parameters, a feature that we now can exploit numerically. We givea short list of examples below. For further examples we refer to Glau (2016a)and Glau (2016b).3. Representation (36) displays the entries of the stiffness matrix as Fourier inte-grals . Moreover, the nodes appear as Fourier variables. As a consequence, FastFourier Transform (FFT) methods can be used to accelerate their simultaneouscomputation.
Expression (3) introduces operators A for L´evy processes L in terms of thecharacteristic triplet ( b, σ, F ). The following examples present the respective symbolsfor some well known L´evy models. Example 4.6 (Symbol in the Black-Scholes (BS) model) . In the univariate Black-Scholes model, determined by the Brownian volatility σ > , the symbol is given by A ( ξ ) = A bs ( ξ ) = iξb + 12 σ ξ , (37) with drift set b to b = r − σ . (38) Example 4.7 (Symbol in the Merton model) . In the Merton model where σ > , λ > , α ∈ R and β > , the symbol computes to A ( ξ ) = A merton ( ξ ) = 12 σ ξ + iξb − λ (cid:16) e − iαξ − β ξ − (cid:17) (39) with drift set to b = r − σ − λ (cid:18) e α + β − (cid:19) , (40) as required by the no-arbitrage condition. Example 4.8 (Symbol in the CGMY model) . In the CGMY model of Carr et al.(2002) with σ > , C > , G ≥ , M ≥ and Y ∈ (1 , , the symbol computes to A ( ξ ) = A cgmy ( ξ ) = iξb + 12 σ ξ − C Γ( − Y ) (cid:2) ( M + iξ ) Y − M Y + ( G − iξ ) Y − G Y (cid:3) , (41) for all ξ ∈ R , with drift b set to b = r − σ − C Γ( − Y ) (cid:2) ( M − Y − M Y + ( G + 1) Y − G Y (cid:3) (42) for martingale pricing. xample 4.9 (Symbol in the NIG model) . With σ > , α > , β ∈ R and δ > such that α > β , the symbol of the NIG model is given by A ( ξ ) = A nig ( ξ ) = 12 σ ξ + iξb − δ (cid:16)p α − β − p α − ( β − iξ ) (cid:17) (43) for all ξ ∈ R with drift given by b = r − σ − δ (cid:16)p α − β − p α − ( β + 1) (cid:17) (44)Implementing (36), we encounter new numerical challenges: From the perturba-tion study in Section 3.3 that we need to evaluate the integrals at high precision.Consider first the Black-Scholes model and choose the piecewise linear hat functionsas basis elements as a toy example. Applying a standard Matlab integration routinewill lead to considerable errors. To understand the effect, let us have a closer lookat the integrands in (36), which we depict in Figure 2. ξ c ϕ ( ξ ) Graph of c ϕ ξ
100 110 120 130 140 150 160 170 180 190 200 c ϕ ( ξ ) × -4 ξ c ϕ ( ξ ) × -6 Figure 2: Graph of c ϕ , the Fourier transform of the hat function ϕ centered overthe origin, evaluated over three subintervals of R + . The mesh is chosen with h = 1.The oscillations and the rather slow decay to zero complicate numerical integrationwith high accuracy requirements considerably when c ϕ is involved.More precisely, we show several integrands of A ∈ R N × N in the representationprovided by (36) of Corollary 4.4 with the symbol given by Example 4.6. Eachintegrand is evaluated for a different value of l − k over three different subintervalstaken from the unbounded integration range. In the Fourier approach of calculating15
10 20 30 40 50 60 70 80 90 100 ξ -505 × -3 Integrand of A bsij j − i = 0 j − i = 5 j − i = 150
100 110 120 130 140 150 160 170 180 190 200 ξ -505 × -4 j − i = 0 j − i = 5 j − i = 150 ξ -505 × -6 j − i = 0 j − i = 5 j − i = 150 Figure 3: The integrand for the Black-Scholes stiffness matrix A kl for several valuesof l − k . The grid of the hat functions spans the interval [ − ,
5] with 150 equidistantlyspaced inner nodes and grid fineness h = 0 . | l − k | .16he stiffness matrix A ∈ R N × N via the respective symbol, the integrands of A kl wouldhave to be numerically integrated for all l − k ∈ {− ( N − , . . . , − , , , . . . , N − } .The larger | j − i | , however, the more severe the numerical challenges for evaluatingthe integrand, as Figure 3 demonstrates. All integrands illustrated therein decayrather slowly. Additionally, oscillations increase in | j − i | . In combination, thesetwo observation seriously threaten a numerically reliable evaluation of the integral.With increasing values of | j − i | , the oscillations of the integrand accelerate and thenumber of necessary supporting points for accurate integration soars.These considerations show us that we need to further investigate the problemto obtain a flexible method to compute the stiffness matrix reliably and with lowcomputational cost. The path that we propose here is to modify the problem insuch a way that the resulting integrands decay much faster so that the domain ofintegration can be chosen considerably smaller and a usual integration routine suchas Matlab’s function quadgk is sufficient. To do so, we first observe that the hatfunctions, which we used in our toy example, are piecewise linear functions. Whilebeing continuous they are not continuously differentiable everywhere and thus lacksmoothness on an elementary level already. This lack of smoothness translates intoa slow decay of their Fourier transform.Therefore, we propose to replace the piecewise linear basis functions by basisfunctions that display considerably higher regularity leading to appealing decayproperties of the integrands in (36). In the following two sections, we present twodifferent approaches to implement such a problem modification. It is well known that convolution with a smooth function has a smoothing effect onthe function that the convolution is applied to. Our first basis function alternativewill therefore be a classic hat function smoothed by convolution. Functions thatqualify for this smoothing by convolution are called mollifiers. In order to choosean appropriate mollifier for our purposes—the fast and accurate computation of theintegrals in (36), the mollifiers need to display two essential features:(1) The Fourier transform of the modified basis function needs to be available.(2) The smoothing effect needs to be steerable through a parameter.As the Fourier transform of the convolution of two functions is the product of thetwo Fourier transformed functions, (2) boils down to the availability of the Fouriertransform of the mollifier. Since the Fourier transform of standard mollifiers isnot available in closed form, we widen the range of the standard mollifiers andallow for non-compact support. More precisely, we call the sequence m = ( m k ) k ∈ N , m k ∈ L ( R d ) for all k ∈ N , a mollifier , if1. m k ≥
0, for all k ∈ N ,2. R R d m k ( x ) d x = 1, and3. for all ̺ > R R d \ B ̺ (0) m k ( x ) d x → k → ∞ .17 -2 -1.5 -1 -0.5 0 0.5 1 1.5 200.20.40.60.81 hat function ε = 0 . ε = 0 . ε = 0 . Figure 4: A comparison between the classic hat function ϕ on a grid with h = 1 andthe mollified hat function ϕ ε = ϕ ∗ m ε Gaussian for several values of ε ∈ { . , . , . } using the Gaussian mollifier of Example 4.10.Feature (1) is often required and we follow the usual construction here. ByProposition and Definition 2.14 in Alt (2011) we can adjust the influence of mollifi-cation by a mollification ε . To this end let m ∈ L ( R d ) with m ≥ , and Z R d m ( x ) d x = 1 . (45)Define m ε = 1 ε d m (cid:16) · ε (cid:17) . (46)Then for each ̺ > R R d m ε ( x ) d x = 1 and R R d \ B ̺ (0) m ε ( x ) d x → ε →
0. Consequently, for each null sequence ( ε k ) k ∈ N the sequence ( m ε k ) k ∈ N is amollifier in the sense of our definition. Example 4.10 (A mollifier based on the Normal distribution) . We present anexample for a mollifier. Define m Gaussian ( x ) = 1 √ π e − x . (47) Then we call ( m ε k Gaussian ) k ∈ N defined according to (46) a Gaussian mollifier . Thecharacteristic function of the Gaussian mollifier is known in closed form, \ m ε Gaussian ( ξ ) = exp (cid:18) − ε ξ (cid:19) , (48) thus exhibiting exponential decay. It is a well known result, that mollified functions f ∗ m k converge to f in L p ( R d ),1 ≤ p < ∞ when k tends to infinity, see for example Satz 2.15 in Alt (2011).Figure 5 displays the integrand in the stiffness matrix of the Black-Scholes model.The integrand is evaluated on three subintervals of the semi-infinite integration18 × -3 ε -Molli fi ed Integrand of A bsij , i − j = 0 × -3 ε = 0 . hε = 0 . h ξ
100 110 120 130 140 150 160 170 180 190 200 × -4 ε -Molli fi ed Integrand of A bsij , i − j = 0 × -7 ε = 0 . hε = 0 . h ξ × -10 ε -Molli fi ed Integrand of A bsij , i − j = 0 × -177 ε = 0 . hε = 0 . h Figure 5: The integrand of A kl , the stiffness matrix of the Black-Scholes model withmollified hat functions as basis functions for the main diagonal entry, l − k = 0.region. The grid setting is identical to the one of Figure 3. Instead of classic hatfunctions their mollified counterparts have been employed as basis functions usingthe Gaussian mollifier of Example 4.10 as smoothing influence. Even with just aslight mollification influence, ε = 0 . h , the the integrand decays noticeably fasterthanks to the mollification. For moderate values of ε = 0 . h the integrand decaysto zero rapidly.Before we test the performance of this approach to modify the Galerkin schemein Section 5 below, we introduce our second approach. Instead of mollification of piecewise linear basis functions, we can choose basis func-tions that display higher regularity itself. We therefore investigate a well-establishedclass of Finite Element basis functions as candidates for our purposes, namely cubicsplines. Spline theory is a well-investigated field that applies to a much broadercontext than we consider here. We refer the reader to Schumaker (2007) for an in-troduction and overview. From our perspective, splines are smooth basis functions.19heir Fourier transform is accessible and the theory of function spaces they span iswell-established. We give the definition of the Irwin-Hall cubic spline that inheritsis name from related probability distribution.We define the univariate
Irwin-Hall spline ϕ : R → R + by ϕ ( x ) = 14 ( x + 2) , − ≤ x < − | x | − x + 4 , − ≤ x < − x ) , ≤ x ≤ , elsewhere (49)for all x ∈ R . The spline ϕ has compact support on [ − ,
2] and is a cubic spline.We use it to define a spline basis:
Definition 4.1 (Spline basis functions on an equidistant grid) . Choose N ∈ N .Assume an equidistant grid Ω = { x , . . . , x N } , x j ∈ R for all j = 1 , . . . , N , withmesh fineness h > . Let ϕ the Irwin-Hall spline of (49) . For j = 1 , . . . , N define ϕ j ( x ) = ϕ (( x − x j ) /h ) , ∀ x ∈ R . We call ϕ j the spline basis function associated to node j . For a given grid Ω = { x , . . . , x N } , x j ∈ R , Definition 4.1 provides the set ofspline basis functions that we also use in our numerical implementation, later. Instandard spline literature, the set of Irwin-Hall spline basis functions is usually en-riched with additional splines associated with the first and the last node of thegrid that provide further flexibility in terms of Dirichlet and Neumann boundaryconditions, see for example Schumaker (2007). We omit the three Irwin-Hall basisfunctions associated with either of the first and the last grid nodes thus implicitlyprescribing Dirichlet, Neumann and second order derivative zero boundary condi-tions. Lemma 4.11 (Fourier Transform of the Irwin-Hall spline) . Let ϕ be the Irwin-Hallcubic spline of (49) . Then its Fourier transform b ϕ is given by b ϕ ( ξ ) = 3 ξ (cos(2 ξ ) − ξ ) + 3) (50) for all ξ ∈ R . The proof of the Lemma follows by elementary calculation. This immediatelygives the following
Corollary 4.12 (Fourier Transform of spline basis functions on an equidistant grid) . Choose N ∈ N . Assume an equidistant grid Ω = { x , . . . , x N } , x j ∈ R for all j = 1 , . . . , N , with mesh fineness h > and let ϕ j be the spline basis functionassociated with node j as defined in Definition 4.1. Its Fourier transform is givenby c ϕ j ( ξ ) = e iξx j h ξ (cos(2 ξh ) − ξh ) + 3) for all ξ ∈ R . c ϕ ( ξ ) Graph of c ϕ hat hat ε spline ξ
100 110 120 130 140 150 160 170 180 190 200 c ϕ ( ξ ) × -4 ξ c ϕ ( ξ ) × -6 Figure 6: Graphs of the Fourier transforms of all basis function candidates presentedin this section, evaluated over three subintervals of R + . The mesh is chosen with h = 1, the mollification parameter is again set to ε = 0 . h .Figure 6 compares the decay behaviour of the Fourier transforms of all threebasis function types. As Figure 2 already illustrated, the Fourier transform of theclassic hat functions exhibits both slow decay rates and oscillatory behaviour. Instark contrast the Fourier transforms of the mollified hats as well as the Fouriertransform of Irwin-Hall splines visually decays to zero instantly. In case of themollified hat functions this is due to the exponential decay of the Fourier transformof the Gaussian mollifier while for splines Corollary 4.12 displays a polynomial decayof order 4. In this regard, both alternatives to the classic hat functions are promisingcandidates for the implementation of the symbol method of Corollary 4.4. In theupcoming Section 5 we put that promise to the test. In this section we implement the pricing PIDEs for plain vanilla call and put optionsand test the two approaches to the symbol method experimentally.
Theorem 5.1 (Feynman-Kac) . Let ( L t ) t ≥ be a (time-homogeneous) L´evy process.Consider the PIDE ∂ t V C,P + A V C,P = 0 , for almost all t ∈ (0 , T ) V C,P (0) = g C,P , (51)21 here A is the operator associated with the symbol of ( L t ) t ≥ . Assume further theassumptions (A1)–(A3) of Eberlein and Glau (2011) to hold. Then (51) possessesa unique weak solution V C,P ∈ W (0 , T ; H α/ η ( R d ) , L η ( R d )) (52) where α > is the Sobolev index of the symbol of ( L t ) t ≥ and η ∈ R d is chosen ac-cording to Theorem 6.1 in Eberlein and Glau (2011). If additionally g C,Pη ∈ L ( R d ) ,then the relation V C,P ( T − t, x ) = E (cid:2) g C,P ( L T − t + x ) (cid:3) (53) holds for all t ∈ [0 , T ] , x ∈ R d .Proof. The result proved in Eberlein and Glau (2011) and follows from Theorem 6.1therein.
Remark 5.2.
Setting g C,P = g C in (51) , the payoff profile of a European calloption, results in V C being the price of a European call option. Analogously, setting g C,P = g P , the payoff profile of a European put option, results in V P being the priceof a European call option. Algorithm 1 summarizes the abstract structure of a general FEM solver based onthe symbol method. By plugging the symbol associated to the model of choice intothe computation of line 9 of the algorithm, the solver instantly adapts to that model.In other words, only one line needs to be specified to obtain a model specific solverfor option pricing. As Examples 4.6, 4.7, 4.8, 4.9 and others emphasize, the symbolexists in analytically (semi–)closed form for many models, indeed. Algorithm 1 thusprovides a very appealing tool for FEM pricing in practice.
As we derive prices of plain vanilla European call and put options, the solution tothe respective pricing PIDE is defined on the whole real line. As a first step towardsa discretization, we want to truncate the domain to bounded interval ( a, b ) and wechoose to implement zero boundary conditions. To do so, we follow the standardprocedure to subtract an appropriate auxiliary function ψ . Having chosen ψ , theresulting modified problem for φ = V C,P − ψ is ∂ t φ ( t, x ) + ( A φ ) ( t, x ) = f ( t, x ) , ∀ ( t, x ) ∈ (0 , T ) × R φ (0 , x ) = g Ψ ( x ) , ∀ x ∈ R , (54)where g Ψ ( x ) = g ( x ) − ψ (0 , x ) for all x ∈ R and the right hand side f is given by f ( t, x ) := − ( ∂ t ψ ( t, x ) + ( A ψ )( t, x )) . The solution V C,P to the original problem (51) can easily be restored by V C,P = φ + ψ . We establish the properties that ψ needs to provide, later, where we willpresent some examples, as well.The right hand side in vector notation is given by F ( t k ) = ( F ( t k ) , . . . , F N ( t k )) ∈ R N for each t k on the time grid with F j ( · ), j = 1 , . . . , N , given by F j = − Z R ( ∂ t ψ ( t, x ) + ( A ψ )( t, x ) + rψ ( t, x )) ϕ j ( x ) d x (55)22 lgorithm 1 A symbol method based FEM solver Choose equidistant space grid x i , i = 1 , . . . , N Choose basis functions ϕ i , i = 1 , . . . , N , with ϕ i ( x ) = ϕ ( x − x i ) for some ϕ Choose equidistant time grid T j , j = 0 , . . . , M procedure Compute Mass Matrix M Derive the mass matrix M ∈ R N × N by M kl = R R ϕ l ( x ) ϕ k ( x ) d x, ∀ k, l = 1 , . . . , N procedure Compute Stiffness Matrix A Derive the stiffness matrix A ∈ R N × N by plugging the symbol A of thechosen model into the following formula and computing A kl = π R R A ( ξ ) e iξ ( x k − x l ) | c ϕ ( ξ ) | d ξ, ∀ k, l = 1 , . . . , N using numerical integration procedure Run Theta Scheme
Choose a function ψ to subtract from the original pricing problem to obtain azero boundary problem and retrieve the respective basis function coefficientvectors ψ k ∈ R N , k = 1 , . . . , M . Consider the suggestions by Lemma 5.3 orLemma 5.4 for plain vanilla European options below. Choose an appropriate basis function coefficient vector V ∈ R N matchingthe initial condition of the transformed problem Derive the right hand side vectors F k ∈ R N , k = 0 , . . . , M . ConsultLemma 5.3 or Lemma 5.4 matching the choice of ψ . Choose θ ∈ [0 ,
1] and run the iterative scheme for k = 0 : ( M − V k +1 = ( M + ∆ t θ A ) − (cid:0) ( M − ∆ t (1 − θ ) A ) V k + F k + θ (cid:1) end procedure Reconstruct Solution to Original Problem
Add previously subtracted right hand side ψ to the solution of the trans-formed problem by computing e V k = V k + ψ k , k = 0 , . . . , M to retrieve the basis function coefficient vectors e V k , k = 0 , . . . , M , to theoriginal pricing problemfor all j = 1 , . . . , N .We observe that the operator A applied to the auxiliary function ψ appears inthe integral. For the same reasons as in the computation of the stiffness matrixentries, we decide to apply the symbol method for the computation of the entries ofthe right hand side F ∈ R N . We pursue these considerations in following section. F First, we need to choose an appropriate auxiliary function ψ . As its purpose is toenable us to truncate the domain and insert zero boundary conditions, we need toinspect the limit behaviour of the price value V C ( x, t ) → , x → −∞ , t ∈ [0 , T ] V C ( x, t ) → e x − Ke − rt , x → + ∞ , t ∈ [0 , T ] (56)23or call options and V P ( x, t ) → Ke − rt − e x , x → −∞ , t ∈ [0 , T ] V P ( x, t ) → , x → + ∞ , t ∈ [0 , T ] (57)for put options. This is the usual way to obtain the auxiliary function. Now, inregard to our specific approach, relying on the Fourier transforms, we identify ad-ditional desirable feature for the auxiliary function. We denote b ψ ( t, z ) := \ ψ ( t, · )( z ).Consider a smooth function ψ : [0 , T ] × R → R such that ψ ( t ) ∈ H α/ η ( R ) for all t ∈ [0 , T ] for some η ∈ R . Then, for the second summand in (55) we have by applyingthe symbol method of Lemma 4.1 that Z R ( A ψ )( t, x ) ϕ j ( x ) d x = 12 π Z R A ( ξ − iη ) b ψ ( t, ξ − iη ) c ϕ j ( ξ + η ) d ξ, (58)where A denotes the symbol of the model. With the above identity, we are able toderive the right hand side ( F j ) j =1 ,...,N of (55) in terms of Fourier transforms by F j = − π Z R (cid:16)d ∂ t ψ ( t, ξ − iη ) + A ( ξ − iη ) b ψ ( t, ξ − iη ) + r b ψ ( t, ξ − iη ) (cid:17) c ϕ j ( ξ + η ) d ξ. This shows that ψ is numerically suitable for the purpose of localizing the pricingPIDE if ψ is quickly evaluable on the region [ a, b ] × [0 , T ] and the integrals deter-mining F j can be numerically evaluated fast for all j = 1 , . . . , N . The first featureallows retransforming the solution of the localized problem into the solution of theoriginal pricing PIDE, while the second grants the fast numerical derivation of theright hand side in equation (54). These considerations lead us to the following list ofdesirable features for the auxiliary function ψ that is required to obey the respectivelimit conditions (56), (57):1. a (semi-)closed expression of the function ψ ,2. a (semi-)closed expression of its Fourier transform b ψ
3. and fast decay of | b ψ ( ξ ) | and | d ∂ t ψ ( ξ ) | for | ξ | → ∞ .The smoother ψ , the faster | b ψ | decays. We thus need different functions ψ to subtractthat not only fulfil the appropriate boundary conditions (56) or (57) but that arealso as smooth as possible.In the following two subsections we will analyze two candidates for ψ that displaythe desired features.A first suggestion for ψ consists in using Black-Scholes prices as functions in x =log( S ) ∈ [ a, b ] and time to maturity t ∈ [0 , T ] for localization of the pricing PIDE.We express the price of a European option with payoff profile g C,P in the Black-Scholes model in terms of (generalized) Fourier transforms and define ψ accordingly,as the following Lemma explains. Lemma 5.3 (Subtracting Black-Scholes prices) . Choose a Black-Scholes volatility σ > . Define ψ to be the associated Black-Scholes price, ψ ( t, x ) = ψ bs ( t, x ) = e − ηx e − rt π Z R e iξx d g C,P ( − ( ξ + iη )) ϕ bs t,σ ( ξ + iη ) d ξ, (59)24 ith ϕ bs t,σ ( z ) = e tA bs ( z ) . We denote A the symbol of the associated operator A , r ≥ the prevailing risk-free interest rate and choose η < − and η > in the case of acall option and η for the put. Then the right hand side F : [0 , T ] → R N computes to F j ( t ) = 12 π Z R (cid:18) (cid:16) A bs − A (cid:17) ( ξ − iη ) (cid:19)d g C,P ( ξ − iη ) exp (cid:16) − t (cid:16) r + A bs ( ξ − iη ) (cid:17)(cid:17) c ϕ j ( ξ + iη ) d ξ (60) for all j = 1 , . . . , N .Proof. In order to derive the right hand side, we need to represent ψ in Fourierterms. Since for call and put options, ψ / ∈ L ( R ), we compute the (generalized)Fourier transform of ψ or the Fourier transform of ψ η = e η · g , respectively. We get ψ η ( t, x ) = e − rt π Z R e − iξx d g C,P ( ξ − iη ) ϕ bs t,σ ( − ( ξ − iη )) d ξ. (61)The integral in (61) is a Fourier (inversion) integral. We read off c ψ η ( t, ξ ) = d g C,P ( ξ − iη ) exp (cid:16) − t (cid:16) r + A bs ( ξ − iη ) (cid:17)(cid:17) , (62)where we used the relation between the characteristic function and the symbol of aprocess. Now, [ ∂∂t ψ η ( t, ξ ) = − (cid:16) r + A bs ( ξ − iη ) (cid:17) c ψ η ( t, ξ ) . (63)Finally, since ψ bs ∈ H α/ η ( R ), we have that Z R ( A ψ bs )( t, x ) ϕ j ( x ) d x = 12 π Z R A ( ξ − iη ) \ ψ bs ( t, · )( ξ − iη ) [ ϕ j − η ( ξ ) d ξ. (64)Collecting our results proves the claim.The candidate ψ = ψ bs matches the desired criteria. It is quickly evaluable,since functions for yielding Black-Scholes prices are implemented in many code li-braries. Also, the integral in (60) is numerically accessible, since the integranddecays fast. Observe that FFT techniques could be employed to computed F j ( t ) forall j = 1 , . . . , N simultaneously. A major disadvantage of choosing ψ = ψ bs , how-ever, lies in the fact that t ∈ [0 , T ] can not be separated from the integrand in (60).Consequently, F j ( t k ), j = 1 , . . . , N , k = 1 , . . . , M , must be numerically evaluatedon each time grid node individually. This results in significant numerical cost. Wetherefore present a second candidate for ψ that avoids this issue. Lemma 5.4 (Subtracting Quasi-Hockey stick multiplied by Normal) . Let σ ψ > .Define ψ C in the call option and ψ P in the put option case by ψ C ( t, x ) = (cid:0) e x − Ke − rt (cid:1) Φ( x ) , ( t, x ) ∈ [0 , T ] × [ a, b ] ,ψ P ( t, x ) = (cid:0) Ke − rt − e x (cid:1) (1 − Φ( x )) , ( t, x ) ∈ [0 , T ] × [ a, b ] , (65)25 here Φ denotes the cumulative distribution function of the normal N (0 , σ ψ ) distri-bution. Furthermore, in the call option case choose η < − and η > in the putoption case. Then, the right hand side F : [0 , T ] → R N computes to F j ( t ) = 12 π Z R (cid:0) A ( ξ − iη ) + r (cid:1) c f N ( ξ − i ( η + 1)) iξ + ( η + 1) c ϕ j ( ξ + iη ) d ξ − e − rt K Z R (cid:0) A ( ξ − iη ) (cid:1) c f N ( ξ − iη ) iξ + η c ϕ j ( ξ + iη ) d ξ ! , (66) for all j = 1 , . . . , N with t ∈ [0 , T ] , where A is the symbol of the associated operator A and with c f N ( ξ ) = exp (cid:18) − σ ψ ξ (cid:19) , the Fourier transform of the normal N (0 , σ ψ ) density.Proof. We consider the call option case first. To derive the expression for F j in (66)we need to compute the Fourier transform of (the appropriately weighted) ψ C . Wechoose η < − t ∈ [0 , T ] arbitrary but fix and compute for K = 1, \ ψ Cη ( t, · )( ξ ) = Z R e iξx e ηx (cid:0) e x − e − rt (cid:1) Φ( x ) d x = Z R e iξx e ( η +1) x Φ( x ) d x − e − rt Z R e iξx e ηx Φ( x ) d x. (67)Integration by parts and l’Hˆopital’s rule that Z R e iξx e ( η +1) x Φ( x ) d x = − iξ + ( η + 1) Z R e i ( ξ − i ( η +1)) x f N ( x ) d x, (68)which can be expressed in terms of the Fourier transform of the normal distributionyielding Z R e iξx e ( η +1) x Φ( x ) d x = − c f N ( ξ − i ( η + 1)) iξ + ( η + 1) . (69)Equivalently, we obtain for the second integral in (67) that Z R e iξx e ηx Φ( x ) d x = − c f N ( ξ − iη ) iξ + η . (70)Assembling these results we find \ ψ Cη ( t, · )( ξ ) = − c f N ( ξ − i ( η + 1)) iξ + ( η + 1) + e − rt c f N ( ξ − iη ) iξ + η . (71)We deduce F j ( t ) = 12 π Z R (cid:0) A ( ξ − iη ) + r (cid:1) c f N ( ξ − i ( η + 1)) iξ + ( η + 1) c ϕ j ( ξ + iη ) d ξ − e − rt Z R A ( ξ − iη ) c f N ( ξ − iη ) iξ + η c ϕ j ( ξ + iη ) d ξ ! (72)26ith c f N ( ξ ) = exp (cid:18) − σ ψ ξ (cid:19) . For the put option case we choose as defined in (65), ψ P ( x, t ) = (cid:0) Ke − rt − e x (cid:1) (1 − Φ( x ))= (cid:0) e x − Ke − rt (cid:1) (Φ( x ) − . (73)Since ∂∂x (Φ( x ) −
1) = ∂∂x Φ( x ) , ∀ x ∈ R , (74)the computations for c ψ Pη follow along the same lines as they do for the call and weget the relation \ ψ Pη ( t, · )( ξ ) = \ ψ Cη ( t, · )( ξ ) , ∀ ( t, ξ ) ∈ [0 , T ] × R , (75)for η set to some η >
0, which proves the claim.
Remark 5.5 (Computational features of ψ C and ψ P ) . While ψ C serves as localizingfunction for the call option case, ψ P can be used in the put option case. Both candi-dates are based on the payoff functions of call and put options but avoid the lack ofdifferentiability with respect to x in x = log( Ke − rt ) for t ∈ [0 , T ] . As a consequence,both ψ C and ψ P are very smooth functions and thus fulfil the requirements collectedabove when σ ψ is chosen small enough. Additionally, the two integrals in (66) donot depend on the time variable t ∈ [0 , T ] and thus need to be computed only once foreach basis function ϕ j . This results in a significant acceleration in computationaltime compared to the suggestion ψ = ψ bs of Lemma 5.3. The previous sections have outlined the necessary consecutive phases in setting up aFinite Element solver for option pricing. In order to obtain a flexible FEM solver forpricing PIDEs in L´evy models, the previous section introduced the symbol methodwhich considers all components of the FEM solver in Fourier space. There, compo-nents are based on the symbol instead of the L´evy measure and become numericallyaccessible. Many examples of asset models for which the associated symbols existin analytically closed form have deemed this alternative approach being worthwhileto pursue. At the same time, however, smoothness of the FEM basis functions be-came a critical issue. In a second step, we therefore investigated two approaches tocombine smoothness and numerical accessibility. Mollified hat functions and splineswere introduced as promising examples to construct a symbol method based FEMsolver with.This section will put that promise to the test. We therefore implemented thesymbol method for both mollified hats and splines. In stark contrast to a directimplementation of (16) and (55), the symbol method enjoys the flexibility of beingable to easily plug in the symbol of any L´evy model for which it is available in ana-lytically closed form. The model restriction of that first approach thus disappears.Having first implemented the method for the Merton model, the extension of the27ode to the NIG and the CGMY model comes with virtually no additional effort. Inthis regard, the method impressively underlines its appeal for applications in prac-tice where the suitability of a model might depend on the asset class it is employedfor. An institution that needs to maintain pricing routines for several asset classeswill thus cherish the flexibility that the symbol method offers, recall Algorithm 1 inthis regard which sketches the implementation of a general, symbol method basedFEM solver that easily adapts to various models.Finally, we conduct an empirical order of convergence study. We consider theunivariate Merton, CGMY and NIG model and investigate the empirical rates ofconvergence for the different implementations as Table 1 summarizes.
Model Symbol Parameter choices Implemented basis functions
Mollified hats SplinesMerton Example 4.7 σ = 0 . α = − . X X β = 0 . λ = 3CGMY Example 4.8 C = 0 . G = 23 . X X M = 27 . Y = 1 . α = 12 . β = − . X X δ = 0 . r = 0 . K = 1 as an example, thus considering the payoff function g ( x ) = max( e x − , . (76)In each study we compute FEM prices for N k basis functions, with N k = 1 + 2 k , k = 4 , . . . , N = 17 basis functions in the most coarse and N = 513 basis functionsin the most granular case. On each grid, the nodes that basis functions are associatedwith are equidistantly spaced from another and the basis functions always span thespace interval Ω = [ − , N time grid =2000 equidistantly spaced time nodes spanning a grid range of two years up untilmaturity, thus covering a time to maturity interval of[ T , T N time ] , with T = 0 and T N time = 2 . (78)28or each k = 4 , . . . ,
9, the resulting price surface constructed by N k basis functionsin space and N time = 2000 grid points in time is computed. A comparison of thesesurfaces is drawn to a price surface of most granular structure based on the sametype of basis function. We call this most granular surface true price surface. It restson N true = N = 1 + 2 = 2049 basis functions in space and N time grid pointsin time spanning the same grid intervals as above, that is Ω = [ − ,
5] in space and[0 ,
2] in time, respectively. The underlying FEM implementation is thus based ondistances h true between grid nodes that basis function are associated with of h mollified hattrue = (5 − ( − / (2 + 2 ) ≈ . ,h splinestrue = (5 − ( − / (4 + 2 ) ≈ . , ∆ t true = 2 / (2000 − ≈ .
001 (79)in space and time, respectively. Note that all space grids are designed in such away that the log-strike log( K ) = 0 is one of the space nodes. For each model andmethod and each k = 4 , . . . ,
9, the (discrete) L error ε L is calculated as ε L ( k ) = vuut ∆ t true · h true · N time X i =1 N true X j =1 (cid:16) P rice true ( i, j ) − P rice k ( i, j ) (cid:17) , wherein P rice true ( i, j ) is the value of the true pricing surface at space node j ∈{ , . . . , } and time node i = 1 , . . . , P rice k ( i, j ) is the respective,linearly interpolated value of the coarser pricing surface with only N k basis functionnodes.Figure 7 summarizes the results of the six studies of empirical order of conver-gence in the Merton, the NIG and the CGMY model in a symbol based implemen-tation once using mollified hats and once using splines as basis functions. In eachimplementation and for all considered models, the (discrete) L error decays expo-nentially with rate 2. The convergence result of Theorem 5.4 by von Petersdorff andSchwab (2003) suggest that this is the best possible rate we can hope for, whichyields the experimental validation of both approaches.29 l og ( ε L ( k )) e rr o r -10-5051015 Merton model (spline)Line with slope = -2 k l og ( ε L ( k )) e rr o r -10-5051015 NIG model (spline)Line with slope = -2 k l og ( ε L ( k )) e rr o r -10-5051015 CGMY model (spline)Line with slope = -2 k l og ( ε L ( k )) e rr o r -10-5051015 Merton model (mhat)Line with slope = -2 k l og ( ε L ( k )) e rr o r -10-5051015 NIG model (mhat)Line with slope = -2 k l og ( ε L ( k )) e rr o r -10-5051015 CGMY model (mhat)Line with slope = -2
Figure 7: Results of the empirical order of convergence study for the Merton, theNIG and the CGMY model using mollified hats (left pictures) and splines (rightpictures) as basis functions. All models are parametrized as stated in Table 1.Additionally, part of a straight line with (absolute) slope of 2 is depicted in eachfigure serving as a comparison. 30
Conclusion and outlook
We have presented a finite element solver that is highly flexible in the model choiceand that maintains numerical feasibility. Invoking the symbol was key. The transi-tion into Fourier space has introduced smoothness as a new requirement to the basisfunctions. We have presented splines and mollified hats as compatible basis func-tions in our approach. Several numerical examples have confirmed the convergencesrates expected by standard theory in both cases. The contribution of the mollifi-cation to the error has not yet been addressed in the literature. One possibility toderive such error estimates lies in adapting the perturbation analysis of von Peters-dorff and Schwab (2003) treating mollification as a numerical perturbation of thestiffness matrix. Let us mention several possible extensions of the approach. First,the implementation naturally extends to time-inhomogeneous L´evy models that weneglected here for notational convenience. Second, combining the symbol methodwith Wavelet basis functions allows for compression techniques that might furtherimprove the overall numerical performance, as Hilber et al. (2013) outlines. Third,the polynomial decay that we observe in our numerical experiments can possibly beimproved to exponential rates by invoking an hp -discontinuous Galerkin scheme, seee.g. Sch¨otzau and Schwab (2006). A Proofs
A.1 Proof of Lemma 3.1
Proof.
We first consider ϕ, ψ ∈ C ∞ ( R ).For F ≡ a jump ( ϕ, ψ ) := − Z R Z R (cid:16) ϕ ( x + y ) − ϕ ( x ) − ϕ ′ ( x ) h ( y ) (cid:17) F (d y ) ψ ( x ) e h η,x i d x, needs to be carefully derived. In order to exploit the identity ϕ ( x + y ) − ϕ ( x ) − yϕ ′ ( x ) = Z y Z z ϕ ′′ ( v ) d v d z we split the integral with respect to the L´evy measure in three parts, set c ( F ) := R | y | < (cid:0) y − h ( y ) (cid:1) F (d y ) − R | y | > h ( y ) F (d y ) and obtain a jump ( ϕ, ψ ) := − Z R Z | y | < Z y Z z ϕ ′′ ( x + v ) d v d zF (d y ) ψ ( x ) e h η,x i d x − c ( F ) Z R ϕ ′ ( x ) ψ ( x ) e h η,x i d x − Z R Z | y | > (cid:0) ϕ ( x + y ) − ϕ ( x ) (cid:1) F (d y ) ψ ( x ) e h η,x i d x. R y R z (cid:12)(cid:12) ϕ ′′ ( v ) (cid:12)(cid:12) d v d z ≤ cy with some constant c > y ∈ [ − , Z R Z | y | < Z y Z z (cid:12)(cid:12) ϕ ′ ( x + v ) (cid:12)(cid:12) d v d zF (d y ) (cid:12)(cid:12) ψ ′ ( x ) + 2 ηψ ( x ) (cid:12)(cid:12) e h η,x i d x ≤ (1 + 2 η ) k ϕ k H η k ψ k H η Z | y | < y F (d y ) (80)we can apply the theorem of Fubini and partial integration to obtain − Z R Z | y | < Z y Z z ϕ ′′ ( x + v ) d v d zF (d y ) ψ ( x ) e h η,x i d x = Z R Z | y | < Z y Z z ϕ ′ ( x + v ) d vF (d y ) (cid:0) ψ ′ ( x ) + 2 ηψ ( x ) (cid:1) e h η,x i d x. This yields the assertion for ϕ, ψ ∈ C ∞ ( R ).Next, we verify that the bilinear form as stated in Lemma 3.1 is well defined for ϕ, ψ ∈ H η ( R ) and is continuous with respect to the norm of H η ( R ). For F ≡ Z R Z | y | > (cid:12)(cid:12) ϕ ( x + y ) − ϕ ( x ) (cid:12)(cid:12) F (d y ) (cid:12)(cid:12) ψ ( x ) (cid:12)(cid:12) e h η,x i d x ≤ F (cid:0) R \ [ − , (cid:1) k ϕ k L η k ψ k L η . Thus a from Lemma 3.1 is a continuous bilinear form on H η ( R ) × H η ( R ) thatcoincides with (9) on the dense subset C ∞ ( R ) × C ∞ ( R ). This proves the assertion. A.2 Proof of Lemma 4.2
Proof.
To prove the assertion, we verify the conditions of Lemma 3 in Glau (2016b),which provides an abstract robustness result for weak solutions. We first observethat the conditions for f n , f, g n , g coincide with those of Lemma 3 in Glau (2016b).Second, we verify conditions (An1)–(An3) of Lemma 3 in Glau (2016b). Thereforewe assign to each u, v ∈ X the coefficients α k ( u ) , α k ( v ) ∈ R for k ≤ N such that u = P Nk =1 α k ( u ) ϕ k and v = P Nk =1 α k ( v ) ϕ k . Thanks to the finite dimensionsionalityof X , there exists a constant e C > u ∈ X , k u k V ≤ N X k =1 (cid:12)(cid:12) α k ( u ) (cid:12)(cid:12) k u k V ≤ C ′ k u k V . (81)Thanks to (27) there exists a sequence 0 < c n → j, k ≤ N , (cid:12)(cid:12) ( a n − a )( ϕ j , ϕ k ) (cid:12)(cid:12) ≤ c n k ϕ j k V k ϕ k k V . (82)Together with assumption (A2) this yields for all j, k ≤ N , (cid:12)(cid:12) a n ( ϕ j , ϕ k ) (cid:12)(cid:12) ≤ C k ϕ k k V k ϕ k k V . (83)32nequalities (83) and (81) together yield for all u, v ∈ X , (cid:12)(cid:12) a n ( u, v ) (cid:12)(cid:12) ≤ N X k =1 N X j =1 (cid:12)(cid:12) α k ( u ) α j ( v ) (cid:12)(cid:12)(cid:12)(cid:12) a n ( ϕ k , ϕ j ) (cid:12)(cid:12) ≤ C N X k =1 N X j =1 (cid:12)(cid:12) α k ( u ) α j ( u ) (cid:12)(cid:12) k ϕ k k V k ϕ k k V ≤ C e C k u k V k v k V , which shows that condition (An1) of Lemma 3 in Glau (2016b) is satisfied.Due to inequalities (82) and (81), we have for all u ∈ X , (cid:12)(cid:12) ( a − a n )( u, u ) (cid:12)(cid:12) ≤ N X k =1 N X j =1 (cid:12)(cid:12) α k ( u ) α j ( u ) (cid:12)(cid:12)(cid:12)(cid:12) a n ( ϕ k , ϕ j ) (cid:12)(cid:12) ≤ c n N X k =1 N X j =1 (cid:12)(cid:12) α k ( u ) α j ( u ) (cid:12)(cid:12) k ϕ j k V k ϕ k k V ≤ c n e C k u k V , which shows assumption (An3) of Lemma 3 in Glau (2016b).Finally, from assumption (A1) and the last inequality for all u ∈ X we obtain a n ( u, u ) ≥ a ( u, u ) − (cid:12)(cid:12) ( a − a n )( u, u ) (cid:12)(cid:12) ≥ G k u k V − G ′ k u k H − c n e C k u k V , which shows that there exists N ∈ N such that condition (An2) of Lemma 3 inGlau (2016b) is satisfied for all n > N . This shows the assertion of Lemma 4.2. References
Alt, H., 2011. Lineare Funktionalanalysis: Eine anwendungsorientierte Einf¨uhrung,6th Edition. Springer Lehrbuch. Springer.Brummelhuis, R., Chan, R., 2014. An RBF scheme for option pricing in exponentialL´evy models. Applied Mathematical Finance 21, 238–269.Carr, P., Geman, H., Madan, D. B., Yor, M., 2002. The fine structure of assetreturns: An empirical investigation. Journal of Business 75 (2), 305–332.Chan, R., Hubbert, S., 2014. Options pricing under the one-dimensional Jump-diffusion model using the radial basis function interpolation scheme. Review ofDerivatives Research 17 (2), 161–189.Cont, R., Voltchkova, E., 2005a. A finite difference scheme for option pricing injump diffusion and exponential L´evy models. SIAM Journal on Numerical Analysis43 (4), 1596–1626. 33ont, R., Voltchkova, E., 2005b. Integro-differential equations for option prices inexponential L´evy models. Finance and Stochastics 9 (3), 299–325.Eberlein, E., Glau, K., 2011. PIDEs for pricing European options in L´evy models –a Fourier approach, technische Universtit¨at M¨unchen.Eberlein, E., Glau, K., 2014. Variational solutions of the pricing PIDEs for Europeanoptions in L´evy models. Applied Mathematical Finance 21 (5-6), 417–450.F¨ollmer, H., 2009. Alles richtig und trotzdem falsch? Mitteilungen der DMV 17,148–154.Glau, K., 2016a. Classification of L´evy processes with parabolic Kolmogorov back-ward equations, forthcoming in SIAM Journal Theory of Probability and its Ap-plications.Glau, K., 2016b. Feynman-Kac type formula for L´evy processes with discontionuouskilling rate, accepted for publication in Finance and Stochastics.Hilber, N., Reich, N., Winter, C., Schwab, C., 2009. Numerical methods for L´evyprocesses. Finance and Stochastics 13 (4), 471–500.Hilber, N., Reichmann, O., Schwab, C., Winter, C., 2013. Computational Methodsfor Quantitative Finance. Springer.Itkin, A., 2015. Efficient solution of backward jump-diffusion PIDEs with splittingand matrix exponentials, forthcoming in Journal of Computational Finance.Itkin, A., Carr, P., 2012. Using Pseudo-Parabolic and Fractional Equations for Op-tion Pricing in Jump Diffusion Models. Computational Economics 40, 63–104.Marazzina, D., Reichmann, O., Schwab, C., 2012. hp-DGFEM for Kolmogorov-Fokker-Planck Equations of Multivariate L´evy Processes. Mathematical Modelsand Methods in Applied Sciences (M3AS) 22 (1).Matache, A.-M., Nitsche, P.-A., Schwab, C., 2005. Wavelet Galerkin pricing of Amer-ican options on L´evy driven assets. Quantitative Finance 5 (4), 403–424.Matache, A.-M., von Petersdorff, T., Schwab, C., 2004. Fast deterministic pricing ofoptions on L´evy driven assets. Mathematical Modelling and Numerical Analysis38 (1), 37–71.Salmi, S., Toivanen, J., von Sydow, L., 2014. An IMEX-scheme for pricing op-tions under stochastic volatility models with jumps. SIAM Journal on ScientificComputing 36 (4), B817–B834.Sch¨otzau, D., Schwab, C., 2006. Time Discretization of Parabolic Problems by the hphp