Urea and sucrose assisted combustion synthesis of LiFePO4/C nano-powder for lithium-ion battery cathode application

1.   Introduction

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$ \acute{e} $vy 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$ \acute{e} $vy-$ \alpha $-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 $ \alpha $-order spatial derivative, where $ \alpha\in (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$ \acute{e} $vy-$ \alpha $-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 $ \mathcal{O}(M^3) $ and storage space of $ \mathcal{O}(M^2) $, 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 $ \mathcal{O}(M^3) $ to $ \mathcal{O}(M\log M) $, and the storage space from $ \mathcal{O}(M^2) $ to $ \mathcal{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.

2.   Mathematical formulation

In this section, the L$ \acute{e} $vy-$ \alpha $-stable process with jumps will be introduced in detail, and the pricing of stock loans underneath will be formulated.

2.1.   The model

The L$ \acute{e} $vy-$ \alpha $-stable process with jumps is defined on a complete filtered probability space $ (\Omega, \mathcal{F}_t, \mathbb{P}) $ with the current time $ t\in[0, T] $. Following the assumptions in ^[9], under this model, the logarithmic price of the underlying, i.e., $ x_t = \log(S_t) $, satisfies the following stochastic differential equation (SDE)

where $ r $, $ D $ and $ t $ are the risk free interest rate, the dividend and the current time respectively, and $ \nu = -\sigma^{\alpha}\sec\frac{\alpha\pi}{2} $ is a convexity adjustment. $ L_t^{\alpha, -1} $ is the maximally skewed L$ \acute{e} $vy-$ \alpha $-stable process with tail index $ \alpha\in (1, 2] $. $ \mathcal{N}_t $ is a Poisson process characterized by the jump intensity $ \xi\geq0 $. $ \{Y_i, i = 1, 2, \ldots\} $ is a sequence of independent and identically distributed hyper-exponential random variables with probability density function

Note that $ \hat{p}_i $ ($ i = 1, 2, ..., m_1 $) and $ \tilde{p}_i $ ($ i = 1, 2, ..., m_2 $) denote the probabilities of the ith positive and negative jumps, respectively. They satisfies $ \sum_{i = 1}^{m_1}\hat{p_i}+\sum_{j = 1}^{m_2}\tilde{p_j} = 1 $. $ \hat{\theta}_i > 1 $ ($ i = 1, ..., m_1 $) are the magnitudes of the upward jumps and $ \tilde{\theta}_j > 0 $ ($ j = 1, ..., m_2 $) are that of the downward random jumps. The average jump size is given by

where $ \mathbb{E}_{\mathbb{Q}} $ is the expectation operator under $ \mathbb{Q} $. And according to the results in References^[11] and ^[23], the asset price process is a martingale under $ \mathbb{Q} $.

2.2.   The FPIDE system

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^{\gamma t} $, where $ K $ is the principal and $ \gamma $ ($ \gamma \geq r $) is the continuously compounded loan interest rate. The payoff of the stock loan at maturity is given by

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 $ S_f $, is the optimal redemption price. According to the non-arbitrage principle, the value of the stock loan at $ t $ in the continuous region is

According to the derivations in ^[9], $ V(x, t) $ satisfies the following FPIDE in the continuous region: for $ x\in (-\infty, x_f) $, $ t\in[0, T) $

where $ x_f = \ln S_f $, $ 1 < \alpha\leq 2 $, and the left-sided Riemann-Liouville fractional derivative $ {}_{-\infty}D^{\alpha}_x $ is defined as

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.,

Therefore, by taking into consideration the boundary condition at $ x = -\infty $ ($ S = 0 $) and the terminal condition at $ T $, the price of the stock loan under the L$ \acute{e} $vy-$ \alpha $-stable process with jumps satisfies

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.

3.   Numerical scheme

Upon establishing a closed pricing system Eq (2.4) for stock loans when the underlying satisfies the L$ \acute{e} $vy-$ \alpha $-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.

3.1.   Coordinate transformation

To simplify the solution process, we use the following new variables

With the transformation details provided in Appendix A, the pricing system Eq (2.4) becomes: for $ z\in(-\infty, z_f] $, $ t\in[0, T] $

where $ a = (r-\nu-\xi\varsigma-D-\gamma) < 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 $ e^{z_f} $. Since the American call option value is always greater than or equal to its payoff value, we directly have $ U(z, t)\geq\max\left(e^z-K, 0\right) $.

3.2.   The penalty method

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

where $ \varepsilon $ is a regularization constant and $ 0 < \varepsilon \ll 1 $, $ H $ is a constant, and $ q(z) $ is defined as $ q(z) = e^z - K $. With the penalty term added on and the truncation of the domain applied, the FPIDE system becomes: for $ z\in[z_{min}, z_{max}] $, $ t\in[0, T] $

where $ \exp(z_{max}) $ and $ \exp(z_{min}) $ denote the maximum and minimum stock prices, respectively. In particular, according to the reasons provided in ^[16] and ^[25], we set $ z_{max} = \ln(3K) $ and $ z_{min} = \ln(0.01) $. In the following, for simplicity purposes, the subscript $ \varepsilon $ of $ U_{\varepsilon}(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 $ [z_{min}, z_{max}] $. 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 $ \varepsilon\rightarrow 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 $ \Delta z = \frac{z_{max}}{M} $ and $ \Delta t = \frac{T}{N} $ for the spatial and time directions, respectively. Denote the nodal value in the $ z $ direction as $ z_j = (j-1)\Delta z $, for $ j = -(M+1), \cdots, -1, 0, 1, \cdots, (M+1) $, and in the $ t $ direction as $ t_i = (i-1)\Delta t $, for $ i = 1, 2, \cdots, (N+1) $.

Next, the Shifted Grünwald-Letnikov formula ^[21] is employed to approximate the left-sided Riemann-Liouville fractional derivative, i.e.,

where $ g_k $ are defined as

In addition, the integral term contained in the governing equation of Eq (3.3) is approximated by the trapezoidal rules ^[9], i.e.,

where

and

Therefore, after the full implicit difference scheme is applied, the FPIDE system Eq (3.3) becomes

with

where $ U_j^i $ denotes the approximation solution of $ U(x, t) $ at $ (x_j, t_i) $.

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 $ \rho_j^M $ in Eq (3.4) and $ R_j $ in Eq (3.5) are bounded and satisfy

Proof. Since $ f_Y(y) $ is the density of a hyper-exponential random variable $ Y $, we have

Therefore, it is clear that

On the other hand, since $ \hat{p}_i\geq0, \hat{\theta}_i > 1 $ and $ \exp(z_{\max}) = 3K $, it is straightforward to obtain

The proof of this lemma is thus completed.

Lemma 3.2. For $ \alpha\in(1, 2] $, $ g_k $ ($ k = 0, 1, \cdots, +\infty $) satisfy

where $ m\geq1 $.

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 $ \Delta t $ satisfies $ \Delta t \leq \frac{1}{|\gamma -r-\xi -1|} $, and the constant $ H $ is bounded below as

where $ \bar{K} = 2[\exp(z_{\max})-K] $, the approximated stock loan values $ \{U_j^i\} $ satisfy

Proof. We shall first show that $ U_j^i\geq q_j $. Define $ u^i_j = U^i_j - q_j $. It is straightforward that $ u^{N+1}_j = U^{N+1}_j-q_j\geq0 $. Substituting $ u_j^i $ into the governing equation contained in Eq (3.3) yields

where

Since $ |\frac{e^{\Delta z}-1}{\Delta z}| \leq \frac{e^{z_{\max}}-1}{z_{\max}}\leq 1 $ and $ \sum^\infty_{k = 0}g_k = 0 $, we have

In addition, according to ^[27], we have $ \sum\limits^\infty_{k = 0}g_k e^{z_{j-k+1}} = e^{z_{j+1}}\sum\limits^\infty_{k = 0}g_k e^{-k\Delta z} $ and $ (1-z)^{\alpha} = \sum\limits^\infty_{k = 0}g_k z^{-k} $ for $ |z|\leq 1 $. It is clear that

Now, let $ \bar{K} = 2[\exp(z_{\max})-K] $, and we have

Therefore, combining with Lemma 3.2, it is straightforward to obtain

Define

and let $ J $ be an index such that $ u_J^{i} = u^i $. Since $ a < 0 $, $ \rho^M_j > 0 $ for all $ j $, and $ g_k > 0 $ for $ k\neq1 $, one can deduce from Eq (3.6) that

After simple algebraic operations, we obtain

On the other hand, according to Eq (3.7) and $ \Delta t\leq\frac{1}{|\gamma-r-\xi -1|} $, it is straightforward that

For the ease of description, let

Define a function $ f(x) $ as

According to the definition of $ f(x) $, it is clear that $ f(u^i)\geq0 $ if $ u^{i+1}\geq 0 $. Since $ f(0) = \Delta t (E_j - H)\leq 0 $, $ f'(x) = A+ \frac{\epsilon \Delta t C}{(x+\epsilon)^2} > 0 $, and $ u^{N+1}_j\geq0 $, we obtain $ u^i\geq0 $, and consequently $ u^i_j\geq0 $. Therefore, $ U^i_j\geq q_j $ is satisfied.

Next, we show that $ U^i_j\geq 0 $. Similarly, define $ U^i = \min_j U^i_j $, and let $ J $ be an index such that $ U^i_J = U^i $. From Eq (3.6), it is clear that

which, after a careful rearrangement, yields

On the other hand, since $ R_j > 0 $ and $ \frac{\Delta t \varepsilon H}{U^i_J-\varepsilon + q_J} > 0 $ because of the fact that $ U_j^i\geq q_j $ for all $ i, j $, we have

Since $ 1-\left(\gamma-r-\xi - 2\sum_{\ell = 0}^M\rho^M_{\ell-j} \right)\Delta t > 0 $, one can deduce that

Since $ U^{N+1}\geq0 $ due to the fact that $ U^{N+1}_j = \max\left[\exp(z_j)-K, 0\right]\geq0 $ for all $ j $, we finally have $ U_j^i\geq0 $. Thus, the proof is thus completed.

We remark that the constraint condition we set on $ \Delta t $, i.e., $ \Delta t \leq \frac{1}{|\gamma -r-\xi -1|} $ is indeed not harsh. This is because the parameters $ \gamma $ and $ r $ satisfy $ \gamma \in (0, 1) $ and $ r \in (0, 1) $, respectively, and moveover the value of $ \xi \leq 10 $ is also reasonable according to the discussions provided in ^[28]. Financially, Theorem 1 ensures that the discrete value $ U_j^i $ cannot fall below the corresponding intrinsic value, a characteristic inherited from the early exercise nature of the stock loan.

4.   The PCGNR method with a circulant pre-conditioner

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 $ \Delta z = \frac{\left(z_{max}-z_{min}\right)}{M} $ and $ z_j = z_{min}+(j-1)\Delta z $, for $ j = 1, 2, \cdots, M+1 $. Define

and

Then the matrix form of Eq (3.6) can now be written as

where $ {\bf U}^i = (U^i_2, U^i_2, ..., U^i_M)^{\tau} $, $ F({\bf U}^i) = (f(U^i_2), f(U^i_3), ..., f(U^i_M))^{\tau} $,

$ \mathcal{I}_{(M-1)\times {(M-1)}} $ is the identity matrix,

$ {\bf R}^i = \left(R_2^i, R_3^i, ..., R_{M}^i\right)^{\tau} $, $ {\bf E}^i = \left(0, 0, ..., \varrho\right)_{M-1}^{\tau} $ and $ \varrho = \left(-\eta_2g_0-\sum_{\ell = 0}^{M-1}\rho_{\ell}\right)U^i_{M+1} $.

After a Newton-type iteration approach is applied, Eq (4.3) becomes

where $ l = 1, 2, \cdots $, the initial guess $ \omega^0 $ is set as $ \omega^0 = \textbf{U}^{i+1} $ for each time level, and $ \delta\omega^l = \omega^l-\omega^{l-1} $. Moreover, $ \mathcal{J}_{\textbf{F}} $ is Jacobian matrix with column vector $ F(\omega^l) $, and $ \kappa \in (0, 1) $ is the damping parameter. We set $ \textbf{U}^i = \omega^l $ once the stopping criterion $ \parallel \omega^l - \omega^{l-1}\parallel_{\infty} \leq tol $ for some $ l $ is reached, where $ tol $ is the stopping tolerance of the iterative method. By denoting

Equation (4.4) can be rewritten as

The most challenging part in solving Eq (4.4) is the high computational cost resulting from the fact that both $ \mathcal{A} $ and $ \mathcal{W} $ are dense matrices. To overcome this difficulty, we first apply the CGNR method ^[30], which is in fact to solve $ [{\bf M}-\mathcal{J}_{\textbf{F}}]^{\tau}{\bf M}\delta \omega^{l} = [{\bf M}-\mathcal{J}_{\textbf{F}}]^{\tau} {\bf b}^{l} $ 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 $ {\bf M}^{\mathcal{T}}{\bf M} $ 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 $ \mathcal{J}_{\textbf{F}} $ is not Toeplitz matrix and we should approximate this matrix as $ a_0 \mathcal{I} $, where $ a_0 $ is the average value of main diagonal elements of matrix $ \mathcal{J}_{\textbf{F}} $. So we structure a Toeplitz matrix as following

Next, the Strang's circulant preconditioner $ s(\mathcal{T}) = [s_{j-k}]_{0\leq j, k < M} $ for matrix $ \mathcal{T} $ is structured as

We use $ \mathcal{P} $ denotes the Strang's circulant preconditioner $ s(\mathcal{T}) = [s_{j-k}]_{0\leq j, k < M} $ to simplify the expression.

Mathematically, after the PCGNR method with a pre-conditioner $ \mathcal{P} $ is applied, Eq (4.4) becomes

The pseudo-code of the PCGNR method is displayed in Algorithm 1. The matrix-vector multiplication only needs $ \mathcal{O}(MlogM) $ operations via the fast Fourier transform (FFT) method.

5.   Numerical experiments and discussions

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 $ \kappa $ and the stopping tolerance $ tol $ are set to be $ 0.2 $ and $ 10^{-6} $, respectively.

5.1.   Performance of the numerical method

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)\geq\max(e^z - K, 0) $ discretely, as shown in Theorem 3.1. In Figure 1, the surfaces of the difference between $ U^i_j $ and $ \max(q_j, 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 $ e^{z_f} $. In Figure 2(a), we display $ U(z, t) $ against the stock price $ e^z $. 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 $ e^{z_f} $ as a function of the time to expiry is plotted. One can clearly observe from this figure that $ e^{z_f} $ 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 $ \alpha = 1.52 $, $ \sigma = 0.2 $, $ r = 0.05 $, $ D = 0.06 $, $ K = 2 $, $ T = 0.2 $, $ \xi = 0.03 $, $ \hat{p} = 0.04 $, $ \gamma = 0.06 $, $ \hat{\theta} = 1.2 $, $ \tilde{\theta} = 0.2, m_1 = 1, m_2 = 1, \hat{p} = 0.5, \tilde{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 = \parallel U-u(z, t)_{\mathrm{ref}}\parallel_2 $, where $ \parallel \cdot\parallel_2 $ is the $ L_2 $ norm for matrix, and $ u(z, t)_{\mathrm{ref}} $ is the benchmark solution was determined directly through matrix operation '$ A\backslash b $' in Matlab with $ M = 2^{11} $ and $ N = 500 $.

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 $ \Delta z = \frac{z_{max}-z_{min}}{2^{11}} $, 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

where $ i $ denotes the $ i $th of Table 2, and $ L_i $ ($ i = 2, 3, 4 $) is the number of time steps displayed in the $ i $th 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.

5.2.   Quantitative analysis

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 $ e^{z_f} $ back to $ V $ and $ S_f $ via $ V(x, t) = e^{\gamma t}U(z, t) $ and $ S_f(t) = e^{\gamma t+z_f} $, 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 $ \xi $. One can observe clearly from this figure that the optimal redemption price increases w.r.t $ \xi $. 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 $ S_f $ w.r.t $ \xi $ holds automatically. This figure also reveals the fact that $ S_f $ 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 $ \hat{p} $. From this figure, it is clear that a larger $ \hat{p} $ results in a lower optimal redemption price. In fact, the logarithmic return of the stock price is decreasing w.r.t $ \hat{p} $, because the return decreases w.r.t $ \xi $ from Eq (2.1) and $ \xi $ increases w.r.t $ \hat{p} $ from Eq (2.2). Therefore, an increasing $ \hat{p} $ would lower the stock price, and thus makes the intermediate American call option $ U(z, t) $ less valuable. Therefore, the optimal exercise price $ e^{z_f} $ of the intermediate American call decreases w.r.t $ \hat{p} $. Since the optimal redemption price of the stock loan $ S_f $ is related to $ z_f $ through $ S_f(t) = e^{\gamma t+z_f} $, it is straightforward that a larger $ \hat{p} $ results in a lower optimal redemption price. Similarly, one could explain the monotonicity of the optimal redemption price w.r.t $ \hat{\theta} $ and $ \tilde{\theta} $. For the length of the paper, we provide those figures in Figure 5 with no detailed explanations.

Finally, we examine how the tail index $ \alpha $ influences the optimal redemption strategy. Several sets of optimal redemption prices with different $ \alpha $ 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 $ \alpha $. Financially, the tail index controls the tail of the distribution of the underlying price, and both tails will be fatter when $ \alpha $ becomes smaller, as suggested in ^[1]. Therefore, as $ \alpha $ becomes smaller, the possibility of larger stock prices increases, and so does the optimal redemption price.

6.   Conclusions

In this paper, the pricing of stock loans under the L$ \acute{e} $vy-$ \alpha $-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.

Acknowledgments

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).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A

Since

we have

Let $ s=y - \gamma t $, and we have

Now, substituting Eqs (A.1)–(A.3) into the pricing system Eq (2.4), the FPIDE system Eq (3.1) can finally be obtained.