Next Article in Journal
A Non-Convex Hybrid Overlapping Group Sparsity Model with Hyper-Laplacian Prior for Multiplicative Noise
Next Article in Special Issue
On Finite-Time Blow-Up Problem for Nonlinear Fractional Reaction Diffusion Equation: Analytical Results and Numerical Simulations
Previous Article in Journal
Fast and Accurate Numerical Algorithm with Performance Assessment for Nonlinear Functional Volterra Equations
Previous Article in Special Issue
Total Value Adjustment of Multi-Asset Derivatives under Multivariate CGMY Processes
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Preconditioned Iterative Method for a Multi-State Time-Fractional Linear Complementary Problem in Option Pricing

1
School of Economics, Guangdong University of Technology, Guangzhou 510520, China
2
Key Laboratory of Digital Economy and Data Governance, Guangdong University of Technology, Guangzhou 510520, China
3
Department of Mathematics, University of Macau, Macau
*
Author to whom correspondence should be addressed.
Fractal Fract. 2023, 7(4), 334; https://doi.org/10.3390/fractalfract7040334
Submission received: 14 March 2023 / Revised: 10 April 2023 / Accepted: 12 April 2023 / Published: 17 April 2023

Abstract

:
Fractional derivatives and regime-switching models are widely used in various fields of finance because they can describe the nonlocal properties of the solutions and the changes in the market status, respectively. The regime-switching time-fractional diffusion equations that combine both advantages are also used in European option pricing; however, to our knowledge, American option pricing based on such models and their numerical methods is yet to be studied. Hence, a fast algorithm for solving the multi-state time-fractional linear complementary problem arising from the regime-switching time-fractional American option pricing models is developed in this paper. To construct the solution strategy, the original problem has been converted into a Hamilton–Jacobi–Bellman equation, and a nonlinear finite difference scheme has been proposed to discretize the problem with stability analysis. A policy-Krylov subspace method is developed to solve the nonlinear scheme. Further, to accelerate the convergence rate of the Krylov method, a tri-diagonal preconditioner is proposed with condition number analysis. Numerical experiments are presented to demonstrate the validity of the proposed nonlinear scheme and the efficiency of the proposed preconditioned policy-Krylov subspace method.

1. Introduction

Option pricing has been a popular research topic in finance, as option instruments are important derivatives in financial markets. The Black–Scholes model [1] is the most famous model for option pricing problems and has been widely used; however, this model has suffered many shortcomings because it is based on geometric Brownian motion and cannot model jumps in the assets’ prices. New models have been proposed to deal with these problems, such as stochastic volatility models [2,3], jump diffusion models [4], variance gamma models [5], tempered stable models [6], self-exciting jump models [7], and Hawkes jump diffusion models [8]. With the Levy processes used in option pricing, new pricing models involving fractional derivatives have emerged.
Option valuation involving fractional derivatives is widely discussed. Cartea and del-Castillo-Negrete find that the barrier option price can be obtained by solving fractional partial differential equations for some particular Levy processes [9], such as finite moment log stable process [10], CGMY process [11], and KoBoL process [6]. Further, spatial fractional derivatives have been used in American options pricing [12], European options pricing [13], and double barrier options pricing [14]; however, these models cannot describe the long-term memory effects of the option value. Hence, the time-fractional derivative started to be used in the option pricing problem, which was introduced by Wyss [15]. After that, two new families of fractional derivatives based on fractional Taylor’s series and modified Riemann–Liouville derivatives were proposed to replace the standard time derivatives to price the options [16]. Additionally, a partial differential equation containing a Caputo fractional derivative was developed to price the European-style options [17]. Subsequently, an option pricing model containing both time and spatial fractional derivatives was proposed in [18].
Although the fractional derivative can describe the nonlocal properties of the assets’ prices and the time-independence properties of the option value, it cannot describe changes in market states, such as bull and bear market transitions. Hence, regime-switching models have been proposed and applied in options pricing [19], time series [20], and credit default swaps [21]. Additionally, the regime-switching option pricing involving fractional derivatives has been developed in [22,23,24], which combines the advantages of both regime-switching and fractional derivatives. For European option pricing, a system of space fractional diffusion equations coupled by Markov terms from the regime-switching finite moment log stable model has been developed in [23,25]. Furthermore, a coupled tempered fractional diffusion equation has been used to value the European options [26,27]. In [28], a system of linear complementary problems (LCPs) involving tempered fractional derivatives generated from the regime-switching CGMY model and KoBoL model has been utilized to price the American options. The above models belong to the regime-switching models with spatial fractional derivatives, and there are relatively few studies of option pricing based on the regime-switching model involving time-fractional derivatives. A regime-switching model containing time-fractional derivatives has been developed in [24] to price the European options, and a fast positivity preserving numerical method has been proposed in [22] to solve it. To the best of our knowledge, no numerical results exist for such a model for pricing the American options. Hence, developing a fast solution strategy for the multi-state time-fractional LCP arising from the regime-switching model with a time-fractional derivative is the main aim of this paper.
As noted above, no study in the literature considers the numerical method for the multi-state time-fractional LCP arising from American option pricing. The development of the fast solution strategy can only refer to the numerical methods for the general American options pricing. The methods for pricing American options can be broadly classified into four categories. The first type of approach focuses on solving the problems via estimating the optimal execution boundary, such as the method proposed in [29]. The second type of approach is mainly converting the problem into a linear problem via designing a semi-implicit scheme or implicit–explicit scheme to solve the problem, such as the L-stable method [19], the linearly implicit predictor–corrector scheme [30], or the IMEX-BDF method [31]. The third type of approach is mainly transforming the problem into a nonlinear equation to solve with iterative methods, such as the preconditioned penalty method [28], the fix-point method [32], the fitted finite volume iterative method [33], the modulus-based matrix splitting iteration method [34], and the modified modulus-based matrix splitting iteration method [35]. The fourth type of approach is designing a iterative method to solve the LCP directly, such as the projected SOR method [36] and the projected algebraic multigrid method [37]. In this paper, we design the fast numerical method following the third type of method because this kind of method can be extended easily to handle other related option pricing problems. However, it needs to solve a nonlinear equation, which requires quite a lot of computational cost. In [38], the authors compare several iterative procedures for solving the nonlinear scheme and find that the performance of the fix-point policy iteration with direct control formulation is a reliable general-purpose method. Inspired by [38], the fast numerical method designed in this paper is based on the policy iteration and Krylov subspace method, which can be expected to be reliable and efficient.
The rest of this article is as follows. In Section 2, a nonlinear finite difference scheme is proposed to discretize the considered problem with stability analysis. In Section 3, a preconditioned policy-Krylov subspace method has been developed with convergence analysis. Numerical experiments are given in Section 4, and conclusions are provided in Section 5.

2. Numerical Scheme and Theoretical Analysis

In this paper, the multi-state time-fractional LCP generated from a regime-switching modified Black–Scholes model with a time-fractional derivative [24,39] for the non-dividend American put option pricing problem is considered. The considered problem can be written as an optimization problem containing N s LCPs coupled with a Markov generator matrix. For  j = 1 , 2 , , N s , the j-th LCP can be written as follows:
L j V j ( S , t ) β j V j ( S , t ) t β j 0 , V j ( S , t ) V * ( S ) , L j V j ( S , t ) β j V j ( S , t ) t β j V j ( S , t ) V * ( S ) = 0 , V j ( S min , t ) = max { K S min , 0 } , t [ 0 , T ) , V j ( S max , t ) = max { K S max , 0 } , t [ 0 , T ) , V j ( S , T ) = max { K S , 0 } , S ( S min , S max ) , V * ( S ) = max { K S , 0 } , S ( S min , S max ) ,
where
L j V j ( S , t ) : = 1 2 σ j 2 S 2 2 V j ( S , t ) S 2 r S V j ( S , t ) S + r V j ( S , t ) i = 1 N s q j , i V i ( S , t ) ,
and β j V j ( S , t ) t β j is the modified right Riemann–Liouville fractional derivative [16] and defined by
β j V j ( S , t ) t β j = 1 Γ ( 1 β j ) t t T V j ( S , ξ ) V j ( S , T ) ( ξ t ) β j d ξ .
In the above model, K is the strike price of the options, and r is the risk-free interest rate. β j is the order of the time-fractional derivative of the j-th state and satisfies 0 < β j < 1 . σ j denotes the volatility when the underlying price is in the j-th regime. S is the asset’s price, and its lower and upper truncated boundaries are S min and S max , respectively. The entry q j , i is the transition rate of the Markov chain.
Inspired by [38], we attempted to devise a fast solution strategy based on policy iteration and direct control formulation. Therefore, the considered LCPs (1) are solved by transforming them into the Hamilton–Jacobi–Bellman (HJB) equations. By defining a new variable U j ( S , t ) = V j ( S , T t ) , the considered LCP (1) can be converted into the following HJB equation:
min { U j ( S , t ) V * ( S ) , 0 C D t β j U j ( S , t ) + L j U j ( S , t ) } = 0 ,
where 0 C D t β j U j ( S , t ) is the Caputo fractional derivative [40], which is defined by
0 C D t β j U j ( S , t ) = 1 Γ ( 1 β j ) 0 t U j ( S , ξ ) ( t ξ ) β j d ξ .
Remark that the transformation between the right-modified Riemann–Liouville fractional derivatives and the Caputo fractional derivatives has been given in [18]. The boundary conditions of HJB Equation (3) are listed by
U j ( S min , t ) = max { K S min , 0 } , U j ( S max , t ) = max { K S max , 0 } .

2.1. Finite Difference Method

A nonlinear finite difference scheme is developed to solve the coupled HJB equation given in (3). As the time-fractional derivative has been transferred into the Caputo fractional derivative, we can use the L 1 approximate [41] to discretize it, which is given by
0 C D t β j U j ( S , t m ) = τ β j Γ ( 2 β j ) [ a 0 ( β j ) U j ( S , t m ) s = 1 m 1 ( a m s 1 ( β j ) a m s ( β j ) ) U k ( S , t s ) a m 1 ( β j ) U j ( S , t 0 ) ] + O ( τ 2 β j ) ,
where
a l ( β j ) = ( l + 1 ) 1 β j l 1 β j , for l 0 .
Furthermore, the central difference method and backward finite difference method are used to discretize the second-order derivative and first-order derivative, respectively, which can be written by
2 U j ( S , t ) S 2 = U j ( S + h , t ) 2 U j ( S , t ) + U j ( S h , t ) h 2 + O ( h 2 ) ,
and
U j ( S , t ) S = U j ( S + h , t ) U j ( S , t ) h + O ( h ) .
Let N and M be the positive integers. The intervals [ S min , S max ] and [ 0 , T ] can be divided into N + 1 and M sub-intervals, respectively. Then, the numerical mesh can be defined by
S i = S min + i h , for i = 0 , 1 , , N + 1 , t m = m τ , for m = 0 , 1 , , M ,
where h = S max S min N + 1 , τ = T / M .
With the finite difference approximate given in (4), (5), and (6), by denoting u j , i ( m ) as the numerical solution of U j ( S i , t m ) , the HJB equation can be discretized. For  U j ( S , t ) C ( 4 , 2 ) ( [ S min , S max ] × [ 0 , T ] ) , and  j = 1 , 2 , , N s , the numerical scheme for HJB Equation (3) can be written by
min u j , i ( m ) v i * , τ β j Γ ( 2 β j ) [ a 0 ( β j ) u j , i ( m ) s = 1 m 1 ( a m s 1 ( β j ) a m s ( β j ) ) u j , i ( s ) a m 1 ( β j ) u j , i ( 0 ) ] + L ˜ u j , i ( m ) = 0 ,
where
L ˜ u j , i ( m ) = 1 2 σ j 2 S i 2 u j , i + 1 ( m ) 2 u j , i ( m ) + u j , i 1 ( m ) h 2 r S i u j , i + 1 ( m ) u j , i ( m ) h + r u j , i ( m ) s = 1 N s q j , s u j , s ( m ) ,
and v i * = V * ( S i ) .

2.2. Matrix Form

To simplify the notation of the matrix form, we define T n ( a , b , c ) as an n-by-n tri-diagonal matrix, which is given by
T n ( a , b , c ) = b c a b c a b .
Then, define the vectors u ( m ) = [ u 1 , 1 ( m ) , u 1 , 2 ( m ) , , u 1 , N ( m ) , u 2 , 1 ( m ) , u 2 , 2 ( m ) , , u N s , N ( m ) ] T , v ˜ * = L ^ v * , where v * = [ v 1 * , v 2 * , , v N * ] T and L ^ is an N s -by-1 column vector with all entries of 1.
With the above notations, the finite difference scheme (7) can be written by the following matrix form:
min { u ( m ) v ˜ * , W u ( m ) h ( m ) + D 1 I N f } = 0 ,
where
W = W ˜ D 1 Q I N , W ˜ = I N S N + r D 1 I N r h D 1 D s 1 A 1 + D 2 D s 2 A 2 .
The entry of Q is q j , i , which is the generator of the Markov chain. Satisfy (1) q j , i 0 , if  i j and (2) q j , j = j i q j , i for each j = 1 , 2 , , N s . The other matrices in (9) are listed by
D s 1 = d i a g { S 1 , S 2 , , S N } , D s 2 = d i a g { S 1 2 , S 2 2 , , S N 2 } , D 1 = d i a g { τ β 1 Γ ( 2 β 1 ) , τ β 2 Γ ( 2 β 2 ) , , τ β N s Γ ( 2 β N s ) } , D 2 = d i a g { τ β 1 σ 1 2 Γ ( 2 β 1 ) 2 h 2 , τ β 2 σ 2 2 Γ ( 2 β 2 ) 2 h 2 , , τ β N s σ N s 2 Γ ( 2 β N s ) 2 h 2 } , A 1 = T N ( 0 , 1 , 1 ) , A 2 = T N ( 1 , 2 , 1 ) .
Additionally, vectors h ( m ) and f in (8) are defined as follows. Vector h ( m ) = [ h 1 , 1 ( m ) , h 1 , 2 ( m ) , , h 1 , N ( m ) , h 2 , 1 ( m ) , h 2 , 2 ( m ) , , h N s , N ( m ) ] T , whose entries are given by
h j , i ( m ) = s = 1 m 1 ( a m s 1 ( β j ) a m s ( β j ) ) u j , i ( s ) + a m 1 ( β j ) u j , i ( 0 ) ,
and f = [ f 1 , 1 , f 1 , 2 , , f 1 , N , f 2 , 1 , f 2 , 2 , , f N s , N ] T with entries
f j , i = 0 , i 1 , N σ j 2 S 1 2 2 h 2 U j ( S min , t m ) , i = 1 , σ j 2 S N 2 2 h 2 U j ( S max , t m ) r S N h U j ( S max , t m ) , i = N .

2.3. Stability Analysis

Although the finite difference method has been developed, its stability property is necessary. With the matrix form of the discrete HJB equation in (8), one lemma can be derived first.
Lemma 1.
For 0 < τ < 1 , the matrix W = [ w p , q ] in (9) is an M-matrix with
| w p , p | s = 1 , s p N N s | w p , s | + 1 + τ β max Γ ( 2 β ^ ) r ,
where β ^ = arg min β j { β 1 , β 2 , , β N s } { Γ ( 2 β j ) } .
Proof. 
By the definition of matrix Q, we know that the row sum of it equals to zero. With the definition of matrix W in (9), and the properties of the Kronnecker product, the matrix r h D 1 D s 1 A 1 + D 2 D s 2 A 2 D 1 Q I N is a weakly diagonal dominant matrix with the positive main diagonal and negative off diagonals. The remaining part of the coefficient matrix W is a positive diagonal matrix I N N s + r D 1 I N . Hence, the matrix W has the positive main diagonal and negative off diagonal with the following property: | w p , p | s = 1 , s p N N s | w p , s | 1 + τ β min Γ ( 2 β ^ ) r , where β ^ = arg min β j { β 1 , β 2 , , β N s } { Γ ( 2 β j ) } .    □
With Lemma 1, we have the following theorem, which ensures the unconditional stability of the nonlinear numerical scheme (7).
Theorem 1.
The nonlinear scheme (7) is unconditionally stable, i.e.,
u ( m ) u ( 0 ) + T β ˜ Γ ( 1 β max ) f ,
where β ˜ = arg max β j { β 1 , β 2 , , β N s } { T β j } .
Proof. 
Let u ( m ) be the solution of the discrete HJB Equation (8) on the m-th time level. Assume that u j ^ ( m ) , i ^ ( m ) ( m ) is an entry of the vector u ( m ) and satisfies
u j ^ ( m ) , i ^ ( m ) ( m ) = max 1 j N s , 1 i N | u j , i ( m ) | = u ( m ) .
Then, we discuss u ( m ) in two cases when the solution of (7) exists.
Case i: For u j ^ ( m ) , i ^ ( m ) ( m ) > v i ^ ( m ) * , we know that u j ^ ( m ) , i ^ ( m ) ( m ) is determined by ( W u ( m ) h ( m ) + D 1 I N f ) ( j ^ ( m ) 1 ) N + i ^ ( m ) = 0 , where ( v ) ( j ^ ( m ) 1 ) N + i ^ ( m ) ) denotes the ( j ^ ( m ) 1 ) N + i ^ ( m ) -th entry of any given vector v . Hence, with Lemma 1, we obtain the following inequalities:
u ( m ) = | u j ^ ( m ) , i ^ ( m ) ( m ) | w ( j ^ ( m ) 1 ) N + i ^ ( m ) , ( j ^ ( m ) 1 ) N + i ^ ( m ) u j ^ ( m ) , i ^ ( m ) ( m ) s ( j ^ ( m ) 1 ) N + i ^ ( m ) w ( j ^ ( m ) 1 ) N + i ^ ( m ) , s u j ^ ( m ) , i ^ ( m ) ( m ) w ( j ^ ( m ) 1 ) N + i ^ ( m ) , ( j ^ ( m ) 1 ) N + i ^ ( m ) u j ^ ( m ) , i ^ ( m ) ( m ) + s ( j ^ ( m ) 1 ) N + i ^ ( m ) w ( j ^ ( m ) 1 ) N + i ^ ( m ) , s u s ( m ) = h j ^ ( m ) , i ^ ( m ) ( m ) τ β j ^ ( m ) Γ ( 2 β j ^ ( m ) ) f j ^ ( m ) , i ^ ( m ) s = 1 m 1 ( a m s 1 ( β j ^ ( m ) ) a m s ( β j ^ ( m ) ) ) u j ^ ( m ) , i ^ ( m ) ( s ) + a m 1 ( β j ^ ( m ) ) u j ^ ( m ) , i ^ ( m ) ( 0 ) + τ β j ^ ( m ) Γ ( 2 β j ^ ( m ) ) f j ^ ( m ) , i ^ ( m ) s = 1 m 1 ( a m s 1 ( β j ^ ( m ) ) a m s ( β j ^ ( m ) ) ) u ( s ) + a m 1 ( β j ^ ( m ) ) u ( 0 ) + τ β j ^ ( m ) Γ ( 2 β j ^ ( m ) ) f .
Case ii: For u j ^ ( m ) , i ^ ( m ) ( m ) = v i ^ ( m ) * , it is straightforward to know that u ( m ) v ˜ * .
Combining two cases yields the following inequality:
u ( m ) = | u j ^ ( m ) , i ^ ( m ) ( m ) | max { v ˜ * , s = 1 m 1 ( a m s 1 ( β j ^ ( m ) ) a m s ( β j ^ ( m ) ) ) u ( s ) + a m 1 ( β j ^ ( m ) ) u ( 0 ) + τ β j ^ ( m ) Γ ( 2 β j ^ ( m ) ) f } .
When m = 1 , since u ( 0 ) = v ˜ * , we have
u ( 1 ) max v ˜ * , a 0 ( β j ^ ( 1 ) ) u ( 0 ) + τ β j ^ ( 1 ) Γ ( 2 β j ^ ( 1 ) ) f u ( 0 ) + t 1 β ˜ Γ ( 1 β max ) f .
When m = 2 , with (10), it holds that
u ( 2 ) max v ˜ * , ( a 0 ( β j ^ ( 2 ) ) a 1 ( β j ^ ( 2 ) ) ) u ( 1 ) + a 1 ( β j ^ ( 2 ) ) u ( 0 ) + τ β j ^ ( 2 ) Γ ( 2 β j ^ ( 2 ) ) f max v ˜ * , ( a 0 ( β j ^ ( 2 ) ) a 1 ( β j ^ ( 2 ) ) ) u ( 1 ) + a 1 ( β j ^ ( 2 ) ) u ( 0 ) + τ β j ^ ( 2 ) Γ ( 2 β j ^ ( 2 ) ) ( 1 β j ^ ( 2 ) ) 2 β j ^ ( 2 ) f u ( 0 ) + t 2 β ˜ Γ ( 1 β max ) f
Assume that the following formula holds when m = 1 , 2 , , l , which is given by
u ( m ) u ( 0 ) + t m β ˜ Γ ( 1 β max ) f .
Then, when m = l + 1 , with (10), it yields that
u ( l + 1 ) max { v ˜ * , s = 1 l ( a l s ( β j ^ ( l + 1 ) ) a l s + 1 ( β j ^ ( l + 1 ) ) ) u ( s ) + a l ( β j ^ ( l + 1 ) ) u ( 0 ) + τ β j ^ ( l + 1 ) Γ ( 2 β j ^ ( l + 1 ) ) f } max { v ˜ * , s = 1 l ( a l s ( β j ^ ( l + 1 ) ) a l s + 1 ( β j ^ ( l + 1 ) ) ) u ( 0 ) + t s β ˜ Γ ( 1 β max ) f + a l ( β j ^ ( l + 1 ) ) u ( 0 ) + t l + 1 β j ^ ( l + 1 ) Γ ( 1 β j ^ ( l + 1 ) ) f } u ( 0 ) + t l + 1 β ˜ Γ ( 1 β max ) f .
Hence, the nonlinear scheme (7) is unconditionally stable, i.e., 
u ( m ) u ( 0 ) + T β ˜ Γ ( 1 β max ) f .
   □

3. Preconditioned Policy-Krylov Subspace Method

In order to solve the HJB Equation (8), the policy iteration method [42] is used as a component in constructing the fast algorithm.

3.1. Policy Iteration Method

The discrete HJB equation can be written by the following component-wise form [42]:
θ i ( u ( m ) v ˜ * ) i + ( 1 θ i ) ( W u ( m ) h ( m ) + D 1 I N f ) i = 0 .
where ( v ) i denotes the i-th entry of any given vector v and 1 i N N s . The vector θ contains the control parameters and its i-th entry is valued by θ i = arg min θ { 0 , 1 } { θ ( u ( m ) v ˜ * ) i + ( 1 θ ) ( W u ( m ) h ( m ) + D 1 I N f ) i } . Then, the policy iteration method can be used to solve (11). The framework of the policy iteration method is presented in Algorithm 1.
Algorithm 1 Framework of the policy iteration method
Let ( u ( m ) ) k be the k-th iteration of the policy iteration method for computing the solution u ( m ) . Let θ k R N N s , M k R N N s × N N s , and b k R N N s . For  1 i N N s , we have
( θ ( m ) ) i k = arg min θ { 0 , 1 } { θ ( ( u ( m ) ) k v ˜ * ) i + ( 1 θ ) ( W u ( m ) h ( m ) + D 1 I N f ) i } , ( M ( m ) ) i k = ( θ ( m ) ) i k ( I N N s ) i + ( 1 ( θ ( m ) ) i k ) ( W ) i , ( b ( m ) ) i k = ( θ ( m ) ) i k ( v * ) i + ( 1 ( θ ( m ) ) i k ) ( h ( m ) D 1 I N f ) i .
Find ( u ( m ) ) k + 1 R N N s , such that
( M ( m ) ) k ( u ( m ) ) k + 1 = ( b ( m ) ) k .
With Algorithm 1, it is known that a linear system needs to be solved in each policy iteration. To ensure the algorithm is reliable, the coefficient matrix ( M ( m ) ) k in (12) should be invertible. Hence, we have the following lemma. Remark that the notation M m , k is used to stand for ( M ( m ) ) k to simplify the statement in the subsequent theoretical analysis.
Lemma 2.
For 0 < β j < 1 , the coefficient matrix M m , k is invertible and satisfies
( M m , k ) 1 1 .
Proof. 
With Lemma 1, we know that matrix W is an M-matrix. Coefficient matrix M m , k is constructed by selecting rows from W and the identity matrix I N N s , thus matrix M m , k is also an M-matrix and invertible. In addition, the entries [ i j , k ] of the identity matrix hold that | i k , k | j k | i j , k | = 1 . Refer to [43], ( M m , k ) 1 1 can be obtained. □
The matrix W is an M-matrix, then the properties of the LCP considered in this paper are similar to the case in [42]. Refer to the proof in [42], the following theorem can be obtained straightforwardly, which ensures that the policy iteration method converges to the solution in finite steps.
Theorem 2.
Let ( u k ) k = 0 be the sequence generated by Algorithm 1. Then,
u k + 1 u k , n N a n d u n = u * , n K
where u * is the solution of (11) and K is an integer independent of u 0 .

3.2. Preconditioned Krylov Subspace Method

In each policy iteration, a linear system (12) is needed to be solved. Based on the structure of the matrix M m , k , a fast Krylov subspace method is used as the linear solver. The matrix-vector product M m , k v can be obtained by
M m , k v = Θ m , k v + ( I N N s Θ m , k ) ( W ˜ v vec ( ( D 1 Q V T ) T ) ) ,
where v is any given vector, and V is an N-by- N s matrix reshaped by v column by column. Θ m , k is the diagonal matrix, and its i-th diagonal entry equals the control parameter ( θ ( m ) ) i k in Algorithm 1. Since matrix W ˜ is a block diagonal matrix with tri-diagonal blocks, with (13), the operation cost of M m , k v is O ( N N s 2 ) , and each iteration of the Krylov subspace method also costs O ( N N s 2 ) .
Although the operation cost of each iteration of the Krylov subspace method is O ( N N s 2 ) , the total cost for solving the linear system (12) depends on the iteration number. If the coefficient matrix M m , k is ill conditioned, the Krylov subspace method converges to the solution slowly. Hence, the preconditioning technique is considered to reduce the iteration number and avoid this situation. However, extra computational costs are needed for the preconditioning technique in each iteration. The preconditioner based on circulant matrices will not be considered because its extra computational cost usually is O ( N N s log N N s ) due to the fast Fourier transformation. In this section, we develop a simple preconditioner with the tri-diagonal structure, which is defined by
P m , k = Θ m , k I N s N + ( I N N s Θ m , k ) ( W ˜ D 1 Q d I N ) ,
where Q d = d i a g { q 1 , 1 , q 2 , 2 , , q N s , N s } . Because the preconditioner P m , k keeps the tri-diagonal structure, the matrix-vector product ( P m , k ) 1 v costs O ( N N s ) based on the Thomas algorithm, where v is any given vector. Although the operation cost for the preconditioned Krylov subspace method is still O ( N N s 2 ) per iteration, some theorems are needed to ensure the proposed preconditioning technique is useful. First, we need to ensure the preconditioner is invertible.
Theorem 3.
The preconditioner P m , k is invertible and satisfies
( P m , k ) 1 1 .
Proof. 
By the definition of the preconditioner P m , k , we know that it is constructed by selecting the row between the matrices W ˜ D 1 Q d I N and I N N s . Refer to the proof of Lemmas 1 and 2, it is easy to know that P m , k is an M-matrix and satisfies ( P m , k ) 1 1 . □
Despite the invertibility of the proposed preconditioner, proved by Theorem 3, the convergence rate of the preconditioned Krylov subspace method is critical. It is well known that the convergence rate of the Krylov subspace method is influenced by the condition number of the coefficient matrix. Therefore, the following theorem is given to prove that the condition number of the coefficient matrix with the proposed preconditioner has an upper boundary.
Theorem 4.
For 0 < β j < 1 and 0 < τ < 1 , the condition number of the preconditioned coefficient matrix ( P m , k ) 1 M m , k is bounded, that is,
c o n d ( ( P m , k ) 1 M m , k ) τ β min max k N s { | q k , k | } + 1 2 .
Proof. 
With Theorem 3, it is straightforward that
( P m , k ) 1 M m , k = ( P m , k ) 1 ( M m , k P m , k ) + I N N s ( P m , k ) 1 ( I N N s Θ ) D 1 ( Q d Q ) I N + 1 τ β min max k N s { | q k , k | } + 1 .
Additionally, with Lemma 2, it yields
( ( P m , k ) 1 M m , k ) 1 = ( M m , k ) 1 P m , k = ( M m , k ) 1 ( P m , k M m , k ) + I N N s ( M m , k ) 1 ( P m , k M m , k ) + 1 ( M m , k ) 1 D 1 ( Q Q d ) I N + 1 τ β min max k N s { | q k , k | } + 1 .
Hence, the condition number of the preconditioned matrix ( P m , k ) 1 M m , k is bounded, which can be written by
c o n d ( ( P m , k ) 1 M m , k ) = ( P m , k ) 1 M m , k ( ( P m , k ) 1 M m , k ) 1 τ β min max k N s { | q k , k | } + 1 2 .
With Theorem 4, we know that the proposed preconditioner is able to ensure that the condition number of the preconditioned matrix is bounded. The upper bound is affected by Q, but τ β min usually being sufficiently small can help us to ensure that the upper bound is not too large. Hence, the proposed preconditioner can be expected to reduce the condition of the coefficient matrix and improve the convergence rate of the Krylov subspace method.
Remark 1.
Although Theorem 4 holds for the left preconditioner, this theorem still holds for the right preconditioner with similar proofs. The numerical tests with the right preconditioning technique are also given in Section 4.

4. Numerical Experiment

In this section, three American option pricing problems with time-fractional derivatives containing two states, four states, and eight states are considered as the numerical examples to test our proposed preconditioned policy-Krylov subspace method. All the numerical experiments were carried out in Matlab 2016a in the workstation with the following configuration: Intel(R) Core(TM) i9-10900K CPU 3.70 GHz and 64 GB RAM. The bi-conjugate gradient stabilized method and the GMRES method are used in this section, and their stopping criterion is 10 10 . The stopping criterion of the policy iteration method is 10 7 . The initial guesses of two iterative methods are both the numerical solution on the previous time step. The restart of the GMRES method is 20 in this section. In addition, to facilitate the calculation, in this section, we set a new symbol N ˜ that satisfies N ˜ = N + 1 .
In the numerical experiment, we assume that the strike price K = 50 , the risk-free interest rate r = 0.05 , and the expiry moment T = 1 , i.e., one year. The truncation boundaries [ S min , S max ] = [ 0 , 100 ] . The other parameters of the three cases are listed by
(a)
N s = 2 , σ = [ 0.2 , 0.4 ] , β = [ 0.8 , 0.95 ] , Q = 2 2 3 3 .
(b)
N s = 4 , σ = [ 0.2 , 0.3 , 0.15 , 0.2 ] , β = [ 0.8 , 0.95 , 0.6 , 0.7 ] ,
Q = 45 20 15 10 30 80 30 20 10 10 60 40 20 10 20 50 .
(c)
N s = 8 , σ = [ 0.2 , 0.3 , 0.15 , 0.2 , 0.5 , 0.2 , 0.4 , 0.3 ] , β = [ 0.8 , 0.95 , 0.6 , 0.7 , 0.8 , 0.3 , 0.2 , 0.6 ] ,
Q = 192 34 16 36 31 23 25 27 32 125 3 37 15 2 29 7 17 21 144 19 9 9 25 14 36 5 4 118 16 14 18 25 7 34 7 13 145 32 21 31 10 24 9 36 5 98 11 3 5 14 16 14 37 1 124 37 5 20 1 4 38 6 7 81 .
The first part of the numerical experiment is to test the convergence rate of the nonlinear scheme (7) and the effectiveness of the proposed solution strategy. The performance of the Krylov subspace method without preconditioning technique is also presented in order to demonstrate the performance of the proposed preconditioner P m , k . Since the exact solution is absent, the numerical solution from a fine grid where N ˜ = 2 13 and M = 2 11 is regarded as the exact solution to compute the truncation errors. Related numerical results from the bi-conjugate gradient stabilized method are presented in Table 1, while the numerical results from the GMRES method are given in Table 2.
In Table 1 and Table 2, the notation `Error’ denotes the infinity norm of the truncation error of the nonlinear scheme (7) solved by the preconditioned policy-Krylov subspace method, and `Rate’ means the convergence rate. The notations `CG’ and `PCG’ stand for the bi-conjugate gradient stabilized method without and with the proposed preconditioner, respectively. Similarly, `GMRES’ and `PGMRES’ denote the GMRES method and the GMRES method with the proposed preconditioner, respectively. The iteration numbers of the policy iteration method with different linear solvers are almost the same. Thus only the average number of policy iterations from the preconditioned Krylov subspace method is given and denoted by `Ite-Out’. The notations `Ite-In’ and `Time’ are the average iteration numbers of the linear solvers and CPU time, respectively. Remark that the left preconditioning technique is used for the preconditioned bi-conjugate gradient stabilized method, and the right preconditioning technique is used for the preconditioned GMRES method.
From Table 1 and Table 2, it can be seen that the convergence rate of the proposed nonlinear scheme is first-order. Note that the errors of the two tables are close because both of these linear solvers solve the same problem with the same stopping criterion. The policy iteration method converges to the solution rapidly, which can be found by `Ite-Out’. Compared with the bi-conjugate gradient stabilized method and GMRES method, we find that the GMRES method takes more iterations for the same problem in our numerical examples. Significantly, the proposed preconditioner can improve the convergence rate of the bi-conjugate gradient stabilized and GMRES methods. As analyzed in Theorem 4, the upper bound on the condition number of Case (c) is higher than Case (a) and (b) due to the number of states being larger. Hence, the bi-conjugate gradient stabilized method and the GMRES method with preconditioner needs more iterations to converge.
To further analyze the effect of the preconditioning technique and verify Theorem 4, the second part of the numerical experiment focuses on testing the eigenvalues and condition numbers of the preconditioned matrix ( P m , k ) 1 M m , k and coefficient matrix M m , k . It is difficult to choose a specific matrix to finish the test since the matrices P m , k and M m , k are determined by m and k. Given the existence of the American option free boundary, there is usually a point θ * such that the matrix satisfies θ = 1 when S < θ * and θ = 0 when S θ * . To facilitate further analysis, we used Case (b) to finish the test and set a particular P and M that satisfied θ * = N ˜ 2 . Then, the eigenvalue distributions of the coefficient matrix M and preconditioned matrix P 1 M when M = N ˜ = 2 10 are given in Figure 1. Additionally, the condition numbers of them when M = N ˜ are given in Figure 2.
From Figure 1a, we know that the maximum of the real part of the eigenvalues of the coefficient matrix is more than 600, which may cause the slow convergence rate of the Krylov subspace method. Conversely, the eigenvalues in Figure 1b are clustered around one after preconditioning, although the gathering radius is close to 0.3 , which helps improve the convergence rate of the Krylov subspace method.
In Figure 2, the condition number of the coefficient matrix M, the condition number of the preconditioned matrix P 1 M , and the upper boundary of the condition number given by Theorem 4 are presented. From Figure 2a, the condition number is lower than 5 with the preconditioning technique. The theoretical upper boundary decreases from 28.65 to 3.33 as M increases. In Figure 2b, the condition number without the preconditioning technique is given to compare with the condition number of P 1 M . One can see that the condition number is increasing with M and the maximum value is 3318.69 when M = 2 11 . With the proposed preconditioning technique, the condition number is kept below 5, which demonstrates that the proposed preconditioner can reduce the condition number of the coefficient matrix M.
In the last part of the numerical experiment, the option values of three cases are presented in Figure 3 when M = N ˜ = 2 12 , respectively. In Figure 3, we can find that multi-regimes cause the option value differences between different regimes.

5. Concluding Remark

This paper considers the multi-state time-fractional LCP for American options pricing and proposes a preconditioned policy-Krylov subspace method to handle this problem. The considered problem has been transformed into HJB equations coupled by the Markov generator matrix, then a nonlinear finite difference scheme is developed to discretize it with unconditional stability guarantees. A preconditioned policy-Krylov subspace method has been proposed to solve the discrete HJB equation with convergence guarantees. The condition number of the preconditioned coefficient matrix has been proven to have an upper boundary. The numerical experiments are presented to demonstrate the efficiency of the proposed nonlinear finite difference scheme and the preconditioned policy-Krylov subspace method.
In addition, the considered problem is based on a non-dividends option pricing model, and the dividend models are yet to be addressed. The problem of pricing American options with dividends, as well as the problem of pricing options with fractional order derivatives in space-time, will be addressed in our future work.

Author Contributions

Conceptualization, S.-L.L.; Methodology, X.C. and X.G.; Software, X.G.; Formal analysis, X.C.; Writing—original draft, X.C.; Writing—review & editing, X.G., S.-L.L. and Y.S.; Supervision, S.-L.L. All authors have read and agreed to the published version of the manuscript.

Funding

The first author is supported by Guangdong Basic and Applied Basic Foundation (Grant No. 2020A1515110991) and the National Natural Science Foundation of China (Grant No. 12101137). The third author is supported by the University of Macau (Grant No. MYRG2020-00208-FST and MYRG2022-00262-FST). The fourth author is supported by Guangdong Basic and Applied Basic Foundation (Grant No. 2022A1515011125) and the National Natural Science Foundation of China (Grant No. 72271064).

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the reviewers for their valuable comments and suggestions that have improved the quality of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Black, F.; Scholes, M. The pricing of options and corporate liabilities. J. Politics Econ. 1973, 81, 637–654. [Google Scholar] [CrossRef]
  2. Stein, E.M.; Stein, J.C. Stock price distributions with stochastic volatility:an analytic approach. Rev. Financ. Stud. 1991, 4, 727–752. [Google Scholar] [CrossRef]
  3. Heston, S.L. A closed-form solution for options with stochastic volatility with applications to bond and currency options. Rev. Financ. Stud. 1993, 6, 327–343. [Google Scholar] [CrossRef]
  4. Kou, S.G. A jump-diffusion model for option pricing. Manag. Sci. 2002, 48, 1086–1101. [Google Scholar] [CrossRef]
  5. Madan, D.B.; Carr, P.P.; Chang, E.C. The variance gamma process and option pricing. Eur. Financ. Rev. 1998, 2, 79–105. [Google Scholar] [CrossRef]
  6. Boyarchenko, S.; Levendorskii, S. Non-Gaussian Merton-Black-Scholes Theory; World Scientific: Singapore, 2002. [Google Scholar]
  7. Fulop, A.; Li, J.; Yu, J. Self-exciting jumps, learning, and asset pricing implications. Rev. Financ. Stud. 2015, 28, 876–912. [Google Scholar] [CrossRef]
  8. Hawkes, A.G. Hawkes jump-diffusions and finance: A brief history and review. Eur. J. Financ. 2022, 28, 627–641. [Google Scholar] [CrossRef]
  9. Cartea, A.; del Castillo-Negrete, D. Fractional diffusion models of option prices in markets with jumps. Physical A 2007, 374, 749–763. [Google Scholar] [CrossRef]
  10. Carr, P.; Wu, L. The finite moment log stable process and option pricing. J. Financ. 2003, 58, 753–777. [Google Scholar] [CrossRef]
  11. Carr, P.; Geman, H.; Madan, D.B.; Yor, M. The fine structure of asset returns: An empirical investigation. J. Bus. 2002, 75, 305–333. [Google Scholar] [CrossRef]
  12. Chen, W.; Wang, S. A penalty method for a fractional order parabolic variational inequality governing American put option valuation. Comput. Math. Appl. 2014, 67, 77–90. [Google Scholar] [CrossRef]
  13. Meng, Q.; Ding, D.; Sheng, Q. Preconditioned iterative methods for fractional diffusion models in finance. Numer. Meth. Part. Differ. Equ. 2014, 31, 1382–1395. [Google Scholar] [CrossRef]
  14. Zhang, H.; Liu, F.; Turner, I.; Chen, S. The numerical simulation of the tempered fractional Black–Scholes equation for European double barrier option. Appl. Math. Model. 2016, 40, 5819–5834. [Google Scholar] [CrossRef]
  15. Wyss, W. The fractional Black-Scholes equation. Fract. Calc. Appl. Anal. 2000, 3, 51–61. [Google Scholar]
  16. Jumarie, G. Derivation and solutions of some fractional Black–Scholes equations in coarse-grained space and time. Application to Merton’s optimal portfolio. Comput. Math. Appl. 2010, 59, 1142–1164. [Google Scholar] [CrossRef]
  17. Cartea, A. Derivation pricing with marked point processes using tick-by-tick data. Quant. Financ. 2013, 13, 111–123. [Google Scholar] [CrossRef]
  18. Zhang, H.; Liu, F.; Chen, S.; Anh, V.; Chen, J. Fast numerical simulation of a new time-space fractional option pricing model governing European call option. Appl. Math. Comput. 2018, 339, 186–198. [Google Scholar] [CrossRef]
  19. Yousur, M.; Khaliq, A.Q.M.; Liu, R.H. Pricing American options under multi-state regime switching with an efficient L-stable method. Int. J. Comput. Math. 2015, 92, 2530–2550. [Google Scholar]
  20. Hamilton, J.D. Analysis of time series subject to changes in regime. J. Econ. 1990, 45, 39–70. [Google Scholar] [CrossRef]
  21. He, X.; Chen, W. A Monte-Carlo based approach for pricing credit default swaps with regime switching. Comput. Math. Appl. 2018, 76, 1758–1766. [Google Scholar] [CrossRef]
  22. Koleva, M.N.; Vulkov, L.G. Fast positivity preserving numerical method for time-fractional regime-switching option pricing problem. In Advanced Computing in Industrial Mathematics; Georgiev, I., Kostadinov, H., Lilkova, E., Eds.; Springer: Cham, Switzerland, 2023; pp. 88–99. [Google Scholar]
  23. Lin, S.; He, X. A regime switching fractional Black–Scholes model and European option pricing. Commun. Nonlinear Sci. Numer. Simul. 2020, 85, 105222. [Google Scholar] [CrossRef]
  24. Saberi, E.; Hejazi, S.R.; Dastranj, E. A new method for option pricing via time-fractional PDE. Asian-Eur. J. Math. 2018, 11, 1850074. [Google Scholar] [CrossRef]
  25. Zhou, Z.; Ma, J.; Gao, X. Convergence of iterative Laplace transform methods for a system of fractional PDEs and PIDEs arising in option pricing. East Asian J. Appl. Math. 2018, 8, 782–808. [Google Scholar] [CrossRef]
  26. Chen, X.; Deng, D.; Lei, S.; Wang, W. An implicit-explicit preconditioned direct method for pricing options under regime-switching tempered fractional partial differential models. Numer. Algorithms 2021, 87, 939–965. [Google Scholar] [CrossRef]
  27. Wang, M.; Wang, C.; Yin, J. A second-order ADI method for pricing options under fractional regime-switching models. Netw. Heterog. Media 2023, 18, 647–663. [Google Scholar] [CrossRef]
  28. Lei, S.; Wang, W.; Chen, X.; Ding, D. A fast preconditioned penalty method for American options pricing under regime-switching tempered fractional diffusion models. J. Sci. Comput. 2018, 75, 1633–1655. [Google Scholar] [CrossRef]
  29. Boyarchenko, S.; Levendorskii, S. American options in regime-switching models. SIAM J. Control Optim. 2009, 48, 1353–1376. [Google Scholar] [CrossRef]
  30. Khaliq, A.K.Q.; 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]
  31. Wang, W.; Chen, Y.; Fang, H. On the variable two-step IMEX BDF method for parabolic integro-differential equations with nonsmooth initial data arising in finance. SIAM J. Numer. Anal. 2019, 57, 1289–1317. [Google Scholar] [CrossRef]
  32. Shi, X.; Yang, L.; Huang, Z. A fixed point method for the linear complementarity problem arising from American option pricing. Acta Math. Appl. Sin. Engl. Ser. 2016, 32, 921–932. [Google Scholar] [CrossRef]
  33. Gan, X.; Yin, J. Pricing American options under regime-switching model with a Crank-Nicolson fitted finite volume method. East Asian Appl. Math. 2020, 10, 499–519. [Google Scholar]
  34. Bai, Z. Modulus-based matrix splitting iteration methods for linear complementarity problems. Numer. Linear Algebra Appl. 2010, 17, 917–933. [Google Scholar] [CrossRef]
  35. Zheng, H.; Vong, S. A modified modulus-based matrix splitting iteration method for solving implicit complementarity problems. Numer. Algorithms 2019, 82, 573–592. [Google Scholar] [CrossRef]
  36. Cryer, C.W. The solution of a quadratic programming problem using systematic overrelaxation. SIAM J. Control Optim. 1971, 9, 385–392. [Google Scholar] [CrossRef]
  37. Toivanen, J.; Oosterlee, C.W. A projected algebraic multigrid method for linear complementarity problems. Numer. Math. Theor. Meth. Appl. 2012, 5, 85–98. [Google Scholar] [CrossRef]
  38. Huang, Y.; Forsyth, P.A.; Labahn, G. Methods for pricing American options under regime switching. SIAM J. Sci. Comput. 2011, 33, 2144–2168. [Google Scholar] [CrossRef]
  39. Chen, C.; Wang, Z.; Yang, Y. A new operator splitting method for American options under fractional Black-Scholes models. Comput. Math. Appl. 2019, 77, 2130–2144. [Google Scholar] [CrossRef]
  40. Podlubny, I. Fractional Differential Equations; Academic Press: New York, NY, USA, 1999. [Google Scholar]
  41. Gao, G.; Sun, Z.; Zhang, H. A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications. J. Comput. Phys. 2014, 259, 33–50. [Google Scholar] [CrossRef]
  42. Reisinger, C.; Witte, J.H. On the use of policy iteration as an easy way of pricing American options. SIAM J. Financ. Math. 2012, 3, 459–478. [Google Scholar] [CrossRef]
  43. Varah, J.M. A lower bound for the smallest singular value of a matrix. Linear Algebra Appl. 1975, 11, 3–5. [Google Scholar] [CrossRef]
Figure 1. Eigenvalues distribution: (a) eigenvalue of coefficient matrix M; (b) eigenvalue of preconditioned matrix P 1 M .
Figure 1. Eigenvalues distribution: (a) eigenvalue of coefficient matrix M; (b) eigenvalue of preconditioned matrix P 1 M .
Fractalfract 07 00334 g001
Figure 2. Condition number of the coefficient matrix: (a) condition number of preconditioned matrix P 1 M ; (b) comparison of condition number of P 1 M and M.
Figure 2. Condition number of the coefficient matrix: (a) condition number of preconditioned matrix P 1 M ; (b) comparison of condition number of P 1 M and M.
Fractalfract 07 00334 g002
Figure 3. Option value of three examples: (a) option value of two regimes; (b) option value of four regimes; (c) option value of eight regimes.
Figure 3. Option value of three examples: (a) option value of two regimes; (b) option value of four regimes; (c) option value of eight regimes.
Fractalfract 07 00334 g003
Table 1. Numerical results from the bi-conjugate gradient stabilized methods.
Table 1. Numerical results from the bi-conjugate gradient stabilized methods.
CGPCG
NMErrorRateIte-OutIte-InTimeIte-InTime
Case (a)
2 8 2 6 3.9728 × 10 3 -2.8677.460.37 s2.560.08 s
2 9 2 7 1.9224 × 10 3 1.04722.82106.051.24 s2.080.21 s
2 10 2 8 9.1278 × 10 4 1.07462.86141.444.96 s1.780.79 s
2 11 2 9 3.9949 × 10 4 1.19212.86184.5523.49 s1.544.40 s
Case (b)
2 8 2 6 8.9047 × 10 3 -2.7362.250.38 s7.920.19 s
2 9 2 7 4.5215 × 10 3 0.97782.7791.381.59 s5.790.51 s
2 10 2 8 2.1635 × 10 3 1.06342.71138.528.79 s4.612.26 s
2 11 2 9 9.3962 × 10 4 1.20322.69208.5744.75 s3.709.30 s
Case (c)
2 8 2 6 1.9890 × 10 3 -3.33427.404.83 s14.840.55 s
2 9 2 7 9.8072 × 10 4 1.02013.32848.7335.33 s11.462.04 s
2 10 2 8 4.8702 × 10 4 1.00993.261705.16242.65 s9.437.15 s
2 11 2 9 2.1722 × 10 4 1.16483.24--7.4130.50 s
Table 2. Numerical results from the GMRES methods.
Table 2. Numerical results from the GMRES methods.
GMRESPGMRES
NMErrorRateIte-OutIte-InTimeIte-InTime
Case (a)
2 8 2 6 3.9724 × 10 3 -2.86156.681.66 s5.040.11 s
2 9 2 7 1.9220 × 10 3 1.04742.82242.876.81 s4.120.26 s
2 10 2 8 9.1239 × 10 4 1.07492.86406.8636.64 s3.330.90 s
2 11 2 9 3.9910 × 10 4 1.19292.86692.43392.76 s3.055.00 s
Case (b)
2 8 2 6 8.9046 × 10 3 -2.73113.891.58 s13.860.30 s
2 9 2 7 4.5214 × 10 3 0.97782.77191.738.20 s10.590.72 s
2 10 2 8 2.1634 × 10 3 1.06352.71372.0499.65 s8.293.07 s
2 11 2 9 9.3955 × 10 4 1.20322.69760.57637.37 s6.5511.78 s
Case (c)
2 8 2 6 1.9890 × 10 3 -3.331411.1136.60 s24.770.98 s
2 9 2 7 9.8069 × 10 4 1.02023.328428.321409.52 s18.564.21 s
2 10 2 8 4.8698 × 10 4 1.00993.2632,658.1217,206.00 s14.5610.57 s
2 11 2 9 2.1718 × 10 4 1.16503.24--11.5539.96 s
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

Chen, X.; Gong, X.; Lei, S.-L.; Sun, Y. A Preconditioned Iterative Method for a Multi-State Time-Fractional Linear Complementary Problem in Option Pricing. Fractal Fract. 2023, 7, 334. https://doi.org/10.3390/fractalfract7040334

AMA Style

Chen X, Gong X, Lei S-L, Sun Y. A Preconditioned Iterative Method for a Multi-State Time-Fractional Linear Complementary Problem in Option Pricing. Fractal and Fractional. 2023; 7(4):334. https://doi.org/10.3390/fractalfract7040334

Chicago/Turabian Style

Chen, Xu, Xinxin Gong, Siu-Long Lei, and Youfa Sun. 2023. "A Preconditioned Iterative Method for a Multi-State Time-Fractional Linear Complementary Problem in Option Pricing" Fractal and Fractional 7, no. 4: 334. https://doi.org/10.3390/fractalfract7040334

APA Style

Chen, X., Gong, X., Lei, S. -L., & Sun, Y. (2023). A Preconditioned Iterative Method for a Multi-State Time-Fractional Linear Complementary Problem in Option Pricing. Fractal and Fractional, 7(4), 334. https://doi.org/10.3390/fractalfract7040334

Article Metrics

Back to TopTop