Numerical Method for a System of PIDEs Arising in American Contingent Claims under FMLS Model with Jump Diffusion and Regime-Switching Process

This paper investigates a numerical method for solving fractional partial integro-diﬀerential equations (FPIDEs) arising in American Contingent Claims, which follow ﬁnite moment log-stable process (FMLS) with jump diﬀusion and regime switching. Mathematically, the prices of American Contingent Claims satisfy a system of d problems with free-boundary values, where d is the number of regimes of the market. In addition, an optimal exercise boundary is needed to setup with each regime. Therefore, a fully implicit scheme based on the penalty term is arranged. In the end, numerical examples are carried out to verify the obtained theoretical results, and the impacts of state variables in our model on the optimal exercise boundary of American Contingent Claims are analyzed.


Introduction
More and more researchers pay much attention to the pricing problem of American Contingent Claims since the pricing model is free-boundary problem and it is of great academic value to solve this kind of model. Hence, we consider the stock loan, which is a kind of famous American contingent claim, in this paper. In fact, Xia and Zhou had studied the stock loan with infinite maturity based on the Black-Scholes framework in 2007 [1]. Following this contribution, many scholars also consider the same topic by the Black-Scholes model under other different conditions. For instance, Chen et al. studied the stock loan with finite maturity under the framework that risk-free interest rate follows the Rednleman-Bartter model without drift term [2]. Cai et al. investigated the value and optimal redemption price of stock loan with infinite and finite maturity under the framework of the hyperexponential jump diffusion model, respectively [3]. Liang et al. provided the formulae of stock loan with automatic termination clause, cap, and margin [4]. Prager and Zhang [5] considered the finite horizon loan valuation under various stock models that include classical aerometric Brownian motion, mean-reverting, and two-state regime switching with both mean-reverting and geometric Brownian motion states.
Most models in literatures use a continuous-time Brownian motion as their primary source of stochasticity. ese stock models include geometric Brownian motion with stochastic interest rate and diffusion with hyperexponential jump. However, even if considering the actual phenomenon, the continuous-time Brownian motion does not fully reflect the stochastic nature of financial markets. First of all, empirical evidence has showed that the distribution of asset return can appear "leptokurtic distribution" (see, e.g., [6][7][8]), which has a higher peak and two heavier tails compared with the standard normal distribution. Secondly, from an economic perspective, regime-switching behavior captures the changing preferences and beliefs of investors concerning stock prices as the state of financial market changes. e presence of regime-switching dynamics in financial market has been well acknowledged. A frequently observed phenomenon is that a transition in a real business cycle, expansion, and contraction usually leads to significant changes of stock returns, interest rates, and financial indices.
ese changes exhibit certain cyclic or periodic patterns. Consequently, there is a need for more realistic models that better capture asymmetric distribution, phenomenon of jumps, and the dynamic changes of risk asset. Naturally, the finite moment log-stable model with jump diffusion and regime switching (FMLSJR) [9] is employed to capture dynamical process of underlying asset price.
Under framework of FMLSJR, the values of stock loan follow a d-dimension coupled partial integral-differential equations (PIDEs). And there are d-free boundaries in the pricing model. Associated with each regime, there is a function of optimal redemption prices and a space-fractional partial integral-differential equation. And for the fractional model, the problem of how to solve it has attracted the interest of many researchers. For instance, Jumarie [10] presented the solutions of time-fractional and space-fractional BCS models about the European option by Fourier transforms and Mittag-Leffler function. Liang et al. [4] used the Laplace transform technique to obtain a completed pricing formula of option in the case of the bifractional Black-Scholes-Merton model. Following this work, Kumar et al. [11] investigated the time-fractional BS European option pricing model and obtained an analytical solution. Under the framework of the FMLS model, Chen et al. [12] considered the analytical solution of European option by using a Fourier integral transform and Fox functions. According to the literature review above, it appears that the favored methods used to consider the analytical solution of the fractional models of European option are via integral transform methods. e solutions obtained by means of these methods usually take the form of a convolution of some special functions, which make it difficult to compute. Moreover, these integral transform methods may not be effective for the American Contingent Claims. Hence, studying the numerical approximate solutions of these models appears to be a very practical and important research objective.
ere are a lot of recent publications investigating the numerical method for the fractional models of European and American Contingent Claims. For example, Cartea and del-Castillo-Negrete [13] employed the shifted Grawld-Letnik scheme to discrete the fractional derivatives and set an implicit scheme to solve the pricing model of exotic option. Marom and Momoniat [14] presented a comparison of numerical solutions of the FMLS, KoBol, and CGMY models. Li [15] considered a difference scheme with firstorder accuracy for a space-fractional BS model. By using a wavelet technique, G. Hariharan [16] gave a numerical scheme based on the wavelet technique to solve the timefractional BS model arising from European option pricing problem. Chen [17] investigated a second-order finite difference method for the one-dimensional FMLS model governing the valuation of European options and a power penalty method for a space-fractional differential linear complementarity problem arising in the valuation of American options on a single asset and then extended the above methods to the corresponding two dimensional models arising in the valuation of European and American options on two assets. e author also gave the detailed numerical analysis of the methods. Zhang et al. [18] considered the European double barrier option in the case of tempered fractional Black-Scholes equation, and they employed a fast biconjugate gradient stabilized method to solve the linear system. Fan et al. [19] set a fully implicit scheme with first-order accuracy based on the penalty function for the stock loan pricing model, and their method can be used to other pricing model of American Contingent Claims. Zhou et al. [9] investigated the iterative Laplace transform methods for system of fractional PDEs and PIDEs arising in option pricing.
However, from the existing literature, it appears that research on the numerical simulation of the coupled fractional differential equations for American Contingent Claims is also relatively limited. So, in this paper, we give a numerical scheme based on penalty term and present a technique to solve the coupled system. ere are two contributions of this paper. First of all, we prove the lower bound of stock loan values at the case of fractional order with jumps and mechanism transfer. Especially, the involvement of mechanism transfer increases the difficulty of proof. Because of the existence of mechanism transfer term, the objective equation in the model is coupled by d partial differential equations. A skill is employed to solve the nonlinear equation in this paper to avoid repeated iteration. Secondly, we present a technique to solve the final coupled system, which makes the numerical results more accurate than the results calculated by the iterative algorithm. e rest of the paper is organized as follows. In Section 2, the completed pricing model is given. And we introduce our numerical method by extending the penalty method and prove the effectiveness of this method in Section 3.
e numerical examples and analysis of parameter impact are displayed in Section 4, and we give the conclusions in the final section.

The Mathematical Model
Let (Ω, F t , P) be a complete filtering probability space, where t ∈ [0, T] and T > 0 is a fixed time expiration. Assume that there exists an equivalent martingale measure Q, under which the dynamics of the logarithmic price of the risky asset x t � log(S t ) are given by the following stochastic differential equation as in [9]: where χ(t) is a continuous-time Markov chain with d− states χ k , k ∈ D: � 1, 2, . . . , d { }. Assume that at each state χ k , the risk-free interest rate r(χ k ) � r k , σ(χ k ) � σ k > 0, and the convexity adjustment Discrete Dynamics in Nature and Society is the maximally skewed Lévy process, and the tail index α ∈ (1, 2). N t is a Poisson process, which is characterized by a jump intensity parameter ξ ≥ 0. Y i , i � 1, 2, . . . is a sequence of independent and identically distributed hyperexponential random variables with probability density function as follows: where p i ≥ 0, i � 1, 2, . . . , m 1 and p i ≥ 0, i � 1, 2, . . . , m 2 are the probabilities of different kinds of positive and negative jumps, respectively. Moreover, it satisfies . . , m 1 , and θ j > 0, j � 1, . . . , m 2 , are magnitude parameters of different types of upwards and downward random jumps, respectively. e average jump size is given as where E Q is the expectation operator under probability measure Q. Let the following be the generator matrix of the Markov chain process: with constants q ij ≥ 0, for all i, j, and According to [9], the value V(k; x, t) of stock loans satisfies the following coupled FPIDEs: which is left-sided Riemann-Liouville fractional derivatives.
So, the boundary conditions of the mathematical model can be obtained as follows ( [2,20,21]): where S f,k (e x f,k ) denotes the optimal redeem price of stock loan contract. For every k ∈ D, the model is a free-boundary problem. Associated with each regime is an optimal exercise and a FPIDE, and it is impossible to obtain the analytic solution; therefore, the numerical method should be considered. Moreover, according to the conditions of stock loan contract, the investor should continue to hold the contract if the redemption value is less than the holding value. erefore, the function V(k; x, t) must satisfy the following inequality: for all x ≤ x f and 0 ≤ t ≤ T.

Model Transformation.
To avoid the effect of time-factor e ct in our model on the numerical results, we first introduce a new variable system as follows: Discrete Dynamics in Nature and Society By variable transformation (the details are displayed in Appendix A), the pricing systems (6) and (8) can be written as where . And, inequality (9) should be rewritten as for all z ≤ z f,k and 0 ≤ t ≤ T.
It is straightforward to obtain that the function U(k; z, t) in system (11)-(15) can be viewed as an American call option with strike price K and free-boundary e z f,k . e analytic solution of U(k; z, t) is laborious and even impossible to achieve; therefore, we next employ a finite difference scheme based on the penalty function to solve it. Now, we extend the penalty method to model (11)-(15) and develop an efficient numerical method via the implementation of a fully implicit scheme. So, we let 0 < ϵ ≪ 1 be a small regularization parameter and C > 0 is a constant. By adding the penalty terms [22], to the corresponding equations in (11), we obtain the following system of nonlinear space-fractional partial integraldifferential equations on a fixed domain: where e z max denotes the maximum stock price, z ∈ (− ∞, z max ], and t ∈ [0, T]. According to the numerical experiments conducted by Kandilarov and Valkov [23], the maximum value of stock is equal to 3 or 4 times of strike price, and we take e z max � 3K in this paper. In addition, we will omit the subscript ϵ of function U ϵ (k; z, t) for convenience.
For the integral term, the trapezoidal formula is used to discrete it at grid point (z j , t i ) (see [9]) as follows: In addition, the left-hand Riemann-Liouville fractional derivative can be approximated by the first-order Grünwald-Letnikov formula as follows (see [12,24]): where g ι are the fractional difference coefficients given as follows: Using the discretization and applying the fully implicit method, it leads to a nonlinear algebraic equation: with the boundary and terminal conditions as follows: For our proposed difference scheme, we have a discrete analogue form of the important property inherited by the stock loan model. Prior to presenting the proof, we need the following lemma of [24].

Lemma 2.
Both the coefficients ρ M j in equation (24) and R i,j k in equation (25) are bounded, and Proof. From the fact that f Y (y) is the density function of hyperexponential random variables Y, it follows that In addition, p i ≥ 0, θ i > 1 and exp(z max ) � 3K or 4K, so it is straightforward to obtain □ Discrete Dynamics in Nature and Society
Proof. First, we prove U i,j for all i, j, and k. It is straightforward to obtain u N+1,j k � U N+1,j k − q j ≥ 0. Hence, by substituting u i,j k into (28) and after simplification, it yields where From the fact ∞ ι�0 g ι e z j− ι+1 � e z j+1 ∞ ι�0 g ι e − ιΔz and according to [14], (1 − z) α � ∞ ι�0 g ι z − ι when |z| ≤ 1. erefore, we obtain that Discrete Dynamics in Nature and Society To sum up, combining Lemma 2, it is straightforward to obtain Define and let (k 0 , j 0 ) be a pair of indices, such that u i � u i,j 0 k 0 . So, from (28), it follows that To reduce the form of inequality (43), we obtain that On the other hand, Hence, it is straightforward to obtain that Define a function So, if u i+1 ≥ 0, we can obtain From it follows that Based on the condition u N+1,j k ≥ 0, we can obtain u i ≥ 0. Naturally, u i,j for all i, j, and k. Next, we prove that U i,j k ≥ 0 for all i, j, and k. Following the idea above, define and let (k 0 , j 0 ) be a pair of indices, such that U i,j 0 k 0 � U i . It follows from (28) that By simple calculation, the following equality holds: From above, the inequalities U j,i k ≥ q j for all i, j, and k are verified. So, In addition, U N+1,j k � max(exp(z j ) − K, 0) ≥ 0 for all j, k. Namely, U N+1 ≥ 0. erefore, it yields that for all i, j, and k, which completes the proof. (59) Here, we take z min � ln(0.01). en, we can check that U(z min , t) � 0. In such a case, we redefine the space step Δz � (z max − z min )/M, and z j � z min +(j − 1)Δz, j � 1, 2, . . . , M + 1. (60) then the matrix form of numerical scheme (28) can be written as follows: where I denotes the (M − 1) × (M − 1) identity matrix.
In fact, system (62) is a coupling nonlinear equation so that it is very hard to obtain the numerical solution. Hence, we let 8 Discrete Dynamics in Nature and Society A � M 1 Δtq 12 I · · · Δq 1 d I where en, (62) can be written as the following nonlinear system: (68) Note that system (68) cannot be solved directly since the penalty function F(U i ) is nonlinear with respect to U i . In the following, Newton iteration approach is employed to solve where l � 1, 2, . . . and κ ∈ (0, 1) is a damping parameter. e initial value ω 0 � U i+1 for each time level is the given initial guess and δω l � ω l − ω l− 1 . J F is a Jacobian matrix with column vector F(ω l ). We choose U i � ω l . e condition ‖ω l − ω l− 1 ‖ ∞ ≤ θ is the stopping criterion, where the positive control θ is a sufficiently small tolerance number. In this paper, we take κ � 0.2, θ � 10 − 4 , and ε � 10 − 5 .

Numerical Examples and Discussion
4.1. Discussion on the Numerical Method. In this part, we provide some numerical examples to support the theoretical results that we have obtained in this paper. e model parameters are set to be r 1 � 0.05, r 2 � 0.04, σ 1 � 0.20, σ 2 � 0.24, D � 0.8, c � 0.6, ξ � 0.01, p � 0.02, and θ � 1.2. All contracts of stock loan have maturity T � 2 years and principal K � 20. First of all, to verify the conclusion that inequality U i,j k ≥ max(q j , 0) holds in eorem 1, we show the mesh surface of U i,j k − max(q j , 0) in Figure 1 for different time and stock price in two states. As shown in Figure 1, it is straightforward to conclude that the present difference scheme conserves the inequality U i,j k − max(q j , 0) ≥ 0 for all i, j, and k. In Figure 2, the values of U(k; x, t) as a function of both stock price S t and time t over the rectangular domain [S min , 60] × [0, 2] under difference economic state are displayed. It can be observed from Figure 2 that our numerical method produces smooth and stable approximation solutions.
In fact, the function U(k; z, t) in models (18)- (21) can be recognized as an American call option with strike K and optimal exercise boundary, exp(z f,k ). en, the values of U(k; z, t) are not less than the payoff function max(e z − K, 0), and the function exp(z f,k ) is increasing with respect to the duration in every state. Figures 3 and 4 show these two facts, respectively. In Ref. [9], for solving the coupled equations resulted by regime switch, the authors proposed the iteration algorithm. Hence, we use this idea to solve system (62), and the consumed CPU time to run our method and the iteration algorithm for different k and M and fixed N � 50 is displayed in Table 1. All numerical computations were carried out by using MATLAB on the Lenovo S5 laptop with configuration as follows: Intel(R) Core(TM) i5-6300HQ, 2.30 GHz. We solve the linear systems by Matlab operator A∖b. In Table 1, CPU − time 1 and CPU − time 2 denote consumed CPU time of the iteration algorithm and our method, respectively. According to the results in Table 1, we can obtain that our method is more effective, and the advantages of our method become more obvious with the increase in k.

Impact of Parameters.
In this subsection, we discuss some interesting implications by considering the two-state regime-switching model in pricing stock loans. Moreover, we compare the result with the optimal redemption prices obtained by using the framework without regime switching.
In Figure 5, there are four different curves of optimal redemption prices under different parameters, and the horizontal axis denotes the time to expiration. e black and blue curves are obtained by the models without regime switching as (Q � 0). In the case of r 1 � 0.05 or r 1 � 0.04, the optimal redemption prices calculated by our model are bigger than those obtained by the model without regime switching, which means the investors should choose a higher price to redeem stock when there is uncertainty in the financial market. However, as the maturity closes to zero, the uncertainty should be disappeared such that the holder should exercise the stock loan contract at the price level of Ke cT as shown in Figure 5.
Discrete Dynamics in Nature and Society As shown in Figure 6, the introduction of regime switching produces higher optimal redemption prices than those calculated by the framework without regime switching (Q � 0).
is phenomenon is the same as an American option, of which the optimal exercise values are monotonically increasing functions of volatility. e difference between blue line and green line (red line and black line) reflects the extra values of having a certain positive probability that the stock should stay at the economy state with higher risk. However, as the maturity of stock loan contract approaches, the expected amount of time spent on other economy state also decreases, which can result in the extra value due to decreasing regime switching. In addition, as time approaches to the expiration, the investor should redeem stock at price level of Ke cT as shown in Figure 6.

Conclusions
In this paper, we investigate the pricing model of stock loan under the FMLS model with jump diffusion process in a d-state regime-switching economy, which was formulated as a d-dimension coupling fractional partial integro-differential equations. We first introduce a penalty term to transform the original model into one with fixed domain in every state. Secondly, a fully nonlinear implicit difference scheme is proposed, and we prove the numerical solution following the constraint V(k, x, t) ≥ max(e x − Ke ct , 0) inherited by the stock loan. Moreover, our numerical examples show an important phenomenon that the uncertainty coming from the regime switching should raise the optimal redemption prices of stock loans. And the results in this paper can be used to do research studies on other types of American Contingent Claims.      In fact, during the process of solving the coupling linear equation, the values of objective function in different economical state are calculated simultaneously rather than obtained through iteration, respectively; therefore, our method is more efficient. Under the same accuracy, our method is more effective than the iteration algorithm as shown the results of simulation in Table 1.
In future research, first of all, one could expand the number of underlying assets, which is also meaningful. Secondly, the investor can terminate the contract by paying liquidated damages. Finally, the fractional derivatives and the coupled term usually result in a full or dense coefficient matrix, which has significant computational and storage requirements; therefore, one could improve the computing speed of our method.