# Laplace Transform Methods for a Free Boundary Problem of Time-Fractional Partial Differential Equation System

Authors: Zhiqiang Zhou; Xuemei Gao
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

We study the pricing of the American options with fractal transmission system under two-state regime switching models. This pricing problem can be formulated as a free boundary problem of time-fractional partial differential equation (FPDE) system. Firstly, applying Laplace transform to the governing FPDEs with respect to the time variable results in second-order ordinary differential equations (ODEs) with two free boundaries. Then, the solutions of ODEs are expressed in an explicit form. Consequently the early exercise boundaries and the values for the American option are recovered using the Gaver-Stehfest formula. Numerical comparisons of the methods with the finite difference methods are carried out to verify the efficiency of the methods.


---

## Body

## 1. Introduction

Markov regime switching models were first introduced by Hamilton [1] and recently have become popular in financial applications including equity options [2–17], bond prices and interest rate derivatives [18–20], portfolio selection [21], and trading rules [22–26]. The Markov regime switching models allow the model parameters (drift and volatility coefficients) to depend on a Markov chain which can reflect the information of the market environments and at the same time preserve the simplicity of the models. However when the model parameters are governed by the Markov chain, the valuation of the options becomes complex.American option pricing is a kind of classical free boundary problems, where the free boundary refers to early exercise boundary or optimal exercise boundary. Laplace transform methods have been developed to solve the free boundary problems arising in American option pricing under geometric Brownian motion (GBM) (see Mallier and Alobaidi [27], Zhu [28], and Zhu and Zhang [29]) and constant elasticity of variance (CEV) (see Wong and Zhao [30]). The essential idea of the Laplace transform methods for solving the American option pricing problems is described as follows. Applying Laplace transform to the governing free boundary partial differential equations (PDEs) with respect to the time variable results in a boundary value problem of second-order ordinary differential equations (ODEs). Unlike the case for the European option pricing, the ODE involves an unknown boundary (in Laplace space) which needs to be solved. If the solution of the ODE with the unknown boundary can be expressed in an explicit form, then a nonlinear algebraic equation for the unknown boundary can be derived. Using a simple solver for the nonlinear algebraic equation, we can get the value of the unknown boundary in Laplace space. Consequently the early exercise boundary and the value for the American option can be achieved using the inverse Laplace transform.Chen et al. [31] studied a predictor-corrector finite difference methods for pricing American options under the FMLS model which is a kind of space-fractional derivative model. Liang et al. [32] used Fourier transform to solve a bifractional Black-Scholes-Merton differential equation. In this paper, we study the pricing of American options with time-fractional model which has essential difference to the space-fractional model. Following the model in [33], we assume that the underlying asset price still follows the classical Brownian motion, but the change in the option price is considered as a fractal transmission system. The price of such American option follows time-fractional partial differential equations (PDEs) with free boundaries. The solution of the time-fractional PDEs with two free boundaries is more challenging than solving the fractional PDEs with fixed boundary for European option pricing in [33] and much more complex than solving single Black-Scholes-Merton differential equation in [32].In this paper, we develop the Laplace transform methods for pricing time-fractional American options under regime switching models. The fundamental difference to the GBM and CEV American option pricing is that the governing free boundary PDEs for the regime switching American option pricing are a coupled system. Therefore the Laplace transform in time for the system of free boundary PDEs leads to a coupled system of ODEs whose analytic solutions take much trouble to the approach. A careful derivation in different cases makes the approach successful. Numerical examples are provided to verify the efficiency of the approach. The generated option value and early exercise boundary are compared with the finite difference method which are usually used as the benchmark methods in the area of option pricing.


## 2. Laplace Transform Methods for American Option Pricing

Let the underlying asset pricesSt follow a two-state regime switching model under risk-neutral measure:(1)dStSt=rYtdt+σYtdWt,where Wt is a standard Brownian motion and Y(t) is a continuous-time Markov chain with two states (y1,y2). Assume that, at each state Y(t)=yk,k=1,2, the interest rates ryk=rk and volatility σyk=σk for k=1,2 are nonnegative constants. Let(2)A=-q1q1q2-q2be the generator matrix of the Markov chain process with qi≥0(i=1,2). Let T be the option expiration date, τ=T-t be the time to maturity, Vi(S,τ) represent the American put option prices with stock price S, and S¯i(τ) be the corresponding optimal exercise boundaries at time τ. Then the American put options satisfy the following system of PDEs with free boundaries, for i=1,2:(3)∂Vi∂τ=LiVi+qiVj,j≠i,S>S¯iτ,τ>0,where Li are the differential operators defined by(4)Li=12σi2S2∂2∂S2+riS∂∂S-ri+qi.However, it is argued that the time derivative ∂Vi/∂τ should be replaced by the fractional derivative ∂αVi/∂τα(0<α≤1) under the assumption that the change in the option price follows a fractal transmission system (see, e.g., [32–34]). So, the time-fractional American put options satisfy the following system of PDEs:(5)∂Vi∂τα=LiVi+qiVj,j≠i,S>S¯iτ,τ>0,with the free boundary conditions and the initial condition:(6)ViS,τ=K-S,0≤S≤S¯iτ,(7)ViS,0=max⁡K-S,0,(8)limS↑∞⁡ViS,τ=0,(9)limS↓S¯iτ⁡ViS,τ=K-S¯iτ,(10)limS↓S¯iτ⁡∂∂SViS,τ=-1,where ∂Vi/∂τα is defined as the Caputo fractional derivative for i=1,2:(11)∂Vi∂τα≡1Γ1-α∫0τ1τ-θα∂ViS,θ∂θdθ,0<α<1,∂Vi∂τ,α=1.Condition (10) is called “high-contact condition” or “smooth-pasting condition.”Forλ>0, we define the Laplace-Carson transform (LCT) of the American put option price Vi(x,τ) as(12)Vi^S,λ=∫0+∞ViS,τλe-λτdτ≕LCViS,τλ,which is essentially the Laplace transform (LT) as(13)LCViS,τλ=λLViS,τλ.Using the Laplace transform formula for the Caputo fractional derivative (see (2.253) in [35]),(14)L∂ViS,τ∂ταλ=λαLViS,τλ-λα-1ViS,0,and the relationship (13), the LCT for ∂Vi/∂τα is found as(15)LC∂ViS,τ∂ταλ=λαLCViS,τλ-λαViS,0.The use of LCT simplifies notations in the later analysis and makes inverse Laplace solutions have a unique singularity λ=0. We note that for the case α=1 formula (15) is still true. Denoting S^i(λ)=LC[S¯i(τ)](λ) and taking LCTs to (3)–(10), we obtain the following ODEs, for i=1,2:(16)12σi2S2∂2V^i∂S2+riS∂V^i∂S-ri+qi+λαV^i+qiV^j+λαmax⁡K-S,0=0,S>S^iλ,(17)V^iS,λ=K-S,0≤S≤S^iλ,(18)limS↑∞⁡V^iS,λ=0,(19)limS↓S^iλ⁡V^iS,λ=K-S^iλ,(20)limS↓S^iλ⁡∂∂SV^iS,λ=-1.The existence of Laplace transform S^i(λ) and Vi^(S,λ) is guaranteed by the continuity and boundedness of S¯i(τ) and Vi(S,τ). To solve (16)–(20), we use the following variable and function transformations:(21)x=ln⁡SK,x^iλ=ln⁡S^iλK,V^iS,λ=KU^ix,λ.Let ai=(1/2)σi2,bi=ri-1/2σi2,ci=ri+qi. Then the calculations yield an equivalent form to (16)–(20) as follows:(22)ai∂2U^i∂x2+bi∂U^i∂x-ci+λαU^i+qiU^j+λαmax⁡1-ex,0=0,x>x^iλ,(23)U^ix,λ=1-ex,-∞≤x≤x^iλ,(24)limx↑∞⁡U^ix,λ=0,(25)limx↓x^iλ⁡U^ix,λ=1-ex^iλ,(26)limx↓x^iλ⁡∂∂xU^ix,λ=-ex^iλ.Compared to (16), the coefficients of the first and second derivatives in (22) become constant which will simplify the calculation. However, (22) includes two unknown variables x^i(i=1,2) and a nonsmooth term max⁡(1-ex,0), which will make the solutions somewhat complicated.From the definition ofS^i(λ) and monotone decreasing property of S¯i(τ) (see, e.g., [10]), we have(27)0<S^iλ=LCS¯τλ<S¯i0=K,and therefore(28)-∞<x^iλ=ln⁡S^iλK<0.Moreover U^i(x,λ) and its first-order derivative are continuous at x=0 (x=0 means that S^i=K). Without loss of generality, we assume x^1(λ)<x^2(λ). So -∞<x^1(λ)<x^2(λ)<0. The solutions of (22)–(26) can be solved separately on four intervals I1=(-∞,x^1),I2=[x^1,x^2),I3=[x^2,0), and I4=[0,+∞).Proposition 1.
The solutions of (22)–(26) on interval I1=(-∞,x^1) are given by(29)U^1x,λ=1-ex,U^2x,λ=1-ex,x∈I1.Proof.
The proof is straightforward by noticing (23).The solutions of ODEs on the other intervals can be written as a general solution to the corresponding homogeneous part plus a particular one satisfying the nonhomogeneous system. More precisely, we discuss the solutions on the different intervals as follows.Proposition 2.
On intervalI2=[x^1,x^2), the solutions of (22)–(26) are given by(30)U^1x,λ=ΦxC1+μλx,U^2x,λ=1-ex,x∈I2,and constant C1 is uniquely determined by the boundary conditions (25) and (26); that is,(31)C1=Φ~-1x^1f1x^1.Here the unspecified notations are defined in (34)–(37) below.Proof.
Forx∈I2, the solution U^2(x,λ)=1-ex and then U^1(x,λ) satisfies the following equation:(32)a1∂2U^1∂x2+b1∂U^1∂x-c1+λαU^1+q1+λα1-ex=0.By solving ODE (32), we get the general solution(33)U^1x,λ=C11ϕ1x+C12ϕ2x+μλx,where μλ(x) is a particular solution and ϕ1(x) and ϕ2(x) are the fundamental increasing and decreasing solutions with explicit forms:(34)μλx=q1+λαc1+λα+q1+λαa1+b1-c1-λαex,(35)ϕ1x=eγ1+x,ϕ2x=eγ1-x,γ1±=-b1±b12+4a1c1+λα2a1.Note that γ1+>0 and γ1-<0. Denoting(36)C1=C11C12,Φx=ϕ1xϕ2xT,μ~λx=μλxddxμλx,(37)Φ~x=ΦxddxΦx,f1x=1-ex-ex-μ~λx,and using (23), we complete the proof.Now we derive the solutions of (22)–(26) on interval I3=[x^2,0). Let y(x)=[y1(x),y2(x),y3(x),y4(x)]T with(38)y1x=U^1x,λ,y2x=∂U^1x,λ∂x,y3x=U^2x,λ,y4x=∂U^2x,λ∂x.Then (22) is equivalent to the following first-order ODE system:(39)y′x=Ayx+gx,where(40)A=0100c1+λαa1-b1a1-q1a100001-q2a20c2+λαa2-b2a2,gx=0λαex-10λαex-1.Lemma 3.
For any real numberλ>0, the homogeneous equation y′(x)=Ayx has real valued fundamental solutions(41)Ψ~x=eν1xξ1,eν2xξ2,eν3xξ3,eν4xξ4,where ν1 and ν2 are different negative real numbers, ν3 and ν4 are different positive real numbers, and ξj(j=1,2,3,4) are real valued column vectors.Proof.
It is well-known thaty′(x)=Ay(x) has fundamental solution matrix as (41), with νj and ξj representing characteristic values and corresponding characteristic vectors. So we only need to prove all νj are real numbers and different from each other. The characteristic polynomial of (39) is(42)Pν=det⁡νI-A=ν-100-c1+λαa1ν+b1a1q1a1000ν-1q2a20-c2+λαa2ν+b2a2=ν2+b1a1ν-c1+λαa1ν2+b2a2ν-c2+λαa2-q1q2a1a2.Therefore,(43)Pν=1a1a2a1ν2+b1ν-c1+λαa2ν2+b2ν-c2+λα-q1q2=1a1a2ν-γ1-ν-γ1+ν-γ2-ν-γ2+-q1q2,where(44)γi±=-bi±bi2+4aici+λα2ai.Since ai=σi/2>0 and ci=ri+qi>0, for any λ≥0 we know that all γi± are real numbers and γi-<0(i=1,2) and γi+>0(i=1,2). Moreover, we have(45)Pγi±=-q1q2a1a2<0,i=1,2,P0=c1+λαc2+λα-q1q2a1a2>0,limν→±∞⁡Pν>0.Then there must exist four real roots νj(j=1,2,3,4) solving characteristic equation P(ν)=0, and we can list the ranges of these values:(46)-∞<ν1<mini⁡γi-,maxi⁡γi-<ν2<0,0<ν3<mini⁡γi+,maxi⁡γi+<ν4<+∞.Note that the roots νj(j=1,2,3,4) can be searched by the secant method on their corresponding ranges as above.Proposition 4.
LetΨ(x) be a submatrix extracting the first and the third rows from Ψ~(x). Then the solutions U^(x,λ)=(U^1,U^2)T of (22)–(26) on interval I3=[x^2(λ),0) are given by(47)U^x,λ=ΨxC2+ηλx,where Ψ~(x) is given by (41) and C2 by (52) below.Proof.
It can be verified that (22) has a particular solution(48)ηλx=η1,λxη2,λx=-c1+λαq1q2-c2+λα-1-λα-λα+δ1q1q2δ2-1λαexλαex.Here, δi=ai+bi-(ci+λ)=-qi-λ for i=1,2. Denote C2=[C21,C22,C23,C24]T. Then the solution to (39) can be written as (49)yx=Ψ~xC2+η~λx,η~λx=η1,λx,ddxη1,λx,η2,λx,ddxη2,λxT.Using the continuity and smoothing-pasting of solutions at x=x^2, we know that C2 satisfies(50)Ψ~x^2C2+η~λx^2=Φ~x^2C1+μ~λx^21-ex^2-ex^2.Substituting C1=Φ~-1(x^1)f1(x^1) into the above equation and denoting(51)f2x^1,x^2=Φ~x^2Φ~-1x^1f1x^1+μ~λx^21-ex^2-ex^2-η~λx^2,we obtain that(52)C2=Ψ~-1x^2f2x^1,x^2.Thus the proof is complete.Proposition 5.
LetΨ^(x) be the submatrix extracting the first and the second columns from Ψ(x), where Ψ(x) is given in Proposition 4. Then the solutions U^(x,λ)=(U^1,U^2)T of (22)–(26) on interval I4=[0,+∞) are given by(53)U^x=Ψ^xC3.with C3 given by (56).Proof.
Sincemax⁡(1-ex,0)=0 for x∈I4=[0,+∞), (22) becomes homogeneous. Furthermore using condition limx↑∞U^i(x,λ)=0 (see (24)), we can write the solution to the homogeneous part of (24) (or homogeneous part of (39)) as(54)U^x=Ψ^xC3,where Ψ^(x) is the submatrix containing the first and the second columns of Ψ(x) and C3=[C31,C32]T. Notations Ψ~,Ψ, and Ψ^ are 4×4, 2×4, and 2×2 matrices, respectively. Since the first-order derivative of U^i(x,λ) is also continuous at x=0, by evaluating the first-order derivatives of (53) and (47) at x=0, we obtain that(55)Ψ′0C2+ηλ′0=Ψ^′0C3,where Ψ′(0)=(dΨ(x)/dx)|x=0 and Ψ^′(0)=(dΨ^(x)/dx)|x=0. Therefore, we obtain(56)C3=Ψ^′0-1Ψ′0C2+ηλ′0.SinceU^i(x,λ) is continuous at x=0, by evaluating (53) and (47) at x=0, we have(57)Ψ0C2+ηλ0=Ψ^0C3.Substituting (52) and (56) into (57) and simplifying the calculation, we obtain a system of nonlinear algebra equations with respect to x^1 and x^2(58)PλΨ~-1x^2f2x^1,x^2+Qλ=0,where constant matrices Pλ and Qλ are defined by(59)Pλ=Ψ^0Ψ^′0-1Ψ′0-Ψ0,Qλ=Ψ^0Ψ^′0-1ηλ′0-ηλ0.We summarize the above discussions in the following theorem.Theorem 6.
For any fixedλ>0, the optimal exercise boundaries x^1(λ) and x^2(λ) can be numerically computed from nonlinear equations (58) and the solutions U^(x,λ)=[U^1(x,λ),U^2(x,λ)]T of (22)–(26) are given by (60)U^1x,λ=1-ex,x∈I1,U^1x,λ=ΦxC1+μλx,x∈I2,U^2x,λ=1-ex,x∈I1∪I2,U^x,λ=ΨxC2+ηλx,x∈I3,U^x,λ=Ψ^xC3,x∈I4,where C1,C2, and C3 are given by formulas (31), (52), and (56), respectively; the fundamental solution matrices Φ(x), Ψ(x), and Ψ^(x) are defined in (37) and Propositions 4 and 5, respectively; the particular solutions μλ(x) and ηλ(x) are given by (34) and (48), respectively.Finally, the original optimal exercise boundariesS¯i(τ) and American put option prices Vi(S,τ) can be expressed in terms of the Laplace inversion:(61)S¯iτ=KeL-1x^iλ/λ,i=1,2,ViS,τ=KL-11λU^iln⁡SK,λ,i=1,2.The numerical approximation for Laplace inversionf(t)=L-1f^(λ) can be computed using the most powerful Gaver-Stehfest algorithm (see, e.g., [36])(62)ft=ln⁡2t∑k=1nCknf^kln⁡2t,with(63)Ckn=-1k+n/2∑j=k+1/2min⁡k,n/2jn/22j!n/2-j!j!j-1!k-j!2j-k!;here n must be taken as even number.


## 3. Numerical Examples

In this section, the Laplace transform method (LTM) is compared with the finite difference method (FDM). We use symbolic manipulation and multiprecision computing provided by Mathematica 9.0.Table1 lists the computational values of American put options under regime switching. Columns entitled “LTM" and “FDM” represent the option values obtained by Laplace transform method and finite different method, respectively. For the LTM the Gaver-Stehfest formula is applied with the number n=6 and for the FDM the number of time nodes M=200 and the number of space nodes N=400. The relative error (RE) is the absolute value of the difference between the computed LTM option values and the FDM values divided by the FDM option values. In the implementation of the LTM, the work precision is set by 30 significant digits, while the numerical values are presented with 5 significant digits.Table 1
Prices of American put option with parametersσ1=0.8,σ2=0.3,r1=0.1,r2=0.05,q1=6,q2=9,K=9,T=1 and different values of α.
α
S
0
Regime 1
Regime 2
LTM
FDM
RE (%)
LTM
FDM
RE (%)
α
=
1
3.00
6.0000
6.0000
0.0000
6.0000
6.0000
0.0000
4.50
4.4924
4.5431
1.1160
4.5278
4.5117
0.3568
6.00
3.3692
3.4138
1.3065
3.3071
3.3502
1.2865
7.50
2.5536
2.5835
1.1573
2.4741
2.5026
1.1388
9.00
1.9512
1.9712
1.0146
1.8630
1.8816
0.9885
10.50
1.5047
1.5177
0.8566
1.4146
1.4265
0.8342
12.00
1.1717
1.1795
0.6613
1.0847
1.0915
0.6230
α=0.7
3.00
6.0000
6.0000
0.0000
6.0000
6.0000
0.0000
4.50
4.5084
4.5436
0.7747
4.5181
4.5149
0.0709
6.00
3.3787
3.3988
0.5914
3.3229
3.3420
0.5715
7.50
2.5258
2.5381
0.4846
2.4442
2.4537
0.3872
9.00
1.9016
1.9089
0.3824
1.7972
1.8037
0.3604
10.50
1.4633
1.4672
0.2658
1.3659
1.3691
0.2337
12.00
1.1517
1.1531
0.1214
1.0668
1.0678
0.0937
α=0.4
3.00
6.0000
6.0000
0.0000
6.0000
6.0000
0.0000
4.50
4.5275
4.5370
0.2094
4.5067
4.5112
0.0998
6.00
3.3630
3.3677
0.1396
3.3102
3.3146
0.1327
7.50
2.4738
2.4749
0.0444
2.3871
2.3877
0.0251
9.00
1.8271
1.8257
0.0767
1.7102
1.7085
0.0995
10.50
1.3915
1.3883
0.2305
1.2882
1.2849
0.2568
12.00
1.0940
1.0891
0.4499
1.0094
1.0046
0.4778
CPU time
1.98 s for LTM; 532.47 s for FDMThe relative errors shown in Table1 are less than 1.5% for all the cases, which illustrates that the LTM is competitively accurate for solving American option pricing problems. Since there are no exact solutions of FPDEs, we cannot judge which one of the two results is more accurate. In general, the relative error less than five percent is acceptable in financial engineering practice.From the last row of Table1, we see that the CPU time of LTM is much less than that of FDM. The computational cost of LTM is uniform at every time t=T-τ and the CPU consumption mainly comes from solving n nonlinear algebraic equations with two variables x^1 and x^2 (see expression (58)). FDM is a time-stepping scheme which is needed to solve M nonlinear systems at time τ=T or t=0. Moreover, for the global dependence of the fractional differential operator, FDM needs to do some additional computation. We observe that the numbers n=6 of the quadrature nodes are much smaller than the numbers M=200 of time steps. In other words, to achieve the similar accuracy, there are much fewer nonlinear systems needed to be solved by the LTM than that by the FDM.Figure1 plots the early exercise boundaries obtained by LTM and FDM with the same values of parameters as in Table 1. It can be seen that the boundaries computed by FDM and LTM are pretty close. Again it shows that the LTM is robust in generating the early exercise boundaries.Figure 1
Early exercise boundaries computed by LTM and FDM withα=0.7 and the same parameters as in Table 1. (a) is for regime 1 and (b) for regime 2.
(a)
(b)


## 4. Conclusions

In this paper, we have developed Laplace transform methods to solve the time-fractional American option pricing under regime switching models. The value of the American option with regime switching is formulated as the solution to a free boundary problem of time-fractional partial differential equation system. The Laplace transform is executed for the time variables and the resulting system PDEs are solved analytically. Consequently a system of nonlinear algebraic equations for the free boundaries is obtained and solved using secant methods. Finally numerical Laplace inversion is applied to recover the early exercise boundaries and the option values. Comparisons between the LTM and the benchmark FDM are made via numerical examples, which shows that the LTM is efficient for pricing American options with regime switching. However, the LTM is still challenging for more complex models like the regime switching models with state-dependent jump diffusions. This will be left for the future studies.



---

*Source: 6917828-2017-07-10.xml*