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

1. IntroductionThe 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 [V.sub.t] 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, [[??].sub.t] denotes the Brownian motion process under risk neutral measure P, [mu], [alpha], [lambda], r (risk-free interest rate) and a (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

[mathematical expression not reproducible], (1)

the price of American-style put option with underlying stochastic volatility [V.sub.t], current time t, maturity T, and strike price K, where [T.sub.0,T] is the optimal stopping set with the Brownian motion filtration for t [member of] [0, T].

Using transformation of time variable [tau] = 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, [tau]) := f(V, T - [tau]) = f(V, t) and [V.sub.f]([tau]), respectively. Then the valuation of American volatility put option can be formulated as a free-boundary problem

[partial derivative]P(V, [tau])/[partial derivative][tau] = AP(V, [tau])) V > [V.sub.f]([tau]), [tau] > 0, (2)

P(V, [tau]) = K - V, 0 [less than or equal to] V [less than or equal to] [V.sub.f] ([tau]) (3)

[partial derivative]P([V.sub.f]([tau]), [tau])/[partial derivative]V = -1, (4)

[mathematical expression not reproducible], (5)

P(V,0) = max (K - V, 0), (6)

where we can verify that

[mathematical expression not reproducible], (7)

where [1.sub.{x}] 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 [V.sub.f](0) are given in the third column of Table 2, where [V.sup.*] is the unique positive root of the following nonlinear algebraic equation:

(r + [lambda] ln V - [alpha]) V - r K = 0, (8)

It is proved by [17, 18] for GBMP model and [18] for MRGP model, MRSRP model, and MRLP model that the free-boundary functions [V.sub.f]([tau]) are continuously differentiable, strictly decreasing and convex on the interval V [member of] [0, +[infinity]), which are confirmed by the simulation results in [19]. Moreover we note that the free-boundary functions [V.sub.f]([tau]) are bounded on the interval ([V.bar], [bar.V]], where [V.bar] = [V.sub.f]([infinity]) are the early exercise boundaries of the corresponding perpetual American volatility put options, whose values satisfy the following ODEs:

A[P.sup.[infinity]] (V) = 0, V > [V.sub.f] ([infinity]), (9)

[P.sup.[infinity]] (V) = K - V, 0 [less than or equal to] V [less than or equal to] [V.sub.f] ([infinity]) (10)

d[P.sup.[infinity]] ([V.sub.f]([infinity]))/dV = -1, (11)

[mathematical expression not reproducible], (12)

where the differential operators A are listed in the third column of Table 1 with replacement of P by [P.sup.[infinity]].

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 [member of] [0, [V.sub.f]([tau])] for [tau] [member of] [0, [infinity]), then the value of American volatility put is given by P(V, [tau]) = K - V, and we can verify that it is equivalent to the following PDE:

[mathematical expression not reproducible]. (13)

Combining (2) and (13) gives a unified form of PDE:

[mathematical expression not reproducible]. (14)

Recall that [V.bar] [member of] [0,[V.sub.f]([tau])] for all [tau] [member of] [0,+[infinity]); the value of American volatility put option equals K -V on the region V [member of] [0, [V.bar]]. Therefore it suffices to solve the PDE (14) with conditions

P([V.bar], [tau]) = K - [V.bar], (15)

[partial derivative]P ([V.bar], [tau])/[partial derivative]V = -1, (16)

[mathematical expression not reproducible], (17)

P(V,0) = max (K - V, 0). (18)

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

[mathematical expression not reproducible], (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

[L.sub.C] [P (V, [tau])] (z) = zL [P (V, [tau])] (z). (20)

Let [[tau].sub.f](V) be the inverse function of [V.sub.f](t). Since [[tau].sub.f](V) is a strictly decreasing function on the interval [mathematical expression not reproducible] is given by

[mathematical expression not reproducible]. (21)

Applying LCT to PDE (14) with respect to the time variable on the fixed region V [member of] [[V.bar], +[infinity]), we obtain that

[mathematical expression not reproducible], (22)

and LCT to conditions (15)-(17) give that

[??] ([V.bar], z) = K - [V.bar], (23)

[partial derivative] [??] ([V.bar], z)/[partial derivative]V = -1, (24)

[mathematical expression not reproducible]. (25)

The ODEs (22)-(25) cannot be solved directly as the unknown [V.bar] and the function [[tau].sub.f] (V) are involved. In next paragraphs, we derive the value [V.bar] and design an approximation for [[tau].sub.f](V).

Assume that [phi](V) is a solution of (9) and satisfies condition (12). Then [P.sup.[infinity]](V) = C[phi](V) with C being a constant is also the solution of (9) with (12). Using (10) and (11), we have

C[phi]([V.sub.f] ([infinity])) = K - [V.sub.f] ([infinity]), (26)

[mathematical expression not reproducible]. (27)

This gives a nonlinear algebraic equation

[mathematical expression not reproducible], (28)

with unknown [V.sub.f]([infinity]). Solving (28), we get the value [V.bar] = [V.sub.f]([infinity]). In the following, we aim to solve [phi](V) exactly that allows us to avoid using the iterative procedure in [16] for approximating [[tau].sub.f](V) and optimizing NCIMs.

Theorem 1. Under GBMP model, the value of [V.sub.f]([infinity]) is given by

[V.sub.f]([infinity]) = K[gamma]/[gamma] - 1, (29)

[mathematical expression not reproducible]. (30)

Proof. The proof is provided by [20].

Theorem 2. Under MRGP model, the value of [P.sup.[infinity]] (V) is given by

[P.sup.[infinity]](V) = C[psi] [r/2[lambda], 1/2, [([alpha], [lambda]V).sup.2]/[lambda][[sigma].sup.2]], (31)

where [psi](*) 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

[mathematical expression not reproducible]. (32)

Using the variable transformation y = [alpha] - [lambda]V, and the function transformation [??](y) = [P.sup.[infinity]](([alpha] - y)/[lambda]), further inserting

[partial derivative][P.sup.[infinity]]/ [partial derivative]V = -[lambda] [partial derivative][??]/ [partial derivative]t, (33)

[mathematical expression not reproducible], (34)

into (32) yields

[mathematical expression not reproducible], (35)

whose basic solutions are just the degenerate hypergeometric functions [mathematical expression not reproducible] (see, e.g., [21]).

Thus the general solution of (32) can be expressed as

[mathematical expression not reproducible], (36)

where [C.sub.1] and [C.sub.2] are constants to be determined.

Noting that [lim.sub.V[right arrow][infinity]][phi][r/2[lambda], 1/2, [([alpha] - [lambda]V).sup.2]/[lambda][[sigma].sup.2]] = +[infinity], from condition (12), we obtain that [C.sub.1] = 0. Thus the general solution of (32) is

[P.sup.[infinity]] (V) = [C.sub.2][psi][ r/2[lambda], 1/2 [([alpha] - [lambda]V).sup.2]/[lambda][[sigma].sup.2]], (37)

where we complete the proof.

Theorem 3. Under MRSRP model, the value of [P.sup.[infinity]](V) is given by

[mathematical expression not reproducible], (38)

where [W.sub.k,m](*) is the Whittaker function (see, e.g., [22]), and k = 1/4- r/2[lambda], m = i/4, C is a constant to be determined.

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

[mathematical expression not reproducible]. (39)

Using the variable transformation y = [lambda]V/[[sigma].sup.2], and the function transformation [??](y) = [P.sup.[infinity]]([[sigma].sup.2]y/[lambda]), further inserting

[partial derivative][P.sup.[infinity]]/ [partial derivative]V = -[lambda]/[[sigma].sup.2] [partial derivative][??]/ [partial derivative]t, (40)

[mathematical expression not reproducible], (41)

into (39) yields

[mathematical expression not reproducible], (42)

further using the function transformation [??](y) = [y.sup.-1/4] [e.sup.y/2]W(y), and inserting

[mathematical expression not reproducible], (43)

[mathematical expression not reproducible], (44)

into (42) yields

[mathematical expression not reproducible], (45)

where

m = 1/4, k = 1/4 - r/2[lambda], (46)

whose basic solutions are just Whittaker functions [M.sub.k,m](y) and [W.sub.k,m](y) (see, e.g., [22]). Therefore the basic solutions of equation (39) are

[mathematical expression not reproducible], (47)

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

[P.sup.[infinity]] (V) = [C.sub.1][phi](V) + [C.sub.2][psi](V), (48)

where [C.sub.1] and [C.sub.2] are constants to be determined.

Noting that [lim.sub.V[right arrow]][psi](V) = +[infinity], from condition (12), we obtain that [C.sub.2] = 0. Thus the general solution of (39) is

[mathematical expression not reproducible], (49)

where we complete the proof.

Theorem 4. Under MRLP model, the value of [P.sup.[infinity]] (V) is given by

[P.sup.[infinity]] (V) = C[psi] [r/2[lambda], 1/2, ([lambda] ln V - [alpha] + [[sigma].sup.2]/2)2/[lambda][[sigma].sup.2]], (50)

where [psi](*) is the degenerate hypergeometric function, and C is a constant to be determined.

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

[mathematical expression not reproducible]. (51)

Using the variable transformation y = X ln V - [alpha] + [[sigma].sup.2]/2, and the function transformation [mathematical expression not reproducible], further inserting

[mathematical expression not reproducible], (52)

[mathematical expression not reproducible], (53)

into (51) yields

[mathematical expression not reproducible], (54)

whose basic solutions are just the degenerate hypergeometric functions [mathematical expression not reproducible]. Thus the general solution of (51) can be expressed as

[mathematical expression not reproducible], (55)

where [C.sub.1] and [C.sub.2] are constants to be determined.

Noting that [mathematical expression not reproducible], from condition (12), we obtain that [C.sub.1] = 0. Thus the general solution of equation (51) is

[mathematical expression not reproducible], (56)

where we complete the proof.

Recall again that [[tau].sub.f](V) is the strictly decreasing function on the interval ([V.bar], [bar.V]] with [[tau].sub.f]([V.bar]+) = [infinity], [[tau].sub.f]([bar.V]) = 0. Thus it is similar to [16], and we construct an approximation of [[tau].sub.f](V)

[mathematical expression not reproducible], (57)

where [l.sub.k], k = i,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 [V.sub.i] = [V.bar] + i[DELTA]V, i = 0,1,..., N with [V.sub.N] = [V.sub.max], where [V.sub.max] is a sufficiently larger number such that [mathematical expression not reproducible] be the approximation of [??]([V.sub.i], z) at mesh point V = [V.sub.i] for any complex value z, i.e., [mathematical expression not reproducible]; then the ODE (22) can be discretized as, for i = i, 2, ..., N - 1,

[mathematical expression not reproducible], (58)

with boundary conditions

[mathematical expression not reproducible], (59)

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

(A - zI) [??] (z) = g (z), (60)

where the coefficient matrix A is defined as

[mathematical expression not reproducible], (61)

moreover the solution vector [??](z) and the constant vector g(z) are defined as, respectively,

[mathematical expression not reproducible], (62)

[mathematical expression not reproducible], (63)

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

[mathematical expression not reproducible]. (64)

Let a(z) be the first row of [(A - zI).sup.-1]; then a(z) is given by the following linear system:

a (z) (A -zI) = [1, 0, 0, ..., 0]. (65)

Thus the first element in the solution vector is given by

[[??].sub.1] (z) = a (z) g (z; [l.sub.1], [l.sub.2]), (66)

where g(z;[l.sub.1], [l.sub.2]) := g(z), and the approximation of [[tau].sub.f]([V.sub.i]) involved in g(z), for i = i, 2, ..., N - i, is given by (57).

We find the parameters [[l.sub.k]], 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 [z.sub.1], [z.sub.2],..., [z.sub.m] (m [greater than or equal to] 2) and find the values of unknown [[l.sub.k]], k = 1, 2 that minimize the sum of residue

[mathematical expression not reproducible]. (67)

Inserting (66) into (67), we can find the values of parameters [[l.sub.k]], k = i,2 by solving the following discrete optimization problem, for given real positive numbers [z.sub.1], [z.sub.2], ..., [z.sub.m] (m [greater than or equal to] 2),

[mathematical expression not reproducible], (68)

for approximation formula (57). After carrying out the optimization procedure (68), then we can find the approximation of the early exercise boundary [V.sub.f]([tau]) as the inverse of the function [[tau].sub.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,i,..., N,

P ([V.sub.i], [tau]) = [L.sup.-1] [[??]([V.sub.i], z)/z], (69)

and apply the cubic spline interpolation to obtain the option value function P(V, [tau]).

In the following, the hyperbolic contour integral methods for Laplace inversion are studied to restore the option prices. Denote

V = ([V.sub.1], [V.sub.2],..., [V.sub.N-1])', (70)

P ([tau]) := (P ([V.sub.1], [tau]), P([V.sub.2], [tau]),..., P ([V.sub.N-1], [tau]))', (71)

[mathematical expression not reproducible]. (72)

And the term g(z) can be rewritten as two parts

g (z) = [g.sub.1] (z) + [g.sub.2] (z), (73)

With

[mathematical expression not reproducible], (74)

[mathematical expression not reproducible], (75)

where

[mathematical expression not reproducible], (76)

and

[mathematical expression not reproducible], (77)

for i = i,2,..., N - 1.

From (76) and (77), we can see that [mathematical expression not reproducible], which guarantee that the following Laplace inversion formula is feasible (see [16]):

[mathematical expression not reproducible], (78)

where the contours [[GAMMA].sub.1] and [[GAMMA].sub.2] should be chosen as two special Hankel contours as in Figure 1 which enclose all singularities of [mathematical expression not reproducible], respectively. From expression (78), we can see that [mathematical expression not reproducible] lie in the left-half plane of the contour [[GAMMA].sub.1] or equivalently the spectrum [LAMBDA]( A) lies in a sectorial region [[summation].sub.[delta]] that will be proved in the following.

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

[[GAMMA].sub.1] :z([zeta]) = v [i + sin (i[zeta] - [theta])], -[infinity] < [zeta] < [infinity], (79)

where parameters v > 0 and [theta] set the width and the asymptotic angle of the hyperbolic contour, respectively.

Substituting the contour (79) into (78) gives

[mathematical expression not reproducible], (80)

where [[zeta].sub.m] = mh for the trapezoidal rule and L is the number of quadrature nodes. The parameters h, v are given by (see, e.g., [23])

[mathematical expression not reproducible], (81)

which ensure that all singularities of the integrand in (80) lie in a sectorial region

[mathematical expression not reproducible], (82)

where [theta] is a free parameter. The choice of the parameters (81) leads to the following predicted convergence rate:

[mathematical expression not reproducible], (83)

where

R([theta]) = [[pi].sup.2] - 2[pi][theta] -2[pi][delta]/[rho]([theta]). (84)

and [parallel]x[parallel] is some vector norm. According to expressions (81)-(84), after the semiangle [delta] of the sectorial region (82) is given, then we can get the optimal [theta] by maximizing the function R([theta]), and thus the optimal parameters h and v can be computed from the formulas (81).

The parameterized hyperbolic contour [[GAMMA].sub.2] on the right-hand side is selected as

[[GAMMA].sub.2] :z([zeta]) = v [i + sin (i[zeta] + [theta])], -[infinity] < [zeta] < [infinity], (85)

where parameters h and v are given by formulas (81) with setting [delta] = 0 and

[mathematical expression not reproducible]. (86)

To analyze the spectrum of the coefficient matrix (61), we follow the idea in [24] and construct a Toeplitz matrix [bar.A] to replace analyzing (61)

[mathematical expression not reproducible], (87)

where the elements in the matrix [bar.A] are defined as

[mathematical expression not reproducible], (88)

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.bar], [V.sub.max]], respectively. With the constructed Toeplitz matrix, we can analyze the spectrum [LAMBDA](A) and estimate the parameters of hyperbolic contour roughly.

Theorem 5. The numerical range W([bar.A]) of [bar.A] defined by (87) lies in a sectorial region [[summation].sub.[delta]] such that

W ([bar.A]) [subset or equal to] [[summation].sub.[delta]] [equivalent to] {z [member of] C : [absolute value of arg (- z)] [less than or equal to] [delta]}, (89)

with

[delta] = arctan ([absolute value of b]/2[square root of ar]]. (90)

And

[LAMBDA] ([bar.A]) [subset or equal to] W' ([bar.A]) [subset or equal to] [[summation].sub.[delta]]. (91)

Proof. Denote f([eta]) the generating function of [bar.A] and [OMEGA](f([eta])) = [f([eta]) | [eta] [member of] (-[pi],[pi])). We shall use the following two conclusions (see, e.g., [8,16]): the numerical range W([bar.A]) is a subset of the closure of convex hull of [OMEGA](f([eta])), i.e., w([bar.A]) [subset or equal to] [bar.conv([OMEGA](f))]; let A, D [member of] [C.sup.nxn] and suppose that D is positive definite; then the spectrum [LAMBDA](DA) is a subset of the angular numerical range of A, i.e., [LAMBDA](DA) [subset or equal to] W' (A)

The generating function can be written as

[mathematical expression not reproducible], (92)

inserting (88) into (92) yields

[mathematical expression not reproducible], (93)

further using (93), we have

[mathematical expression not reproducible], (94)

which shows that [OMEGA](f([eta])) [subset or equal to] [[summation].sub.[delta]]. Apparently the sectorial region Zs is a closed convex hull; thus we have

W([bar.A]) [subset or equal to] [bar.conv ([OMEGA] (f))] [subset or equal to] [[summation].sub.[delta]]. (95)

Since W'([bar.A]) is also a sectorial region whose angle is just the angular opening of the smallest angular sector that includes W([bar.A]), thus we have

[LAMBDA](A) [subset or equal to] W' ([bar.A]) [subset or equal to] [[summation].sub.[delta]], (96)

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])

RE := [absolute value of ([P.sub.NCIM] - [P.sub.FDM])]/[P.sub.FDM], (97)

In the test, we set [V.sub.max] = 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 [z.sub.k] = k (k = 1, 2, ..., m) with m = 4 (see expression (68)). The Matlab command "fminsearch" is used to search the approximate parameters [l.sub.k], k = 1, 2, with initial values [l.sub.k] = 1, k = 1,2, besides we use the technique in [16] to get the optimal [V.bar] 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 [L.sup.2]-norm, for the given volatility model. The [L.sup.2]-norm is defined as

[mathematical expression not reproducible]. (98)

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 [mu] = 0.4, [sigma] = 0.5, r = 0.05, and K =i. 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.

4.2. Implementations for MRGP Model. In this test, the model parameters are taken as [alpha] = 1.6, [sigma] = 1, [lambda] = 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.

4.3. Implementations for MRSRP Model. In this test, the model parameters are taken as [sigma] = 1.4, [lambda] = 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.

4.4. Implementations for MRLP Model. In this test, the model parameters are taken as [alpha] = 0.6, [lambda] = 0.i, [sigma] = 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.

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.

https://doi.org/10.1155/2018/1838521

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

[1] M. Brenner and D. Galai, "New financial instruments for hedging changes in volatility," Financial Analysts Journal, vol. 45, no. 4, pp. 61-65, 1989.

[2] R. E. Whaley, "Derivatives on market volatility: hedging tools long overdue," The Journal of Derivatives, vol. 1, no. 1, pp. 71-84, 1993.

[3] A. Grunbichler and F. A. Longstaff, "Valuing futures and options on volatility," Journal of Banking & Finance, vol. 20, no. 6, pp. 985-1001, 1996.

[4] A. Sepp, "VIX option pricing in a jump-diffusion model," Risk, pp. 84-89, 2008.

[5] P. Carr and R. Lee, "Volatility Derivatives," Annual Review of Financial Economics, vol. 1, pp. 1-20, 2009.

[6] 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.

[7] 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.

[8] 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.

[9] Y.-N. Lin and C.-H. Chang, "VIX option pricing," Journal of Futures Markets, vol. 29, no. 6, pp. 523-543, 2009.

[10] 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.

[11] 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.

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

[13] 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.

[14] Z. Wang and R. T. Daigler, "The performance of VIX option pricing models: empirical evidence beyond simulation," Journal ofFutures Markets, vol. 31, no. 3, pp. 251-281, 2011.

[15] J. Detemple and C. Osakwe, "The valuation of volatility options," European Finance Review, vol. 4, pp. 21-50, 2000.

[16] 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.

[17] 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.

[18] H. K. Liu, "Convexity of the exercise boundary of the American-style put," http://cn.arxiv.org/pdf/1304.5337v3.

[19] J. Ma, W. Li, and X. Han, "Stochastic lattice models for valuation of volatility options," Economic Modelling, vol. 47, pp. 93-104, 2015.

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

[21] 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.

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

[23] J. A. Weideman, "Parabolic and hyperbolic contours for computing the Bromwich integral," Mathematics of Computation, vol. 76, no. 259, pp. 1341-1356, 2007.

[24] 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.

[25] 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.

[26] 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.

Yong Chen (iD) (1) and Jianjun Ma (2)

(1) School of Economic Mathematics, Southwestern University of Finance and Economics, Chengdu, 611130, China

(2) 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

Caption: Figure 1: An example of Hankel contours.

Caption: Figure 2: Early exercise boundaries (left) under GBMP model with T/2000 for FDM and L =i4 for NCIM, NCIM errors in log-scale versus L (right).

Caption: Figure 3: Early exercise boundaries (left) under MRGP model with T/2000 for FDM and L = 6 for NCIM, NCIM errors in log-scale versus L (right).

Caption: Figure 4: Early exercise boundaries (left) under MRSRP model with T/2000 for FDM and L =i0 for NCIM, NCIM errors in log-scale versus L (right).

Caption: Figure 5: Early exercise boundaries (left) under MRLP model with T/2000 for FDM and L =i4 for NCIM, NCIM errors in log-scale versus L (right).

Caption: Table 1: Volatility models and infinitesimal generators.

Caption: Table 3: Difference operators.

Caption: Table 4: Coefficients.

Table 2: The values of A(K - V) and [V.sub.f](0). Model si(K--V) [V.sub.f] (0) GBMP (r - 2[mu] - K [[sigma].sup.2])V - rK MRGP (r + [lambda]) V - min (K, [alpha] + rK/ [alpha]- rK [lambda] + r) MRSRP (r + 2 [lambda])V - min (K, [sigma].sup.2] + rK/ [[sigma].sup.2] -rK 2[lambda] + r) MRLP (r + [lambda] ln V - min (K, [V.sup.*]) [alpha]) V - rK Table 5: Comparison results (option values, relative errors and CPU time). FDM [V.sub.0] T/250 T/500 T/1000 T/2000 value 0.5 0.50000 0.50000 0.50000 0.50000 1 0.19484 0.19489 0.19491 0.19493 1.5 0.10127 0.10132 0.10135 0.10136 CPU times (s) 315 453 732 1335 NCIM [V.sub.0] L=14 RE (%) value 0.5 0.50000 0.00000 1 0.19441 0.26676 1.5 0.10113 0.22691 CPU times (s) 186 -- Table 6: Comparison results (option values, relative errors and CPU time). FDM T/250 T/500 T/1000 T/2000 value 0.5 0.69600 0.69624 0.69636 0.69642 1 0.58613 0.58649 0.58668 0.58677 1.5 0.52618 0.52663 0.52685 0.52696 CPU times (s) 613 785 1141 1808 NCIM L=6 value RE (%) 0.5 0.69152 0.70217 1 0.58206 0.80270 1.5 0.52280 0.78755 CPU times (s) 87 -- Table 7: Comparison results (option values, relative errors and CPU time). FDM n T/250 T/500 T/1000 T/2000 value 4 2.46906 2.46997 2.47044 2.47068 5 2.15034 2.15138 2.15189 2.15216 6 1.87964 1.88074 1.88129 1.88158 CPU times (s) 551 716 1045 1729 NCIM n L=10 RE (%) value 4 2.45192 0.75931 5 2.13482 0.80570 6 1.86568 0.84503 CPU times (s) 95 -- Table 8: Comparison results (option values, relative errors and CPU time). FDM [V.sub.0] T/250 T/500 T/1000 T/2000 value 4 1.94291 1.94363 1.94399 1.94418 5 1.63948 1.64031 1.64073 1.64095 6 1.40529 1.40617 1.40662 1.40684 CPU times (s) 431 582 1519 NCIM [V.sub.0] L=14 RE (%) value 4 1.93598 0.42177 5 1.62844 0.76236 6 1.39302 0.98234 CPU times (s) 128 --

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Chen, Yong; Ma, Jianjun |

Publication: | Discrete Dynamics in Nature and Society |

Date: | Jan 1, 2018 |

Words: | 5870 |

Previous Article: | Multiperiod Telser's Safety-First Portfolio Selection with Regime Switching. |

Next Article: | Research on Cascading Failure Model of Urban Regional Traffic Network under Random Attacks. |

Topics: |