
In this paper, the pricing of stock loans when the underlying follows a Lˊevy-α-stable process with jumps is considered. Under this complicated model, the stock loan value satisfies a fractional-partial-integro-differential equation (FPIDE) with a free boundary. The difficulties in solving the resulting FPIDE system are caused by the non-localness of the fractional-integro differential operator, together with the nonlinearity resulting from the early exercise opportunity of stock loans. Despite so many difficulties, we have managed to propose a preconditioned conjugate gradient normal residual (PCGNR) method to price efficiently the stock loan under such a complicated model. In the proposed approach, the moving pricing domain is successfully dealt with by introducing a penalty term, however, the semi-globalness of the fractional-integro operator is elegantly handled by the PCGNR method together with the fast Fourier transform (FFT) technique. Remarkably, we show both theoretically and numerically that the solution determined from the fixed domain problem by the current method is always above the intrinsic value of the corresponding option. Numerical experiments suggest the accuracy and advantage of the current approach over other methods that can be compared. Based on the numerical results, a quantitative discussion on the impacts of key parameters is also provided.
Citation: Congyin Fan, Wenting Chen, Bing Feng. Pricing stock loans under the Lˊevy-α-stable process with jumps[J]. Networks and Heterogeneous Media, 2023, 18(1): 191-211. doi: 10.3934/nhm.2023007
[1] | Bing Feng, Congyin Fan . American call option pricing under the KoBoL model with Poisson jumps. Networks and Heterogeneous Media, 2025, 20(1): 143-164. doi: 10.3934/nhm.2025009 |
[2] | Guillermo Reyes, Juan-Luis Vázquez . The Cauchy problem for the inhomogeneous porous medium equation. Networks and Heterogeneous Media, 2006, 1(2): 337-351. doi: 10.3934/nhm.2006.1.337 |
[3] | Yongqiang Zhao, Yanbin Tang . Approximation of solutions to integro-differential time fractional wave equations in $ L^{p}- $space. Networks and Heterogeneous Media, 2023, 18(3): 1024-1058. doi: 10.3934/nhm.2023045 |
[4] | Diandian Huang, Xin Huang, Tingting Qin, Yongtao Zhou . A transformed $ L1 $ Legendre-Galerkin spectral method for time fractional Fokker-Planck equations. Networks and Heterogeneous Media, 2023, 18(2): 799-812. doi: 10.3934/nhm.2023034 |
[5] | Kexin Li, Hu Chen, Shusen Xie . Error estimate of L1-ADI scheme for two-dimensional multi-term time fractional diffusion equation. Networks and Heterogeneous Media, 2023, 18(4): 1454-1470. doi: 10.3934/nhm.2023064 |
[6] | Hirofumi Notsu, Masato Kimura . Symmetry and positive definiteness of the tensor-valued spring constant derived from P1-FEM for the equations of linear elasticity. Networks and Heterogeneous Media, 2014, 9(4): 617-634. doi: 10.3934/nhm.2014.9.617 |
[7] | Nils Svanstedt . Multiscale stochastic homogenization of monotone operators. Networks and Heterogeneous Media, 2007, 2(1): 181-192. doi: 10.3934/nhm.2007.2.181 |
[8] | Henri Berestycki, Guillemette Chapuisat . Traveling fronts guided by the environment for reaction-diffusion equations. Networks and Heterogeneous Media, 2013, 8(1): 79-114. doi: 10.3934/nhm.2013.8.79 |
[9] | Ciro D’Apice, Umberto De Maio, T. A. Mel'nyk . Asymptotic analysis of a perturbed parabolic problem in a thick junction of type 3:2:2. Networks and Heterogeneous Media, 2007, 2(2): 255-277. doi: 10.3934/nhm.2007.2.255 |
[10] |
Linglong Du .
Long time behavior for the visco-elastic damped wave equation in |
In this paper, the pricing of stock loans when the underlying follows a Lˊevy-α-stable process with jumps is considered. Under this complicated model, the stock loan value satisfies a fractional-partial-integro-differential equation (FPIDE) with a free boundary. The difficulties in solving the resulting FPIDE system are caused by the non-localness of the fractional-integro differential operator, together with the nonlinearity resulting from the early exercise opportunity of stock loans. Despite so many difficulties, we have managed to propose a preconditioned conjugate gradient normal residual (PCGNR) method to price efficiently the stock loan under such a complicated model. In the proposed approach, the moving pricing domain is successfully dealt with by introducing a penalty term, however, the semi-globalness of the fractional-integro operator is elegantly handled by the PCGNR method together with the fast Fourier transform (FFT) technique. Remarkably, we show both theoretically and numerically that the solution determined from the fixed domain problem by the current method is always above the intrinsic value of the corresponding option. Numerical experiments suggest the accuracy and advantage of the current approach over other methods that can be compared. Based on the numerical results, a quantitative discussion on the impacts of key parameters is also provided.
A stock loan is a contract that uses shares of preferred stocks as collateral to secure a loan from another party. There are two positions in a stock loan contract. One is a stock-owing investor who delivers his stocks to a financial institution providing stock loan service. The other is the service provider who gives the investor the right allowing him/her to redeem the stocks at any valid time by repaying the principal and loan interest. As mentioned in [1], stock loans can not only transfer the risk of holding stocks to financial institutions but also overcome the selling restrictions and establish market liquidity. Mathematically, the pricing of stock loans is equivalent to the valuation of American calls with strike prices dependent on the time. In recent years, the trading volumes of stock loans are increasing dramatically, and the pricing of this kind of derivatives has thus received great attention in both industrial and academic areas.
In the literature, the pricing of stock loans was first attempted with the assumption that their maturities are infinite. For example, Xia and Zhou (2007) derived a closed-form analytical solution for the price of stock loans under the Black-Scholes (BS) model [2]. Their work was then explored by Liang et al. (2010) by adding additional features into the formulation, such as the automatic termination clause, cap and margin [3]. Cai and Zhang (2020) proposed a novel and unified framework for the valuation of stock loans with infinite maturity under general regime-switching exponential lˊevy models [4]. However, some researchers pointed out that though the maturities of stock loans are usually quite long, they can never be overlooked by assuming that they are infinite. As a result, they concentrated on pricing stock loans with finite maturities. In this case, it is almost impossible to derive closed-form analytical solutions, and most of the work is done by an approximation method. Examples in this category include using the binomial tree method [5] and the Laplace transform method [6] to price stock loans with different dividend distributions, the Lagrange finite element method to solve the price of stock loans with accumulative dividends [7], the asymptotic expansion method to price stock loans under a fast mean-reverting stochastic volatility model [8], and the projected successive-over-relaxation (PSOR) method to price stock loans in an incomplete market.
However, in most of the work mentioned above, the authors underestimate the probability of large underlying price changes over small time steps, a phenomenon often observed in financial markets. To incorporate the features of both Geometric Brownian motion and jumps, the Lˊevy-α-stable process with jumps is proposed in [9], which could degenerate to the standard Brownian motion, Poisson process and compound Poisson processes by specific parameter settings. Financially, this model takes into consideration the "asymmetric distribution" and "jump" of asset prices, and is able to account for the overall movements of stock prices and significant price changes due to market turbulence. Mathematically, the non-localness resulting from the abnormal diffusion with jumps leads to a FPIDE. In comparison with the BS system with jumps, the current pricing system replaces the second-order spatial derivative in the original system by an α-order spatial derivative, where α∈(1,2]. Financially, the use of the fractional-integro operator implies that the price of the derivative depends on the information of the portfolio over a range of underlying values rather than some localized information.
In the quantitative finance area, the application of fractional calculus attracts interest from a number of researchers. For example, under a modified BS equation with a time-fractional derivative, Chen et al. derived closed-form solutions for double barrier options [10]. Regarding the BS equation with spatial derivatives, Carr and Wu (2002) established the finite moment log-stable (FMLS) model based on market observations, and further investigated its performance empirically [11]. Cartea and del-Castillo-Negrete (2005) solves numerically the price of barrier options under the FMLS model using a finite difference method [12]. Chen et al. considered analytically the pricing of European-style options under different spatial fractional diffusion models, such as the FMLS model[13], the CGMY model[14] and the KoBol model[15]. They also solved the pricing of American options under the FMLS model-based on predictor-corrector framework [16]. Zhou et al. (2018) used the Laplace transform method to solve a series of FPDES and FPIDEs arising in the option pricing field [9].
In this paper, we consider the pricing of stock loans with finite maturity under the Lˊevy-α-stable process with jumps. The moving pricing domain caused by the early exercise nature of stock loans is firstly fixed by adding a small and continuous penalty term. An implicit finite difference scheme is then applied to discretize the resulting FPIDE system. However, due to the non-localness of the fractional-integro differential operator, a full coefficient matrix with Toeplitz structure is obtained, which requires the computational cost in the order of O(M3) and storage space of O(M2), where M is the number of nodal points in the spatial direction [17,18]. Through the application of the preconditioned conjugate gradient normal residual (PCGNR) method with a Strang's circulant pre-conditioner [19] and the fast Fourier transform (FFT) technique, the computational cost reduces significantly from O(M3) to O(MlogM), and the storage space from O(M2) to O(M). Hence the circulant preconditioning technique is used to investigate the option pricing under the framework of fractional diffusion models [20,21,22].
The rest of the paper is organized as follows. In Section 2, the mathematical formulation for the stock loan is provided. The numerical method and corresponding theoretical analysis are presented in Section 3. In Section 4, the PCGNR method is applied. Numerical results and discussions are presented in Section 5. Concluding remarks are given in the last section.
In this section, the Lˊevy-α-stable process with jumps will be introduced in detail, and the pricing of stock loans underneath will be formulated.
The Lˊevy-α-stable process with jumps is defined on a complete filtered probability space (Ω,Ft,P) with the current time t∈[0,T]. Following the assumptions in [9], under this model, the logarithmic price of the underlying, i.e., xt=log(St), satisfies the following stochastic differential equation (SDE)
dxt=(r−D−ν−ξς)dt+σdLα,−1t+d(Nt∑i=1Yi), | (2.1) |
where r, D and t are the risk free interest rate, the dividend and the current time respectively, and ν=−σαsecαπ2 is a convexity adjustment. Lα,−1t is the maximally skewed Lˊevy-α-stable process with tail index α∈(1,2]. Nt is a Poisson process characterized by the jump intensity ξ≥0. {Yi,i=1,2,…} is a sequence of independent and identically distributed hyper-exponential random variables with probability density function
fY(y)=m1∑i=1ˆpiˆθie−ˆθiy1{y≥0}+m2∑j=1˜pj˜θje−˜θjy1{y≤0}. |
Note that ˆpi (i=1,2,...,m1) and ˜pi (i=1,2,...,m2) denote the probabilities of the ith positive and negative jumps, respectively. They satisfies ∑m1i=1^pi+∑m2j=1~pj=1. ˆθi>1 (i=1,...,m1) are the magnitudes of the upward jumps and ˜θj>0 (j=1,...,m2) are that of the downward random jumps. The average jump size is given by
ς=EQ[exp(Y1)−1]=m1∑i=1ˆpiˆθiˆθi−1+m2∑j=1˜pj˜θjˆθj+1−1, | (2.2) |
where EQ is the expectation operator under Q. And according to the results in References[11] and [23], the asset price process is a martingale under Q.
Now, we turn to formulate mathematically the pricing of stock loans under the proposed model. Mathematically, as pointed out in [2], a stock loan can be regarded as an American call option with time-dependent strike price Keγt, where K is the principal and γ (γ≥r) is the continuously compounded loan interest rate. The payoff of the stock loan at maturity is given by
Π(xT,T;α)=max(ex−KeγT,0). |
In addition, similar to American calls, there is a particular value dividing the underlying into two regions at every valid time. One is called the exercise region in which it is optimal for the investor to redeem the stocks. The other is the continuous region in which it is better for the investor to hold the contract. This particular value of the underlying, denoted by Sf, is the optimal redemption price. According to the non-arbitrage principle, the value of the stock loan at t in the continuous region is
V(x,t)=e−r(T−t)EQ[Π(xT,T)|Ft]. |
According to the derivations in [9], V(x,t) satisfies the following FPIDE in the continuous region: for x∈(−∞,xf), t∈[0,T)
∂V(x,t)∂t+(r−D−ν−ξς)∂V(x,t)∂x+ξ∫+∞−∞V(x+y,t)fY(y)dy+ν−∞DαxV(x,t)=(r+ξ)V(x,t), |
where xf=lnSf, 1<α≤2, and the left-sided Riemann-Liouville fractional derivative −∞Dαx is defined as
−∞DαxV(x,t)=1Γ(2−α)∂2∂x2∫x−∞V(z,t)(x−z)α−1dz. | (2.3) |
In addition to the governing Eq (2.3), appropriate boundary conditions are also required in order to actually determine the stock loan price. Similar to [1], we assume the continuity of the stock loan price and its delta across the free boundary to ensure the optimality of holding the stock loan in the continuous region, i.e.,
V(xf,t)=exf−Keγt,∂V(xf,t)∂x=Sf=exf. |
Therefore, by taking into consideration the boundary condition at x=−∞ (S=0) and the terminal condition at T, the price of the stock loan under the Lˊevy-α-stable process with jumps satisfies
{∂V(x,t)∂t+(r−D−ν−ξς)∂V(x,t)∂x+ξ∫+∞−∞V(x+y,t)fY(y)dy+ν−∞DαxV(x,t)=(r+ξ)V(x,t),limx→−∞V(x,t)=0,V(xf,t)=exf−Keγt,∂V(xf,t)∂x=exf,V(x,T)=max(ex−KeγT,0). | (2.4) |
We remark that the above FPIDE system is much more difficult to solve than the corresponding BS case with jumps, with the main difficulty resulting from the non-localness of the fractional-integro differential operator. In the following section, a new numerical scheme is proposed to solve for Eq (2.4) efficiently.
Upon establishing a closed pricing system Eq (2.4) for stock loans when the underlying satisfies the Lˊevy-α-stable process with jumps, a new numerical approach is proposed to solve the nonlinear FPIDE system Eq (2.4) efficiently. This will be illustrated in detail in this section.
To simplify the solution process, we use the following new variables
z=x−γt,U(z,t)=e−γtV(x,t),zf=xf−γt. |
With the transformation details provided in Appendix A, the pricing system Eq (2.4) becomes: for z∈(−∞,zf], t∈[0,T]
{∂U(z,t)∂t+a∂U(z,t)∂z+ξ∫+∞−∞U(z+y,t)fY(y)dy+ν−∞DαzU(z,t)=(r+ξ−γ)U(z,t),U(zf,t)=ezf−K,∂U(zf,t)∂z=ezf,limz→−∞U(z,t)=0,U(z,T)=max(ez−K,0), | (3.1) |
where a=(r−ν−ξς−D−γ)<0.
Clearly, the undetermined function U(z,t) can now be viewed as the price of an American call with strike K and a free boundary ezf. Since the American call option value is always greater than or equal to its payoff value, we directly have U(z,t)≥max(ez−K,0).
In fact, the FPIDE system Eq (3.1) is nonlinear, with the nonlinearity arising from the early exercise opportunity of the stock loan. As described in the Reference [24], if one would like to simply solve a partial differential equation that automatically fulfills the extra requirement. An equation that approximates this property fairly well can be derived by adding a nonlinear penalty term to governing. And our pricing mathematical model can be solved by this method, therefore we define
εHUε(z,t)+ε−q(z), | (3.2) |
where ε is a regularization constant and 0<ε≪1, H is a constant, and q(z) is defined as q(z)=ez−K. With the penalty term added on and the truncation of the domain applied, the FPIDE system becomes: for z∈[zmin,zmax], t∈[0,T]
{∂Uε(z,t)∂t+a∂Uε(z,t)∂z+ξ∫+∞−∞Uε(z+y,t)fY(y)dy+ν−∞DαzUε(z,t)+εHUε+ε−q=(r+ξ−γ)Uε(z,t),Uε(zmin,t)=0,Uε(z,T)=max(ez−K,0),Uε(zmax,t)=ezmax−K, | (3.3) |
where exp(zmax) and exp(zmin) denote the maximum and minimum stock prices, respectively. In particular, according to the reasons provided in [16] and [25], we set zmax=ln(3K) and zmin=ln(0.01). In the following, for simplicity purposes, the subscript ε of Uε(z,t) is omitted unless otherwise stated. We remark that though the above FPIDE system Eq (3.3) is nonlinear, the moving boundary disappears and the solution domain now becomes fixed as [zmin,zmax]. Moreover, the new system Eq (3.3) is closely related to the original system Eq (3.1) in the sense that the solution of the former approaches that of the latter once ε→0. In addition, the solution of Eq (3.3) is also greater than the corresponding intrinsic value in a discrete sense. This issue will be further explored both theoretically and numerically in the following work, and is also one of the innovations of the current paper.
Now, place respectively 2M+3 and N+1 uniform grids in the z and t directions. The grid sizes are defined as Δz=zmaxM and Δt=TN for the spatial and time directions, respectively. Denote the nodal value in the z direction as zj=(j−1)Δz, for j=−(M+1),⋯,−1,0,1,⋯,(M+1), and in the t direction as ti=(i−1)Δt, for i=1,2,⋯,(N+1).
Next, the Shifted Grünwald-Letnikov formula [21] is employed to approximate the left-sided Riemann-Liouville fractional derivative, i.e.,
−∞DαxUij≈1(Δx)α∞∑k=0gkUij−k+1, |
where gk are defined as
gk=(−1)k(αk),where(αk)=α(α−1)⋯(α−k+1)k!. |
In addition, the integral term contained in the governing equation of Eq (3.3) is approximated by the trapezoidal rules [9], i.e.,
∫+∞−∞U(zj+y,ti)fY(y)dy≈M∑ℓ=0ρMℓ−j[Uiℓ+Uiℓ+1]+Rj, |
where
ρMj=12∫(j+1)ΔzjΔzfY(y)dy={12∑m1ℓ=1ˆpℓ(e−ˆθℓjΔz−e−ˆθℓ(j+1)Δz),j≥0,12∑m2ℓ=1˜pℓ(e˜θℓ(j+1)Δz−e˜θℓjΔz),j≤0, | (3.4) |
and
Rj=∫+∞zM+1−zj(ezj+y−Kerti)fY(y)dy,=ezjm1∑ℓ=1ˆpℓˆθℓˆθℓ−1e(1−ˆθℓ)(zmax−zj)−Kertim1∑ℓ=1ˆpℓe−ˆθℓ(zmax−xj),=(ezmax−Kerti)m1∑ℓ=1ˆpℓˆθℓe−ˆθℓ(zmax−xj). | (3.5) |
Therefore, after the full implicit difference scheme is applied, the FPIDE system Eq (3.3) becomes
Ui+1j−UijΔt+aUij−Uij−1Δz+ν(Δz)α∞∑k=0gkUij+1−k+M∑ℓ=0ρMℓ−j[Uiℓ+Uiℓ+1]+Rj+ϵHUij+ϵ−qj+(γ−ξ−r)Uij=0, | (3.6) |
with
Ui−(M+1)=0,UiM+1=exmax−K,UN+1j=max(exj−K,0), |
where Uij denotes the approximation solution of U(x,t) at (xj,ti).
With the discretization in hand, we shall show that the solution of the difference Eq (3.6) is always greater than or equal to the corresponding payoff value, which is an important property inherited from that of the stock loan. Prior to display the proof, two lemmas need to be proved first.
Lemma 3.1. Both the coefficients ρMj in Eq (3.4) and Rj in Eq (3.5) are bounded and satisfy
M∑−(M1+1)ρMj≤12,Rj≤ezmax−K. | (3.7) |
Proof. Since fY(y) is the density of a hyper-exponential random variable Y, we have
0≤fY(y)≤1. |
Therefore, it is clear that
M∑j=−∞ρMj=M∑j=−∞12∫(j+1)ΔzjΔzfY(y)dy≤12∫+∞−∞fY(y)dy=12. |
On the other hand, since ˆpi≥0,ˆθi>1 and exp(zmax)=3K, it is straightforward to obtain
Rj=(ezmax−Kerti)m1∑ℓ=1ˆpℓˆθℓe−ˆθℓ(zmax−xj)≤(ezmax−Kerti)m1∑ℓ=1ˆpℓ≤ezmax−K. |
The proof of this lemma is thus completed.
Lemma 3.2. For α∈(1,2], gk (k=0,1,⋯,+∞) satisfy
g0=1,g1=−α,0⩽⋯≤g3≤g2≤1,∞∑k=0gk=0,m∑k=0gk<0, |
where m≥1.
Proof. The proof of this lemma is the same as the one provided in [25] and [26], and is thus omitted. Readers can refer to these two references for interest.
Theorem 3.1. If Δt satisfies Δt≤1|γ−r−ξ−1|, and the constant H is bounded below as
H≥∣γ−r−ξ∣ezmax+∣a∣ezmax+ν(ezmax−1)α(zmax)α+ˉK, |
where ˉK=2[exp(zmax)−K], the approximated stock loan values {Uij} satisfy
Uij≥max(qj,0),j=−(M+1),⋯,−2,−1,0,⋯,(M+1),i=1,2,⋯,(N+1). |
Proof. We shall first show that Uij≥qj. Define uij=Uij−qj. It is straightforward that uN+1j=UN+1j−qj≥0. Substituting uij into the governing equation contained in Eq (3.3) yields
[1−(γ−r−ξ)Δt−aΔtΔx]uij=ui+1j−aΔtΔzuij+νΔt(Δz)α∞∑k=0gkuij−k+1+εHuij+ε−qj+M∑ℓ=0ρMℓ−j[uiℓ+uiℓ+1]Δt−ΔtEj | (3.8) |
where
Ej=−(γ−r−ξ)ezj−aΔz(ezj−1−ezj)−ν(Δz)α∞∑k=0gk(ej−k+1−K)−M∑ℓ=0ρMℓ−j[ezℓ+ezℓ+1−2K]−Rj. |
Since |eΔz−1Δz|≤ezmax−1zmax≤1 and ∑∞k=0gk=0, we have
∣Ej∣≤|γ−r−ξ|ezmax+|a|ezmax+Rj+|M∑ℓ=0ρMℓ−j[ezℓ+ezℓ+1−2K]|+|ν(Δz)α∞∑k=0gkej−k+1|. |
In addition, according to [27], we have ∞∑k=0gkezj−k+1=ezj+1∞∑k=0gke−kΔz and (1−z)α=∞∑k=0gkz−k for |z|≤1. It is clear that
|ν(Δz)α∞∑k=0gke−kΔt|=|ν(1−e−Δx)α(Δz)α|≤ν(ezmax−1)α(zmax)α. |
Now, let ˉK=2[exp(zmax)−K], and we have
|M∑ℓ=0ρMℓ−j[ezℓ+ezℓ+1−2K]|≤|M∑ℓ=0ρMℓ−jˉK|≤ˉKM∑ℓ=0ρMℓ−j≤12ˉK. |
Therefore, combining with Lemma 3.2, it is straightforward to obtain
|Ej|≤|γ−r−ξ|ezmax+|a|ezmaxezmax−1ezmax+ν(ezmax−1)α(zmax)α+ˉK. |
Define
ui=minjuij, |
and let J be an index such that uiJ=ui. Since a<0, ρMj>0 for all j, and gk>0 for k≠1, one can deduce from Eq (3.6) that
[1−(γ−r−ξ)Δt−aΔtΔz]ui≥ui+1J−aΔtΔzui+νΔt(Δz)α∞∑k=0gkui+εHuj+ε−qj+2M∑ℓ=0ρMℓ−juiΔt−ΔtEj. |
After simple algebraic operations, we obtain
[1−(γ−r−ξ−2M∑ℓ=0ρMℓ−j)Δt]ui−εHui+ε+ΔtE(j)≥ui+1J≥ui+1. |
On the other hand, according to Eq (3.7) and Δt≤1|γ−r−ξ−1|, it is straightforward that
1−(γ−r−ξ−2M∑ℓ=0ρMℓ−j)Δt≥0. |
For the ease of description, let
A=[1−(γ−r−ξ−2M∑ℓ=0ρMℓ−j)Δt]. |
Define a function f(x) as
f(x)=Ax−εΔtCx+ε+ΔtEj. |
According to the definition of f(x), it is clear that f(ui)≥0 if ui+1≥0. Since f(0)=Δt(Ej−H)≤0, f′(x)=A+ϵΔtC(x+ϵ)2>0, and uN+1j≥0, we obtain ui≥0, and consequently uij≥0. Therefore, Uij≥qj is satisfied.
Next, we show that Uij≥0. Similarly, define Ui=minjUij, and let J be an index such that UiJ=Ui. From Eq (3.6), it is clear that
[1−(γ−r−ξ)Δt−aΔtΔx]Ui≥Ui+1J−aΔtΔzUi+νΔt(Δz)α∞∑k=0gkUi+εHUi+ε−qj+2M∑ℓ=0ρMℓ−jUiΔt, |
which, after a careful rearrangement, yields
[1−(γ−r−ξ−2M∑ℓ=0ρMℓ−j)Δt]Ui≥εΔtCUiJ+ε−qJ+ΔtRJ+Ui+1. |
On the other hand, since Rj>0 and ΔtεHUiJ−ε+qJ>0 because of the fact that Uij≥qj for all i,j, we have
[1−(γ−r−ξ−2M∑ℓ=0ρMℓ−j)Δt]Ui≥Ui+1. |
Since 1−(γ−r−ξ−2∑Mℓ=0ρMℓ−j)Δt>0, one can deduce that
[1−(γ−r−ξ−2M∑ℓ=0ρMℓ−j)Δt]N+1−iUi≥UN+1. |
Since UN+1≥0 due to the fact that UN+1j=max[exp(zj)−K,0]≥0 for all j, we finally have Uij≥0. Thus, the proof is thus completed.
We remark that the constraint condition we set on Δt, i.e., Δt≤1|γ−r−ξ−1| is indeed not harsh. This is because the parameters γ and r satisfy γ∈(0,1) and r∈(0,1), respectively, and moveover the value of ξ≤10 is also reasonable according to the discussions provided in [28]. Financially, Theorem 1 ensures that the discrete value Uij cannot fall below the corresponding intrinsic value, a characteristic inherited from the early exercise nature of the stock loan.
Since the discrete system Eq (3.6) is still nonlinear, the Newton iteration method is used. However, the non-localness of the fractional-integro differential operator results in coefficient matrix in a dense form. Further exploration of the structure of the coefficient matrix is thus needed to enhance the computational efficiency while decreasing the storage space[29], which is the main concern of the current subsection.
For illustration purposes, we reset the index j and refine the grid size in the z direction as Δz=(zmax−zmin)M and zj=zmin+(j−1)Δz, for j=1,2,⋯,M+1. Define
η1=aΔtΔz,η2=−Δtν(Δz)α,η3=1−Δt(γ−ξ−r)−η1, | (4.1) |
and
WMl=ρMl+ρMl−1,l=0,±1,±2,...,±(M+1). | (4.2) |
Then the matrix form of Eq (3.6) can now be written as
(η3I(M−1)×(M−1)+η1B+η2A−ΔtW)Ui−F(Ui)=Ui+1−Ei−ΔtRi. | (4.3) |
where Ui=(Ui2,Ui2,...,UiM)τ, F(Ui)=(f(Ui2),f(Ui3),...,f(UiM))τ,
f(Uij)=ΔtϵHUij+ϵ−qj, |
I(M−1)×(M−1) is the identity matrix,
W=[WM0WM1WM2⋯WMM−2WMM−1WM−1WM0WM1⋯WMM−3WMM−2WM−2WM−1WM0⋯WMM−4WMM−3⋮⋱⋱⋱⋮⋮WM2−MWM3−MWM4−M⋯WM0WM1WM1−MWM2−MWM3−M⋯WM−1WM0], |
B=[0⋯001⋯00⋮⋮⋮⋮0⋯000⋯10]M−1×M−1,A=[g1g00⋯0g2g10⋯0g3g2g0⋯0⋮⋱⋱⋱⋮gM−2gM−3⋯g1g0gM−1gM−2⋯g2g1], |
Ri=(Ri2,Ri3,...,RiM)τ, Ei=(0,0,...,ϱ)τM−1 and ϱ=(−η2g0−∑M−1ℓ=0ρℓ)UiM+1.
After a Newton-type iteration approach is applied, Eq (4.3) becomes
[η3I+η1B+η2A−ΔtW−JF(ωl−1)]δωl=Ui+1−Ei−ΔtRi+F(ωl−1)−[η3I+η1B+η2A−ΔtW]ωl−1ωl=ωl−1+κ(ωl−ωl−1), | (4.4) |
where l=1,2,⋯, the initial guess ω0 is set as ω0=Ui+1 for each time level, and δωl=ωl−ωl−1. Moreover, JF is Jacobian matrix with column vector F(ωl), and κ∈(0,1) is the damping parameter. We set Ui=ωl once the stopping criterion ∥ωl−ωl−1∥∞≤tol for some l is reached, where tol is the stopping tolerance of the iterative method. By denoting
M=η3I+η1B+η2A−ΔtW,bl=Ui+1−Ei−ΔtRi+F(ωl−1)−[η3I+η1B+η2A−ΔtW]δωl−1, |
Equation (4.4) can be rewritten as
[M−JF(ωl−1)](δωl)=bl. | (4.5) |
The most challenging part in solving Eq (4.4) is the high computational cost resulting from the fact that both A and W are dense matrices. To overcome this difficulty, we first apply the CGNR method [30], which is in fact to solve [M−JF]τMδωl=[M−JF]τbl instead of Eq (4.5).
However, by noticing that the convergence rate of the CGNR method is still quite low due to the fact that the conditional number of the matrix MTM is large, a pre-conditioner technique is applied to accelerate the convergence rate of the CGNR method. It is straightforward to find that the matrix JF is not Toeplitz matrix and we should approximate this matrix as a0I, where a0 is the average value of main diagonal elements of matrix JF. So we structure a Toeplitz matrix as following
T=M−a0I. |
Next, the Strang's circulant preconditioner s(T)=[sj−k]0≤j,k<M for matrix T is structured as
sj={Tj,0≤j<M/2,0,j=M/2ifMiseven,andj=(M+1)/2if M is odd,Tj−M,M/2<j<M,Tj+M,0<−j<M. |
We use P denotes the Strang's circulant preconditioner s(T)=[sj−k]0≤j,k<M to simplify the expression.
Mathematically, after the PCGNR method with a pre-conditioner P is applied, Eq (4.4) becomes
[(P)−1(M−JF)]T[(P)−1(M−JF)]δωl=[(P)−1(M−JF)]T(P)−1bl, |
The pseudo-code of the PCGNR method is displayed in Algorithm 1. The matrix-vector multiplication only needs O(MlogM) operations via the fast Fourier transform (FFT) method.
Algorithm 1. PCGNR method for solving (M−JF)(δωl)=bl with a pre-conditioner P. |
Given the initial guess x0, and a stopping tolerance tol. |
Compute r0=[P−1(bl−(M−JF))x0], z0=[(P)−1(M−JF)]τr0,p0=z0,mr=||r0||22. |
For i=0,1,..., |
wi=[(P)−1(M−JF)]τpi, |
αi=||zi||22/||wi||22, |
xi+1=xi+αipi, |
ri+1=ri−αiwi, |
zi+1=[(P)−1(M−JF)]Tri+1, |
βi=||zi+1||22/||zi||22, |
pi+1=zi+1+βipi, |
res=||ri+1||22. |
If res/mr<tol, stop; |
otherwise, set δωl=xi+1. |
End for |
This section investigates the performance of the newly proposed numerical scheme and the impact of the introduction of the fractional diffusion with jumps in the optimal redemption price of the stock loan. All numerical computations in this section are carried out by Matlab2020 on a Lenovo S5 machine. In addition, the damping parameter κ and the stopping tolerance tol are set to be 0.2 and 10−6, respectively.
Before showing the performance of the current scheme, we examine whether or not the numerical solution preserves the basic properties of stock loans. This could be viewed as a necessary condition for the reliability of the proposed approach. Theoretically, our numerical solution should satisfy U(x,t)≥max(ez−K,0) discretely, as shown in Theorem 3.1. In Figure 1, the surfaces of the difference between Uij and max(qj,0) with different parameter settings are displayed, which implies that the inequality is indeed preserved. On the other hand, as pointed out in previous sections, U(z,t) contained in Eq (3.1) is equivalent to an American call option with strike price K and optimal exercise price ezf. In Figure 2(a), we display U(z,t) against the stock price ez. From this figure, it is clear that the current numerical solution satisfies well the smooth pasty condition across the free boundary. In addition, in Figure 2(b), the optimal exercise price ezf as a function of the time to expiry is plotted. One can clearly observe from this figure that ezf is increasing with respect to (w.r.t) the time to expiry, which is another important property of American calls. All these suggest that the current numerical solution is indeed reliable, and no additional artificial errors are brought in when implementing the proposed method on a computer.
To further investigate the performance of the current method, we compare the computational efficiency of the Gaussian elimination (GE), the CGNR method, and the PCGNR method, as shown in Table 1. The parameters adopted for computing this table are α=1.52, σ=0.2, r=0.05, D=0.06, K=2, T=0.2, ξ=0.03, ˆp=0.04, γ=0.06, ˆθ=1.2, ˜θ=0.2,m1=1,m2=1,ˆp=0.5,˜p=0.5. Moreover, in this table, Ite-In denotes the average iteration number required in each time step, and the error is defined as Err=∥U−u(z,t)ref∥2, where ∥⋅∥2 is the L2 norm for matrix, and u(z,t)ref is the benchmark solution was determined directly through matrix operation 'A∖b' in Matlab with M=211 and N=500.
GE | CGNR | PCGNR | ||||||
M | Time(s) | Err | Ite−In | Time(s) | Err | Ite−In | Time(s) | Err |
25+1 | 33.6213 | 0.0636 | 49.7512 | 0.8722 | 0.0708 | 4.8703 | 0.7531 | 0.0901 |
26+1 | 137.7118 | 0.0249 | 49.7512 | 0.9842 | 0.0290 | 4.8702 | 0.8753 | 0.0352 |
27+1 | 611.6442 | 0.0108 | 45.8657 | 2.8870 | 0.0120 | 5.5608 | 1.1039 | 0.0101 |
28+1 | 2469.0242 | 0.0044 | 44.4080 | 7.3217 | 0.0058 | 6.2370 | 4.5179 | 0.0046 |
29+1 | 18258.5597 | 0.0019 | 44.3532 | 43.3345 | 0.0028 | 6.2371 | 16.9851 | 0.0020 |
210+1 | ** | ** | 45.8301 | 223.1818 | 0.0013 | 6.8123 | 18.9329 | 0.0016 |
One can clearly observe from Table 1 that for a fixed No. of nodal points, the total CPU times required by the CGNR and PCGNR to produce the same level of error are significantly less than that of the GE. Furthermore, the average inner iteration numbers required by the PCGNR method are the least. These suggest the superiority of the PCGNR method in computational efficiency over the GE and CGNR method.
Besides the computational efficiency, the determination of the order of convergence of a particular numerical method is also equally important. Theoretically, from the adopted discretization scheme, the current method should be first-order and second-order convergent in the time direction and spatial directions, respectively. We now turn to investigate this issue numerically. To check the order of convergence in t direction, we fix the size in the x direction to be fairly small as Δz=zmax−zmin211, and increase the grid number in the t direction from 100 to 400. The results are displayed in Table 2. In this table, the notation 'Order' denotes the convergence order and is defined as
Orderi=ln(Erri)−ln(Erri−1)ln(Li)−ln(Li−1), |
Number of time steps | Error | Order |
100 | 0.2005 | – |
200 | 0.1226 | 0.7098 |
300 | 0.0886 | 0.8015 |
400 | 0.0691 | 0.8653 |
where i denotes the ith of Table 2, and Li (i=2,3,4) is the number of time steps displayed in the ith line of this table. From this table, it is clear that our scheme is first-order convergent in the time direction. Similarly, the convergence order in the spatial direction is also examined, and the results are displayed in Table 3 From this table, one can conclude that the convergence order in the x direction of our method is around 1.5. The slight loss of convergence order in the spatial direction may result from the approximation method adopted for approximating the integral and factional derivative. Note that the parameters adopted to compute Table 2 and Table 3 are the same as those used in Table 1.
Number of spatial steps | Error | Order |
33 | 0.0927 | – |
65 | 0.0381 | 1.2828 |
128 | 0.0126 | 1.5964 |
257 | 0.0048 | 1.3923 |
With the validation of the proposed numerical approach, it suffices for us to examine the impacts of the introduction of the factional diffusion with jumps on the optimal redemption strategy of stock loans. For illustration purposes, we sometimes need to convert U and ezf back to V and Sf via V(x,t)=eγtU(z,t) and Sf(t)=eγt+zf, respectively.
We first investigate how the discrete jumps affect the optimal redemption strategy of the stock loan. Depicted in Figure 3 is the optimal redemption price as a function of the time to expiry with different jump intensity ξ. One can observe clearly from this figure that the optimal redemption price increases w.r.t ξ. Financially, a larger jump intensity means that the stock price would change more often, and the stock loan contract will be more valuable because it contains more risks now. According to the smooth pasty condition across the free boundary, the monotonicity of Sf w.r.t ξ holds automatically. This figure also reveals the fact that Sf is not monotonic w.r.t the time to expiry. This is not surprising, and could also be plausibly explained by the same reasons provided in [1].
In Figure 4, the optimal redemption price is plotted against the time to expiry with different probabilities of positive jumps ˆp. From this figure, it is clear that a larger ˆp results in a lower optimal redemption price. In fact, the logarithmic return of the stock price is decreasing w.r.t ˆp, because the return decreases w.r.t ξ from Eq (2.1) and ξ increases w.r.t ˆp from Eq (2.2). Therefore, an increasing ˆp would lower the stock price, and thus makes the intermediate American call option U(z,t) less valuable. Therefore, the optimal exercise price ezf of the intermediate American call decreases w.r.t ˆp. Since the optimal redemption price of the stock loan Sf is related to zf through Sf(t)=eγt+zf, it is straightforward that a larger ˆp results in a lower optimal redemption price. Similarly, one could explain the monotonicity of the optimal redemption price w.r.t ˆθ and ˜θ. For the length of the paper, we provide those figures in Figure 5 with no detailed explanations.
Finally, we examine how the tail index α influences the optimal redemption strategy. Several sets of optimal redemption prices with different α values are computed and displayed in Figure 6. From the curves in the figure, it is straightforward to find that the optimal redemption price is monotonically decreasing w.r.t the tail index α. Financially, the tail index controls the tail of the distribution of the underlying price, and both tails will be fatter when α becomes smaller, as suggested in [1]. Therefore, as α becomes smaller, the possibility of larger stock prices increases, and so does the optimal redemption price.
In this paper, the pricing of stock loans under the Lˊevy-α-stable process with jumps is investigated. The valuation problem is formulated as solving a nonlinear FPIDE system with a free boundary denoting the optimal redemption price. By adding a penalty term, the solution domain becomes fixed, and the resulting system is then solved by a PCGNR method. Numerical experiments suggest that the current method is indeed reliable and does have advantages in computational speed and storage space. Based on the numerical results, some future research directions are expected. Firstly, it is feasible to apply the current method to solve the pricing of more complicated stock loans under the fractional diffusions with jumps, such as stock loans with automatic termination clauses, caps and margins. Secondly, it is also promising to extend the current method to the pricing of other American-style derivatives under the current model. Lastly, to further enhance the efficiency of the current method, some other linearization techniques could be explored and adopted together with the PCGNR method.
This work is supported by the Ministry of Education of Humanities and Social Science Project of China No. 21YJAZH005, Research Base of Humanities and Social Sciences outside Jiangsu Universities "Research Center of Southern Jiangsu Capital Market" (2017ZSJD020), Natural science research project of Guizhou Education Department with grant KY[2021]276 and Natural science research project of Guizhou University of commerce (American Contingent Claims valuation with regime-switching under the Framework of Machine Learning).
The authors declare there is no conflict of interest.
Since
z=x−γt,U(z,t)=e−γtV(x,t),zf,k=xf,k−γt, |
we have
∂V(x,t)∂t=γeγtU(z,t)+eγt∂U(z,t)∂t−γeγt∂U(k;z,t)∂z, | (A.1) |
∂V(x,t)∂x=eγt∂U(z,t)∂z,−∞DαxV(x,t)=1Γ(2−α)d2dx2∫x−∞V(y,t)(x−y)α−1dy=1Γ(2−α)d2dz2∫z+γt−∞V(y,t)(z+γt−y)α−1dy. | (A.2) |
Let s=y−γt, and we have
−∞DαxV(x,t)=eγtΓ(2−α)d2dz2∫z−∞e−γtV(s+γt,t)(z−s)α−1ds=eγtΓ(2−α)d2dz2∫z−∞U(s,t)(z−s)α−1ds,=eγt−∞DαzU(z,t). | (A.3) |
Now, substituting Eqs (A.1)–(A.3) into the pricing system Eq (2.4), the FPIDE system Eq (3.1) can finally be obtained.
[1] |
W. Chen, L. Xu, S. P. Zhu, Stock loan valuation under a stochstic interest rate model, Comput. Math. Appl., 70 (2015), 1757–1771. https://doi.org/10.1016/j.camwa.2015.07.019 doi: 10.1016/j.camwa.2015.07.019
![]() |
[2] |
J. Xia, X. Y. Zhou, Stock loans, Math. Financ., 17 (2007), 307–317. https://doi.org/10.1111/j.1467-9965.2006.00305.x doi: 10.1111/j.1467-9965.2006.00305.x
![]() |
[3] |
Z. Liang, W. Wu, S. Jiang, Stock loan with automatic termination clause, cap and margin, Comput. Math. Appl., 60 (2010), 3160–3176. https://doi.org/10.1016/j.camwa.2010.10.021 doi: 10.1016/j.camwa.2010.10.021
![]() |
[4] |
N. Cai, W. Zhang, Regime classification and stock loan valuation, Oper. Res., 68 (2020), 965–983. https://doi.org/10.1287/opre.2019.1934 doi: 10.1287/opre.2019.1934
![]() |
[5] |
M. Dai, Z. Q. Xu, Optimal redeeming strategy of stock loans with finite maturity, Math. Financ., 21 (2011), 775–793. https://doi.org/10.1111/j.1467-9965.2010.00449.x doi: 10.1111/j.1467-9965.2010.00449.x
![]() |
[6] |
X. Lu, E. R. M. Putri, Semi-analytic valuation of stock loans with finite maturity, Commun. Nonlinear. Sci., 27 (2015), 206–215. https://doi.org/10.1016/j.cnsns.2015.03.007 doi: 10.1016/j.cnsns.2015.03.007
![]() |
[7] |
D. Prager, Q. Zhang, Stock loan valuation under a regime-switching model with mean-reverting and finite maturity, J. Syst. Sci. Complex., 23 (2010), 572–583. https://doi.org/10.1007/s11424-010-0146-7 doi: 10.1007/s11424-010-0146-7
![]() |
[8] |
T. W. Wong, H. Y. Wong, Stochastic volatility asymptotics of stock loans: valuation and optimal stopping, J. Math. Anal. Appl., 394 (2012), 337–346. https://doi.org/10.1016/j.jmaa.2012.04.067 doi: 10.1016/j.jmaa.2012.04.067
![]() |
[9] |
Z. Zhou, J. Ma, X. Gao, Convergence of iterative Laplace transform methods for a system of fractional PDEs and PIDEs arising in option pricing, E. Asian. J. Appl. Math., 8 (2018), 782–808. https://doi.org/10.4208/eajam.130218.290618 doi: 10.4208/eajam.130218.290618
![]() |
[10] |
W. Chen, X. Xu, S. P. Zhu, Analytically pricing double barrier options based on a time-fractional Black-Scholes equation, Comput. Math. Appl., 69 (2015), 1407–1419. https://doi.org/10.1016/j.camwa.2015.03.025 doi: 10.1016/j.camwa.2015.03.025
![]() |
[11] |
P. Carr, H. Geman, D. B. Madan, Y. Marc, The fine structure of asset returns: an empirical investigation, J. Bus., 75 (2002), 305–332. https://doi.org/10.1086/338705 doi: 10.1086/338705
![]() |
[12] |
A. Cartea, D. del-Castillo-Negrete, Fractional diffusion models of option prices in markets with jumps, Physica A, 374 (2007), 749–763. https://doi.org/10.1016/j.physa.2006.08.071 doi: 10.1016/j.physa.2006.08.071
![]() |
[13] |
W. Chen, X. Xu, S. P. Zhu, Analytically pricing European-style options under the modified Black-Scholes equation with a spatial-fractional derivative, Q. Appl. Math., 72 (2014), 597–611. https://doi.org/10.1090/S0033-569X-2014-01373-2 doi: 10.1090/S0033-569X-2014-01373-2
![]() |
[14] |
W. Chen, M. Du, X. Xu, An explicit closed-form analytical solution for European options under the CGMY model, Commun. Nonlinear. Sci., 42 (2017), 285–297. https://doi.org/10.1016/j.cnsns.2016.05.026 doi: 10.1016/j.cnsns.2016.05.026
![]() |
[15] |
W. Chen, S. Lin, Option pricing under the KoBol model, Anziam. J., 60 (2018), 175–190. https://doi.org/10.1017/S1446181118000196 doi: 10.1017/S1446181118000196
![]() |
[16] |
W. Chen, X. Xu, S. P. Zhu, A predictor-corrector approach for pricing American options under the finite moment log-stable model, Appl. Numer. Math., 97 (2015), 15–29. https://doi.org/10.1016/j.apnum.2015.06.004 doi: 10.1016/j.apnum.2015.06.004
![]() |
[17] |
X. M. Gu, T. Z. Huang, C. C. Ji, B. Carpentieri, A. A. Alikhanov, Fast iterative method with a second-order implicit difference scheme for time-space fractional convection-diffusion equation. J. Sci. Comput., 72 (2017): 957–985. https://doi.org/10.1007/s10915-017-0388-9 doi: 10.1007/s10915-017-0388-9
![]() |
[18] |
X. M. Gu, Y. L. Zhao, X. L. Zhao, B. Carpentieri, Y. Y. Huang, A note on parallel preconditioning for the all-at-once solution of Riesz fractional diffusion equations, Numer. Math-Theory. Me., 14 (2021): 893–919. https://doi.org/10.48550/arXiv.2003.07020 doi: 10.48550/arXiv.2003.07020
![]() |
[19] |
X. Chen, D. Ding, S. L. Lei, W. Wang, A fast preconditioned iterative method for two-dimensional options pricing under fractional differential models. Comput. Math. Appl., 79 (2020): 440–456. https://doi.org/10.1016/j.camwa.2019.07.010 doi: 10.1016/j.camwa.2019.07.010
![]() |
[20] |
S. L. Lei, W. Wang, X. Chen, D. Ding, A fast preconditioned penalty method for American options pricing under regime-switching tempered fractional diffusion models, J. Sci. Comput., 75 (2018): 1633–1655. https://doi.org/10.1007/s10915-017-0602-9 doi: 10.1007/s10915-017-0602-9
![]() |
[21] |
W. Wang, X. Chen, D. Ding, S. L. Lei, Circulant preconditioning technique for barrier options pricing under fractional diffusion models, Int. J. Comput. Math., 92 (2015): 2596–2614. http://dx.doi.org/10.1080/00207160.2015.1077948 doi: 10.1080/00207160.2015.1077948
![]() |
[22] |
X. Chen, W. Wang, D. Ding, S. L. Lei, A fast preconditioned policy iteration method for solving the tempered fractional HJB equation governing American options valuation, Comput. Math. Appl., 73 (2017): 1932–1944. http://dx.doi.org/10.1016/j.camwa.2017.02.040 doi: 10.1016/j.camwa.2017.02.040
![]() |
[23] |
S. G. Kou, A jump-diffusion model for option pricing, Manage. Sci., 48 (2002), 1086–1101. https://doi.org/10.1287/mnsc.48.8.1086.166 doi: 10.1287/mnsc.48.8.1086.166
![]() |
[24] |
B. F. Nielsen, O. Skavhaug, A. Tveito, Penalty and front-fixing methods for the numerical solution of American option problems, J. Comput. Financ., 4 (2002), 69–97. https://doi.org/10.21314/jcf.2002.084 doi: 10.21314/jcf.2002.084
![]() |
[25] |
S. Chen, F. Liu, X. Jiang, I. Turner, V. Anh, A fast semi-implicit difference method for a nonlinear two-sided space-fractional diffusion equation with variable diffusivity coefficients, Appl. Math. Comput., 257 (2015), 591–601. https://doi.org/10.1016/j.amc.2014.08.031 doi: 10.1016/j.amc.2014.08.031
![]() |
[26] |
H. Hejazi, T. Moroney, F. Liu, Stability and convergence of a finite volume method for the space fractional advection-dispersion equation, J. Comput. Appl. Math., 255 (2014), 684–697. https://doi.org/10.1016/j.cam.2013.06.039 doi: 10.1016/j.cam.2013.06.039
![]() |
[27] |
O. Marom, E. Momoniat, A comparison of numerical solutions of fractional diffusion models in finance, Nonlinear. Anal-Real., 10 (2009), 3435–3442. https://doi.org/10.1016/j.nonrwa.2008.10.066 doi: 10.1016/j.nonrwa.2008.10.066
![]() |
[28] |
N. Cai, S. G. Kou, Option pricing under a mixed-exponential jump diffusion model, Manage. Sci., 57 (2011), 2067–2081. https://doi.org/10.1287/mnsc.1110.1393 doi: 10.1287/mnsc.1110.1393
![]() |
[29] |
Y. L. Zhao, P. Y. Zhu, X. M. Gu, X. L. Zhao, H. Y. Jian, A preconditioning technique for all-at-once system from the nonlinear tempered fractional diffusion equation, J. Sci. Comput., 83 (2020): 10. https://doi.org/10.1007/s10915-020-01193-1 doi: 10.1007/s10915-020-01193-1
![]() |
[30] |
S. L. Lei, H. W. Sun, A circulant preconditioner for fractional diffusion equations, J. Comput. Phys., 242 (2013), 715–725. https://doi.org/10.1016/j.jcp.2013.02.025 doi: 10.1016/j.jcp.2013.02.025
![]() |
1. | Xu Chen, Xin-Xin Gong, Youfa Sun, Siu-Long Lei, A Preconditioned Policy–Krylov Subspace Method for Fractional Partial Integro-Differential HJB Equations in Finance, 2024, 8, 2504-3110, 316, 10.3390/fractalfract8060316 | |
2. | Beng Feng, Congyin Fan, American call option pricing under the KoBoL model with Poisson jumps, 2025, 20, 1556-1801, 143, 10.3934/nhm.2025009 |
GE | CGNR | PCGNR | ||||||
M | Time(s) | Err | Ite−In | Time(s) | Err | Ite−In | Time(s) | Err |
25+1 | 33.6213 | 0.0636 | 49.7512 | 0.8722 | 0.0708 | 4.8703 | 0.7531 | 0.0901 |
26+1 | 137.7118 | 0.0249 | 49.7512 | 0.9842 | 0.0290 | 4.8702 | 0.8753 | 0.0352 |
27+1 | 611.6442 | 0.0108 | 45.8657 | 2.8870 | 0.0120 | 5.5608 | 1.1039 | 0.0101 |
28+1 | 2469.0242 | 0.0044 | 44.4080 | 7.3217 | 0.0058 | 6.2370 | 4.5179 | 0.0046 |
29+1 | 18258.5597 | 0.0019 | 44.3532 | 43.3345 | 0.0028 | 6.2371 | 16.9851 | 0.0020 |
210+1 | ** | ** | 45.8301 | 223.1818 | 0.0013 | 6.8123 | 18.9329 | 0.0016 |
Number of time steps | Error | Order |
100 | 0.2005 | – |
200 | 0.1226 | 0.7098 |
300 | 0.0886 | 0.8015 |
400 | 0.0691 | 0.8653 |
Number of spatial steps | Error | Order |
33 | 0.0927 | – |
65 | 0.0381 | 1.2828 |
128 | 0.0126 | 1.5964 |
257 | 0.0048 | 1.3923 |
GE | CGNR | PCGNR | ||||||
M | Time(s) | Err | Ite−In | Time(s) | Err | Ite−In | Time(s) | Err |
25+1 | 33.6213 | 0.0636 | 49.7512 | 0.8722 | 0.0708 | 4.8703 | 0.7531 | 0.0901 |
26+1 | 137.7118 | 0.0249 | 49.7512 | 0.9842 | 0.0290 | 4.8702 | 0.8753 | 0.0352 |
27+1 | 611.6442 | 0.0108 | 45.8657 | 2.8870 | 0.0120 | 5.5608 | 1.1039 | 0.0101 |
28+1 | 2469.0242 | 0.0044 | 44.4080 | 7.3217 | 0.0058 | 6.2370 | 4.5179 | 0.0046 |
29+1 | 18258.5597 | 0.0019 | 44.3532 | 43.3345 | 0.0028 | 6.2371 | 16.9851 | 0.0020 |
210+1 | ** | ** | 45.8301 | 223.1818 | 0.0013 | 6.8123 | 18.9329 | 0.0016 |
Number of time steps | Error | Order |
100 | 0.2005 | – |
200 | 0.1226 | 0.7098 |
300 | 0.0886 | 0.8015 |
400 | 0.0691 | 0.8653 |
Number of spatial steps | Error | Order |
33 | 0.0927 | – |
65 | 0.0381 | 1.2828 |
128 | 0.0126 | 1.5964 |
257 | 0.0048 | 1.3923 |