# Numerical Contour Integral Methods for Free-Boundary Partial Differential Equations Arising in American Volatility Options Pricing

Authors: Yong Chen; Jianjun Ma
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

The aim of this paper is to study the numerical contour integral methods (NCIMs) for solving free-boundary partial differential equations (PDEs) from American volatility options pricing. Firstly, the governing free-boundary PDEs are modified as a unified form of PDEs on the fixed space region; then performing Laplace-Carson transform (LCT) leads to ordinary differential equations (ODEs) which involve the unknown inverse functions of free boundaries. Secondly, the inverse free-boundary functions are approximated and optimized by solving of the free-boundary values of the perpetual American volatility options. Finally, the ODEs are solved by the finite difference methods (FDMs), and the results are restored via the numerical Laplace inversion. Numerical results confirm that the NCIMs outperform the FDMs for solving free-boundary PDEs in regard to the accuracy and computational time.


---

## Body

## 1. Introduction

The pricing problems of European-style options written on volatility are widely studied in the literature by deriving the explicit formulas [1–6], by using the simulation approaches [7, 8], by using the PDE approaches [5, 9–13], and by the empirical analysis [14]. However, the valuation of American-style volatility options is more challenging due to the early exercise (free boundary) property. Detemple and Osakwe [15] derive the early exercise premium (EEP) representations of American-style volatility options under the general stochastic volatility models: geometric Brownian motion process (GBMP), mean-reverting Gaussian process (MRGP), mean-reverting square-root process (MRSRP), and the mean reverting log process (MRLP). The EEP representations result in complex integral equations; however the integral equations cannot be analytically solved and the iterative methods for finding the numerical solutions of the integral equations somehow lack convincing accuracy and reliability.This paper studies the numerical contour integral methods (NCIMs) for solving free-boundary partial differential equations (PDEs) from American volatility options written on the volatility whose price follows four well-known models: GBMP, MRGP, MRSRP, and MRLP. The original idea of NCIMs is developed by Zhou, Ma, and Sun [16] for solving free-boundary problems of space-fractional diffusion equations; then Ma and Zhu [8] prove the convergence rates of such methods under the regime-switching European option pricing, and it can be described as follows. By approximating the inverse function of free boundary, the free-boundary PDE could be written as a unified form of PDE defined on a fixed space region; then applying Laplace transform to the unified PDE with respect to time variable results in a boundary value problem of ODE; finally the finite difference method (FDM) is adopted for solving the aimed ODE, and the computational results are restored via the numerical Laplace inversion. However, it should be pointed out that the original paper needs another consumable iterative algorithm for approximating and optimizing the inverse function of free boundary, and this paper has no such procedure by means of solving of the free-boundary values of the perpetual American volatility options, which can optimize the approximated inverse free-boundary functions.The remaining parts of this paper are arranged as follows: In Section2, the free-boundary PDEs governing the American volatility option price are presented and some accompanied properties are introduced. In Section 3, the NCIMs are studied for solving the aimed free-boundary PDEs and some essential theorems are studied. In Section 4, numerical examples are carried out to verify the effectiveness of NCIMs. Conclusions are given in the final section.


## 2. Model Descriptions

LetVt denote the volatility of the underlying asset. Then we consider the following four stochastic volatility models [15] listed in the second column of Table 1. Meanwhile the infinitesimal generators are listed in the third column of Table 1.Table 1
Volatility models and infinitesimal generators.
Model Volatility process A P ( V , τ ) GBMP d V t = V t ( 2 μ + σ 2 ) d t + 2 σ d Z ~ t 2 σ 2 V 2 ∂ 2 P ∂ V 2 + ( 2 μ + σ 2 ) V ∂ P ∂ V - r P MRGP d V t = ( α - λ V t ) d t + σ d Z ~ t σ 2 2 ∂ 2 P ∂ V 2 + ( α - λ V ) ∂ P ∂ V - r P MRSRP d V t = ( σ 2 - 2 λ V t ) d t + 2 σ V t d Z ~ t 2 σ 2 V ∂ 2 P ∂ V 2 + ( σ 2 - 2 λ V ) ∂ P ∂ V - r P MRLP d ln ⁡ ( V t ) = α - λ ln ⁡ ( V t ) d t + σ d Z ~ t σ 2 V 2 2 ∂ 2 P ∂ V 2 + α - λ ln ⁡ V V ∂ P ∂ V - r PIn Table1, Z~t denotes the Brownian motion process under risk neutral measure P, μ, α, λ, r (risk-free interest rate) and σ (volatility of volatility) are constants.In the following, we describe the free-boundary partial differential equations from American volatility options. Due to the limitation of length, the put options are considered in this paper and the call options are studied similarly without essential difficulties. Denote by(1)fV,t=supτ∈T0,T⁡Ee-rτ-tmax⁡K-Vτ,0∣Vt=V,the price of American-style put option with underlying stochastic volatility Vt, current time t, maturity T, and strike price K, where T0,T is the optimal stopping set with the Brownian motion filtration for t∈[0,T].Using transformation of time variableτ=T-t (namely, the time to maturity), then the option price and the early exercise boundary (namely, the free-boundary function) are denoted by P(V,τ)≔f(V,T-τ)=f(V,t) and Vf(τ), respectively. Then the valuation of American volatility put option can be formulated as a free-boundary problem(2)∂PV,τ∂τ=APV,τ,V>Vfτ,τ>0,(3)PV,τ=K-V,0≤V≤Vfτ(4)∂PVfτ,τ∂V=-1,(5)limV→∞⁡PV,τ=0,(6)PV,0=max⁡K-V,0,where we can verify that(7)APV,τ1V≤Vfτ=AK-V,where 1{·} represents the indicator function and the values of A(K-V) for different models are listed in the second column of Table 2. Moreover the early exercise boundaries at initial Vf(0) are given in the third column of Table 2, where V∗ is the unique positive root of the following nonlinear algebraic equation:(8)r+λln⁡V-αV-rK=0,Table 2
The values ofA(K-V) and Vf(0).
Model A ( K - V ) V f ( 0 ) GBMP ( r - 2 μ - σ 2 ) V - r K K MRGP ( r + λ ) V - α - r K min ⁡ K , α + r K λ + r MRSRP ( r + 2 λ ) V - σ 2 - r K min ⁡ K , σ 2 + r K 2 λ + r MRLP ( r + λ ln ⁡ V - α ) V - r K min ⁡ K , V ∗It is proved by [17, 18] for GBMP model and [18] for MRGP model, MRSRP model, and MRLP model that the free-boundary functions Vf(τ) are continuously differentiable, strictly decreasing and convex on the interval V∈[0,+∞), which are confirmed by the simulation results in [19]. Moreover we note that the free-boundary functions Vf(τ) are bounded on the interval (V_,V¯], where V_=Vf(∞) are the early exercise boundaries of the corresponding perpetual American volatility put options, whose values satisfy the following ODEs:(9)AP∞V=0,V>Vf∞,(10)P∞V=K-V,0≤V≤Vf∞(11)dP∞Vf∞dV=-1,(12)limV→∞⁡P∞V=0,where the differential operators A are listed in the third column of Table 1 with replacement of P by P∞.


## 3. Numerical Contour Integral Methods for Free-Boundary PDEs from American Volatility Options Pricing

The aim of this section is to study the numerical contour integral methods (NCIMs) for solving the free-boundary PDEs (2)-(6).From the PDEs (2)-(6), if the volatility falls into the exercise region V∈[0,Vf(τ)] for τ∈[0,∞), then the value of American volatility put is given by P(V,τ)=K-V, and we can verify that it is equivalent to the following PDE:(13)∂PV,τ∂τ=APV,τ-AK-V,0≤V≤Vfτ,τ>0.Combining (2) and (13) gives a unified form of PDE:(14)∂PV,τ∂τ=APV,τ-AK-V1V≤Vfτ,V≥0,τ>0.Recall thatV_∈0,Vf(τ) for all τ∈[0,+∞); the value of American volatility put option equals K-V on the region V∈[0,V_]. Therefore it suffices to solve the PDE (14) with conditions(15)PV_,τ=K-V_,(16)∂PV_,τ∂V=-1,(17)limV→∞PV,τ=0,(18)PV,0=maxK-V,0.Now we focus on the Laplace-Carson transform (LCT) methods for solving the unified PDE (14). For any complex value z, the LCT is defined as(19)P^V,z=∫0+∞PV,τze-zτdτ≔LCPV,τz,which is essentially the same as the Laplace transform (LT), and the reason for using the LCT is to simplify the notations in later analysis. The relationship between LCT and LT is given by(20)LCPV,τz=zLPV,τz.Letτf(V) be the inverse function of Vf(τ). Since τf(V) is a strictly decreasing function on the interval (V_,V¯] with τf(V_+)=∞, τf(V¯)=0, then the LCT of 1V≤Vf(τ) is given by(21)LC1V≤Vfτ=z∫0+∞e-zτ1V≤Vfτdτ=z∫0+∞e-zτ1τ≤τfVdτ=z∫0τfVe-zτdτ=1-e-zτfV.Applying LCT to PDE (14) with respect to the time variable on the fixed region V∈[V_,+∞), we obtain that(22)AP^V,z-zP^V,z-AK-V1-e-zτfV+zmax⁡K-V,0=0,V≥V_,and LCT to conditions (15)-(17) give that(23)P^V_,z=K-V_,(24)∂P^V_,z∂V=-1,(25)limV→∞⁡P^V,z=0.The ODEs (22)-(25) cannot be solved directly as the unknown V_ and the function τf(V) are involved. In next paragraphs, we derive the value V_ and design an approximation for τf(V).Assume thatϕ(V) is a solution of (9) and satisfies condition (12). Then P∞(V)=Cϕ(V) with C being a constant is also the solution of (9) with (12). Using (10) and (11), we have(26)CϕVf∞=K-Vf∞,(27)CdϕVdV∣V=Vf∞=-1.This gives a nonlinear algebraic equation(28)ϕVf∞=Vf∞-KdϕVdV∣V=Vf∞,with unknown Vf(∞). Solving (28), we get the value V_=Vf(∞). In the following, we aim to solve ϕ(V) exactly that allows us to avoid using the iterative procedure in [16] for approximating τf(V) and optimizing NCIMs.Theorem 1.
Under GBMP model, the value ofVf(∞) is given by(29)Vf∞=Kγγ-1,(30)γ=-2μ-σ2-2μ-σ22+8σ2r4σ2.Proof.
The proof is provided by [20].Theorem 2.
Under MRGP model, the value ofP∞(V) is given by(31)P∞V=Cψr2λ,12,α-λV2λσ2,where ψ(·) is the degenerate hypergeometric function (see, e.g., [21]), and C is a constant to be determined.Proof.
Under MRGP model, the ODE (9) can be rewritten as(32)σ22∂2P∞∂V2+α-λV∂P∞∂V-rP∞=0.
Using the variable transformationy=α-λV, and the function transformation P~(y)=P∞(α-y)/λ, further inserting(33)∂P∞∂V=-λ∂P~∂y,(34)∂2P∞∂V2=λ2∂2P~∂y2,into (32) yields(35)λ2σ22∂2P~∂y2-λy∂P~∂y-rP~=0,whose basic solutions are just the degenerate hypergeometric functions ϕr/2λ,1/2,y2/λσ2 and ψr/2λ,1/2,y2/λσ2 (see, e.g., [21]).
Thus the general solution of (32) can be expressed as(36)P∞V=C1ϕr2λ,12,α-λV2λσ2+C2ψr2λ,12,α-λV2λσ2,where C1 and C2 are constants to be determined.
Noting thatlimV→∞ϕr/2λ,1/2,α-λV2/λσ2=+∞, from condition (12), we obtain that C1=0. Thus the general solution of (32) is(37)P∞V=C2ψr2λ,12,α-λV2λσ2,where we complete the proof.Theorem 3.
Under MRSRP model, the value ofP∞(V) is given by(38)P∞V=CeλV/2σ2λVσ2-1/4Wk,mλVσ2,where Wk,m(·) is the Whittaker function (see, e.g., [22]), and k=1/4-r/2λ,m=1/4, C is a constant to be determined.Proof.
Under MRSRP model, the ODE (9) can be rewritten as(39)2σ2V∂2P∞∂V2+σ2-2λV∂P∞∂V-rP∞=0.
Using the variable transformationy=λV/σ2, and the function transformation P~(y)=P∞(σ2y/λ), further inserting(40)∂P∞∂V=λσ2∂P~∂y,(41)∂2P∞∂V2=λ2σ4∂2P~∂y2,into (39) yields(42)2λy∂2P~∂y2+λ1-2y∂P~∂y-rP~=0,further using the function transformation P~(y)=y-1/4ey/2W(y), and inserting(43)∂P~∂y=y-1/4ey/212-14yWy+W′y,(44)∂2P~∂y2=y-1/4ey/2516y2-14y+14Wy+1-12yW′y+W′′y,into (42) yields(45)d2Wydy2+-14+ky+1/4-m2y2Wy=0,where(46)m=14,k=14-r2λ,whose basic solutions are just Whittaker functions Mk,m(y) and Wk,m(y) (see, e.g., [22]). Therefore the basic solutions of equation (39) are(47)ϕV=λVσ2-1/4eλV/2σ2Wk,mλVσ2,ψV=λVσ2-1/4eλV/2σ2Mk,mλVσ2,and the general solution of (39) can be expressed as(48)P∞V=C1ϕV+C2ψV,where C1 and C2 are constants to be determined.
Noting thatlimV→∞ψ(V)=+∞, from condition (12), we obtain that C2=0. Thus the general solution of (39) is(49)P∞V=C1eλV/2σ2λVσ2-1/4Wk,mλVσ2,where we complete the proof.Theorem 4.
Under MRLP model, the value ofP∞(V) is given by(50)P∞V=Cψr2λ,12,λln⁡V-α+σ2/22λσ2,where ψ(·) is the degenerate hypergeometric function, and C is a constant to be determined.Proof.
Under MRLP model, the ODE (9) can be rewritten as(51)σ2V22∂2P∞∂V2+α-λln⁡VV∂P∞∂V-rP∞=0.
Using the variable transformationy=λln⁡V-α+σ2/2, and the function transformation P~(y)=P∞ey-σ2/2+α/λ, further inserting(52)∂P∞∂V=λey-σ2/2+α/λ∂P~∂y,(53)∂2P∞∂V2=λ2e2y-σ2/2+α/λ∂2P~∂y2-λe2y-σ2/2+α/λ∂P~∂y,into (51) yields(54)σ2λ22∂2P~∂y2-λy∂P~∂y-rP~=0,whose basic solutions are just the degenerate hypergeometric functions ϕr/2λ,1/2,y2/λσ2 and ψr/2λ,1/2,y2/λσ2. Thus the general solution of (51) can be expressed as(55)P∞V=C1ϕr2λ,12,λln⁡V-α+σ2/22λσ2+C2ψr2λ,12,λln⁡V-α+σ2/22λσ2,where C1 and C2 are constants to be determined.
Noting thatlimV→∞ϕr/2λ,1/2,λln⁡V-α+σ2/22/λσ2=+∞, from condition (12), we obtain that C1=0. Thus the general solution of equation (51) is(56)P∞V=C2ψr2λ,12,λln⁡V-α+σ2/22λσ2,where we complete the proof.Recall again thatτf(V) is the strictly decreasing function on the interval (V_,V¯] with τf(V_+)=∞, τf(V¯)=0. Thus it is similar to [16], and we construct an approximation of τf(V)(57)τfV=0,V≥V¯,l1-ln⁡V-V_V¯-V_l2,V_<V<V¯,+∞,V=V_,where lk,k=1,2 are positive numbers to be determined by the following optimization.Now we adopt the finite difference method (FDM) to solve the ODE (22) with boundary conditions (23)-(25). We define uniform mesh Vi=V_+iΔV,i=0,1,…,N with VN=Vmax, where Vmax is a sufficiently larger number such that P^(Vmax,z)≈0. Let P^i be the approximation of P^(Vi,z) at mesh point V=Vi for any complex value z, i.e., P^i(z)≈P^(Vi,z); then the ODE (22) can be discretized as, for i=1,2,…,N-1,(58)A~P^iz-zP^iz=AK-Vi1-e-zτfVi-zmax⁡K-Vi,0,with boundary conditions (59)P^0=K-V_,P^N=0,where the difference operators A~ are defined in the second column of Table 3.Table 3
Difference operators.
Model A ~ P ^ i ( z ) GBMP 2 σ 2 V i 2 P ^ i + 1 - 2 P ^ i + P ^ i - 1 Δ V 2 + ( 2 μ + σ 2 ) V i P ^ i + 1 - P ^ i - 1 2 Δ V - r P ^ i MRGP σ 2 2 P ^ i + 1 - 2 P ^ i + P ^ i - 1 Δ V 2 + ( α - λ V i ) P ^ i + 1 - P ^ i - 1 2 Δ V - r P ^ i MRSRP 2 σ 2 V i P ^ i + 1 - 2 P ^ i + P ^ i - 1 Δ V 2 + ( σ 2 - 2 λ V i ) P ^ i + 1 - P ^ i - 1 2 Δ V - r P ^ i MRLP σ 2 V i 2 2 P ^ i + 1 - 2 P ^ i + P ^ i - 1 Δ V 2 + ( α - λ ln ⁡ V i ) V i P ^ i + 1 - P ^ i - 1 2 Δ V - r P ^ iUsing (59), the linear system (58) can be expressed as in the matrix form(60)A-zIP^z=gz,where the coefficient matrix A is defined as(61)A=b1c10⋯0a2b2c2⋯0⋮⋱⋱⋱⋮0…aN-2bN-2cN-20…0aN-1bN-1,moreover the solution vector P^(z) and the constant vector g(z) are defined as, respectively,(62)P^z=P^1z,P^2z,⋯,P^N-2z,P^N-1z′,(63)gz=g1z-a1P^0,g2z,⋯,gN-2z,gN-1z-cN-1P^N′,where elements in the matrix A are listed in Table 4, and(64)giz=AK-Vi1-e-zτfVi-zmax⁡K-Vi,0.Table 4
Coefficients.
Model a i b i c i GBMP 2 σ 2 V i 2 Δ V 2 - ( 2 μ + σ 2 ) V i 2 Δ V - 4 σ 2 V i 2 Δ V 2 - r 2 σ 2 V i 2 Δ V 2 + ( 2 μ + σ 2 ) V i 2 Δ V MRGP σ 2 2 Δ V 2 - α - λ V i 2 Δ V - σ 2 Δ V 2 - r σ 2 2 Δ V 2 + α - λ V i 2 Δ V MRSRP 2 σ 2 V i Δ V 2 - σ 2 - 2 λ V i 2 Δ V - 4 σ 2 V i Δ V 2 - r 2 σ 2 V i Δ V 2 + σ 2 - 2 λ V i 2 Δ V MRLP σ 2 V i 2 2 Δ V 2 - ( α - λ ln ⁡ V i ) V i 2 Δ V - σ 2 V i 2 Δ V 2 - r σ 2 V i 2 2 Δ V 2 + ( α - λ ln ⁡ V i ) V i 2 Δ VLeta(z) be the first row of (A-zI)-1; then a(z) is given by the following linear system:(65)azA-zI=1,0,0,…,0.Thus the first element in the solution vector is given by(66)P^1z=azgz;l1,l2,where g(z;l1,l2)≔g(z), and the approximation of τf(Vi) involved in g(z), for i=1,2,…,N-1, is given by (57).We find the parameterslk,k=1,2 such that the condition (24) holds for all real values of z>0. A practicable way given by Zhou, Ma, and Sun [16] is to consider a sequence of z1,z2,…,zm(m≥2) and find the values of unknown lk,k=1,2 that minimize the sum of residue(67)∑j=1m∂P^V_,zj∂V--1≈∑j=1mP^V1,zj-K-V_ΔV+1.Inserting (66) into (67), we can find the values of parameters {lk},k=1,2 by solving the following discrete optimization problem, for given real positive numbers z1,z2,…,zm(m≥2),(68)minl1,l2≥0∑j=1mazjgzj;l1,l2-K-V_ΔV+1,for approximation formula (57). After carrying out the optimization procedure (68), then we can find the approximation of the early exercise boundary Vf(τ) as the inverse of the function τf(V).After solving the linear system (60) and only by recognizing the relationship between LCT and LT (20), we can use Laplace inversion to recover the nodal values, for i=0,1,…,N,(69)PVi,τ=L-1P^Vi,zz,and apply the cubic spline interpolation to obtain the option value function P(V,τ).In the following, the hyperbolic contour integral methods for Laplace inversion are studied to restore the option prices. Denote(70)V=V1,V2,…,VN-1′,(71)Pτ≔PV1,τ,PV2,τ,…,PVN-1,τ′,(72)P^z≔P^V1,z,P^V2,z,…,P^VN-1,z′.And the term g(z) can be rewritten as two parts(73)gz=g1z+g2z,With(74)g1z=g1,1z-a1P^0,g1,2z,…,g1,N-2z,g1,N-1z-cN-1P^N′,(75)g2z=g2,1z,g2,2z,…,g2,N-2z,g2,N-1z′,where(76)g1,iz=-AK-Vie-zτfVi1τ>τfVi+AK-Vi-zmax⁡K-Vi,0,and(77)g2,iz=-AK-Vie-zτfVi1τ≤τfVi,for i=1,2,…,N-1.From (76) and (77), we can see that ezτg1(z)/z→0 as Re(z)→-∞ and ezτg2(z)/z→0 as Re(z)→+∞, which guarantee that the following Laplace inversion formula is feasible (see [16]):(78)Pτ=P1τ+P2τ=12πi∫Γ1ezτP^1zzdz+12πi∫Γ2ezτP^2zzdz,=12πi∫Γ1ezτA-zI-1g1zzdz+12πi∫Γ2ezτA-zI-1g2zzdz,where the contours Γ1 and Γ2 should be chosen as two special Hankel contours as in Figure 1 which enclose all singularities of P^1(z)/z and P^1(z)/z, respectively. From expression (78), we can see that ezτP^1(z)/z→0 as Re(z)→-∞ and ezτP^2(z)/z→0 as Re(z)→+∞ provided that all singularities of (A-zI)-1g1(z)/z and (A-zI)-1g2(z)/z lie in the left-half plane of the contour Γ1 or equivalently the spectrum Λ(A) lies in a sectorial region Σδ that will be proved in the following.Figure 1
An example of Hankel contours.Following [16, 23], the parameterized hyperbolic contour Γ1 on the left-hand side is selected as(79)Γ1:zζ=ν1+sin⁡iζ-θ,-∞<ζ<∞,where parameters ν>0 and θ set the width and the asymptotic angle of the hyperbolic contour, respectively.Substituting the contour (79) into (78) gives(80)P1τ=12πi∫Γ1ezτP^1zzdz=12πi∫-∞+∞ezζτz′ζP^1zζzζdζ≈h2πi∑m=-LLezζmτz′ζmP^1zζmzζm,where ζm=mh for the trapezoidal rule and L is the number of quadrature nodes. The parameters h, ν are given by (see, e.g., [23])(81)h=ρθL,ν=4πθ-π2+2πδρθLτ,ρθ=cosh-1⁡2θ4θ-π+2δsin⁡θ,which ensure that all singularities of the integrand in (80) lie in a sectorial region(82)Σδ=z∈C:arg⁡-z≤δ,δ∈0,π2,where θ is a free parameter. The choice of the parameters (81) leads to the following predicted convergence rate:(83)P1τ-h2πi∑m=-LLezζmτz′ζmP^1zζmzζm≤Oe-RθL,where(84)Rθ=π2-2πθ-2πδρθ,and · is some vector norm. According to expressions (81)-(84), after the semiangle δ of the sectorial region (82) is given, then we can get the optimal θ by maximizing the function R(θ), and thus the optimal parameters h and ν can be computed from the formulas (81).The parameterized hyperbolic contourΓ2 on the right-hand side is selected as(85)Γ2:zζ=ν1+sin⁡iζ+θ,-∞<ζ<∞,where parameters h and ν are given by formulas (81) with setting δ=0 and(86)P2τ≈h2πi∑m=-LLezζmτz′ζmP^2zζmzζm.To analyze the spectrum of the coefficient matrix (61), we follow the idea in [24] and construct a Toeplitz matrix A¯ to replace analyzing (61)(87)A¯=b¯c¯0⋯0a¯b¯c¯⋯0⋮⋱⋱⋱⋮0…a¯b¯c¯0⋯0a¯b¯,where the elements in the matrix A¯ are defined as(88)a¯=aΔV2-b2ΔV,b¯=-2aΔV2-r,c¯=aΔV2+b2ΔV,where a and b are the mean values of the coefficients of diffusion term and convection term of the studied PDEs on the truncated domain [V_,Vmax], respectively. With the constructed Toeplitz matrix, we can analyze the spectrum Λ(A) and estimate the parameters of hyperbolic contour roughly.Theorem 5.
The numerical rangeW(A¯) of A¯ defined by (87) lies in a sectorial region Σδ such that(89)WA¯⊆Σδ≡z∈C:arg⁡-z≤δ,with(90)δ=arctan⁡b2ar.And(91)ΛA¯⊆W′A¯⊆Σδ.Proof.
Denotef(η) the generating function of A¯ and Ωf(η)=f(η)∣η∈(-π,π). We shall use the following two conclusions (see, e.g., [8, 16]): the numerical range W(A¯) is a subset of the closure of convex hull of Ωf(η), i.e., W(A¯)⊆conv(Ω(f))¯; let A,D∈Cn×n and suppose that D is positive definite; then the spectrum Λ(DA) is a subset of the angular numerical range of A, i.e., ΛDA⊆W′(A).
The generating function can be written as(92)fη=b¯ei0η+c¯e-iη+a¯eiη=b¯+a¯+c¯cos⁡η+ia¯-c¯sin⁡η,inserting (88) into (92) yields(93)fη=2aΔV2cos⁡η-1-r-ibΔVsin⁡η,η∈-π,0∪0,π,further using (93), we have(94)Im⁡fηRe⁡fη=b/ΔVsin⁡ηr+1-cos⁡η2a/ΔV2=b/ΔV2sin⁡η/2cos⁡η/2r+2sin2⁡η/22a/ΔV2=b/ΔVcos⁡η/2r/2sin⁡η/2+2sin⁡η/2a/ΔV2≤b2ar,which shows that Ωf(η)⊆Σδ. Apparently the sectorial region Σδ is a closed convex hull; thus we have(95)WA¯⊆convΩf¯⊆Σδ.Since W′(A¯) is also a sectorial region whose angle is just the angular opening of the smallest angular sector that includes W(A¯), thus we have(96)ΛA¯⊆W′A¯⊆Σδ,where we complete the proof.With Theorem5, the Laplace inversion formulas (80) and (86) can be properly used.


## 4. Numerical Examples

In this section, we use NCIMs for solving the free-boundary PDEs (2)-(6) under four stochastic volatility processes: GBMP, MRGP, MRSRP, and MRLP. The results are compared with that by the FDMs [25]; moreover the relative error (RE) criterion is adopted to measure the computational accuracy. Relative error is defined by (see, e.g., [26])(97)RE≔PNCIM-PFDMPFDM,In the test, we setVmax=4K and T=2, the space mesh partition number N=4000 for FDMs and NCIMs. Moreover the FDMs, when time mesh sizes equal T/2000, are taken as the benchmark methods. The positive values zk=k(k=1,2,…,m) with m=4 (see expression (68)). The Matlab command “fminsearch” is used to search the approximate parameters lk, k=1,2, with initial values lk=1,k=1,2, besides we use the technique in [16] to get the optimal V_ for MRGP and MRLP models since it achieves more accurate results. Moreover we compute the option prices in one Table, plot one figure including two subfigures, i.e., the early exercise boundaries (left), NCIMs errors versus L (right) in the sense of L2-norm, for the given volatility model. The L2-norm is defined as(98)ξL2≔1N+1∑i=0Nξi21/2.The codes are run in MATLAB R2014a on a PC with the configuration: AMD, CPU A10-9600P@2.40 GHz, and 24.0 GB RAM.


### 4.1. Implementations for GBMP Model

In this test, the model parameters are taken asμ=0.4, σ=0.5, r=0.05, and K=1. The numerics in Table 5 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.25% in most cases, whereas the CPU time for NCIM is less than the FDM, and the reasons for NCIM achieving good performance have two aspects: on the one hand, the NCIM has the exponential-order convergence rate with respect to the number of the quadrature nodes L for the hyperbolic contour integral (see formula (83)); on the other hand, there is no iterative procedure for numerical finding the free-boundary value of the corresponding perpetual American volatility put; thus it is not like [16]. In Figure 2 (left), we can see the boundary obtained by NCIM is very close to that by FDM. Moreover in Figure 2 (right), we can see that the error of the NCIM with respect to L roughly behaves proportional to the theoretical ones exp⁡(-0.2313L), which means that the error exponentially decays with respect to L.Table 5
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=14 value value RE (%) 0.5 0.50000 0.50000 0.50000 0.50000 0.50000 0.00000 1 0.19484 0.19489 0.19491 0.19493 0.19441 0.26676 1.5 0.10127 0.10132 0.10135 0.10136 0.10113 0.22691  CPU times (s) 315 453 732 1335 186 –Figure 2
Early exercise boundaries (left) under GBMP model withT/2000 for FDM and L=14 for NCIM, NCIM errors in log-scale versus L (right).


### 4.2. Implementations for MRGP Model

In this test, the model parameters are taken asα=1.6, σ=1, λ=2, r=0.05, and K=1. The numerics in Table 6 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.8% in most cases, whereas the CPU time for NCIM is much less than the FDM. In Figure 3 (left), we can see the boundary obtained by NCIM almost matches that by FDM; moreover the realized convergence rate of NCIM roughly behaves proportional to the theoretical ones from the right subfigure.Table 6
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=6 value value RE (%) 0.5 0.69600 0.69624 0.69636 0.69642 0.69152 0.70217 1 0.58613 0.58649 0.58668 0.58677 0.58206 0.80270 1.5 0.52618 0.52663 0.52685 0.52696 0.52280 0.78755  CPU times (s) 613 785 1141 1808 87 –Figure 3
Early exercise boundaries (left) under MRGP model withT/2000 for FDM and L=6 for NCIM, NCIM errors in log-scale versus L (right).


### 4.3. Implementations for MRSRP Model

In this test, the model parameters are taken asσ=1.4, λ=0.001, r=0.05, and K=5. The numerics in Table 7 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.8% in most cases, whereas the CPU time for NCIM is less than the FDM. In Figure 4 (left), we can see the boundary obtained by NCIM is very close to that by FDM; moreover the realized convergence rate of NCIM roughly behaves proportional to the theoretical ones from the right subfigure.Table 7
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=10 value value RE (%) 4 2.46906 2.46997 2.47044 2.47068 2.45192 0.75931 5 2.15034 2.15138 2.15189 2.15216 2.13482 0.80570 6 1.87964 1.88074 1.88129 1.88158 1.86568 0.84503  CPU times (s) 551 716 1045 1729 95 –Figure 4
Early exercise boundaries (left) under MRSRP model withT/2000 for FDM and L=10 for NCIM, NCIM errors in log-scale versus L (right).


### 4.4. Implementations for MRLP Model

In this test, the model parameters are taken asα=0.6, λ=0.1, σ=1, r=0.05, and K=5. The numerics in Table 8 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.9% in most cases, whereas the CPU time for NCIM is less than the FDM. In Figure 5 (left), we can see the boundary obtained by NCIM is very close to that by FDM; moreover the realized convergence rate of NCIM roughly behaves proportional to the theoretical ones from the right subfigure.Table 8
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=14 value value RE (%) 4 1.94291 1.94363 1.94399 1.94418 1.93598 0.42177 5 1.63948 1.64031 1.64073 1.64095 1.62844 0.76236 6 1.40529 1.40617 1.40662 1.40684 1.39302 0.98234  CPU times (s) 431 582 1519  128 –Figure 5
Early exercise boundaries (left) under MRLP model withT/2000 for FDM and L=14 for NCIM, NCIM errors in log-scale versus L (right).


## 4.1. Implementations for GBMP Model

In this test, the model parameters are taken asμ=0.4, σ=0.5, r=0.05, and K=1. The numerics in Table 5 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.25% in most cases, whereas the CPU time for NCIM is less than the FDM, and the reasons for NCIM achieving good performance have two aspects: on the one hand, the NCIM has the exponential-order convergence rate with respect to the number of the quadrature nodes L for the hyperbolic contour integral (see formula (83)); on the other hand, there is no iterative procedure for numerical finding the free-boundary value of the corresponding perpetual American volatility put; thus it is not like [16]. In Figure 2 (left), we can see the boundary obtained by NCIM is very close to that by FDM. Moreover in Figure 2 (right), we can see that the error of the NCIM with respect to L roughly behaves proportional to the theoretical ones exp⁡(-0.2313L), which means that the error exponentially decays with respect to L.Table 5
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=14 value value RE (%) 0.5 0.50000 0.50000 0.50000 0.50000 0.50000 0.00000 1 0.19484 0.19489 0.19491 0.19493 0.19441 0.26676 1.5 0.10127 0.10132 0.10135 0.10136 0.10113 0.22691  CPU times (s) 315 453 732 1335 186 –Figure 2
Early exercise boundaries (left) under GBMP model withT/2000 for FDM and L=14 for NCIM, NCIM errors in log-scale versus L (right).


## 4.2. Implementations for MRGP Model

In this test, the model parameters are taken asα=1.6, σ=1, λ=2, r=0.05, and K=1. The numerics in Table 6 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.8% in most cases, whereas the CPU time for NCIM is much less than the FDM. In Figure 3 (left), we can see the boundary obtained by NCIM almost matches that by FDM; moreover the realized convergence rate of NCIM roughly behaves proportional to the theoretical ones from the right subfigure.Table 6
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=6 value value RE (%) 0.5 0.69600 0.69624 0.69636 0.69642 0.69152 0.70217 1 0.58613 0.58649 0.58668 0.58677 0.58206 0.80270 1.5 0.52618 0.52663 0.52685 0.52696 0.52280 0.78755  CPU times (s) 613 785 1141 1808 87 –Figure 3
Early exercise boundaries (left) under MRGP model withT/2000 for FDM and L=6 for NCIM, NCIM errors in log-scale versus L (right).


## 4.3. Implementations for MRSRP Model

In this test, the model parameters are taken asσ=1.4, λ=0.001, r=0.05, and K=5. The numerics in Table 7 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.8% in most cases, whereas the CPU time for NCIM is less than the FDM. In Figure 4 (left), we can see the boundary obtained by NCIM is very close to that by FDM; moreover the realized convergence rate of NCIM roughly behaves proportional to the theoretical ones from the right subfigure.Table 7
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=10 value value RE (%) 4 2.46906 2.46997 2.47044 2.47068 2.45192 0.75931 5 2.15034 2.15138 2.15189 2.15216 2.13482 0.80570 6 1.87964 1.88074 1.88129 1.88158 1.86568 0.84503  CPU times (s) 551 716 1045 1729 95 –Figure 4
Early exercise boundaries (left) under MRSRP model withT/2000 for FDM and L=10 for NCIM, NCIM errors in log-scale versus L (right).


## 4.4. Implementations for MRLP Model

In this test, the model parameters are taken asα=0.6, λ=0.1, σ=1, r=0.05, and K=5. The numerics in Table 8 show that the prices obtained by NCIM and FDM benchmark are very close and the relative errors are less than 0.9% in most cases, whereas the CPU time for NCIM is less than the FDM. In Figure 5 (left), we can see the boundary obtained by NCIM is very close to that by FDM; moreover the realized convergence rate of NCIM roughly behaves proportional to the theoretical ones from the right subfigure.Table 8
Comparison results (option values, relative errors and CPU time).
V 0 FDM NCIM T/250 T/500 T/1000 T/2000 L=14 value value RE (%) 4 1.94291 1.94363 1.94399 1.94418 1.93598 0.42177 5 1.63948 1.64031 1.64073 1.64095 1.62844 0.76236 6 1.40529 1.40617 1.40662 1.40684 1.39302 0.98234  CPU times (s) 431 582 1519  128 –Figure 5
Early exercise boundaries (left) under MRLP model withT/2000 for FDM and L=14 for NCIM, NCIM errors in log-scale versus L (right).


## 5. Conclusions

This paper studied the NCIMs for solving free-boundary PDEs from American volatility options pricing. The free-boundary PDEs could be written as a unified form of PDEs on a fixed space region, and performing the Laplace transform gave the parameterized ODEs with unknown inverse function of free boundary, then the value of early exercise boundary of the perpetual American volatility put was solved which enabled the inverse function of free boundary to be well approximated, finally the FDM was adopted to solve the ODEs, and the results were restored via numerical Laplace inversion. Numerical comparisons of the NCIMs with the FDMs were conducted to verify the effectiveness of the method.



---

*Source: 1838521-2018-11-04.xml*