Discrete Dynamics in Nature and Society
Volume 2018 (2018), Article ID 1838521, 12 pages
https://doi.org/10.1155/2018/1838521

Research Article

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

Yong Chen¹ and Jianjun Ma²

¹School of Economic Mathematics, Southwestern University of Finance and Economics, Chengdu, 611130, China
²International Business School, Sichuan International Studies University, Chongqing, 400031, China

Correspondence should be addressed to Yong Chen; chenyong168999@163.com

Received 15 May 2018; Revised 15 October 2018; Accepted 16 October 2018; Published 4 November 2018

Academic Editor: Yong Zhou

Copyright © 2018 Yong Chen and Jianjun Ma. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

Table 1: Volatility models and infinitesimal generators.

Model
Volatility process

GBMP

MRGP

MRSRP

MRLP

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 by

the 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 problem

where we can verify that

where 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:

Table 2: The values of and .

Model

GBMP

MRGP

MRSRP

MRLP

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 as

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 , , then the LCT of 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., [21]), and is a constant to be determined.

Proof. Under MRGP model, the ODE (9) can be rewritten as

Using the variable transformation , and the function transformation , further inserting

into (32) yields

whose basic solutions are just the degenerate hypergeometric functions and (see, e.g., [21]).
Thus the general solution of (32) can be expressed as

where and are constants to be determined.
Noting that , from condition (12), we obtain that . Thus the general solution of (32) is

where we complete the proof.

Theorem 3. Under MRSRP model, the value of is given by

where 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 as

Using the variable transformation , and the function transformation , further inserting

into (39) yields

further using the function transformation , and inserting

into (42) yields

where

whose basic solutions are just Whittaker functions and (see, e.g., [22]). Therefore the basic solutions of equation (39) are

and the general solution of (39) can be expressed as

where and are constants to be determined.
Noting that , from condition (12), we obtain that . Thus the general solution of (39) is

where we complete the proof.

Theorem 4. Under MRLP model, the value of is given by

where is the degenerate hypergeometric function, and is a constant to be determined.

Proof. Under MRLP model, the ODE (9) can be rewritten as

Using the variable transformation , and the function transformation , further inserting

into (51) yields

whose basic solutions are just the degenerate hypergeometric functions and . Thus the general solution of (51) can be expressed as

where and are constants to be determined.
Noting that , from condition (12), we obtain that . Thus the general solution of equation (51) is

where 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.

Table 3: Difference operators.

Model

GBMP

MRGP

MRSRP

MRLP

Using (59), the linear system (58) can be expressed as in the matrix form

where the coefficient matrix is defined as

moreover the solution vector and the constant vector are defined as, respectively,

where elements in the matrix are listed in Table 4, and

Table 4: Coefficients.

Model

GBMP

MRGP

MRSRP

MRLP

Let be the first row of ; then is given by the following linear system:

Thus the first element in the solution vector is given by

where , 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. Denote

And the term can be rewritten as two parts

With

where

and

for .

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.

Figure 1: An example of Hankel contours.

Following [16, 23], the parameterized hyperbolic contour on the left-hand side is selected as

where parameters 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 on the right-hand side is selected as

where 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 as

where 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 that

with

And

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 as

inserting (88) into (92) yields

further using (93), we have

which shows that . Apparently the sectorial region is a closed convex hull; thus we have

Since is also a sectorial region whose angle is just the angular opening of the smallest angular sector that includes , thus we have

where 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 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 , , , 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 .

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

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 with for FDM and for NCIM, NCIM errors in log-scale versus (right).

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.

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

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 with for FDM and for NCIM, NCIM errors in log-scale versus (right).

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.

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

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 with for FDM and for NCIM, NCIM errors in log-scale versus (right).

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.

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

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 with for FDM and for NCIM, NCIM errors in log-scale versus (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.

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.

References

M. Brenner and D. Galai, “New financial instruments for hedging changes in volatility,” Financial Analysts Journal, vol. 45, no. 4, pp. 61–65, 1989.

R. E. Whaley, “Derivatives on market volatility: hedging tools long overdue,” The Journal of Derivatives, vol. 1, no. 1, pp. 71–84, 1993.

A. Grünbichler and F. A. Longstaff, “Valuing futures and options on volatility,” Journal of Banking & Finance, vol. 20, no. 6, pp. 985–1001, 1996.

A. Sepp, “VIX option pricing in a jump-diffusion model,” Risk, pp. 84–89, 2008.

P. Carr and R. Lee, “Volatility Derivatives,” Annual Review of Financial Economics, vol. 1, pp. 1–20, 2009.

S.-P. Zhu and G.-H. Lian, “An analytical formula for VIX futures and its applications,” Journal of Futures Markets, vol. 32, no. 2, pp. 166–190, 2012.

D. Psychoyios and G. Skiadopoulos, “Volatility options: hedging effectiveness, pricing, and model error,” Journal of Futures Markets, vol. 26, no. 1, pp. 1–31, 2006.

J. Ma and Z. Zhou, “Convergence analysis of iterative Laplace transform methods for the coupled PDEs from regime-switching option pricing,” Journal of Scientific Computing, vol. 75, no. 3, pp. 1656–1674, 2018.

Y.-N. Lin and C.-H. Chang, “VIX option pricing,” Journal of Futures Markets, vol. 29, no. 6, pp. 523–543, 2009.

D. Psychoyios, G. Dotsis, and R. N. Markellos, “A jump diffusion model for VIX volatility options and futures,” Review of Quantitative Finance and Accounting, vol. 35, no. 3, pp. 245–269, 2010.

Y.-N. Lin and C.-H. Chang, “Consistent modeling of S&P 500 and VIX derivatives,” Journal of Economic Dynamics and Control (JEDC), vol. 34, no. 11, pp. 2302–2319, 2010.

Y.-N. Lin and C.-H. Chang, “Rejoinder to a remark on Lin and Chang's paper 'Consistent modeling of S&amp;P 500 and VIX derivatives',” Journal of Economic Dynamics and Control (JEDC), vol. 36, no. 5, pp. 716–718, 2012.

G. D. Graziano and L. Torricelli, “Target volatility option pricing,” International Journal of Theoretical and Applied Fiance, vol. 15, pp. 1250005-1–1250005-17, 2012.

Z. Wang and R. T. Daigler, “The performance of VIX option pricing models: empirical evidence beyond simulation,” Journal of Futures Markets, vol. 31, no. 3, pp. 251–281, 2011.

J. Detemple and C. Osakwe, “The valuation of volatility options,” European Finance Review, vol. 4, pp. 21–50, 2000.

Z. Q. Zhou, J. T. Ma, and H.-W. Sun, “Fast Laplace transform methods for free-boundary problems of fractional diffusion equations,” Journal of Scientific Computing, vol. 74, no. 1, pp. 49–69, 2018.

X. Chen and J. Chadam, “A mathematical analysis of the optimal exercise boundary for American put options,” SIAM Journal on Mathematical Analysis, vol. 38, no. 5, pp. 1613–1641, 2006/07.

H. K. Liu, “Convexity of the exercise boundary of the American-style put,” http://cn.arxiv.org/pdf/1304.5337v3.

J. Ma, W. Li, and X. Han, “Stochastic lattice models for valuation of volatility options,” Economic Modelling, vol. 47, pp. 93–104, 2015.

L. Jiang, Mathematical Modeling and Methods of Option Pricing, World Scientific Publishing, Singapore, 2005.

A. D. Polyanin and V. F. Zaitsev, Handbook of Exact Solutions for Ordinary Differential Equations, Chapman and Hall/CRC Press, Boca Raton, Fla, USA, 2003.

M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover Publications, New York, NY, USA, 1965.

J. A. Weideman, “Parabolic and hyperbolic contours for computing the Bromwich integral,” Mathematics of Computation, vol. 76, no. 259, pp. 1341–1356, 2007.

H.-K. Pang and H.-W. Sun, “Fast numerical contour integral method for fractional diffusion equations,” Journal of Scientific Computing, vol. 66, no. 1, pp. 41–66, 2016.

A. Chockalingam and K. Muthuraman, “An approximate moving boundary method for American option pricing,” European Journal of Operational Research, vol. 240, no. 2, pp. 431–438, 2015.

Z. Zhou and X. Gao, “Laplace transform methods for a free boundary problem of time-fractional partial differential equation system,” Discrete Dynamics in Nature and Society, vol. 2017, Article ID 6917828, 9 pages, 2017.