A high-order finite difference method for option valuation

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A high-order finite difference method for option valuation Mehzabeen Jumanah Dilloo, Désiré Yannick Tangman * Department of Mathematics, University of Mauritius, Réduit, Mauritius

article

info

Article history: Received 11 April 2016 Received in revised form 26 April 2017 Accepted 6 May 2017 Available online xxxx Keywords: High-order scheme Local mesh refinement Exponential time integration Merton’s jump–diffusion model Heston’s stochastic volatility model Nonlinear Black–Scholes equation

a b s t r a c t In this paper, we propose the use of an efficient high-order finite difference algorithm to price options under several pricing models including the Black–Scholes model, the Merton’s jump–diffusion model, the Heston’s stochastic volatility model and the nonlinear transaction costs or illiquidity models. We apply a local mesh refinement strategy at the points of singularity usually found in the payoff of most financial derivatives to improve the accuracy and restore the rate of convergence of a non-uniform high-order five-point stencil finite difference scheme. For linear models, the time-stepping is dealt with by using an exponential time integration scheme with Carathéodory–Fejér approximations to efficiently evaluate the product of a matrix exponential with a vector of option prices. Nonlinear Black–Scholes equations are solved using an efficient iterative scheme coupled with a Richardson extrapolation. Our numerical experiments clearly demonstrate the highorder accuracy of the proposed finite difference method, making the latter a method of choice for solving both linear and nonlinear partial differential equations in financial engineering problems. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction In financial engineering, several methods have been proposed in literature to improve the accuracy of numerical option pricing algorithms. Low-order finite difference schemes for solving partial differential equations (PDEs) have been enhanced through the common use of extrapolation techniques [1]. Unfortunately, extrapolation is highly dependent on the smooth convergence of the numerical methods underhand and can give inaccurate results if the convergence rate is too erratic. In [2–4], the authors have used spectral methods to improve the accuracy of European style options under one-dimensional models. However, these methods require the solutions of dense linear systems in RN ×N which have O(N 3 ) complexity, where N represents the number of nodes used to discretise the spatial domain. Although such schemes give accurate results for European options, they are less efficient for American options because spectral convergence is usually not achieved due to the difficulty in accurately capturing the free boundary location [5]. On the other hand, [6–8] have developed high-order compact (HOC) schemes which have linear complexity when computing the option prices. However, the derivation of these schemes involves the differentiation and intricate manipulations of the PDE at hand in order to replace the higher order derivatives in the leading truncation error terms [6]. Such an approach will be generally more difficult to implement for PDEs with variable and time dependent coefficients such as those that arise in the nonlinear Black–Scholes equations. Moreover, it is well known that the point of singularity that exists at the strike price hampers the order of convergence of these HOC discretisations. So far, one of the best remedial approaches was suggested in [7] where a local mesh refinement strategy is applied by adding four nodes in the vicinity of that point of singularity at each refinement stage.

*

Corresponding author. E-mail address: [email protected] (D.Y. Tangman).

http://dx.doi.org/10.1016/j.camwa.2017.05.006 0898-1221/© 2017 Elsevier Ltd. All rights reserved.

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

2

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

When the assumptions of no transaction costs and a perfectly liquid market are relaxed, the Black–Scholes pricing PDE becomes nonlinear. Such problems do not admit closed-form solutions even for vanilla options, and since the nonlinear term is found in the diffusion coefficient, the problem becomes stiff. Numerically, this requires stability restrictions to be imposed on the time step if an explicit scheme is used to discretise this nonlinear term [9]. In [10], the second-order derivative within the nonlinear coefficient is approximated over twice the spatial mesh size used to discretise the derivatives with linear coefficients. Although this strategy helps to improve the stability of the implicit–explicit type of schemes, it does not completely alleviate the problem. More recently, splitting methods based on locally one-dimensional (LOD) backward Euler (implicit) method [11] and on LOD Crank–Nicolson method [12] have been suggested to solve the nonlinear Black– Scholes equations. Even though explicit formulae could be obtain while implementing the stable first-order implicit and second-order Crank–Nicolson schemes, restrictive positivity-preserving conditions on the time step are required to obtain convergent solutions to the exact solution. We also want to mention that [13,14] have proposed high-order schemes for the numerical solution of the nonlinear Black–Scholes equation, but since their methods did not cater for the singularities in the payoff functions, only a lower order convergence could be observed. Stochastic volatility models have been developed in order to produce option prices that better describe the volatility smiles and skews of real market data, and one of the most popular stochastic volatility model is the Heston’s model [15]. Several numerical methods, namely finite element-finite volume methods [16,17], spectral element approximations [18] and method of line approaches [19] among others, have been proposed to price options under this model. The resulting twodimensionality of the problem and the presence of second-order mixed spatial partial derivatives complicate the application of high-order schemes in the pricing of financial derivatives. Indeed, in [20], the derivation of a HOC scheme for the Heston’s model imposes the mesh size along each spatial direction to be equal, which can be quite restrictive. More recently, [5,21] have used the grid transformation of [22] which concentrates grid nodes at the strike price to build high-order schemes for stochastic volatility models. But this requires the evaluation of the Jacobian and the Hessian of transformations which are not very practical for multiple points of singularity. In this paper, we propose to use the non-uniform discretisation due to Bowen and Smith [23] for option pricing since such discretisation confers the following advantages over classical spatial discretisation as it:

• grants us with added flexibility to construct our grid such that refinement points can be easily added in the vicinity of several singular points,

• helps to reduce the number of computational nodes used, whereby far-field boundary conditions [24] can be naturally implemented using a coarse grid to extend our computational domain without the use of complex artificial boundary conditions [25–27], and this significantly improves the efficiency and practicality of the PDE approach, • results in banded pentadiagonal linear systems which have the same linear computational complexity as for tridiagonal linear systems, • gives high-order convergent prices for European options, path dependent options such as barrier options and American options, • gives high-order convergence for the nonlinear Black–Scholes and the Heston’s two-dimensional problems, and these of course, represent the main contributions of the present work. We also show how to solve the integro term in the Merton’s jump–diffusion model in O(N) computations per time step improving on the O(N log N) in [7,28]. Recently, [29] have used a second-order finite difference discretisation in space with a discontinuous Galerkin finite element in time to solve Merton’s partial integro-differential equation. However, the multigrid algorithm used to solve the dense linear systems is less efficient compared to our proposed scheme which solves only sparse linear systems. Our non-uniform pentadiagonal scheme can be easily combined with an exact in time exponential time integration scheme [2,3,30,31] to generate highly accurate option prices. An efficient iterative scheme [32] coupled with a repeated extrapolation technique [33] is further proposed to solve the nonlinear problems. The paper is structured as follows. In Section 2, we describe the different option pricing problems and models that we later consider in our numerical experiments. In Section 3, we give the high-order non-uniform finite difference approximations used in space and we explain the local mesh refinement strategy used to restore the high-order convergence rate. We then describe the time-stepping schemes for both linear and nonlinear problems, and we briefly discuss the stability analysis of the proposed scheme. Section 4 numerically demonstrates the efficiency and accuracy of the new scheme for option pricing, and finally, we conclude in Section 5 where we also give some promising scope for future work. 2. Option pricing We state here the different pricing problems and models that we shall study in this paper. Under the classical Black– Scholes model, we aim at pricing various types of options which have payoffs involving different points of singularity. This is to demonstrate the versatility of our approach. Then, we consider various types of nonlinear Black–Scholes PDEs arising from popular transaction costs and illiquidity models. Further, we show how pricing can be extended to other realistic models like the Merton’s jump–diffusion and Heston’s stochastic volatility models. Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

3

2.1. The Black–Scholes model Under the Black–Scholes model [34], options denoted by V (t ; S) satisfy the celebrated Black–Scholes PDE 1 ∂ 2V ∂V ∂V + σ 2 S 2 2 + rS − rV = 0, ∂t 2 ∂S ∂S where S is the asset price, σ is the variance of the return on the asset and r is the risk free interest rate.

(1)

Under this model, European and digital call options satisfy the Black–Scholes PDE (1) with terminal conditions V (T ; S) = max(S − E , 0),

(2)

and V (T ; S) = B H(S − E), respectively. Both payoffs admit a singularity at the strike E. H(·) is the Heaviside function and B is the notional amount paid if S > E at expiry. Barrier options are path dependent options that come into existence or are terminated when the asset price S breaches a certain barrier level H. In this paper, we consider an up-and-out barrier option which satisfies terminal condition (2) and an added condition that V (t ; S ≥ H) = 0 where t refers to the discrete or continuous monitoring dates. We further consider a butterfly spread option whose payoff has three points of singularity and is given as V (T ; S) = max(S − E1 , 0) − 2 max(S − E2 , 0) + max(S − E3 , 0), where E1 < E2 < E3 such that E2 = (E1 + E3 )/2. We refer the reader to [4,35] for more details about these options. Except for the continuous up-and-out barrier option which satisfies homogeneous Dirichlet boundary condition at the barrier H, Neumann boundary conditions can be implemented as

∂ 2V ≈ 0 as S → 0 and S → +∞, ∂ S2

(3)

for all the options discussed so far since they have linear payoff functions in S (i.e. ∂ V (T ; S)/∂ S ≈ constant) as S → 0 and S → +∞ [36]. Closed-form solutions under the Black–Scholes model [35] are also available for these options and we shall later use them to assess the efficiency of our numerical method. 2.1.1. American put option pricing problem There exists no easily computable analytical price for American options due to their distinctive early exercise feature. The valuation of such a financial contract first necessitates the solution of a time dependent free boundary location, and this problem can be posed as the free boundary value problem: 1 ∂ 2V ∂V ∂V + σ 2 S 2 2 + rS − rV = 0, ∂t 2 ∂S ∂S V (T ; S) = max(E − S , 0), V (t ; Sf (t)) = E − Sf (t), ∂ 2V ≈ 0, ∂ S2

S ∈ [Sf (t), +∞) and t ∈ [0, T ], S ≥ 0, t ∈ [0, T ], as S → +∞,

where Sf (t) is the optimal early exercise boundary. ˆ discrete monitoring dates will converge to the It is well known that a Bermudan option which can be early exercised at M ˆ → ∞. Away from monitoring dates, we solve the usual European style options price of an American option in the limit M and at monitoring dates tmˆ , we simply impose the payoff constraint condition V (tmˆ ; S) = max(V (tmˆ ; S), E − S),

ˆ. ˆ = 1, 2, . . ., M m

(4)

Since the first-order convergence of Bermudan options to American options with respect to the number of monitoring dates ˆ has been well established, Richardson extrapolation can be applied to find the value of American options from Bermudan M options [37–39]. 2.2. Nonlinear Black–Scholes PDEs Relaxing the assumptions that there are no transaction costs and that the market is liquid, the linear Black–Scholes equation (1) is substituted by nonlinear pricing equations where the nonlinearity occurs in the coefficient of the diffusion term. The general form of the pricing PDE can be formulated as

( ) 2 1 2 ∂ V ∂ 2V ∂V ∂V 2∂ V + ˆ σ S, , S + rS − rV = 0, 2 2 ∂t 2 ∂S ∂S ∂S ∂S

(5)

where the function ˆ σ 2 is different for various popular transaction costs and illiquidity models. Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

4

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

Under the Leland’s transaction cost model,

ˆ σ =σ 2

2

∂ 2V 1 + Le sign ∂ S2

(

(

))

,

where σ is the original volatility in the Black–Scholes model and the Leland’s number is

√

kˆ √ . π σ δt 2

Le =

Here, kˆ represents the round trip transaction cost per unit dollar of the transaction and δ t is the interval between successive revisions of the portfolio [10,40,41]. For the risk adjusted pricing methodology [42,43],

⎛

C 2M ∂ 2V S ˆ σ 2 = σ 2 ⎝1 + 3 2π ∂ S 2

(

) 13

⎞ ⎠,

where M ≥ 0 denotes the transaction cost measure and C ≥ 0 is the risk premium measure. Barles and Soner [44] proposed a more complex transaction cost model where

( ( )) ∂ 2V ˆ σ 2 = σ 2 1 + ψ exp[r(T − t)]a2 S 2 2 . ∂S Here, the parameter a depends on the risk aversion factor and the number of options to be sold, and the function ψ (·) admits an implicit solution given in [14]. We also consider Frey’s illiquidity model [45] for which

σ2

ˆ σ2 = (

1 − α S ∂∂ SV2 2

)2 ,

and α , which usually takes small values [46], is the parameter modelling the liquidity of the market. Therefore, one may use the Taylor’s series approximation about α = 0 to have [47]

( ) ∂ 2V ˆ σ 2 = σ 2 1 + 2α S 2 . ∂S 2.3. The Merton’s jump–diffusion model This model [48] caters for the infrequent jumps that have been observed in stock prices. The price dynamics is dSt St

= (r − λJ ζ )dt + σ dWt + (Y − 1)dNt ,

where Nt is a Poisson process with arrival intensity λJ and Y − 1 is an impulse function making S jump to SY . Both Nt and Y are known to be uncorrelated to the Wiener process Wt . Under this model, Y is taken from a log-normal distribution with mean µJ and variance σJ2 . Therefore, the expectation of Y − 1, denoted by ζ , is given by

[ ζ = exp µJ +

σJ2

] − 1.

2

Considering the transformation x = ln(S /E), the computation of the value of a contingent claim V (t ; x) satisfies the partial integro-differential equation (PIDE)

( ) ∫ ∂V 1 2 ∂ 2V 1 2 ∂V + σ + r − σ − λJ ζ − (r + λJ )V + λJ V (t ; x + y)f (y)dy = 0, ∂t 2 ∂ x2 2 ∂x ℜ

(6)

where 1

[

f (y) = √ exp − 2π σJ

(y − µJ )2 2σJ2

] ,

and the boundary conditions in (3) are converted to

∂ 2V ∂V ≈ ∂ x2 ∂x

as x → −∞ and x → +∞.

(7)

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

5

2.4. The Heston’s stochastic volatility model Since the Black–Scholes assumption of constant volatility has been empirically shown to be flawed, stochastic volatility models have been proposed and one of the most popular models in literature is the Heston’s stochastic volatility model [15] which is specified as follows dSt St

= rdt +

√ vt dWt ,

with dvt = κ (θ − vt )dt + σv

√ vt dWtv ,

where the variance of the asset price is denoted by vt , the long term mean is θ , the mean reversion rate is κ , the volatility of the variance is σv and the correlation coefficient between the two Wiener processes Wt and Wtv is ρ . Under such a two-dimensional model, it can be shown that a contingent claim V (t ; S , v ) satisfies the following Garman’s type PDE

∂V ∂ 2V 1 ∂ 2V 1 ∂ 2V ∂V ∂V + v S 2 2 + σv2 v 2 + ρσv v S + rS + κ (θ − v ) − rV = 0. ∂t 2 ∂S 2 ∂v ∂ S ∂v ∂S ∂v The payoffs for exotic and path dependent options are similar as discussed earlier except that they now hold for all values of vt . The determination of appropriate boundary conditions along v under the Heston’s model remains an open problem. However, Dirichlet boundary conditions have been discussed in [49] where the authors also determine that if κθ −σv2 /2 ≥ 0, there is no need to implement boundary condition at the degenerate boundary v → 0. Following [5,20,50], one can also opt for Neumann boundary conditions which read as

∂V ≈ 0 as v → 0 and v → +∞. ∂v

(8)

We stress here that the Dirichlet boundary conditions obtained in [49] for different types of options as v → +∞ are derived from the Neumann boundary condition (8) which is consistent with market observations [50]. However, as we shall later see in our numerical experiments, care must be taken before imposing the far-field boundary condition (8) on a truncated domain. 3. Finite difference discretisation In this section, we outline the derivations of the non-uniform high-order approximations to the first-order and secondorder partial derivatives. Also, we describe the exponential time integration scheme used, and finally we demonstrate numerically that the proposed scheme is unconditionally stable. 3.1. Local mesh refinement strategy We start by constructing a uniform grid for our spatial domain, say the x-domain given by

{

Ω1 x = xi ∈ ℜ : xi = xmin + i1x, i = 0, . . . , n, 1x =

xmax − xmin n

}

,

and using the Taylor’s series expansion, the fourth-order approximations to the partial derivatives on Ω1x work out to be

∂V Vi−2 − 8Vi−1 + 8Vi+1 − Vi+2 1x4 (5) ≈ + V (xi ) + O(higher order terms), ∂x 121x 30

(9)

∂ 2V −Vi−2 + 16Vi−1 − 30Vi + 16Vi+1 − Vi+2 1x4 (6) ≈ + V (xi ) + O(higher order terms), 2 2 ∂x 121x 90 but as reported in [7], only a second-order convergence is observed due to the points of singularity that exist in most payoff functions. The remedial approach is to concentrate grid nodes at these points of singularity so as to minimise the magnitude of the error at these locations and hence, restore the high-order rate of convergence [2,21,22]. Currently, one of the best strategies was proposed in [7] whereby four nodes are added at each refinement stage. In our work, we propose to add only two nodes at each refinement stage. The strategy is illustrated in Fig. 1 as we double the value of n and increase the number of refinement stages R. We have assumed that there is a singularity at the node xi = 100 such that at each refinement stage, we add two additional nodes at (xi + xi±1 )/2. Thus, with R refinement stages, the total number of nodes now becomes N = n + 2R + 1. Here, we point out that for a discrete barrier option, the same strategy is also applied at the barrier level H but for continuously monitored barrier options, only nodes in the continuation region are added such that only one additional node is added near H at each refinement stage. Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

6

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

Fig. 1. Discretised domain after two steps of the local mesh refinement strategy as we double the value of n and increase the refinement stages R. We assume a point of singularity is located at 100.

3.2. Non-uniform high-order spatial discretisations Once the grid has been refined, it is no longer uniform and we cannot use the high-order spatial discretisations (9) defined on a uniform grid. Therefore, we next outline the derivations of the high-order approximations on non-uniform grids due to Bowen and Smith [23]. Let Dd,˜n [V ] denote the n˜ -point stencil approximation to the dth derivative of V at the point xi , where d < n˜ . Consequently, the following expression holds

⏐ n˜ ∑ ∂ d V ⏐⏐ ≈ D [ V ] = wk V (xk ), d ,˜ n ⏐ ∂ xd x=xi k=1

and in order to find the weights wk , k = 1 to n˜ , in the above linear combination, we use the first n˜ terms in the Taylor’s series expansion of V (xk ) about xi to obtain Dd,˜n [V ] =

n˜ ∑ k=1

⎡

( )⎤ n˜ −1 ∑ αkj ⎣wk Dj,˜n [V ] ⎦ , j! j=0

where αk = xk − xi . Then, by comparing the coefficients, the weights of the dth derivative of V satisfy the following set of linear equations 1

⎛

1

⎜ α1 ⎜ ⎜ .. ⎜ . ⎜ ⎜ ⎝ α1n˜ −1 ( ) n˜ − 1 !

α2 .. . α2n˜ −1 ) n˜ − 1 !

(

··· ··· .. . ···

1

⎞

⎛ ⎞ ⎛ ⎞ δ0,d ⎟ w1 αn˜ ⎟ w ⎟ ⎜ 2 ⎟ ⎜ δ1,d ⎟ .. ⎟⎜ . ⎟ = ⎜ . ⎟, . ⎟⎝ . ⎠ ⎝ . ⎠ . . ˜n−1 ⎟ ⎠ w αn˜ δn˜ −1,d n˜ ( ) n˜ − 1 !

where

δk,d

1, 0,

{

if k = d; otherwise.

Since we are interested in the case n˜ = 5, the weights can be easily solved by using the Mathematica code shown in Fig. 2. The approximations in Fig. 2 are used for i = 2 to N −3, where N is the total number of nodes after R refinement stages have been performed. Changing the reference point xi to x0 , x1 , xN −2 and xN −1 , the one-sided and the noncentral approximations Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

7

Fig. 2. Mathematica code to obtain the weights for D1,5 and D2,5 . wk is given by the kth row in each column.

can be calculated by simply using

{0, x1 − x0 , x2 − x0 , x3 − x0 , x4 − x0 }, {x0 − x1 , 0, x2 − x1 , x3 − x1 , x4 − x1 }, {αi−2 , αi−1 , αi , αi+1 , αi+2 } = { xN −5 − xN −2 , xN −4 − xN −2 , xN −3 − xN −2 , 0, xN −1 − xN −2 }, ⎪ ⎩ {xN −5 − xN −1 , xN −4 − xN −1 , xN −3 − xN −1 , xN −2 − xN −1 , 0},

⎧ ⎪ ⎨

at x0 ; at x1 ; at xN −2 ; at xN −1 ,

in Fig. 2. Also, the corresponding errors derived in [23] are expressed as

∂V e4 (5) e1 e4 − e5 (6) (xi ) = − V (xi ) − V (xi ) + O(higher order terms), (10) ∂x 120 720 2 ∂ V e3 (5) e1 e3 − e4 (6) D2,5 [V ] − (xi ) = V (xi ) + V (xi ) + O(higher order terms), (11) ∂ x2 60 360 where ej , j = 1 to 5, are formed by adding together all distinct products of j distinct objects from the set {αi−2 , αi−1 , αi , αi+1 , αi+2 }. For instance, D1,5 [V ] −

e4 = αi−2 αi−1 αi αi+1 + αi−2 αi−1 αi αi+2 + αi−2 αi−1 αi+1 αi+2 + αi−2 αi αi+1 αi+2 + αi−1 αi αi+1 αi+2 . Clearly, on a uniform grid, ej = 0 for j odd and with e4 = 41x4 , we recover the fourth-order discretisations (9) for uniform grids. Hence, coupling the above scheme with the local mesh refinement strategy, we expect a high-order accuracy in space. An explanation for the high-order rate of convergence based on the non-uniform grid in Fig. 1 is provided in the Appendix. 3.3. Integro discretisation under the Merton’s jump–diffusion model For the discretisation of the integral term in (6), we let z = x + y such that

∫

∫

Ix (t) =

V (t ; z)f (z − x)dz .

V (t ; x + y)f (y)dy = ℜ

ℜ

The approximation to the above integral can be obtained by using the high-order Simpson’s composite rule combined with a fast Fourier transform (FFT) technique for efficiently computing the dense matrix–vector multiplication in O(N log N) complexity [7]. Since this integral term is non-local, the solution depends on an additional boundary correction term that is dependent on the option being priced. For more details about the implementation, the reader is referred to [7,28]. Another approach to work out the integral has been proposed in [51] where the authors have found that Ix (t) is the exact solution to the problem

σJ2 ∂ 2U 2µJ ∂ U ∂U = + , −∞ ≤ x ≤ ∞, 0 ≤ τ ≤ , 2 ∂τ ∂ x2 2 σJ ∂ x

(12)

at τ = σJ2 /2 subject to the initial condition U(0; x) = V (t ; x) and the same boundary conditions (7). After discretising the spatial derivatives in the above PDE, we come up with an ordinary differential equation (ODE) of the form dU dτ

= Dx2 U +

2µJ

σJ2

Dx1 U ,

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

8

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

where Dx1 U and Dx2 U are the non-uniform high-order approximations to the first-order and second-order derivatives of U with respect to x respectively, and U is a vector which approximates U on the spatial computational domain. Using separation of variables, the solution to the above ODE can be easily solved using a one-step exponential time integration to give

[ Ix (t) = exp

σJ2 2

] Dx2

+µ

x J D1

V (t ; x).

(13)

This approach is versatile and does not depend on the options being priced. Therefore, there is no need to use the asymptotic behaviour of the option values at far-field in order to compute the correction terms needed when the PIDE (6) is localised [7,28]. However, computing the matrix exponential in (13) is a computational task with O(N 3 ) complexity if the expm Matlab in-built function is used [52,53]. We shall next describe how this task can be achieved in O(N) complexity which is also an improvement over the O(N log N) complexity. 3.4. Exponential time integration scheme with Carathéodory–Fejér approximations (ETI-CF) Once the spatial derivatives in the pricing PDEs are approximated using high-order discretisations, we remain to approximate the temporal derivative which also affects the general rate of convergence and the computational efficiency of a numerical method. Here, we apply a time reversed transformation, τ = T − t, to transform our problem into the initial boundary value problem dV

= LV + λJ Ix (τ ), (14) dτ where Ix (τ ) is given in (13) and λJ = 0 for the Black–Scholes and Heston’s models. L represents the discretised spatial operator given by 1

σ 2 Is2 Ds2 + rIs Ds1 − rDs0 , 2 under the Black–Scholes model, by L=

L=

1 2

(15)

( ) 1 σ 2 Dx2 + r − σ 2 − λJ ζ Dx1 − (r + λJ )Dx0 , 2

for the Merton’s jump–diffusion model and by 1 Iv ⊗ (Is2 Ds2 ) + σv2 (Iv Dv2 ) ⊗ Ds0 + ρσv (Iv Dv1 ) ⊗ (Is Ds1 ) + rDv0 ⊗ (Is Ds1 ) + [Iκ (θ −v ) Dv1 ] ⊗ Ds0 − rDv0 ⊗ Ds0 , 2 2 under the Heston’s stochastic volatility model, where ⊗ is the Kronecker product. Here, Is represents a diagonal matrix containing the grid nodes in S, and in our case, Ds0 is an identity matrix on S. We have similar easily understandable notations in the x-domain and v -domain. Solution to ODE (14) can be obtained by separation of variables. Under the Black–Scholes and Heston’s models, we yield L=

1

V (τ + 1τ ) = exp[A]V (τ ),

(16)

where matrix A = L1τ for a time step 1τ = T /m with m being the number of time steps. For European style options under linear models, only one time step (m = 1) is required since the exponential time integration (ETI) advances the unconstrained solution exactly in time such that no temporal errors exist. However, for American or nonlinear problems with time dependent coefficients, there will exist a temporal error which we shall efficiently eliminate using a Richardson extrapolation technique. Following the work of [54], the matrix exponential and vector product can be efficiently approximated by exp[A]V (τ ) ≈

η ∑ (

cj A − zj D0

)−1

V (τ ) ,

(17)

j=1

where the residues cj and poles zj can be computed based on the singular value decomposition of a Hankel matrix [54]. The approximation error is as small as O(9.28903−η ), and since the option values are real, the residues and poles come in complex conjugate pairs [3] such that only η/2 linear solves are necessary to find (16). However, in the case of the Merton’s jump–diffusion model, we note that the integral term Ix (τ ) = exp[B]V (τ ; x), where the matrix B = (σJ2 /2)Dx2 + µJ Dx1 , destroys the sparsity structure of the right hand side of (14). Consequently, we propose to manipulate the code in [54] to evaluate V (τ + 1τ ) = exp[A + λJ 1τ exp(B)]V (τ ). Since A and B are made up of the same matrices Dx1 Dx2 ≈

∂ ∂x

(

∂ ∂ x2 2

∂ ∂ = 2 ∂ x3 ∂x 3

) =

2

(

∂ ∂x

)

(18) Dx1

and

Dx2 ,

and given that

≈ Dx2 Dx1 ,

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

9

(18) can be rewritten as V (τ + 1τ ) = exp[λJ 1τ exp(B)] exp(A)V (τ ) .

[

]

V 1 = exp(A)V (τ ) can be obtained from (17) where cj and zj are computed using the code produced in [54]. Then, the code line, F = exp(scl ∗ (t − 1)./(t + 1 + 1e − 16)), in [54] is then modified as F = exp[λJ 1τ (exp(scl ∗ (t − 1)./(t + 1 + 1e −16)))]−1 to obtain V 2 = exp[λJ 1τ exp(B)]V 1 − V 1 , so that V (τ + 1τ ) = V 2 + V 1 . The same formula in (17) is used where B and V 1 are used instead of A and V (τ ). As such, we are able to price options under the Merton’s jump–diffusion model in O(N) complexity since only η sparse linear solves are performed over one time step. 3.5. An iterative method to solve nonlinear Black–Scholes equations In order to solve our nonlinear pricing problems, we first discretise the nonlinear PDE (5) in space using the non-uniform five-point stencil approximations. This results in an ODE of the form dV dτ

= LV + LN (V )V ,

(19)

where L is the discretised differential operator under the Black–Scholes model, and LN (V ) represents the nonlinear discretised operator given as

⎧ ( s ) s 1 2 ⎪ ⎪ ⎪ σ Is2 Le sign D2 V D2 , ⎪ 2 ⎪ ⎪ ( 2 ) 13 ⎪ ⎪ 3 ⎪ ⎨ σ 2 I 2 C M I Ds V Ds , s 2 s 2 2π LN (V ) = 2 ⎪ )] [ ( 1 ⎪ 2 ⎪ σ Is2 ψ exp[r τ ]a2 Is2 Ds2 V Ds2 , ⎪ ⎪ ⎪2 ⎪ ⎪ ) ⎪1 2 ( ⎩ σ Is2 2Is α Ds2 V Ds2 , 2

for the Leland’s model; for the risk adjusted pricing methodology; for the Barles and Soner’s model; for the Frey’s model.

If LN (V ) was known and was independent of τ , then the value of a European option at τ + 1τ could be calculated from V τ +1τ = eL1τ +LN (V )1τ V τ , where V τ denotes the price of the European call option at time τ . However, since LN (V ) is unknown, we propose to use an iterative scheme to solve (19). Such iterative techniques have been employed to solve the nonlinear PDEs arising from the 0 application of the penalty method for pricing American options [32,55]. Assuming that V τ +1τ = V τ , the algorithm of our iterative scheme over each time step can be described as follows for a maximum number of iterations denoted by Mitr. while error>tolerance and q ≤ Mitr do q

q+1

V τ +∆τ = eL∆τ +LN (V τ +∆τ )∆τV τ q+1

q

error = V τ +∆τ − V τ +∆τ 





∞

q = q+1 end while q+1

In order to gain in computational efficiency, we use the Carathéodory–Fejér approximations to evaluate V τ +1τ . Here, we point out that the temporal discretisation is very important to yield correct values since this approach can be verified to give only a first-order accuracy with respect to time. Therefore, to recover a high-order convergence, we employ the Richardson extrapolation technique [33] for which the tableau when doubling m at each stage k is given as Tk,1 = Vmk , Tk,j =

k = 1 , . . . , s,

2j−1 Tk,j−1 − Tk−1,j−1 2j−1 − 1

(20)

,

for k = 2, . . . , s

and

j = 2, . . . , k.

3.6. Stability analysis We briefly conduct a stability analysis when our proposed approach is used to price options under the various models discussed before. Consider a semi-discretised system of the form dV dτ

= LV (τ ),

(21)

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

10

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

Table 1 Table of parameters.

σ = 0.2, r = 0.05, T = 1, E = 100. σ = 0.2, r = 0.05, T = 1, E = 100, λJ = 1, µJ = 0, σJ = 0.2. r = 0.05, T = 0.25, E = 100, κ = 2.0, ρ = −0.5, θ = 0.1, σv = 0.1. kˆ = 0.01, δ t = 0.01, C = 2, M = 0.01, a = 0.02, α = 0.05.

Black–Scholes model Merton’s model Heston’s model Nonlinear Black–Scholes models

(a) Merton’s model.

(b) Heston’s model.

Fig. 3. Logarithmic norm µ[L] against N for linear models.

where here, for simplicity of illustration, L also involves the integral operator of the Merton’s model. For (21) to be stable, we need to show whether there exists rigorous estimates on the exponential of the matrix τ L of the type

∥ exp[τ L]∥ ≤ K exp[τ ω].

(22)

If ω and K are of moderate size, then (22) implies that errors will propagate only in a mild fashion. If ω ≤ 0 and K = 1, then the semi-discretised system satisfies a strong stability property (contractivity) [56]. In this paper, we consider the logarithmic norm in the Euclidean space which is given by

{ } ) 1( µ[L] = max λ : λ is the eigenvalue of L + LT , 2

to find upper bounds for contractivity. The following theorem from [56] allows us to relate the exponential of the matrix τ L to the logarithmic norm. Theorem. Let L = (Lij ) ∈ ℜN ×N and ω ∈ ℜ. Then

µ[L] ≤ ω ⇔ ∥ exp[τ L]∥ ≤ exp[τ ω],

for all τ ≥ 0.

So, the upper bounds on µ[L] should be negative (ω ≤ 0) to prove strong stability. Exhaustive derivations are involved for non-uniform finite difference approximations and we refer interested readers to [56,57]. Consequently, it is easier to show numerically that µ[L] ≤ 0. For the parameters given in Table 1, Fig. 3 shows that the logarithmic norms of the operators in space are less than zero and decrease monotonically as N] increases. For nonlinear Black–Scholes models, we follow the [ same procedures to demonstrate in Fig. 4 that µ L + LN (V ) ≤ 0, where L(V ) is calculated using the values of V from the last iteration in the last time step. Therefore, we can conclude that the proposed scheme is unconditionally stable. 3.7. Construction of the non-uniform computational grid Truncating the semi-infinite S-domain [0, +∞] under the Black–Scholes and Heston’s models, and the infinite x-domain [−∞, +∞] under the Merton’s model is the first task to be performed. If the domain is truncated at large values, the computational efficiency is lowered as more nodes are required to achieve satisfactory accuracy while truncating the domain over a small range is known to affect the global rate of convergence. We use the truncation range

[ ] √ √ √ √ [alx , blx ] = (c1 + x0 ) − l c2 + c4 , (c1 + x0 ) + l c2 + c4 , Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

(a) Leland’s model.

(b) Risk adjusted pricing methodology.

(c) Barles and Soner’s model.

(d) Frey’s model.

11

Fig. 4. Logarithmic norm µ L + LN (V ) against N for nonlinear Black–Scholes models.

[

]

as suggested in [58] with l depending on a user-defined tolerance level and cn denoting the nth cumulant of xt = ln(St /E). The cumulants can be obtained using cn =

1 ∂ n (t φ (u)) ⏐ i

∂ un

⏐ ⏐ ⏐

,

u=0

where t φ (u) is the exponent of the characteristic function and i is the unit imaginary number. We refer to [58] for the detailed computations of these cumulants. Since we are working with the non-uniform approximations and based on our intensive numerical experiments, we propose to use a uniform grid 1x over the range [a3x , b3x ] and an appropriate coarser uniform grid, say 51x, for the remaining grid nodes on the interval [a7x , b7x ]. Similar partially uniform grid can be also constructed using 1S and 51S along the S-domain. An analogous approach is adopted while truncating the v -domain [0, +∞] under the Heston’s model. Since values of v > 0.5 are not economically significant [5], limiting the domain at the small value vmax = 0.5 is computationally favourable. However, applying the boundary condition (8) as v → +∞ at vmax = 0.5 which is not large enough can result in relatively larger errors, especially at vmax itself. Consequently, we propose to use a uniform grid with 1v over the range [0, vmax ] and extend the domain to sufficiently large vmaxext with a much larger mesh size. We shall assess the efficiency of this strategy in the next section. 4. Numerical experiments In this section, we illustrate the performance of the proposed approach which consists of using the non-uniform highorder discretisations in space with the local mesh refinement at the singularities coupled with the ETI-CF scheme with η = 12. We use the parameters given in Table 1 unless stated otherwise. All our experiments are done on an Intel(R) Core(TM)i5-3230M CPU @2.60 GHz processor with 4GB RAM using Matlab R13. We first price European call options under the Black–Scholes model where we compare the results of our approach with those from two methods. The first comparison method consists of the central second-order discretisations in space combined with the ETI-CF temporal scheme (ETI-CF2). The second one is composed of the spatial HOC scheme in [6] with a local mesh Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

12

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

Fig. 5. Log–log plot of error at S0 = 100 against N when pricing European call options under the Black–Scholes model. Table 2 Infinity norm of the errors and convergence rates when pricing discrete and continuous barrier options under the Black–Scholes model using our nonuniform high-order ETI-CF scheme. R

1 2 3 4 5

1S 2 1 0.5 0.25 0.125

Discrete barrier

Continuous barrier

N

Error

Order

Time (s)

N

Error

Order

Time (s)

144 288 572 1136 2260

4.4960e−3 2.6875e−4 1.6433e−5 1.0460e−6 6.4732e−8

– 4.0643 4.0316 3.9737 4.0142

0.0935 0.0954 0.0999 0.1088 0.1382

63 126 249 492 975

1.4482e−3 9.0777e−5 5.6732e−6 3.5472e−7 2.2187e−8

– 3.9958 4.0001 3.9994 3.9989

0.0087 0.0090 0.0093 0.0097 0.0100

refinement strategy [7] coupled with an implicit Richardson extrapolation time-stepping scheme (HOC IM-RE). In this case, we use the first-order implicit scheme with 10 initial time steps for 4 Richardson extrapolation stages where we double the number of time steps at each stage to eliminate the temporal errors. The slopes of the lines in Fig. 5 show that a high-order convergence is attained only for our method and the HOC IM-RE scheme while the ETI-CF2 scheme converges more slowly as expected. However, we can also observe that the HOC IM-RE approach which is based on a uniform grid uses more points than our method to reach the same level of accuracy. This demonstrates a gain in computational efficiency by the proposed scheme on which an adaptive grid can be formulated. Next, we price continuous and discrete up-and-out barrier options for 12 monitoring dates under the Black–Scholes model where the barrier level is set at H = 120. Here, we apply our refinement strategy at the strike E and the barrier level H. In Table 2, we report the infinity norm of the errors computed over the whole of our computational domain for both discrete and continuous barriers. Reference prices for continuous barriers are computed in [35] and for the discretely monitored barriers, we use the very accurate COS method [39]. It can be observed that a high-order rate of convergence is achieved with accuracy up to 8 decimal places reached in less than one second. The extra time required for discrete barrier options is due to the fact that the 6 linear systems need to be solved over each monitoring date. We also calculate the option price sensitivities commonly referred to as the Greeks. For the Delta (1) and Gamma (Γ ) which are approximated as

∆=

∂V ≈ Ds1 V , ∂S

and Γ =

∂ 2V ≈ Ds2 V , ∂ S2

respectively, we consider the digital and butterfly options since the singular points that exist in their payoff functions have notoriously been known to produce spurious oscillations for schemes which are not L-stable [59]. For B = 100, E1 = 90, E2 = 100, and E3 = 110, we record the infinity norm of the errors and convergence rates for the 1 and Γ of these options in Table 3. The results show that the kinks – at E for the digital option, and at E1 , E2 and E3 for the butterfly option – do not reduce the high-order rate of convergence for the Greek parameters. We next solve the nonlinear Black–Scholes PDEs using our approach combined with the iterative method for a tolerance of 1e-10 and Mitr = 5. We compare the results of our approach with different existing techniques when pricing a European call option. The first comparative method (Ankudinova & Ehrhardt method) [10] consists of treating the coefficient of the secondorder derivative in the nonlinear term explicitly and discretising the latter with a central second-order approximation with a larger mesh size 21S. Then, linear terms are discretised using central second-order approximations on a mesh size 1S, Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

13

Table 3 Infinity norm of errors and convergence rates of 1 and Γ for digital and butterfly options under the Black–Scholes model using our non-uniform high-order ETI-CF scheme. Digital option N

50 100 198 392

Butterfly option

1

Γ

1

N

Error

Order

Error

Order

2.5918e−3 1.6637e−4 1.0372e−5 6.5556e−7

– 3.9614 4.0037 3.9838

1.6667e−4 1.0830e−5 6.6343e−7 4.1716e−8

– 3.9439 4.0289 3.9913

96 198 396 786

Γ

Error

Order

Error

Order

2.1929e−4 1.4172e−5 8.8827e−7 5.5552e−8

– 3.9517 3.9959 3.9991

2.3061e−5 1.4981e−6 9.4007e−8 5.8791e−9

– 3.9443 3.9942 3.9991

Table 4 Convergence of the price of a European call option at S0 = 100 when solving the nonlinear problems under the Leland’s model and the risk adjusted pricing methodology (RAPM). For the Ankudinova & Ehrhardt method, we use m = 2n and m = 4n respectively. For RAPM, our method requires an initial number of 10 time steps for s = 5. Leland’s model Our method

Ankudinova & Ehrhardt method

N

Error

Order

Time (s)

n

Error

Order

Time (s)

52 102 200 394

1.8603e−3 1.1109e−4 6.7751e−6 3.9501e−7

– 4.0657 4.0354 4.1003

0.0257 0.0263 0.0384 0.0423

80 160 320 640

5.7449e−2 1.4307e−2 3.6288e−3 9.8260e−4

– 2.0056 1.9791 1.8848

0.0453 0.1008 0.1719 0.4957

5.8182e−2 1.4132e−2 3.1646e−3 4.9750e−4

– 2.0416 2.1589 2.6692

0.0929 0.1628 0.4250 1.3463

Reference N = 780

11.827536

Risk adjusted pricing methodology 52 102 200 394

2.2901e−3 2.8711e−4 1.5876e−5 5.1627e−7

Reference N = 780

– 2.9957 4.1767 4.9426

0.6263 0.8236 1.1783 1.7991

80 160 320 640 12.725266

and the Crank–Nicolson scheme is used to advance the solution in time. We also compare our results with the locally one-dimensional methods namely the LOD backward Euler method in [11] and the LOD Crank–Nicolson method in [12]. The method in [10] is chosen for comparison because of its accuracy and those in [11,12] are chosen for their stability. Using the parameters in Table 1, we tabulate our results for the Leland’s model and the( risk adjusted ) pricing methodology in Table 4. For the Leland’s model where the nonlinear term involves the simple term sign ∂ 2 V /∂ S 2 , our intensive numerical experiments reveal that approximations over one time step is sufficient to give high-order convergent results. We analyse the convergence of the error using reference prices from our approach with N = 780 for both the Leland’s model and the risk adjusted pricing methodology since we believe that our scheme gives better accuracy. Clearly, our pentadiagonal iterative scheme exhibits a high-order convergence compared to the Ankudinova & Ehrhardt method which demonstrates only a second-order convergence. Moreover, we have observed that the latter requires stability restrictions on the time step because of the explicit treatment of the nonlinear diffusion terms. For instance, if we were to use the same number of time steps and spatial steps (m = n) for the risk adjusted pricing methodology, the value would drop down to 12.7236 when n = 2560 after being 12.7250 for n = 1280. Indeed, approximating the nonlinear terms with larger mesh sizes improves the stability of the scheme but does not solve the instability problem completely. For the Barles and Soner’s model, we compare our values with those from the Ankudinova & Ehrhardt method and the LOD method in [11]. For the Frey’s model, we set r = 0 to be able to compare our results with the LOD method in [12]. Table 5 illustrates that our approach converges faster for both models. The Barles and Soner’s model takes relatively more time to run compared to the other models because of the need to interpolate values for the function ψ (·) at each time step. For Frey’s model, we can indeed observe that the LOD Crank–Nicolson method is a better method compared to the LOD backward Euler method. However, the results show that our approach is still the most efficient method. We notice that a large number of time steps is needed for the LOD methods because the splitting method used induces a splitting error which depends on the size of the time step. We have solved the linear Black–Scholes PDE with the LOD methods, and if m is chosen too small, convergence to another solution other than the exact solution was observed. Consequently, our approach proves to be computationally efficient and stable for solving nonlinear PDEs. We also use our approach to price European options under the Merton’s jump–diffusion model. In Table 6, we clearly see that without a local mesh refinement at the strike E, only an error of about e−4 is obtained compared to an error of e−7 when using 5 refinement stages, that is, with only 10 additional grid nodes. This strategy is clearly shown to restore the high-order rate of convergence. Moreover, we compare the results obtained when the integral term of the PIDE (6) is treated differently. Using our non-uniform high-order approach in space, we first evaluate the matrix exponential in (13) using the one-step ETI-CF approach in Section 3.4 (ETI-CF). Secondly, we implement the high-order composite Simpson’s rule-FFT approach in [7] Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

14

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

Table 5 Convergence of the price of a European call option at S0 = 100 when solving the nonlinear Black–Scholes PDEs resulting from the Barles and Soner’s model and the Frey’s model. For the Ankudinova & Ehrhardt method, we use m = 8n while for our method we extrapolate with s = 5 and s = 6 respectively. Barles and Soner’s model Our method

Ankudinova & Ehrhardt method

LOD backward Euler method

N

Error

Time (s)

n

Error

Time (s)

n

m

Error

Time (s)

26 52 102 200

1.7250e−2 2.8294e−3 2.9710e−4 1.1064e−5

1.8456 1.9825 2.2494 2.6143

40 80 160 320

2.4416e−1 6.1639e−2 1.5022e−2 3.1022e−3

0.4985 0.5701 1.2194 2.5247

32 64 128 256

1600 6400 25600 102400

2.5498e−1 5.6024e−2 1.7811e−2 8.3881e−3

2.3280 8.2473 32.778 153.58

Reference N = 394

12.662399

Frey’s model Our method

LOD backward Euler method

LOD Crank–Nicolson method

N

Error

Time (s)

n

m

Error

Time (s)

n

m

Error

Time (s)

208 402 788

1.5336e−4 1.3552e−5 3.6852e−7

2.1267 3.1966 5.7022

128 256 512

25600 102400 409600

1.9791e−2 7.5643e−3 4.6700e−3

15.1872 68.2701 336.719

128 256 512

1600 6400 25600

1.4897e−2 2.8749e−3 1.0025e−4

1.3065 2.3933 5.4616

Reference N = 1558

8.891599

Table 6 Infinity norm of errors and convergence rates when pricing European call options under Merton’s jump–diffusion model.

1x 0.1 0.05 0.025 0.0125 0.00625

Our method with no LMR

Our method with LMR

n

Error

Order

R

N

Error

Order

30 60 122 244 488

1.2758e−1 3.0867e−2 7.6551e−3 1.9114e−3 4.7762e−4

– 2.0473 2.0116 2.0018 2.0007

1 2 3 4 5

32 64 128 252 498

3.1543e−2 2.1939e−3 1.4499e−4 9.2338e−6 5.8682e−7

– 3.8457 3.9195 3.9729 3.9759

(a) Log–log plot of error against N.

(b) Log–log plot of error against time.

Fig. 6. Error graphs when pricing European call options under the Merton’s jump–diffusion model at S0 = 100.

and use a first-order implicit–explicit Richardson extrapolation time-stepping scheme with initially 40 time steps for 4 extrapolation stages as we double the number of time steps (SIMP-FFT-IMEX-RE). Here, the composite Simpson’s rule-FFT approach in [7] is used for comparison since it is the most efficient method for treating the integral term under the Merton’s model in O(N log N) complexity. To the best of our knowledge, no numerical PDE based algorithms have been able to evaluate such non-local integral in less than O(N log N) since the complexity to solve (12) for each time step in [51] is O(mN). We note, from Fig. 6a, that a high-order convergence is reached for both approaches with the ETI-CF scheme using less number of nodes than the SIMP-FFT-IMEX-RE one. This is because the SIMP-FFT-IMEX-RE scheme needs a uniform grid over the local domain to efficiently evaluate the integral component in (6). Moreover, the error plot against the computational time in Fig. 6b indicates clearly that the ETI-CF method is far more efficient having O(N) complexity compared to the O(N log N) complexity of the SIMP-FFT-IMEX-RE approach. Indeed, one can observe that the time taken by the ETI-CF is less than that of the SIMP-FFT-IMEX-RE approach to reach a similar level of accuracy. Next, we price American put options under the Black–Scholes and Merton’s jump–diffusion models by extrapolating Bermudan option prices. The parameters T = 0.5 and λJ = 0.1 are used for these numerical experiments. In order to Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

15

Table 7 Relative errors and convergence rates when pricing American put options at S0 = 100 under the Black–Scholes and Merton’s jump–diffusion models with reference prices given by 4.65568477 and 4.84116836 respectively using the COS method in [39]. Extrapolating Bermudan options with an initial of 16 monitoring dates for 6 extrapolation stages R

0 1 2 3 4

Black–Scholes model

Merton’s jump–diffusion model

1S

N

Error

Order

Time (s)

1x

N

Error

Order

Time (s)

5 2.5 1.25 0.625 0.3125

63 85 138 236 479

5.4856e−2 2.7150e−3 1.7563e−4 1.2544e−5 9.5516e−7

– 4.3366 3.9503 3.8075 3.7151

0.8000 0.9063 1.1226 1.5853 2.5251

0.1 0.05 0.025 0.0125 0.00625

12 30 73 154 303

2.4816e−1 1.1851e−2 7.7410e−4 6.4345e−5 2.9717e−6

– 4.3882 3.9364 3.5886 4.4365

1.0091 1.1155 1.3504 1.9399 3.1174

Fig. 7. Free boundary location for the American put option under the Black–Scholes model.

price American options, we first run our previously used code for European options and enforce the payoff constraint (4) at ˆ to locate the optimal exercise boundaries Sf (tmˆ ) under Bermudan early exercise ˆ = 1, 2, . . . , M, each monitoring date tmˆ , m rights. Then, the grid is constructed such that local refinement is performed at these optimal locations. For more efficiency, the lower boundary value of our computational grid can be taken relatively close to Sf (tMˆ ). Under the Black–Scholes model, Fig. 7 illustrates a plot of the optimal exercise boundary computed by our method, and it is seen to be in good agreement with the optimal free boundary locations of the Goodman & Ostrov method [60] and of the COS method [39]. The methods proposed in [39,60] are well-known to give accurate results for the free boundary location of American options. Thus, we can see that the value of Smin can be chosen as 80. Indeed, the free boundary function being a monotonic function, it suffices to know that under the Black–Scholes model [61] Smin ≈ Sf (∞) =

µ E, µ−1

where µ =

( ) √( )2 − r − 21 σ 2 − r − 12 σ 2 + 2σ 2 r σ2

.

ˆ k = 8(2k ), k = 1, . . . , s, the corresponding value of an American option can Having the Bermudan option values VMˆ for M k

ˆ s ). In be obtained from Ts,s using the Richardson extrapolation tableau (20) with a computational complexity of about O(N M Table 7, we can observe, for s = 6, that highly accurate option prices are obtained for American options due to the high-order rate of convergence, and the algorithm requires only a few seconds to execute. The reference prices are again taken from the highly accurate COS method in [39]. Finally, we conduct numerical experiments for pricing European call options under the two-dimensional Heston’s stochastic volatility model. First, we investigate the fall in the error at vmax = 0.5 using two approaches. The first one involves truncating the v -domain at the small but economically large enough value of vmax = 0.5 while the second one (our approach) consists of truncating the v -domain at a larger value of vmaxext . From our intensive numerical experiments, we have seen that vmaxext ≈ 2 is large enough so as not to have significant influences on the solution in the domain of interest [0, vmax ]. For computational efficiency, we again exploit the non-uniformity of our proposed scheme and use a coarse grid in the extended v -domain. In [5], the far-field boundary condition (8) was imposed at vmax = 0.5 but as we can see from Fig. 8, the errors are relatively large and stagnate when the v -domain is truncated at the too small value of vmax = 0.5. However, for our approach where only a few additional nodes have been used to implement the boundary condition (8) at vmaxext = 2, the errors are seen to drop with a high-order rate of convergence. Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

16

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

Fig. 8. Log–log plot of errors at S0 = 100 and vmax = 0.5 against N. vmax = 0.5 and vmaxext = 2. Refinement is carried out at S0 = 100 and v0 = 0.2, and the reference solution is taken from [58]. Table 8 Error and convergence rates when pricing European call options at S0 = 100 and v0 under the Heston’s model. The reference solution is 5.1993497562 for v0 = 0.04 and 8.9913266587 for v0 = 0.2 [58]. For the Ballestra method, we take n initial number of time steps for a two-stage Richardson extrapolation, and the number of steps in v -domain is n/4. Our method R

2 3 4 5 6

N

52 102 200 394 780

Ballestra & Cecere method

v0 = 0.04

v0 = 0.2

Time (s)

Error

Order

Error

Order

3.4936e−3 2.7304e−4 1.8562e−5 1.1705e−6 9.4680e−8

– 3.6775 3.8787 3.9872 3.6279

2.4121e−3 1.4781e−4 9.0757e−6 5.5731e−7 3.0454e−8

– 4.0285 4.0256 4.0255 4.1938

0.5770 1.3707 2.9078 6.2121 14.054

n

20 40 80 160 320

v0 = 0.04

v0 = 0.2

Time (s)

Error

Order

Error

Order

1.1678e−2 1.7392e−3 7.2250e−5 3.2205e−6 6.0193e−7

– 2.7473 4.5893 4.4876 2.2751

3.5784e−2 8.6869e−4 5.0841e−4 4.9572e−4 4.9231e−4

– 5.3644 0.7728 0.0365 0.0098

0.0915 0.3651 2.5774 24.068 378.35

We also compare our results with the recently proposed method in [5] (Ballestra & Cecere method) which consists of using grid transformations with Chebyshev approximations that are known to give spectral convergence. Further, the authors employed an operator splitting method in space to avoid working with excessively large dense matrices and implemented an implicit–explicit time-stepping scheme combined with a two-stage Richardson extrapolation. From Table 8, we can observe that our method is computationally more efficient while for the Ballestra & Cecere method, we can see that the error stagnates at about e−4 for v0 = 0.2. While the transformation used in [5] can easily concentrate grid nodes at the strike price E, it is more difficult to do so at v0 = 0.2 since such transformation is almost linear for small values of v . However, the Chebyshev points used in [5] naturally concentrate grid nodes at boundaries such that we can observe, from Table 8, a high accuracy when v0 = 0.04. Nonetheless, our method is still superior where better accuracy is obtained more rapidly for all economically significant values of v . We have used only 56 + 2R nodes in v with R refinement stages at v0 . 5. Conclusion In this paper, we have proposed the use of a non-uniform high-order spatial discretisation for option pricing. Such discretisation can easily adapt itself to the partially non-uniform grid that results from the refinement at the points of singularity to restore the high-order convergence rate and the coarser grid construction which accounts for the implementation of farfield boundary conditions. We combine this high-order spatial discretisation with an exponential time integration scheme efficiently implemented using Carathéodory–Fejér approximations. This gives a highly accurate technique which exhibits a high-order rate of convergence for pricing exotic and American options under the Black–Scholes and Merton’s jump– diffusion models. Moreover, using an iterative scheme with a Richardson extrapolation time-stepping scheme, we have also been able to solve several nonlinear Black–Scholes option pricing problems with a higher order accuracy. This new approach is also shown to efficiently price options under the two-dimensional model of Heston. Indeed, this scheme can be easily adapted to price many other financial instruments like interest rate derivatives [62,63] as well as to solve the computational challenges under modern financial models such as regime-switching models with jumps. We shall report on these findings in a future work. Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

17

Acknowledgements We wish to thank the anonymous referees for their insightful and constructive comments which have significantly improved the presentation of our work. The research work of M. J. Dilloo is supported by a postgraduate research scholarship from the University of Mauritius under the grant number Q0083. Appendix. Explaining the high-order convergence of our non-uniform spatial discretisation In this section, we describe how a high-order convergence is obtained with our technique using the non-uniform grid in Fig. 1. First, we consider a uniform grid with n + 1 nodes where we let the x-domain be discretised with 1x. At the point of singularity, it is well-known that we have a reduced order of convergence of O(1x2 ). As we halve 1x (double n) with an additional refinement stage, we can observe from Fig. 1 that the point of singularity at x = 100 lies on a uniform grid with 1x/4 and 1x/16 respectively. Hence, the ratio of convergence can be obtained as

( 1x )2 1x2 4 ( 1x )2 = ( 1x )2 = 16. 4

16

Here, we stress on the fact that the errors in option pricing are usually dominant around the points of singularity. Consequently, by applying a refinement strategy at those points, we can expect to yield a high-order rate of convergence when reporting the errors at the kinks or the errors in the infinity norm over all the grid points including the refinement points. Indeed, our discretised spatial operator consists of the matrices D1,5 [V ] and D2,5 [V ] which satisfy (10) and (11) respectively, with e3 = αi−2 αi−1 αi+1 + αi−2 αi−1 αi+2 + αi−2 αi+1 αi+2 + αi−1 αi+1 αi+2 , e4 = αi−2 αi−1 αi+1 αi+2 , where αk = xk − xi , and since xi is a grid node, αi = 0. Here, we argue that our local mesh refinement strategy gives rise to partially uniform grids, that is, the construction of the computational grid is not completely arbitrary. Therefore, in our case, αk is proportional to 1x which halves as we double n. Thus, D2,5 [V ] is only fourth-order on a uniform grid (e3 = 0) but third order on a non-uniform grid (e3 ̸ = 0) compared to D1,5 [V ] which always gives a fourth-order convergence (e4 ∝ 1x4 ). Regions around the kink are smooth but involve grid nodes that lie on non-uniform grids. To illustrate how a high-order rate of convergence can be expected for D2,5 [V ], we consider such a point, say x = 95, which lies on a non-uniform grid at both the first and second stage refinement as shown in Fig. 1. It can be verified that e3 = 31x3 /16 and e3 = 1x3 /128 at the first and second stage refinement respectively. As such the ratio of convergence can be calculated as 31x3 128 16 1x3

= 24,

which implies a higher order convergence under such a local mesh refinement strategy. We also note that the node x = 90 which lies on a non-uniform grid at the first stage refinement will again be found on a uniform grid as we continue to double the number of nodes after R ≥ 2 such that e3 = 0. Similarly, x = 95 itself will be on a uniform grid for the case when R ≥ 3. It can also be argued that except for the additional points of refinement, we could have enforced the discretisation to be on a uniform grid as in [7,64]. For example, x = 95 is not on a uniform grid if we use {85, 90, 95, 97.5, 100} but is on a uniform grid if one uses {85, 90, 95, 100, 105} which forms part of our computational domain. This is a rather simple implementation issue that we do not deem worthy to pursue, and we prefer to keep the discretisation simple since points previously on non-uniform grids, like x = 95, will eventually form part of a uniform grid as we continue doubling the value of n. Since we have combined this non-uniform approximation in space with an exact in time stable ETI scheme, we can therefore expect to yield a high-order rate of convergence as demonstrated in all our numerical experiments. References [1] L.V. Ballestra, Repeated spatial extrapolation: An extraordinary efficient approach for option pricing, J. Comput. Appl. Math. 256 (2014) 83–91. [2] D.Y. Tangman, A. Gopaul, M. Bhuruth, Exponential time integration and Chebychev discretisation schemes for fast pricing of options, Appl. Numer. Math. 58 (2008) 1309–1319. [3] D.Y. Tangman, A.A.I. Peer, N. Rambeerich, M. Bhuruth, Fast simplified approaches to Asian option pricing, J. Comput. Finance 14 (2011) 3–36. [4] E. Pindza, K.C. Patidar, E. Ngounda, Implicit-explicit predictor–corrector methods combined with improved spectral methods for pricing European style vanilla and exotic options, Electron. Trans. Numer. Anal. 40 (2013) 268–293. [5] L.V. Ballestra, L. Cecere, A fast numerical method to price American options under the Bates model, Comput. Math. Appl. 72 (2016) 1305–1319. [6] W.F. Spotz, G.F. Carey, Extension of high-order compact schemes to time-dependent problems, Numer. Methods Partial Differential Equations 17 (2001) 657–672. [7] S.T. Lee, H. Sun, Fourth-order compact scheme with local mesh refinement for option pricing in jump-diffusion model, Numer. Methods Partial Differential Equations 28 (2012) 1079–1098.

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

18

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

[8] D.Y. Tangman, A. Gopaul, M. Bhuruth, Numerical pricing of options using high-order compact finite difference schemes, J. Comput. Appl. Math. 218 (2008) 270–280. [9] S. Mashayekhi, J. Hugger, Finite difference schemes for a nonlinear Black–Scholes model with transaction cost and volatility risk, Acta Math. Univ. Comenian. 84 (2015) 255–266. [10] J. Ankudinova, M. Ehrhardt, On the numerical solution of nonlinear Black–Scholes equations, Comput. Math. Appl. 56 (2008) 799–812. [11] J. Guo, W. Wang, An unconditionally stable, positivity-preserving splitting scheme for nonlinear Black–Scholes equations with transaction costs. Sci. World J. Volume 2014. [12] J. Guo, W. Wang, On the numerical solution of nonlinear option pricing equation in illiquid markets, Comput. Math. Appl. 69 (2015) 117–133. [13] B. Düring, M. Fournié, High order compact finite difference schemes for a nonlinear Black–Scholes equation, Int. J. Theor. Appl. Finance 6 (2003) 767–789. [14] R. Company, E. Navarro, J. Pintos, E. Ponsoda, Numerical solution of linear and nonlinear Black–Scholes option pricing equations, Comput. Math. Appl. 56 (2008) 813–821. [15] S. Heston, A closed-form solution for options with stochastic volatility with applications to bond and currency options, Rev. Financ. Stud. 6 (1993) 327–343. [16] R. Zvan, P.A. Forsyth, K.R. Vetzal, Penalty methods for American options with stochastic volatility, J. Comput. Appl. Math. 91 (1998) 199–218. [17] L.V. Ballestra, C. Sgarra, The evaluation of American options in a stochastic volatility model with jumps: An efficient finite element approach, Comput. Math. Appl. 60 (2010) 1571–1590. [18] W. Zhu, D.A. Kopriva, A spectral element approximation to price European options with one asset and stochastic volatility, J. Sci. Comput. 42 (2010) 426–446. [19] C. Chiarella, B. Kang, G.H. Meyer, The evaluation of barrier option prices under stochastic volatility, Comput. Math. Appl. 64 (2012) 2034–2048. [20] B. Düring, M. Fournié, High-order compact finite difference scheme for option pricing in stochastic volatility models, J. Comput. Appl. Math. 236 (2012) 4462–4473. [21] B. Düring, M. Fournié, C. Heuer, High-order compact finite difference schemes for option pricing in stochastic volatility models on non-uniform grids, J. Comput. Appl. Math. 271 (2014) 247–266. [22] C.W. Oosterlee, C.C.W. Leentvaar, X. Huang, Accurate American option pricing by grid stretching and high-order finite differences, Working papers, DIAM, Delft University of Technology, The Netherlands, 2005. [23] M.K. Bowen, R. Smith, Derivative formulae and errors for non-uniformly spaced points, Proceedings of the Royal Society A: Math. Phys. Eng. Sci. 461 (2005) 1975–1997. [24] R. Kangro, R. Nicolaides, Far field boundary conditions for Black–Scholes equations, SIAM J. Numer. Anal. 38 (2000) 1357–1368. [25] M. Ehrhardt, R.E. Mickens, A fast, stable and accurate numerical method for the Black–Scholes equation of American options, Int. J. Theor. Appl. Finance 11 (2008) 471–501. [26] H. Han, X. Wu, A fast numerical method for the Black–Scholes equation of American options, SIAM J. Numer. Anal. 41 (2003) 2081–2095. [27] H.Y. Wong, J. Zhao, An artificial boundary method for American option pricing under the CEV model, SIAM J. Numer. Anal. 46 (2008) 2183–2209. [28] A. Almendral, C.W. Oosterlee, Numerical valuation of options with jumps in the underlying, Appl. Numer. Math. 53 (2005) 1–18. [29] W. Wang, Y. Chen, Fast numerical valuation of options with jump under Merton’s model, J. Comput. Appl. Math. 318 (2017) 79–92. [30] Z. Cen, A. Le, A. Xu, Exponential time integration and second-order difference scheme for a generalized Black–Scholes equation, J. Appl. Math. 2012 (2012). http://dx.doi.org/10.1155/2012/796814. [31] N. Rambeerich, D.Y. Tangman, A. Gopaul, M. Bhuruth, Exponential time integration for fast finite element solutions of some financial engineering problems, J. Comput. Appl. Math. 224 (2009) 668–678. [32] P.A. Forsyth, K.R. Vetzal, Quadratic convergence for valuing American options using a penalty method, SIAM J. Sci. Comput. 23 (2002) 2095–2122. [33] L. Feng, V. Linetsky, Pricing options in jump-diffusion models: An extrapolation approach, Oper. Res. 56 (2008) 304–325. [34] F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Polit. Econ. 81 (1973) 637–654. [35] E.G. Haug, The Complete Guide To Option Pricing Formulas, second ed., McGraw-Hill, 2007. [36] H. Windcliff, P.A. Forsyth, K.R. Vetzal, Analysis of the stability of the linear boundary condition for the Black–Scholes equation, J. Comput. Finance 8 (2004) 65–92. [37] C.C. Chang, S.L. Chung, R.C. Stapleton, Richardson extrapolation technique for pricing American-style options, J. Futures Mark. 27 (2007) 791–817. [38] R. Lord, F. Fang, F. Bervoets, C.W. Oosterlee, A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes, SIAM J. Sci. Comput. 30 (2008) 1678–1705. [39] F. Fang, C.W. Oosterlee, Pricing early-exercise and discrete barrier options by Fourier–Cosine series expansions, Numer. Math. 114 (2009) 27–62. [40] H.E. Leland, Option pricing and replication with transactions costs, J. Finance 40 (1985) 1283–1301. [41] B. Bensaid, J.P. Lesne, H. Pagès, J. Scheinkman, Derivative asset pricing with transaction costs, Math. Finance 2 (1992) 63–86. [42] M. Jandačka, D. Ševčovič, On the risk-adjusted pricing-methodology-based valuation of vanilla options and explanation of the volatility smile, J. Appl. Math. 3 (2005) 235–258. [43] P. Wilmott, T. Hoggard, A.E. Whalley, Hedging option portfolios in the presence of transaction costs, Adv. Futures Opt. Res. 7 (1994) 21–35. [44] G. Barles, H.M. Soner, Option pricing with transaction costs and a nonlinear Black–Scholes equation, Finance Stoch. 2 (1998) 369–397. [45] R. Frey, P. Patie, Risk management for derivatives in illiquid markets: A simulation study, in: Advances in Finance and Stochastics, Springer, 2002, pp. 137–159. [46] U. Çetin, R. Jarrow, P. Protter, M. Warachka, Pricing options in an extended Black–Scholes economy with illiquidity: Theory and empirical evidence, Rev. Financ. Stud. 19 (2006) 493–529. [47] R. Frey, U. Polte, Nonlinear Black–Scholes equations in finance: Associated control problems and properties of solutions, SIAM J. Control Optim. 49 (2011) 185–204. [48] R.C. Merton, Option pricing when underlying stock returns are discontinuous, J. Financ. Econ. 3 (1976) 125–144. [49] S.P. Zhu, W.T. Chen, A predictor–corrector scheme based on the ADI method for pricing American puts with stochastic volatility, Comput. Math. Appl. 62 (2011) 1–26. [50] N. Clarke, K. Parrott, Multigrid for American option pricing with stochastic volatility, Appl. Math. Finance 6 (1999) 177–195. [51] P. Carr, A. Mayo, On the numerical evaluation of option prices in jump diffusion processes, Eur. J. Finance 13 (2007) 353–372. [52] N.J. Higham, The scaling and squaring method for the matrix exponential revisited, SIAM J. Matrix Anal. Appl. 26 (2005) 1179–1193. [53] C. Moler, C.V. Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Rev. 45 (2003) 3–49. [54] N.L. Trefethen, J.A.C. Weideman, T. Schmelzer, Talbot quadratures and rational approximations, BIT Numer. Math. 46 (2006) 653–670. [55] Y. d’Halluin, P.A. Forsyth, G. Labahn, A penalty method for American options with jump diffusion processes, Numer. Math. 97 (2004) 321–352. [56] K.J. in’t Hout, K. Volders, Stability of central finite difference schemes on non-uniform grids for the Black–Scholes equation, Appl. Numer. Math. 59 (2009) 2593–2609. [57] K.J. in’t Hout, K. Volders, Stability of central finite difference schemes for the Heston PDE, Numer. Algorithms 60 (2012) 115–133.

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.

M.J. Dilloo, D.Y. Tangman / Computers and Mathematics with Applications (

)

–

19

[58] F. Fang, C.W. Oosterlee, A novel pricing method for European options based on Fourier–Cosine series expansions, SIAM J. Sci. Comput. 31 (2008) 826–848. [59] D.M. Pooley, K.R. Vetzal, P.A. Forsyth, Convergence remedies for non-smooth payoffs in option pricing, J. Comput. Finance 6 (2003) 25–40. [60] J. Goodman, D.N. Ostrov, On the early exercise boundary of the American put option, SIAM J. Appl. Math. 62 (2002) 1823–1835. [61] Y.K. Kwok, Mathematical Models of Financial Derivatives, second ed., 2008, Chapter 5. [62] R.K. Coonjobeharry, D.Y. Tangman, M. Bhuruth, A two-factor jump-diffusion model for pricing convertible bonds with default risk, Int. J. Theor. Appl. Finance 19 (2016). http://dx.doi.org/10.1142/S0219024916500461. [63] D.Y. Tangman, N. Thakoor, K. Dookhitram, M. Bhuruth, Fast approximations of bond option prices under CKLS models, Finance Res. Lett. 8 (2011) 206–212. [64] N. Thakoor, D.Y. Tangman, M. Bhuruth, Fast valuation of CEV American options, Wilmott Mag. (2015) 54–61.

Please cite this article in press as: M.J. Dilloo, D.Y. Tangman, A high-order finite difference method for option valuation, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.05.006.