Next Article in Journal
On the Realization of Exact Upper Bounds of the Best Approximations on the Classes H1,1 by Favard Sums
Next Article in Special Issue
Application of the Triple Sumudu Decomposition Method for Solving 1+1 and 2+1-Dimensional Boussinesq Equations
Previous Article in Journal
Infinite Series Concerning Tails of Riemann Zeta Values
Previous Article in Special Issue
Exact and Approximate Solutions for Linear and Nonlinear Partial Differential Equations via Laplace Residual Power Series Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Efficiency of Some Predictor–Corrector Methods with Fourth-Order Compact Scheme for a System of Free Boundary Options

1
Department of Mathematics and Statistics, University of Calgary, Calgary, AB T2N 1N4, Canada
2
Mathematics and Statistics, Louisiana Tech University, Ruston, LA 71272, USA
*
Author to whom correspondence should be addressed.
Axioms 2023, 12(8), 762; https://doi.org/10.3390/axioms12080762
Submission received: 26 June 2023 / Revised: 17 July 2023 / Accepted: 31 July 2023 / Published: 2 August 2023
(This article belongs to the Special Issue Differential Equations and Related Topics)

Abstract

:
The trade-off between numerical accuracy and computational cost is always an important factor to consider when pricing options numerically, due to the inherent irregularity and existence of non-linearity in many models. In this work, we first present fast and accurate (1,2) and (2,2) predictor–corrector methods with a fourth-order compact finite difference scheme for pricing coupled system of the non-linear free boundary option pricing problem consisting of the option value and delta sensitivity. To predict the optimal exercise boundary, we set up a high-order boundary scheme, which is strategically derived using a combination of the fourth-order Robin boundary scheme and the fourth-order compact finite difference scheme near boundary. Furthermore, we implement a three-step high-order correction scheme for computing interior values of the option value and delta sensitivity. The discrete matrix system of this correction scheme has a tri-diagonal structure and strictly diagonal dominance. This nice feature allows for the implementation of the Thomas algorithm, thereby enabling fast computation. The optimal exercise boundary value is also corrected in each of the three correction steps with the derived Robin boundary scheme. Our implementations are fast on both coarse and very refined grids and provide highly accurate numerical approximations. Moreover, we recover a reasonable convergence rate. Further extensions to high-order predictor two-step corrector schemes are elaborated.

1. Mathematical Model

Predictor–corrector numerical methodologies have been implemented for solving American-style option pricing partial differential equation (PDE) models, as seen in the past literature. Solution frameworks with the predictor–corrector schemes have been implemented for solving these models in the form of variational inequality [1], with a penalty term [2,3] or as a free boundary problem [4,5,6,7]. The first two formulations introduce some constraints that remove the free boundary, known as the optimal exercise boundary. Here instead, because of our interest in computing the optimal exercise boundary, option value, and delta sensitivity with much precision, we consider the American-style options as a free boundary problem. If we consider non-dividend-paying put options V ( t , S ) written on an underlying asset with a price S ( t ) , strike price K, and a time to maturity T, then the free boundary partial differential equation can be written as
V t ( t , S ) + σ 2 2 V S S ( t , S ) + r S V S ( t , S ) r V ( t , S ) = 0 ,
where the asset price S ( t ) is driven by geometric Brownian motion W ( t ) , given as
d S ( t ) = μ S ( t ) d t + σ S ( t ) d W ( t ) ,
where σ represents the volatility and r is the interest rate. Under risk neutral probability Q , we have
d S ( t ) = r S ( t ) d t + σ S ( t ) d W ( t ) .
Let τ = T t ; then, the governing differential equation in (1) can be changed to the free boundary PDE as follows:
V t ( t , S ) σ 2 2 V S S ( t , S ) r S V S ( t , S ) + r V ( t , S ) = 0 , S > s f ( τ ) ;
V ( τ , S ) = K S , V S ( τ , S ) = 1 , S < s f ( τ ) .
The boundary and initial conditions are assumed to be
V ( τ , s f ( τ ) ) = K s f ( τ ) , V S ( τ , s f ( τ ) ) = 1 ;
V ( τ , ) = 0 , V S ( τ , ) = 0 ;
V ( 0 , S ) = ( K S ) + .
To fix the free boundary and to further remove the convective term, we implement the Landau transformation [8] and further take the derivative as follows:
x = ln S s f ( τ ) , V ( τ , S ) = U ( τ , x ) , W ( τ , x ) = U x ( τ , x ) .
The transformation above results in a system of non-linear American option pricing models consisting of option values and delta sensitivity as follows:
U τ ( τ , x ) σ 2 2 U x x ( τ , x ) ξ ( τ ) U x ( τ , x ) + r U ( τ , x ) = 0 , x > 0 ;
W τ ( τ , x ) σ 2 2 W x x ( τ , x ) ξ ( τ ) U x x ( τ , x ) + r W ( τ , x ) = 0 , x > 0 ;
ξ ( τ ) = r + s f ( τ ) s f ( τ ) σ 2 2 ,
U ( τ , x ) = K e x s f ( τ ) , W ( τ , x ) = e x s f ( τ ) , x < 0 ,
U ( 0 , x ) = 0 , U ( τ , 0 ) = K s f ( τ ) , U ( τ , ) = 0 , τ > 0 ;
W ( 0 , x ) = 0 , W ( τ , 0 ) = s f ( τ ) , W ( τ , ) = 0 , τ > 0 .
Note that, for the second derivative of the option value, we have the following implied boundary conditions.
U x x ( τ , 0 ) = 2 r K σ 2 s f ( τ ) , U x x ( τ , ) = 0 .
To the best of our knowledge, only a few authors have solved the American-style options as a free boundary problem with predictor–corrector schemes and their extensions based on a front-fixing approach. Zhu and Zhang [6] used a second-order central difference scheme with the (1,2) Euler-CN scheme to solve the fixed free boundary American option model. Zhu and Chen [7] further extended the method in [6] for solving stochastic volatility with the ADI time integration scheme. Chen et al. [4] solved the American options under the finite moment log-stable model with the predictor–corrector method in [6]. Hajipour and Malek [5] used both BDF3 for the predictor and corrector schemes and recovered the second-order convergence rate result. The above-mentioned methods give at most a second-order convergence rate in space and do not consider a system of pricing consisting of option value and Greeks. In this study, we propose a kind of fast and accurate method of explicit–implicit predictor–corrector schemes for solving systems of American-style option pricing models as a free boundary problem by using fourth-order compact finite difference and the discrete solution of the option value and delta sensitivity. Here, the systems of free boundary PDEs are used to solve the option value, delta sensitivity, and optimal exercise boundary simultaneously. In addition, we further explore and study a suite of high-order predictor–corrector schemes and understand their relative performance for solving the above-mentioned PDEs. To the best of our knowledge, we are the first to propose predictor–corrector methods with a fourth-order compact finite difference scheme for solving systems of fixed free boundary option pricing PDEs that simultaneously approximate the option value, delta sensitivity, and the optimal exercise boundary.
The remaining sections of this work are organized as follows: In Section 2, we describe order (1,2) Euler-CN and (2,2) Leapfrog-CN predictor–corrector methods for solving the system of PDEs in (10) and (11) using a fourth-order compact finite difference scheme. In Section 3, we explore the high-order predictor–corrector schemes. We verify the performance of these schemes with some of the existing methods and present our results and findings in Section 4. We then conclude our study in Section 5.

2. Order (1,2) and (2,2) Predictor–Corrector Compact Differencing

In the section, we present a suite of predictor–corrector fourth-order compact finite difference schemes for solving a system of fixed free boundary option pricing PDEs consisting of the option value and delta sensitivity. Our space-time grid is defined as follows:
x i + 1 x i = h = x M M , τ n + 1 τ n = k = T N , x [ 0 , x M ] , τ [ 0 , T ] .
We denote the numerical approximations of the option value, delta sensitivity, and the optimal exercise boundary at grid points x i and τ n as u ( τ n , x i ) , w ( τ n , x i ) , and s f ( τ n ) , respectively. x M is denoted the far right artificial boundary, which replaces the infinite right boundary. Moreover, the far right boundary value for the option value and delta sensitivity is as follows:
u ( τ n , x M ) = 0 , w ( τ n , x M ) = 0 .

2.1. Order (1,2) Euler-CN PC Method with Fourth-Order Compact Scheme

Here, we set up a high-order scheme for predicting the optimal exercise boundary, which is based on the combination of the fourth-order scheme near boundary and a third-order Robin-style boundary scheme. To this end, we present the two high-order boundary and near-boundary schemes by considering the following Lemma for the third-order Robin boundary scheme.
Lemma 1.
Assume that  u ( τ n , x i ) C 1 , 2 ( ( 0 , T ] , [ 0 , x M ] ) . It holds that
7 2 u ( τ n , x 0 ) + 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) = 3 h u x ( τ n , x 0 ) + h 2 u x x ( τ n , x 0 ) + O ( h 4 ) .
Proof. 
Let x i + j x i = j h . The Taylor series expansion around x 0 = 0 gives
u ( τ n , x 1 ) = u ( τ n , x 0 ) + h u x ( τ n , x 0 ) + h 2 2 u x x ( τ n , x 0 ) + h 3 6 u x x x ( τ n , x 0 ) + O ( h 4 ) ,
u ( τ n , x 2 ) = u ( τ n , x 0 ) + 2 h u x ( τ n , x 0 ) + 4 h 2 2 u x x ( τ n , x 0 ) + 4 h 3 3 u x x x ( τ n , x 0 ) + O ( h 4 ) .
Multiplying (20) by 8 and then subtracting from (21), we obtain (19).    □
It is important to mention, as we will see below, that the manner for which we implement (19) enables us to obtain a fourth-order approximation. To obtain the optimal exercise boundary scheme, we construct the following fourth-order compact scheme near boundary [9,10,11]:
u x x ( τ n , x i ) = 2 h 2 u ( τ n , x i 1 ) 2 u ( τ n , x i ) + u ( τ n , x i + 1 ) 1 2 h 2 u x ( τ n , x i + 1 ) u x ( τ n , x i 1 ) + O ( h 4 ) = U ( τ n , x i ) + O ( h 4 ) , i = 1 , 2 .
If we consider value matching, smooth pasting, and the implied second derivative boundary conditions as follows,
u ( τ n , x 0 ) = K s f ( τ n ) , w ( τ n , x 0 ) = u x ( τ n , x 0 ) = s f ( τ n ) , u x x ( τ n , x 0 ) = 2 r K σ 2 s f ( τ n ) ,
we obtain from Lemma 1 that
4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) = 7 2 + 2 r h 2 σ 2 K 7 2 + 3 h + h 2 s f ( τ n ) + O ( h 4 ) .
Based on the forward Euler scheme and the fixed free boundary American options representing the option value in (10), we obtain
σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) = σ 2 2 4 u x x ( τ n , x 1 ) 1 2 u x x ( τ n , x 2 ) = 4 u τ ( τ n , x 1 ) 1 2 u τ ( τ n , x 2 ) 1 s f ( τ n ) d s f ( τ n ) d τ 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) r σ 2 2 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + r 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) = 4 u ( τ n + 1 , x 1 ) u ( τ n , x 1 ) k 1 2 u ( τ n + 1 , x 2 ) u ( τ n , x 2 ) k 1 s f ( τ n ) s f ( τ n + 1 ) s f ( τ n ) k 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) r σ 2 2 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + r 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) + O ( h 4 + k ) .
Multiplying it by k, we then obtain
σ 2 k 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) = 4 u ( τ n + 1 , x 1 ) u ( τ n , x 1 ) 1 2 u ( τ n + 1 , x 2 ) u ( τ n , x 2 ) s f ( τ n + 1 ) s f ( τ n ) s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) r σ 2 2 k 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + r k 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) + O ( h 4 + k 2 ) .
Note that
4 u ( τ , x 1 ) 1 2 u ( τ , x 2 ) = 7 2 + 2 r h 2 σ 2 K 7 2 + 3 h + h 2 s f ( τ ) .
Denote
a 1 = 7 2 + 2 r h 2 σ 2 , a 2 = 7 2 + 3 h + h 2 .
We may simplify (26) to
σ 2 k 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) = a 1 K a 2 s f ( τ n + 1 ) + ( r k 1 ) [ a 1 K a 2 s f ( τ n ) ] s f ( τ n + 1 ) s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) r σ 2 2 k 1 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + O ( h 4 + k 2 ) .
Denote
a 3 ( τ n ) = ( r k 1 ) [ a 1 K a 2 s f ( τ n ) ] , a 4 ( τ n ) = r σ 2 2 k 1 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
The optimal exercise boundary is then predicted from the Euler scheme and fourth-order difference scheme as follows:
s f ( p ) ( τ n + 1 ) = H 1 ( τ n ) s f ( τ n ) + O ( h 4 + k 2 ) ,
H 1 ( τ n ) = σ 2 k 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) a 3 ( τ n ) + a 4 ( τ n ) a 1 K a 2 s f ( τ n ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
Remark 1.
At this point, it is worth acknowledging the work of Egorova et al. [12], where they formulated a similar equation in (32) with the second-order central finite difference scheme to solve the American option pricing model with the second-order explicit method. They further implemented another similar approach for solving the regime-switching American option pricing problem [13]. Here, our formulation ensures that the present scheme admits fourth-order accuracy, and we only use (32) to predict the optimal exercise boundary and then use (27) to correct the optimal exercise after each implicit iteration step.
Furthermore, the time-dependent coefficient associated with the transformed American option pricing model is then predicted as follows:
ξ p ( τ n + 1 2 ) = r + 2 s f ( p ) ( τ n + 1 ) s f ( τ n ) k s f ( p ) ( τ n + 1 ) + s f ( τ n ) σ 2 2 ,
which is a part of the discrete system when correcting the option value, delta sensitivity, and optimal exercise boundary. It is worth mentioning that because we present two systems of fixed free boundary PDEs consisting of the option value and its first derivative known as the delta sensitivity, we can easily use the discrete solution of the delta sensitivity for the prediction of the optimal exercise boundary. Note that the left boundary value of the option value, which is used to formulate the first step of the correction scheme, can be obtained from the predicted value of the optimal exercise boundary as follows:
u ( τ n + 1 , x 0 ) = K s f ( p ) ( τ n + 1 ) .
Below, we present the well-known fourth-order compact finite difference scheme for the interior approximations of the option value and delta sensitivity with the implicit Crank–Nicholson (CN) time integration scheme:
u x x ( τ n , x i 1 ) + 10 u x x ( τ n , x i ) + u x x ( τ n , x i + 1 ) = 12 h 2 u ( τ n , x i 1 ) 2 u ( τ n , x i ) + u ( τ n , x i + 1 ) + O ( h 4 ) , i = 2 , 3 M 1 ,
Implementing this high-order interior scheme for the two systems of transformed and nonlinear PDEs in (10) and (11), we obtain a scheme for the option value as follows:
σ 2 2 u x x ( τ n + 1 2 , x i 1 ) + 10 u x x ( τ n + 1 2 , x i ) + u x x ( τ n + 1 2 , x i + 1 ) = u τ ( τ n + 1 2 , x i 1 ) + 10 u τ ( τ n + 1 2 , x i ) + u τ ( τ n + 1 2 , x i + 1 ) ξ c ( τ n + 1 2 ) w ( τ n + 1 2 , x i 1 ) + 10 w ( τ n + 1 2 , x i ) + w ( τ n + 1 2 , x i + 1 ) + r u ( τ n + 1 2 , x i 1 ) + 10 u ( τ n + 1 2 , x i ) + u ( τ n + 1 2 , x i + 1 ) = 6 σ 2 h 2 u ( τ n + 1 2 , x i 1 ) 2 u ( τ n + 1 2 , x i ) + u ( τ n + 1 2 , x i + 1 ) .
Applying the CN method in time for (36), we then obtain the discrete equation as follows:
u τ ( τ n + 1 , x i 1 ) u τ ( τ n , x i 1 ) 12 k + 10 u τ ( τ n + 1 , x i ) u τ ( τ n , x i ) 12 k + u τ ( τ n + 1 , x i + 1 ) u τ ( τ n , x i + 1 ) 12 k = σ 2 4 h 2 u ( τ n + 1 , x i 1 ) 2 u ( τ n + 1 , x i ) + u ( τ n + 1 , x i + 1 ) + σ 2 4 h 2 u ( τ n , x i 1 ) 2 u ( τ n , x i ) + u ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) 24 [ w ( τ n + 1 , x i 1 ) + w ( τ n , x i 1 ) ] + 10 ξ c ( τ n + 1 2 ) 24 [ w ( τ n + 1 , x i ) + w ( τ n , x i ) ] + ξ c ( τ n + 1 2 ) 24 [ w ( τ n + 1 , x i + 1 ) + w ( τ n , x i + 1 ) ] r 24 [ u ( τ n + 1 , x i 1 ) + u ( τ n , x i 1 ) ] 10 r 24 [ u ( τ n + 1 , x i ) + u ( τ n , x i ) ] r 24 [ u ( τ n + 1 , x i + 1 ) + u ( τ n , x i + 1 ) ] .
Simplifying further, we obtain a discrete equation as follows:
1 12 σ 2 k 4 h 2 + r k 24 u ( τ n + 1 , x i 1 ) + 10 12 + σ 2 k 2 h 2 + 10 r k 24 u ( τ n + 1 , x i ) + 1 12 σ 2 k 4 h 2 + r k 24 u ( τ n + 1 , x i + 1 ) = 1 σ 2 k 4 h 2 r k 24 u ( τ n , x i 1 ) + 1 σ 2 k 2 h 2 10 r k 24 u ( τ n , x i ) + 1 + σ 2 k 4 h 2 r k 24 u ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) k 24 [ w ( τ n + 1 , x i 1 ) + w ( τ n , x i 1 ) ] + 10 ξ c ( τ n + 1 2 ) k 24 [ w ( τ n + 1 , x i ) + w ( τ n , x i ) ] + ξ c ( τ n + 1 2 ) k 24 [ w ( τ n + 1 , x i + 1 ) + w ( τ n , x i + 1 ) ] .
Denote
a i , i = 10 12 + σ 2 k 2 h 2 + 10 r k 24 , a i , i 1 = a i , i + 1 = 1 12 σ 2 k 4 h 2 + r k 24 ;
b i , i = 10 12 σ 2 k 2 h 2 10 r k 24 , b i , i 1 = b i , i + 1 = 1 12 + σ 2 k 4 h 2 r k 24 .
Hence, (38) can be re-written as
a i , i 1 u ( τ n + 1 , x i 1 ) + a i , i u ( τ n + 1 , x i ) + a i , i + 1 u ( τ n + 1 , x i + 1 ) = b i , i 1 u ( τ n , x i 1 ) + b i , i u ( τ n , x i ) + b i , i + 1 u ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) k 24 [ w ( τ n + 1 , x i 1 ) + w ( τ n , x i 1 ) ] + 10 ξ c ( τ n + 1 2 ) k 24 [ w ( τ n + 1 , x i ) + w ( τ n , x i ) ] + ξ c ( τ n + 1 2 ) k 24 [ w ( τ n + 1 , x i + 1 ) + w ( τ n , x i + 1 ) ] .
It can be seen that (41) presents a tridiagonal discrete matrix system. Moreover, it is easy to see that
| a i , i | > | a i , i 1 + a i , i + 1 | .
Thus, our discrete system admits strictly diagonal dominance and can easily be solved using the Gauss–Seidel (GS) iteration method or Thomas algorithm. Moreover, we use the previous solution of w ( τ n + 1 , x i ) when solving (41) using the Thomas algorithm. That is, we use w ( τ n + 1 l , x i ) when solving for u ( τ n + 1 l + 1 , x i ) .
Once u n + 1 is obtained, the optimal exercise boundary is then corrected based on (43):
s f ( c ) ( τ n + 1 l ) = M ( τ n + 1 ) a 1 K a 2 + O ( h 4 ) , l = 1 , 2 , 3 ;
where
M ( τ n + 1 ) = 4 u ( τ n + 1 , x 1 ) 1 2 u ( τ n + 1 , x 2 ) , a 1 = 7 2 + 2 r h 2 σ 2 , a 2 = 7 2 + 3 h + h 2 .
Here, l represent the number of iterations carried out at each time-level. In this research work, we only consider three iterations at each time-level based on our implemented corrector scheme since we have observed that three iterations with our scheme yield a very good approximation. Moreover, the optimal exercise boundary is re-corrected whenever the option’s value is approximated in each iteration step based on (43). For brevity, we remove the l notation in the optimal exercise boundary. The boundary value of the delta sensitivity is then computed from s f ( c ) ( τ n + 1 ) as follows:
u ( τ n + 1 , x 0 ) = K s f ( c ) ( τ n + 1 ) , w ( τ n + 1 , x 0 ) = s f ( c ) ( τ n + 1 ) .
The time-dependent coefficient is corrected as follows:
ξ c ( τ n + 1 2 ) = r + 2 s f ( c ) ( τ n + 1 ) s f ( τ n ) k s f ( c ) ( τ n + 1 ) + s f ( τ n ) σ 2 2
Similarly for the delta sensitivity, we obtain the discrete system as follows:
σ 2 2 w x x ( τ n + 1 2 , x i 1 ) + 10 w x x ( τ n + 1 2 , x i ) + w x x ( τ n + 1 2 , x i + 1 ) = w τ ( τ n + 1 2 , x i 1 ) + 10 w τ ( τ n + 1 2 , x i ) + w τ ( τ n + 1 2 , x i + 1 ) ξ c ( τ n + 1 2 ) h 2 u ( τ n + 1 2 , x i 1 ) 2 u ( τ n + 1 2 , x i ) + u ( τ n + 1 2 , x i + 1 ) + r w ( τ n + 1 2 , x i 1 ) + 10 w ( τ n + 1 2 , x i ) + w ( τ n + 1 2 , x i + 1 ) = 6 σ 2 h 2 w ( τ n + 1 2 , x i 1 ) 2 w ( τ n + 1 2 , x i ) + w ( τ n + 1 2 , x i + 1 ) .
Note that
ξ c ( τ n + 1 2 ) 12 u x x ( τ n + 1 2 , x i 1 ) + 10 u x x ( τ n + 1 2 , x i ) + u x x ( τ n + 1 2 , x i + 1 ) = ξ c ( τ n + 1 2 ) h 2 u ( τ n + 1 2 , x i 1 ) 2 u ( τ n + 1 2 , x i ) + u ( τ n + 1 2 , x i + 1 ) .
With the CN method, we then obtain the discrete equation for the delta sensitivity:
w τ ( τ n + 1 , x i 1 ) w τ ( τ n , x i 1 ) 12 k + 10 w τ ( τ n + 1 , x i ) w τ ( τ n , x i ) 12 k + w τ ( τ n + 1 , x i + 1 ) w τ ( τ n , x i + 1 ) 12 k = σ 2 4 h 2 w ( τ n + 1 , x i 1 ) 2 w ( τ n + 1 , x i ) + w ( τ n + 1 , x i + 1 ) + σ 2 h 2 w ( τ n , x i 1 ) 2 w ( τ n , x i ) + w ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) k 2 h 2 [ u ( τ n + 1 , x i 1 ) + u ( τ n , x i 1 ) ] ξ c ( τ n + 1 2 ) k h 2 [ u ( τ n + 1 , x i ) + u ( τ n , x i ) ] ξ c ( τ n + 1 2 ) k 2 h 2 [ u ( τ n + 1 , x i + 1 ) + u ( τ n , x i + 1 ) ] r k 24 [ w ( τ n + 1 , x i 1 + w ( τ n , x i 1 ) ] 10 r k 24 [ w ( τ n + 1 , x i ) + w ( τ n , x i ) ] r k 24 [ w ( τ n + 1 , x i + 1 ) + w ( τ n , x i + 1 ) ] .
Simplifying further, we obtain a discrete equation for the delta sensitivity:
1 12 σ 2 k 4 h 2 + r k 24 w ( τ n + 1 , x i 1 ) + 10 12 + σ 2 k 2 h 2 + 10 r k 24 w ( τ n + 1 , x i ) + 1 12 σ 2 k 4 h 2 + r k 24 w ( τ n + 1 , x i + 1 ) = 1 + σ 2 k 4 h 2 r k 24 w ( τ n , x i 1 ) + 1 σ 2 k 2 h 2 10 r k 24 w ( τ n , x i ) + 1 + σ 2 k 4 h 2 r k 24 w ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) k 2 h 2 [ u ( τ n + 1 , x i 1 ) + u ( τ n , x i 1 ) ] ξ c ( τ n + 1 2 ) k h 2 [ u ( τ n + 1 , x i ) + u ( τ n , x i ) ] + ξ c ( τ n + 1 2 ) k 2 h 2 [ u ( τ n + 1 , x i + 1 ) + u ( τ n , x i + 1 ) ] .
Similar to our description for the option value, we obtain a scheme for the delta sensitivity as follows:
a i , i 1 w ( τ n + 1 , x i 1 ) + a i , i w ( τ n + 1 , x i ) + a i , i + 1 w ( τ n + 1 , x i + 1 ) = b i , i 1 w ( τ n , x i 1 ) + b i , i w ( τ n , x i ) + b i , i + 1 w ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) k h 2 [ u ( τ n + 1 , x i 1 ) + u ( τ n , x i 1 ) ] ξ c ( τ n + 1 2 ) k 2 h 2 [ u ( τ n + 1 , x i ) + u ( τ n , x i ) ] + ξ c ( τ n + 1 2 ) k h 2 [ u ( τ n + 1 , x i + 1 ) + u ( τ n , x i + 1 ) ] ,
which can be solved using the Thomas algorithm. In summary, our low-order (1,2) predictor–three-step-corrector method with the fourth-order compact finite difference scheme enables us to approximate the optimal exercise boundary simultaneously with the option value and delta sensitivity by solving the set of discrete equations in (31), (41), and (51).

2.2. Stability Analysis

Here, we carry out a stability analysis of the order (1,2) predictor–corrector fourth-order compact finite difference scheme. Because we use a one-step predictor scheme and three implicit correction steps at each time-level, it is observable that the influence of the correction scheme will become dominant. Moreover, the optimal exercise boundary, left boundary values and the time dependent coefficient in the model are predicted and further corrected at each stage of the three implicit correction step based on the optimal exercise boundary equation and the scheme presented in (31) and (27), respectively. To this end, we first investigate the stability of the high-order correction scheme with three iterative steps as given below.
u ( τ n + 1 , x i 1 ) u ( τ n , x i 1 ) 12 k + 10 u ( τ n + 1 , x i ) u ( τ n , x i ) 12 k + u ( τ n + 1 , x i + 1 ) u ( τ n , x i + 1 ) 12 k = σ 2 4 h 2 u ( τ n + 1 , x i 1 ) 2 u ( τ n + 1 , x i ) + u ( τ n + 1 , x i + 1 ) + σ 2 4 h 2 u ( τ n , x i 1 ) 2 u ( τ n , x i ) + u ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) 24 [ w ( τ n + 1 , x i 1 + w ( τ n , x i 1 ) ] 10 ξ c ( τ n + 1 2 ) 24 [ w ( τ n + 1 , x i ) + w ( τ n , x i ) ] + ξ c ( τ n + 1 2 ) 24 [ w ( τ n + 1 , x i + 1 ) + w ( τ n , x i + 1 ) ] r 24 [ u ( τ n + 1 , x i 1 + u ( τ n , x i 1 ) ] 10 r 24 [ u ( τ n + 1 , x i ) + u ( τ n , x i ) ] r 24 [ u ( τ n + 1 , x i + 1 ) + u ( τ n , x i + 1 ) ] .
w ( τ n + 1 , x i 1 ) w ( τ n , x i 1 ) 12 k + 10 w ( τ n + 1 , x i ) w ( τ n , x i ) 12 k + w ( τ n + 1 , x i + 1 ) w ( τ n , x i + 1 ) 12 k = σ 2 4 h 2 w ( τ n + 1 , x i 1 ) 2 w ( τ n + 1 , x i ) + w ( τ n + 1 , x i + 1 ) + σ 2 h 2 w ( τ n , x i 1 ) 2 w ( τ n , x i ) + w ( τ n , x i + 1 ) + ξ c ( τ n + 1 2 ) k 2 h 2 [ u ( τ n + 1 , x i 1 ) + u ( τ n , x i 1 ) ] ξ c ( τ n + 1 2 ) k h 2 [ u ( τ n + 1 , x i ) + u ( τ n , x i ) ] + ξ c ( τ n + 1 2 ) k 2 h 2 [ w ( τ n + 1 , x i + 1 ) + w ( τ n , x i + 1 ) ] r k 24 [ w ( τ n + 1 , x i 1 ) + w ( τ n , x i 1 ) ] 10 r k 24 [ w ( τ n + 1 , x i ) + w ( τ n , x i ) ] r k 24 [ w ( τ n + 1 , x i + 1 ) + w ( τ n , x i + 1 ) ] .
Our analysis follows the matrix form of the von Neumann stability analysis and we let the Fourier modes
u i n = υ n e i j β h , w i n = ψ n e i j β h , j = 1 .
Note that
ξ ( τ n + 1 2 ) = r + 2 s f ( τ n + 1 ) s f ( τ n ) k s f ( τ n + 1 ) + s f ( τ n ) σ 2 2 .
Let μ = σ 2 k 4 h 2 > 0 . Substituting (54) into (52) and (51), we obtain
υ n + 1 1 1 3 sin 2 β h 2 + 4 μ sin 2 β h 2 + r k r k 2 sin 2 β h 2 υ n 1 1 3 sin 2 β h 2 4 μ sin 2 β h 2 + r k r k 2 sin 2 β h 2 ξ ( τ n + 1 2 ) ψ n + 1 1 2 1 6 sin 2 β h 2 + ψ n 1 2 1 6 sin 2 β h 2 = 0 ,
ψ n + 1 1 1 3 sin 2 β h 2 + 4 μ sin 2 β h 2 + r k r k 2 sin 2 β h 2 ψ n 1 1 3 sin 2 β h 2 4 μ sin 2 β h 2 + r k r k 2 sin 2 β h 2 ξ ( τ n + 1 2 ) 2 ξ ( τ n + 1 2 ) h 2 υ n + 1 sin 2 β h 2 2 ξ ( τ n + 1 2 ) h 2 υ n sin 2 β h 2 = 0 .
If we denote
ϕ 1 = 1 1 3 sin 2 β h 2 + 4 μ sin 2 β h 2 + r k r k 2 sin 2 β h 2 ,
ϕ 2 = 1 1 3 sin 2 β h 2 4 μ sin 2 β h 2 + r k r k 2 sin 2 β h 2 ,
ϕ 3 ( τ n ) = ξ ( τ n + 1 2 ) 1 2 1 6 sin 2 β h 2 , ϕ 4 ( τ n ) = 2 ξ ( τ n + 1 2 ) h 2 sin 2 β h 2 ;
then we obtain two systems of equations as follows:
ϕ 1 υ n + 1 ϕ 2 υ n = ϕ 3 ( τ n ) ψ n + 1 + ϕ 3 ( τ n ) ψ n ,
ϕ 1 ψ n + 1 ϕ 2 υ n = ϕ 4 ( τ n ) υ n + 1 + ϕ 4 ( τ n ) υ n .
Presenting (61) and (62) in matrix form, we obtain
ϕ 1 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 υ n + 1 ψ n + 1 = ϕ 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 2 υ n ψ n .
Denote
A n = ϕ 1 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 , B n = ϕ 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 2 ;
Thus, we have
υ n + 1 ψ n + 1 = A n 1 B n υ n ψ n .
Here, it is important to observe that the two matrices A n and B n are not constant due to the time-dependent coefficients. We have to ensure that A n is invertible. This indicates that the determinant of A n should not be zero at any time-level n. Note that
0 sin 2 β h 2 1 , ϕ 3 ( τ n ) > 0 ;
ϕ 3 ( τ n ) ϕ 4 ( τ n ) = 1 2 1 6 sin 2 β h 2 2 ξ 2 ( τ n + 1 2 ) h 2 sin 2 β h 2 0 .
Furthermore, because μ > 0 , we have ϕ 1 ϕ 2 , implying
| A n | = ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) > 0 , | B n | = ϕ 2 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) > 0 .
Hence, we can invert both matrices A n and B n even though both matrices are time-dependent, which gives
υ n + 1 ψ n + 1 = A n 1 B n υ n ψ n = 1 ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 ϕ 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 2 υ n ψ n .
Simplifying further, we obtain
υ n + 1 ψ n + 1 = ϕ 1 ϕ 2 + ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 ϕ 3 ( τ n ) + ϕ 2 ϕ 3 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 2 ϕ 4 ( τ n ) + ϕ 1 ϕ 4 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 3 ( τ n ) ϕ 4 ( τ n ) + ϕ 1 ϕ 2 ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) υ n ψ n = C n υ n ψ n ,
where C n is the amplification matrix. Let λ n represent the eigenvalue of the matrix C n . To confirm the unconditional stability of the coupled implicit discrete system, we need to show that | λ n | 1 . It can be seen that λ n satisfies
λ n 2 2 λ n ϕ 1 ϕ 2 + ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) + [ ϕ 1 ϕ 2 + ϕ 3 ( τ n ) ϕ 4 ( τ n ) ] 2 [ ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ] 2 [ ϕ 2 ϕ 4 ( τ n ) + ϕ 1 ϕ 4 ( τ n ) ] [ ϕ 1 ϕ 3 ( τ n ) + ϕ 2 ϕ 3 ( τ n ) ] [ ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ] 2 = 0 ,
in which the solutions are
λ n ( 1 , 2 ) = ϕ 1 ϕ 2 + ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ± ( ϕ 1 + ϕ 2 ) ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) .
Note that we obtain ϕ 3 ( τ n ) ϕ 4 ( τ n ) 0 and the complex solutions
λ n ( 1 , 2 ) = ϕ 1 ϕ 2 + ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ± j ( ϕ 1 + ϕ 2 ) ϕ 3 ( τ n ) ϕ 4 ( τ n ) ϕ 1 2 ϕ 3 ( τ n ) ϕ 4 ( τ n ) ,
where j = 1 . For simplicity, let κ ( τ n ) = ϕ 3 ( τ n ) ϕ 4 ( τ n ) . Hence, we obtain
| λ n ( 1 , 2 ) | = [ ϕ 1 ϕ 2 κ ( τ n ) ] 2 [ ϕ 1 2 + κ ( τ n ) ] 2 + ( ϕ 1 + ϕ 2 ) 2 κ ( τ n ) [ ϕ 1 2 + κ ( τ n ) ] 2 = [ ϕ 1 2 + κ ( τ n ) ] [ ϕ 2 2 + κ ( τ n ) ] [ ϕ 1 2 + κ ( τ n ) ] 2 .
Simplifying further, we then obtain that
| λ n ( 1 , 2 ) | = ϕ 2 2 + κ ( τ n ) ϕ 1 2 + κ ( τ n ) 1 , n .
Hence, the interior implicit three-step correction scheme based on the second-order CN time integration method and fourth-order compact finite difference scheme is unconditionally stable.
Next, we navigate the stability of the boundary Euler scheme for predicting the optimal exercise boundary, left boundary values of the option value and the delta sensitivity, and the time-dependent coefficient of the convective term. To this end, we recall the optimal exercise boundary predictor equation:
s f ( p ) ( τ n + 1 ) = H 1 ( τ n ) s f ( τ n ) ,
where
H 1 ( τ n ) = σ 2 k 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) a 3 ( τ n ) + a 4 ( τ n ) a 1 K a 2 s f ( τ n ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
To ensure the stability of (71), we need
| H 1 ( τ n ) | = | σ 2 k 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) a 3 ( τ n ) + a 4 ( τ n ) a 1 K a 2 s f ( τ n ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) | < 1 ,
which implies
1 < σ 2 k 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) a 3 ( τ n ) + a 4 ( τ n ) a 1 K a 2 s f ( τ n ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) < 1 .
From (30), we can further deduce the term in (74) as follows:
a 3 ( τ n ) + a 4 ( τ n ) a 1 K = ( r k 1 ) [ a 1 K a 2 s f ( τ n ) ] + ν k 1 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) a 1 K = r k 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) + ν k 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) ,
where ν = r σ 2 2 . For simplicity, let
a 5 ( τ n ) = a 2 s f ( τ n ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) = 3 h + h 2 s f ( τ n ) 7 2 w ( τ n , x 0 ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
Note that the delta sensitivity is monotonically increasing and non-positive with
s f ( τ n ) = w ( τ n , x 0 ) w ( τ n , x i ) 0 , n , i 0 .
Moreover, the optimal exercise boundary is monotonically decreasing with
0 < s f ( ) s f ( τ n ) s f ( 0 ) , s f ( 0 ) = K , s f ( ) = γ γ + 1 K , γ = 2 r σ 2 .
We refer the reader to the work of Kwok [14], Musiela [15], and Zhang et al. [16] for details. Thus, we have
7 2 w ( τ n , 0 ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) > 0 , ( 3 h + h 2 ) s f ( τ n ) > 0 , h > 0 , n ,
and a 5 ( τ n ) 0 . Denote
Q u ( τ n ) = 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) ,
Q w ( τ n ) = 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
Then, (74) becomes
1 < 1 k σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) a 5 ( τ n ) < 1 .
Note that
1 k σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) a 5 ( τ n ) < 1 ,
implying that
k σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) a 5 ( τ n ) > 0 .
To ensure that k > 0 , we have to confirm
σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) a 5 ( τ n ) > 0 .
From (76)–(79), we have confirmed that a 5 ( τ n ) > 0 . Moreover, we see that
σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) = 4 u τ ( τ n , x 1 ) 1 2 u τ ( τ n , x 2 ) Q w ( τ n ) s f ( τ n ) s f ( τ n ) + O ( h 4 ) = a 2 + Q w ( τ n ) s f ( τ n ) s f ( τ n ) + O ( h 4 ) = a 2 s f ( τ n ) + Q w ( τ n ) s f ( τ n ) s f ( τ n ) + O ( h 4 ) .
Since the first derivative of the optimal exercise boundary s f ( τ n ) is non-positive τ n , when h is very small, we obtain
σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) = a 2 s f ( τ n ) + Q w ( τ n ) s f ( τ n ) s f ( τ n ) > 0 ,
and
σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) a 5 ( τ n ) > 0 .
Next, we verify the right bound as given below
1 < 1 k σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) a 5 ( τ n ) .
Simplifying further, we then obtain
k < 2 a 5 ( τ n ) σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) .
Considering the non-positivity of the delta sensitivity coupled with an extrapolated Taylor series expansion, we can further obtain
4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) = 7 2 w ( τ n , x 0 + h w x ( τ n , x 0 ) + O ( h 2 ) .
Considering the left boundary value of the delta sensitivity and implied second derivative left boundary condition, we further obtain
4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) = 7 2 s f ( τ n ) + 6 r h σ 2 K 3 h s f ( τ n ) + O ( h 2 ) .
Hence, with (76) and a very small h, we can obtain
a 5 ( τ n ) = 7 2 + 3 h + h 2 s f ( τ n ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) = 7 2 + 3 h + h 2 s f ( τ n ) 7 2 s f ( τ n ) + 6 r h σ 2 K 3 h s f ( τ n ) = h 2 s f ( τ n ) + 6 r h σ 2 K h 2 + 6 r h σ 2 K .
Substituting to (90), we then obtain
k < 2 h 2 + 6 r h σ 2 K σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r Q u ( τ n ) ν Q w ( τ n ) .
Note that
σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) = σ 2 2 h 2 8 u ( τ n , x 0 ) 17 u ( τ n , x 1 ) + 10 u ( τ n , x 2 ) u ( τ n , x 3 ) + σ 2 2 h 2 2 h w ( τ n , x 0 ) h 4 w ( τ n , x 1 ) 2 h w ( τ n , x 2 ) + h 4 w ( τ n , x 3 ) .
Let
P u ( τ n ) = 8 u ( τ n , x 0 ) 17 u ( τ n , x 1 ) + 10 u ( τ n , x 2 ) u ( τ n , x 3 ) ,
P w ( τ n ) = 2 w ( τ n , x 0 ) 1 4 w ( τ n , x 1 ) 2 w ( τ n , x 2 ) + 1 4 w ( τ n , x 3 ) .
Thus,
k < 4 h 2 + 6 r h σ 2 K P u ( τ n ) + h P w ( τ n ) + 2 r σ 2 h 2 Q u ( τ n ) 2 σ 2 h 2 ν Q w ( τ n ) h 2 σ 2 = C ( h , τ n ) h 2 σ 2 .
Here, we need to further obtain at least a reasonable upper bound for C ( h , τ n ) in (98) given as
C ( h , τ n ) = 4 h 2 + 6 r h σ 2 K P u ( τ n ) + h P w ( τ n ) + 2 r σ 2 h 2 Q u ( τ n ) 2 σ 2 h 2 ν Q w ( τ n ) .
To this end, we further simplify the denominator of (99). Let x i x 0 = i h . If we derive Taylor series expansion around x 0 , we further obtain the following:
u ( τ n , x i ) = u ( τ n , x 0 ) + i h u x ( τ n , x 0 ) + i 2 h 2 2 u x x ( τ n , x 0 ) + O ( h 3 )
h w ( τ n , x i ) = h w ( τ n , x 0 ) + i h 2 u x x ( τ n , x 0 ) + O ( h 3 )
If we consider (100), (101), left boundary values of the option value and delta sensitivity, and the implied left boundary value of the second derivative of the option value, we can further simplify the terms in P u ( τ n ) and P w ( τ n ) given in (96) and (97) as follows:
17 u ( τ n , x 1 ) = 17 + 17 r h 2 σ 2 K + 17 + 17 h + 17 h 2 2 s f ( τ n ) + O ( h 3 ) ,
10 u ( τ n , x 2 ) = 10 + 40 r h 2 σ 2 K 10 + 20 h + 20 h 2 s f ( τ n ) + O ( h 3 ) ,
u ( τ n , x 3 ) = 1 + 9 r h 2 σ 2 K + 1 + 3 h + 9 h 2 2 s f ( τ n ) + O ( h 3 ) .
h 4 w ( τ n , x 1 ) = r h 2 2 σ 2 K + 1 4 h + h 2 s f ( τ n ) + O ( h 3 ) ,
2 h w ( τ n , x 2 ) = 8 r h 2 σ 2 K + 2 h + 4 h 2 s f ( τ n ) + O ( h 3 ) ,
h 4 w ( τ n , x 3 ) = 3 r h 2 2 σ 2 K 1 4 h + 3 h 2 s f ( τ n ) + O ( h 3 ) ,
Note that
8 u ( τ n , x 0 ) = 8 K 8 s f ( τ n ) , 2 h w ( τ n , x 0 ) = 2 h s f ( τ n ) .
Substituting them into (96) and (97), we obtain
P u ( τ n ) = 14 r h 2 σ 2 K 7 h 2 s f ( τ n ) + O ( h 3 ) ,
h P w ( τ n ) = 7 r h 2 σ 2 K + 7 h 2 2 s f ( τ n ) + O ( h 3 ) .
Furthermore, for a very small h, we have
2 r σ 2 h 2 Q u ( τ n ) = 2 r σ 2 h 2 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) = 2 r σ 2 h 2 7 2 + 2 r h 2 σ 2 K 7 2 + 3 h + h 2 s f ( τ n ) = 4 r 2 h 4 σ 4 K + 7 r h 2 σ 2 K 7 r h 2 σ 2 s f ( τ n ) 6 r h 3 σ 2 s f ( τ n ) 2 r h 4 σ 2 s f ( τ n ) 4 r 2 h 4 σ 4 K 6 r h 3 σ 2 K 2 r h 4 σ 2 K ,
2 h 2 ν σ 2 Q w ( τ n ) = 2 h 2 ν σ 2 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) = 2 h 2 r σ 2 2 σ 2 6 r h σ 2 K 7 2 + 3 h s f ( τ n ) .
With further simplification of (112), we obtain
2 h 2 ν σ 2 Q w ( τ n ) = 12 h 3 r 2 σ 4 K + 7 r h 2 σ 2 s f ( τ n ) + 6 r h 3 σ 2 s f ( τ n ) + 6 r h 3 σ 2 K 7 h 2 2 s f ( τ n ) 3 h 3 s f ( τ n ) .
Hence, we have
2 r σ 2 h 2 Q u ( τ n ) 2 h 2 ν σ 2 Q w ( τ n ) 4 r 2 h 4 σ 4 K 2 r h 4 σ 2 K 12 h 3 r 2 σ 4 K + 7 r β h 2 σ 2 K + 6 β r h 3 σ 2 K 7 h 2 2 K 3 h 3 K .
Here, β = γ γ + 1 , where β is defined in (78). For a very small h, we further obtain
P u ( τ n ) + h P w ( τ n ) = 7 r h 2 σ 2 K 7 h 2 2 s f ( τ n ) 7 r h 2 σ 2 7 h 2 2 K .
Thus, we have
P u ( τ n ) + h P w ( τ n ) + 2 r σ 2 h 2 Q u ( τ n ) 2 h 2 ν σ 2 Q w ( τ n ) 7 h 2 σ 2 2 r β σ 2 K + 3 h 3 σ 4 2 β r σ 2 4 r 2 σ 4 K + 2 r h 4 σ 4 2 r σ 2 K ,
and
C ( h , τ n ) 4 σ 2 6 r + σ 2 h 7 h σ 2 2 r β σ 2 + 3 h 2 2 β r σ 2 4 r 2 σ 4 + 2 r h 3 2 r σ 2 .
To conclude our proof, we need to ensure that
7 h σ 2 2 r β σ 2 + 3 h 2 2 β r σ 2 4 r 2 σ 4 + 2 r h 3 2 r σ 2 > 0 ,
which implies that
2 r β > max σ 2 , 4 r 2 + σ 4 σ 2 .
Hence, we conclude that if the assumption in (119) holds, then
k < 4 6 r + σ 2 h 7 h σ 2 2 r β σ 2 + 3 h 2 2 β r σ 2 4 r 2 σ 4 + 2 r h 3 2 r σ 2 h 2 .
Remark 2.
For the remainder of the low- and high-order predictor–corrector schemes considered in this section and the following section, we will not take into account the stability analysis. We limit our stability analysis to the (1,2) Euler-CN predictor–corrector compact finite difference scheme described above.

2.3. Order (2,2) Leapfrog-CN PC Method with Fourth-Order Compact Scheme

For the Leapfrog predictor scheme, we obtain
σ 2 k 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) = 4 u ( τ n + 1 , x 1 ) u ( τ n 1 , x 1 ) 1 2 u ( τ n + 1 , x 2 ) u ( τ n 1 , x 2 ) s f ( τ n + 1 ) s f ( τ n 1 ) s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) 2 ν k 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + 2 r k 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) + O ( h 4 + k 3 ) .
Simplifying further, we obtain
σ 2 k 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) = a 1 K a 2 s f ( τ n + 1 ) + ( 2 r k 1 ) a 1 K + 2 a 2 r k s f ( τ n ) a 2 s f ( τ n 1 ) s f ( τ n + 1 ) s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) 2 ν k s f ( τ n 1 ) s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + O ( h 4 + k 3 ) .
Let
a 6 ( τ n | τ n 1 ) = a 1 K + ( 2 r k 1 ) a 1 K + 2 a 2 r k s f ( τ n ) a 2 s f ( τ n 1 ) ,
a 7 ( τ n | τ n 1 ) = 2 ν k s f ( τ n 1 ) s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
Subsequently, the optimal exercise boundary is then predicted using the second-order Leapfrog time integration method and fourth-order difference scheme as follows:
s f ( p ) ( τ n + 1 ) = H 2 ( τ n | τ n 1 ) s f ( τ n ) + O ( h 4 + k 3 ) ,
H 2 ( τ n | τ n 1 ) = σ 2 k 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) a 6 ( τ n | τ n 1 ) + a 7 ( τ n | τ n 1 ) a 2 s f ( τ n ) + 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
Remark 3.
In our numerical experiment and study, we initialize the (2,2) Leapfrog-CN predictor–corrector fourth-order compact scheme with the (1,2) Euler-CN predictor–corrector fourth-order compact scheme.

3. Order (3,3) and (4,4) Predictor–Corrector Compact Differencing

In this subsection, we extend our research study to high-order predictor–corrector methods for pricing free boundary options using fourth-order compact finite difference scheme. To this end, we present a suite of order (3,3) and (4,4) predictor–corrector schemes. Here, our choices of (4,4) predictor–corrector schemes exclude Adam Bashforth (AB) and Adam–Moulton (AM) methods because of the irregularity that is inherent in the locality of τ = 0 and x = 0 and the singularity of the s f ( 0 ) therein. Thus, we want to avoid, as much as possible, using more previous derivative (in time) of the optimal exercise boundary, option value, and delta sensitivity for predicting and correcting the present value, especially when the latter is close to the terminal point. Implementing the above-mentioned high-order schemes (AB and AM) adaptively can be more suitable for allowing optimal selection of the time step at each time-level. However, this is not within the scope of our research exploration of (3,3) and (4,4) predictor–corrector schemes.

3.1. Order (3,3) PC Method with Fourth-Order Compact Scheme

To formulate order (3,3) predictor–corrector methods with a fourth-order compact finite difference scheme, we first derive an equation that computes the first derivative of the optimal exercise boundary, which can be used to set up the optimal exercise boundary predictor scheme. To this end, we take the derivative of (24) and obtain
4 u τ ( τ n , x 1 ) 1 2 u τ ( τ n , x 2 ) = 7 2 + 3 h s f ( τ n ) + O ( h 4 ) .
Furthermore, considering the PDE governing the fixed free boundary American options, we obtain
σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) = σ 2 2 4 u x x ( τ n , x 1 ) 1 2 u x x ( τ n , x 2 ) = 4 u τ ( τ n , x 1 ) 1 2 u τ ( τ n , x 2 ) 1 s f ( τ n ) s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) ν 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + r 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) = 7 2 + 3 h + h 2 + 1 s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) s f ( τ n ) ν 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) + r 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) .
Further simplification reveals that
s f ( τ n ) = F ( τ n ) G ( τ n ) ,
with
F ( τ n ) = σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) r 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) + ν 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) ,
G ( τ n ) = 7 2 + 3 h + h 2 + 1 s f ( τ n ) 4 w ( τ n , x 1 ) 1 2 w ( τ n , x 2 ) .
For the third-order predictor scheme, as mentioned earlier, we avoid schemes that require more of the previous derivatives of the optimal exercise boundary, like the high-order Adam–Bashforth methods. Here, we consider two third-order predictor time integration methods. The first predictor scheme requires only the first-order derivative of the optimal exercise boundary at level n. To this end, we present the following Lemma.
Lemma 2.
Assume that  s f ( τ n ) C 1 ( ( 0 , T ] ) . It holds that
s f ( p ) ( τ n + 1 ) = 3 2 s f ( τ n ) + 3 s f ( τ n 1 ) 1 2 s f ( τ n 2 ) + 3 k s f ( p ) ( τ n ) + O ( k 4 ) .
Proof. 
Let τ n + m τ n = m k . The Taylor series expansion around τ n gives
s f ( τ n + 1 ) = s f ( τ n ) + k s f ( τ n ) + k 2 2 s f ( τ n ) + k 3 6 s f ( τ n ) + O ( k 4 ) ,
s f ( τ n 1 ) = s f ( τ n ) k s f ( τ n ) + k 2 2 s f ( τ n ) k 3 6 s f ( τ n ) + O ( k 4 ) ,
s f ( τ n 2 ) = s f ( τ n ) 2 k s f ( τ n ) + 2 k 2 s f ( τ n ) 4 k 3 3 s f ( τ n ) + O ( k 4 ) .
Subtracting (133) from (134), we obtain
s f ( τ n + 1 ) s f ( τ n 1 ) = 2 k s f ( τ n ) + k 3 3 s f ( τ n ) + O ( h 4 ) .
Multiplying (134) by 4 and subtracting it from (135), we obtain
4 s f ( τ n 1 ) s f ( τ n 2 ) = 3 s f ( τ n ) 2 k s f ( τ n ) + 2 k 3 3 s f ( τ n ) + O ( k 4 ) .
Multiplying (136) by 2, we further obtain
2 s f ( τ n + 1 ) 2 s f ( τ n 1 ) = 4 k s f ( τ n ) + 2 k 3 3 s f ( τ n ) + O ( k 4 ) .
Subtracting (137) from (138) with some rearrangement of result, we then obtain (132).    □
Here, we label the third-order predictor scheme in the above Lemma as P 3 a . The second third-order predictor scheme requires only two previous derivative of the optimal exercise boundary and is given as follows:
s f ( p ) ( τ n + 1 ) = 4 s f ( τ n ) + 5 s f ( τ n 1 ) + 4 k s f ( p ) ( τ n ) + 2 k s f ( τ n 1 ) + O ( k 4 ) .
One may see the work of Adam [17,18] and Carpenter et al. [19,20] on how to derive this scheme and its extensions. The third-order predictor scheme is labelled as P 3 b . Here,
s f ( p ) ( τ n ) = F ( τ n ) G ( τ n ) , s f ( τ n 1 ) = F ( τ n 1 ) G ( τ n 1 ) ,
F ( τ n ) = σ 2 2 4 U ( τ n , x 1 ) 1 2 U ( τ n , x 2 ) + r 4 u ( τ n , x 1 ) 1 2 u ( τ n , x 2 ) + ν 4 u x ( τ n , x 1 ) 1 2 u x ( τ n , x 2 ) ,
G ( τ n ) = 7 2 + 3 h + h 2 + 1 s f ( τ n ) 4 u x ( τ n , x 1 ) 1 2 u x ( τ n , x 2 ) ,
ξ n + 1 ( p ) = r + s f ( p ) ( τ n ) s f ( p ) ( τ n + 1 ) σ 2 2 .
To implement our high-order corrector scheme, we need to compute the boundary value of the option value from the predicted optimal exercise boundary as follows:
u ( τ n + 1 , x 0 ) = K s f ( p ) ( τ n + 1 ) .
For the third- and fourth-order corrector schemes, we construct the near-boundary scheme to take into account the associated derivative term in the time integration method. To this end, for the option value, we use a one-sided fourth-order combined compact finite difference [9,10,11] as follows:
14 u x x ( τ n , x i ) 5 u x x ( τ n , x i + 1 ) + 4 u x x ( τ n , x i + 2 ) u x x ( τ n , x i + 3 ) = 2 h 2 u ( τ n , x i 1 ) 2 u ( τ n , x i ) + u ( τ n , x i + 1 ) + O ( h 4 ) , i = 1 .
Furthermore, we use a different fourth-order near-boundary combined compact finite difference scheme for the delta sensitivity by considering the following Lemma.
Lemma 3.
Assume that  u ( τ n , x i ) C 1 , 3 ( ( 0 , T ] , ( 0 , x M ] ) . It holds that
u x x x ( τ n , x i ) = 15 2 h 3 u ( τ n , x i + 1 ) u ( τ n , x i 1 ) ] 3 2 h 2 u x ( τ n , x i 1 ) + 8 u x ( τ n , x i ) + u x ( τ n , x i + 1 ) + O ( h 4 ) , i = 1 .
Proof. 
Let x i + j x i = j h . The Taylor series expansions give
u ( τ n , x i 1 ) u ( τ n , x i + 1 ) = 2 h u x ( τ n , x i ) + h 3 3 u x x x ( τ n , x i ) + h 5 60 u x x x x x ( τ n , x i ) + O ( h 7 ) ,
u ( τ n , x i 1 ) + u ( τ n , x i + 1 ) = 2 u ( τ n , x i ) + h 2 u x x ( τ n , x i ) + h 4 12 u x x x x ( τ n , x i ) + O ( h 6 ) .
Taking derivative of (148) with respect to x and multiplying it by h 5 , we obtain
h 5 u x ( τ n , x i 1 ) + h 5 u x ( τ n , x i + 1 ) = 2 h 5 u x ( τ n , x i ) + h 3 5 u x x x ( τ n , x i ) + h 60 u x x x x x ( τ n , x i ) + O ( h 7 ) .
Subtracting (147) from (149) and re-arranging the result, we then obtain (146).    □
Remark 4.
We observed that the fourth-order combined compact scheme derived in the above Lemma has been recently presented and used in the work of Abrahamsen and Fornberg [21] and as such, we acknowledge and reference their work here.
The high-order scheme in (146) enables us to approximate the near boundary value of delta sensitivity using the discrete value of the option value. Moreover, fewer nodal points of the higher derivatives are used when deriving the near-boundary scheme for the delta sensitivity, as shown in (146). It is worth mentioning that for the delta sensitivity, the near-boundary fourth-order combined compact scheme in (146) is very suitable and provides more accurate result when compared with (145). However, for the option value, we still use (145) when accounting for the near-boundary scheme. For the rest of the interior points for the option value and delta sensitivity, we use the same fourth-order compact scheme presented in (35). Hence, we obtain two systems as follows:
A u u x x = B u u + b u , A w w x x = B w w + b w ,
where A u , B u , b u , A w , B w , b w are given as follows:
A u = 14 5 4 1 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 ,
A w = 10 0 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 1 0 0 0 0 0 0 1 10 ;
B u = 1 h 2 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 ,
B w = 1 h 2 120 15 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 12 0 0 0 0 0 0 12 24 ,
b u = 12 h 2 u 0 n 0 0 0 0 0 0 0 0 , b w = 75 h 3 u 2 n u 0 n 15 h 2 w 0 n 0 0 0 0 0 0 0 0 .
To further show the importance of the presented fourth-order near-boundary scheme for the delta sensitivity, we compute the condition number of A u and A w using the cond function in MATLAB. We observe that the condition number of matrix A w is 1.49997, while the condition number of matrix A u is 2.02913. It is seen that the condition number of A w is less than the condition number of A w , hence validating the benefit of our presented fourth-order near-boundary scheme for the delta sensitivity. With the discrete matrix system above, we then obtain the semi-discrete system for the option value as given below:
u τ n = σ 2 2 A u 1 B u u n + b u n + ξ n w n r u n + O ( h 4 ) .
For corrections using a high-order implicit scheme, we consider two methods that require fewer derivatives of the option value and delta sensitivity. The first third-order implicit corrector scheme is the BDF3 method.
u n + 1 = 18 11 u n 9 11 u n 1 + 2 11 u n 2 + 6 11 k u τ n + 1 + O ( k 3 + h 4 ) .
The second third-order corrector scheme requires only two derivative terms at τ n + 1 and τ n as
u n + 1 = 4 5 u n + 1 5 u n 1 + 2 5 k u τ n + 1 + 4 5 k u τ n + O ( k 3 + h 4 ) .
One may see the work of Adam [17,18] and Carpenter et al. [19,20] on how to derive the scheme in (153) and its extensions. We label it as C 3 . The last third-order predictor scheme we considered is the third-order Adam–Moulton corrector scheme as follows:
u n + 1 = u n + 5 12 k u τ n + 1 + 8 12 k u τ n 1 12 k u τ n 1 + O ( k 3 + h 4 ) .
We label (166) as A M 3 . The optimal exercise boundary, its first derivative, and the non-constant coefficient of the first derivative of the option value are then corrected as follows:
s f ( c ) ( τ n + 1 l ) = M ( τ n + 1 ) a 1 K a 2 + O ( h 4 ) , l = 1 , 2 .
a 1 = 7 2 + 2 r h 2 σ 2 , a 2 = 7 2 + 3 h + h 2 ;
u ( τ n + 1 , x 0 ) = K s f ( c ) ( τ n + 1 ) ,
M ( τ n + 1 ) = 4 u ( τ n + 1 , x 1 ) 1 2 u ( τ n + 1 , x 2 ) ,
s f ( τ n + 1 ) = M τ ( τ n + 1 ) a 2 + O ( h 4 ) ,
M τ ( τ n + 1 ) = 4 u τ ( τ n + 1 , x 1 ) 1 2 u τ ( τ n + 1 , x 2 ) ,
ξ n + 1 ( c ) = r + s f ( τ n + 1 ) s f ( τ n + 1 ) σ 2 2 .
The corrected ξ n + 1 ( c ) and s f ( c ) ( τ n + 1 l ) are then used to obtain the discrete solution of the delta sensitivity:
w ( τ n + 1 , x 0 ) = s f ( c ) ( τ n + 1 ) ,
w τ n = σ 2 2 A w 1 B w w n + b w n + ξ n + 1 ( c ) A u 1 B u u n + b u n r w n + O ( h 4 ) .
For the BDF3, we obtain
w n + 1 = 18 11 w n 9 11 w n 1 + 2 11 w n 2 + 6 11 k w τ n + 1 + O ( k 3 + h 4 ) .
Similarly, for the second third-order corrector scheme
w n + 1 = 4 5 w n + 1 5 w n 1 + 2 5 k w τ n + 1 + 4 5 k w τ n + O ( k 3 + h 4 ) .
The last third-order predictor scheme is the third-order Adam–Moulton corrector scheme.
w n + 1 = w n + 5 12 k w τ n + 1 + 8 12 k w τ n 1 12 k w τ n 1 + O ( k 3 + h 4 ) .
Remark 5.
For the order (3,3) predictor–corrector fourth-order compact finite difference shown above, including the order (4,4) we will present below, we use only two correction iterative steps to correct the option value, delta sensitivity, and optimal exercise boundary, unlike the (1,2) and (2,2) schemes where we implemented three correction iterative steps.

3.2. Order (4,4) PC Method with Fourth-Order Compact Scheme

For the fourth-order predictor scheme, we consider the scheme below [19,20]:
s f ( p ) ( τ n + 1 ) = 9 s f ( τ n ) + 9 s f ( τ n 1 ) + s f ( τ n 2 ) + 6 k s f ( p ) ( τ n ) + 6 k s f ( τ n 1 ) + O ( k 5 ) .
We label it as P 4 . For the fourth-order correction time integration scheme, we consider the BDF4 scheme as follows:
u n + 1 = 48 35 u n 36 25 u n 1 + 16 25 u n 2 3 25 u n 3 + 12 25 k u τ n + 1 + O ( k 4 + h 4 ) ,
and another fourth-order time integration method that require fewer previous first derivative term as follows [22]:
u n + 1 = 9 17 u n + 9 17 u n 1 1 17 u n 2 + 6 17 k u τ n + 1 + 18 17 k u τ n + O ( k 4 + h 4 ) .
We label the correction scheme in (169) as C 4 . The choice of BDF4 and C 4 as the corrector scheme over fourth-order AB and AM is mainly to avoid using more of the previous first derivatives (in time) of the option value close to τ = 0 . This is similar to our reason for deriving a new third-order predictor scheme for this work, as presented in (132).
Similarly, for the delta sensitivity, we obtain
w n + 1 = 48 35 w n 36 25 w n 1 + 16 25 w n 2 3 25 w n 3 + 12 25 k w τ n + 1 + O ( k 4 + h 4 ) ,
w n + 1 = 9 17 w n + 9 17 w n 1 1 17 w n 2 + 6 17 k w τ n + 1 + 18 17 k w τ n + O ( k 4 + h 4 ) .
For the initial procedure, we use the SSPRK3 method for both order (3,3) and (4,4) predictor–corrector fourth-order compact finite difference schemes.
Remark 6.
It is worth mentioning that we can further recombine our predictor–corrector methods to include order (2,3) and (3,4). Moreover, an adaptive implementation can further be generated with (1,2), (2,3), and (3,4) predictor–corrector schemes. However, for brevity, our description will not include order (2,3) and (3,4) predictor–corrector compact finite difference schemes.

4. Numerical Results

The experiments were performed on a computer with a 12th Gen Intel(R) Core(TM) i7-12700H at 2.30 GHz on a 64-bit Windows 11 operating system. Furthermore, MATLAB Programming Language was used for numerical experiments and visualization. We present our findings and results in the following subsections.

4.1. Order (1,2) and (2,2) Predictor–Corrector Schemes

In this subsection, we first consider (1,2) and (2,2) predictor–corrector fourth-order compact finite difference schemes for solving a system of free boundary American put option pricing problem. Here, we would like to verify the performance of the presented PC-compact finite difference scheme in terms of solution accuracy, computational runtime, and convergence rate. Furthermore, we want to compare their results with the existing methods.
To this end, we first consider the example in the work of Zhu and Zhang [6] with the parameters listed in Table 1.
Our focus here is to compare the values of the optimal exercise boundary obtained based on our method and the method of Zhu and Zhang [6]. s f ( T ) = 76.163220 is the benchmark value, which was obtained from the highly accurate sixth-order compact finite difference scheme [23] with step size of h = 0.01 . The solution profiles with the (1,2) predictor–corrector fourth-order compact finite difference scheme were plotted in Figure 1 and Figure 2. The values and computational runtime in seconds are listed in Table 2 and Table 3, respectively.
In Figure 1 and Figure 2, we observe the smoothness in the numerical approximation of the optimal exercise boundary, option value, and delta sensitivity. Moreover, we observe that the results of order (1,2) and (2,2) predictor–corrector fourth-order compact finite difference schemes agree well with the result obtained from the method of Zhu and Zhang [6]. It is very important to observe from Table 2 and Table 3 that with small number of grid intervals ( h = 0.025 , k = h 2 ), the order (1,2) and (2,2) predictor–corrector fourth-order compact finite difference schemes provide a solution accuracy much closer to the result of the sixth-order compact finite difference up to 6-digit. Moreover, our (1,2) and (2,2) predictor–corrector schemes are very fast in computation even though our system of PDEs simultaneously solve the option value and delta sensitivity.
Here, it is worth mentioning that we define the runtime in seconds as the time required to run the whole code that approximates the option value, delta sensitivity, and the optimal exercise boundary simultaneously, excluding initialization and visualization. Initialization and visualization take insignificant computational runtime. However, it could be slightly significant here because we achieve reasonable accuracy with a very coarse grid and very little computational runtime. The main reason that enables the order (1,2) Euler-CN and (2,2) Leapfrog-CN predictor–corrector scheme to be very fast is the implementation of the Thomas algorithm in the three-iterative step of the corrector scheme.
To fully describe the computational speed we achieved with order (1,2) Euler-CN and (2,2) Leapfrog-CN predictor–corrector compact finite difference schemes, we computed the optimal exercise boundary, option value, and delta sensitivity simultaneously using a very small-time step k = 4 × 10 6 with varying step sizes. We then compared the results with the one obtained from the sixth-order compact finite difference scheme. The numerical solution of the optimal exercise boundary and the computational runtime with these methods are listed in Table 4 and Table 5.
The importance of our order (1,2) and (2,2) predictor–corrector fourth-order compact finite difference schemes is that if we fix a very small-time step, as shown in Table 5, their computational runtime grows almost insignificantly and linearly as the step size in the space grid is exceedingly reduced. We observe this nice feature across all our numerical experiment with the above low-order (1,2) and (2,2) predictor–corrector fourth-order compact finite difference schemes. It is easy to further observe in Table 4 and Table 5 that both predictor–corrector schemes outperform the sixth-order compact scheme method in terms of computational runtime from h = 0.006 . It is also interesting to see that when h = 0.003 , the results of the (1,2) Euler-CN and (2,2) Leapfrog-CN predictor–corrector schemes are the same as that of the sixth-order compact scheme up to seven and eight digits even though our (1,2) and (2,2) PC schemes are implemented with a fourth-order compact scheme, signaling that our implementation shows a highly accurate result with a very little computational runtime. Moreover, the computational runtime for the sixth-order compact finite difference scheme is significantly high with a very refine space and time grid. Notwithstanding, we have already achieved seven-digit accuracy with the sixth-order compact finite difference method with h = 0.01 and a runtime of 9 seconds, which is the main advantage of the sixth-order compact finite difference scheme over our (1,2) and (2,2) predictor–corrector schemes.
To better understand the importance of the three-step correction scheme we proposed based on the order (1,2) Euler-CN and order (2,2) Leapfrog-CN fourth-order compact finite difference scheme, we compared the computational results of the optimal exercise boundary based on one, two, three, and four iteration steps with the implicit scheme, as described in Table 6. From Table 6, we observe that as the number of iteration increases from 1 to 3, there is a substantial improvement in the solution accuracy. However, beyond a number of iterations of 3, the improvement in the solution accuracy is insignificant. Hence, there is no need to waste more computational effort, which is the reason why we use three iteration steps.
In the second experiment, we compared the results of the option value and delta sensitivity with the existing methods. To this end, we considered the example in the work of Tangman et al. [24] with the parameter given in Table 7. For numerical comparison, we used a benchmark value obtained from the Binomial method with 15,001 steps [25]. We used a step size of h = 0.01 and a time step of k = h 2 . The numerical results are listed in Table 8 and Table 9. For this example, we display the plot profiles for the (2,2) predictor–corrector fourth-order compact finite difference scheme in Figure 3 and Figure 4. The curves of the optimal exercise boundary, option value, and delta sensitivity admit sufficient smoothness. Moreover, we observe a close match between the result obtain from the low-order (1,2) and (2,2) predictor–corrector schemes as compared with the well-known methods, including the benchmark from the Binomial method [25].
In the last experiment, our interest was to compute the convergence rates of the present predictor–corrector schemes. To this end, we considered parameters in the second example. The only difference is that we used different time to maturity, given as T = 0.125 . Here, we fixed the time step k = 10 6 and then varied the step size, as shown in Table 10 and Table 11, where the maximum errors and convergence rate are displayed. The convergence rate results for the low-order (1,2) and (2,2) predictor–corrector schemes are reasonable. Moreover, further improvement can be achieved if the singularity associated with the first derivative of the optimal exercise boundary and inherent discontinuity associated with the model are carefully improved. In general, our proposed low-order (1,2) and (2,2) predictor–corrector fourth-order compact finite schemes which simultaneously solve the optimal exercise boundary, option value, and delta sensitivity are very fast and provide reasonable solution accuracy on both coarse and very refined grids.

4.2. Order (3,3) and (4,4) Predictor–Corrector Schemes

Here, we investigate the performance of the presented (3,3) and (4,4) predictor–corrector integration methods with fourth-order compact finite difference scheme for solving the American option pricing model as a fixed free boundary problem. Consider the example in the work of Bunch and Johnson [28] with the following parameters displayed in Table 12.
The benchmark solution as presented in the work of Bunch and Johnson [28] was obtained from the Binomial tree method with 10,000 steps. Our numerical approximations with this high-order predictor–corrector schemes are presented in Table 13 and Table 14. From Table 13 and Table 14, we observe reasonable solution accuracy across the high-order predictor–corrector schemes. Moreover, for the order (3,3) predictor–corrector fourth-order compact finite difference schemes, we observe that the BDF3 scheme has the lowest accuracy and P 3 b - C 3 has the best accuracy as compared with the benchmark value. For the (4,4) predictor–corrector scheme, P 4 - C 4 performed better than P 4 -BDF4. However, the discrepancy in results among the high-order predictor–corrector compact finite difference schemes is not substantial. Finally, we plotted the curves of the optimal exercise boundary, option value, and delta sensitivity obtained from the P 4 - C 4 predictor–corrector fourth-order compact finite difference schemes in Figure 5, Figure 6 and Figure 7.

5. Conclusions

We have explored the efficiency of a suite of predictor–corrector schemes with fourth-order compact finite difference for solving a system of American option pricing PDEs consisting of the option value and delta sensitivity as a free boundary problem. The numerical results shows that the low-order (1,2) Euler-CN and (2,2) Leapfrog-CN predictor–corrector compact finite difference schemes have great performance in terms of computational speed and accuracy. Moreover, (1,2) Euler-CN and (2,2) Leapfrog-CN are very fast in both coarse and very refined grids and also present highly accurate solutions with coarse grids and little computation runtime. This is largely possible due to the implementation of the Thomas algorithm and fourth-order compact finite difference schemes. Moreover, the speed of our computation is further enhanced by using only three iterative steps at each time-level with the second-order CN and fourth-order compact finite difference correction schemes. We also obtain reasonable convergence rate results with these two implementations. For extension purposes, we have further explored the suite of high-order predictor–corrector schemes and have obtained reasonable solution accuracy. In general, we can validate that our implementations are very efficient and could be extendable to solving other free boundary option pricing problems.
In our future work, we hope to further enhance the performance of our proposed suite of low- and high-order predictor–corrector fourth-order compact finite difference schemes by improving the smoothness and irregularity inherent in the transformed model such that we can recover numerical convergence rate that is consistently in good agreement with the theoretical one. Furthermore, we hope in our future work to extend this method and its variants to non-standard and exotic options including stochastic volatility models.

Author Contributions

Conceptualization, C.N. and W.D.; methodology, C.N.; software, C.N.; validation, C.N.; formal analysis, C.N.; investigation, C.N.; resources, C.N.; data curation, C.N.; writing—original draft preparation, C.N.; writing—review and editing, W.D.; visualization, C.N.; supervision, W.D.; project administration, W.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The first author is funded in part by an NSERC Discovery Grant.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kazmi, K. An IMEX predictor–corrector method for pricing options under regime-switching jump-diffusion models. Int. J. Comput. Math. 2019, 96, 1137–1157. [Google Scholar] [CrossRef]
  2. Kalantari, R.; Shahmorad, S.; Ahmadian, D. The stability analysis of predictor–corrector method in solving American option pricing model. Comput. Econ. 2016, 47, 255–274. [Google Scholar] [CrossRef]
  3. Khaliq, A.Q.M.; Voss, D.A.; Kazmi, S.H.K. A linearly implicit predictor–corrector scheme for pricing American options using a penalty method approach. J. Bank. Financ. 2006, 30, 489–502. [Google Scholar] [CrossRef]
  4. Chen, W.; Xu, X.; Zhu, S.P. A predictor–corrector approach for pricing American options under the finite moment log-stable model. Appl. Numer. Math. 2015, 97, 15–29. [Google Scholar] [CrossRef]
  5. Hajipour, M.; Malek, A. Efficient high-order numerical methods for pricing of options. Comput. Econ. 2015, 45, 31–47. [Google Scholar] [CrossRef]
  6. Zhu, S.P.; Zhang, J. A new predictor–corrector scheme for valuing American puts. Appl. Math. Comput. 2011, 217, 4439–4452. [Google Scholar] [CrossRef]
  7. Zhu, S.P.; Chen, W.T. A predictor–corrector scheme based on the ADI method for pricing American puts with stochastic volatility. Comput. Math. Appl. 2011, 62, 1–26. [Google Scholar] [CrossRef] [Green Version]
  8. Wu, L.; Kwok, Y.K. A front-fixing finite difference method for the valuation of American options. J. Financ. Eng. 1997, 6, 83–97. [Google Scholar]
  9. Zhao, J.; Corless, R.M. Compact finite difference method for integro-differential equations. Appl. Math. Comput. 2006, 177, 271–288. [Google Scholar] [CrossRef]
  10. Zhao, J. Highly accurate compact mixed methods for two point boundary value problems. Appl. Math. Comput. 2007, 188, 1402–1418. [Google Scholar] [CrossRef]
  11. Zhao, J.; Davison, M.; Corless, R.M. Compact finite difference method for American option pricing. J. Comput. Appl. Math. 2007, 206, 306–321. [Google Scholar] [CrossRef] [Green Version]
  12. Egorova, V.N.; Jódar, L. Solving American option pricing models by the front fixing method: Numerical analysis and computing. In Abstract and Applied Analysis; Hindawi: London, UK, 2014; Volume 2014. [Google Scholar]
  13. Egorova, V.N.; Company, R.; Jódar, L. A new efficient numerical method for solving American option under regime switching model. Comput. Math. Appl. 2016, 71, 224–237. [Google Scholar] [CrossRef]
  14. Kwok, Y.K. Mathematical Models of Financial Derivatives; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  15. Musiela, M.; Rutkowski, M. Martingale Methods in Financial Modelling; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  16. Zhang, K.; Song, H.; Li, J. Front-fixing FEMs for the pricing of American options based on a PML technique. Appl. Anal. 2015, 94, 903–931. [Google Scholar] [CrossRef]
  17. Adam, Y. A Hermitian finite difference method for the solution of parabolic equations. Comput. Math. Appl. 1975, 1, 393–406. [Google Scholar] [CrossRef] [Green Version]
  18. Adam, Y. Highly accurate compact implicit methods and boundary conditions. J. Comput. Phys. 1977, 24, 10–22. [Google Scholar] [CrossRef]
  19. Carpenter, M.H.; Gottlieb, D.; Abarbanel, S. The stability of numerical boundary treatments for compact high-order finite-difference schemes. J. Comput. Phys. 1993, 108, 272–295. [Google Scholar] [CrossRef] [Green Version]
  20. Carpenter, M.H.; Gottlieb, D.; Abarbanel, S. Stable and accurate boundary treatments for compact, high-order finite-difference schemes. Appl. Numer. Math. 1993, 12, 55–87. [Google Scholar] [CrossRef]
  21. Abrahamsen, D.; Fornberg, B. Solving the Korteweg-de Vries equation with Hermite-based finite differences. Appl. Math. Comput. 2021, 401, 126101. [Google Scholar] [CrossRef]
  22. Li, M.; Tang, T. A compact fourth-order finite difference scheme for unsteady viscous incompressible flows. J. Sci. Comput. 2001, 16, 29–45. [Google Scholar] [CrossRef]
  23. Nwankwo, C.; Dai, W. Sixth-order compact differencing with staggered boundary schemes and 3(2) Bogacki-Shampine pairs for pricing free-boundary options. arXiv 2022, arXiv:2207.14379. [Google Scholar]
  24. Tangman, D.Y.; Gopaul, A.; Bhuruth, M. A fast high-order finite difference algorithm for pricing American options. J. Comput. Appl. Math. 2008, 222, 17–29. [Google Scholar] [CrossRef] [Green Version]
  25. Leisen, D.P.; Reimer, M. Binomial models for option valuation-examining and improving convergence. Appl. Math. Financ. 1996, 3, 319–346. [Google Scholar] [CrossRef]
  26. Ikonen, S.; Toivanen, J. Operator splitting methods for American option pricing. Appl. Math. Lett. 2004, 17, 809–814. [Google Scholar] [CrossRef] [Green Version]
  27. Brennan, M.J.; Schwartz, E.S. The valuation of American put options. J. Financ. 1977, 32, 449–462. [Google Scholar] [CrossRef]
  28. Bunch, D.S.; Johnson, H. The American put option and its critical stock price. J. Financ. 2000, 55, 2333–2356. [Google Scholar] [CrossRef]
Figure 1. Optimal exercise boundary and option value with order (1,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ).
Figure 1. Optimal exercise boundary and option value with order (1,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ).
Axioms 12 00762 g001
Figure 2. Delta sensitivity with order (1,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ).
Figure 2. Delta sensitivity with order (1,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ).
Axioms 12 00762 g002
Figure 3. Optimal exercise boundary and option value with order (2,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ). s f ( T ) = 76.284497.
Figure 3. Optimal exercise boundary and option value with order (2,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ). s f ( T ) = 76.284497.
Axioms 12 00762 g003
Figure 4. Delta sensitivity with order (2,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ).
Figure 4. Delta sensitivity with order (2,2) predictor–corrector scheme ( h = 0.01 , k = h 2 ).
Axioms 12 00762 g004
Figure 5. Optimal exercise boundary with P 4 - C 4 and varying strike price.
Figure 5. Optimal exercise boundary with P 4 - C 4 and varying strike price.
Axioms 12 00762 g005
Figure 6. Option value with P 4 - C 4 and varying strike price.
Figure 6. Option value with P 4 - C 4 and varying strike price.
Axioms 12 00762 g006
Figure 7. Delta sensitivity with P 4 - C 4 and varying strike price.
Figure 7. Delta sensitivity with P 4 - C 4 and varying strike price.
Axioms 12 00762 g007
Table 1. Option parameters.
Table 1. Option parameters.
  KTr σ x M
  1001.000.100.303.00
Table 2. Optimal exercise boundary with the predictor–corrector schemes (No. of grid interval N = 120 , k = h 2 ).
Table 2. Optimal exercise boundary with the predictor–corrector schemes (No. of grid interval N = 120 , k = h 2 ).
PC(1,2) Euler-CN(2,2) Leapfrog-CNZhu and Zhang [6] (N = 400)
s f ( T ) 76.16322776.16315176.163742
Table 3. Computational runtime(s) of the predictor–corrector schemes (No. of grid interval N = 120 ).
Table 3. Computational runtime(s) of the predictor–corrector schemes (No. of grid interval N = 120 ).
PC(1,2) Euler-CN(2,2) Leapfrog-CN
Total Runtime (s)0.2 s0.2 s
Table 4. Computational runtime(s) and with k = 4 × 10 6 and varying step sizes.
Table 4. Computational runtime(s) and with k = 4 × 10 6 and varying step sizes.
  Methods h = 0.01 h = 0.006 h = 0.003
  (1,2) Euler-CN24 s32 s56 s
  (2,2) Leapfrog-CN21 s26 s51 s
  Sixth-order compact scheme9 s180 s4785 s
Table 5. Optimal exercise boundary with k = 4 × 10 6 and varying step sizes.
Table 5. Optimal exercise boundary with k = 4 × 10 6 and varying step sizes.
  Methods h = 0.01 h = 0.006 h = 0.003
  (1,2) Euler-CN76.16307876.16318876.163226
  (2,2) Leapfrog-CN76.16305976.16318276.163225
  Sixth-order compact scheme76.16322076.16322676.163226
Table 6. Optimal exercise boundary with different number of iterations using Euler-CN predictor–corrector compact scheme ( h = 0.01 ).
Table 6. Optimal exercise boundary with different number of iterations using Euler-CN predictor–corrector compact scheme ( h = 0.01 ).
No. of Iteration1234
s f ( T ) 76.16384176.16338876.16327676.163268
Table 7. Option parameters.
Table 7. Option parameters.
  KTr σ x M
  1003.000.050.203.00
Table 8. Option value with predictor–corrector schemes. WK is the method of Wu and Kwok [8]. OP is the operator splitting method [26]. BS1 and BS2 is the method of Brennan and Schwartz [27].
Table 8. Option value with predictor–corrector schemes. WK is the method of Wu and Kwok [8]. OP is the operator splitting method [26]. BS1 and BS2 is the method of Brennan and Schwartz [27].
SEuler-CNLeapfrog-CNBinomialWKOPBS-1BS-2
8020.282020.282020.279720.282520.279520.278520.2783
9013.307713.307713.307513.311713.307413.304713.3072
1008.71078.71078.71068.71358.71048.70708.7102
1105.68265.68265.68255.68675.68245.67915.6822
1203.69653.69653.69643.70013.69633.69353.6961
Table 9. Delta sensitivity with predictor–corrector schemes. WK is the method of Wu and Kwok [8]. OP is the operator splitting method [26]. BS1 and BS2 is the method of Brennan and Schwartz [27].
Table 9. Delta sensitivity with predictor–corrector schemes. WK is the method of Wu and Kwok [8]. OP is the operator splitting method [26]. BS1 and BS2 is the method of Brennan and Schwartz [27].
SEuler-CNLeapfrog-CNBinomialWKOPBS-1BS-2
80−0.8537−0.8537−0.8536−0.8508−0.8536−0.8539−0.8536
90−0.5619−0.5619−0.5619−0.5600−0.5619−0.5621−0.5619
100−0.3706−0.3706−0.3706−0.3694−0.3706−0.3707−0.3706
110−0.2436−0.2436−0.2436−0.2429−0.2436−0.2436−0.2436
120−0.1594−0.1594−0.1594−0.1589−0.1594−0.1593−0.1594
Table 10. Convergence rate result of the option value with (1,2) Euler-CN compact scheme.
Table 10. Convergence rate result of the option value with (1,2) Euler-CN compact scheme.
  hMaximum ErrorsConvergence Rate
  0.2
  0.13.636750127435
  0.050.5417349218612.747
  0.0250.1696922272371.675
  0.01250.0266272376172.672
  0.006250.0019699274623.757
Table 11. Convergence rate result of the option value with (2,2) Leapfrog-CN compact scheme.
Table 11. Convergence rate result of the option value with (2,2) Leapfrog-CN compact scheme.
  hMaximum ErrorsConvergence Rate
  0.2
  0.13.662480988725
  0.050.5940740514362.624
  0.0250.1531107022381.956
  0.01250.0324270921142.239
  0.006250.0015850821854.355
Table 12. Option parameters.
Table 12. Option parameters.
     S 0 Tr σ x M
  404.000.04880.303.00
Table 13. Option value with the (3,3) and(4,4) predictor–corrector schemes ( h = 0.01 , k = h 2 ). The benchmark value is the Binomial method (BM) with 10,000 steps.
Table 13. Option value with the (3,3) and(4,4) predictor–corrector schemes ( h = 0.01 , k = h 2 ). The benchmark value is the Binomial method (BM) with 10,000 steps.
K P 3 a -BDF3 P 3 b -BDF3 P 3 a - C 3 P 3 b - C 3 P 3 b -M3 P 4 -BDF4 P 4 - C 4 BM
350.69780.69780.69780.69770.69760.69770.69750.6975
402.48372.48362.48322.48322.48332.48332.48302.4825
455.70735.70725.70675.70685.70685.70695.70655.7056
Time (s)5.95.85.75.75.85.95.5None
Table 14. Delta sensitivity with the (3,3) and (4,4) predictor–corrector schemes ( h = 0.01 , k = h 2 ). The benchmark value is the Binomial method (BM) with 10,000 steps.
Table 14. Delta sensitivity with the (3,3) and (4,4) predictor–corrector schemes ( h = 0.01 , k = h 2 ). The benchmark value is the Binomial method (BM) with 10,000 steps.
K P 3 a -BDF3 P 3 b -BDF3 P 3 a - C 3 P 3 b - C 3 P 3 b -M3 P 4 -BDF4 P 4 - C 4 BM
35−0.1740−0.1740−0.1740−0.1740−0.1740−0.1740−0.1740−0.1741
40−0.4418−0.4418−0.4419−0.4419−0.4419−0.4419−0.4419−0.4420
45−0.7262−0.7262−0.7263−0.7263−0.7263−0.7263−0.7264−0.7266
Time (s)5.95.85.75.75.85.95.5None
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Nwankwo, C.; Dai, W. Efficiency of Some Predictor–Corrector Methods with Fourth-Order Compact Scheme for a System of Free Boundary Options. Axioms 2023, 12, 762. https://doi.org/10.3390/axioms12080762

AMA Style

Nwankwo C, Dai W. Efficiency of Some Predictor–Corrector Methods with Fourth-Order Compact Scheme for a System of Free Boundary Options. Axioms. 2023; 12(8):762. https://doi.org/10.3390/axioms12080762

Chicago/Turabian Style

Nwankwo, Chinonso, and Weizhong Dai. 2023. "Efficiency of Some Predictor–Corrector Methods with Fourth-Order Compact Scheme for a System of Free Boundary Options" Axioms 12, no. 8: 762. https://doi.org/10.3390/axioms12080762

APA Style

Nwankwo, C., & Dai, W. (2023). Efficiency of Some Predictor–Corrector Methods with Fourth-Order Compact Scheme for a System of Free Boundary Options. Axioms, 12(8), 762. https://doi.org/10.3390/axioms12080762

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop