An Adaptive and Explicit Fourth Order Runge-Kutta-Fehlberg Method Coupled with Compact Finite Differencing for Pricing American Put Options
11 An Adaptive and Explicit Fourth Order Runge-Kutta-Fehlberg Method Coupled with Compact Finite Differencing for Pricing American Put Options.
Chinonso I. Nwankwo a, *, Weizhong Dai b a Department of Mathematics , Statistics, and Computer Science, University of Illinois at Chicago, Chicago, IL 60607, USA b Department of Mathematics and Statistics, Louisiana Tech University, Ruston LA 71272, USA * Corresponding author, [email protected], [email protected]
Abstract
We propose an adaptive and explicit fourth-order Runge-Kutta-Fehlberg method coupled with a fourth-order compact scheme to solve the American put options problem. First, the free boundary problem is converted into a system of partial differential equations with a fixed domain by using logarithm transformation and taking additional derivatives. With the addition of an intermediate function with a fixed free boundary, a quadratic formula is derived to compute the velocity of the optimal exercise boundary analytically. Furthermore, we implement an extrapolation method to ensure that at least, a third-order accuracy in space is maintained at the boundary point when computing the optimal exercise boundary from its derivative. As such, it enables us to employ fourth-order spatial and temporal discretization with Dirichlet boundary conditions for obtaining the numerical solution of the asset option, option Greeks, and the optimal exercise boundary. The advantage of the Runge-Kutta-Fehlberg method is based on error control and the adjustment of the time step to maintain the error at a certain threshold. By comparing with some existing methods in the numerical experiment, it shows that the present method has a better performance in terms of computational speed and provides a more accurate solution.
Keywords:
American put options, logarithmic transformation, optimal exercise boundary, compact finite difference method, Runge-Kutta-Fehlberg method, fixed free boundary Introduction
American style option, written on an asset ๐ ๐ก with the strike price, ๐พ and expiration time ๐ differs from the European option due to the early (optimal) exercise boundary which leads to a free boundary problem Let ๐(๐, ๐ ) denote the option price, ๐ ๐ (๐) represent the optimal exercise boundary and ๐ = ๐ โ ๐ก . Then,
๐(๐, ๐) satisfies the coupled free boundary value problem: โ ๐๐(๐, ๐)๐๐ + 12 ๐ ๐ ๐ ๐(๐, ๐)๐๐ + ๐๐ ๐๐(๐, ๐)๐๐ โ ๐๐(๐, ๐) = 0, for ๐ > ๐ ๐ (๐), (1๐) ๐(๐, ๐) = ๐พ โ ๐, for ๐ < ๐ ๐ (๐). (1๐) Here, the initial and boundary conditions are given as:
๐(๐, 0) = ๐๐๐ฅ( ๐พ โ ๐, 0), ๐ ๐ (0) = ๐พ; (1๐) ๐(๐ ๐ , ๐) = ๐พ โ ๐ ๐ (๐), ๐(0, ๐) = ๐พ, ๐(โ, ๐) = 0, ๐๐๐ ๐(๐ ๐ , ๐) = โ1, (1๐) This early exercise boundary presents an advantage and a challenge in the valuation of American options and solving (1), respectively. In terms of advantage, it provides a possibility to exercise the options early, however, we have some level of complexity. This is because the early exercise boundary and the American option values are simultaneously obtained (McKean, 1965) when solving the free boundary problem. It is well known that due to this complexity, there is no closed-form or analytical formula for evaluation of the American option. Hence, numerical, semi-numerical, and analytical approximation present a choice for solving (1). Several numerical methods have been proposed for solving the American options problem with the front-fixing approach. In particular, the second-order explicit and implicit finite difference schemes have been used (Company et al., 2014; Company et al., 2016; Nielsen et al., 2001; Wu and Kwok, 1997; Ballestra, 2018). Moreover, the degeneracy that occurs in the front-fixing method of Wu and Kwok (1998) was further pointed out in the work of Kim et al. (2013). The finite element method has also been implemented for solving American options based on the front-fixing approach (Holmes and Yang, 2008; Zhang et al., 2014; Song et al. 2017). Holmes and Yang (2008) implemented the Crank-Nicholson method and Zhang et al. (2014) and Song et al. (2017) used the first-order backward method based on the temporal discretization. Over the decades, embedded high order Runge-Kutta adaptive methods have been developed (Fehlberg, 1969; Verner, 1978, Dormand and Prince, 1980; Cash and Karp, 1990; Papageorgiou and Tsitouras, 1996; Tsitouras, 1998; Macdougall and Verner, 2002; Verner, 2010) and implemented in several works of literature (Simos, 1993; Simos and Papakaliatakis, 1998, Burden et al., 2016) for solving heat problems, stochastic wave and Schrodinger equations (Wilkie and Cetinbas 2005; Treblay and Carrighton, 2004), harmonic oscillator problem (Hoover et al. 2016), thin-film model (Romeo et al. 2008) and Lotka Voltera prey-predator model (Paul et al., 2016). The adaptive Runge-Kutta method, which is more effective than the classical fourth Runge-Kutta method, is based on the control of error estimation which results in time step adjustment and optimal selection of time step at each time level. The variation in time step by optimal selection of the time step at each time level provides some computational benefits. In terms of cost, it enables the selection of large time steps in a local region with sufficient smoothness and small-time steps in a local region with large variation and discontinuity (William and Saul, 1992). The motivation of this work is to implement an adaptive and explicit fourth-order Runge-Kutta-Fehlberg time integration method coupled with a fourth-order compact scheme for solving the American put options problem based on the front-fixing approach. To the best of our knowledge, we are the first to implement this combination for solving the American option problem. In short, by implementing a method of extrapolation, we improve on the techniques of Kim et al. (2013) which avoid the degeneracy that occurs near the optimal exercise boundary to obtain the velocity of the latter analytically with high order accuracy in space. We then apply adaptive fourth-order Runge-Kutta-Fehlberg methods to compute the optimal exercise boundary with high order accuracy in time. Coupled with a compact finite difference scheme for spatial discretization, more accurate numerical solutions of the option price, and option Greeks are obtained with fast computation. The rest of the paper is organized as follows. In section 2, we discuss the various transformation involved in our method. In section 3, we employ a compact scheme in the spatial discretization and adaptive Runge-Kutta-Fehlberg method for temporal discretization. In section 4, we investigate and compare the numerical performance of an adaptive Runge-Kutta-Fehlberg with the classical Runge-Kutta and other existing methods and conclude the paper in section 5. Transformations and Free Boundary Analysis
Front-Fixing Logarithmic Transformation
Here, we first employ logarithmic transformation (Egorova et. al., 2016; Wu and Kwok, 1997) to fix the free boundary by using the following relation ๐ฅ = ln ๐๐ ๐ (๐) = ln ๐ โ ln ๐ ๐ (๐) , ๐(๐ฅ, ๐) = ๐(๐, ๐). (2) Applying it to equation (1), we then have ๐๐(๐ฅ, ๐)๐๐ โ 12 ๐ ๐ ๐(๐ฅ, ๐)๐๐ฅ โ (๐ โฒ๐ ๐ ๐ + ๐ โ ๐ ๐(๐ฅ, ๐) = ๐พ โ ๐ = ๐พ โ ๐ ๐ (๐)๐ ๐ฅ , for ๐ฅ < 0; (3๐) where the initial condition (1) is changed to ๐(๐ฅ, 0) = max(๐พ โ ๐พ๐ ๐ฅ , 0) = 0, ๐ฅ โฅ 0, ๐ ๐ (0) = ๐พ. (3๐) By letting ๐ฅ โ 0 โ , we obtain from (3c) that ๐(0, ๐) = ๐พ โ ๐ ๐ (๐) . Thus, together with (1d), we obtain the boundary condition for (3a) as
๐(0, ๐) = ๐พ โ ๐ ๐ (๐), ๐(โ, ๐) = 0. (3๐) Taking further derivative to remove the first-order derivative in (3a), we obtain a system of coupled partial differential equations consisting of the asset, delta, gamma, and speed options as follows: ๐๐(๐ฅ, ๐)๐๐ โ 12 ๐ ๐ ๐(๐ฅ, ๐)๐๐ฅ โ (๐ โฒ ๐ (๐)๐ ๐ (๐) + ๐ โ ๐ ๐๐(๐ฅ, ๐)๐๐ โ 12 ๐ ๐ ๐(๐ฅ, ๐)๐๐ฅ โ (๐ โฒ ๐ (๐)๐ ๐ (๐) + ๐ โ ๐ ๐(๐ฅ, ๐)๐๐ฅ + ๐๐(๐ฅ, ๐) = 0, (4๐) ๐๐(๐ฅ, ๐) ๐๐ โ 12 ๐ ๐ ๐(๐ฅ, ๐)๐๐ฅ โ (๐ โฒ ๐ (๐)๐ ๐ (๐) + ๐ โ ๐ ๐(๐ฅ, ๐)๐๐ฅ + ๐๐(๐ฅ, ๐) = 0, (4๐) ๐๐(๐ฅ, ๐)๐๐ โ 12 ๐ ๐ ๐(๐ฅ, ๐)๐๐ฅ โ (๐ โฒ ๐ (๐)๐ ๐ (๐) + ๐ โ ๐ ๐(๐ฅ, ๐)๐๐ฅ + ๐๐(๐ฅ, ๐) = 0, (4๐) where ๐ฅ > 0, and the initial and boundary conditions for ๐(๐ฅ, ๐) , ๐(๐ฅ, ๐)
๐(๐ฅ, ๐), and
๐(๐ฅ, ๐) are given as : ๐(๐ฅ, 0) = 0, ๐(๐ฅ, 0) = 0, ๐(๐ฅ, 0) = 0, ๐(๐ฅ, 0) = 0, ๐ ๐ (0) = ๐พ, ๐ฅ โฅ 0; (4๐) ๐(0, ๐) = ๐พ โ ๐ ๐ (๐), ๐(0, ๐) = ๐(0, ๐) = ๐(0, ๐) = โ๐ ๐ (๐); (4๐) ๐(โ, ๐) = 0, ๐(โ, ๐) = 0, ๐(โ, ๐) = 0, ๐(โ, ๐) = 0. (4๐)
Transformed Function with the Fixed Free Boundary
Due to degeneracy that occurs near the optimal exercise boundary, we adopt the idea in the work of Kim et al. (2013, 2017) that implements an intermediate (square root) function to avoid such degeneracy. The transformed function is of the form
๐(๐ฅ, ๐) = โ๐(๐ฅ, ๐) โ ๐พ + ๐ ๐ฅ ๐ ๐ (๐), ๐(๐ฅ, ๐) = ๐ (๐ฅ, ๐) + ๐พ โ ๐ ๐ฅ ๐ ๐ (๐), (5๐) with ๐(๐ฅ, ๐) {= 0, ๐ฅ โ [ln ๐ ๐ (โ) โ ln ๐ ๐ (0)],> 0, ๐ฅ โ (0, โ). (5๐) Here , ๐ ๐ (โ) is the asymptotically optimal exercise boundary given as follows: ๐ ๐ (โ) = ๐พ๐พ + 1 ๐พ, ๐พ = 2๐๐ . (5๐) By computing the higher derivatives of
๐(๐ฅ, ๐) and
๐(๐ฅ, ๐) when ๐ฅ = 0 using (1) and (5), Kim et al. (2013) obtained the derivative of the optimal exercise boundary by taking the Taylor expansion of
๐(๐ฅ, ๐) near the optimal exercise boundary with up to fourth-order accuracy as follows:
๐(๐ฅฬ , ๐) = ๐(0, ๐) + ๐ฅฬ ๐ โฒ (0, ๐) + ๐ฅฬ โฒโฒ (0, ๐) + ๐ฅฬ โฒโฒโฒ (0, ๐) + ๐(๐ฅฬ ). (6) Here,
๐(0, ๐) = 0, ๐ โฒ (0, ๐) = โ๐๐๐ , ๐ โฒโฒ (0, ๐) = โ 2๐ ๐ โ๐๐3๐ ; (7๐) ๐ โฒโฒโฒ (0, ๐) = 2๐ ๐2 โ๐๐3๐ + ๐โ๐๐2๐ . (7๐) Note that ๐ ๐ = (๐ + 1๐ ๐ ๐๐ ๐ ๐๐ โ ๐ and ๐ฅฬ โช ๐ฅ is arbitrary and very close to the optimal exercise boundary. Substituting (7a)-(7c) into (6), we obtain a quadratic equation with respect to ๐๐ ๐ / ๐๐. Thus, the velocity of the optimal exercise boundary is presented in quadratic form as follows: ๐๐ ๐ ๐๐ = โ๐ โ โ๐ โ 4๐๐2๐ , (8) where ๐( ๐ฅ ฬ , ๐, ๐ ๐ ) = ๐ ๐ฅ ฬ ๐ ๐2 , ๐ ( ๐ฅ ฬ , ๐, ๐ ๐ ) = โ ๐ ๐ฅ ฬ ๐ ๐ + 2 ๐๐ ๐ฅ ฬ ๐ ๐ ; (9๐) ๐(๐ฅฬ , ๐, ๐ ๐ ) = โ๐ + ๐ ๐ฅฬ ๐ โ ๐๐ ๐ฅฬ + ๐๐ ๐ฅฬ + ๐๐ ๐ฅฬ , ๐ = โ ๐๐พ, ๐ = ๐ โ ๐ . (9๐) Because of the ๐ฅฬ associated with ๐( ๐ฅ ฬ , ๐, ๐ ๐ ), we observed that the order of accuracy in space might not be at least, ๐(๐ฅฬ ) at the boundary when computing the optimal exercise boundary from its derivative. It is well known that for a stable scheme, a third-order boundary condition is consistent with a fourth-order interior scheme (Adam, 1977). To achieve high order accuracy in space, we present an improvement to (8) through a method of extrapolation to increase the truncation error of ๐(๐ฅฬ , ๐) , up to sixth order accuracy by considering the following lemma:
Lemma . Assume
๐(๐ฅฬ , ๐) โ ๐ถ ๐+3 [0, ๐ฟ] , then it holds ๐ ๐(0, ๐) + ๐ ๐(๐ฅฬ , ๐) + โฏ + ๐ ๐ ๐(๐๐ฅฬ , ๐) = ๐ ๐ฅฬ ๐ โฒ (0, ๐) + ๐ ๐ฅฬ ๐ โฒโฒ (0, ๐) + ๐ ๐ฅฬ ๐ โฒโฒโฒ (0, ๐) + ๐(๐ฅฬ ๐+3 ), (10๐) which gives โฒ (0, ๐) + 99๐ฅฬ โฒโฒ (0, ๐) + 9๐ฅฬ (3) (0, ๐) + ๐(๐ฅฬ ). (10๐) Proof . Applying the Taylor expansion at , we obtain ๐(๐ฅฬ , ๐) = ๐(0, ๐) + ๐ฅฬ ๐ โฒ (0, ๐) + ๐ฅฬ โฒโฒ (0, ๐) + ๐ฅฬ (3) (0, ๐) + ๐ฅฬ
24 ๐ (4) (0, ๐) + ๐ฅฬ
120 ๐ (5) (0, ๐)+ ๐(๐ฅฬ ), (11๐) ๐(2๐ฅฬ , ๐) = ๐(0, ๐) + 2๐ฅฬ ๐ โฒ (0, ๐) + 2๐ฅฬ ๐ โฒโฒ (0, ๐) + 4๐ฅฬ (3) (0, ๐) + 2๐ฅฬ (4) (0, ๐) + 4๐ฅฬ
15 ๐ (5) (0, ๐)+ ๐(๐ฅฬ ), (11๐) ๐(3๐ฅฬ , ๐) = ๐(0, ๐) + 3๐ฅฬ ๐ โฒ (0, ๐) + 9๐ฅฬ โฒโฒ (0, ๐) + 9๐ฅฬ (3) (0, ๐) + 27๐ฅฬ (4) (0, ๐) + 81๐ฅฬ
40 ๐ (5) (0, ๐)+ ๐(๐ฅฬ ). (11๐) Multiplying (11a) by 16 and subtracting from (11b), we obtain โฒ (0, ๐) + 6๐ฅฬ ๐ โฒโฒ (0, ๐) + 4๐ฅฬ (3) (0, ๐) โ 2๐ฅฬ
15 ๐ (5) (0, ๐)+ ๐(๐ฅฬ ). (11๐) Multiplying (11b) by 81/16 and subtracting from (11c), we obtain โฒ (0, ๐) + 45๐ฅฬ โฒโฒ (0, ๐) + 9๐ฅฬ (3) (0, ๐) โ 27๐ฅฬ
40 ๐ (5) (0, ๐)+ ๐(๐ฅฬ ). (11๐) Multiplying (11d) by 81/16 and subtracting from (11e), we obtain (10b). Substituting (7) into (10b), we have ๐ โ๐๐๐ฅฬ + 3๐ ๐2 โ๐๐๐ฅฬ ๐ + 9๐โ๐๐๐ฅฬ + ๐(๐ฅฬ ). (11๐) The new quadratic form of the derivative of the optimal exercise boundary is then given from (11f) as follows: (๐๐ ๐ ๐๐ ) + ๐ ๐๐ ๐ ๐๐ + ๐ = 0, (12๐) where ๐๐ ๐ ๐๐ = โ๐ โ โ๐ โ 4๐2 , ๐(๐ฅฬ , ๐, ๐ ๐ ) = (2 โ 11๐ ๐ , (12๐) ๐(๐ฅฬ , ๐, ๐ ๐ ) = 3๐ ๐ ๐2 ๐ฅฬ [โ81๐(๐ฅฬ , ๐) + 818 ๐(2๐ฅฬ , ๐) โ ๐(3๐ฅฬ , ๐)]+ 3๐ ๐ ๐2 ๐ฅฬ ( ๐ฅฬ 4๐ โ ๐ฅฬ + ๐ฅฬ ๐ + ๐ฅฬ ). (12๐) With (12), we can ensure at least, a third-order approximation in space at the boundary. To approximate the optimal exercise boundary with high order accuracy in time, we then implement an adaptive fourth-order Runge-Kutta-Fehlberg method for temporal discretization which is detailed in the following section. Numerical method
We solve the discretized system of PDEs that consists of the asset, delta, gamma, and speed options in a uniform space grid and non-uniform adaptive time grid [0, โ) ร [0 ๐] . We replace the infinite domain with an estimated boundary ๐ฅ ๐ (Egorova et al., 2016; Kangro and Nicolaides, 2000; Toivanen, 2010,) . Let ๐ represent the node points in the grid and M represent the numbers of grid points, respectively, then we have ๐ฅ ๐ = ๐โ, โ = ๐ฅ ๐ ๐ , ๐ โ [0, M], (13)
Here, the numerical solutions of the asset options, option Greeks, and optimal exercise boundary are represented as (๐ข) ๐๐ , (๐ค) ๐๐ , (๐ฆ) ๐๐ , (๐ง) ๐๐ , and ๐ ๐๐ . Fourth-Order Compact Finite Difference Scheme
We employ a compact finite difference scheme (Zhang and Wang, 2012; Bhatt and Khaliq, 2015) for the spatial discretization of our model. For the interior points, we use the compact scheme discretization as follows: ๐ โฒโฒ (๐ฅ ๐โ1 ) + 10๐ โฒโฒ (๐ฅ ๐ ) + ๐ โฒโฒ (๐ฅ ๐+1 ) = 12โ [๐(๐ฅ ๐โ1 ) โ 2๐(๐ฅ ๐ ) + ๐(๐ฅ ๐+1 )] + ๐(โ ). (14๐) For ๐ = 1 and ๐ = ๐ โ 1, we employ a one-sided formula as follows: โฒโฒ (๐ฅ ) โ 5๐ โฒโฒ (๐ฅ ) + 4๐ โฒโฒ (๐ฅ ) โ ๐ โฒโฒ (๐ฅ ) = 12โ [๐(๐ฅ ) โ 2๐(๐ฅ ) + ๐(๐ฅ )] + ๐(โ ). (14๐) โฒโฒ (๐ฅ ๐โ1 ) โ 5๐ โฒโฒ (๐ฅ ๐โ2 ) + 4๐ โฒโฒ (๐ฅ ๐โ3 ) โ ๐ โฒโฒ (๐ฅ ๐โ4 )= 12โ [๐(๐ฅ ๐โ2 ) โ 2๐(๐ฅ
๐โ1 ) + ๐(๐ฅ ๐ )] + ๐(โ ). (14๐) The matrix-vector form is as follows:
๐ด = 12โ [ โ2 1 0 โฏ 01 โ 2 1 โฎ 1 โ 2 1 1 โ 2 1 0 โฑ โฑ โฑ 0 1 โ 2 1 โฎ 1 โ 2 10 โฏ 0 1 โ 2] ๐โ1ร๐โ1 โฒ ๐ต = [ 14 โ 5 4 โ 1 0 โฏ 0 1 10 1 โฎ 1 10 1 1 10 1 0 โฑ โฑ โฑ 0 โฎ 1 10 1 0 โฏ 0 โ 1 4 โ 5 14 ]
๐โ1ร๐โ1โฒ , ๐ ๐ข = 12โ [ ๐ข ๐ = 0] ๐โ1ร1 ; ๐ ๐ค = 12โ [ ๐ค ๐ = 0 ] ๐โ1ร1 , ๐ ๐ฆ = 12โ [ ๐ฆ ๐ = 0 ] ๐โ1ร1 , ๐ ๐ง = 12โ [ ๐ง ๐ = 0 ] ๐โ1ร1 . (14๐)
Hence, ๐ โฒโฒ = ๐ต โ๐ (๐ด๐ + ๐ ๐ข ), ๐ โฒโฒ = ๐ต โ๐ (๐ด๐ + ๐ ๐ค ), (15๐) ๐ โฒโฒ = ๐ต โ๐ (๐ด๐ + ๐ ๐ฆ ), ๐ โฒโฒ = ๐ต โ๐ (๐ด๐ + ๐ ๐ง ). (15๐) Substituting (14) into (4), we recast our partial differential equations in the form of a system of ordinary differential equations as follows: ๐๐๐๐ = ๐(๐, ๐), ๐๐๐๐ = ๐(๐, ๐), ๐๐๐๐ = ๐(๐, ๐), ๐๐๐๐ = ๐(๐, ๐), (16) where ๐(๐, ๐) = ๐ โ๐ (๐ด๐ + ๐ ๐ข ) + ๐ ๐ ๐ โ ๐๐, (17๐) ๐(๐, ๐) = ๐ โ๐ (๐ด๐ + ๐ ๐ค ) + ๐ ๐ ๐ต โ๐ (๐ด๐ + ๐ ๐ข ) โ ๐๐, (17๐) ๐(๐, ๐) = ๐ โ๐ (๐ด๐ + ๐ ๐ฆ ) + ๐ ๐ ๐ต โ๐ (๐ด๐ + ๐ ๐ค ) + โ๐๐, (17๐) ๐(๐, ๐) = ๐ โ๐ (๐ด๐ + ๐ ๐ง ) + ๐ ๐ ๐ต โ๐ (๐ด๐ + ๐ ๐ฆ ) + โ๐๐. (17๐) We would like to point out some flexibility in this work based on the explicit approach. The rate of change of the optimal exercise boundary is independent of the higher derivatives (delta, gamma, and speed options). By computing the optimal exercise boundary first, we could implement a Dirichlet boundary condition based on (4f). Moreover, the numerical solutions of the asset and delta options with optimal exercise boundary as a coupled system are independent of the higher derivatives (gamma and speed option). The choice of including the gamma and speed options in the coupled system is to approximate them with high order accuracy. Furthermore, it is important to mention that if the choice is to obtain the numerical solutions of the asset option and optimal exercise boundary only, we can further introduce a compact discretization of the first derivative to accommodate such possibility as follows: ๐ โฒ (๐ฅ ๐โ1 ) + 4๐ โฒ (๐ฅ ๐ ) + ๐ โฒ (๐ฅ ๐+1 ) = 3โ [๐(๐ฅ ๐+1 ) โ ๐(๐ฅ ๐โ1 )] + ๐(โ ). (18๐) For ๐ = 1 and ๐ = ๐ โ 1, we employ a one-sided formula as follows (Zhang and Wang, 2012; Bhatt and Khaliq, 2015): โฒ (๐ฅ ) + ๐ โฒ (๐ฅ ) = 1โ [โ 1112 ๐(๐ฅ ) โ 4๐(๐ฅ ) + 6๐(๐ฅ ) โ 43 ๐(๐ฅ ) + 14 ๐(๐ฅ )] + ๐(โ ). (18๐) โฒ (๐ฅ ๐โ1 ) + ๐ โฒ (๐ฅ ๐โ2 )= 1โ [1112 ๐(๐ฅ ๐ ) โ 4๐(๐ฅ ๐โ1 ) + 6๐(๐ฅ
๐โ2 ) โ 43 ๐(๐ฅ
๐โ3 ) โ 14 ๐(๐ฅ
๐โ4 )] + ๐(โ ). (18๐) The matrix-vector form is as follows:
๐ท = [ 4 1 0 โฏ 0 1 4 1 โฎ 1 4 1 1 4 1 0 โฑ โฑ โฑ 0 โฎ 1 4 10 โฏ 0 1 4 ]
๐โ1ร๐โ1 , ๐ฬ ๐ข = 1112โ [ โ๐ข ๐ = 0] ๐โ1ร1 , (18๐) ๐ถ = 3โ [ โ 43 2 โ 49 112 0 โฏ 0โ1 0 1 โฎ โ1 0 1 โ1 0 1 0 โฑ โฑ โฑ 0 โ1 0 1 โฎ โ1 0 10 โฏ โ 112 49 โ 2 43]
๐โ1ร๐โ1 โฒ (18๐) where ๐ โฒ = ๐ท โ๐ (๐ถ๐ + ๐ฬ ๐ข ). Implementing it in (3), we then have ๐๐๐๐ = ๐ โ๐ (๐ด๐ + ๐ ๐ข ) + ๐ ๐ ๐ท โ๐ (๐ถ๐ + ๐ฬ ๐ข ) โ ๐๐, (18๐) which can be used to obtain the numerical solutions of the optimal exercise boundary and the asset option. Adaptive and Classical Fourth-Order Time Integrators
Adaptive Runge-Kutta-Fehlberg Method:
By recasting our system of discretized partial differential equations in the form of ordinary differential equations, we then present (16)-(17) in explicit form as follows: ๐๐ ๐ ๐๐ = ๐ โ๐ (๐ด๐ ๐ + ๐ ๐ข๐ ) + ๐ ๐ ๐ ๐ โ ๐๐ ๐ , (19๐) ๐๐ ๐ ๐๐ = ๐ โ๐ (๐ด๐ ๐ + ๐ ๐ค๐ ) + ๐ ๐ ๐ต โ๐ (๐ด๐ ๐ + ๐ ๐ข๐ ) โ ๐๐ ๐ , (19๐) ๐๐ ๐ ๐๐ = ๐ โ๐ (๐ด๐ ๐ + ๐ ๐ฆ๐ ) + ๐ ๐ ๐ต โ๐ (๐ด๐ ๐ + ๐ ๐ค๐ ) โ ๐๐ ๐ , (19๐) ๐๐ ๐ ๐๐ = ๐ โ๐ (๐ด๐ ๐ + ๐ ๐ง๐ ) + ๐ ๐ ๐ต โ๐ (๐ด๐ ๐ + ๐ ๐ฆ๐ ) โ ๐๐ ๐ . (19๐) In this work, we implement a fourth-order adaptive Runge-Kutta-Fehlberg method (Fehlberg, 1969) based on the coefficients of Cash and Karp (1990). Runge-Kutta-Fehlberg method uses a fifth-order Runge-Kutta method to estimate the local truncation error of the fourth-order Runge-Kutta method (Burden et al., 2016). With a given tolerance, the optimal time step is obtained for each time level. For brevity, we only describe function which follows from (19) for computing the new values and error of the asset option from the RKF as follows: The fourth-order and fifth-order Runge-Kutta methods are given as 1 ๐ ๐+1 = ๐ ๐ + ( 37378 ๐ณ ๐ข1 + 250621 ๐ณ ๐ข3 + 125594 ๐ณ ๐ข4 + 5121771 ๐ณ ๐ข6 ), (20๐) ๐ฬ ๐+1 = ๐ ๐ + ( 282527648 ๐ณ ๐ข1 + 1857548384 ๐ณ ๐ข3 + 1352555296 ๐ณ ๐ข4 + 27714336 ๐ณ ๐ข5 + 14 ๐ณ ๐ข6 ), (20๐) respectively, and the error estimated as ๐ ๐ข = โ๐ฬ ๐+1 โ ๐ ๐+1 โ โ < ๐, (20๐) where ๐ณ ๐ข1 = ๐(๐ ๐ , ๐ ๐ )๐, ๐ณ ๐ข2 = ๐ (๐ ๐ + 15 ๐ณ ๐ข1 , ๐ ๐ + 15 ๐ณ ๐ค1 ) ๐; (20๐) ๐ณ ๐ข3 = ๐ (๐ ๐ + 340 ๐ณ ๐ข1 + 940 ๐ณ ๐ข2 , ๐ ๐ + 340 ๐ณ ๐ค1 + 940 ๐ณ ๐ค2 ) ๐, (20๐) ๐ณ ๐ข4 = ๐ (๐ ๐ + 310 ๐ณ ๐ข1 โ 910 ๐ณ ๐ข2 + 65 ๐ณ ๐ข3 , ๐ ๐ + 310 ๐ณ ๐ค1 โ 910 ๐ณ ๐ค2 + 65 ๐ณ ๐ค3 ) ๐, (20๐) ๐ณ ๐ข5 = ๐ (๐ ๐ โ 1154 ๐ณ ๐ข1 + 52 ๐ณ ๐ข2 โ 7027 ๐ณ ๐ข3 + 3527 ๐ณ ๐ข4 , ๐ ๐ โ 1154 ๐ณ ๐ค1 + 52 ๐ณ ๐ค2 โ 7027 ๐ณ ๐ค3 + 3527 ๐ณ ๐ค4 ) ๐, (20๐) ๐ณ ๐ข6 = ๐ (๐ ๐ + 163155296 ๐ณ ๐ข1 + 175512 ๐ณ ๐ข2 + 57513824 ๐ณ ๐ข3 + 44275110592 ๐ณ ๐ข4 + 2534096 ๐ณ ๐ข5 , ๐ ๐ + 163155296 ๐ณ ๐ค1 + 175512 ๐ณ ๐ค2 + 57513824 ๐ณ ๐ค3 + 44275110592 ๐ณ ๐ค4 + 2534096 ๐ณ ๐ค5 ) ๐. (20โ) Similarly, the mathematical formulation in (20) also follows for computing the option Greeks. Hence, for the sake of brevity, we skip them. Here, ๐ represents the time step . If the condition in (20c) fails based on an arbitrary ๐ , an optimal parameter is determined, from which a new time step is calculated until an optimal time step that satisfied (20c) is obtained. Moreover, if the condition in (20c) is satisfied, a new time step is also estimated which will be used in the next time level. The calculation is done (William and Saul, 1992; Clayton et al., 2019) as follows: ๐ ๐๐๐ค = {0.9๐ ๐๐๐ (๐๐๐/ ๐ ๐ข ) , ๐ โค ๐ ๐ข , ๐๐๐ (๐๐๐/ ๐ ๐ข ) , ๐ > ๐ ๐ข . (21) Classical Runge-Kutta Method:
Here, as a comparison and to calculate the convergent rate with a constant time step, we employ a fourth-order explicit Runge-Kutta method (RK4) for temporal discretization. We fully describe the procedure for solving (19) using the Runge-Kutta method as follows: ๐น ๐ข1 = ๐(๐ ๐ , ๐ ๐ )๐, ๐น ๐ค1 = ๐(๐ ๐ , ๐ ๐ )๐; (22๐) ๐น ๐ฆ1 = ๐(๐ ๐ , ๐ ๐ )๐, ๐น ๐ง1 = ๐(๐ ๐ , ๐ ๐ )๐; (22๐) ๐น ๐ข2 = ๐ (๐ ๐ + 12 ๐น ๐ข1 , ๐ ๐ + 12 ๐น ๐ค1 ) ๐, ๐น ๐ค2 = ๐ (๐ ๐ + 12 ๐น ๐ค1 , ๐ ๐ + 12 ๐น ๐ข1 ) ๐; (22๐) ๐น ๐ฆ2 = ๐ (๐ ๐ + 12 ๐น ๐ฆ1 , ๐ ๐ + 12 ๐น ๐ค1 ) ๐, ๐น ๐ง2 = ๐ (๐ ๐ + 12 ๐น ๐ง1 , ๐ ๐ + 12 ๐น ๐ฆ1 ) ๐; (22๐) ๐น ๐ข3 = ๐ (๐ ๐ + 12 ๐น ๐ข2 , ๐ ๐ + 12 ๐น ๐ค2 ) ๐, ๐น ๐ค3 = ๐ (๐ ๐ + 12 ๐น ๐ค2 , ๐ ๐ + 12 ๐น ๐ข2 ) ๐; (22๐) ๐น ๐ฆ3 = ๐ (๐ ๐ + 12 ๐น ๐ฆ2 , ๐ ๐ + 12 ๐น ๐ค2 ) ๐, ๐น ๐ง3 = ๐ (๐ ๐ + 12 ๐น ๐ง2 , ๐ ๐ + 12 ๐น ๐ฆ2 ) ๐; (22๐) ๐น ๐ข4 = ๐(๐ ๐ + ๐น ๐ข3 , ๐ ๐ + ๐น ๐ค3 )๐, ๐น ๐ค4 = ๐(๐ ๐ + ๐น ๐ค3 , ๐ ๐ + ๐น ๐ข3 )๐; (22๐) ๐น ๐ฆ4 = ๐(๐ ๐ + ๐น ๐ฆ1 , ๐ ๐ + ๐น ๐ค1 )๐, ๐น ๐ง4 = ๐(๐ ๐ + ๐น ๐ง1 , ๐ ๐ + ๐น ๐ฆ1 )๐; (22โ) ๐ ๐+1 = ๐ ๐ + ๐6 (๐น ๐ข1 + 2๐น ๐ข2 + 2๐น ๐ข3 + ๐น ๐ข4 ), ๐ ๐+1 = ๐ ๐ + ๐6 (๐น ๐ค1 + 2๐น ๐ค2 + 2๐น ๐ค3 + ๐น ๐ค4 ); (23๐) ๐ ๐+1 = ๐ ๐ + ๐6 (๐น ๐ฆ1 + 2๐น ๐ฆ2 + 2๐น ๐ฆ3 + ๐น ๐ฆ4 ), ๐ ๐+1 = ๐ ๐ + ๐6 (๐น ๐ง1 + 2๐น ๐ง2 + 2๐น ๐ง3 + ๐น ๐ง4 ). (23๐) Approximation of the Optimal Exercise Boundary : Because of the explicit nature of our proposed method, we need to approximate the optimal exercise boundary before computing the asset option and option Greeks. To achieve this, we discretize (12b) using both adaptive and classical RK4 methods.
Let ๐๐ ๐๐ ๐๐ = ๐(๐ ๐๐ , ๐ข ๐ฅฬ ๐ ) = โ๐ ๐ โ โ(๐ ๐ ) โ 4๐ ๐ with ๐ ๐ = (2 โ 11๐ ๐๐ , (24๐) ๐ ๐ = 3๐ (๐ ๐๐ ) ๐ฅฬ [โ81๐ ๐ฅฬ ๐ + 818 ๐ โ ๐ ]+ 3๐ (๐ ๐๐ ) ๐ฅฬ ( ๐ฅฬ 4๐ โ ๐ฅฬ + ๐ฅฬ ๐ + ๐ฅฬ ). (24๐) For the adaptive Runge-Kutta-Fehlberg method, the fourth order Runge-Kutta method ๐ ๐๐+1 = ๐ ๐๐ + ( 37378 ๐ ๐ ๐ + 250621 ๐ ๐ ๐ + 125594 ๐ ๐ ๐ + 5121771 ๐ ๐ ๐ ), (25๐)
3 is computed simultaneously with the fifth-order Runge-Kutta method ๐ ฬ ๐๐+1 = ๐ ๐๐ + ( 282527648 ๐ ๐ ๐ + 1857548384 ๐ ๐ ๐ + 1352555296 ๐ ๐ ๐ + 27714336 ๐ ๐ ๐ + 14 ๐ ๐ ๐ ), (25๐) and the error estimated as ๐ ๐ ๐ = |๐ ๐๐+1 โ ๐ ฬ ๐๐+1 | < ๐, (25๐) where ๐ ๐ ๐ = ๐(๐ ๐๐ , ๐ข ๐ฅฬ ๐ )๐, ๐ ๐ ๐ = ๐ (๐ ๐๐ + 15 ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐; (25๐) ๐ ๐ ๐ = ๐ (๐ ๐๐ + 340 ๐ ๐ ๐ + 940 ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐, (25๐) ๐ ๐ ๐ = ๐ (๐ ๐๐ + 310 ๐ ๐ ๐ โ 910 ๐ ๐ ๐ + 65 ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐, (25๐) ๐ ๐ ๐ = ๐ (๐ ๐๐ โ 1154 ๐ ๐ ๐ + 52 ๐ ๐ ๐ โ 7027 ๐ ๐ ๐ + 3527 ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐, (25๐) ๐ ๐ ๐ = ๐ (๐ ๐๐ + 163155296 ๐ ๐ ๐ + 175512 ๐ ๐ ๐ + 57513824 ๐ ๐ ๐ + 44275110592 ๐ ๐ ๐ + 2534096 ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐. (25โ) For the classical Runge-Kutta method, we compute as follows: ๐ ๐ ๐ = ๐(๐ ๐๐ , ๐ข ๐ฅฬ ๐ )๐, ๐ ๐ ๐ = ๐ (๐ ๐๐ + 12 ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐, ๐ ๐ ๐ = ๐ (๐ ๐๐ + 12 ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐; (26๐) ๐ ๐ ๐ = ๐ (๐ ๐๐ + ๐ ๐ ๐ , ๐ข ๐ฅฬ ๐ ) ๐, ๐ ๐๐+1 = ๐ ๐๐ + 16 (๐ ๐ ๐ + 2๐ ๐ ๐ + 2๐ ๐ ๐ + ๐ ๐ ๐ ). (26๐) Here, we choose ๐ฅฬ = 2โ in our numerical experiment.
Computational Procedure using Adaptive RKF Time Integrator
In this section, we describe the implementation and algorithm for computing the asset, delta, gamma, and speed options using the adaptive Runge-Kutta-Fehlberg methods based on the coefficients of
Cash and Karp (1990). It is worth noting that in this work, we restrict the error estimate only with the asset option. That is, we use only ๐ ๐ข to confirm to optimal time step. When approximating our numerical solutions, there is a threshold for ๐ above which the optimal exercise boundary, when computed from the quadratic equation, will give a complex value. We first adapt our code to find a maximum ๐ that guarantees a real value for the optimal exercise boundary before 4 proceeding to find and obtain the optimal time step and numerical approximation(s) at each time level, respectively. Algorithms for obtaining the numerical solutions of the optimal exercise boundary, asset option, and the option Greeks using the fourth-order adaptive Runge-Kutta methods are described below. Algorithm.
Algorithm for the Runge-Kutta-Fehlberg method (RKF). initialize ๐ก = 0, โ, ๐, ๐, ๐ด, ๐ต , and ๐๐๐ โณ The initial choice of ๐ is arbitrary and independent of โ initialize ๐ ๐๐ , ๐ ๐ , ๐ ๐ , ๐ ๐ and ๐ ๐ while ๐ก < ๐ if ๐ก + ๐ > ๐ ๐ = ๐ โ ๐ก endif 7. while true compute ๐ ๐๐+1 โณ based on (25) if ๐ ๐๐+1 is a real value, break โณ obtain a maximum ๐ that guarantee real value for ๐ ๐๐+1 else ๐ = ๐๐ โณ if the initial choice of ๐ = โ endif 12. end while 13. compute ๐ ๐ข , ๐ ๐ค , ๐ ๐ฆ , and ๐ ๐ง compute ๐ ๐+1 , ๐ฬ ๐+1 , ๐ ๐+1 , ๐ ๐+1 , and ๐ ๐+1 โณ based on (22) and (23) compute ๐ ๐ข = โ ๐ ฬ ๐+1 โ ๐ ๐+1 โ โ , if ๐ ๐ข < ๐๐๐ set ๐ ๐ = ๐ ๐+1 , ๐ ๐ = ๐ ๐+1 , ๐ ๐ = ๐ ๐+1 , ๐ ๐ = ๐ ๐+1 and ๐ ๐๐ = ๐ ๐๐+1 set ๐ฟ ๐ข = 0.9(๐๐๐/๐ ๐ข ) and ๐ = ๐ฟ ๐ข ๐ โณ based on (24) ๐ก = ๐ก + ๐ else 21. set ๐ฟ ๐ข = 0.9(๐๐๐/๐ ๐ข ) and ๐ = ๐ฟ ๐ข ๐ โณ based on (24) endif 23. repeat Numerical Experiment and Discussion
In this section, the numerical performance of the proposed method is investigated and validated using two examples and further compared with the existing results. The numerical experiment was carried out on the mesh with a uniform grid size.
Example 1.
Consider the example provided in the work of Zhu (2006). The following data are presented
๐พ = 100, ๐ = 1, ๐ = 10%, ๐ = 30%. (27)
In this example, we focus on comparing the values of the optimal exercise boundary. We compared the results of the adaptive Runge-Kutta-Fehlberg method (RKF) coupled with a finite compact scheme (FCS-5 RKF) with that of the method of Zhu (2006), the numerical method of Wu and Kwok (1997), and the classical Runge-Kutta method with a finite compact scheme which we label as FCS-RK4. The results were listed in Table 1. The plots of the asset option, option Greeks, and optimal exercise boundary were displayed in Figs. 1. From Table 1, one can observe that our numerical approximations of the optimal exercise boundary for the FCS-RK4 and FCS-RKF are very close to the analytical approximation of Zhu (2006). Furthermore, we further observe from Fig. 1 that the FCS-RK4 method does not approximate the higher derivative (i.e., speed option) accurately when compared with the FCS-RKF.
Fig. 1a . Asset option, option Greeks and optimal exercise boundary with FCS-RKF ( โ = 0.01, ๐ =๐, ๐ = 10 โ8 ). Fig. 1b . Asset option, option Greeks, and optimal exercise boundary with FCS-RK4 ( โ = 0.01, ๐ = 0.0001 ๐ = ๐ ). Table 1.
Comparison of the optimal exercise boundary for example 1 ( ๐ = ๐, โ = 0.01 ). Optimal Exercise Boundary
Zhu (2007) Wu and Kwok (1997) FCS-RK4 FCS-RKF 76.11
Example 2.
Consider the example in the work of Kim et al. (2013). We compare our result with Kim et al. (2013), the moving boundary method (Muthuraman, 2008), and the Binomial method (Cox et al., 1979) which is used as the benchmark result. Here, we use (18) to compute only the option values and the optimal exercise boundary. The following data are considered 7
๐พ = 100, ๐ = 0.5, ๐ = 5%, ๐ = 20%. (28)
The results were listed in Table 2. In Tables 3 and 4, we present the total CPU time(s), optimal exercise boundary value, and the global minimum and maximum optimal time step of the FCS-RKF based on varying tolerance ๐ and step size. The plot of adaptive optimal time step selection for each time level based on varying tolerance ๐ and step size was further displayed in Fig. 2. From Table 2, one can easily see the better performance of the FCS-RKF. With ๐ = 10 โ8 and from โ =0.0125, the results obtained from the FCS-RKF are very close to the one obtained from the binomial method that serves as a benchmark in this example. From Tables 3 and 4, we observe the dependent of the optimal time selection on the tolerance ๐ and fixed step size in direct proportion. From Fig. 2, we observe the concentration of small oscillation near the payoff ( ๐ = 0 ) as the tolerance and the step size is reduced. This is expected because of the variation due to the discontinuity in the first derivative of the asset option that normally occurs at the payoff. Hence, a small varying time step is needed in the region near the payoff. Furthermore, it is important to observe from Fig. 2 that the adaptive method could be very useful in detecting unknown locations of discontinuity or rapid variation in systems (Gear and รsterby, 1984; Dieci and Lopez, 2012) when it is implemented efficiently. Table 2.
Comparison of the asset option in example 2.
S Binomial MBM Kim et al. FCS-RK4 80 20.0000 20.0000 20.0000 20.0000 90 10.6661 10.6680 10.6661 10.6680 100 4.6556 4.6504 4.6549 4.6588 110 1.6681 1.6629 1.6686 1.6704 120 0.4976 0.4993 0.4985 0.4988 S FCS-RKF ( ๐ = 10 โ8 ) โ =
80 20.0000 20.0000 20.0000 90 10.6653 10.6660 10.6661 100 4.6539 4.6554 4.6555 110 1.6666 1.6678 1.6679 120 0.4968 0.4974 0.4975
Table 3.
Performance of the FCS-RKF based on fixed step size ( โ = 0.01 ) and varying tolerance. S ๐ = 10 โ3 ๐ = 10 โ5 ๐ = 10 โ8 CPU time(s) 23.61 57.57 1667.53 ๐ ๐ (๐) Table 4.
Performance of the FCS-RKF based on fixed tolerance ๐ = 10 โ8 and varying step size. S โ = 0.1 โ = 0.05 โ = 0.01 CPU time(s) 3.93 19.97 1667.53 ๐ ๐ (๐) Fig. 2a.
Optimal time step selection for each time level using FCS-RKF with a fixed โ = 0.01 . Fig. 2b.
Optimal time step selection for each time level using FCS-RKF with a fixed ๐ = 10 โ8 . 9 To check the maximum errors and convergent rates in space, we used (12), example 1, and FCS-RK4 and selected a constant time step ๐ = 1 ร 10 โ6 with varying step size โ = 0.2, 0.1, 0.05, 0.025, 0.0125 and . We then compared the convergence rate in space based on the extrapolation method we implemented in (12) with the method of Kim et al. (2013) in (8). We further check the maximum errors and the convergent rates of the delta option based on the method of extrapolation with a constant time step ๐ =1 ร 10 โ6 and varying step size โ = 0.2, 0.1, 0.05, 0.025 , and . The results are displayed in Table 5. For the asset option, the method of Kim et al. (2013) provides up to second order in space as seen in Table 5. However, by improving it with a method of extrapolation, we obtained a high order accuracy in space for both the asset and delta options which is in close agreement with the theoretical convergent rate. Table 5.
Maximum errors and convergent rates in space of the asset and delta options with ๐ = 1 ร 10 โ6 . Asset option Kim et al. (2013) Our method โ maximum error convergent rate maximum error convergent rate ~ ~ ~ ~ โ1 ~ โ0 ~ โ1 โ0.236 โ2 โ1 โ3 โ2 โ4 โ2 โ5 Delta Option (Our method) maximum error convergent rate ~ ~ โ0 ~ โ1 โ2 โ3 โ5 Conclusion
We have proposed an adaptive and explicit fourth-order Runge-Kutta-Fehlberg method with a fourth-order compact scheme for pricing American options. By implementing logarithmic transformation, taking a further derivative, improving the method of Kim et al. (2013) by adopting a method of extrapolation, we then obtain an analytical formula for the velocity of the optimal exercise boundary with high order accuracy in space. We further recast the free boundary problem to a system of coupled ordinary differential equations, employ compact finite difference for spatial discretization with Dirichlet boundary condition and implement an adaptive fourth-order Runge-Kutta-Fehlberg method for temporal 0 discretization. This enables us to approximate the optimal exercise boundary, options value, and option Greeks in the set of coupled ODEs with high order accuracy both in space and in time. Furthermore, we check the convergent rate of our numerical method with the FCS-RK4 method and confirm that the numerical convergent rate is in close agreement with the theoretical convergent rate. By further comparing the result from the FCS-RKF method with the existing methods including the classical Runge-Kutta (FCS-RK4), we then validate the superiority of the adaptive method.
References Adam, Y. (1977). Highly accurate compact implicit methods and boundary conditions. Journal of Computational Physics, 24, 10-22. 2.
Ballestra, L. V. (2018). Fast and accurate calculation of American option prices. Decisions in Economics and Finance, 41 (399-426). 3.
Bhatt, H. P., and Khaliq, A. Q. M. (2016). Fourth-order compact schemes for the numerical simulation of coupled Burgersโ equation. Computer Physics Communications, 200 (117-138). 4.
Burden, R. L., Faires, D. J., and Burden, A. M. (2016). Numerical Analysis. Cengage Learning, Boston, MA. 5.
Cash R. J., and Karp, A. H. (1990). A variable order Runge-Kutta for initial value problems with rapidly varying right-hand sides. ACM Transaction on Mathematical Software, 16 (201-222). 6.
Clayton, S. L., Lemma, M., and Chowdhury, A. (2019). Numerical solutions of nonlinear ordinary differential equations by using adaptive Runge-Kutta method. Journal of Advances in Mathematics, 16 (147-154). 7.
Company, R., Egorova, V.N., and Jรณdar, L. (2014). A positive, stable, and consistent front-fixing numerical scheme for American options. Springer International Publishing AG 2016. 8.
Company, R., Egorova, V.N., and Jรณdar, L. (2014). Solving American option pricing models by the front fixing method: numerical analysis and computing. Abstract and Applied Analysis. 9.
Cox, J. C., Ross, S. A., and Rubinstein, M. (1979). Option pricing: a simplified approach. Journal of Financial Economics, 7 (229โ263). 10.
Dieci, L. and Lopez, L. (2012). A survey of numerical methods for IVPs of ODEs with discontinuous right-hand side. Journal of Computational and Applied Mathematics, 236 (3967โ3991). 11.
Dormand, J. R., and Prince, J. P. (1980). A family of embedded Rung-Kutta formulae. Journal of Computational and Applied Mathematics, 6. 1 12.
Egorova, V. N., Company, R., and Jรณdar, L. (2016). A new efficient numerical method for solving American option under regime switching model. Computers and Mathematics with Applications, 71 (224โ237). 13.
Fehlberg, E. (1969). Low-order classical Runge-Kutta formulas with step size control and their application to some heat transfer problems. NASA Technical Report 315. 14.
Holmes, A. D., and Yang, H. (2008). A front-fixing finite element method for the valuation of American options. SIAM Journal of Scientific Computing, 30 (2158-2180). 15.
Gear, W. C. and รsterby, O. (1984). Solving ordinary differential equations with discontinuities. ACM Transaction in Mathematical Software, 10 (23-44). 16.
Hoover W. G., Sprot, J. C., and Hoover C. G. (2016). Adaptive Runge-Kutta integration for stiff systems: Comparing Nose and Nose-Hoovers dynamics for the harmonic oscillator. American Journal of Physics, 84. 17.
Kangro, R., and Nicolaides, R. (2000). Far field boundary conditions for Black--Scholes equations. SIAM Journal on Numerical Analysis, 38 (4), 1357โ1368. 18.
Kim, B. J., Ma, Y., and Choe, H. J. (2013). A simple numerical method for pricing an American put option. Journal of Applied Mathematics. 19.
Kim, B. J., Ma, Y., and Choe, H. J. (2017). Optimal exercise boundary via intermediate function with jump risk. Japan Journal of Industrial and Applied Mathematics, 34 (779-792). 20.
Macdougall, T., and Verner, J. H. (2002). Global error estimators for 7, 8 Runge-Kutta pairs. Numerical Algorithm, 31 (215-231). 21.
McKean, JR., H. P. (1965). A free boundary problem for the heat equation arising from a problem in mathematical economics. Industrial Management Review, 6 (32-39). 22.
Muthuraman, K. (2008). A moving boundary approach to American option pricing. Journal of Economics, Dynamics, and Control, 32 (3520โ3537). 23.
Nielsen, B. F., Skavhaug O., and Tveito (2002). A penalty and front-fixing methods for the numerical solution of American option problems. Journal of Computational Finance, 5 (4). 24.
Papageorgiou, G., and Tsitouras, C. (1996). Continuous extensions to high order Runge-Kutta methods. International Journal of Computer Mathematics, 65 (273-291). 25.
Paul, S., Mondal, S. P., and Bhattacharya, P. (2016). Numerical solution of Lotka Volterra prey predator model by using Runge-Kutta-Fehlberg method and Laplace Adomain decomposition method. Alexandria Engineering Journal. 2 26.
Romeo, A., Finocchio, G., Carpentieri, M., Torres. L., Consolo, G., and Azzerboni, B. (2008). A numerical solution of the magnetization reversal modeling in a permalloy thin film using fifth order Runge-Kutta method with adaptive step size control. Physica, 403 (464-468). 27.
Simos, T. E. (1993). A Runge-Kutta Fehlberg method with phase-lag of order infinity for initial-value problems with oscillation solution. Computers and Mathematics with Application, 25 (95-101). 28.
Simos, T. E., and Papakaliatakis, G. (1998). Modified Runge-Kutta Verner methods for the numerical solution of initial and boundary-value problems with engineering application. Applied Mathematical Modelling, 22 (657-670). 29.
Song, H., Zhang, K., and Li, Y. (2017). Finite element and discontinuous Galerkin methods with perfect matched layers for American options. Numerical Mathematics Theory Methods and Application, 10 (829-521). 30.
Toivanen, J. (2010). Finite difference methods for early exercise options. Encyclopedia of Quantitative Finance. 31.
Tremblay, J. C., and Carrington Jr., T. (2004). Using preconditioned adaptive step size Runge-Kutta methods for solving the time-dependent Schrodinger equation. Journal of Chemical Physics, 121. 32.
Tsitouras, C. (1998). A parameter study of explicit Runge-Kutta pairs of orders 6(5). Applied Mathematics Letter, 11 (65-69). 33.
Verner, J. H. (1978). Explicit Runge-Kutta methods with estimates of the local truncation error. Siam Journal of Numerical Analysis, 15 (772โ790). 34.
Verner, J. H. (2010). Numerically optimal Runge-Kutta pairs with interpolants. Numerical Algorithm, 53 (383-396). 35.
William, H. P., and Saul, A. T. (1992). Adaptive stepsize Runge-Kutta Integration. Computer in Physics, 8. 36.
Wilkie, J., and Cetinbas, M. 2005. Variable-stepsize Runge-Kutta for stochastic Schrodinger equations. Physics Letters A, 337 (166-182). 37.
Wu, L., and Kwok, Y.K. (1997). A front-fixing method for the valuation of American options. J. Finance. Eng. 6 (2), 83โ97. 38.
Zhang, P., and Wang, J. (2012). A predictor-corrector compact finite difference scheme for Burgersโ equation. Applied Mathematics and Computation, 219 (892-898). 39.