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

Authors: Congyin Fan; Peimin Chen
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

This paper investigates a numerical method for solving fractional partial integro-differential equations (FPIDEs) arising in American Contingent Claims, which follow finite moment log-stable process (FMLS) with jump diffusion and regime switching. Mathematically, the prices of American Contingent Claims satisfy a system ofd 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.


---

## Body

## 1. 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. These 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–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. The 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. These 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 ad-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. The 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.There 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 first-order 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 time-fractional 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. The 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.There 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 byd 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.The rest of the paper is organized as follows. InSection 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. The numerical examples and analysis of parameter impact are displayed in Section 4, and we give the conclusions in the final section.


## 2. The Mathematical Model

LetΩ,ℱt,ℙ 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 ℚ, under which the dynamics of the logarithmic price of the risky asset xt=logSt are given by the following stochastic differential equation as in [9]:(1)dxt=rχt−νχt−ξςdt+σχtdLtα,−1+d∑i=1NtYi,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=rk, σχk=σk>0, and the convexity adjustment νχk=νk=−σkαsecαπ/2. The term Ltα,−1 is the maximally skewed Lévy process, and the tail index α∈1, 2. Nt is a Poisson process, which is characterized by a jump intensity parameter ξ≥0. Yi,i=1,2,… is a sequence of independent and identically distributed hyperexponential random variables with probability density function as follows:(2)fYy=∑i=1m1p^iθ^ie−θ^iy1y≥0+∑j=1m2p˜jθ˜je−θ˜jy1y≤0,where p^i≥0, i=1,2,…,m1 and p˜i≥0, i=1,2,…,m2 are the probabilities of different kinds of positive and negative jumps, respectively. Moreover, it satisfies ∑i=1m1p^i+∑j=1m2p˜j=1. Similarly, the parameters θ^i>1, i=1,…,m1, and θ˜j>0, j=1,…,m2, are magnitude parameters of different types of upwards and downward random jumps, respectively. The average jump size is given as(3)ς=EℚexpY1−1=∑i=1m1p^iθ^iθ^i−1+∑j=1m2p˜jθ˜jθ^j+1−1,where Eℚ is the expectation operator under probability measure ℚ. Let the following be the generator matrix of the Markov chain process:(4)Q=−q11q12q13⋯q1dq21−q22q23⋯q2d⋮⋱⋱⋱⋮qd−1,1qd−1,2⋯−qd−1,d−1qd−1,dqd1qd2⋯qd,d−1−qd,d,with constants qij≥0, for all i,j, and(5)∑j=1,j≠kdqkj=qkk,k∈D.According to [9], the value Vk;x,t of stock loans satisfies the following coupled FPIDEs:(6)∂Vk;x,t∂t+rk−νk−D−ξς∂Vk;x,t∂x+ξ∫−∞+∞Vk;x+y,tfYydy+νkDxα−∞Vk;x,t+∑l=1,l≠kdqklVl;x,t=rk+qkk+ξVk;x,t,where 1<α<2, x∈−∞,xf. And,(7)Dxα−∞Vk;x,t=1Γ2−α∂2∂x2∫−∞xVk;z,tx−zα−1dz,which is left-sided Riemann–Liouville fractional derivatives.So, the boundary conditions of the mathematical model can be obtained as follows ([2, 20, 21]):(8)limx⟶−∞Vk;x,t=0,Vk;xf,k,t=exf,k−Keγt,∂Vk;xf,k,t∂x=Sf,k=exf,k,Vk;x,T=maxex−KeγT,0,where Sf,kexf,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. Therefore, the function Vk;x,t must satisfy the following inequality:(9)Vk;x,t≥maxex−Keγt,0,for all x≤xf and 0≤t≤T.


## 3. Fully Implicit Scheme Based on Penalty Method

### 3.1. Model Transformation

To avoid the effect of time-factoreγt in our model on the numerical results, we first introduce a new variable system as follows:(10)z=x−γt,Uk;z,t=eγtVk;x,t,zf,k=xf,k−γt.By variable transformation (the details are displayed in Appendix A), the pricing systems (6) and (8) can be written as(11)∂Uk;z,t∂t+ak∂Uk;z,t∂z+ξ∫−∞+∞Uk;z+y,tfYydy+νDzα−∞Uk;z,t+∑l=1,l≠kdqklUl;x,t=bkUk;z,t,(12)Uk;zf,t=expzf,k−K,(13)∂Uk;f,k,t∂z=expzf,k,(14)limz⟶−∞Uk;z,t=0,(15)Uk;x,T=maxez−K,0,where ak=rk−νk−D−γ−ξς and bk=rk+qkk+ξ−γ. And, inequality (9) should be rewritten as(16)Uk;z,t≥maxez−K,0,for all z≤zf,k and 0≤t≤T.It is straightforward to obtain that the functionUk;z,t in system (11)–(15) can be viewed as an American call option with strike price K and free-boundary ezf,k. The analytic solution of Uk;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],(17)εCUϵk;z,t+ϵ−qz,k∈D,qz=ez−K,to the corresponding equations in (11), we obtain the following system of nonlinear space-fractional partial integral-differential equations on a fixed domain:(18)∂Uϵk;z,t∂t+ak∂Uϵk;z,t∂z+ξ∫−∞+∞Uϵk;z+y,tfYydy+∑l=1,l≠kdqklUϵl;z,t+νkDzα−∞Uϵk;z,t+ϵCUϵk;z,t+ϵ−qz=bkUϵk;z,t,(19)Uϵk;zmax,t=ezmax−K,(20)limz⟶−∞Uϵk;z,t=0,(21)Uϵk;z,T=maxez−K,0,where ezmax denotes the maximum stock price, z∈−∞,zmax, 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 ezmax=3K in this paper. In addition, we will omit the subscript ϵ of function Uϵk;z,t for convenience.


### 3.2. Implicit Scheme

We first takeΔz>0 as spatial step, such that MΔz=zmax, where M is a positive integer and then place N+1 uniform grids in the t direction, namely, Δt=T/N. That is,(22)xj=j−1Δz,ti=i−1Δt,where j=∞,…,−1,0,1,…,M+1; i=1,2,…,N+1; and Uki,j denotes the value of function Uk;xj,τi at point xj,τi and k state.For the integral term, the trapezoidal formula is used to discrete it at grid pointzj,ti (see [9]) as follows:(23)∫−∞+∞Uk;zj+y,tifYydy=∑ι=0Mρι−jMUkι,i+Ukι+1,i+Rki,j,where(24)ρjM=12∫jΔzj+1ΔzfYydy=12∑ℓ=1m1p^ℓe−θ^ℓjΔz−e−θ^ℓj+1Δz,j≥0,12∑ℓ=1m2p˜ℓeθ˜ℓj+1Δz−eθ˜ℓjΔz,j≤0,(25)Rki,j=∫zM+1−zj+∞ezj+y−KerktifYydy=exj∑ℓ=1m1p^ℓθ^ℓθ^ℓ−1e1−θ^ℓzmax−zj−Kerkti∑ℓ=1m1p^ℓe−θ^ℓxmax−xj=ezmax−Kerkti∑ℓ=1m1p^ℓθ^ℓe−θ^ℓxmax−xj.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]):(26)Dzα−∞Uki,j=1Δzα∑ι=0∞gιUki,j−ι+1,where gι are the fractional difference coefficients given as follows:(27)gι=−1ιαι,αι=αα−1⋯α−ι+1ι!.Using the discretization and applying the fully implicit method, it leads to a nonlinear algebraic equation:(28)Uki+1,j−Uki,jΔt+akUki,j−Uki,j−1Δz+νΔzα∑ι=0∞gιUki,j+1−ι+ϵHUki,j+ϵ−qj+∑ℓ=0Mρℓ−jMUki,ℓ+Uki,ℓ+1+∑l=1,l≠kdqklUli,j+Rki,j−bkUki,j=0,with the boundary and terminal conditions as follows:(29)limj⟶−∞Uki,j=0,Uki,M+1=exmax−K,UkN+1,j=maxexj−K,0.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 1.
The coefficientsgι satisfy(30)g0=1,g1=−α,0≤⋯≤g3≤g2≤1,∑ι=0∞gι=0,for 1<α<2.Lemma 2.
Both the coefficientsρjM in equation (24) and Rki,j in equation (25) are bounded, and(31)∑j=−∞MρjM≤12,Rki,j≤ezmax−K.Proof.
From the fact thatfYy is the density function of hyperexponential random variables Y, it follows that(32)∑j=−∞MρjM=∑j=−∞M12∫jΔzj+1ΔzfYydy≤12∫−∞+∞fYydy=12.
In addition,p^i≥0,θ^i>1 and expzmax=3K or 4K, so it is straightforward to obtain(33)Rj,ki=ezmax−Kerkti∑ℓ=1m1p^ℓθ^ℓe−θ^ℓzmax−xj≤ezmax−Kerkti∑ℓ=1m1p^ℓ≤ezmax−K.Theorem 1.
If the time stepΔt≤1/qkk−bk−2∑ℓ=0Mρℓ−j0M and(34)C≥bkezmax+akezmaxezmax−1zmax+4ezmax−K+νezmax−1αzmaxα,then the approximate stock loan values Uki,j generated by the scheme (28) satisfy(35)Uki,j≥maxqj,0,for j=−∞,…,−2,−1,0,…,M+1, i=N+1,N,…,1, and k∈D.Proof.
First, we proveUki,j≥qj. Let(36)uki,j=Uki,j−qj,for all i,j, and k. It is straightforward to obtain ukN+1,j=UkN+1,j−qj≥0. Hence, by substituting uki,j into (28) and after simplification, it yields(37)1+Δtbk−akΔtΔzuki,j=uki+1,j+νΔtΔzα∑ι∞gιuki,j−ι+1+ϵΔtCuki,j+ϵ−qj−aΔtΔzuki,j−1+∑ℓ=0Mρℓ−jMuki,ℓ+uki,ℓ+1Δt+Δt∑l=1Mqkluki,j−ΔtEj,k,where(38)Ej,k=bkqj−akΔzqj−1−qj−νΔzα∑ι=0∞gιqj−ι+1−∑ℓ=0Mρℓ−jMqℓ+qℓ+1−∑l=1Mqklqj−Rki,j=bezj−akΔzezj−1−ezj−νΔzα∑ι=0∞gιqj−ι+1−K−∑ℓ=0Mρℓ−jMezℓ+ezℓ+1−2K−qkkqj−Rki,j.
FromeΔz−1/Δz≤ezmax−1/zmax and ∑ι=0∞gι=0, we have that(39)|Ej,k≤bkezj+νΔzα∑k=0∞gkej−k+1−K+aΔzezj−1−ezj+∑ℓ=0Mρℓ−jMezℓ+ezℓ+1−2K+Rj,ki+ezj−Kqkk≤bkezmax+akezmaxezmax−1zmax+Rj,ki+3ezmax−K+νΔzα∑k=0∞gkej−k+1.
From the fact∑ι=0∞gιezj−ι+1=ezj+1∑ι=0∞gιe−ιΔz and according to [14], 1−zα=∑ι=0∞gιz−ι when z≤1. Therefore, we obtain that(40)νΔzα∑ι=0∞gιe−ιΔt=ν1−e−ΔxαΔzα≤νezmax−1αzmaxα.
To sum up, combining Lemma2, it is straightforward to obtain(41)Ej,k≤bkezmax+akezmaxezmax−1zmax+4ezmax−K+νezmax−1αzmaxα.
Define(42)ui=minj,kuki,j,and let k0,j0 be a pair of indices, such that ui=uk0i,j0. So, from (28), it follows that(43)1+Δtbk−akΔtΔzui≥uk0i+1,j0+νΔtΔzα∑ι∞gιui+ϵΔtCui+ϵ−akΔtΔzui+∑ℓ=0Mρℓ−jMui+uiΔt+Δt∑l=1Mqklui−ΔtEj0,k0.
To reduce the form of inequality (43), we obtain that(44)1−qkk−bk−2∑ℓ=0Mρℓ−j0MΔtui−ϵΔtCui+ϵ+ΔtEj0,k0≥uk0i+1,j0≥ui+1.
On the other hand,(45)Δt≤1qkk−bk−2∑ℓ=0Mρℓ−j0M.
Hence, it is straightforward to obtain that(46)1−qkk−bk−2∑ℓ=0Mρℓ−j0MΔt≥0.
Let(47)A=1−γ−rk−ξ−2∑ℓ=0Mρℓ−jMΔt.
Define a function(48)fx=Ax−εΔtCx+ε+ΔtEj,k.
So, ifui+1≥0, we can obtain(49)fui≥0.
From(50)f0=ΔtEj−C≤0,it follows that(51)f′x=A+ϵΔtCx+ϵ2>0.
Based on the conditionukN+1,j≥0, we can obtain ui≥0. Naturally, uki,j≥0. Consequently,(52)Uki,j≥qj,for all i,j, and k.
Next, we prove thatUki,j≥0 for all i,j, and k. Following the idea above, define(53)Ui=minj,kUki,j,and let k0,j0 be a pair of indices, such that Uk0i,j0=Ui. It follows from (28) that(54)1+bkΔt−akΔtΔzUi≥Uk0i+1,j0+νΔtΔzα∑ι∞gιUi+ϵΔtCUi+ϵ−qj−akΔtΔzUi+∑ℓ=0Mρℓ−jMUi+UiΔt+Δt∑l=1MqklUi+Rki,j.
By simple calculation, the following equality holds:(55)1−qkk−bk−2∑ℓ=0Mρℓ−jMΔtUi−ϵΔtCUi+ϵ−qj0≥Uk0i+1,j0≥Ui+1.
From above, the inequalitiesUkj,i≥qj for all i,j,andk are verified. So,(56)1−qkk−bk−2∑ℓ=0Mρℓ−jMΔtUi≥Ui+1.
Notice that1−qkk−bk−2∑ℓ=0Mρℓ−jMΔt≥0, then(57)Ui≥Ui+1.
In addition,UkN+1,j=maxexpzj−K,0≥0 for all j,k. Namely, UN+1≥0. Therefore, it yields that(58)Uki,j≥0,for all i,j, and k, which completes the proof.


### 3.3. Implementation of Numerical Method

In practice, to implement the numerical scheme (28), the semi-infinite domain −∞,zmax×0,T must be truncated into a finite domain:(59)zmin,zmax×0,T.Here, we takezmin=ln0.01. Then, we can check that Uzmin,t=0. In such a case, we redefine the space step Δz=zmax−zmin/M, and(60)zj=zmin+j−1Δz,j=1,2,…,M+1.Let(61)ηk=akΔtΔz,βk=1−akΔtΔz+Δtbk,λk=−ΔtνkΔzα,WjM=ρjM+ρj−1M,j=0,±1,±2,…,±M−1,then the matrix form of numerical scheme (28) can be written as follows:(62)βkℐ+ℬk+Ak−ΔtWUki−FUki−Δt∑l=1,l≠kdqlkUki=Uki+1−Eki−ΔtRki,k=1,2,…,d,where Uki=Uki,2,Uki,3,…,Uki,Mτ, FUki=fUki,2,fUki,3,…,fUki,M, and(63)fUki,j=ΔtϵCUki,j+ϵ−qj.ℐ denotes the M−1×M−1 identity matrix.(64)W=W0MW1M⋯WM−2MWM−1MW−1MW0M⋯WM−3MWM−2MW−2MW−1M⋯WM−4MWM−3M⋮⋱⋱⋮⋮W2−MMW3−MM⋯W0MW1MW1−MMW2−MM⋯W−1MW0M,ℬk=ηk0⋯001⋯00⋮⋮⋮⋮0⋯000⋯10,Ak=λkg1g0⋯0g2g1⋯0g3g2⋯0⋮⋯⋱⋮gM−1gM−2⋯g1.Rki=Rki,2,Rki,3,…,Rki,Mτ, Ei=0,0,…,ϱM−1τ, and(65)ϱ=−η2g0−Δt∑ℓ=0M−1ρℓ−ΔtqkkUM+1i.In fact, system (62) is a coupling nonlinear equation so that it is very hard to obtain the numerical solution. Hence, we let(66)Ui=U1i,U2i,…,Udiτ,A=ℳ1Δtq12ℐ⋯Δq1dℐΔtq21ℐℳ2⋯Δtq2dℐ⋮⋮⋮⋮Δtqd1ℐΔtqd2ℐ⋯ℳd,where(67)ℳk=βkℐ+ℬk+Ak−ΔtW,k=1,2,…,d,ℱUi=FU1i,FU2i,…,FUdiτ,Ei=E1i,E2i,…,Ediτ,Ri=R1i,R2i,…,Rdiτ.Then, (62) can be written as the following nonlinear system:(68)AUi−ℱUi=Ui+1−Ei−ΔtRi.Note that system (68) cannot be solved directly since the penalty function ℱUi is nonlinear with respect to Ui. In the following, Newton iteration approach is employed to solve(69)A−Jℱωl−1δωl=Ui+1−Ei−ΔtRi+ℱωl−1−Aωl−1,ωl=ωl−1+κωl−ωl−1,where l=1,2,… and κ∈0,1 is a damping parameter. The initial value ω0=Ui+1 for each time level is the given initial guess and δωl=ωl−ωl−1. Jℱ is a Jacobian matrix with column vector ℱωl. We choose Ui=ωl. The 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.


## 3.1. Model Transformation

To avoid the effect of time-factoreγt in our model on the numerical results, we first introduce a new variable system as follows:(10)z=x−γt,Uk;z,t=eγtVk;x,t,zf,k=xf,k−γt.By variable transformation (the details are displayed in Appendix A), the pricing systems (6) and (8) can be written as(11)∂Uk;z,t∂t+ak∂Uk;z,t∂z+ξ∫−∞+∞Uk;z+y,tfYydy+νDzα−∞Uk;z,t+∑l=1,l≠kdqklUl;x,t=bkUk;z,t,(12)Uk;zf,t=expzf,k−K,(13)∂Uk;f,k,t∂z=expzf,k,(14)limz⟶−∞Uk;z,t=0,(15)Uk;x,T=maxez−K,0,where ak=rk−νk−D−γ−ξς and bk=rk+qkk+ξ−γ. And, inequality (9) should be rewritten as(16)Uk;z,t≥maxez−K,0,for all z≤zf,k and 0≤t≤T.It is straightforward to obtain that the functionUk;z,t in system (11)–(15) can be viewed as an American call option with strike price K and free-boundary ezf,k. The analytic solution of Uk;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],(17)εCUϵk;z,t+ϵ−qz,k∈D,qz=ez−K,to the corresponding equations in (11), we obtain the following system of nonlinear space-fractional partial integral-differential equations on a fixed domain:(18)∂Uϵk;z,t∂t+ak∂Uϵk;z,t∂z+ξ∫−∞+∞Uϵk;z+y,tfYydy+∑l=1,l≠kdqklUϵl;z,t+νkDzα−∞Uϵk;z,t+ϵCUϵk;z,t+ϵ−qz=bkUϵk;z,t,(19)Uϵk;zmax,t=ezmax−K,(20)limz⟶−∞Uϵk;z,t=0,(21)Uϵk;z,T=maxez−K,0,where ezmax denotes the maximum stock price, z∈−∞,zmax, 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 ezmax=3K in this paper. In addition, we will omit the subscript ϵ of function Uϵk;z,t for convenience.


## 3.2. Implicit Scheme

We first takeΔz>0 as spatial step, such that MΔz=zmax, where M is a positive integer and then place N+1 uniform grids in the t direction, namely, Δt=T/N. That is,(22)xj=j−1Δz,ti=i−1Δt,where j=∞,…,−1,0,1,…,M+1; i=1,2,…,N+1; and Uki,j denotes the value of function Uk;xj,τi at point xj,τi and k state.For the integral term, the trapezoidal formula is used to discrete it at grid pointzj,ti (see [9]) as follows:(23)∫−∞+∞Uk;zj+y,tifYydy=∑ι=0Mρι−jMUkι,i+Ukι+1,i+Rki,j,where(24)ρjM=12∫jΔzj+1ΔzfYydy=12∑ℓ=1m1p^ℓe−θ^ℓjΔz−e−θ^ℓj+1Δz,j≥0,12∑ℓ=1m2p˜ℓeθ˜ℓj+1Δz−eθ˜ℓjΔz,j≤0,(25)Rki,j=∫zM+1−zj+∞ezj+y−KerktifYydy=exj∑ℓ=1m1p^ℓθ^ℓθ^ℓ−1e1−θ^ℓzmax−zj−Kerkti∑ℓ=1m1p^ℓe−θ^ℓxmax−xj=ezmax−Kerkti∑ℓ=1m1p^ℓθ^ℓe−θ^ℓxmax−xj.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]):(26)Dzα−∞Uki,j=1Δzα∑ι=0∞gιUki,j−ι+1,where gι are the fractional difference coefficients given as follows:(27)gι=−1ιαι,αι=αα−1⋯α−ι+1ι!.Using the discretization and applying the fully implicit method, it leads to a nonlinear algebraic equation:(28)Uki+1,j−Uki,jΔt+akUki,j−Uki,j−1Δz+νΔzα∑ι=0∞gιUki,j+1−ι+ϵHUki,j+ϵ−qj+∑ℓ=0Mρℓ−jMUki,ℓ+Uki,ℓ+1+∑l=1,l≠kdqklUli,j+Rki,j−bkUki,j=0,with the boundary and terminal conditions as follows:(29)limj⟶−∞Uki,j=0,Uki,M+1=exmax−K,UkN+1,j=maxexj−K,0.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 1.
The coefficientsgι satisfy(30)g0=1,g1=−α,0≤⋯≤g3≤g2≤1,∑ι=0∞gι=0,for 1<α<2.Lemma 2.
Both the coefficientsρjM in equation (24) and Rki,j in equation (25) are bounded, and(31)∑j=−∞MρjM≤12,Rki,j≤ezmax−K.Proof.
From the fact thatfYy is the density function of hyperexponential random variables Y, it follows that(32)∑j=−∞MρjM=∑j=−∞M12∫jΔzj+1ΔzfYydy≤12∫−∞+∞fYydy=12.
In addition,p^i≥0,θ^i>1 and expzmax=3K or 4K, so it is straightforward to obtain(33)Rj,ki=ezmax−Kerkti∑ℓ=1m1p^ℓθ^ℓe−θ^ℓzmax−xj≤ezmax−Kerkti∑ℓ=1m1p^ℓ≤ezmax−K.Theorem 1.
If the time stepΔt≤1/qkk−bk−2∑ℓ=0Mρℓ−j0M and(34)C≥bkezmax+akezmaxezmax−1zmax+4ezmax−K+νezmax−1αzmaxα,then the approximate stock loan values Uki,j generated by the scheme (28) satisfy(35)Uki,j≥maxqj,0,for j=−∞,…,−2,−1,0,…,M+1, i=N+1,N,…,1, and k∈D.Proof.
First, we proveUki,j≥qj. Let(36)uki,j=Uki,j−qj,for all i,j, and k. It is straightforward to obtain ukN+1,j=UkN+1,j−qj≥0. Hence, by substituting uki,j into (28) and after simplification, it yields(37)1+Δtbk−akΔtΔzuki,j=uki+1,j+νΔtΔzα∑ι∞gιuki,j−ι+1+ϵΔtCuki,j+ϵ−qj−aΔtΔzuki,j−1+∑ℓ=0Mρℓ−jMuki,ℓ+uki,ℓ+1Δt+Δt∑l=1Mqkluki,j−ΔtEj,k,where(38)Ej,k=bkqj−akΔzqj−1−qj−νΔzα∑ι=0∞gιqj−ι+1−∑ℓ=0Mρℓ−jMqℓ+qℓ+1−∑l=1Mqklqj−Rki,j=bezj−akΔzezj−1−ezj−νΔzα∑ι=0∞gιqj−ι+1−K−∑ℓ=0Mρℓ−jMezℓ+ezℓ+1−2K−qkkqj−Rki,j.
FromeΔz−1/Δz≤ezmax−1/zmax and ∑ι=0∞gι=0, we have that(39)|Ej,k≤bkezj+νΔzα∑k=0∞gkej−k+1−K+aΔzezj−1−ezj+∑ℓ=0Mρℓ−jMezℓ+ezℓ+1−2K+Rj,ki+ezj−Kqkk≤bkezmax+akezmaxezmax−1zmax+Rj,ki+3ezmax−K+νΔzα∑k=0∞gkej−k+1.
From the fact∑ι=0∞gιezj−ι+1=ezj+1∑ι=0∞gιe−ιΔz and according to [14], 1−zα=∑ι=0∞gιz−ι when z≤1. Therefore, we obtain that(40)νΔzα∑ι=0∞gιe−ιΔt=ν1−e−ΔxαΔzα≤νezmax−1αzmaxα.
To sum up, combining Lemma2, it is straightforward to obtain(41)Ej,k≤bkezmax+akezmaxezmax−1zmax+4ezmax−K+νezmax−1αzmaxα.
Define(42)ui=minj,kuki,j,and let k0,j0 be a pair of indices, such that ui=uk0i,j0. So, from (28), it follows that(43)1+Δtbk−akΔtΔzui≥uk0i+1,j0+νΔtΔzα∑ι∞gιui+ϵΔtCui+ϵ−akΔtΔzui+∑ℓ=0Mρℓ−jMui+uiΔt+Δt∑l=1Mqklui−ΔtEj0,k0.
To reduce the form of inequality (43), we obtain that(44)1−qkk−bk−2∑ℓ=0Mρℓ−j0MΔtui−ϵΔtCui+ϵ+ΔtEj0,k0≥uk0i+1,j0≥ui+1.
On the other hand,(45)Δt≤1qkk−bk−2∑ℓ=0Mρℓ−j0M.
Hence, it is straightforward to obtain that(46)1−qkk−bk−2∑ℓ=0Mρℓ−j0MΔt≥0.
Let(47)A=1−γ−rk−ξ−2∑ℓ=0Mρℓ−jMΔt.
Define a function(48)fx=Ax−εΔtCx+ε+ΔtEj,k.
So, ifui+1≥0, we can obtain(49)fui≥0.
From(50)f0=ΔtEj−C≤0,it follows that(51)f′x=A+ϵΔtCx+ϵ2>0.
Based on the conditionukN+1,j≥0, we can obtain ui≥0. Naturally, uki,j≥0. Consequently,(52)Uki,j≥qj,for all i,j, and k.
Next, we prove thatUki,j≥0 for all i,j, and k. Following the idea above, define(53)Ui=minj,kUki,j,and let k0,j0 be a pair of indices, such that Uk0i,j0=Ui. It follows from (28) that(54)1+bkΔt−akΔtΔzUi≥Uk0i+1,j0+νΔtΔzα∑ι∞gιUi+ϵΔtCUi+ϵ−qj−akΔtΔzUi+∑ℓ=0Mρℓ−jMUi+UiΔt+Δt∑l=1MqklUi+Rki,j.
By simple calculation, the following equality holds:(55)1−qkk−bk−2∑ℓ=0Mρℓ−jMΔtUi−ϵΔtCUi+ϵ−qj0≥Uk0i+1,j0≥Ui+1.
From above, the inequalitiesUkj,i≥qj for all i,j,andk are verified. So,(56)1−qkk−bk−2∑ℓ=0Mρℓ−jMΔtUi≥Ui+1.
Notice that1−qkk−bk−2∑ℓ=0Mρℓ−jMΔt≥0, then(57)Ui≥Ui+1.
In addition,UkN+1,j=maxexpzj−K,0≥0 for all j,k. Namely, UN+1≥0. Therefore, it yields that(58)Uki,j≥0,for all i,j, and k, which completes the proof.


## 3.3. Implementation of Numerical Method

In practice, to implement the numerical scheme (28), the semi-infinite domain −∞,zmax×0,T must be truncated into a finite domain:(59)zmin,zmax×0,T.Here, we takezmin=ln0.01. Then, we can check that Uzmin,t=0. In such a case, we redefine the space step Δz=zmax−zmin/M, and(60)zj=zmin+j−1Δz,j=1,2,…,M+1.Let(61)ηk=akΔtΔz,βk=1−akΔtΔz+Δtbk,λk=−ΔtνkΔzα,WjM=ρjM+ρj−1M,j=0,±1,±2,…,±M−1,then the matrix form of numerical scheme (28) can be written as follows:(62)βkℐ+ℬk+Ak−ΔtWUki−FUki−Δt∑l=1,l≠kdqlkUki=Uki+1−Eki−ΔtRki,k=1,2,…,d,where Uki=Uki,2,Uki,3,…,Uki,Mτ, FUki=fUki,2,fUki,3,…,fUki,M, and(63)fUki,j=ΔtϵCUki,j+ϵ−qj.ℐ denotes the M−1×M−1 identity matrix.(64)W=W0MW1M⋯WM−2MWM−1MW−1MW0M⋯WM−3MWM−2MW−2MW−1M⋯WM−4MWM−3M⋮⋱⋱⋮⋮W2−MMW3−MM⋯W0MW1MW1−MMW2−MM⋯W−1MW0M,ℬk=ηk0⋯001⋯00⋮⋮⋮⋮0⋯000⋯10,Ak=λkg1g0⋯0g2g1⋯0g3g2⋯0⋮⋯⋱⋮gM−1gM−2⋯g1.Rki=Rki,2,Rki,3,…,Rki,Mτ, Ei=0,0,…,ϱM−1τ, and(65)ϱ=−η2g0−Δt∑ℓ=0M−1ρℓ−ΔtqkkUM+1i.In fact, system (62) is a coupling nonlinear equation so that it is very hard to obtain the numerical solution. Hence, we let(66)Ui=U1i,U2i,…,Udiτ,A=ℳ1Δtq12ℐ⋯Δq1dℐΔtq21ℐℳ2⋯Δtq2dℐ⋮⋮⋮⋮Δtqd1ℐΔtqd2ℐ⋯ℳd,where(67)ℳk=βkℐ+ℬk+Ak−ΔtW,k=1,2,…,d,ℱUi=FU1i,FU2i,…,FUdiτ,Ei=E1i,E2i,…,Ediτ,Ri=R1i,R2i,…,Rdiτ.Then, (62) can be written as the following nonlinear system:(68)AUi−ℱUi=Ui+1−Ei−ΔtRi.Note that system (68) cannot be solved directly since the penalty function ℱUi is nonlinear with respect to Ui. In the following, Newton iteration approach is employed to solve(69)A−Jℱωl−1δωl=Ui+1−Ei−ΔtRi+ℱωl−1−Aωl−1,ωl=ωl−1+κωl−ωl−1,where l=1,2,… and κ∈0,1 is a damping parameter. The initial value ω0=Ui+1 for each time level is the given initial guess and δωl=ωl−ωl−1. Jℱ is a Jacobian matrix with column vector ℱωl. We choose Ui=ωl. The 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.


## 4. 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.The model parameters are set to ber1=0.05, r2=0.04, σ1=0.20, σ2=0.24, D=0.8, γ=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 Uki,j≥maxqj,0 holds in Theorem 1, we show the mesh surface of Uki,j−maxqj,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 Uki,j−maxqj,0≥0 for all i,j, and k. In Figure 2, the values of Uk;x,t as a function of both stock price St and time t over the rectangular domain Smin,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.Figure 1
Mesh surface ofUki,j−maxqj,0 with M=26andN=100: (a) State 1 and (b) State 2.
(a)(b)Figure 2
Mesh surface of functionUk;x,t with M=26andN=100.
(a)(b)In fact, the functionUk;z,t in models (18)–(21) can be recognized as an American call option with strike K and optimal exercise boundary, expzf,k. Then, the values of Uk;z,t are not less than the payoff function maxez−K,0, and the function expzf,k is increasing with respect to the duration in every state. Figures 3 and 4 show these two facts, respectively.Figure 3
The values of functionUk;z,t with M=26andN=100.Figure 4
Optimal exercise priceexpzf,k under parameters with M=210andN=100.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 kandM 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−time1 and CPU−time2 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.Table 1
The comparison between our method and the method by using the iteration algorithm.
MCPU−time1CPU−time2k=2271.20101.2299283.96723.51002918.560115.690021087.956461.9984k=3278.23161.92462839.91094.298129205.120119.4887210513.890681.3718k=42715.78019.48712879.310215.711029429.364835.93142101012.3243121.6666


### 4.2. 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 Figure5, there are four different curves of optimal redemption prices under different parameters, and the horizontal axis denotes the time to expiration. The black and blue curves are obtained by the models without regime switching as Q=0. In the case of r1=0.05 or r1=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γT as shown in Figure 5.Figure 5
Optimal redemption pricesexpzf,k with σ1=σ2=0.24,D=0.8,γ=0.6,ξ=0.01,p^=0.02,θ^=1.2,T=2,andK=20.As shown in Figure6, the introduction of regime switching produces higher optimal redemption prices than those calculated by the framework without regime switching Q=0. This phenomenon is the same as an American option, of which the optimal exercise values are monotonically increasing functions of volatility. The 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γT as shown in Figure 6.Figure 6
Optimal redemption pricesexpzf,k with r1=r2=0.04,D=0.8,γ=0.6,ξ=0.01,p^=0.02,θ^=1.2,T=2,andK=20.


## 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.The model parameters are set to ber1=0.05, r2=0.04, σ1=0.20, σ2=0.24, D=0.8, γ=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 Uki,j≥maxqj,0 holds in Theorem 1, we show the mesh surface of Uki,j−maxqj,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 Uki,j−maxqj,0≥0 for all i,j, and k. In Figure 2, the values of Uk;x,t as a function of both stock price St and time t over the rectangular domain Smin,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.Figure 1
Mesh surface ofUki,j−maxqj,0 with M=26andN=100: (a) State 1 and (b) State 2.
(a)(b)Figure 2
Mesh surface of functionUk;x,t with M=26andN=100.
(a)(b)In fact, the functionUk;z,t in models (18)–(21) can be recognized as an American call option with strike K and optimal exercise boundary, expzf,k. Then, the values of Uk;z,t are not less than the payoff function maxez−K,0, and the function expzf,k is increasing with respect to the duration in every state. Figures 3 and 4 show these two facts, respectively.Figure 3
The values of functionUk;z,t with M=26andN=100.Figure 4
Optimal exercise priceexpzf,k under parameters with M=210andN=100.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 kandM 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−time1 and CPU−time2 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.Table 1
The comparison between our method and the method by using the iteration algorithm.
MCPU−time1CPU−time2k=2271.20101.2299283.96723.51002918.560115.690021087.956461.9984k=3278.23161.92462839.91094.298129205.120119.4887210513.890681.3718k=42715.78019.48712879.310215.711029429.364835.93142101012.3243121.6666


## 4.2. 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 Figure5, there are four different curves of optimal redemption prices under different parameters, and the horizontal axis denotes the time to expiration. The black and blue curves are obtained by the models without regime switching as Q=0. In the case of r1=0.05 or r1=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γT as shown in Figure 5.Figure 5
Optimal redemption pricesexpzf,k with σ1=σ2=0.24,D=0.8,γ=0.6,ξ=0.01,p^=0.02,θ^=1.2,T=2,andK=20.As shown in Figure6, the introduction of regime switching produces higher optimal redemption prices than those calculated by the framework without regime switching Q=0. This phenomenon is the same as an American option, of which the optimal exercise values are monotonically increasing functions of volatility. The 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γT as shown in Figure 6.Figure 6
Optimal redemption pricesexpzf,k with r1=r2=0.04,D=0.8,γ=0.6,ξ=0.01,p^=0.02,θ^=1.2,T=2,andK=20.


## 5. Conclusions

In this paper, we investigate the pricing model of stock loan under the FMLS model with jump diffusion process in ad-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 Vk,x,t≥maxex−Keγt,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 Table1.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.



---

*Source: 8818876-2021-07-09.xml*