A Preconditioned Policy–Krylov Subspace Method for Fractional Partial Integro-Differential HJB Equations in Finance

: To better simulate the prices of underlying assets and improve the accuracy of pricing financial derivatives, an increasing number of new models are being proposed. Among them, the Lévy process with jumps has received increasing attention because of its capacity to model sudden movements in asset prices. This paper explores the Hamilton–Jacobi–Bellman (HJB) equation with a fractional derivative and an integro-differential operator, which arise in the valuation of American options and stock loans based on the Lévy-α -stable process with jumps model. We design a fast solution strategy that includes the policy iteration method, Krylov subspace method, and banded preconditioner, aiming to solve this equation rapidly. To solve the resulting HJB equation, a finite difference method including an upwind scheme, shifted Grünwald approximation, and trapezoidal method is developed with stability and convergence analysis. Then, an algorithmic framework involving the policy iteration method and the Krylov subspace method is employed. To improve the performance of the above solver, a banded preconditioner is proposed with condition number analysis. Finally, two examples, sugar option pricing and stock loan valuation, are provided to illustrate the effectiveness of the considered model and the efficiency of the proposed preconditioned policy–Krylov subspace method.


Introduction
The option-pricing problem has always been a hot research topic in finance, as options are important tools for investment and hedging risk.Among a variety of option products, American options have gained popularity in the market because they can be exercised before maturity.This feature makes them the predominant type of options sold in the market.Concurrently, stock loans have become commonplace in the current financial market, and determining how much money can be loaned against the current value of stocks has also emerged as a significant area of research.Given that both of these financial instruments can be exercised (or repaid) at any time, their prices (or loan amounts) can be determined by solving the following linear complementarity problem (LCP): where L is a linear operator containing partial derivatives, U(x, t) represents the value, and U * (x, t) signifies the constraint.
To solve the LCP (1), it is necessary to understand the definition of the operator L, which depends on the assumed underlying asset price behavior.According to the well-known Black-Scholes model [1], the underlying asset price is modeled as follows: where S t is the assets' price, W t represents the standard Brownian motion, ν is the drift rate of the asset's return, and: LU(S, t) := − ∂U(S, t) ∂t with r being the interest rate, σ being the volatility, and S being the asset's price.However, this model cannot account for many empirical observations of asset prices, such as sudden price movements.Hence, new models have been proposed to handle these issues, including stochastic volatility models [2,3], jump diffusion models [4], self-exciting jump models [5], Hawkes jump diffusion models [6], mixed fractional Brownian model [7], two-factor nonaffine stochastic volatility model [8], and sub-fractional Brownian model [9,10].The model based on the Lévy process, notable for its ability to model price jumps and transform into a fractional diffusion equation [11], is among these.With different density functions, different models have been proposed, such as the KoBoL model [12], the CGMY model [13], and the FMLS model [14].The Lévy process with jumps has also been proposed to more accurately describe the price of underlying assets [15,16].In this paper, this stochastic process is used to model the underlying assets, defined by where x t = ln S t and t denotes the current time, with ν = − 1 2 σ α sec( απ 2 ) being the convexity adjustment.The variable L α,−1 t represents the Lévy-α-stable process with maximum skewness, where the tail index α falls within the range of (1,2).N t denotes a Poisson process characterized by the jump intensity λ ≥ 0. {Y i , i = 1, 2 . ..} is a sequence of independent and identically distributed random variables, and ξ = e u J + σ 2 J 2 − 1, where the parameters u J and σ J represent the expectation and standard deviation of the jumps, respectively.Then, the operator LU(x, t) includes both a fractional derivative and an integral operator (for more details, refer to the following section).By solving the LCP in (1) with this LU(x, t) under specific boundary and initial conditions, we can determine the price of American options or the value of stock loans.
To address the LCP arising from the valuation of American options and stock loans, various numerical algorithms have been developed.These can be categorized into five primary strategies.The first focuses on identifying the optimal execution boundary to resolve the challenges, as demonstrated by the approach presented in [17].The second strategy transforms the problem into a linear framework through either semi-implicit or implicitexplicit schemes, such as the L-stable method [18], the linearly implicit predictor-corrector scheme [19], and the IMEX BDF method [20].The third approach converts the problem into a nonlinear equation, utilizing iterative methods for solution acquisition, including the fixed-point method [21], the preconditioned penalty method [22], and the modulus-based matrix splitting iteration method [23].The fourth approach involves devising iterative methods for direct LCP resolution, exemplified by the projected SOR method [24] and the projected algebraic multigrid method [25].The last category of algorithms consists of the recently popular deep learning algorithms, such as the neural network method [26] and the physics-informed neural network method [27].For the LCP discussed in this paper, a Laplace transform method is introduced for rapid resolution by avoiding time marching [15].Moreover, a fast solution strategy that combines Newton's method with the preconditioned conjugate gradient normal residual method is proposed in [16], benefiting from fast Fourier transformation acceleration and achieving an operational cost of O(N log N) per inner iteration.This paper introduces an alternative approach, translating the LCP into a Hamilton-Jacobi-Bellman (HJB) equation, and devising a fast algorithm for its resolution, complete with theoretical assurances.
This study aims to devise a fast algorithm for solving the HJB equation, which includes a fractional derivative and an integral operator derived from the LCP based on the Lévyα-stable process with jumps model.To tackle this equation, a nonlinear finite difference method is developed, accompanied by stability and convergence analyses.Subsequently, the policy iteration method and the Krylov subspace method are applied to solve the resultant finite difference scheme.A banded preconditioner is also proposed to enhance the convergence speed of the internal iterative method, supported by theoretical analysis.These steps enable the development of a preconditioned policy-Krylov subspace method for solving the fractional partial integro-differential HJB equations in finance, ensuring an efficient and theoretically sound solution.
The rest of this article is organized as follows: The description of the equation considered in this paper is illustrated in Section 2. A nonlinear finite difference scheme is proposed in Section 3 to discretize the HJB equation with the unconditional stability guarantee and first-order convergence rate.In Section 4, a fast algorithm framework incorporating the policy iteration method and the Krylov subspace method is introduced.In Section 5, a banded preconditioner is proposed with theoretical analysis about the condition number of the preconditioned matrix.Numerical experiments are given in Section 6, and conclusions are drawn in Section 7.

Description of the Equation
As described in Section 1, the price of the American options and the stock loans can be obtained by solving the LCP in (1).Referring to [28], solving the LCP can be achieved by converting it into an HJB equation.Then the pricing model based on the Lévy-α-stable process with jumps considered in this paper can be solved by solving the following HJB equation, which is defined as follows: where x = ln S, and S is the asset's price.In the above equation, LU(x, t) is defined by where D α − U(x, t) is the left Riemann-Liouville fractional derivative [29], and the definition of the parameters r, ν, λ, ξ has been given in Section 1.The boundary conditions and initial condition are determined by the function Ψ(x, t), which has different meanings in different problems.This paper considers the problems of pricing American options and stock loans; hence, the definition of Ψ(x, t) is given as follows: Ψ(x, t) = max(0, e x − K), for American option pricing, max{0, e x − Ke γt }, for stock loan pricing, for x ∈ (−∞, ∞) and t ∈ [0, T].For the American option pricing, K is the strike price.For stock loan pricing, K is the principal value, and γ is the interest rate of the loan [30].
The main difficulty in solving the HJB Equation ( 2), apart from it being a nonlinear equation, also lies in how to handle the Riemann-Liouville fractional derivative and the integral operator.First is the Riemann-Liouville fractional derivative, which is defined as follows: Because it is a global operator, using the shifted Grünwald approximation [29] will result in a dense coefficient matrix after discretization, thus making the solution challenging.
On the other hand, the integro-differential part +∞ −∞ U(x + y, t) f (y)dy also poses computational difficulties.If the trapezoidal method is used for discretization [15], it will also result in a dense coefficient matrix.Therefore, in summary, when designing a fast algorithm to solve the HJB equation, it is necessary to address the computational challenges posed by a fractional derivative and an integro-differential operator.In addition, differing from Refs.[15,16], the probability density function f (z) of Y i is described by the following Gaussian distribution formula [31]:

Finite Difference Method with Theoretical Analysis
As mentioned above, since LU(x, t) in the HJB Equation ( 2) involves both a fractional derivative and an integral operator, finding the analytical solution to the HJB equation is challenging.Consequently, numerical methods become the primary approach for solving the aforementioned HJB equation.In this section, a finite difference method is developed to discretize the equation, which includes the shifted Grünwald approximation, an upwind scheme and the trapezoidal formulas.

Finite Difference Method
Let N and M be positive integers.Divide the spatial interval [x L , x R ] and temporal interval [0, T] into N + 1 and M sub-intervals, respectively.That is, With above finite difference mesh, the left fractional derivative D α − U(x, t) can be discretized using the shifted Grünwald approximation [29], expressed as for all 1 < α < 2. The sequence {g For the advection term, the upwind scheme [32] is applied, as given by In addition, the integral part in (3) can be approximated by the trapezoidal formulas.To simplify the notation, we denote Ω( û, σ, x) as where dz denotes the cumulative distribution function with expectation û and standard deviation σ, evaluated at the values in x.
To simplify the analysis, we assume that the truncation domain is sufficiently large; i.e., the left and right boundary conditions for the American call option are 0 and e x − K, respectively.As mentioned before, the value outside the domain satisfies that U(z, t) = Ψ(z, t) for z ∈ (−∞, x L ) ∪ (x R , ∞).Similar to the approach in [33], with the notation given in (7), when x = x i , we have and define ρ j as follows, Additionally, η i varies based on the financial product in question.For the American option pricing, η i is defined as follows: Similarly, by assuming that the truncation domain for the stock loan is sufficiently large, then the left and right boundary conditions are set to 0 and e x − Ke γt , respectively.Thus, the operator +∞ −∞ U(x i + z, t m ) f (z)dz, used in the stock loan pricing, can be discretized according to (8), and η i is given by 6) and ( 8), the HJB Equation ( 2) can be discretized as where and

Matrix Form
To simplify the notation of the matrix form, we use T n to denote an n-by-n Toeplitz matrix, which is defined by Then, the matrix form of the numerical scheme (11) can be written as where and I N is an N by N identity matrix; other matrices in ( 14) can be expressed by n−1 , . . ., g (α) 2 ; g where To simplify the notation, we assume that the truncation domain is large, then the specific form of the entry f m i is given as follows:

Stability and Convergence Analysis
Although the numerical scheme has been developed, theorems about the stability and convergence are still needed to ensure the numerical solution can be obtained correctly.For theoretical analysis, the following proposition is introduced.
Proposition 1 ([34]).The sequence {g (α) k } defined in (5) has the following properties: The above proposition provides the properties of the shifted Grünwald approximation.Additionally, one more lemma is needed to analyze the properties of the integral operator, which is given as follows: Lemma 1.With the density function f (z), the entries of the matrix S satisfy Proof.By the definition of s j , we have On the other hand, s j is equal to the probability of a variable X that follows a normal distribution, in the range of (j − 1)h to (j + 1)h, that is, Pr((j − 1)h < X < (j + 1)h).Therefore, its value is definitely greater than 0 for all j.
In the following analysis, the properties of the M-matrix play an important role.Therefore, the definition of the M-matrix is provided first.A matrix is called an M-matrix if it is a diagonally dominant matrix, with all entries on the main diagonal being positive and all off-diagonal entries being negative [35].To prove the stability and convergence of the numerical scheme, combining Proposition 1 and Lemma 1, the properties of matrix W are analyzed.Lemma 2. For 1 < α < 2, the matrix W = [w i,j ] N i,j=1 is an M-matrix, and satisfies Proof.By the definition of D and G, and with Proposition 1, it is straightforward to understand that the matrices τc 1 h D and τc 2 h α G are the weakly diagonally dominant matrices which have positive diagonals and non-positive off-diagonals.Let us denote a matrix W as W = (1 − τc 3 )I N + τc 4 S.Then, we obtain that Based on the definition of W in (15), and given that τs j < 0 and |τ|c 4 > |τ|c 4 s 0 , W is a diagonally dominant matrix with positive diagonals and negative off-diagonals.Since the sum of matrices W, τc 1 h D and τc 2 h α G constitutes the matrix W, we know that W is an M-matrix and satisfies |w To prove that the numerical scheme is stable, we assume that u m i and u m i are both the numerical solutions of the scheme (13).Denote that T , then we have the following theorem.

Proof.
Assume that the i 0 -th entry of the vector We can then proceed with the proof as follows: since u m i and u m i are both the numerical solutions of the scheme, we are given the following two equations: By assuming the existence of the solution, the proof can be discussed in four cases.Case 1: Assume that in the i 0 -th row of the vector u m and u m , two solutions satisfy Then we have the following equation: With Lemma 2, this yields Case 2: In the i 0 -th row of the vector u m and u m , two solutions satisfy . Then, we have On the other hand, this leads to With Lemma 2, we have Case 3: In the i 0 -th row of the vector u m and u m , two solutions satisfy . Then, we have and Case 4: In the i 0 -th row of the vector u m and u m , two solutions satisfy The following inequality can be obtained: By combing the above four cases, we know that the proposed scheme ( 13) is unconditionally stable, that is, Theorem 2. The proposed nonlinear scheme ( 13) is convergent with first order, meaning that where m = 1, 2, . . ., M, and C is a positive constant.
Proof.Similar to the proof of Theorem 1, by denoting R m i = C(h + |τ|), we know that the solutions u m i and U(x i , t m ) satisfy the following two equations: respectively.Denote that e m i = u m i − U(x i , t m ) and R max = max i=1,...,N;m=1,...,M {|R m i |}; then we can start the analysis.Similarly to the proof of Theorem 1, by assuming that . ., N, the proof can be completed by dividing into four cases.Case 1: Two solutions satisfy . Then we have the following equation: With Lemma 2, we have Case 2: Two solutions satisfy ∑ N j=1 w i 0 ,j u m j − u m−1 i 0 Then, On the other hand, this leads to With Lemma 2, the following inequality can be obtained.
Case 3: Two solutions satisfy ∑ N j=1 w i 0 ,j u m j − u m−1 i 0 Then, Case 4: Two solutions satisfy ∑ N j=1 w i 0 ,j u m j − u m−1 i 0 The following inequality can be obtained: By combing the above four cases, we have With Theorems 1 and 2, the proposed numerical scheme has been shown to be stable and convergent.However, given its nature as a nonlinear finite difference scheme, a numerical algorithm is required to address this problem.Consequently, a preconditioned policy-Krylov subspace method is developed in the subsequent section.

Fast Policy-Krylov Subspace Iterative Method
Inspired by [28,36], this paper employs the policy iteration method to solve the discrete HJB equation.Within each policy iteration, a linear system must be solved.The coefficient matrix W being a Toeplitz matrix allows for the consideration of the Krylov subspace iterative method, accelerated with the fast Fourier transformation (FFT), as the inner iterative method for solving this linear system.

Policy Iteration Method
The component-wise form based on the HJB Equation ( 13) can be written as, where (v) i represents the i-th entry of any given vector v when 1 ≤ i ≤ N. ϕ is a control parameter that can take two values, 0 or 1.It is defined that the sequence {u m+1,k } will converge to the solution u m+1 , where k = 0, 1, 2 . . . .Then, the framework of the policy iteration method [28], starting with the initial numerical value u m+1,0 , is presented in Algorithm 1. Similar to the above, let (V) i represent the i-th row for any given matrix V when 1 ≤ i ≤ N.

Algorithm 1 Policy iteration method
Since u m+1,k is the k-th iteration of the policy iteration method for computing the solution for 1 ≤ i ≤ N, where and After introducing the structure of the policy iteration method, an additional lemma is needed to prove that the exact solution can be obtained in a finite number of steps.

Lemma 3.
The coefficient matrix M m,k is invertible and satisfies Proof.From Lemma 2, we know that the matrix W is an M-matrix The coefficient matrix M m,k is composed by selecting rows from the identity matrix I N and the matrix W. Since both I N and W are M-matrices, M m,k is also an M-matrix and satisfies Referring to Theorem 1 in [37], it can also be proven that ∥(M m,k ) −1 ∥ ∞ ≤ 1.
Based on the properties outlined in Lemma 3 and referring to [28], it can be proven that the policy iteration method can be shown to converge to the solution after a finite number of steps.

Fast Krylov Subspace Method
With the policy iteration method, the discrete HJB equation can be solved.However, a linear system exists within each policy iteration that needs to be solved.With Algorithm 1, the coefficient matrix can be expressed as follows: where Φ m,k represents a diagonal matrix, and its i-th diagonal entry is the control parameter ϕ m,k i .Specifically, it is defined as Φ m,k = diag(ϕ m,k 1 , ϕ m,k 2 , . . ., ϕ m,k N ).Similar to the case in [38], the fast Krylov subspace method can be used to solve the linear system, given the Toeplitz structure of the matrix W. Hence, in this paper, the fast policy-Krylov subspace iterative method is utilized to solve the discrete HJB equation.It is noteworthy that the operational cost of each iteration of the Krylov subspace method is only O(N log N), where N is the matrix size.

Preconditioning Technique
Although the preconditioning method proposed in [38] can be utilized to enhance the efficiency of the Krylov subspace method discussed in this paper, our ongoing pursuit of faster and more stable algorithms remains paramount.Consequently, we reference the preconditioning approach in Ref. [39], designing a new banded preconditioner in this section.

Banded Preconditioner
Define the matrix P m,k l as a banded matrix that approximates the coefficient matrix M m,k with 2l − 1 bands.In this paper, we consider the situation where l ≥ 2. The composition of the matrix P m,k l is structured as: 20) where 0 , 0, . . ., 0), S T = T N (0, 0, . . ., 0, s 1−l , s 2−l , . . ., , s −1 ; s 0 ; s 1 , . . ., s l−2 , s l−1 , 0, . . ., 0), After proposing the band preconditioner, it becomes necessary to further analyze the properties, such as the invertibility of the preconditioning matrix and the condition number of the preconditioned matrix.To simplify the subsequent analysis, we drop the indices m and k; that is, we use P −1 l M to represent (P m,k l ) −1 M m,k .

Properties of the Preconditioned Matrix
To ensure the proposed preconditioner is feasible, a theorem is required to prove that P l is invertible.The theorem is stated as follows: Theorem 3. The proposed preconditioner P l is invertible and satisfies Proof.Based on the structure of the matrix P l as described in (20), and a proof similar to that in Lemma 3, it can be easily demonstrated that P l is an M-matrix and satisfies With the above theorem, we can ascertain that the preconditioner P l is feasible.Since the preconditioner P l is a banded matrix, the computational cost of the LU factorization of P l is O(l 2 N), and the cost for the matrix-vector product P −1 l v is O(lN) for any given vector v.With the proposed preconditioning technique, the operational cost of the Krylov subspace iterative method remains O(N log N) per iteration when l ≪ N.However, the computational cost of the preconditioned Krylov subspace method still depends on the number of iterations.In order to prove that the preconditioner proposed in this paper can improve the convergence rate of the Krylov subspace method, the following analysis will examine the condition number of the preconditioned matrix, defined as follows: If the coefficient matrix is ill-conditioned, that is, if the condition number is large, the fast Krylov subspace method converges to the solution slowly.Hence, in the following part, we will prove that the condition number κ ∞ (P −1 l M) has an upper bound.Then, a relatively fast convergence rate of the preconditioned Krylov subspace method can be expected.The following lemmas are provided to support further analysis before analyzing the condition number: Lemma 4 ( [40]).For α ∈ (1, 2), we have with C being a positive constant.
Lemma 5.If u J = 0 and σ J = , where C ′ is a positive constant.In the more general case, when σ J is not specified, one obtains: Proof.When u J = 0, σ J = √ C ′ h, with Chebyshev's inequality, we have Next, we proceed to prove the scenario where σ J is unspecified.Based on the properties of the Q-function, we derive the following equation.
To analyze the condition number of the preconditioned matrix, an additional assumption is needed: there exists a positive constant C such that With the above assumption, it will be shown that there exists an upper bound on the condition number of the preconditioned matrix P −1 l M. Thus, the relatively fast convergence rate of the preconditioned Krylov subspace method can be expected.

Numerical Experiments
In this section, two examples are given to demonstrate the accuracy of the considered model and the effectiveness of the proposed preconditioned iterative method.All numerical experiments are carried out using Matlab 2023a with the following configuration: an Intel(R) Core(TM) i9-10900K CPU at 3.70 GHz and 64 GB RAM.The stopping criterion for the policy iteration method is set to 10 −9 , while the stopping criterion for the Krylov subspace method is 10 −12 .The generalized minimal residual (GMRES) method [41] is chosen to represent the Krylov subspace method.The initial guess of the GMRES method is the zero vector, whereas the initial guess of the policy iteration method is the numerical solution from the previous time step.The GMRES method is restarted every 20 iterations.To simplify the notation in the following section, let N = N + 1.

American Call Option
In the first example, we consider an American option pricing problem.To test the accuracy of the considered model for the option pricing problem in the real financial market, the considered model, i.e., the FMLSJ model, is used to price the sugar call option contract SR403 on 23 November 2023, in the Zhengzhou Commodity Exchange.Additionally, the Black-Scholes (BS) model and the finite moment log stable (FMLS) model are used as comparative models for pricing this option.In this test, r = 2.561%, derived from the annual Shibor rate on the trading day, and T = 56/179, x l = ln(1), x r = ln(20000), the strike prices K = [6500, 6600, 6700, 6800, 6900, 7000, 7100], and the option prices in the market are C = [300.5, 219.5, 167, 117, 76.5, 49, 33.5].Note that the right preconditioner is used in these numerical experiments.
Denote that Θ represents the set of the parameters of the models; V model,i and V market,i are the i-th values of the option from the model and the market, respectively.The particle swarm optimization algorithm is then used to obtain the model's parameters, with the objective function being The specific values of these parameters and the mean square error (MSE) from different model pricing are shown in Table 1.From Table 1, it is evident that the MSE of the FMLSJ model is the smallest, being only about one-third of the error compared to the other two models.To further distinguish the performance of different models in terms of pricing errors, Figure 1 is provided.Figure 1a offers a comparison of the three models against the actual option prices.From the figure, it can be observed that the prices calculated via the FMLS model and the BS model are relatively similar, which explains why their MSEs are close to each other.In contrast, the FMLSJ model's prices are closer to the actual prices.Figure 1b presents the relative errors under different option prices.This figure also demonstrates that the errors of the FMLS model and the BS model are very similar.Meanwhile, the relative errors of the FMLSJ model are all below 7%, indicating superior performance compared to the first two models.
In the second numerical test of this example, the parameters listed in Table 1, obtained from the option contract SR403, are used to assess the performance of the proposed fast solution strategy.The other parameters are set as follows: r = 2.561%, T = 1, x max = ln(100), x min = ln(0.1),K = 50.Additionally, to evaluate the effectiveness of the banded preconditioner used in this paper, the preconditioner based on the circulant approximation from [38] is employed for comparison.Specifically, this preconditioner is defined as: where s(.) denotes Strang's approximation.For simplicity, the following text will use the notation GMRES to refer to the GMRES algorithm without preconditioning, SGMRES to denote the GMRES algorithm with the circulant preconditioner (21), and BGMRES(l) to indicate the GMRES method preconditioned by the banded preconditioner with 2l − 1 bands.Additionally, in the following tables, "Error" represents the infinite norm of the error between the numerical solution and the exact solution, "Rate" represents the convergence rate, "Iter-Out" indicates the average number of iterations of the policy iteration method on each time level, and "Iter-In" refers to the average number of iterations of the Krylov subspace method on each policy iteration.Note that the numerical solution on a fine grid with N = 2 16 and M = 2 12 is regarded as the exact solution since its analytical solution is unknown.
From Table 2, it is evident that since the outer iterations all employ the policy iteration method, the number of outer iterations remains the same, as do the errors and convergence order.Meanwhile, Table 3 demonstrates that compared to circulant preconditioning, the band preconditioning method performs better.Additionally, when comparing different values of l, BGMRES(4) achieves good results with the lowest CPU time, whereas l = 7 and l = 13 do not improve the number of iterations in this example and lead to higher CPU times due to the increased computational effort required for obtaining its inversion.Therefore, l = 4 is found to be sufficient to achieve good results, both in terms of the number of iterations and CPU time.Additionally, to verify the effectiveness of Theorem 4, the condition numbers of the preconditioning matrices for different l will be presented.However, it should be noted that the coefficient matrix in the policy iteration method changes with each time step and iteration.To facilitate verification, the parameters of the model from Table 2 are utilized, with N = 2 14 , M = 2 7 , and Φ = diag(ϕ 1 , ϕ 2 , . . ., ϕ N ), where Then, the condition number of the coefficient matrix M is 365.1499.The condition numbers of the preconditioned coefficient matrix for various l values are depicted in Figure 2. As stated in Theorem 4, it is observed from Figure 2 that as l increases, the condition number gradually decreases and more closely approaches 1.Given that the original condition number reaches 365.1499, it can be seen that the banded preconditioner effectively reduces the condition number of the coefficient matrix, and when l = 2, the condition number is already very close to 1.

Stock Loan
In this part, we consider a stock loan problem to test the performance of the proposed preconditioned policy-Krylov subspace method.The parameters are set as follows: σ = 0.2, λ = 0.01, u J = 0.8, σ J = 0.01, γ = 0.1, r = 0.05, T = 1.Experience from Example 1 indicates that a smaller value of l can achieve satisfactory results, so in this example, the value of l is uniformly set to 4. Additionally, in this example, the value of α will be set to 1.2, 1.5, and 1.8, respectively, to assess the algorithm's performance with different fractional orders.Similar to Example 1, the numerical solution on a fine grid with N = 2 15 and M = 2 11 is regarded as the exact solution.
For α = 1.2, 1.5, and 1.8, the comparison results of different iterative methods are presented in Tables 4, 5 and 6, respectively.The errors in these tables are derived from the results of the BGMRES(4) method.It is observed that the convergence rate is first order for all results.Compared with the results of the iterative method without preconditioning, it is evident that as α increases, more iterations are required.However, both the circulant preconditioner and the proposed banded preconditioner can reduce the number of iterations, with the band preconditioning technique demonstrating significantly better performance in terms of both the number of iterations and CPU time.To further illustrate the effectiveness of the model and the impact of its parameters, Figure 3 is presented.The same parameters from this subsection are used, specified as follows: u J = 0.8, σ J = 0.01, γ = 0.1, r = 0.05, T = 1.Any missing parameters are provided in the figure captions.Additionally, all models are solved using the BGMRES(4) method with settings of N = 2 12 and M = 2 10 .The figures demonstrate that the stock loan prices, similar to American call options, are positioned above the pay-off function.Figure 3a shows that as α increases, the curve approaches the pay-off function, aligning with results from Ref. [42]. Figure 3b illustrates that as volatility σ increases, the price of the stock loan becomes higher.Figure 3c indicates that as λ, the jump frequency, increases, price fluctuations become more frequent, thereby enhancing the value of the stock loan.These figures confirm that the solutions obtained from the considered model meet our expectations for the curve of stock loan valuation.

Conclusions
In this paper, we consider a fractional partial integro-differential HJB equation resulting from the American option pricing and stock loan pricing based on the Lévy-α-stable process with jumps model, and propose a preconditioned policy-Krylov subspace method to solve it.A finite difference scheme is developed to discretize the HJB equation, with stability and first-order convergence analysis.The coefficient matrix generated from the fractional derivative and the integro-differential operator is an M-matrix and possesses the Toeplitz structure, leading us to use the policy iteration method and the fast Krylov subspace method as the outer and inner iterative methods, respectively.To accelerate the convergence rate of the Krylov subspace method, a banded preconditioner is proposed.It is proven that the condition number of the preconditioned matrix has an upper boundary affected by the bandwidth.Finally, numerical examples of an American options pricing problem and a stock loan valuation problem are provided to demonstrate the effectiveness of the proposed fast solution strategy.

Figure 1 .
Figure 1.Comparisons between three models.(a) Option price comparison.(b) Pricing errors of three models.

Figure 2 .
Figure 2. Condition number of the preconditioned matrix.

Table 1 .
Parameters and MSE from different models.

Table 2 .
Comparisons of errors and outer iteration numbers.

Table 3 .
Comparisons between different linear solvers.