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

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.


Introduction
The pricing problems of European-style options written on volatility are widely studied in the literature by deriving the explicit formulas [1][2][3][4][5][6], by using the simulation approaches [7,8], by using the PDE approaches [5,[9][10][11][12][13], and by the empirical analysis [14].However, the valuation of Americanstyle 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 freeboundary functions.
The remaining parts of this paper are arranged as follows: In Section 2, 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.

Model Descriptions
Let   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.
In Table 1, Z denotes the Brownian motion process under risk neutral measure P, , , ,  (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 the price of American-style put option with underlying stochastic volatility   , current time , maturity , and strike price , where T 0, is the optimal stopping set with the Brownian motion filtration for  ∈ [0, ].
Using transformation of time variable  =  −  (namely, the time to maturity), then the option price and the early exercise boundary (namely, the free-boundary function) are denoted by (, ) fl (,  − ) = (, ) and   (), respectively.Then the valuation of American volatility put option can be formulated as a free-boundary problem (, 0) = max ( − , 0) , where we can verify that where 1 {⋅} represents the indicator function and the values of A(−) for different models are listed in the second column of Table 2. Moreover the early exercise boundaries at initial   (0) are given in the third column of Table 2, where  * is the unique positive root of the following nonlinear algebraic equation: It is proved by [17,18] for GBMP model and [18] for MRGP model, MRSRP model, and MRLP model that the free-boundary functions   () are continuously differentiable, strictly decreasing and convex on the interval  ∈ [0, +∞), which are confirmed by the simulation results in [19].Moreover we note that the free-boundary functions   () are bounded on the interval (, ], where  =   (∞) are the early exercise boundaries of the corresponding perpetual American volatility put options, whose values satisfy the following ODEs: lim where the differential operators A are listed in the third column of Table 1 with replacement of  by  ∞ .

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).
Now we focus on the Laplace-Carson transform (LCT) methods for solving the unified PDE (14).For any complex value , the LCT is defined as (19) 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 Let   () be the inverse function of   ().Since   () is a strictly decreasing function on the interval (, ] with   (+) = ∞,   () = 0, then the LCT of 1 {≤  ()} is given by Applying LCT to PDE (14) with respect to the time variable on the fixed region  ∈ [, +∞), we obtain that and LCT to conditions ( 15)- (17) give that The ODEs ( 22)-( 25) cannot be solved directly as the unknown  and the function   () are involved.In next paragraphs, we derive the value  and design an approximation for   ().
Assume that () is a solution of ( 9) and satisfies condition (12).Then  ∞ () = () with  being a constant is also the solution of ( 9) with (12).Using (10) and (11), we have This gives a nonlinear algebraic equation with unknown   (∞).Solving (28), we get the value  =   (∞).In the following, we aim to solve () exactly that allows us to avoid using the iterative procedure in [16] for approximating   () and optimizing NCIMs.
Theorem 1.Under GBMP model, the value of   (∞) is given by Proof.The proof is provided by [20].
Theorem 2. Under MRGP model, the value of  ∞ () is given by where (⋅) is the degenerate hypergeometric function (see, e.g., [ ]), and  is a constant to be determined.
Theorem 4.Under MRLP model, the value of  ∞ () is given by where (⋅) is the degenerate hypergeometric function, and  is a constant to be determined.
Recall again that   () is the strictly decreasing function on the interval (, ] with   (+) = ∞,   () = 0. Thus it is similar to [16], and we construct an approximation of   () where ℓ  ,  = 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   =  + Δ,  = 0, 1, . . .,  with   =  max , where  max is a sufficiently larger number such that P( max , ) ≈ 0. Let P be the approximation of P(  , ) at mesh point  =   for any complex value , i.e., P () ≈ P(  , ); then the ODE (22) where the difference operators Ã are defined in the second column of Table 3.Using (59), the linear system (58) can be expressed as in the matrix form where the coefficient matrix A is defined as moreover the solution vector P() and the constant vector g() are defined as, respectively, where elements in the matrix A are listed in Table 4, and Let a() be the first row of (A − I) −1 ; then a() is given by the following linear system: Thus the first element in the solution vector is given by where g(; ℓ 1 , ℓ 2 ) fl g(), and the approximation of   (  ) involved in g(), for  = 1, 2, . . .,  − 1, is given by (57).We find the parameters {ℓ  },  = 1, 2 such that the condition (24) holds for all real values of  > 0. A practicable way given by Zhou, Ma, and Sun [16] for approximation formula (57).After carrying out the optimization procedure (68), then we can find the approximation of the early exercise boundary   () as the inverse of the function   ().
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  = 0, 1, . . ., , and apply the cubic spline interpolation to obtain the option value function (, ).
In the following, the hyperbolic contour integral methods for Laplace inversion are studied to restore the option prices.Denote And the term g() can be rewritten as two parts With where and for  = 1, 2, . . .,  − 1.From ( 76) and (77), we can see that   |g 1 ()/| → 0 as Re() → −∞ and   |g 2 ()/| → 0 as Re() → +∞, which guarantee that the following Laplace inversion formula is feasible (see [16]): where the contours Γ 1 and Γ 2 should be chosen as two special Hankel contours as in Figure 1 which enclose all singularities of P1 ()/ and P1 ()/, respectively.From expression (78), we can see that   P1 ()/z → 0 as Re() → −∞ and   P2 ()/ → 0 as Re() → +∞ provided that all singularities of (A − I) −1 g 1 ()/ and (A − I) −1 g 2 ()/ 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.Following [16,23], the parameterized hyperbolic contour Γ 1 on the left-hand side is selected as where parameters ] > 0 and  set the width and the asymptotic angle of the hyperbolic contour, respectively.Substituting the contour (79) into (78) gives where   = ℎ for the trapezoidal rule and  is the number of quadrature nodes.The parameters ℎ, ] are given by (see, e.g., [23]) which ensure that all singularities of the integrand in (80) lie in a sectorial region where  is a free parameter.The choice of the parameters (81) leads to the following predicted convergence rate: where 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 (), and thus the optimal parameters ℎ and ] can be computed from the formulas (81).The parameterized hyperbolic contour Γ 2 on the righthand side is selected as where parameters ℎ and ] are given by formulas (81) with setting  = 0 and 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) where the elements in the matrix A are defined as where  and  are the mean values of the coefficients of diffusion term and convection term of the studied PDEs on the truncated domain [,  max ], respectively.With the constructed Toeplitz matrix, we can analyze the spectrum Λ(A) and estimate the parameters of hyperbolic contour roughly.
Theorem 5. e numerical range W(A) of A defined by ( ) lies in a sectorial region Σ  such that with ) . And Proof.Denote () the generating function of A and Ω(()) = {() |  ∈ (−, )}.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 Ω(()), i.e., W(A) ⊆ conv(Ω()); let A, D ∈ C × 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 inserting (88) into (92) yields which shows that Ω(()) ⊆ Σ  .Apparently the sectorial region Σ  is a closed convex hull; thus we have 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 where we complete the proof.

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]) In the test, we set  max = 4 and  = 2, the space mesh partition number  = 4000 for FDMs and NCIMs.Moreover the FDMs, when time mesh sizes equal /2000, are taken as the benchmark methods.The positive values   =  ( = 1, 2, . . ., ) with  = 4 (see expression (68)).The Matlab command "fminsearch" is used to search the approximate parameters ℓ  ,  = 1, 2, with initial values ℓ  = 1,  = 1, 2, besides we use the technique in [16] to get the optimal  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  (right) in the sense of  2 -norm, for the given volatility model.The  2 -norm is defined as The codes are run in MATLAB R2014a on a PC with the configuration: AMD, CPU A10-9600P@2.40GHz, and 24.0 GB RAM.
. .Implementations for GBMP Model.In this test, the model parameters are taken as  = 0.4,  = 0.5,  = 0.05, and  = 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  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  roughly behaves proportional to the theoretical ones exp(−0.2313),which means that the error exponentially decays with respect to .
. .Implementations for MRGP Model.In this test, the model parameters are taken as  = 1.6,  = 1,  = 2,  = 0.05, and  = 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.
. .Implementations for MRSRP Model.In this test, the model parameters are taken as  = 1.4, = 0.001,  = 0.05, and  = 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.
. .Implementations for MRLP Model.In this test, the model parameters are taken as  = 0.6,  = 0.1,  = 1,  = 0.05, and  = 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.

Conclusions
This paper studied the NCIMs for solving free-boundary PDEs from American volatility options pricing.The freeboundary 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.

Table 5 :
Comparison results (option values, relative errors and CPU time).

Table 6 :
Comparison results (option values, relative errors and CPU time).

Table 7 :
Comparison results (option values, relative errors and CPU time).

Table 8 :
Comparison results (option values, relative errors and CPU time).