A FINITE DIFFERENCE METHOD FOR PRICING EUROPEAN AND AMERICAN OPTIONS UNDER A GEOMETRIC PROCESS

. In this paper we develop a numerical approach to a fractional-order diﬀerential Linear Complementarity Problem (LCP) arising in pricing Euro- pean and American options under a geometric L´evy process. The LCP is ﬁrst approximated by a nonlinear penalty fractional Black-Scholes (fBS) equation. We then propose a ﬁnite diﬀerence scheme for the penalty fBS equation. We show that both the continuous and the discretized fBS equations are uniquely solvable and establish the convergence of the numerical solution to the viscosity solution of the penalty fBS equation by proving the consistency, stability and monotonicity of the numerical scheme. We also show that the discretization has the 2nd-order truncation error in both the spatial and time mesh sizes. Numerical results are presented to demonstrate the accuracy and usefulness of the numerical method for pricing both European and American options under the geometric L´evy process.


1.
Introduction. Pricing financial options has attracted much attention from both mathematicians and financial engineers in the last decade. A financial option is a contract that gives its owner the right, not obligation, to buy (call option) or to sell (put option) a fixed quantity of a stock at a fixed price (strike price) on (European type) or before (American type) a given expiry date. In a complete market, Black & Scholes [2] demonstrated that the price of a European option on a stock, whose price follows geometric Brownian motion with constant drift and volatility, satisfies a second order partial differential equation, known as the Black-Scholes (BS) equation (or model), with proper boundary and terminal conditions. One major shortcoming of the BS model is that the Gaussian shocks used in the model underestimate the probability that stock prices exhibit large movements or jumps over small time steps as illustrated by empirical data in financial market. To overcome this problem, we assume that the underlying stock price S t of an option follows, as proposed in [5], the following a geometric Lévy process d(ln S t ) = (r − v)dt + dL t (1.1)

WEN CHEN AND SONG WANG
with the solution where T is a future known time, r is the risk-free interest rate, v a convexity adjustment so that the expected value of S T becomes E[S T ] = e r(T −t) S t , and dL t is the increment of a Lévy process under the equivalent martingale measure (EMM). Boyarchenko and Levendorskii [4] proposed the use of a modified Lévy-stable (LS) (Lévy-α-stable) process to model the dynamics of securities. This modification introduces a damping effect in the tails of the LS distribution, which are known as KoBoL processes. Carr, Geman, Madan and Yor [5] proposed the CGMY process including both positive and negative jumps. In this paper, we are concerned with options based on finite moment log-stable (FMLS) processes. In [6], the authors show that a classical hedging portfolio can be substantially improved by employing 'non-local' or fractional differential operators. Since over a time step ∆t, the stock price S t can diffuse or jump to a value S t+∆t far away from S t , the localized information becomes less relevant. The fractional derivative weighs information of the portfolio over a range of values of the underlying stock [6] rather than localized information. When the Brownian motion component is replaced by a Lévy process, the Black-Scholes equation becomes a partial integro-differential equation (PIDE). In [6], by Fourier transform, the PIDE is written as a fractional partial differential equation. In what follows, we refer it to as a fractional Black-Scholes (fBS) equation. Fractional partial differential equations (fPDEs), as generalizations of classical integer-order partial differential equations, are increasingly used to model problems with memory in many areas such as finance and fluid flows. Fractional spatial derivatives are used to model anomalous diffusion or dispersion in which a particle spreads at a rate inconsistent with the classical Brownian motion. Since closed-form solutions to fPDEs of practical significance can rarely be found, various numerical techniques have been proposed for fPDEs. Fractional derivatives can be represented in different forms such as those of Riemann-Liouville (RL) and Grüwald-Letnikov (GL) [19]. Most existing discretization methods have been developed for fPDEs in GL form (cf., for example, [21,22,7,18,26]). Unlike a European option whose value is determined by an fBS equation, the value of an American option under the Lévy process is governed by a linear differential complementarity problem involving the fBS operator. Various penalty methods have been developed for solving complementarity problems in both infinite and dimensions [25,16,13,28,14,29,17,8]. In this work, we develop a numerical method for the fractional differential linear complementarity problem (LCP), or the variational inequality, arising from pricing American options under the Lévy process. We first approximate the LCP by a nonlinear fBS equation using the linear penalty approach used in [30,1,29]. We then develop a finite difference scheme based on a numerical quadrature rule for the spatial integral term and Crank-Nicolson timestepping scheme for the penalized nonlinear fBS equation which contains the fBS governing European option valuation as a special case. The truncation error of this discretization is shown to be of 2nd-order in both space and time. We will show the solution to the discretized system converges to the exact viscosity solution of the penalized fBS equation by proving that the numerical scheme is consistent, stable and monotone. Numerical results will be presented to demonstrate the accuracy and usefulness of the numerical scheme using some model fPDEs and fBS equations.
The organization of this paper is as follows. In Section 2, we will first give a brief account of the continuous LCP governing the American option valuation and apply the penalty method to the LCP to yield a penalized nonlinear fBS equation. We then discuss briely the unique solvability of the penalized fBS equation. In Section 3, we develop a discretization scheme for the fractional derivative and a full discretization scheme for the penalized fBS equation. The consistency, stability, and monotonicity of the numerical method are proved in Section 4. In Section 5, we first use a model fPDE to demonstrate that our numerical method is 2nd-order accurate in both space and time. We then present numerical results on European and American call and put options to show that the method produces practically useful results.
2. The continuous problem and its unique solvability.
2.1. The continuous problem. A time-dependent random variable X t is a Lévy process, if and only if it has independent and stationary increments with the following log-characteristic function in Lévy-Khintchine representation and Ψ(ξ) is the characteristic exponent of the Lévy process which is a combination of a drift component, a Brownian motion component and a jump component. These three components are determined by the Lévy-Khintchine triplet (m, σ 2 , W ). If the Lévy measure is of the form W (dx) = w(x)dx, w(x) is then called the Lévy density. For an LS process, the Lévy density is given by where D > 0, p, q ∈ [−1, 1] and p + q = 1 satisfying 0 < α ≤ 2. The characteristic exponent of the LS process is The parameters α and σ are respectively the stability index and scaling parameter. The parameter s := p − q is called the skewness parameter satisfying −1 ≤ s ≤ 1, and m is a location parameter. When s = 1 (resp. s = −1) the random variable X is maximally skewed to the left (resp. right). When α = 2 and s = 0, it becomes Gaussian case. A particular feature of the FMLS process is that it only exhibits downwards jumps, while upwards movements have continuous paths. The characteristic exponent of the LS process with s = −1, is where ν := 1 2 σ α sec απ 2 is the convexity adjustment of the random walk. For a given function u(x), one form of the α-th derivative of u is for x > x 0 , where x 0 is a given real number and Γ(·) denotes the Gamma function. When u(x 0 ) = 0 and u (x 0 ) = 0, it reduces to the Caputo's representation of the fractional partial derivative. It is shown in [6] that the value U of a European option written on a stock, whose price S follows (1.1) with L t = Ψ F M LS defined in (2.1), is determined by the following fBS equation: with the boundary and terminal conditions: satisfying the compatibility conditions U 0 (T ) = U * (x min ) and U 1 (T ) = U * (x max ), where x = ln S, x min << 0 and x max > 0 are two constants representing the lower and upper bounds for x, and In (2.3), r is the risk-free rate and α ∈ (1, 2) is the order of fractional derivative, and U 0 (t), U 1 (t) and U * (x) are known functions. For vanilla options, U * (x) = max{e x − K, 0} for a call and U * (x) = max{K − e x , 0} for a put, where K denotes the strike price of the option. Note that the original spatial solution domain is (−∞, ∞). In computation, we truncate this infinite domain by the lower and upper bounds x min and x max , as done in (2.3)-(2.4). Thus, we assume that x min and x max are chosen such that x min << 0 and e xmax >> K. Note that the use of (2.2) requires U (x min , t) = 0 and U x (x min , t) = 0 in order to avoid the singularity at x min . Both of these can be achieved, up to a truncation error, by transforming (2.3) into an fBS equation satisfying the homogeneous Dirichlet boundary condition.
Let F (x, t) be the function defined by Clearly, F (x, t) satisfies the boundary conditions (2.4) and (2.9) (up to a truncation error). Also, F is an exponential function of x and thus it is invariant under the 1st and α-th order differentiation operations with respect to x. Taking LF from both sides of (2.3) and introducing a new variable V ( where f (x, t) = LF . The boundary and terminal conditions (2.4)-(2.5) then become From the definitions of u and F and x = ln S, we have Ke xmin e xmax − e xmin + U S (S(x min ), t)e xmin = 0, (2.9) since U S (S(x), t) is bounded on (−∞, x max ). Thus, V x (x min , t) → 0 exponentially as x min → −∞. Therefore, from (2.2), (2.7) and (2.9), we see that the fractional derivative involved in LV now becomes the following Caputo's form: up to a truncation error when x min << 0. While the price of an European option is governed by (2.3)-(2.5), it is well known that the price of an American option is determined by the following a linear complementarity problem [24]: for (x, t) ∈ I × [0, T ) satisfying the boundary and terminal conditions (2.7)-(2.8).
In [24], the authors proposed and analyzed a power penalty method for (2.12) for the case that the Black-Scholes operator is the 2nd order differential operator, i.e., α = 2 in (2.3). The penalty method proposed in [24] is extended to (2.12) in [9]. In this work, we use the linear penalty method to solve (2.12) as used in [1,29], i.e., we approximate (2.12) by the the following penalized fBS equation: with the boundary and initial conditions (2.7)-(2.8), where λ > 1 is a penalty constant and [z] + = max{0, z} for any function z. The convergence of V λ to V for the case that α ∈ (1, 2) and λ > 1 is in [9]. In the present work, we will concentrate on the construction and the convergence of a discretization scheme for (2.13). Note that (2.13) contains (2.6) as a special case when λ = 0. Therefore, in what follows, our discussion will be focused on (2.13) unless mentioned otherwise.

2.2.
The variational problem and its solvability. We now consider the unique solvability of (2.13). First, we formulate it as a variational problem, and then we show that the variational problem has a unique solution. We start this discussion by introducing some function spaces.
For an open set Ω ⊂ R and 1 ≤ p ≤ ∞, we let L p (Ω) = v : ( Ω |v(x)| p dx) 1/p < ∞} denote the space of all p-power integrable functions on Ω equipped with the usual L p -norm · L p (Ω) , and (·, ·) denote the usual inner product. For any γ ∈ (0, 1], we let H γ (R) := v : v and −∞ D γ x v ∈ L 2 (R) | · | γ and · γ are two functionals defined respectively as for any u ∈ H γ (R). Then it is easy to show that | · | γ and · γ are semi-norm and norm on H γ (R) respectively. It has been shown in [12] that H γ (R) equipped with · γ is a Sobolev space. Similarly to the above definition of fractional Sobolev space, we also define the Sobolev space of functions having a support on the finite interval I = (x min , x max ) given by In what follows, we also use ·, · to denote the duality paring between H γ 0 (I) and its dual space H −γ 0 (I). Using the notation defined above, we pose the following variational problem: with the boundary and initial conditions (2.7)-(2.8): where A(·, ·; t) is a bilinear form defined by Using integration by parts, it is easy to verify that Problem (2.1) is the variational problem of (2.13) with (2.6)-(2.8) (cf. [12]). It has also been shown in [12] that the bilinear form A(·, ·; t) is coercive and continuous, as given in the following lemma: Lemma 2.1. There exist positive constants C 1 and C 2 such that for any v, w ∈ H α/2 0 (I), Using this lemma, we have the following theorem. The proof of this theorem is just a re-statement of that of Theorem 3.1 of [24] which shows that the nonlinear form on the RHS of (2.14), (i.e., the nonlinear operator on the RHS of (2.13)) is strongly monotone and continuous, based on Lemma 2.1. For brevity, we omit this discussion.
3. Discretization of (2.13) . We now consider the discretization of the fractional partial differential equation (2.13). Several efficient and accurate discretization scheme have been proposed for linear, nonlinear and penalized 2nd-order Black-Scholes equations [27,23,1,17,15]. However, to our best knowledge, there is essentially no work on the numerical approximation of the penalized fBS equation. In this section, we propose a discretization method for (2.13).

3.1.
Discretization of the α-th derivative. Let the interval I = (x min , x max ) be divided into M sub-intervals with mesh nodes where h = (x max − x min )/M . For clarity, we omit the variable t in this subsection. When 1 < α < 2, from (2.10) we have that, for any i ∈ {1, 2, ..., M }, To discretize I i k , we first rewrite it as In the above, we used a truncated Taylor expansion for V xx (y) − V xx (x k ). The derivatives V xx (x k ) and V xxx (x k ) are then approximated respectively by the following finite differences: where V i denotes an approximation to V (x i ) for any feasible i. The two integrals on the RHS of (3.3) can be evaluated exactly, and using (3.1), it is easy to show that these integrals are given by Therefore, I i k can then be approximated by For P i k and Q i k , we have the following lemma: The proof is trivial and thus it is omitted.
We now show f 2 (k) is also strictly increasing. Differentiating f 2 (k), we have It remains to prove f 3 is strictly increasing. Differentiating f 3 (k), we have (2) The proof of this is trivial and we omit it.
(3) For the finite difference scheme in (3.9), the approximation of the α-th derivative of V (x) = 1 becomes exact, i.e., where ∆t = T /N . Using the central differencing for the first derivative in space, Crank-Nicolson time stepping method and the scheme (3.9) for the α-th derivative, we construct the following discretization scheme for (2.3): with the boundary and terminal conditions: for i = 1, 2, ..., M − 1 and j = 0, 2, ..., N − 1, where V j i denotes an approximation to V (x i , t j ), f k i = f (x i , t k ) for k = j and j + 1, and Let µ = −b ∆t Γ(2−α)h α and η = a ∆t 2h . We rewrite equation (3.14) as for j = 0, 1, ..., N − 1. This system can further be written as the following matrix equation: where I denotes the (M − 1) × (M − 1) identity matrix, Clearly, G and B arise respectively from the discretization of the α-th derivative and the first derivative ∂V /∂x. It is trivial to show that the elements c ij are given by We comment that C is a Toeplitz matrix as the elements in each diagonal of C is a constant. Also, (3.16) is a nonlinear system for V j+1 since the diagonal matrix D is a nonlinear and non-smooth function of V j+1 . Note that D is monotonically increasing in V j+1 , in practice a Newton's or non-smooth Newton's method can be used for solving (3.16) numerically. 4. Consistency, stability and monotonicity of the scheme. In this section, we show that the solution to (3.14) converges to the viscosity solution to (2.13) by proving that the numerical scheme proposed in the previous section is consistent, stable and monotone. We start this discussion with the following theorem: Theorem 4.1. The finite difference scheme for (2.3), defined by (3.14), is consistent, with a truncation error of order O(∆t 2 + h 2 ).
Proof. In what follows, we let C denote a generic positive constant, independent of h and ∆t. We first consider the truncation error in the approximation D α h to D α at x i for any i = 1, 2, ..., M − 1. From (3.2) and (3.3) we have that, for any function V (x) sufficiently smooth on I, where E i denotes the following remainder in the approximation of V (y)−V (x k−1 ) by a truncated Taylor's expansion: for a ξ k ∈ (x k−1 , x k ). From this equality we have, for i = 1, 2, ..., M − 1, where ||V (4) || ∞ denotes the maximum norm of V (4) on I. Since where the definition of S k is obvious. Using the expansion we can easily show that

Thus,
i k=1 S k can be written as 3) The first term on the RHS of (4.3) has the following upper bound: For the Gamma function, when 0 < a < 1 and n is a positive integer, we have Γ(n + a) Γ(n + 1) < n a−1 from which, it can be shown that Using this inequality and the properties of the Gamma function, we have The third term on the RHS of (4.3) can be estimated as follows: Replacing the three terms on the RHS of (4.3) with their respective upper bounds (4.4)-(4.6) and combining the resulting estimate with (4.2) and (4.1), we have |V (4.8) From (4.7), (4.8) and (3.9), we see that the truncation error in the discretization of the α-th derivative is bounded by as the last integral in the above expression exists. Finally, it is well-known that Crank-Nicolson time stepping scheme, the central differencing and the mid-point quadrature rule used in (3.14) are all at least 2nd-order accurate on uniform meshes. Combining this fact with (4.9), we have Therefore, the discretization is consistent.
Theorem 4.2. The finite difference scheme defined by (3.14) is unconditionally stable.
Proof. Let us first consider the case that λ = 0 in (2.13), or d i = 0 in (3.14) for all i.
We use the discrete Fourier transform to prove the stability of the Crank-Nicolson method. Using µ and η introduced in Subsection 3.2, we rewrite (3.14) as for i = 1, 2, ..., M − 1 and any feasible j, wheref j i = (f j i + f j+1 i )/2. This system has the matrix form (3.16) and from the definition (3.17) we see that all the coefficient matrices in (3.16) are Toeplitz matrices. Thus, each of the terms in (3.16) can be written as convolution of a coefficient vector with a finite support, one of (· · · , 0, V n 1 , ..., V n M −1 , 0, · · · ) and (· · · , 0,f j 1 , ...,f j M −1 , 0, · · · ) for n = j and j +1. Applying the discrete Fourier transform to the system, or equivalently replacing V n m andf n m in (4.10) with W n e mξhi andF n e mξhi respectively for all admissible m and n, we have g k e ξ(−k+1)hi + r∆t where ξ ∈ [− π h , π h ], and W n andF n are respectively the discrete Fourier transform of V n andf n for n = j and j + 1. Dividing both sides by e iξhi and rearranging the resulting equation, we get To prove the stability, it suffices to show that for a constant L > 0, independent of h and ∆t. This is obvious if one can prove that A ≥ 0.

For any function
denote the discrete L 2 -norm of u. Then, using the properties of the discrete Fourier and its inverse transforms (particularly Parseval's equality) we have from which we have (recallL is a generic positive constant) Therefore, the numerical method is unconditionally is stable when λ = 0. In the case that λ > 0, from the definition of d i we have that, for any feasible i and k, where ρ k i = 0 or 1. Therefore, (3.14) can be written as 13) whereρ j i = (ρ i i + ρ j+1 i )/2 = 0, 1/2 or 1. Comparing (4.13) with (4.10) we see that (4.13) is in the same form as (4.10) with r∆t replaced with ∆(r + λρ k i ) for k = j or j + 1 and ∆f j i with ∆t f j i +ρ j i λV * i . All of these terms are of the order O(∆t). Thus, following the same analysis presented above for λ = 0, we have that the scheme is also stable when λ > 0. Therefore, we have proved the theorem.
We now show that the numerical scheme is monotone. Proof. For notation simplicity, we first consider the case that λ = 0. Let where f (β) := (1 + 2.5 × 2 β − 3 × 3 β + 4 β ). It now suffices to show that f (β) > 0 for β ∈ (1, 2). Since f (2) = 0, we need only to show that f (β) is strictly decreasing for β ∈ (1, 2). Differentiating f gives It is easy to show, even graphically, that f (β) < 0 for all β ∈ [1, 2]. Therefore, We now use the above result to prove the monotonicity of F j+1 i . When ∆t ≤ 2 r , we have from the definition of F j+1 i that, for any ε > 0 and feasible i and j, When λ > 0, it is easy to see d(V k i ) is monotonically increasing in V k i for any feasible k. Therefore, both d(V j i ) and d(V j+1 i ) are monotone and thus the scheme is also monotone.
Combining Theorems 4.1, 4.2 and 4.3 ,we have the following convergence result. In fact, conventionally, Theorems 4.1 and 4.2 already imply the convergence of our numerical scheme. Barles and Souganidis showed in [3] that any finite difference scheme for a general nonlinear 2nd-order PDE which is locally consistent, stable and monotone generates a solution converging uniformly on a compact subset of [0, T ] × R to the unique viscosity solution of the PDE. In [11] and [10], Cont and Tankov extended this result to partial integro-differential equations (PIDEs). Since (2.3) is an PIDE, Theorem 4.4 is just a consequence of the results established in [3,11,10].
We also comment that though the theoretical results in this section have been established for Crank-Nicolson's time-stepping scheme, they hold true for a general two-level time-stepping scheme with a splitting parameter θ ∈ [0.5, 1]. However, when θ ∈ (0.5, 1], the truncation error in Theorem 4.1 is of the order O(∆t + h 2 ) instead of O(∆t 2 + h 2 ). For brevity, we omit this discussion.

Numerical experiments.
In this section, we first use an example with a known exact solution to demonstrate the rate of convergence of our scheme. We then show the usefulness and practicality of the method by applying it to several European and American option pricing problems. All the computations have been performed in double precision under MATLAB environment.

Example 1. Fractional diffusion equation with non-homogeneous boundary conditions:
∂u(x, t) ∂t The exact solution to the above problem is u(x, t) = te 2x . Note for this test we have u(−5, t) ≈ 0 and u x (−5, t) = 0, and so, we may straightforwardly apply our numerical scheme to this test without the transformation used in Section 2. This problem is solved using a sequence of meshes h k = ∆t k = 1 5 × 2 −k for k = 0, 1, ..., 5. For each k, the following discrete maximum norm is computed: where U = (U j i ) denotes the numerical solution. These computed errors, along with computed rates of convergence log 2 (E k+1 /E k ), for k = 0, 1, ..., 5 are listed in Table  5.1, from which we see that the rates of convergence of our method are of order O(∆t 2 + h 2 ), coinciding with the truncation error established in Theorem 4.1. For comparison, we have also solved the problem using a combination of the Crank-Nicolson time-stepping scheme and two popular existing spatial finite difference methods proposed respectively in [20] and [18]. These two existing methods are denoted as LG and L2 respectively. The computed errors E LG k 's and E L2 k 's and the rates of convergence for LG and L2 are also listed in Table 5.1, from which it is clear that both of the existing methods are 1st-order accurate, one order lower than our method.
To solve this problem, we choose a mesh with h = 0.02 and ∆t = 1/52. The numerical solution from our method for α = 1.5 is plotted in Figure 5.1 against t and the original independent variable S = e x . From the figure we see that the method is numerically very stable. To see the influence of α on the option price, we solve the problem for α = 1.3, 1.5, 1.7 and 2, and plot in Figure 5.2, the difference between the values from the fBS equation (i.e., α < 2) and the BS equation (α = 2), denoted respectively C f BS (x, t) and C BS (x, t), at t = 0. From this figure we see that the call price increases as α decreases when S is greater than a critical value. This is likely because when α is close to 1, the solution to the fBS equation exhibits jump (or convection) nature, while when α is close to 2, it is of mainly diffusion nature. As a result, an option on a stock of jump nature is more expensive than one of diffusion nature, similar to the case that an option price is an increasing function of the volatility.
The problem has been solved using the same mesh as that for Example 2, and the solution for α = 1.5 is depicted in Figure 5  This problem is solved by our numerical scheme on the uniform partition of the solution domain (ln(0.1), ln(100)) × (0, 1) in (x, t) with M = 100 and N = 104. The penalty parameter is chosen to be λ = 10 10 . The difference between the computed value of this American option with α = 1.5 and the lower bound U * is plotted in, Figure 5.5. To see the difference between the American and European put options, we have also plot the difference between the computed European option value from Example 3 and U * in Figure 5.5. From the figure, we see that the American option is more expensive than its European counterpart. Also, it can be seen that the value  of the American option is bounded below by U * , while the value of the European option falls below U * in a sub-region of the solution domain.
Finally, we plot in Figure 5.6 the differences between the American put values from the fBS model and that from the BS model at the cross-section t = 0. From the figure we see that the value of the American option is a decreasing function of α when S is greater than a critical value, as observed in the Examples 2 and 3.
6. Conclusion. In this paper, we constructed and analyzed a novel 2nd-order finite difference method for the penalized fractional Black-Scholes equation governing European and American option pricing. We have proved the convergence of numerical method by showing that the method is consistent, stable and monotone. In  particular, we have shown that the truncation error of the scheme is of 2nd-order as opposed to the 1st-order truncation errors for the existing numerical methods for the fBS equation. Numerical experiments have been carried out to verify the theoretical findings. The numerical results show that our method is 2nd-order accurate and gives practically useful and correct results when it is used for pricing European and American options.