A new method for evaluating options based on multiquadric RBF-FD method

Applied Mathematics and Computation 308 (2017) 130–141

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A new method for evaluating options based on multiquadric RBF-FD method Ahmad Golbabai∗, Ehsan Mohebianfar School of Mathematics, Iran University of Science & Technology, Narmak, Tehran, Iran

a r t i c l e

i n f o

MSC: 97N50 91G80 Keywords: Local meshless method Radial basis function Black–Scholes model Unconditional stability

a b s t r a c t In this paper, a new local meshless approach based on radial basis functions (RBFs) is presented to price the options under the Black–Scholes model. The global RBF approximations derived from the conventional global collocation method usually lead to ill-conditioned matrices. Employing the idea of local approximants of the finite difference (FD) method and combining it with the radial basis function (RBF) method can result in a local meshless approach such as RBF-FD. It removes the difficulty of ill-conditionness of the original method. The new proposed approach is unconditionally stable as it is shown by VonNeumann stability analysis. It is fast and produces high accurate results as shown in numerical experiments. Moreover, we took into account the variation of shape parameter and analyzed numerically the behavior of the RBF-FD method. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Options are among the most important financial derivatives. They are usually divided in two main categories: call and put options. A call option is the right to buy a particular asset for an agreed amount (strike price) at a specified time in the future (maturity). One would exercise the option at maturity if the stock is above the strike price and would not otherwise. If we use S to mean the stock price and K the strike price, then at maturity the option is worth,

gc (S ) = max(S − K, 0 ),

(1)

which is called the payoff function. On the other hand, a put option is the right to sell a particular asset for a strike price at a maturity. The holder of a put option wants the stock price to fall so that he can sell the asset for more than it is worth. The payoff function for a put option is,

g p (S ) = max(K − S, 0 ).

(2)

Now the option is only exercised if the stock falls below the strike price. The option is called an Eruopean option if it can be exercised only at maturity, and called American option if it can be exercised not only at maturity but also at any time prior to maturity. The options were first traded on 1973 where the Chicago board options exchange (CBOE) first created standardized and listed options. Initially, there were just calls while puts were not even introduced until 1977 [1,2].

∗

Corresponding author. E-mail addresses: [email protected] (A. Golbabai), [email protected] (E. Mohebianfar).

http://dx.doi.org/10.1016/j.amc.2017.03.019 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

131

If one own a call (put) option he wants the stock to rise (fall) as much as possible. The higher (lower) the stock price the greater his profit. But our decision to buy it will depend on how much it costs. The options are valuable and how much the options must be priced is one of the hot issues in the financial markets [1,2]. Black and Scholes considered the option pricing and published their seminal work in 1973 [3]. They presented a mathematical model for calculating the fair price of option. According to the Black–Scholes model, the price of the underlying asset S, satisfies under the risk neutral measure [1,2] the following stochastic differential equation,

ds = rSdt + σ SdW,

(3)

where r is the (constant) interest rate, σ is the (constant) volatility and W is a Wiener standard process. Let V(S, t) denote the price of a option on the underlying asset with maturity T and strike price K. For S ∈ (0, +∞ ) and t ∈ [0, T), the function V(S, t) satisfies the following one-dimensional linear parabolic partial differential equation called Black–Scholes equation,

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

(4)

1.1. European option The Black–Scholes Eq. (4) leads to an European call problem if one add the following boundary and final conditions,

V ( 0, t ) = 0,

V (S, t ) ∼ S

as

S → +∞,

V (S, T ) = gc (S ),

(5) (6)

where the payoff function gc (S) is defined as (1). It also leads to an European put problem if one considers the boundary and final conditions as follows,

V (0, t ) = Ke−r (T −t ) ,

lim V (S, t ) = 0,

S→+∞

V (S, T ) = g p (S ),

(7) (8)

where the payoff function gp (S) is defined as (2). The exact closed-form solution of the European case is provided and can be seen in [1,2,4]. 1.2. American option As mentioned before, American option can be exercised any time prior to maturity with optimal exercise stock value S f = S f (t ) which can be obtained by finding the intersection of both curves V(S, t) and payoff function. It makes American option a free boundary problem because a priori the location of the boundary Sf is unknown [7]. For more information, one can refere to [7]. So then, American option problem does not have an exact closed-form solution due to the unknown free boundary Sf (t) [1,4,7]. Hence, a variety of numerical approaches have been spread to solve this problem. They all suffer from free boudary difficulty. Some authors have been overpassed it by applying penalty method together with the finite volume or RBF methods for both one-dimensional and multi-dimensional problems [5,6]. To resolve this difficulty in our numerical results, we impose the boundary condition for American options as follows [4,7],

V (S, t ) = max{V (S, T ), V (S, t )},

(9)

where the V(S, T) is the payoff value given by the final conditions (6) and (8) for both call and put options respectively. Note that in the physical point of view, the difference between pricing the American options and the European options is the propagation process as an effect of the moving of the unknown free boundary Sf (t). This places an additional restriction at any time t on the solution that its value must be at least the payoff V(S, T) [8,9]. 1.3. Barrier option Options with the barrier feature known as Barrier options, which can take either European or American forms, are considered to be the simplest types of path dependent options. Barrier options distinctive feature is that the payoff depends not only on the final price of the underlying asset, but also on whether the asset price has breached some barrier level during the life of the option. An out-barrier option is one where the option is nullified prior to maturity if the underlying asset price touches the barrier level. The holder of the option may be compensated by a rebate payment for the cancellation of the option. An in-barrier option is one where the option only comes in existence if the asset price crosses the barrier level, though the holder has already paid the option premium up front. When the barrier is upstream with respect to the asset price, the barrier option is called an up-option; otherwise, it is called a down-option. One can identify four types of European barrier options for each types of call and put cases such as up-and-out, down-and-out, up-and-in and down-and-in [10]. As it is clear, the European barrier options have features similar to ordinary European options, except that when the asset price crosses the specified barrier level [10]. For example, a down-and-out put option is the same put option, except that it

132

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141 Table 1 Some well-known RBFs. Name

φ (r, c)

Gaussian (GA) Multiquadric (MQ) Inverse multiquadric (IMQ)

r exp √ (− c2 ) 2 + r2 c √ ( c2 + r2 )−1 2

becomes nullified when the asset price S falls below the downstream out-barrier level B. So in computational experiments, it is enough that one takes V(S, t) = 0 when S < B and solve an ordinary European put option otherwise. Some conventional approaches for solving the Black–Scholes equation are the finite difference method (FDM), the finite element method (FEM), the finite volume method (FVM), the fast Fourier transform [5,6,11–17], the binomial and triniomial trees [18–21] and asymptotic analytic approach [22]. As well as, authors in [23] used the variational iteration and homotopy perturbation methods for solving the Black–Scholes equation. Recently, some meshless methods based on radial basis functions (RBFs) have been presented. Authors in [4,24–26] have considered the global RBF method for solving the onedimensional Black–Scholes equation which apply the Kansa collocation method. They all suffer of dense RBF matrices with high condition numbers. Amani et al. [27] proposed local weak form meshless techniques based on the radial point interpolation (RPI) method and local boundary integral equation (LBIE) method. Moreover, some efforts have been done for solving two-dimensional case [28–30]. In [31], authors have proposed a local RBF method for pricing American options under Merton’s jump diffusion model. In this paper, we develop a new local RBF scheme based on FDM (RBF-FD) to solve the one-dimensional Black–Scholes equation. It results in more accurate solutions rather than global RBF scheme and overcomes the difficulty of dense and ill-conditioned RBF matrices. Moreover, we show that new scheme is unconditionally stable. The rest of this paper is as follows. In Section 2, the RBF-FD method is introduced briefly. It is employed to solve the Black–Scholes Eq. (4) in Section 3. An unconditional stability based on the Von-Neumann analysis is presented for the new proposed approach in Section 4. The numerical results and their interpretations are given in Section 5 where the new proposed scheme is compared with the global RBF scheme and some numerical efforts are done to consider numerically the effects of varying the shape parameter. Finally, we conclude this paper in Section 6. 2. The RBF approach based on FD scheme Some limitations in mesh generation, remeshing and making the approximation scheme in customary mesh-based methods such as FEM and FVM tend the general interests to use the meshless methods that remove the limitations of classical mesh-based methods. In meshless methods, only a cloud of points without any information about nodal connections is used to discretize the domain. One of the most favorite meshless methods is the RBF method which is constructed by radial kernels. They are (conditionally) positive definite, rotationally and translationally invariant. These features make its application straightway specially for approximating the solution of problems with high dimensions. The RBFs contain two useful specifications: a set of scattered centers with possibility of selecting their positions and existence of a free positive parameter known as the shape parameter [32]. Some of the well-known RBFs are mentioned in Table 1 where r is the Euclidean distance between any two points x, y ∈ Rd , i.e. r = x − y2 and c is the shape parameter [33–35]. The shape parameter, c, controls the shape of the basis functions and interchanges the error and stability of interpolation process. This behavior is manifested as Uncertainty Principle [36] which refers to the fact that the accuracy and well conditionness cannot be occured at the same time for a RBF approximant. In general, the RBF method builds the interpolant S(x) for approximating the function u(x) by a linear combination of translates of a single radial function on a set of centers XC = {xc1 , ..., xcN } ⊆ Rd [33,34]. So, for data values ucj = u(xcj ), j = 1, . . . , N, at given set of centers XC ,

u (x ) S (x ) =

N 

α j φ j ( x, c ),

(10)

j=1

where φ j (x, c ) = φ (x − xcj 2 , c ) for j = 1, . . . , N is a RBF corresponding to jth center. The unknown coefficients {α j }Nj=1 will be determined by collocating (10) in the same centers so that S(xci ) = uci , i = 1, . . . , N. It leads to a nonsingular algebraic linear system [33–35]. Kansa applied linear partial differential operator L on (10) to approximate Lu(x) in desired node xi which can be selected as a center [37,38] so that,

Lu(xi )

N 

α j Lφ j ( x i , c ) .

(11)

j=1

In the notation Lφ j (xi , c), it is understood that L is applied to the first variable, then evaluated. The relation (11) is a global approximation, i.e. for approximating of L at one node, all nodes in the domain are used. Many scientists have used this

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

133

Table 2 RBF-FD coefficients for the first and second derivatives. Node

First derivative N=3

Si Si Si Si Si

− 2h −h +h + 2h

Second derivative N=5



2

− 21h 1 + 2hc2 0   2 1 1 + 2hc2 2h



 1 +

1 12h − 32h

0 2 3h



1

 

8h2 c2 2 + 2ch2 2



1 + 2ch2   2 1 + 8ch2

− 121h

N=3 1 h2



N=5 2



1 + hc2   2 − h22 1 + hc2   2 1 1 + hc2 h2



2



h − 121h2 1 + 74 7 c2  4 37h2 1 + 3h2  14c2  h2 − 25h2 1 + 74 35c2  2 4 h 1 + 37 3h2  14c2 h2  − 121h2 1 + 74 7c 2

approach successfully where their results have showed advantages of RBF method for solving partial diffrential problems rather than traditional mesh-based methods [37–42]. The RBF matrices obtained from this global collocation method are very dense when a global RBF such as MQ is employed. It results in ill-conditioned RBF matrices specially as the dimension of matrices is increased. To overpass this issue and improve the accuracy of solution, one may apply the schemes proposed by authors [43] such as LU decomposition schemes, replacement of global solvers by block partitioning, matrix preconditioners and knot adaptivity that minimizes the total number of knots required in a simulation problem as well as the combination of these schemes can produce more accurate results [43]. On the other hand, one can employ a variable shape parameter strategy that assigns different shape parameters to each center. The effect of the variable shape parameter strategies in reducing condition number is considered in many papers [32,39,40,44–46]. One of the best recently idea to overcome this difficulty is to localize the collocation method. It produces sparse system matrices using local approximations. One can apply the finite diffrence (FD) method for this problem. It uses only some nodes surrounding of xi , called stencil, to approximate the derivetives of function u(x) at xi . Suppose  = {x1 , ..., xN } ⊆ Rd (i ) be a partition of domain and S I = {x1(i ) , ..., xN } ⊆  be a stencil corresponding to xi including NI nodes so that xi ∈ S I and I

NI ≤ N. So, for linear partial differential operator L we can use the local approximant,

Lu(xi )

NI  ω (j i) u(x(ji) ),

(12)

j=1 I where unknown weights {ω (j i ) } j=1 are usually computed using the polynomial interpolation [47,48]. This approach becomes very restricted in high dimensions because it is possible only for some kinds of structured nodes . It severely limits the geometric flexibility of the method [49]. This limitation can be eliminated by using the RBFs for calculating the unknown NI weights {ω (j i ) } j=1 called RBF-FD approach [50–52].

N

I I Using a set of RBFs as {φ j (x, c )} j=1 centered at S I for determining the unknown weights {ω (j i ) } j=1 in (12),

N

L φk ( x i , c ) =

NI  ω (j i) φ j (xk , c ),

N

k = 1, ..., NI ,

(13)

j=1

it leads to the following algebraic linear system,

ωI = [L]I ,

(14)

where  is the NI × NI interpolant matrix with entries,

φk j = φ j (xk , c ), k, j = 1, ..., NI ,

(15) (i ) NI

ωI is the NI × 1 coefficient vector including the weight coefficients {ω j } j=1 , named RBF-FD coefficients, and [L]I is the NI × 1 right-hand side vector including the values Lφ k (xi , c) for k = 1, ..., NI . The interpolant matrix  is nonsingular [33–35], therefore one can obtain the weights vector ωI ,

ωI = −1 [L]I .

(16)

In fact, the RBF-FD method completely works like the same FD method except in the way of calculating the weights I {ω (j i) }Nj=1 , so it can considered as an improved FD method. Bayona et al. [53] have provided some formulas for RBF-FD coefficients when the operator L is the first or the second order derivative operator. They have computed them for MQ function and some stencils including NI equispaced nodes using Maple software. Their results are presented in Table 2 for NI = 3, 5 in the limit c h. In this paper, we use their results for specifying the RBF-FD coefficients in numerical experiments. Some authors have employed this approach successfully for solving different problems and carried out some advantages of the RBF-FD rather than RBF and FD approaches [54–56]. It is well known that the local method lacks the spectral accuracy of the global RBF method [53]. Some efforts recently have been done to detect the dependence of the error with distance between nodes h, shape parameter c and number of centers N. Ding et al. [57] obtained an error estimate experimentally.

134

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

Authors in [58,59] showed that in the one dimensional case RBF-FD differentiation is equivalent to standard FD method when c → ∞. Authors in [52] introduced new local scheme as RBF-HFD scheme using the Hermite RBF interpolation method. Bayona et al. [53] also have considered the convergence properties of the RBF-FD formulas on uniform and scattered nodes and have analyzed the dependencies of the error with nodal distance h, shape parameter c and number of centers N. 3. Application of the RBF-FD approach to price the options under the Black–Scholes model In this section, we present a new RBF-FD scheme for solving the Black–Scholes Eq. (4). The scheme is proposed for the put options, however, it is valid for the call options. For the sake of simplicity, let us make the change of variable τ = T − t which allows us to rewrite the relations (4), (7) and (8) as follows,

∂ V (S, τ ) = LV (S, τ ). ∂τ

(17)

where,

LV (S, τ ) =

1 2 2 ∂ 2V (S, τ ) ∂ V (S, τ ) σ S + rS − rV (S, τ ), 2 ∂S ∂ S2

(18)

with modified boundary conditions,

V (0, τ ) = Ke−rτ ,

lim V (S, τ ) = 0,

(19)

S→+∞

and initial condition,

V (S, 0 ) = g p (S ),

(20)

where gp (S) is given by (2). 3.1. Time discretization At the first, we discretize (17) in time. Let Vn (S) denote the approximation of V(S, τ n ) for n = 0, 1, ..., M − 1, so that τn = n τ where τ is the uniform time step size. Note that the VM (S) is the desired option price which is obtained for time step τ M . We approximate the time derivative ∂ V(S, τ )∂ τ using the forward finite difference approximation as follows,

∂ V (S, τ ) V n+1 (S ) − V n (S ) = . ∂τ τ

(21)

To discretize the operartor (18) in time, the following time scheme known as θ -method is employed,

LV (S, τ ) = θ LV n+1 (S ) + (1 − θ )LV n (S ).

(22)

By substituting (21) and (22) in (17) one can obtain,

V n+1 (S ) − V n (S )

τ

= θ LV n+1 (S ) + (1 − θ )LV n (S ).

(23)

3.2. Spatial discretization Assume S ∈ [0, Smax ] and  = {S1 , ..., SN } be an uniform partition of it. In our numerical experiments, we select Smax (i ) sufficiently large to satisfy the right-end boundary condition (19). Furthermore S I = {S1(i ) , ..., SN } ⊆  be an uniform stencil corresponding to Si so that i = 2, ..., N − 1 and NI ≤ N. By collocating (23) at Si we have,

V n+1 (Si ) − V n (Si )

τ

where for k = n, n + 1,

LV k (Si ) =

= θ LV n+1 (Si ) + (1 − θ )LV n (Si ),



i = 2, ..., N − 1,

I

(24)



1 2 ∂ 2V k (S )  ∂ V k (S )  σ ( Si )2  + rSi(i)  − rV k (Si ). 2 2 ∂ S S=Si ∂S S=Si

(25)

The first and second order derivatives in (25) at Si can be approximated using a linear combinations of values of Vk (S) at S I so that,

 NI  ∂ V k (S )  = α (k,i)V k (S(ji) ),  ∂ S S=Si j=1 j

(26)

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

135

 NI  ∂ 2V k (S )  β (j k,i)V k (S(ji) ).  = 2 ∂S S=Si j=1

(27)

I I The unknown weights coefficients {α (j k,i ) } j=1 and {β (j k,i ) } j=1 are the RBF-FD coefficients corresponding to the first and sec-

N

N

N

I ond derivatives respectively. They can be determined by using RBFs centered at S I , i.e. {φ j (S, c )} j=1 . Using (26) and (27) in (25) and substituting it in (24) for time steps k = n, n + 1 result in the following time semi-discretization equation,

NI 



NI    n  ζ jn+1 (Si )V n+1 (S(ji) ) + ηin+1 (Si )V n+1 (Si ) = − ζ j (Si )V n (S(ji) ) − κin (Si )V n (Si ),

j=1 S(ji ) =Si

(28)

j=1 S(ji ) =Si

where for j = 1, ..., i − 1, i + 1, ...NI ,

1 2

ζ jn+1 (Si ) = θ σ 2 τ (Si )2 β (j n+1,i) + θ r τ Si α (j n+1,i) , 1 2

ζ jn (Si ) = (1 − θ )σ 2 τ (Si )2 β (j n,i) + (1 − θ )r τ Si α (j n,i) ,

(29) (30)

and for j = i,

1 2

ηin+1 (Si ) = θ σ 2 τ (Si )2 βi(n+1,i) + θ r τ Si αi(n+1,i) − (1 + θ r τ ),

(31)

  1 κin (Si ) = (1 − θ )σ 2 τ (Si )2 βi(n,i) + (1 − θ )r τ Si αi(n,i) + 1 − (1 − θ )r τ .

(32)

2

3.2.1. Crank–Nicolson scheme The time semi-discretization Eq. (28) for θ = 12 leads to the famous Crank–Nicolson (CN) scheme. In this case, for i = 2, ..., N − 1, Eq. (28) tends to the following matrix system,

AVn+1 = − ( BVn + C ),

(33)

where A and B are the (N − 2 ) × (N − 2 ) RBF-FD matrices which are NI -diagonal because of selecting uniform stencils S I . Their diagonal entries are given with (31) and (32) and their off-diagonal entries are specified using (29) and (30) respectively. C = [C1 , 0, ..., 0, CN−2 ]T is the boundary vector of dimension N − 2 so that,

C1 = ζ1n (S1 )V n (S1 ) + ζ1n+1 (S1 )V n+1 (S1 ),

(34)

CN−2 = ζNn (SN )V n (SN ) + ζNn+1 (SN )V n+1 (SN ),

(35)

where Vk (S1 ) and Vk (SN ) for k = n, n + 1 can be determined using the boundary conditions (19). Furthermore, both Vn and Vn+1 are column vectors of dimension N − 2 including the option price for Si , i = 2, ..., N − 1 and time steps n and n + 1 respectively. 3.2.2. Backward Euler scheme Using θ = 1 leads to the backward Euler (BE) scheme and for i = 2, ..., N − 1 Eq. (28) tends to the the following matrix system,

AVn+1 = − ( Vn + C ),

(36)

where the matrices A and the column vectors Vn and Vn+1 are defined as before. In this case, the entries C1 and CN−2 of the boundary vector C turn to,

C1 = ζ1n+1 (S1 )V n+1 (S1 ),

(37)

CN−2 = ζNn+1 (SN )V n+1 (SN ).

(38)

3.3. The initial step The option price vector VM can be determined using the CN scheme by solving (33) or the BE scheme by (36) M times for n = 0, ..., M − 1. As initial step one can use the initial condition (20) to calculate the initial vector V0 = [V 0 (S2 ), ..., V 0 (SN−1 )]T given for n = 0. Note that the system matrices A, B and C are calculated only one time in both schemes. They are constant for any time step n.

136

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

4. Stability analysis In this section, we prove that the proposed RBF-FD time semi-discretization (28) is unconditionally stable based on VonNeumann analysis when the CN scheme (θ = 12 ) or the BE scheme (θ = 1) is employed. For simplicity the stencils with three uniform nodes are used. For stencils with five uniform nodes, however, the stability achieved with some modifications. Let the stencils with three uniform nodes be used so that S I = {Si−1 , Si , Si+1 }, it turns (28) to,

ζin−1+1 (Si )V n+1 (Si−1 ) + ηin+1 (Si )V n+1 (Si ) + ζin+1+1 (Si )V n+1 (Si+1 ) = −ζin−1 (Si )V n (Si−1 ) − κin (Si )V n (Si ) − ζin+1 (Si )V n (Si+1 ), (39) where for time steps n, n + 1 and n = 0, 1, ..., M − 1, the coefficients

ηin+1 (Si )

ζ jn+1 (Si ) and ζ jn (Si ) for j = i − 1, i + 1 can be κin (Si ) can be determined using the relations

by (29) and (30) respectively, as well as the coefficient and (32), respectively. Assume that the solutions of (39) for k = n, n + 1 are as follows,

V k ( S j ) = ξ k e ϕ j ,

j = i − 1, i, i + 1,

specified (31) and

(40)

where  is the imaginary unit and ϕ is real. By substituting (40) in (39) one can obtain,

ζin−1+1 (Si )ξ n+1 eϕ (i−1) + ηin+1 (Si )ξ n+1 eϕ i + ζin+1+1 (Si )ξ n+1 eϕ (i+1) = −ζin−1 (Si )ξ n eϕ (i−1) − κin (Si )ξ n eϕ i − ζin+1 (Si )ξ n eϕ (i+1) .

(41)

With some simplification,

−ζin−1 (Si )e−ϕ − κin (Si ) − ζin+1 (Si )eϕ

ξ=

ζin−1+1 (Si )e−ϕ + ηin+1 (Si ) + ζin+1+1 (Si )eϕ

.

(42)

It leads to,

X − Y1 , Z + Y2

ξ= where,

 n   ζi+1 (Si ) + ζin−1 (Si ) cos ϕ + κin (Si ) ,

(44)

 n+1   ζi+1 (Si ) + ζin−1+1 (Si ) cos ϕ + ηin+1 (Si ) ,

(45)

X =− Z=

(43)

and,

Y1 = Y2 =



 ζin+1 (Si ) − ζin−1 (Si ) sin ϕ ,



 ζin+1+1 (Si ) − ζin−1+1 (Si ) sin ϕ .

(46) (47)

Using the RBF-FD coefficients in the Table 2 for uniform stencils, (43) results in,

 2 X 2 + Y12 ξ  = . 2 2 Z + Y2

(48)

This implies |ξ | < 1 for both CN and BE schemes, which shows that the proposed RBF-FD time semi-discretization (28) is unconditionally stable. 5. Numerical experiments In this section, we apply the new proposed RBF-FD method with both CN (θ = 12 ) and BE (θ = 1) schemes and stencils with three equispaced nodes i.e. NI = 3 to evaluate some options of different kinds. For this issue, some benchmarks of literature are employed from [4] that authors have presented a global RBF method. Although it has better results rather than conventional methods as FDM and FEM but suffers from dense and ill-conditioned RBF matrices with high condition numbers. This matter raises the computational cost and dependency of the accuracy of method to the value of shape parameter. The RBF-FD removes these problems but note that in all numerical experiments the value of shape parameter is selected sufficiently large to satisfy the limit condition c h due to the RBF-FD coefficients formulas in Table 2. All numerical results carried out for the MQ function using Matlab software with PC Laptop with an Intel(R) Core(TM)2 Duo CPU T6400 2 GHz 2 GB RAM.

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

137

Table 3 European problem parameters. Strike price

Volatility

Interest rate

Maturity

Maximum price

K = 10

σ = 0.2

r = 0.05

T = 0.5

Smax = 30

Table 4 Comparison the results of both RBF and RBF-FD methods for the European problem. N

(L2 )RBF

(L2 )RBF −F D CN

BE

CN

BE

41 61 81

3.9e − 03 1.5e − 03 7.0e − 03

2.1040e − 04 1.8166e − 04 3.0590e − 05

7.2375e − 04 7.0978e − 04 4.2307e − 04

0.011 0.015 0.024

0.006 0.009 0.011

CPUTime(s)

Table 5 Efficiency of the RBF-FD method using the CN scheme for the European problem. M

N

L∞

L2

ROC

CPUTime (s)

64 128 256 512 1024

64 128 256 512 1024

7.6196e − 03 5.3735e − 04 5.5813e − 05 5.4101e − 06 6.4431e − 07

2.4115e − 04 1.2640e − 05 7.1623e − 06 3.8866e − 07 2.2141e − 07

— 3.82 3.27 3.37 3.07

0.026 0.057 0.182 4.239 35.620

5.1. European case To show the efficiency of the new proposed RBF-FD method, we employ it for solving the European put problem with parameters given in Table 3 which has been deriven from [4]. All of comparisons are based on the root mean square and the maximum norms as follows,



1 L2 = N−2

N−1  i=2



2

M (S ) , Vexact (Si ) − Vapp. i



(49)



M L∞ = maxi=2,...,N−1 Vexact (Si ) − Vapp. (Si ),

(50)

M (S ) are the approximation of option values produced using the RBF-FD where Vexact (Si ) are the exact solutions [1,2] and Vapp. i or the RBF methods.

5.1.1. Example 1 In this example the RBF-FD method is compared with the RBF method for constant M = 5 and different N. The L2 norm for both approaches is presented in the Table 4 for a range of stock price as in [4]. As shown in this table, the new method with both CN and BE schemes is in good agreement with the exact solution and very fast because of including diagonal matrices. Moreover, it produces more accurate results rather than the RBF approach so that for N = 81 the root mean square error using the CN scheme is equal to 3.0590e − 05. 5.1.2. Example 2 Aim of this example is to show the high accuracy of the new RBF-FD approach. As well as, the rate of convergence (ROC) for numerical results is calculated using the following relation,

ROC = log2

Max. Error f or previous row

Max. Error f or current row

.

(51)

In the Tables 5 and 6, the results based on L2 and L∞ norms, the ROC and the CPU time are presented for the CN and BE schemes respectively when the number of time discretization steps, M, and the number of space discretization steps, N, are increased at the same time up to 1024. Following the numerical findings in the different errors, convinces us that the error of the proposed approach decreases very rapidly as the number of time and space discretization steps increase. Although both schemes are good performance but in the CN scheme, the produced results are a bit better than the BE scheme. Such that for M = N = 1024 the maximum and root mean square errors produced using the CN scheme are L∞ = 6.4431e − 07 and L2 = 2.2141e − 07 which are very significant. The ROC obtained approximately about 3.5 and 3 in the Table 5 and Table 6, respectively which verifies high accuracy of the new method. As well as the reported CPU time in two cases are very low such that for M = N = 1024 it is equal to 35.620 and 22.802, respectively. These CPU times illustrate that the new approach

138

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141 Table 6 Efficiency of the RBF-FD method using the BE scheme for the European problem. M

N

L∞

L2

ROC

CPUTime (s)

64 128 256 512 1024

64 128 256 512 1024

8.7501e − 03 9.2239e − 04 1.2427e − 04 2.4340e − 05 4.0041e − 06

2.6979e − 04 2.2469e − 05 1.0909e − 05 1.8424e − 06 7.4006e − 07

— 3.25 2.89 2.35 2.60

0.011 0.027 0.146 2.552 22.802

Table 7 American problem parameters. Strike price

Volatility

Interest rate

Maturity

Maximum price

K = 100

σ = 0.3

r = 0.1

T =1

Smax = e6

Table 8 Efficiency of the RBF-FD method for the American problem. S

Binomial n=10 0 0

RBF

RBF-FD CN

BE

80 100 120 140 160 180 200

20.2689 8.3348 3.2059 1.1789 0.4231 0.1502 0.0529

20.2777 8.3338 3.2108 1.1829 0.4254 0.1523 0.0570

20.2369 8.3187 3.1983 1.1755 0.4217 0.1500 0.0533

20.2163 8.2905 3.1818 1.1704 0.4222 0.1519 0.0511

Table 9 Barrier problem parameters. Strike price

Volatility

Interest rate

Maturity

Maximum price

K = 35

σ = 0.1

r = 0.1

T =1

Smax = 70

is very fast. Overall, as already pointed out, the presented numerical results show that the RBF-FD method is stable, fast and produces reasonable and accurate results. 5.2. American case As an American option, we apply the new RBF-FD method for solving the put problem with parameters given in Table 7 which has been deriven from [4,60]. To satisfy the optimal exercise for the valuation of the American put option, we use the update procedure as mentioned in [4]. We update the elements of Vn in the valuation of the European put option at each time level τ n for n = 1, ..., M by,

V n (Si ) = max{K − S(i ), V n (Si )},

i = 2, ..., N − 1.

(52)

5.2.1. Example 3 As mentioned before, an analytical solution for the American option is not possible. So we compare its values with the Binomial method [61] and RBF method [4] for a range of stock price S and constant time discretization steps so that M = 100. For the RBF method N = 101 is selected. We specify N = 161 so that the option values for the RBF-FD method obtain for the same stock price considered in [4] for both Binomial and RBF methods. According to the numerical results in the Table 8, the new method with both CN and BE schemes is comparable and in good agreement with the Binomial and RBF methods. 5.3. Barrier case The new proposed approach can be applied for all kinds of barrier options. So we employ a down-and-out barrier option with parameters given in the Table 9. 5.3.1. Example 4 For some stock price S and constant M = N = 701 the results are presented in the Table 10 which show the outcomes from the RBF-FD approach can be considered as good approximations of the exact ones. In Tables 11 and 12, the presented

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

139

Table 10 Results of both RBF and RBF-FD methods for the Barrier problem. S

Exact

RBF-FD CN

BE

30 32 34 36 38 40

0 0 0 0.1596 0.0467 0.0117

0 0 0 0.1595 0.0467 0.0117

0 0 0 0.1583 0.0459 0.0114

Table 11 Efficiency of the RBF-FD method using the CN scheme for the Barrier problem. N

M

L∞

L2

ROC

CPUTime (s)

64 128 256 512 1024

64 128 256 512 1024

7.9282e − 03 9.6230e − 04 8.7962e − 05 9.0519e − 06 1.0074e − 06

3.2556e − 04 5.5262e − 05 9.6017e − 06 1.6844e − 06 2.9669e − 07

— 3.04 3.45 3.28 3.17

0.036 0.058 0.206 4.731 37.892

Table 12 Efficiency of the RBF-FD method using the BE scheme for the Barrier problem. N

M

L∞

L2

ROC

CPUTime (s)

64 128 256 512 1024

64 128 256 512 1024

2.1060e − 03 3.0082e − 04 5.2654e − 05 9.4607e − 06 1.3919e − 06

8.7952e − 05 2.1880e − 05 6.8759e − 06 1.5008e − 06 9.9273e − 07

— 2.81 2.51 2.23 2.76

0.008 0.028 0.148 2.591 22.806

results illustrate that the new approach using both CN and BE schemes works well for the Barrier problem so that it is accurate, fast and reliable. 5.4. Effect of the shape parameter As mentioned before, the existence of the free shape parameter can be considered as an advantage of the RBFs. Without any additional computational cost or any change in other parameters of the method, varying the value of the shape parameter can cause more accurate results. Although this issue should be done with more carefully in the conventional RBF method due to Uncertainty Principle, but it can be done with more confidence in the new RBF-FD method because of the spars matrices. 5.4.1. Example 5 In the Fig. 1, the RBF-FD is compared with the standard FD approach for the same European case with parameters given in the Table 3. As it can be seen, the proposed approach has better performance when the shape parameter is sufficiently large to satisfy the coefficients in the Table 2 and it leads to the standard FD approach as c → ∞. Moreover, the maximum error can be optimized with an optimal shape parameter. 6. Conclusion The RBF-FD method as a new local meshless method is employed for evaluating different options under the Black– Scholes model such as the European, American and Barrier options. It goes through the difficulty of the original collocation approach and leads to well-posed diagonal matrices. The new approach can be used with the CN and BE schemes which its unconditional stability is shown using Von-Neumann analysis in those two cases. As test examples, we applied our new approach for some put options. We compared our results for the European case with the outcomes of the RBF method, for the American case with the Binomial and the RBF methods and for the Barrier case with the exact solutions. These comparisons showed that the new proposed method is competitive with traditional ones and it does not suffer from the difficulty of the global RBF method. It produces high accurate results with very low CPU time. Moreover, it has significant ROC which is calculated numerically. As well as, we considered the effect of varying

140

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

Standard FD RBF−FD(CN) RBF−FD(BE)

−2

Max. Error

10

−3

10

1

2

3

4 5 Shape parameter

6

7

8

Fig. 1. Comparison of the RBF-FD and the standard FD methods for M = N = 100.

the shape parameter briefly as an example. The numerical experiments show that the outcomes of the method can be optimized with an optimal shape parameter. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

P. Wilmott, Derivatives: The Theory and Practice of Financial Engineering, J. Wiley, 1998. J. Hull, Options, Futures, other Derivatives (seventh ed.). University of Toronto: Prentice Hall, 2002., F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Polit. Econ. 81 (1973) 637–659. Y.C. Hon, X. Mao, A radial basis function method for solving options pricing models, Financial Eng. 8 (1999) 31–49. P.A. Forsyth, K.R. Vetzal, Quadratic convergence for valuing American options using a penalty method, SIAM J. Sci. Comput. 23 (2002) 2095–2122. A.Q.M Khaliq, D.A Voss, S.H.K. Kazmi, Adaptive h-methods for pricing american options, J. Comput. Appl. Math. 222 (2008) 210–227. R. Seydel, Tools for Computational Finance, (third ed.). Springer, Berlin, 2012, Z. Wu, Y.C. Hon, Convergence error estimate in solving free boundary diffusion problem by radial basis functions method, Eng. Anal. Bound. Elem. 27 (2003) 73–79. K.H.M Khabir, Numerical singular perturbation approaches based on spline approximation methods for solving problems in computational finance, 2011, PhD Thesis, University of the Western Cape. Y.K. Kwok, Mathematical Models of Financial Derivatives, Springer Berlin, Heidelberg, 2008. M. Brennan, E. Schwartz, Finite difference methods and jump processes arising in the pricing of contingent claim: a synthesis, J. Financial Quantitative Anal. 13 (1978) 461–474. R. Zvan, P.A. Forsyth, K. Vetzal, A finite volume approach for contingent claims valuation, IMA J. Num. Anal 21 (2001) 703–731. X. Wu, W. Kong, A highly accurate linearized method for free boundary problems, Comput. Math. Appl. 50 (2005) 1241–1250. J. Zhao, M. Davison, R.M. Corless, Compact finite difference method for American option pricing, J. Comput. Appl. Math. 206 (2007) 306–321. N. BF, O. Skavhaug, A. Tveito, Penalty methods for the numerical solution of American multi-asset option problems, J. Comput. Appl. Math. 222 (2008) 3–16. M. Yousuf, K. AQM, B. Kleefeld, The numerical approximation of nonlinear Black–Scholes model for exotic path-dependent American options with transaction cost, Int. J. Comput. Math. 89 (2012) 1239–1254. L.V. Ballestra, L. Cecere, A numerical method to compute the volatility of the fractional Brownian motion implied by american options, Int. J. Comput. Math. 26 (2013) 203–220. J.C. Cox, S.A. Ross, M. Rubinstein, Option pricing: a simplified approach, J. Financial Econ. 7 (1979) 229–263. M. Broadie, J. Detemple, American option valuation: new bounds, approximations, and a comparison of existing methods, Rev. Financial Stud. 9 (1996) 1211–1250. M. Gaudenzi, F. Pressacco, An efficient binomial method for pricing american put options, Decis. Econ. Finance 4 (2003) 1–17. S. Chung, C.C. Chang, R.C. Stapleton, Richardson extrapolation techniques for the pricing of American-style options, J. Futures Market. 27 (2007) 791–817. S.M.M. Kazemi, M. Dehghan, A.F. Bastani, Asymptotic expansion of solutions to the Black–Scholes equation arising from american option pricing near the expiry, J. Comput. Appl. Math 311 (2017) 11–37. M. Dehghan, S. Pourghanbar, Solution of the Black–Scholes equation for pricing of barrier option. Z. Naturforschung a, 2011, 66a, 289–296 Y.C. Hon, A quasi-radial basis functions method for American options pricing, Comput. Math. Appl. 43 (2002) 513–524. M.D. Marcozzi, S. Choi, C.S. Chen, On the use of boundary conditions for variational formulations arising in financial mathematics, Appl. Math. Comput. 124 (2003) 197–214. A. Golbabai, D. Ahmadian, M. Milev, Radial basis functions with application to finance: American put option under jump diffusion, Math. Comput. Model. 55 (2012) 1354–1362. J. Amani Rad, P. Kourosh, S. Abbasbandy, Local weak form meshless techniques based on the radial point interpolation (RPI) method and local boundary integral equation (LBIE) method to evaluate european and american options, Commun. Nonlinear Sci. Numer. Simul. 22 (2015) 1178–1200. A. Khaliq, G. Fasshauer, D. Voss, Using meshfree approximation for multi-asset American option problems, J. Chin. Inst. Eng. 27 (2004) 563–571.

A. Golbabai, E. Mohebianfar / Applied Mathematics and Computation 308 (2017) 130–141

141

[29] U. Petterssona, E. Larssona, G. Marcussonb, J. Perssonc, Improved radial basis function methods for multi-dimensional option pricing, J. Comput. Appl. Math. 222 (2008) 82–93. [30] L.V. Ballestra, G. Pacelli, Pricing European and american options with two stochastic factors: a highly efficient radial basis function approach, J. Econ. Dyn. Control 37 (2013) 1142–1167. [31] A.A.E.F Saib, D.Y. Tangman, M. Bhuruth, A new radial basis functions method for pricing American options under Merton’s jump-diffusion model, Int. J. Comput. Math. 89 (2012) 1164–1185. [32] A. Golbabai, E. Mohebianfar, H. Rabiei, On the new variable shape parameter strategies for radial basis functions, Comput. Appl. Math. 34 (2015) 691–704. [33] M.D. Buhmann, Radial Basis Functions: Theory, Implementation, Cambridge University Press, University of Gissen, 2004. [34] H. Wendland, Scattered Data Approximation. Cambridge University Press, 2005, [35] S.A. Sarra, E.J. Kansa, Multiquadric Radial Basis Function Approximation Methods for the Numerical Solution of Partial Differential Equations. Tech Science Press, 2009 [36] R. Schaback, Improved error bounds for scattered data interpolation by radial basis functions, Math. Comput. 68 (1999) 201–206. [37] E.J. Kansa, Multiquadrics-a scattered data approximation scheme with applications to computational fluid dynamics-I surface approximations and partial derivative estimates, Comput. Math. Appl. 19 (1990) 127–145. [38] E.J. Kansa, Multiquadrics-a scattered data approximation scheme with applications to computational fluid dynamics-II solutions to parabolic, hyperbolic and elliptic partial differential equations, Comput. Math. Appl. 19 (1990) 147–161. [39] A. Golbabai, H. Rabiei, Hybrid shape parameter strategy for the RBF approximation of vibrating systems, Int. J. Comput. Math. 89 (2012) 2410–2427. [40] A. Golbabai, H. Rabiei, A meshfree method based on radial basis functions for the eigenvalues of transient stokes equations, Eng. Anal. Bound. Elem. 36 (2012) 1555–1559. [41] M. Dehghan, M. Abbaszadeh, A. Mohebbi, A meshless technique based on the local radial basis functions collocation method for solving parabolic– parabolic Patlak-Keller-Segel chemotaxis model, Eng. Anal. Bound. Elem. 56 (2015) 129–144. [42] M. Dehghan, V. Mohammadi, The numerical simulation of the phase field crystal (PFC) and modified phase field crystal (MPFC) models via global and local meshless methods, Comput. Methods Appl. Mech. Eng. 298 (2016) 453–484. [43] E.J. Kansa, Y.C. Hon, Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations, Comput. Math. Appl. 39 (20 0 0) 123–137. [44] R.E. Carlson, T.A. Foley, The parameter r2 in multiquadric interpolation, Comput. Math. Appl. 21 (1991) 29–42. [45] E.J. Kansa, R.E. Carlson, Improved accuracy of multiquadric interpolation using variable shape parameters, Comput. Math. Appl. 24 (1992) 99–120. [46] S.A. Sarra, D. Sturgill, A random variable shape parameter strategy for radial basis function approximation methods, Eng. Anal. Bound. Elem. 33 (2009) 1239–1245. [47] B. Fornberg, Calculation of weights in finite difference formulas, SIAM Rev., 1998, 40, 685–691 [48] M. Dehghan, A computational study of the one-dimensional parabolic equation subject to nonclassical boundary specifications, Numer. Methods Part Differ. Equ. 22 (1) (2006) 220–257. [49] W. Cheney, An introduction to approximation theory (2nd ed.). AMS Cheslea Publishing: American Mathematical Society, New York, 20 0 0 [50] A.I. Tolstykh, D.A. Shirobokov, On using radial basis functions in a finite-difference mode with applications to elasticity problems, Comput. Mech. 33 (2003) 68–79. [51] C. Shu, H. Ding, K.S. Yeo, Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations, Comput Meth Appl Mech Eng 192 (2003) 941–954. [52] G.B. Wright, B. Fornberg, Scattered node compact finite difference-type formulas generated from radial basis functions, J. Comput. Phys. 212 (2006) 99–123. [53] V. Bayona, M. Moscoso, M. Carretero, RBF-FD formulas and convergence properties, J. Comput. Phys. 229 (2010) 8281–8295. [54] J.G. Wang, G.R. Liu, A point interpolation meshless method based on radial basis functions, Int. J. Numer.Methods Eng. 54 (2002) 1623–1648. [55] S. Chantasiriwan, Investigation of the use of radial basis functions in local collocation method for solving diffusion problems, Int. Commun. Heat Mass Transf. 31 (2004) 1095–1104. [56] B. Sarler, R. Vertnik, Meshfree explicit local radial basis function collocation method for diffusion problems, Comput. Math. Appl. 51 (2006) 1260–1282. [57] H. Ding, C. Shu, D.B. Tang, Error estimates of local multiquadric-based differential quadrature (LMQDQ) method through numerical experiments, Int. J. Numer. Methods Eng. 63 (2005) 1513–1529. [58] T.A. Driscoll, B. Fornberg, Interpolation in the limit of increasingly flat radial basis functions, Comput. Math. Appl. 43 (2002) 413–422. [59] B. Fornberg, G.B. Wright, E. Larsson, Some observations regarding interpolants in the limit of flat radial basis functions, Comput. Math. Appl. 47 (2004) 37–55. [60] L. Wu, Y.K. Kwok, A front-fixing finite difference method for the valuation of American options, J. Financial Eng. 6 (1997) 83–95. [61] D.J. Higham, Nine ways to implement the binomial method for option valuation in MATLAB, SIAM Rev., 2002, 44, 4, 661–677