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.
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 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.
2. 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, denotes the Brownian motion process under risk neutral measure , , , , (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 bythe price of American-style put option with underlying stochastic volatility , current time , maturity , and strike price , where is the optimal stopping set with the Brownian motion filtration for .
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 and , respectively. Then the valuation of American volatility put option can be formulated as a free-boundary problemwhere we can verify thatwhere represents the indicator function and the values of for different models are listed in the second column of Table 2. Moreover the early exercise boundaries at initial 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 , 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:where the differential operators are listed in the third column of Table 1 with replacement of by .
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 for , then the value of American volatility put is given by , and we can verify that it is equivalent to the following PDE:
Combining (2) and (13) gives a unified form of PDE:
Recall that for all ; the value of American volatility put option equals on the region . Therefore it suffices to solve the PDE (14) with conditions
Now we focus on the Laplace-Carson transform (LCT) methods for solving the unified PDE (14). For any complex value , the LCT is defined aswhich 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 , , then the LCT of is given by
Applying LCT to PDE (14) with respect to the time variable on the fixed region , we obtain thatand LCT to conditions (15)-(17) give thatThe 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 haveThis gives a nonlinear algebraic equationwith 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 bywhere is the degenerate hypergeometric function (see, e.g., [21]), and is a constant to be determined.
Proof. Under MRGP model, the ODE (9) can be rewritten asUsing the variable transformation , and the function transformation , further insertinginto (32) yieldswhose basic solutions are just the degenerate hypergeometric functions and (see, e.g., [21]).
Thus the general solution of (32) can be expressed aswhere and are constants to be determined.
Noting that , from condition (12), we obtain that . Thus the general solution of (32) iswhere we complete the proof.
Theorem 3. Under MRSRP model, the value of is given bywhere is the Whittaker function (see, e.g., [22]), and , is a constant to be determined.
Proof. Under MRSRP model, the ODE (9) can be rewritten asUsing the variable transformation , and the function transformation , further insertinginto (39) yieldsfurther using the function transformation , and insertinginto (42) yieldswherewhose basic solutions are just Whittaker functions and (see, e.g., [22]). Therefore the basic solutions of equation (39) areand the general solution of (39) can be expressed aswhere and are constants to be determined.
Noting that , from condition (12), we obtain that . Thus the general solution of (39) iswhere we complete the proof.
Theorem 4. Under MRLP model, the value of is given bywhere is the degenerate hypergeometric function, and is a constant to be determined.
Proof. Under MRLP model, the ODE (9) can be rewritten asUsing the variable transformation , and the function transformation , further insertinginto (51) yieldswhose basic solutions are just the degenerate hypergeometric functions and . Thus the general solution of (51) can be expressed aswhere and are constants to be determined.
Noting that , from condition (12), we obtain that . Thus the general solution of equation (51) iswhere we complete the proof.
Recall again that is the strictly decreasing function on the interval with , . Thus it is similar to [16], and we construct an approximation of where 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 with , where is a sufficiently larger number such that . Let be the approximation of at mesh point for any complex value , i.e., ; then the ODE (22) can be discretized as, for ,with boundary conditions 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 formwhere the coefficient matrix is defined asmoreover the solution vector and the constant vector are defined as, respectively,where elements in the matrix are listed in Table 4, and
Let be the first row of ; then is given by the following linear system:Thus the first element in the solution vector is given bywhere , and the approximation of involved in , for , is given by (57).
We find the parameters such that the condition (24) holds for all real values of . A practicable way given by Zhou, Ma, and Sun [16] is to consider a sequence of and find the values of unknown that minimize the sum of residue
Inserting (66) into (67), we can find the values of parameters by solving the following discrete optimization problem, for given real positive numbers ,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 ,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. DenoteAnd the term can be rewritten as two partsWithwhereandfor .
From (76) and (77), we can see that as and as , which guarantee that the following Laplace inversion formula is feasible (see [16]):where the contours and should be chosen as two special Hankel contours as in Figure 1 which enclose all singularities of and , respectively. From expression (78), we can see that as and as provided that all singularities of and lie in the left-half plane of the contour or equivalently the spectrum lies in a sectorial region that will be proved in the following.
Following [16, 23], the parameterized hyperbolic contour on the left-hand side is selected aswhere parameters and set the width and the asymptotic angle of the hyperbolic contour, respectively.
Substituting the contour (79) into (78) giveswhere 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 regionwhere is a free parameter. The choice of the parameters (81) leads to the following predicted convergence rate:whereand 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 on the right-hand side is selected aswhere parameters and are given by formulas (81) with setting and
To analyze the spectrum of the coefficient matrix (61), we follow the idea in [24] and construct a Toeplitz matrix to replace analyzing (61)where the elements in the matrix are defined aswhere and are the mean values of the coefficients of diffusion term and convection term of the studied PDEs on the truncated domain , respectively. With the constructed Toeplitz matrix, we can analyze the spectrum and estimate the parameters of hyperbolic contour roughly.
Theorem 5. The numerical range of defined by (87) lies in a sectorial region such thatwithAnd
Proof. Denote the generating function of and . We shall use the following two conclusions (see, e.g., [8, 16]): the numerical range is a subset of the closure of convex hull of , i.e., ; let and suppose that is positive definite; then the spectrum is a subset of the angular numerical range of , i.e., .
The generating function can be written asinserting (88) into (92) yieldsfurther using (93), we havewhich shows that . Apparently the sectorial region is a closed convex hull; thus we haveSince is also a sectorial region whose angle is just the angular opening of the smallest angular sector that includes , thus we havewhere we complete the proof.
With Theorem 5, 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])
In the test, we set and , the space mesh partition number for FDMs and NCIMs. Moreover the FDMs, when time mesh sizes equal , are taken as the benchmark methods. The positive values with (see expression (68)). The Matlab command “fminsearch” is used to search the approximate parameters , , with initial values , 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 -norm, for the given volatility model. The -norm is defined as
The codes are run in MATLAB R2014a on a PC with the configuration: AMD, CPU [email protected] GHz, and 24.0 GB RAM.
4.1. Implementations for GBMP Model
In this test, the model parameters are taken as , , , and . 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 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 , which means that the error exponentially decays with respect to .
4.2. Implementations for MRGP Model
In this test, the model parameters are taken as , , , , and . 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 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.
4.3. Implementations for MRSRP Model
In this test, the model parameters are taken as , , , and . 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 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.
4.4. Implementations for MRLP Model
In this test, the model parameters are taken as , , , , and . 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 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.
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.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.