Journal of Applied Mathematics

Journal of Applied Mathematics / 2012 / Article

Research Article | Open Access

Volume 2012 |Article ID 796814 | https://doi.org/10.1155/2012/796814

Zhongdi Cen, Anbo Le, Aimin Xu, "Exponential Time Integration and Second-Order Difference Scheme for a Generalized Black-Scholes Equation", Journal of Applied Mathematics, vol. 2012, Article ID 796814, 12 pages, 2012. https://doi.org/10.1155/2012/796814

Exponential Time Integration and Second-Order Difference Scheme for a Generalized Black-Scholes Equation

Academic Editor: Kai Diethelm
Received20 Sep 2011
Revised05 Dec 2011
Accepted05 Dec 2011
Published19 Feb 2012

Abstract

We apply an exponential time integration scheme combined with a central difference scheme on a piecewise uniform mesh with respect to the spatial variable to evaluate a generalized Black-Scholes equation. We show that the scheme is second-order convergent for both time and spatial variables. It is proved that the scheme is unconditionally stable. Numerical results support the theoretical results.

1. Introduction

The pricing and hedging of derivative securities, also known as contingent claims, is a subject of much practical importance. One basic type of derivative is an option. The owner of a call option has the right but not the obligation to purchase an underlying asset (such as a stock) for a specified price (called the exercise price or strike price) on or before a expiry date. A put option is similar except the owner of such a contract has the right but not the obligation to sell. Options which can be exercised only on the expiry date are called European, whereas options which can be exercised any time up to and including the expiry date are classified as American. It was shown by Black and Scholes [1] that these option prices satisfy a second-order partial differential equation with respect to the time horizon 𝑡 and the underlying asset price 𝑥. This equation is now known as the Black-Scholes equation and can be solved exactly when the coefficients are constant or space-independent. However, in many practical situations, numerical solutions are normally sought. Therefore, efficient and accurate numerical algorithms are essential for solving this problem accurately.

The Black-Scholes differential operator at 𝑥=0 is degenerate. A common and widely used approach by many authors dealing with finite difference/volume/element methods for the Black-Scholes equation is to apply an Euler transformation to remove the degeneracy of the differential operator when the parameters of the Black-Scholes equation are constant or space independent, see for example [2–7]. As a result of the Euler transformation, the transformed interval becomes (−∞,∞). However, the truncation on the left-hand side of the domain to artificially remove the degeneracy may cause computational errors. Furthermore, the uniform mesh on the transformed interval will lead to the originally grid points concentrating around 𝑥=0 inappropriately. Moreover, when a problem is space-dependent, this transformation is impossible, and thus the Black-Scholes equation in the original form needs to be solved.

It is well known that when using the standard finite difference method to solve those problems involving the convection-diffusion operator, such as the Black-Scholes differential operator, numerical difficulty can be caused. The main reason is that when the volatility or the asset price is small, the Black-Scholes differential operator becomes a convection-dominated operator. Hence, the implicit Euler scheme with central spatial difference method may lead to nonphysical oscillations in the computed solution. The implicit Euler scheme with upwind spatial difference method does not have this disadvantage, but this difference scheme is only first-order convergent [8]. Recently, a stable fitted finite volume method [9] is employed for the discretization of the Black-Scholes equation. But it is also first-order convergent. In [10, 11] numerical methods of option pricing models are studied by applying a standard finite volume method to obtain a difference scheme. Their numerical schemes use central difference for a given mesh, but switch to upstream weighting for a small number of nodes, which are second-order spatial convergent. In this paper we change the grid spacing to a piecewise uniform mesh which is constructed so that central difference is used everywhere.

For time discretization, explicit schemes are easy to implement but suffer from stability problems. Some well-known second-order implicit schemes, such as Crank-Nicolson method, are prone to spurious oscillations unless the time step size is no more than twice the maximum time step size for an explicit method (see Zvan et al. [10, 11]). Although the fully implicit backward Euler method may be used to accurately solve the Black-Scholes PDE due to its strong stability properties, it is only first-order accurate in time. Exponential time integration has gained importance following the work of Cox and Matthews [12] and with recent developments in efficient methods for computing the matrix exponential [13–16], this time evolution method is likely to be a popular choice for solving large semidiscrete systems arising in various numerical computations.

In [17, 18], we have presented robust difference schemes for the Black-Scholes equation, which is based on the implicit Euler method for time discretization and a central difference method for spatial discretization. In this paper, we apply an exponential time integration scheme combined with a central difference scheme on a piecewise uniform mesh with respect to the spatial variable. We show that our scheme is second-order convergent for both time and spatial variables, while in [17, 18] the convergence for the time variable is only first order. It is proved that the scheme is unconditionally stable. Numerical results support this conclusion.

The rest of the paper is organized as follows. In the next section we discuss the continuous model of the Black-Scholes equations. The discretization method is described in Section 3. It is shown that the finite difference scheme is second-order convergent with respect to both spatial and time variables. In Section 4 we prove that the difference scheme is unconditionally stable. Finally numerical examples are presented in Section 5.

2. The Continuous Problem

We consider the following generalized Black-Scholes equation𝜕𝑣−1𝜕𝑡2ğœŽ2(𝑥,𝑡)𝑥2𝜕2𝑣𝜕𝑥2−𝑟(𝑡)𝑥𝜕𝑣𝜕𝑥+𝑟(𝑡)𝑣=0,(𝑥,𝑡)∈ℝ+×(0,𝑇),(2.1)𝑣(𝑥,0)=max(𝑥−𝐾,0),𝑥∈ℝ+[].,(2.2)𝑣(0,𝑡)=0,𝑡∈0,𝑇(2.3) Here 𝑣(𝑥,𝑡) is the European call option price at asset price 𝑥 and at time to maturity 𝑡, 𝐾 is the exercise price, 𝑇 is the maturity, 𝑟(𝑡) is the risk-free interest rate, and ğœŽ(𝑥,𝑡) represents the volatility function of underlying asset. Here we assume that ğœŽ2≥𝛼>0, 𝛽∗≥𝑟≥𝛽>0. When ğœŽ and 𝑟 are constant functions, it becomes the classical Black-Scholes model.

Even though as 𝑥 goes to zero the generalized Black-Scholes operator (2.1) is degenerate, the existence and uniqueness of a solution of (2.1)–(2.3) is well known. Due to the fact that the initial condition is not smooth, the finite difference scheme may not converge to the exact solution.

We first modify the model as follows. Define 𝜋𝜀(𝑦) asğœ‹ğœ€âŽ§âŽªâŽ¨âŽªâŽ©ğ‘(𝑦)=𝑦,𝑦≥𝜀,0+𝑐1𝑦+⋯+𝑐9𝑦9,−𝜀<𝑦<𝜀,0,𝑦≤−𝜀,(2.4) where 0<𝜀≪1 is a transition parameter and 𝜋𝜀(𝑦) is a function which smooths out the original max(𝑦,0) around 𝑦=0. This requires that 𝜋𝜀(𝑦) satisfies𝜋𝜀(−𝜀)=ğœ‹î…žğœ€(−𝜀)=ğœ‹ğœ€î…žî…ž(−𝜀)=ğœ‹ğœ€î…žî…žî…ž(−𝜀)=𝜋𝜀(4)𝜋(−𝜀)=0,𝜀(𝜀)=𝜀,ğœ‹î…žğœ€(𝜀)=1,ğœ‹ğœ€î…žî…ž(𝜀)=ğœ‹ğœ€î…žî…žî…ž(𝜀)=𝜋𝜀(4)(𝜀)=0.(2.5) Using these ten conditions we can easily find that𝑐0=35256𝜀,𝑐1=12,𝑐2=3564𝜀,𝑐4=−35128𝜀3,𝑐6=764𝜀5,𝑐85=−256𝜀7,𝑐3=𝑐5=𝑐7=𝑐9=0.(2.6) Replacing max(𝑥−𝐾,0) in the initial condition (2.2) by the fourth-order smooth function 𝜋𝜀(𝑥−𝐾) we obtain𝜕𝑤−1𝜕𝑡2ğœŽ2(𝑥,𝑡)𝑥2𝜕2𝑤𝜕𝑥2−𝑟(𝑡)𝑥𝜕𝑤𝜕𝑥+𝑟(𝑡)𝑤=0,(𝑥,𝑡)∈ℝ+×(0,𝑇),(2.7)𝑤(𝑥,0)=𝜋𝜀(𝑥−𝐾),𝑥∈ℝ+[].,(2.8)𝑤(0,𝑡)=0,𝑡∈0,𝑇(2.9) The existence and uniqueness of a solution of (2.7)–(2.9) can be found in [19], which also gives the following result: for (𝑥,𝑡)∈[0,+∞)×[0,𝑇],||||‖‖𝜋𝑤(𝑥,𝑡)−𝑣(𝑥,𝑡)≤𝐶𝜀‖‖(𝑥−𝐾)−max(𝑥−𝐾,0)ğ¿âˆž,(2.10) where 𝐶 is a positive constant.

In order to apply the numerical method we need to truncate the infinite domain [0,+∞)×[0,𝑇] into Ω=[0,𝑆max]×[0,𝑇], where 𝑆max is suitably chosen positive number. Based on Wilmott et al.’s estimate [20] that the upper bound of the asset price 𝑆max is typically three or four times the strike price, it is reasonable for us to set 𝑆max=4𝐾. Thus we consider the following problem:𝜕𝑢−1𝜕𝑡2ğœŽ2(𝑥,𝑡)𝑥2𝜕2𝑢𝜕𝑥2−𝑟(𝑡)𝑥𝜕𝑢𝜕𝑥+𝑟(𝑡)𝑢=0,(𝑥,𝑡)∈Ω,(2.11)𝑢(𝑥,0)=𝜋𝜀(𝑥−𝐾),𝑥∈0,𝑆max[]𝑢𝑆,(2.12)𝑢(0,𝑡)=0,𝑡∈0,𝑇,(2.13)max,𝑡=𝑆max−𝐾𝑒−∫𝑇𝑇−𝑡𝑟(𝑠)d𝑠[].,𝑡∈0,𝑇(2.14) The existence and uniqueness of a solution of (2.11)–(2.14) can be also found in [19], which also gives the following result:||||−𝑆𝑢(𝑥,𝑡)−𝑤(𝑥,𝑡)≤𝐾explnmax/𝑥22minÎ©ğœŽ2(𝑇−𝑡),(𝑥,𝑡)∈Ω.(2.15) It follows from (2.10) and (2.15) that we can make the solution of our modified model (2.11)–(2.14) close to that of the original model (2.1)–(2.3) by choosing sufficiently small 𝜀 and sufficiently large 𝑆max. In the remaining of this paper we will consider the model (2.11)–(2.14).

3. Discretization Scheme

The exponential time differencing scheme is an approach used for the numerical solution of a wide range of PDEs that involve spatial variables as well as a time variable. The spatial discretization results in an approximating system of ODEs. For an inhomogeneous linear PDE, this method leads to an inhomogeneous linear system of ODEs in time whose solution satisfies a two-term recurrence relation involving the matrix exponential, where the matrix is determined by the form of spatial discretization applied such as finite difference, finite element, or spectral approach. In this paper we use a central difference scheme on a piecewise uniform mesh for approximating the spatial derivatives.

The use of central difference scheme on the uniform mesh may produces nonphysical oscillations in the computed solution. To overcome this oscillation we use a piecewise uniform mesh Ω𝑁+1 on the space interval [0,𝑆max]:𝑥𝑖=âŽ§âŽªâŽªâŽªâŽªâŽªâŽ¨âŽªâŽªâŽªâŽªâŽªâŽ©â„Žî‚¸ğ›¼â„Žğ‘–=1,1+𝛽∗𝑁(𝑖−1)𝑖=2,…,4𝑁−1,𝐾𝑖=4,𝑁𝐾+𝜀𝑖=4𝑆+1,𝐾+𝜀+max−𝐾−𝜀𝑁3𝑁/4−1(𝐼−𝑁/4−1)𝑖=4+2,…,𝑁,(3.1) whereℎ=𝐾−𝜀1+𝛼/𝛽∗.(𝑁/4−2)(3.2) Here we have used a refined mesh at the region near 𝑥=𝐾 for treating the nonsmoothness of the payoff function. It is easy to see that the mesh sizes â„Žğ‘–=𝑥𝑖−𝑥𝑖−1 satisfyâ„Žğ‘–=âŽ§âŽªâŽªâŽªâŽªâŽ¨âŽªâŽªâŽªâŽªâŽ©ğ›¼â„Žğ‘–=1,ğ›½âˆ—ğ‘â„Žğ‘–=2,…,4𝑁−1,𝜀𝑖=4,𝑁4𝑆+1,max−𝐾−𝜀𝑁3𝑁/4−1𝑖=4+2,…,𝑁.(3.3)

We discretize the generalized Black-Scholes operator using a central difference scheme on the above piecewise uniform mesh:𝐿𝑁𝑈𝑖(𝑡)=𝑑𝑈𝑖(𝑡)âˆ’ğœŽğ‘‘ğ‘¡2𝑖(𝑡)𝑥2ğ‘–â„Žğ‘–+â„Žğ‘–+1𝑈𝑖+1(𝑡)−𝑈𝑖(𝑡)â„Žğ‘–+1−𝑈𝑖(𝑡)−𝑈𝑖−1(𝑡)â„Žğ‘–î‚¶âˆ’ğ‘Ÿ(𝑡)𝑥𝑖𝑈𝑖+1(𝑡)−𝑈𝑖−1(𝑡)â„Žğ‘–+â„Žğ‘–+1+𝑟(𝑡)𝑈𝑖(𝑡)(3.4) for 𝑖=1,…,𝑁−1. This discretization leads to an initial value problem of the form𝑑𝐔𝑑𝑡=𝐀(𝑡)𝐔(𝑡)+𝐟(𝑡),𝐔(0)=𝜋𝜀(𝐱−𝐊),(3.5) where 𝐔(𝑡)=(𝑈1(𝑡),…,𝑈𝑁−1(𝑡))𝑇, the matrix 𝐀(𝑡) of order (𝑁−1) is given byâŽ¡âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ£ğ‘ğ€(𝑡)=1𝑐1ğ‘Žâ‹¯02𝑏2𝑐2â‹¯â‹®â‹®ğ‘Ž3𝑏3𝑐3â‹¯â‹®â‹®â‹®â‹®â‹®â‹¯ğ‘Žğ‘âˆ’2𝑏𝑁−2𝑐𝑁−20â‹¯ğ‘Žğ‘âˆ’1𝑏𝑁−1⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦(3.6) withğ‘Žğ‘–ğœŽ(𝑡)=2𝑖(𝑡)𝑥2ğ‘–î€·â„Žğ‘–+â„Žğ‘–+1î€¸â„Žğ‘–âˆ’ğ‘Ÿ(𝑡)ğ‘¥ğ‘–â„Žğ‘–+â„Žğ‘–+1,ğ‘ğ‘–ğœŽ(𝑡)=−2𝑖(𝑡)𝑥2ğ‘–â„Žğ‘–â„Žğ‘–+1𝑐−𝑟(𝑡),ğ‘–ğœŽ(𝑡)=2𝑖(𝑡)𝑥2ğ‘–î€·â„Žğ‘–+â„Žğ‘–+1î€¸â„Žğ‘–+1+𝑟(𝑡)ğ‘¥ğ‘–â„Žğ‘–+â„Žğ‘–+1for𝑖=1,…,𝑁−1.(3.7) The vectors 𝐟(𝑡) and 𝜋𝜀(𝐱−𝐊) are the corresponding boundary and initial conditions:âŽ›âŽœâŽœâŽœâŽœâŽœâŽœâŽœâŽœâŽğ‘ŽğŸ(𝑡)=1𝑈00⋮0𝑐(𝑡)𝑁−1ğ‘ˆğ‘âŽžâŽŸâŽŸâŽŸâŽŸâŽŸâŽŸâŽŸâŽŸâŽ (𝑡),𝜋𝜀(âŽ›âŽœâŽœâŽœâŽœâŽœâŽœâŽğœ‹ğ±âˆ’ğŠ)=𝜀𝑥1𝜋−𝐾𝜀𝑥2⋮𝜋−𝐾𝜀𝑥𝑁−1î€¸âŽžâŽŸâŽŸâŽŸâŽŸâŽŸâŽŸâŽ âˆ’ğ¾.(3.8)

Lemma 3.1. For each 𝑡 the matrix 𝐀(𝑡) is strictly diagonally dominant.

Proof. It is easy to see that ğ‘Žğ‘–ğœŽ(𝑡)>2𝑖(𝑡)𝑥1ğ‘¥ğ‘–î€·â„Žğ‘–+â„Žğ‘–+1î€¸â„Žğ‘–âˆ’ğ‘Ÿ(𝑡)ğ‘¥ğ‘–â„Žğ‘–+â„Žğ‘–+1≥𝛼𝑥1âˆ’ğ›½âˆ—â„Žğ‘–î€¸ğ‘¥ğ‘–î€·â„Žğ‘–+â„Žğ‘–+1î€¸â„Žğ‘–=î€·âˆ’ğ›¼â„Ž+𝛽∗𝛼/ğ›½âˆ—î€¸â„Žî€¸ğ‘¥ğ‘–î€·â„Žğ‘–+â„Žğ‘–+1î€¸â„Žğ‘–ğ‘=0,2≤𝑖<4ğ‘Žğ‘–î€·(𝑡)â‰¥ğ›¼ğ‘¥ğ‘–âˆ’ğ›½âˆ—â„Žğ‘–î€¸ğ‘¥ğ‘–î€·â„Žğ‘–+â„Žğ‘–+1î€¸â„Žğ‘–ğ‘>0,4≤𝑖≤𝑁−1,(3.9) for sufficiently large 𝑁. Clearly, 𝑏𝑖(𝑡)<0for1≤𝑖≤𝑁−1,𝑐𝑖𝑏(𝑡)>0for1≤𝑖≤𝑁−2,1(𝑡)+𝑐1ğ‘Ž(𝑡)<0,𝑖(𝑡)+𝑏𝑖(𝑡)+ğ‘ğ‘–ğ‘Ž(𝑡)<0,2≤𝑖≤𝑁−2,𝑁−1(𝑡)+𝑏𝑁−1(𝑡)<0.(3.10) Hence we verify that the matrix 𝐀(𝑡) is strictly diagonally dominant.

Solving the system of (3.5) gives𝐔(𝑡)=𝑒∫𝑡0𝐀(𝑠)d𝑠𝐔(0)+𝑡0𝑒−∫𝑠0𝐀(𝑦)d𝑦,𝐟(𝑠)d𝑠(3.11) which satisfies the following recurrence relation:𝐔(𝑡+𝑙)=𝑒∫𝑡𝑡+𝑙𝐀(𝑠)d𝑠𝐔(𝑡)+𝑡𝑡+𝑙𝑒−∫𝑠𝑡𝐀(𝑦)d𝑦𝐟(𝑠)d𝑠,𝑡=0,𝑙,2𝑙,…,(3.12) in which 𝑙 is a constant time step in the discretization of the time variable 𝑡≥0, at the points 𝑡𝑗=𝑗𝑙 (𝑗=0,1,2,…,𝑀).

Using a numerical integration rule for evaluating the quadrature in (3.12) gives the following second-order formula:𝐔(𝑡+𝑙)≈𝑒𝑙𝐀(𝑡)𝐔(𝑡)+𝑡𝑡+𝑙𝑒−(𝑠−𝑡)𝐀(𝑡)𝐟(𝑠)d𝑠=𝑒𝑙𝐀(𝑡)𝐔(𝑡)+𝑡𝑡+𝑙𝑒(𝑡+𝑙−𝑠)𝐀(𝑡)𝐟(𝑠)d𝑠≈𝑒𝑙𝐀(𝑡)𝐔(𝑡)+𝑡𝑡+𝑙𝑒(𝑡+𝑙−𝑠)𝐀(𝑡)𝑙−1(𝑡+𝑙−𝑠)𝐟(𝑡)+𝑙−1(𝑠−𝑡)𝐟(𝑡+𝑙)d𝑠=𝑒𝑙𝐀(𝑡)𝐔(𝑡)+(𝑙𝐀(𝑡))−1𝑙𝑒𝑙𝐀(𝑡)+𝐀−1(𝑡)−𝐀−1(𝑡)𝑒𝑙𝐀(𝑡)𝐟(𝑡)+(𝑙𝐀(𝑡))−1𝐀−1(𝑡)𝑒𝑙𝐀(𝑡)−𝑙𝐼−𝐀−1(𝑡)𝐟(𝑡+𝑙),𝑡=0,𝑙,2𝑙,….(3.13)

Now the problem is how to approximate 𝑒𝑙𝐀(𝑡) to get numerical solution. A good approximation to 𝑒𝑧 is the (𝑏,𝑑) Padé approximation which has the form [21, 22]𝑒𝑧≈𝑅𝑏,𝑑𝑃(𝑧)=𝑑(𝑧)𝑄𝑏,(𝑧)(3.14) where 𝑃𝑑(𝑧) and 𝑄𝑏(𝑧) are the polynomials of degrees 𝑑 and 𝑏, respectively, with real coefficients, in each of which the constant term is unity.

The numerical method to be employed here is based on the use of the following second-order rational approximation:𝑒𝑧≈1+(1−𝑐)𝑧1−𝑐𝑧+(𝑐−1/2)𝑧2.(3.15) So we have1𝑅(𝑙𝐀(𝑡))=𝐼−𝑐𝑙𝐀(𝑡)+𝑐−2𝑙2𝐀2(𝑡)−1[]𝐼+(1−𝑐)𝑙𝐀(𝑡),(3.16) for the matrix exponentials in our recurrence relation (3.13).

Using (3.13) and (3.16) we get𝑙𝑈(𝑡+𝑙)=𝑅(𝑙𝐀(𝑡))𝑈(𝑡)+2[]𝑉(𝑙𝐀(𝑡))𝐟(𝑡)+𝑊(𝑙𝐀(𝑡))𝐟(𝑡+𝑙)(3.17) for 𝑡=0,𝑙,2𝑙,…,𝑀, where1𝑉(𝑙𝐀(𝑡))=𝐼−𝑐𝑙𝐀(𝑡)+𝑐−2𝑙2𝐀2(𝑡)−1,𝑊1(𝑙𝐀(𝑡))=𝐼−𝑐𝑙𝐀(𝑡)+𝑐−2𝑙2𝐀2(𝑡)−11𝐼−2𝑐−2.𝑙𝐀(𝑡)(3.18)

It can be seen that the truncation error of the difference scheme (3.17) with √1/2<𝑐<2−2 is 𝑂(ℎ2+𝑙2) (see [21], e.g.).

4. Stability Analysis

Lemma 4.1. For each 𝑡 the real part of each nonzero eigenvalue of 𝐀(𝑡) is negative.

Proof. From Lemma 3.1 we can obtain that for each 𝑡 the matrix −𝐀(𝑡) is an M-matrix, thus the real part of each nonzero eigenvalue of −𝐀(𝑡) is positive (see [23, Theorem  3.1]). Hence the real part of each nonzero eigenvalue of 𝐀(𝑡) is negative.

Theorem 4.2. The difference scheme (3.17) is unconditionally stable.

Proof. Let 𝜆𝑖 (𝑖=1,2,…,𝑁−1) be eigenvalues of matrix 𝐀(𝑡) for a fixed 𝑡, then 1+(1−𝑐)𝑙𝜆𝑖1−𝑐𝑙𝜆𝑖+(𝑐−1/2)𝑙2𝜆2𝑖(4.1) are eigenvalues of matrix 1𝐼−𝑐𝑙𝐀(𝑡)+𝑐−2𝑙2𝐀2(𝑡)−1[].𝐼+(1−𝑐)𝑙𝐀(𝑡)(4.2) Let 𝜆𝑖 is the conjugate complex of 𝜆𝑖. It is easy to check 1+(1−𝑐)𝑙𝜆𝑖⋅1+(1−𝑐)𝑙𝜆𝑖−1−𝑐𝑙𝜆𝑖+1𝑐−2𝑙2𝜆2𝑖⋅1−𝑐𝑙𝜆𝑖+1𝑐−2𝑙2𝜆2𝑖<0,(4.3) where we have used Lemma 4.1. Hence we have ||||1+(1−𝑐)𝑙𝜆𝑖1−𝑐𝑙𝜆𝑖+(𝑐−1/2)𝑙2𝜆2𝑖||||<1.(4.4) From this we complete the proof.

5. Numerical Experiments

In this section we verify experimentally the theoretical results obtained in the preceding section. Errors and convergence rates for the numerical scheme are presented for two test problems.

Test 1
European call option with parameters: ğœŽ=0.2+0.2(1−𝑡)((𝑥/25−1.2)2/((𝑥/25)2+1.44))), 𝑟=0.06, 𝑇=1, 𝐾=25 and 𝑆max=100.

Test 2
European call option with parameters: ğœŽ=0.2[1+0.1(1−𝑡)(𝑥/(1+𝑥))], 𝑟=0.06, 𝑇=1, 𝐾=25 and 𝑆max=100.

For Tests 1 and 2 we choose √𝜀=0.0001,𝑐=(5/2−2)/2 of the rational approximation (3.17). The computed option value 𝑈 with 𝑁=𝑀=64 are depicted in Figures 1 and 2, respectively.

The exact solutions of our test problems are not available. We use the approximated solution obtained by the implicit Euler method in [17] with 𝑁=2048, 𝑀=2048 as the exact solution since we know this method converges. Because we only know "the exact solution" on mesh points, we use the linear interpolation to get solutions at other points. In this paper 𝑈(𝑥,𝑡) denotes “the exact solution” which is a linear interpolation of the approximated solution 𝑈2048,2048 obtained by the implicit Euler method. We measure the accuracy in the discrete maximum norm𝑒𝑁,𝑀=max𝑖,𝑗||𝑈𝑁,𝑀𝑖𝑗−𝑈𝑥𝑖,𝑡𝑗||,(5.1) and the convergence rate𝑟𝑁,𝑀=log2𝑒𝑁,𝑀𝑒2𝑁,2𝑀(5.2) for the exponential time integration method and the implicit Euler method in [17]. The error estimates and convergence rates in our computed solutions of Tests 1 and 2 from both methods are listed in Tables 1 and 2, respectively.


MN Exponential methodImplicit Euler method
Error 𝑒 𝑁 , 𝑀 Rate 𝑟 𝑁 , 𝑀 Error 𝑒 𝑁 , 𝑀 Rate 𝑟 𝑁 , 𝑀

16 64 1 . 2 5 3 5 𝑒 − 1 — 1 . 7 8 1 7 𝑒 − 1 —
32 128 2 . 9 2 6 8 𝑒 − 2 2.099 8 . 9 5 6 7 𝑒 − 2 0.992
64 256 1 . 5 7 2 5 𝑒 − 2 0.896 4 . 6 8 2 2 𝑒 − 2 0.936


MNExponential methodImplicit Euler method
Error 𝑒 𝑁 , 𝑀 Rate 𝑟 𝑁 , 𝑀 Error 𝑒 𝑁 , 𝑀 Rate 𝑟 𝑁 , 𝑀

16 64 1 . 0 7 1 6 𝑒 − 1 — 1 . 7 4 2 8 𝑒 − 1 —
32 128 2 . 4 7 1 6 𝑒 − 2 2.116 8 . 9 5 0 4 𝑒 − 2 0.961
64 256 1 . 5 8 1 0 𝑒 − 2 0.645 4 . 7 7 8 0 𝑒 − 2 0.906

From the figures it is seen that the numerical solutions by our method are nonoscillatory. We compare the accuracy of the exponential time integration method with the implicit Euler method in [17]. From Tables 1 and 2 we can see the exponential time integration method converges more rapidly than the implicit Euler method in [17], but for the same 𝑀 and 𝑁 the exponential time integration method require more computer time than the implicit Euler method. We also can see that 𝑟𝑁,𝑀 of the exponential time integration method is close to 2, even though this rate is not reached for large 𝑀 and 𝑁, the reason is that the “exact solution” is not really the exact solution. Therefore the numerical results support the convergence estimate of Theorem 4.2.

Acknowledgments

The authors would like to thank the anonymous referee for several suggestions for the improvement of this paper. The work was supported by Zhejiang Province Natural Science Foundation (Grant nos. Y6100021, Y6110310) of China and Ningbo Municipal Natural Science Foundation (Grant no. 2010A610099) of China.

References

  1. F. Black and M. S. Scholes, “The pricing of options and corporate liabilities,” Journal of Political Economy, vol. 81, pp. 637–654, 1973. View at: Google Scholar
  2. G. Courtadon, “A more accurate finite difference approximation for the valuation of options,” Journal of Financial and Quantitative Analysis, vol. 17, pp. 697–703, 1982. View at: Google Scholar
  3. M. Giles and R. Carter, “Convergence analysis of Crank-Nicolson and Rannacher time-marching,” Journal of Computational Finance, vol. 94, pp. 89–112, 2006. View at: Google Scholar
  4. D. M. Pooley, K. R. Vetzal, and P. A. Forsyth, “Convergence remedies for non-smooth payoffs in option pricing,” Journal of Computational Finance, vol. 6, no. 4, pp. 25–40, 2003. View at: Google Scholar
  5. E. Schwartz, “The valuation of warrants: implementing a new approach,” The Journal of Financial Economics, vol. 13, pp. 79–93, 1977. View at: Google Scholar
  6. D. Y. Tangman, A. Gopaul, and M. Bhuruth, “Numerical pricing of options using high-order compact finite difference schemes,” Journal of Computational and Applied Mathematics, vol. 218, no. 2, pp. 270–280, 2008. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  7. J. Zhao, M. Davison, and R. M. Corless, “Compact finite difference method for American option pricing,” Journal of Computational and Applied Mathematics, vol. 206, no. 1, pp. 306–321, 2007. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  8. C. Vázquez, “An upwind numerical approach for an American and European option pricing model,” Applied Mathematics and Computation, vol. 97, no. 2-3, pp. 273–286, 1998. View at: Publisher Site | Google Scholar
  9. S. Wang, “A novel fitted finite volume method for the Black-Scholes equation governing option pricing,” IMA Journal of Numerical Analysis, vol. 24, no. 4, pp. 699–720, 2004. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  10. R. Zvan, P. Forsyth, and K. Vetzal, “Robust numerical methods for PDE models of Asian options,” The Journal of Financial Economics, vol. 1, no. 2, pp. 39–78, 1998. View at: Google Scholar
  11. R. Zvan, P. A. Forsyth, and K. R. Vetzal, “A finite volume approach for contingent claims valuation,” IMA Journal of Numerical Analysis, vol. 21, no. 3, pp. 703–731, 2001. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  12. S. M. Cox and P. C. Matthews, “Exponential time differencing for stiff systems,” Journal of Computational Physics, vol. 176, no. 2, pp. 430–455, 2002. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  13. N. Rambeerich, D. Y. Tangman, A. Gopaul, and M. Bhuruth, “Exponential time integration for fast finite element solutions of some financial engineering problems,” Journal of Computational and Applied Mathematics, vol. 224, no. 2, pp. 668–678, 2009. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  14. D. Y. Tangman, A. Gopaul, and M. Bhuruth, “Exponential time integration and Chebychev discretisation schemes for fast pricing of options,” Applied Numerical Mathematics, vol. 58, no. 9, pp. 1309–1319, 2008. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  15. Z. F. Tian and P. X. Yu, “A high-order exponential scheme for solving 1D unsteady convection-diffusion equations,” Journal of Computational and Applied Mathematics, vol. 235, no. 8, pp. 2477–2491, 2011. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  16. M. Yousuf, “Efficient L-stable method for parabolic problems with application to pricing American options under stochastic volatility,” Applied Mathematics and Computation, vol. 213, no. 1, pp. 121–136, 2009. View at: Publisher Site | Google Scholar
  17. Z. Cen and A. Le, “A robust and accurate finite difference method for a generalized Black-Scholes equation,” Journal of Computational and Applied Mathematics, vol. 235, no. 13, pp. 3728–2733, 2011. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  18. Z. Cen, A. Le, and A. Xu, “A second-order difference scheme for the penalized Black-Scholes equation governing American put option pricing,” Computational Economics. In press. View at: Publisher Site | Google Scholar
  19. C.-K. Cho, T. Kim, and Y. Kwon, “Estimation of local volatilities in a generalized Black-Scholes model,” Applied Mathematics and Computation, vol. 162, no. 3, pp. 1135–1149, 2005. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  20. P. Wilmott, J. Dewynne, and S. Howison, Option Pricing: Mathematical Models and Computation, Oxford Financial Press, Oxford, UK, 1993.
  21. M. Dehghan, “Numerical techniques for a parabolic equation subject to an overspecified boundary condition,” Applied Mathematics and Computation, vol. 132, no. 2-3, pp. 299–313, 2002. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  22. G. D. Smith, Numerical Solution of Partial Differential Equations (Finite Difference Method), The Clarendon Press Oxford University Press, New York, NY, USA, 3rd edition, 1985.
  23. G. Poole and T. Boullion, “A survey on M-matrices,” SIAM Review, vol. 16, pp. 419–427, 1974. View at: Publisher Site | Google Scholar

Copyright © 2012 Zhongdi Cen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder
Views1455
Downloads633
Citations

Related articles

We are committed to sharing findings related to COVID-19 as quickly as possible. We will be providing unlimited waivers of publication charges for accepted research articles as well as case reports and case series related to COVID-19. Review articles are excluded from this waiver policy. Sign up here as a reviewer to help fast-track new submissions.