A radial basis function — Hermite finite difference approach to tackle cash-or-nothing and asset-or-nothing options

Journal of Computational and Applied Mathematics 368 (2020) 112523

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

A radial basis function — Hermite finite difference approach to tackle cash-or-nothing and asset-or-nothing options Yin Yang a , Fazlollah Soleymani b , Mahdiar Barfeie c , Emran Tohidi d ,

∗

a

Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Key Laboratory of Intelligent Computing & Information Processing of Ministry of Education, School of Mathematics and Computational Science, Xiangtan University, Xiangtan 411105, Hunan, China b Department of Mathematics, Institute for Advanced Studies in Basic Sciences (IASBS), Zanjan 45137–66731, Iran c Department of Mathematics and Computer Science, Sirjan University of Technology, Sirjan, Iran d Department of Mathematics, Kosar University of Bojnord, P.O. Box 94156–15458, Bojnord, Iran

article

info

Article history: Received 10 January 2019 Received in revised form 1 August 2019 MSC: 41A10 91G20 65M20 Keywords: RBF-HFD Fourth-order convergence Cash-or-nothing Asset-or-nothing Initial condition

a b s t r a c t A novel local meshfree method is proposed for simulating cash-or-nothing and assetor-nothing options, at which the initial condition is discontinuous. The novelty of the scheme is in the use of Hermite radial basis function (RBF) interpolation on set of local nodes. We also construct the RBF generated Hermite finite difference (HFD) formulas and subsequently obtain their associated differentiation matrices (DMs). The developed accuracy equipped with a θ -shift algorithm for smoothing the discontinuous initial conditions are illustrated with a series of financial tests. © 2019 Elsevier B.V. All rights reserved.

1. Preliminaries and background Mostly, models coming from practical applications in finance or engineering do not possess an exact solution or it is costly to be investigated and implemented, [1,2]. Hence, relying on computational algorithms for such purposes is unavoidable. Some of the most important milestones for solving governing equations in a modeling framework, are the finite difference (FD), finite element (FE), finite volume (FV) and spectral schemes, see [3–6] for some background. In such schemes, the nodes or elements distribution is not changing in a grid. A different and novel family of meshless (also known as meshfree) computational schemes has drawn the attention of researchers in this field. Kernel-based FD schemes show a localized meshless scheme, that turned out to be useful for some specific problems as well as high dimensional problems, see [7,8]. In this context, radial basis functions (RBFs) are famous functions as the kernels. The background notion of RBF interpolation was presented for the first time to estimate in irregular surfaces on scattered data sets, which was then discussed as an effective interpolation approach via a comparative investigation with several other techniques, see the discussions in [9]. Thereafter, many variants of RBFs have been studied and proposed and used to different scientific and engineering problems such as: scattered data interpolation, computational solution of PDEs and partial integro-differential equations (PIDEs), interested reader may refer to the book [10,11]. ∗ Corresponding author. E-mail addresses: [email protected], [email protected], [email protected] (E. Tohidi). https://doi.org/10.1016/j.cam.2019.112523 0377-0427/© 2019 Elsevier B.V. All rights reserved.

2

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

Technically speaking, in terms of interpolation for the function u(x) in M-dimensional space, the following linear combination should be considered [12]: u(x) ≃ Pn,O u(x) =

n ∑

ck φ (x − xk ) .

(1.1)

k=1

Enforcing the interpolation conditions: Pn,O u(xi ) = u(xi ),

i = 1, 2, . . . , n,

(1.2)

yields an n linear equation system having n unknowns {ck }. And so the estimates are obtained. Novel advances have made it doable to improve the grid-based restrictions of FD schemes by considering the kernelbased meshfree setting using RBFs. This procedure gives a meshfree implementation, and is called as the RBF-generated FD (RBF-FD) approach, see [13] for further discussions and background. To illustrate further, in the RBF-FD scheme for solving equations the discretization is calculated by filling several differentiation matrices (DMs) and assembling them into a single large, but sparse, coefficient matrix. The methodology of computing the weighting coefficients for RBF-FD approach was discussed in detail in [14]. 1.1. RBF-Hermite FD (RBF-HFD) method A recently developed member of the family of localized RBF meshfree methods is the RBF-HFD scheme. Recall that an application of the Hermite approach for solving ordinary differential equations (ODEs) and partial differential equations (PDEs) was discussed in the book [15]. Noting that the RBF-HFD is considered to be as an extension over the RBF-FD method, [16], and this makes it a new field of research as well. The objective of RBF-HFD scheme is to tackle a PDE computationally via calculating more accurate estimations to the differential operators in the given PDE. This developed accuracy is attained via additional information from the PDE itself, rather than increasing the stencil size, as is the common procedure to improve the accuracy with the RBF-FD procedures. A thorough discussion along with some application for this new topic has recently been provided in [17]. In this paper, we propose our RBF-HFD approach based on the multiquadric RBF defined as follows [18]:

φ (∥x − xi ∥2 ) =

√

c 2 + ∥x − xi ∥22 ,

i = 1, 2, . . . , m,

(1.3)

where c is the shape parameter. 1.2. Binary options Take into consideration the following Black–Scholes model in pricing a European vanilla option [19,20]:

∂ V (τ , S) 1 ∂ 2 V (τ , S) ∂ V (τ , S) = σ 2S2 + (r − q)S − rV (τ , S), (1.4) ∂τ 2 ∂ S2 ∂S where V (τ , S) is the option value, τ = T − t, T is the time of maturity, q is the dividend, r is the risk-free rate of interest and E is the price of strike. It is well-known that the pricing PDE (1.4), is accompanied with the following initial condition in terms of a call-type option: V (0, S) = (S − E )+ , and for a put-type option by: V (0, S) = (E − S )+ . A binary option (also known as digital option where the initial condition being the function of Heaviside) is an option of exotic nature, at which the payoff is either a constant amount of money or nothing at all, [21]. The two major kinds of such options are the cash-or-nothing and the asset-or-nothing options. The former gives a constant amount of cash if the option expires in-the-money, while the latter provides the value of the underlying security. Consider one asset digital (call) option with the following payoff for cash-or-nothing options [22]: V (0, S) =

{

1 0

if S > E , otherwise.

(1.5)

The cash-or-nothing option pays 1 dollar (or occasionally E dollars) at expiration once the option is in-the-money and 0 otherwise. The asset-or-nothing digital option is the same as cash-or-nothing, except for the fact that payoff is equal to the underlying price. The following is the initial condition for asset-or-nothing option: V (0, S) =

{

S 0

if S > E , otherwise.

(1.6)

The initial conditions for binary put options can similarly be defined. However, we note that the inequality > in (1.5)–(1.6) is sometimes defined as ≥ in the definitions, but we follow (1.5)–(1.6) in this work. The plot of these payoffs is given in Fig. 1.1.

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

Fig. 1.1.

3

The payoffs for cash-or-nothing binary option (left) and the asset-or-nothing binary option (right), when E = 40$ and smax = 3E.

To tackle binary options, the convergence rate of lattice schemes (viz., binomial or trinomial trees) has been analyzed in [23]. It was discussed that non-smooth payoffs reduce the expected convergence speed for (call) options. Since these methods could be considered as explicit FD methods, one could observe similar convergence speed difficulties as well. Pricing (1.4) is defined in an unbounded domain as follows:

{(τ , S) ∈ (0, T ] × [0, +∞)},

(1.7)

subject to an initial condition. To solve (1.4) computationally, we truncate (1.7) into a finite-sized computational one as follows: (τ , S) ∈ (0, T ] × [0, smax ] ≡ (0, T ] × Ω .

(1.8)

Following the discussions given in [24], it can be stated that the price of a European option based on (1.4) is estimated exponentially well by that of an associated barrier option (in log of the barrier), and hence a solution exists after truncating the domain. For further financial discussions, one may refer to [25–27]. 1.3. Motivation Basically, pricing options encounter with some discontinuous points in either the payoff condition (such as (1.5) and (1.6)) or its derivatives. This non-smooth criterion introduces quantization error [20], and could lead to serious degradation in the order of convergence of computational methods. Another note is poor approximations of the solution derivatives, i.e., Greeks, (delta, gamma and theta valuations) even though the option prices seem to be true. Option Greeks measure different factors that affect the price of an option contract. Here, we illustrate the key Greeks: 2 Delta ( ∂∂VS ) which gauges the likelihood that an option you are considering will expire in-the-money; Gamma ( ∂∂ SV2 ) that estimates how much the Delta will change when the stock price changes; and Theta (− ∂∂τV ) which gets a feel for how much value your option might lose each day as it approaches expiration. Armed with Greeks, an option trader can make more informed decisions about which options to trade, and when to trade them. Motivated by this fact for European vanilla options as well as for cash-or-nothing and asset-or-nothing options, a new scheme based on the RBF-HFD is constructed and proposed in this work. A simple yet efficient θ -shift strategy is also presented to handle the discontinuity presenting in the initial value. 1.4. Contributions The objectives of this work are given in what follows:

• To the best of authors knowledge, the weighting coefficients for the multiquadric RBF generated HFD method is computed for both the 1st and 2nd derivatives of a given sufficiently smooth function.

• It is proved that these approximation formulas possess a high speed of convergence four, which is quite accurate and have been obtained by applying only three adjacent points of grid. This high fourth order is valid for interpolation issues though for solving our practical financial problems and due to discontinuity in the domain (payoff), it may lose its accuracy. Nonetheless, the overall scheme based on such a weighting coefficients can compete to the well-known FD and RBF-FD schemes for tackling the similar problem. • Spatial discretization by the fourth-order RBF-HFD formulas is accompanied by a fast time-stepping solver to advance along time in order to tackle digital European options quickly. • A θ -shift strategy is imposed to smooth out the discontinuity of initial conditions and preserve the high rate of convergence in practice.

4

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

1.5. Organization The remaining parts of this work are organized as follows. Section 2 is dedicated to calculate the weighting coefficients for both the 1st and 2nd derivatives of a given function based on the methodology of RBF-HFD. The convergence of these approximation formulas is investigated and it is proved that they have fourth-order convergence rates only by applying three equally-spaced nodes. The application of the approximation formulas for spatially discretizing the Black–Scholes model under the binary options is treated in Section 3. The DMs are sparse and banded. It is shown how we must pursue the semi-discretization and get a set of coupled linear ODEs under appropriate initial/boundary condition. A time-stepping strategy for advancing along time as well as a simple yet effective way for smoothing the initial condition to preserve the convergence rate for digital options is pointed out as well. Computational aspects of the proposed new scheme are illustrated in Section 4 by a concentration on hedging parameters [28], viz, Greeks, which are important for understanding the sensitivity of option pricing. Final comments are given in Section 5 with some outlines for future studies. 2. Weights for the RBF-HFD method The importance of weights in the RBF-HFD approximations is that by considering only three adjacent equally-spaced nodes, one may get the high rate of convergence four for approximating 1st and 2nd derivatives. This is impossible in the FD and the RBF-FD approximations and the order four can be deduced, only if at least five nodes be taken into account. The notion of Hermite here allows us to consider not only the function values at each adjacent nodes, but also the derivative values in the adjacent nodes, [29]. Following this methodology, one is able to retrieve a fourth-order convergence order. Let us consider first the set of equally-spaced points {x1 , x2 , . . . , xn } for discretizing the computational domain [0, smax ], and then a stencil having three nodes as comes next:

{xi−1 , xi , xi+1 } = {xi − h, xi , xi + h},

h > 0.

(2.1)

Now, consider the formulation of the RBF-HFD approximation as follows: f ′ (xi ) ≃ fˆ ′ (xi ) = αi−1 f (xi−1 ) + αi f (xi ) + αi+1 f (xi+1 )

+ βi−1 f ′ (xi−1 ) + βi+1 f ′ (xi+1 ),

(2.2)

2 ≤ i ≤ n − 1.

The following theorem best addresses the error analysis of (2.2) for estimating the 1st derivative of a given function. Theorem 2.1. The RBF-HFD approximation for the 1st derivative of a sufficiently smooth function f using (2.2) is of fourth-order convergence rate (for interpolation), when the weights are defined by:

αi−1 = −αi+1 = −

h

−

2c 2

3 4h

,

(2.3)

αi = 0,

(2.4)

βi−1 = βi+1 = −

h2 2c 2

1

− .

(2.5)

4

Proof. To prove this, we substitute the function f by (1.3) at the stencil nodes to get the following linear system of coupled equations:

√ √ 2hβi−1 h = c 2 + 4h2 αi−1 + c 2 + h2 αi + √ + c αi+1 , √ 2 2 2 c + h2 √c + 4h ( 2 ) c c 2 + h2 αi + c + h2 αi−1 + c 2 αi+1 + h2 αi+1 + hβi+1 = hβi−1 , ( 2 ) √ c + 4h2 αi+1 + 2hβi+1 h = + c 2 + h2 αi + c αi−1 , √ √ 2 2 c 2 + h2 ( 2 c2 )+ 4h 2 h c + h αi + c 2hαi−1 c 2 βi−1 βi+1 , + = ( √ ( ) )3/2 + 3 / 2 2 2 2 2 2 2 c c + 4h c +h c + 4h c2

(

c 2 + h2

)3/2 = √

hαi

c2

+

2hαi+1

h2

+√

c2

+

4h2

+(

c 2 βi+1

c 2 + 4h2

)3/2 +

βi−1 c

.

(2.6) (2.7) (2.8) (2.9)

(2.10)

Solving (2.6)–(2.10) by considering c ≫ h and some simplifications, we get (2.3)–(2.5). This provides the following error equation to estimate the 1st derivative of the function:

( (3) ) ( ) f (xi ) 1 (5) ε (xi ) = − − f x h4 + O h5 , ( ) i 2 3c

120

(2.11)

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

5

wherein ε (xi ) = fˆ ′ (xi ) − f ′ (xi ). The error equation (2.11) proves a fourth-order convergence of the proposed approximation on equally-spaced grids. Noting that to ease up the derivation and understanding that how such high rates can be obtained by three nodes in this methodology, a program in the Wolfram Language Mathematica [30] is provided in the Appendix to confirm this rate of convergence and derivative. This completes the proof. □ One may remove the influence of the free shape parameter c by computing its limit when c → +∞. Hence, we obtain 3

αi−1 = −αi+1 = −

4h

1

, αi = 0, βi−1 = βi+1 = − .

(2.12)

4

To tackle the 2nd derivative case applying the RBF-HFD methodology, we consider (2.1) and write down: f ′′ (xi ) ≃ fˆ ′′ (xi ) = λi−1 f (xi−1 ) + λi f (xi ) + λi+1 f (xi+1 )

+ µi−1 f ′′ (xi−1 ) + µi+1 f ′′ (xi+1 ),

2 ≤ i ≤ n − 1.

(2.13)

Now, we give the weighting coefficients and their associated error equation in the following theorem. Theorem 2.2. The RBF-HFD formula for estimating the 2nd derivative of a sufficiently smooth function f using (2.2) is of fourth-order convergence speed (for interpolation), when the weighting coefficients are given by: 18159h2

333

6

λi−1 = λi+1 = − + + 2, 2 6125c 4 350c 5h ( ) 2 3885 4900 3 12106h − − c4 c2 h2 λi = , 6125 333h2

µi−1 = µi+1 = −

700c 2

−

1 10

(2.14)

(2.15)

.

(2.16)

Proof. Similar to the proof of Theorem 2.1, one may obtain the weighting coefficients (2.14)–(2.16) as long as c ≫ h by constructing the following set of five coupled linear equations at the stencil nodes: c2

(

c 2 + h2

)3/2 =

√

c 2 + 4h2 λi−1 +

µi+1 c 2 µi−1 , c 2 + h2 λi + ( )3/2 + c λi+1 + 2 2 c c + 4h

√

)3/2

( )2 ( )2 λi + c 2 + h2 λi−1 + c 2 + h2 λi+1 + c 2 µi−1 + c 2 µi+1 = , ( )3/2 c c 2 + h2 √ √ c 2 µi+1 µi−1 c2 , )3/2 = c 2 + h2 λi + c 2 + 4h2 λi+1 + ( )3/2 + c λi−1 + ( c c 2 + h2 c 2 + 4h2 ( )) ( λi 1 2µi−1 2λi−1 2 +√ + 3c ( = 0, h √ )5/2 − ( )5/2 c 2 + 4h2 c 2 + h2 c 2 + h2 c 2 + 4h2 ) ⎛ (( ⎞ ) 2 2 2 2 2 c + 4h λ − 3c µ 2 i + 1 i + 1 3c h λ i ⎠. −( +√ )5/2 = h ⎝ ( )5/2 c 2 + h2 c 2 + h2 c 2 + 4h2 1

(

c c 2 + h2

(2.17)

(2.18)

(2.19)

(2.20)

(2.21)

Using a computer algebra system [31], the solution to the system (2.17)–(2.20) can be found as (2.14)–(2.16). Writing Taylor expansions up to order two for the weighting coefficients and putting in (2.13) yields:

( ε (xi ) = −

3 6475c 2 f (4) (xi ) + 48424f ′′ (xi )

(

)

49000c 4

−

)

1 200

f

(6)

( ) (xi ) h4 + O h5 ,

(2.22)

wherein ε (xi ) = fˆ ′′ (xi ) − f ′′ (xi ). The error equation (2.22) reveals a fourth-order convergence speed for the proposed RBF-HFD approximation. □ In this case, the RBF-HFD formula could be simplified in the limit case when c → +∞, and we obtain:

λi−1 = λi+1 =

6 5h2

, λi = −

12 5h2

, µi−1 = µi+1 = −

1 10

.

(2.23)

Noticing that RBF-HFD methods are therefore basically more numerically effective than conventional RBF-FD schemes, as they can obtain higher accuracy and resolution for the same stencil size.

6

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

2.1. Shape parameters Two approaches for scaling the shape parameter are basically taken into account: (i). The first one is inversely proportional to the node distance h. (ii). The second one is to apply a fixed value by increasing the number of discretization nodes. The former selection keeps the ill-conditioning of the interpolation matrix constant, but imposes a stationary error of interpolation that does not decrease to zero in the limit as h goes to zero. Along these lines, a fixed c provides convergence by varying the number of discretization nodes, but the linear systems to calculate the weighting coefficients become ill-conditioned for small h, and convergence is lost because of the round-off errors. Hence, we rely on the first approach to propose c adaptively. Since for Vanilla European (call/put) options, the initial condition is only non-smooth, we may propose the following way for choosing the shape parameter: c

for European vanilla options

= 10h.

(2.24)

But for the binary options, at which the initial condition contains discontinuity, not large enough values for c make the Greeks to have oscillations or inaccuracies. Hence, we choose: c

for digital options

= 250h.

(2.25)

The relations (2.24)–(2.25) vary by increasing the number of discretization points. This means that for any number of discretization points, we obtain an h and then according to the type of the option we use the formulas (2.24)–(2.25) to get the shape parameters. Please see the last table of this paper for further understanding. In fact, these selections are in harmony with the theoretical discussions given in Theorems 2.1–2.2 by fulfilling the fact that c ≫ h. Besides, this value is changing by increasing the number of discretization nodes in the working domain. The performance of our method does not much depended to the choices (2.24)–(2.25), since we have considered the assumption c ≫ h. As such, considering larger values for c removes the effect of RBF approach and tend to regular HDF approximations while lower values for c might increase the accuracy but provide some instability in the numerical results as well. This phenomenon is known as the trade-off balance in this field, [10]. 3. Application in pricing (1.4) In Section 2, we have constructed the weighting coefficients to solve a PDE problem. The RBF-HFD formulas can be applied via semi-discretization approach to discretize the model (1.4) and obtain a set of coupled ODEs ready for the incorporation of time-stepping solvers for finding the final numerical solution. 3.1. Spatial discretizations The differential operator of (1.4) is discretized using DMs consisting of the weighting coefficients of Section 2. Accordingly, we may come up with the following DM including the weights for using the multiquadric RBF-HFD scheme for the convection term: i − j = 1, j − i = 1, otherwise,

⎧ 3 ⎨ − 2ch2 − 4h [ ] h 3 M2 = αi,j n×n = + 4h ⎩ 2c 2 0

(3.1)

and

[ ] M1 = βi,j n×n =

{

2

− 2ch 2 −

|i − j| = 1, otherwise.

1 4

0

(3.2)

The weights (2.6)–(2.8) are useful for the rows 2 to n − 1, while for the first and last rows of the DMs, the weights using a stencil having only two nodes should be constructed based on the FD methodology. We get the following weights: 1

1

h

h

α1,1 = αn,n−1 = − , α1,2 = αn,n =

−

h 2c 2

, β1,2 = βn,n−1 = 0.

(3.3)

Similarly, the DM associated with the diffusion term in (1.4) can be filled as comes next:

N2 = λi,j

[

] n×n

=

1 h2

⎧ ⎪ ⎪ ⎪ ⎨

(

2 3 12106h − 3885 − 4900 4 2 2 c

c

)

h

i = j,

6125

− ⎪ ⎪ ⎪ ⎩ 0

18159h2 6125c 4

+

333 350c 2

+

6 5h2

|i − j| = 1, otherwise,

(3.4)

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

7

and

[ ] N1 = µi,j n×n =

2

{

− 333h − 700c 2

|i − j| = 1, otherwise.

1 10

0

(3.5)

For its first and last rows, we have

λ1,1 = λn,n−1 =

−1 c2

, λ1,2 = λn,n =

2 c2

, µ1,2 = µn,n−1 = 0.

(3.6)

Gathering all the weights coming from the convection, diffusion and reaction terms of L in (1.4) together leads to the following system matrix: An×n =

1 2

σ 2 S2 (H1 N2 ) + (r − q)S(G1 M2 ) − rB,

(3.7)

wherein S H1 G1 B

diag (s1 , s2 , . . . , sn ) , I − M1 , I − N1 , G1 H1 ,

= = = =

(3.8)

and I is the n × n identity matrix. Here H1 , M1 , G1 , N1 are tri-diagonal matrices while S and I are diagonal. It is necessary to note that G1 and H1 satisfy G1 H1 = H1 G1 . Now, we get the following set of coupled homogeneous ODEs under the initial condition defined for European vanilla/binary options with appropriate side conditions imposed directly into the system matrix: BV˙ (τ ) = AV (τ ).

(3.9)

Note that the boundary conditions should be imposed only on A.

3.2. Time-stepping method The n × n matrices A and B are both penta-diagonal banded sparse matrices. We avoid multiplying both sides of (3.9) by the inverse of B so as to avoid encountering fully dense matrices upon implementation of the proposed procedure. Accordingly, we build our time-stepping solver on the sparse system (3.9) to keep the sparsity of the matrices when advancing in time. Upon discretizing in temporal variable τ , many schemes are existing, for example refer to [32]. Explicit methods are basically straightforward to implement, but suffer from stability problems. Implicit schemes are unconditionally stable, but are occasionally time-consuming because of solving nonlinear system of algebraic equations per step, see for example [33] and the references therein. Let us now denote vι as a numerical solution to the true value V (τι ) and choose k + 1 temporal nodes with the step size ∆τ = Tk . Now, we write the general s-stage Runge–Kutta method [34] for the approximate solution to an initial value problem at τι+1 = τι + ∆τ , (0 ≤ ι ≤ k) by: gi

=

vι + ∆τ

s ∑

kj ai,j ,

j=1

ki

=

vι+1

=

f (∆τ ci + τι , gi ) , vι + ∆τ

s ∑

(3.10)

bi ki ,

i=1

wherein Bf (·, ·) = AV (τ ). It is generally assumed that the row-sum conditions hold: ci =

s ∑

ai,j ,

i = 1, 2, . . . , s.

(3.11)

j=1

These conditions effectively determine the points in time, at which the function is sampled. Explicit Runge–Kutta methods are a special case where the matrix Λ = [ai,j ] to be strictly lower triangular. Thus, we may obtain a sixth-order

8

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

explicit Runge–Kutta method (RK6) as follows [35]:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ Λ=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

⎞

0

0

0

0

0

0

0

1

0

0

0

0

0

0 ⎟

3 8 8 27 3(3p1 −7) 392 −3(17p1 +77) 392 (7p1 +22) 12

1 8 2 27

0

0

0

0

0

8 27 −6(p1 −7) 49 (8p1 ) 49 2(7p1 −5) 9

0

0

0

0

3(p1 −21) 392 3(121p1 +21) 1960 −7(3p1 −2) 20

0

0

0 0

(p1 −7) 49

(−p1 −7)

−

49 2 3

(p1 +6) 5

0

−7(9p1 +49)

−7(p1 −7)

90

18

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(3.12)

0

with b = (9/180, 0, 64/180, 0, 49/180, 49/180, 9/180), C = (1, 1/2, 2/3, (7 − p1 )/14, (7 + p1 )/14, 1), and p1 = 211/2 . Notice that a consequence of explicitness is c1 = 0 in (3.11), so that the function is sampled at the beginning of the current integration step. Here, the sixth-order time-stepping solver consists of seven stages and reaches sixth order of convergence. The sixth order says that the error of local truncation is on the order of O(∆τ 7 ), while the total accumulated error is on the order of O(∆τ 6 ), [36, chapter 8]. A sixth-order explicit Runge–Kutta is suggested. The reasoning why such a high order scheme is relevant in this case is not only to get the advantage of explicit schemes by being straightforward in marching along time, but also to have larger domains for the choice of ∆τ . Furthermore, the reason of choosing such a higher order scheme is in the fact that the considered problem has a kink or discontinuity at the strike price in the payoff and since computational errors occurred at the beginning can remain until the end of the time stepping procedure. As such, we rely on a sixth-order solver to tackle this issue. Theorem 3.1. If the set of coupled homogeneous ODEs (3.9) satisfies the Lipschitz condition, then the temporal-solver (3.10)–(3.12) is conditionally stable along time. Proof. To find a stability condition, we proceed as follows, [37]. B is known as a mass matrix. This matrix is diagonalizable and positive definite (since G1 and H1 are positive definite matrices), which make all its eigenvalues to be positive. These could help to build a time-solver on (3.9) without the requirement of directly inverting B. The largest/smallest eigenvalues of B can be represented by: λmax (B) = λmax (G1 ) × λmax (H1 ) and similarly for λmin (B). On the other hand, the Toeplitz structure of our tri-diagonal matrices G1 and H1 , allows us to get the largest/smallest eigenvalues as follows [38]:

λmax (H1 ) λmax (G1 )

= =

1 + 2 cos 1 + 2 cos

√ ( π ) ( h2 (π )

√(

n

√ λmin (H1 )

=

λmin (G1 ) Therefore,

( 1−

1 2

√

=

1−

1 2

cos

1 − 2 cos

+

2c 2

n

(π )

333h2 700c 2

1 4

+

)2

, )2

1 10

c ≫ h,

,

(3.13)

(c 2 +2h2 )2 ,

c ≫ h,

4

n c √ ( π ) ( 333h2 700c 2

n

+

c ≫ h.

1 10

)2

,

c ≫ h.

)( √ ) ( √ 2 1 2 π h2 1 2 π 2 and λmin (B) = = 2 ( 333h + ) cos( ) + 1 ( + ) cos( ) + 1 10 n 4 n 700c 2 2c 2 ) ( ) √ 2 1 2 cos( πn ) × 1 − 2 ( 333h + 10 ) cos( πn ) . Incorporating the time-stepping solver (3.10) on the system 700c 2

λmax (B)

(c 2 +2h2 )2 c4

of ODEs (3.9) yields: vι+1 =

(

(∆τ B−1 A)2 (∆τ B−1 A)3 (∆τ B−1 A)4 2! 3! 4! (∆τ B−1 A)5 (∆τ B−1 A)6 (∆τ B−1 A)7 vι 5! 6! 2160

I + ∆τ B−1 A +

+

+

+

+

)

−

(3.14)

.

Thus, the numerical stability is reduced to:

⏐ ⏐ 2 3 4 5 6 7⏐ ⏐ ⏐1 + ∆τ ωi + (∆τ ωi ) + (∆τ ωi ) + (∆τ ωi ) + (∆τ ωi ) + (∆τ ωi ) − (∆τ ωi ) ⏐ ≤ 1, ⏐ 2 6 24 120 720 2126 ⏐

(3.15)

which is due to (3.14) for any ωi as the eigenvalue of B−1 A. Now by considering: 1

ωmax = (

√ 1−

1 2

(

) cos

2 c 2 +2h2

c4

)( (π ) n

√(

1−2

333h2 700c 2

+

1 10

)2

) λmax (A) , cos

(π ) n

(3.16)

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

Fig. 3.1.

9

Regions for eigenvalues in solving (3.9) using various temporal step sizes labeled above of each image.

the stability condition can now be represented as follows:

⏐ ⏐ ⏐1 + ∆τ ωmax +

(∆τ ωmax )2 2

ωmax ) + (∆τ 120 + 5

+

(∆τ ωmax )6 720

(∆τ ωmax )3 6

+ (∆τ ω24max ) ⏐ ωmax )7 ⏐ − (∆τ2126 ⏐ ≤ 1.

4

(3.17)

Or equivalently to the simplified form below:

⏐ ⏐ ⏐ 7.86071 ⏐ ⏐. ∆τ ≤ ⏐⏐ Re (ωmax ) ⏐

(3.18)

Note that the negative semi-definiteness of A makes all its eigenvalues to have negative real parts. Thus, the proposed scheme has numerical stability if the temporal step size ∆τ satisfies (3.17). The proof is ended. □ To illustrate how the stability region (3.17) works in practice by decreasing the temporal step-size ∆τ , we plot Fig. 3.1 at which blue areas are shown in the left side of the complex plane. The regions show where the eigenvalues of B−1 A (specially the maximum one based on (3.16)) must be located. The regions are getting enlarged rapidly by choosing smaller ∆τ . 3.3. Smoothing the initial condition Smoothing has long been a famous technique to attain stable convergence and in several cases preserve the optimal convergence speed in the presence of non-smoothness of the payoff. In finance, well-resulted techniques are averaging

10

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

the initial condition [39], smoothing in the neighborhood of the discontinuity in [40] by polynomial interpolation, and θ -shifting, [22]. In Algorithm 1, we provide a variant of θ -shift algorithm to smooth out the payoff. This is pursued by suggesting that shifting the grid such that discontinuities occur midway between grid points could improve accuracy. Recalling that authors in [20] pointed out that having the strike price take place midway among mesh nodes basically improves the efficiency of the FD-type schemes. Algorithm 1 A variant of θ -shift algorithm to smooth out the payoff.

• • • • •

Distribute the set of discretization nodes on the computational interval. Compute the payoff based on the contract conditions. Identify the first node after the discontinuity. Replace its value by a θ shifted value. The new shifted payoff can be applied as the modified initial conditions in order to preserve the convergence rate.

4. Computational aspects Numerical results and comparisons are provided for at-the-money options applying our proposed method (PM), the standard second-order uniform FD method with explicit Euler’s method (FDEM), the fourth-order FD method with Euler’s method for time-stepping (FD4), another implementation of the fourth-order uniform FD method but with the sixthorder time solver (3.10) named by FD4-2, and the Monte-Carlo simulation, in the programming package Mathematica 11.0, see [41, chapters 26–27]. Besides, we denote PM0 as the proposed method PM but without shifting strategy which is much useful for binary/digital options. This is pointed out that there are no comparisons with the results of the benchmark work [42], since the main attention in our work is to rely on digital options. The application of our new RBF-HFD weights for tackling digital options is novel. Three tests are considered in what follows: 1. A put option using the parameters: E = 100$, T = 1 year, σ = 0.4, r = 0.05, q = 0.6, smax = 3E .

(4.1)

A reference solution obtained via the exact solution as [20]: Vref (T , E) ≃ 41.3456074011, Gamma ≃ −0.4829545819, Delta ≃ 0.002744605, Theta ≃ −26.690915906. 2. A call option considering the following settings [22] for a cash-or-nothing option: E = 40$, T = 1 year, σ = 0.3, r = 0.05, q = 0, smax = 3E .

(4.2)

Here the reference value is: Vref (T , E) ≃ 0.48193918005 extracting from the Monte-Carlo simulation, and Gamma ≃ 0.03161943213, Delta ≃ −0.000834401461, Theta ≃ 0.02093498959. 3. A call option using the same settings as in Option 2, but for an asset-or-nothing option. Here the reference price is: Vref (T , E) ≃ 24.9700691162, and Gamma ≃ 1.889029461203, Delta ≃ −0.00175653872, Theta ≃ −2.40313231362. In Option 2, the reference solutions are obtained via the built-in command in the programming package Mathematica [43]:

ClearAll["Global‘*"]; sigma = 0.3; r = 0.05; q = 0; T = 1.; e = 40.; FinancialDerivative[{"BinaryCash", "European", "Call"}, {"StrikePrice" -> e, "Expiration" -> T}, {"InterestRate" -> r, "Volatility" -> sigma, "CurrentPrice" -> e, "Dividend" -> q}, {"Value", "Greeks"}] SetAccuracy[%, 12] For Option 3, we have obtained the reference solution, similarly. The boundary conditions for European put option are defined by: V (0, S) =E exp {−r τ }, lim V (S , τ ) =0,

S →∞

(4.3) (4.4)

and for European call option is V (0, S) =0,

(4.5)

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

11

Table 4.1 Convergence results for option 1. Method

n

∆τ

Re(ωmax )

V (T , E)

AE

FDEM

17 33 65 129 257

0.01 0.005 0.001 0.00025 0.0001

−47.6 −240.1 −1091.4 −4698.0 −19621.3

41.1463 41.3275 41.3386 41.3439 41.3455

1.99 × 10−1 1.80 × 10−2 6.98 × 10−3 1.66 × 10−3 9.25 × 10−5

3.53 × 10−3 8.35 × 10−4 2.09 × 10−4 5.29 × 10−5 1.37 × 10−5

1.45 × 10−4 1.58 × 10−5 4.90 × 10−6 1.18 × 10−6 7.67 × 10−8

2.90 × 10−1 6.73 × 10−2 1.65 × 10−2 4.01 × 10−3 1.17 × 10−3

0.01 0.02 0.03 0.12 0.50

FD4

17 33 65 129 257

0.01 0.005 0.001 0.00025 0.00005

−52.9 −293.6 −1393.2 −6123.5 −25846.6

41.4370 41.3969 41.3558 41.3482 41.3461

9.14 × 10−2 5.13 × 10−2 1.02 × 10−2 2.64 × 10−3 5.43 × 10−4

2.82 × 10−4 2.64 × 10−4 7.10 × 10−5 1.76 × 10−5 3.47 × 10−6

9.64 × 10−5 4.03 × 10−5 9.29 × 10−6 2.37 × 10−6 5.59 × 10−7

5.14 × 10−2 7.22 × 10−3 1.27 × 10−3 1.89 × 10−4 2.93 × 10−5

0.02 0.03 0.05 0.19 1.01

FD4–2

17 33 65 129

0.025 0.00625 0.00125 0.0004

−52.9 −293.6 −1393.2 −6123.5

41.3483 41.3528 41.3470 41.3460

2.70 × 10−3 7.23 × 10−3 1.45 × 10−3 4.37 × 10−4

2.07 × 10−4 2.16 × 10−4 6.12 × 10−5 1.52 × 10−5

4.77 × 10−5 1.17 × 10−2 4.24 × 10−6 1.10 × 10−6

1.42 × 10−2 1.55 × 10−5 2.50 × 10−3 7.55 × 10−4

0.03 0.07 0.10 0.35

PM

17 33 65 129

0.025 0.00625 0.00125 0.0004

−78.7 −378.8 −1686.8 −7201.5

41.3810 41.3552 41.3468 41.3455

3.54 × 10−2 9.60 × 10−3 1.25 × 10−3 3.61 × 10−5

6.98 × 10−3 1.55 × 10−3 3.61 × 10−4 7.32 × 10−5

1.32 × 10−4 3.98 × 10−5 1.01 × 10−5 2.32 × 10−6

7.30 × 10−2 1.69 × 10−2 2.63 × 10−3 4.16 × 10−4

0.05 0.06 0.12 0.42

Price

AE

Delta

AE

Gamma

AE

Theta

Time

Table 4.2 Convergence results for option 2. n

∆τ

Re(ωmax )

V (T , E)

AE

9 17 33 65 129 257

0.01 0.005 0.0025 0.001 0.0005 0.0001

−3.9 −26.4 −134.7 −613.5 −2642.3 −11036.6

0.5595 0.4373 0.5020 0.4720 0.4868 0.4794

7.75 × 10−2 4.46 × 10−2 2.00 × 10−2 9.88 × 10−3 4.95 × 10−3 2.46 × 10−3

9.30 × 10−3 1.86 × 10−3 4.77 × 10−4 1.08 × 10−4 1.63 × 10−5 9.59 × 10−6

1.85 × 10−4 5.89 × 10−4 9.04 × 10−5 8.07 × 10−5 3.09 × 10−5 1.79 × 10−5

5.07 × 10−2 3.26 × 10−2 1.11 × 10−2 5.49 × 10−3 2.71 × 10−3 1.37 × 10−3

0.01 0.01 0.02 0.03 0.07 0.49

9 17 33 65 129 257

0.01 0.005 0.0025 0.001 0.0005 0.00004

−2.6 −28.6 −164.0 −782.5 −3443.3 −14537.6

0.5482 0.4390 0.5013 0.4719 0.4868 0.4794

6.62 × 10−2 4.29 × 10−2 1.94 × 10−2 9.96 × 10−3 4.92 × 10−3 2.47 × 10−3

5.86 × 10−3 2.30 × 10−4 1.07 × 10−4 1.69 × 10−5 6.31 × 10−6 4.52 × 10−6

3.11 × 10−4 5.43 × 10−4 5.33 × 10−5 8.86 × 10−5 2.85 × 10−5 1.85 × 10−5

4.61 × 10−2 2.70 × 10−2 1.06 × 10−2 5.50 × 10−3 2.69 × 10−3 1.37 × 10−3

0.02 0.03 0.04 0.05 0.12 1.02

FD4–2

9 17 33 65 129

0.25 0.05 0.01 0.0025 0.000625

−2.6 −28.6 −164.0 −782.5 −3443.3

0.4880 0.4388 0.5013 0.4719 0.4868

6.07 × 10−3 4.30 × 10−2 1.94 × 10−2 9.96 × 10−3 4.92 × 10−3

1.10 × 10−2 3.02 × 10−4 1.39 × 10−4 2.97 × 10−5 3.46 × 10−8

5.29 × 10−4 5.47 × 10−4 5.44 × 10−5 8.84 × 10−5 2.87 × 10−5

5.60 × 10−2 2.72 × 10−2 1.06 × 10−2 5.51 × 10−3 2.69 × 10−3

0.02 0.05 0.07 0.09 0.26

PM

9 17 33 65 129

0.25 0.05 0.01 0.0025 0.000625

−8.2 −44.2 −210.0 −932.1 −3975.7

0.5173 0.4197 0.4839 0.4618 0.4813

3.54 × 10−2 6.21 × 10−2 1.96 × 10−3 2.00 × 10−2 6.13 × 10−4

1.12 × 10−3 1.21 × 10−3 2.19 × 10−4 7.50 × 10−5 1.46 × 10−5

2.73 × 10−4 4.17 × 10−4 7.78 × 10−6 1.42 × 10−4 5.50 × 10−6

3.57 × 10−5 3.50 × 10−2 1.06 × 10−3 1.10 × 10−2 3.35 × 10−4

0.02 0.05 0.07 0.13 0.44

Method

FDEM

FD4

lim V (S , τ ) =smax − E .

S →∞

Price

AE

Delta

AE

Gamma

AE

Theta

Time

(4.6)

Here, to impose the boundary conditions, we replace the first and last rows of A by zero, which stands for differentiating the boundary conditions in terms of τ with a scaling factor, see [34] to observe a background on this issue. This way of incorporating the artificial side conditions asymptotically converge to the boundaries and is fruitful for avoiding reaching inhomogeneous set of coupled ODEs or set of ODEs with time-varying system matrix. An interpolation inside the programming package is also imposed to find the required option values at any positions of the domain. We also set ‘‘PrecisionGoal− > 5,AccuracyGoal− > 5’’ in our codes to speed up the processes as much as possible, [30, chapter 1.11]. This implies that when imposing the stopping terminations and rounding the numbers inside our written codes for solving system (3.9), accuracies up to five decimals is enough. These would terminate the iterative computation sooner than it otherwise would, saving computer time. Time is reported in seconds. Absolute error (AE) which is the absolute value difference of the true and approximate solutions is used as a measure for error study. The numerical results are given in Tables 4.1–4.3 for Options 1–3, respectively. For binary options in Tables 4.2–4.3, the smoothing strategy of θ -shift is applied to preserve the rate of convergence. In Option 2, θ = 0.8 and in Option 3, θ = 7.

12

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

Table 4.3 Convergence results for option 3. Method

n

∆τ

Re(ωmax )

V (T , E)

AE

FDEM

17 33 65 129 257

0.01 0.005 0.002 0.0005 0.0001

−26.4 −134.7 −613.5 −2642.3 −11036.6

23.1679 25.7697 24.5744 25.1681 24.8712

1.80 × 100 7.99 × 10−1 3.95 × 10−1 1.98 × 10−1 9.87 × 10−2

7.81 × 10−2 1.83 × 10−2 4.11 × 10−3 6.92 × 10−4 3.78 × 10−4

2.33 × 10−2 3.69 × 10−3 3.20 × 10−3 1.25 × 10−3 7.10 × 10−4

1.30 × 100 4.46 × 10−1 2.19 × 10−1 1.09 × 10−1 5.43 × 10−2

0.01 0.02 0.03 0.07 0.50

FD4

17 33 65 129 257

0.01 0.005 0.002 0.0005 0.0001

−28.6 −164.0 −782.5 −3443.3 −14537.6

23.3055 25.7622 24.5760 25.1680 24.8713

1.66 × 100 7.92 × 10−1 3.94 × 10−1 1.97 × 10−1 9.87 × 10−2

9.44 × 10−3 2.72 × 10−3 2.19 × 10−4 2.71 × 10−4 1.36 × 10−4

2.03 × 10−2 2.52 × 10−3 3.43 × 10−3 1.17 × 10−3 7.28 × 10−4

1.04 × 100 4.33 × 10−1 2.17 × 10−1 1.08 × 10−1 5.42 × 10−2

0.02 0.03 0.04 0.11 0.56

FD4–2

17 33 65 129

0.0625 0.01 0.0025 0.000625

−28.6 −164.0 −782.5 −3443.3

23.2895 25.7604 24.5743 25.1677

1.68 × 100 7.90 × 10−1 3.95 × 10−1 1.97 × 10−1

1.51 × 10−2 5.24 × 10−3 1.20 × 10−3 2.54 × 10−5

2.07 × 10−2 2.55 × 10−3 3.44 × 10−3 1.17 × 10−3

1.06 4.34 × 10−1 2.18 × 10−1 1.08 × 10−1

0.02 0.06 0.09 0.24

PM0

17 33 65 129

0.0625 0.01 0.0025 0.000625

−44.2 −210.0 −932.1 −3975.7

24.2037 25.9851 24.6296 25.1817

7.66 × 10−1 1.01 × 100 3.40 × 10−1 2.11 × 10−1

1.56 × 10−3 1.80 × 10−4 1.27 × 10−3 1.17 × 10−4

4.05 × 10−3 6.58 × 10−3 2.30 × 10−3 1.46 × 10−3

5.28 × 10−1 5.56 × 10−1 1.87 × 10−1 1.16 × 10−1

0.08 0.09 0.15 040

PM

17 33 65 129

0.0625 0.01 0.0025 0.000625

−44.2 −210.0 −932.1 −3975.7

22.7509 25.1778 24.2278 24.9757

2.21 × 100 2.07 × 10−1 7.42 × 10−1 5.66 × 10−3

2.39 × 10−2 3.27 × 10−3 1.23 × 10−3 2.30 × 10−4

1.29 × 10−2 1.22 × 10−3 5.10 × 10−3 3.28 × 10−5

1.26 × 100 1.14 × 10−1 4.07 × 10−1 3.25 × 10−3

0.05 0.07 0.16 0.40

Price

AE

Delta

AE

Gamma

AE

Theta

Time

Table 4.4 Selecting the shape parameter in option 3 for PM. n

c

9 17 33 65 129 256

3750.000 1875.000 937.500 468.750 234.375 117.188

for digital option

Due to the kink in payoff higher order schemes such as FD4 cannot achieve their higher rates. This can be seen in Table 4.1 for FD4 in comparisons to PM. The real part of the largest eigenvalue of the matrix B−1 A denoted by Re(ωmax ) in absolute value sense is given for both compared methods to check the stability of the DMs constructed by the FD and RBF-HFD methods. In the implementation of the FD method, B is the identity matrix. Results show that the RBF-HFD and the FD schemes are almost similar for European vanilla options based on Table 4.1, but the computed option prices for digital options are much better for the cases, at which the payoff contains discontinuity. Convergence is restored for our proposed scheme time-stepping by smoothing the initial data. Moreover, accuracy is improved in contrast to the FD results. To also illustrate how the shape parameters can be chosen adaptively based on the number of nodes and the type of option, in Table 4.4, the results for the choice of c by increasing the number of discretization points are given. To illustrate the positivity of the option price, checking the results at boundaries as well as to observe correct behavior for hedging parameters, we plot these quantities in Fig. 4.1. Recall that the numerical valuation of digital options, and notably the approximation of their Greeks, delta, gamma and theta, require careful attention. We see that the error of the current RBF-HFD method is smaller than the error of the FD solution. However, the errors of both the FD and the RBF-HFD schemes on the grids are close for the first test. As such, the final usefulness of the proposed approach does not rely only based on the RBF-HFD approach in approximating the spatial variables, but also to the sixth-order time marcher as well as the θ -shifting strategy to smooth out the discontinuous. We end this section by noticing that the choice of an optimal θ in our shifting strategy to smooth out the discontinuous payoffs relies on both the computational interval and the number of discretization points. Accordingly, investigation for such an optimal choice is under investigation in our team. 5. Conclusions For high dimensional problems, models on irregular domain as well as when high order of accuracies are needed, the RBF-FD and the RBF-HFD formulas are attractive choices in contrast to the global RBFs methods. The purpose of this

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

13

Fig. 4.1. Computational solution of Option 3 by PM with n = 129, (first row-left). The (discrete) initial condition by shifting (first row-right). Solution at τ = T in (second row-left), and Delta in (second row-right). Gamma in (third row-left), and Theta in (third row-right).

paper was to contribute in finding the weighting coefficients of the RBF-HFD formulas based on the multiquadric RBF. After proving the high convergence orders of the constructed formulas, they have been used for pricing under binary options with discontinuous initial conditions. Discontinuities in the function of payoff introduce inaccuracies for computational methods, when pricing financial contracts. Particularly, large errors may take place in the approximation of the hedging parameters. While both cashor-nothing and asset-or-nothing options may be applied in analytical asset pricing, they are prone to fraud in their applications and thus removed by regulators in several jurisdictions as a form of gambling. Hence, their accurate pricing and their associated Greeks are of importance. A strategy for preserving the convergence order for non-smooth and discontinuous payoffs was also given. Moreover, several computational problems have been worked out by illustrating their option price graphs as well as their Greeks. The results have been compared by the second and fourth orders FD and Monte-Carlo simulations and showed to be useful. In this paper, we have done the calculations of different methods in a same environment by a laptop having Windows 7 Ultimate, Intel(R) Core(TM) i5-2430M CPU 2.40 GHz processor, HDD internal memory and 16.00 GB of RAM.

14

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

Further investigations could be done in extending these results for American binary options or for other important models such as the one in [44] in forthcoming works. Acknowledgments This work was supported by National Natural Science Foundation of China Project (11671342, 11771369, 11931003), Key Project of Hunan Provincial Department of Education (17A210, 16C0139), and Project of Scientific Research Fund of Hunan Provincial Science and Technology Department (2019YZ3003, 2018JJ2374, 2018WK4006). We are thankful to two anonymous referees whose suggestions and comments helped us to improve this work. Appendix In this appendix a simple yet effective Mathematica program is furnished to ease up the understanding of the procedure of computing the weighting coefficients in the approximation procedure of the first derivative of the function using only three nodes based on the methodology of the RBF-HFD. The fourth order of convergence can be obtained by implementing the last lines of the code.

ClearAll["Global‘*"] $Assumptions = {c > 0, h > 0, r > 0}; phi[r_] := (Sqrt[c^2 + r^2]) MQ1 = D[phi[r], r] // Simplify; MQ2 = D[phi[r], {r, 2}] // Simplify; u1[t_] := Subscript[a, i - 1] f[t - h] + Subscript[a, i] f[t] + Subscript[a, i + 1] f[t + h] + Subscript[b, i - 1] f’[t - h] + Subscript[b, i + 1] f’[t + h] eq1 = ((MQ1 /. {r -> -h}) == {(phi[r] /. {r -> -2 h}), (phi[ r] /. {r -> -h}), (phi[r] /. {r -> 0}) , (MQ1 /. {r -> -2 h}), (MQ1 /. {r -> 0})}.{Subscript[a, i - 1], Subscript[a, i], Subscript[a, i + 1], Subscript[b, i - 1], Subscript[b, i + 1]}) // Simplify eq2 = ((MQ1 /. {r -> 0}) == {(phi[r] /. {r -> -h}), (phi[r] /. {r -> 0}), (phi[ r] /. {r -> h}) , (MQ1 /. {r -> -h}), (MQ1 /. {r -> h})}.{Subscript[a, i - 1], Subscript[a, i], Subscript[a, i + 1], Subscript[b, i - 1], Subscript[b, i + 1]}) // Simplify eq3 = ((MQ1 /. {r -> h}) == {(phi[r] /. {r -> 0}), (phi[r] /. {r -> h}), (phi[ r] /. {r -> 2 h}) , (MQ1 /. {r -> 0}), (MQ1 /. {r -> 2 h})}.{Subscript[a, i - 1], Subscript[a, i], Subscript[a, i + 1], Subscript[b, i - 1], Subscript[b, i + 1]}) // Simplify eq4 = ((MQ2 /. {r -> -h}) == {(MQ1 /. {r -> -2 h}), (MQ1 /. {r -> \ -h}), (MQ1 /. {r -> 0}) , (MQ2 /. {r -> -2 h}), (MQ2 /. {r -> 0})}.{Subscript[a, i - 1], Subscript[a, i], Subscript[a, i + 1], Subscript[b, i - 1], Subscript[b, i + 1]}) // Simplify eq5 = ((MQ2 /. {r -> h}) == {(MQ1 /. {r -> 0}), (MQ1 /. {r -> h}), (MQ1 /. {r -> 2 h}) , (MQ2 /. {r -> 0}), (MQ2 /. {r -> 2 h})}.{Subscript[a, i - 1], Subscript[a, i], Subscript[a, i + 1], Subscript[b, i - 1], Subscript[b, i + 1]}) // Simplify {b1, A1} = Simplify@CoefficientArrays[{eq1, eq2, eq3, eq4, eq5}, {Subscript[a, i - 1], Subscript[a, i], Subscript[a, i + 1], Subscript[b, i - 1], Subscript[b, i + 1]}]; sol = Simplify@(Inverse[A1].(-b1)); taylororder = 2; sol2 = (Normal@Series[Simplify@(sol /. Thread[h -> e c]),

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

15

{e, 0, taylororder}] /. e -> h/c) // FullSimplify (Series[u1[t] /. {Subscript[a, i - 1] -> sol2[[1]], Subscript[a, i] -> sol2[[2]], Subscript[a, i + 1] -> sol2[[3]], Subscript[b, i - 1] -> sol2[[4]], Subscript[b, i + 1] -> sol2[[5]]}, {h, 0, 4}] // FullSimplify) /. t -> Subscript[x, i] References [1] H. Albrecher, A. Binder, V. Lautscham, P. Mayer, Introduction to Quantitative Methods for Financial Markets, Springer, Basel, 2013. [2] A. Basso, S. Funari, Introducing weights restrictions in data envelopment analysis models for mutual funds, Mathematics 6 (2018) 1–24, Art. ID: 164. [3] A. Itkin, Pricing Derivatives under Lévy Models: Modern Finite–Difference and Pseudo–Differential Operators Approach, Birkhäuser, Basel, 2017. [4] V.A. Kovtunenko, K. Kunisch, Revisiting generalized FEM: a Petrov–Galerkin enrichment based FEM interpolation for Helmholtz problem, Calcolo 55 (2018) 1–29, Art. ID: 38. [5] Y. Yang, Jacobi spectral Galerkin methods for fractional integro–differential equations, Calcolo 52 (2015) 519–542. [6] Y. Yang, Y. Huang, Y. Zhou, Numerical solutions for solving time fractional Fokker–Planck equations based on spectral collocation methods, J. Comput. Appl. Math. 339 (2018) 389–404. [7] S. Milovanović, L. von Sydow, Radial basis function generated finite differences for option pricing problems, Comput. Math. Appl. 75 (2018) 1462–1481. [8] L. von Sydow, S. Milovanović, E. Larsson, K. in ’t Hout, M. Wiktorsson, C.W. Oosterlee, V. Shcherbakov, M. Wyns, A. Leitao, S. Jain, T. Haentjens, J. Waldén, BENCHOP – SLV: the BENCHmarking project in option pricing – stochastic and local volatility problems, Int. J. Comput. Math. (2019) 1–14. [9] R.L. Hardy, Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res. 76 (1971) 1905–1915. [10] G.E. Fasshauer, Meshfree methods, in: Handbook of Theoretical and Computational Nanotechnology, Vol. 27, American Scientific Publishers, 2006, pp. 33–97. [11] M.N. Oqielat, O. Ogilat, Application of Gaussion radial basis function with cubic polynomial for modelling leaf surface, J. Math. Anal. 9 (2018) 78–87. [12] P. Assari, M. Dehghan, A meshless discrete collocation method for the numerical solution of singular–logarithmic boundary integral equations utilizing radial basis functions, Appl. Math. Comput. 315 (2017) 424–444. [13] I. Tolstykh, On using RBF–based differencing formulas for unstructured and mixed structured–unstructured grid calculations, in: Proc. 16th IMACS World Congress, Vol. 228, 2000, pp. 4606–4624. [14] G.B. Wright, Radial Basis Function Interpolation: Numerical and Analytical Developments (Ph.D. thesis), University of Colorado Boulder, 2003. [15] L. Collatz, The Numerical Treatment of Differential Equations, third ed., Springer, Berlin, 1966. [16] I.A. Tolstykh, A.D. Shirobokov, On using radial basis functions in a ‘‘finite difference mode’’ with applications to elasticity problems, Comput. Mech. 33 (2003) 68–79. [17] E. Lehto, V. Shankar, G.B. Wright, A radial basis function (RBF) compact finite difference (FD) scheme for reaction–diffusion equations on surfaces, SIAM J. Sci. Comput. 39 (2017) A2129–A2151. [18] V. Bayona, M. Moscoso, M. Kindelan, Gaussian RBF–FD weights and its corresponding local truncation errors, Eng. Anal. Bound. Elem. 36 (2012) 1361–1369. [19] L.J. Höök, G. Ludvigsson, L. von Sydow, The Kolmogorov forward fractional partial differential equation for the CGMY–process with applications in option pricing, Comput. Math. Appl. 76 (2018) 2330–2344. [20] D. Tavella, C. Randall, Pricing Financial Instruments: The Finite Difference Method, Wiley, New York, 2000. [21] K. Back, A Course in Derivative Securities Introduction to Theory and Computation, Springer–Verlag, Berlin, 2005. [22] D.M. Pooley, K.R. Vetzaly, P.A. Forsyth, Convergence remedies for non–smooth payoffs in option pricing, J. Comput. Finance 6 (2003) 25–40. [23] S. Heston, G. Zhou, On the rate of convergence of discrete–time contingent claims, Math. Finance 10 (2000) 53–75. [24] N.C.-H. Leung, C.C. Christara, D.-M. Dang, Partial differential equation pricing of contingent claims under stochastic correlation, SIAM J. Sci. Comput. 40 (2018) B1–B31. [25] L. Li, D. Taguchi, On a positivity preserving numerical scheme for jump–extended CIR process: the alpha–stable case, BIT (2019) http: //dx.doi.org/10.1007/s10543-019-00753-8. [26] S. Stojanovic, Computational Financial Mathematics using Mathematica: Optimal Trading in Stocks and Options, Springer Science + Business Media, LLC, New York, 2003. [27] R. Valkov, Convergence of a finite volume element method for a generalized Black–Scholes equation transformed on finite interval, Numer. Algorithms 68 (2015) 61–80. [28] K.-H. Kim, M.-G. Sin, Efficient hedging in general Black–Scholes model, Appl. Math. Finance 3 (2014) 1–9. [29] G.B. Wright, B. Fornberg, Scattered node compact finite difference–type formulas generated from radial basis functions, J. Comput. Phys. 212 (2006) 99–123. [30] M. Trott, The Mathematica Guidebook for Programming, Springer, New York, NY, USA, 2004. [31] J.G. Sánchez León, Mathematica Beyond Mathematics: The Wolfram Language in the Real World, Taylor & Francis Group, Boca Raton, 2017. [32] M. Sofroniou, G. Spaletta, Extrapolation methods in Mathematica, J. Numer. Anal. Indust. Appl. Math. 3 (2008) 105–121. [33] K. Eriksson, C. Johnson, A. Logg, Explicit time–stepping for stiff ODEs, SIAM J. Sci. Comput. 25 (2003) 1142–1157. [34] M. Sofroniou, R. Knapp, Advanced Numerical Differential Equation Solving in Mathematica, in: Wolfram Mathematica, Tutorial Collection, USA, 2008. [35] H.A. Luther, An explicit sixth–order Runge–Kutta formula, Math. Comp. 22 (1967) 434–436. [36] A.Ü. Keskin, Ordinary Differential Equations for Engineers, Springer International Publishing, AG, 2019. [37] L.F. Shampine, Numerical Solution of Ordinary Differential Equations, Chapman and Hall, 1994. [38] S. Noschese, L. Pasquini, L. Reichel, Tridiagonal Toeplitz matrices: Properties and novel applications, Numer. Linear Algebra Appl. 20 (2013) 302–326. [39] C.C. Christara, N.C.-H. Leung, Analysis of quantization error in financial pricing via finite difference methods, SIAM J. Numer. Anal. 56 (2018) 1731–1757. [40] Z. Cen, A. Le, A robust and accurate finite difference method for a generalized Black–Scholes equation, J. Comput. Appl. Math. 235 (2011) 3728–3733.

16

Y. Yang, F. Soleymani, M. Barfeie et al. / Journal of Computational and Applied Mathematics 368 (2020) 112523

[41] H. Ruskeepää, Mathematica Navigator, third ed., Academic Press, USA, 2009. [42] L. von Sydow, L.J. Höök, E. Larsson, E. Lindström, S. Milovanović, J. Persson, V. Shcherbakov, Y. Shpolyanskiy, S. Sirén, J. Toivanen, J. Waldén, M. Wiktorsson, J. Levesley, J. Li, C.W. Oosterlee, M.J. Ruijter, A. Toropov, Y. Zhao, BENCHOP – The BENCHmarking project in option pricing, Int. J. Comput. Math. 92 (2015) 2361–2379. [43] N.L. Georgakopoulos, Illustrating Finance Policy with Mathematica, Springer International Publishing, UK, 2018. [44] A. Itkin, P. Carr, Using pseudo–parabolic and fractional equations for option pricing in jump diffusion models, Comput. Econ. 40 (2012) 63–104.