Tsukuba Economics Working Papers
No. 2009-008

An Approximation of European Option Prices
under General Diffusion Processes
byHideyuki Takamizawa
April 2009This version: February 2011

UNIVERSITY OF TSUKUBA
Department of Economics
1-1-1 Tennodai
Tsukuba, Ibaraki 305-8571
JAPANAn Approximation of European Option Prices
under General Diﬀusion Processes

Hideyuki Takamizawa
Graduate School of Humanities and Social Sciences,
University of Tsukuba
Tsukuba Ibaraki 305-8571, JAPANE-mail: takamiza[atmark]social.tsukuba.ac.jp

First Draft: Jun 2008 / This Draft: February 2011

Abstract
This study proposes an approximation of European option prices under arbitrary diﬀusion processes of the spot price. The key is to approximate the characteristic functionof the log spot price process as the solution to ordinary diﬀerential equations. The option price is then obtained by the inverse Fourier transform. Numerical experiments,using a model that has the constant elasticity of volatility specification in both the
spot price and volatility processes, confirm reasonable accuracy of the approximation,
except when the volatility process exhibits high variation.

Keywords: Approximation, Conditional moment, European option, Fourier transform,
Stochastic volatility.

11

Introduction

The Fourier Transform (FT) approach to the pricing of options, originally introduced
by Heston (1993), allows for the extension of a class of models under which an eﬃcientcomputation of option prices is possible. This is because the FT approach works when
the characteristic function of a (log) spot price process has a closed form even though
the density function does not. This weaker requirement leads to the development of more
realistic models of the spot price process. Examples include aﬃne jump-diﬀusion models[Bakshi et al. (1997), Bates (1996, 2000), Duﬃe et al. (2000), Pan (2002), and Scott(1997)]; variance-gamma models [Madan and Seneta (1990), and Madan et al. (1998)];
the normal inverse Gaussian model [Barndorﬀ-Nielsen (1998)]; the finite moment log-stablemodel [Carr and Wu (2003a)]; the CGMY model [Carr et al. (2002)]; and discrete-time
GARCH models [Heston and Nandi (2000)].
Nevertheless, it seems that other realistic models are worth pursuing. One notable
example is a model that accommodates both stochastic volatility (SV) and constant elasticity of volatility (CEV), as proposed by Jones (2003), and Melino and Turnbull (1990)
in the option-pricing literature, and by Andersen and Lund (1997), Brenner et al. (1996),
and Gallant and Tauchen (1998) in the time-series literature. This volatility specification
can generate skewness and leptokurtosis in the conditional distribution of the spot price,which are fundamental for many securities in capturing not only time-series properties ofthe spot price but also cross-sectional properties of the option prices observed as impliedvolatility smiles or smirks. Furthermore, in pricing options written on an asset that issensitive to interest rates, the stochastic risk-free rate needs to be incorporated, possiblywith the instantaneous correlation between the spot price and the risk-free rate. In theabove cases, however, even the characteristic function of the spot price process is unavailable in closed form, which makes it diﬃcult to achieve an eﬃcient computation of option
prices based on these models.
This study proposes an approximation of European option prices under general diffusion processes of the spot price. The key is to approximate the characteristic functionof the log price process using a method originally proposed by Shoji (2002) and recently
applied to bond pricing by Takamizawa and Shoji (2009). The method approximates conditional moments of diﬀusion processes as the solution to ordinary diﬀerential equations.
Since the characteristic function is basically given by the expectation, this method can beapplicable. Once an analytical expression of the characteristic function is available, theoption price is obtained by Fourier inversion.
Using an SV model that has the CEV specification in both the spot price and volatility
processes, the accuracy of the approximation is examined by numerical experiments, where
benchmark option prices are computed by the Monte Carlo (MC) method. The numericalresults are summarized as follows. (i) The third-order approximation generally achieves2high accuracy, except when the volatility process exhibits high variation. (ii) When themean reversion of the volatility process is relatively fast, the second-order approximation
performs reasonably well. Actually, the higher accuracy of the third-order approximationis achieved at the cost of computational complexity. This study also explores the possibilityof reducing this complexity without much reducing the accuracy.
Section 2 explains the method of approximating conditional moments, whose application within the FT approach is presented in Section 3. Section 4 examines the accuracyof the approximation. Section 5 identifies conditional third moments of the underlying
processes that significantly aﬀect option prices, which is aimed at reducing computationalburden. Section 6 provides concluding remarks.

2

An Approximation Formula of Conditional Moments

A method of approximating conditional moments, originally proposed by Shoji (2002), is
widely applicable to the computation of up to conditional n-th moments, if they exist,
for a d-dimensional diﬀusion process. An important feature of the method is that all the
moments considered are computed simultaneously as the solution to a system of ordinary
diﬀerential equations. Here, a more specific explanation of the method is presented, taking
the application to the pricing of options into consideration. First is the case with (n, d) =(1, 2), i.e., the conditional first moments of a two-dimensional process, followed by thecases with up to (n, d) = (3, 3), i.e., the conditional third moments of a three-dimensional
process. As seen below, n can be considered as the order of approximation.

2.1

Conditional first moments of a two-dimensional process

Let Xt = (x1,t x2,t )� be a two-dimensional stochastic process, which evolves according to
the following stochastic diﬀerential equation (SDE):
dx1,t = f1 (Xt )dt + ξ1 (Xt )� dWt ,(1)

dx2,t = f2 (Xt )dt + ξ2 (Xt )� dWt ,(2)

where Wt is the two-dimensional Brownian motion, and the drift and diﬀusion functions,fi and ξi (i = 1, 2), satisfy certain technical conditions for the solutions to Eqs.(1) and(2) to exist for an arbitrary X0 . It is also assumed that fi and gij = ξi� ξj (i, j = 1, 2)
are appropriately smooth with respect to Xt . Specifically, according to Shoji (2002), they
belong to C n+1 -class in computing up to conditional n-th moments by the method, which
is needed to prove the convergence of approximate moments to the true ones.
Let Ψs,t be a vector consisting of the first moments of (an increment of) Xt conditioned
on time s ≤ t:

Ψs,t = Es (x1,t − x1,s3

x2,t − x2,s )� ,

where Es stands for the conditional expectation. By integrating Eq.(1) and taking theconditional expectation,Es [x1,t − x1,s ] = Es

�� t
s

�

f1 (Xu )du .

(3)

By applying the Taylor expansion to f1 (Xu ) around Xs up to the first order and substituting this into Eq.(3),Es [x1,t − x1,s ] = f1 (Xs )(t − s)+Es�� t �
s

(1,0)
(0,1)
f1 (Xs )(x1,u − x1,s ) + f1 (Xs )(x2,u − x2,s )��

du + R1 ,

(4)

i+j

where f (i,j) = ∂ i fj , and R1 is a residual term. Eq.(4) can be expressed in a vector form
∂x1 ∂x2

as

�

(1,0)

Es [x1,t − x1,s ] = f1 (Xs )(t − s) + f1
Similarly,

�

(1,0)

Es [x2,t − x2,s ] = f2 (Xs )(t − s) + f2

(0,1)

(Xs )

f1

(Xs )

f2

(0,1)

(Xs )

�� t

Ψs,u du + R1 .

(5)

(Xs )

�� t

Ψs,u du + R2 .

(6)

s

s

Expressing Eqs.(5) and (6) together in a vector form leads toΨs,t = A(Xs )� tsΨs,u du + b(Xs )(t − s) + R ,

(7)

where


A(Xs ) = (1,0)(Xs ) f1

(1,0)

(Xs ) f2

f1
f2

(0,1)

(Xs )

(0,1)

(Xs )



 ,

and R = (R1 R2 )� . Eq.(7) can be developed toΨs,t =� ts

b(Xs ) = 

f1 (Xs )f2 (Xs )eA(Xs )(t−u) b(Xs )du + R̂ ,



 ,

(8)

where R̂ is a residual vector. If, in addition, A is invertible,
Ψs,t = A−1 (Xs ){eA(Xs )(t−s) − I}b(Xs ) + R̂ .

(9)

As seen below, Eqs.(7)–(9) hold for any (n, d) with modification to A(Xs ) and b(Xs ).
Omitting the residual vector leads to the approximation. According to Shoji (2002), both
R and R̂ have order of O((t−s)(n+3)/2 ) when up to conditional n-th moments are computed.
Thus, n can be considered as the order of approximation. In this particular case, it isO((t − s)2 ) for n = 1. Indeed, the conditional moments computed by the method are

more accurate than those computed by the conventional Euler method. To illustrate this,
suppose a one-dimensional lognormal process Xt , which evolves according to the following
1 2

SDE: dXt = aXt dt + bXt dWt . We know that Es [Xt − Xs ] = Xs {e(a− 2 b )(t−s) − 1}. On the4other hand, it is approximated by the Euler method as Es [Xt − Xs ] = aXs (t − s), and by
the proposed method as Es [Xt − Xs ] = Xs {ea(t−s) − 1} by substituting A(Xs ) = a andb(Xs ) = aXs into Eq.(9) (without the residual term). Obviously, the latter is closer to thetrue value. When b is large, however, the approximation error of the proposed method
may not be negligible. A more important message from this simple example, therefore,
is that even though the focus is on the computation of lower-order moments, it is more
appropriate to include information on higher-order moments.

2.2

Up to conditional second moments of a two-dimensional process

The moment vector, Ψs,t , consists of up to the conditional second moments. It is here afive-dimensional vector as�Ψs,t = Es x1,t − x1,s

x2,t − x2,s

(x1,t − x1,s )2(x2,t − x2,s )2��(x1,t − x1,s )(x2,t − x2,s ) .

By applying the Taylor expansion to fi (Xu ) (i = 1, 2) around Xs up to the second
order,fi (Xu ) = fi (Xs )
(1,0)

+fi

(0,1)

(Xs )(x1,u − x1,s ) + fi

1 (2,0)(Xs )(x2,u − x2,s ) + fi (Xs )(x1,u − x1,s )22

1 (0,2)(1,1)+ fi (Xs )(x2,u − x2,s )2 + fi (Xs )(x1,u − x1,s )(x2,u − x2,s ) + ei ,2(10)where ei is a residual term. By substituting Eq.(10) into the right-hand side of Es [xi,t −�

xi,s ] = Es [ st fi (Xu )du],Es [xi,t − xi,s ] = fi (t − s)�
1 (2,0)
(1,0)
(0,1)
+ fi
fi
f
2 i

1 (0,2)
f
2 i

(1,1)
fi

�� t
s

Ψs,u du + Ri ,(11)

where Xs is abbreviated for notational convenience.
Next, by applying the Ito formula to (x1,t − x1,s )2 and taking the conditional expecta-

tion,

2

Es [(x1,t − x1,s ) ] = Es

�� t
s

�

{2f1 (Xu )(x1,u − x1,s ) + g11 (Xu )}du ,(12)

where g11 = ξ1� ξ1 . By applying the Taylor expansion to f1 (Xu ) and g11 (Xu ) around Xs
up to the first and second orders, respectively, the integrand of Eq.(12) becomes
2f1 (Xu )(x1,u − x1,s ) + g11 (Xu )

(1,0)

(0,1)

= g11 (Xs ) + {2f1 (Xs ) + g11 (Xs )}(x1,u − x1,s ) + g11 (Xs )(x2,u − x2,s )
1 (2,0)1 (0,2)(1,0)+{2f1 (Xs ) + g11 (Xs )}(x1,u − x1,s )2 + g11 (Xs )(x2,u − x2,s )2
2
2
(0,1)
(1,1)
+{2f1 (Xs ) + g11 (Xs )}(x1,u − x1,s )(x2,u − x2,s ) + e11 ,

5

(13)

where g (i,j) is defined analogously with f (i,j) and e11 is a residual term. By substituting
Eq.(13) into Eq.(12),Es [(x1,t − x1,s )2 ] = g11 (t − s)
+
×

�

(1,0)
2f1 + g11

� t

(0,1)
g11

1 (2,0)
(1,0)
2f1
+ g11
2

1 (0,2)
g
2 11

(0,1)
(1,1)
2f1
+ g11

�

Ψs,u du + R11 .

s

(14)

A similar manipulation is applied to Es [(x2,t − x2,s )2 ] and Es [(x1,t − x1,s )(x2,t − x2,s )].

Expressing the resulting equations together in a vector form leads to Eq.(7), where


A =

(1,0)

(0,1)

f1


(1,0)

f2


 2f + g (1,0)
1

11

(1,0)

g22

(1,0)

f2 + g12b� = (f1

f2

g11

f1

(0,1)

f2

(0,1)

g11

(0,1)

2f2 + g22

(0,1)

f1 + g12g22

1 (2,0)2 f1
1 (2,0)
2 f2
(1,0)
(2,0)
2f1
+ 12 g11
1 (2,0)
2 g22
(1,0)
(2,0)
f2
+ 12 g12

1 (0,2)2 f1
1 (0,2)2 f2
1 (0,2)
2 g11
(0,1)
(0,2)
2f2
+ 12 g22
(0,1)
(0,2)
f1
+ 12 g12



(1,1)

f1

(1,1)

f2
(0,1)

+ g11

(1,0)

+ g22

2f1
2f2
(1,0)

f1

(1,1)

(0,1)

+ f2

(1,1)
(1,1)

+ g12

g12 ) .

Note that the residual terms, Ri and Rij (i, j = 1, 2), contain (the conditional expectation of) derivatives of fi higher than the first order and derivatives of gij higher than thesecond order. Then, if fi and gij are linear and quadratic, respectively, there is no residual
term. In other words, the conditional moments computed by the method are exact. Evenin this case, it may be beneficial to use this method when the derivation of the closed form
conditional moments is demanding.

2.3

Up to conditional third moments of a three-dimensional process

Since the derivation is basically the same, only the resulting A(Xs ) and b(Xs ) are presented
with a brief explanation of notations. Let Xt = (x1,t x2,t x3,t )� be a three-dimensional
stochastic process. The SDE for x3,t is specified analogously, where the drift and diﬀusion
functions are denoted as f3 (Xt ) and ξ3 (Xt ), respectively, with gij = ξi� ξj (i, j = 1, 2, 3).
Let ψpqr,s,t = Es [(x1,t − x1,s )p (x2,t − x2,s )q (x3,t − x3,s )r ], and a moment vector consistingof up to the conditional second moments can be expressed asΨn=2s,t = (ψ100,s,t

ψ010,s,t

ψ001,s,t

ψ200,s,t

ψ020,s,t

ψ002,s,t

ψ110,s,tψ101,s,t

ψ011,s,t )� ,(15)

where the order of approximation, n = 2, is indicated in the upper suﬃx here. A moment
n=2
vector for n = 3, Ψn=3
s,t , is then obtained by augmenting Ψs,t with a vector consisting of

the conditional third moments:
n=2 �
Ψn=3
s,t = (Ψs,tψ210,s,t

ψ300,s,tψ201,s,t

ψ030,s,t

ψ003,s,t

ψ120,s,t6

ψ111,s,tψ102,s,t

ψ021,s,t

ψ012,s,t )� , (16)





 ,





Table 1 presents the elements of A(Xs ) and b(Xs ) for Ψn=3
s,t . The corresponding matrix
and vector for Ψn=2
s,t are the 9 × 9 upper-left submatrix of A(Xs ) and the 9 × 1 subvectorof b(Xs ) for Ψn=3
s,t , respectively. As seen later, when a two-factor model is considered inwhich the spot price and volatility are stochastic, Table 1 is used for the reference. Whichn to choose depends on the extent to which accuracy one requires. This is examined bynumerical experiments in Section 4. When a three-factor model is considered in which
the risk-free rate is also stochastic or the volatility is driven by two stochastic factors,
the conditional moments of a four-dimensional process are required. The derivation ofA(Xs ) and b(Xs ) for a four-dimensional process is a straightforward extension of thelower-dimensional processes.

3

Application of the Approximation Formula

To clarify how the proposed method is utilized within the FT approach to the pricing ofoptions, the outline of this approach is first provided, followed by the application of theapproximation formula.

3.1

Fourier transform of the option value

Let St be the spot price at time t. The price of a European put option at time 0 with amaturity date t and a strike price K is given by�t

P0 (K, t) = E0 [e− 0 ru du (K − St )1{St <K} ] .

(17)

The price of a call option is immediately obtained by the put-call parity. St and rt may
depend on some state variables, however, they are suppressed at present.
Let st and k be logarithm of St and K, respectively: st = ln St and k = ln K. Define
�t

(18)

P0 (k, t) = ek G0 (k; t, 0) − G0 (k; t, 1) .

(19)

G0 (k; t, a) = E0 [e− 0 ru du east 1{st <k} ] ,and the put option price is given byThus, the pricing of European options is equivalent to the evaluation of G0 (k; t, a).
Following Duﬃe et al. (2000), the Fourier transform of G0 (k; t, a) is given by
H0 (θ; t, a) =

�

e

iθk

dG0 (k; t, a) = E0 [eDefine

�t

−

�t
0

ru du (a+iθ)stφ0 (α; t) = E0 [e− 0 ru du eαst ] ,

7

e

].

(20)

(21)

and H0 (θ; t, a) = φ0 (a + iθ; t). The option value is then recovered by the inverse Fouriertransform as
1
1
G0 (k; t, a) = φ0 (a; t) −
2
π

� ∞
0

1
Im[e−iθk φ0 (a + iθ; t)]dθ ,θ(22)where Im(x) = β for x = α + βi. The feasibility of this FT approach depends on theavailability of an analytical expression of φ0 (α; t).

3.2

Approximation to φ0 (α; t)Suppose the risk-neutral dynamics of the log spot price, st , is given by the following SDE:
dst = µ(st , Yt )dt + σ(st , Yt )dWt ,(23)

where Wt is the Brownian motion under the risk-neutral measure, and the drift anddiﬀusion functions possibly depend on a state vector, Yt . By the absence of arbitrage,µ(st , Yt ) = r(Yt ) − 12 σ 2 (st , Yt ) holds, where the instantaneous risk-free rate, r(Yt ), mayalso depend on (some elements of) Yt .
Now, define

�t

(24)

dZt = µz (st , Yt , Zt )dt + σz (st , Yt , Zt )dWt ,(25)

Zt = e− 0 ru du eαst ,and φ0 (α; t) is the conditional first moment of Zt , i.e., φ0 (α; t) = E0 [Zt ]. By applying theIto formula to Zt ,

where
�

�

1
−r(Yt ) + αµ(st , Yt ) + α2 σ 2 (st , Yt ) Zt ,2
σz (st , Yt , Zt ) = ασ(st , Yt )Zt .

µz (st , Yt , Zt ) =

(26)
(27)

Based on the SDE for Zt together with those for st and Yt , the elements of A and b inEq.(7) are determined.

4

Accuracy of the Approximation

4.1

Setup

The MC method is employed to obtain benchmark prices of European put options, where
the number of repetitions is set to one million with antithetic variates. In generating
sample paths, the SDEs under the risk-neutral measure are discretized by the Euler methodwith a step size of 1/1, 000. The accuracy is then measured by the relative pricing error,
pap /pmc − 1, where pap and pmc denote the approximate and MC prices, respectively.

8The following SV model is considered, which has the CEV specification in both the
spot price and volatility processes, and thus is labeled SV-CEV:
�

�

�
dSt
St γ1 −1
= rt dt + Vt
dW1,t ,
St
S0
dVt = κ2 (θ2 − Vt )dt + σ2 Vtγ2 dW2,t ,

(28)
(29)

where Wi,t (i = 1, 2) are the Brownian motions under the risk-neutral measure withEt [dW1,t dW2,t ] = ρ12 dt. The SV-CEV model nests both the Heston (1993) and Hull andWhite (HW) (1987) models: the former is obtained by γ1 = 1 and γ2 = 0.5, whereas thelatter is by γ1 = γ2 = 1 and θ2 = 0.
Note that the diﬀusion term of the instantaneous rate of return on St depends on
the power of St /S0 , i.e., the normalized price. This specification is convenient in thefollowing reason. For small t, which corresponds to an option maturity, St /S0 is likely
to be around unity, and so is (St /S0 )γ1 −1 with γ1 �= 1. Then, Vt can be interpreted as(nearly) the instantaneous variance of the return, as is the case of γ1 = 1. This helps
select appropriate parameter values of the SDE for Vt together with its initial value, V0 .
In applying the approximation formula, the SDE for the log spot price, st = ln St , isactually required. By the Ito formula,
�

�

�
1
dst = rt − Vt e2(γ1 −1)(st −s0 ) dt + Vt e(γ1 −1)(st −s0 ) dW1,t .
2(30)

As for the instantaneous risk-free rate, both cases of deterministic and stochastic rates
are considered. Since the accuracy results are similar between the two cases, only theformer is reported while the latter is available upon request.
The input values for option prices are as follows: S0 = 1; K = 0.9, 1.0, 1.1; t = 0.5(half a year); V0 = 0.152 ; and rt = 0 (constant). These values are basically determinedfollowing HW (1987, Table II).
Before reporting the numerical results, it is noted that the accuracy of the proposed
method improves by the logarithmic transformation of the volatility process, i.e., vt = ln Vt .
This can be explained as follows. Looking at Table 1, where (st , Vt , Zt ) is substituted for(x1,t , x2,t , x3,t ), E0 [Zt −Z0 ] is in the third row of the moment vector, Ψ0,t . Then, by Eq.(7),the third row of the matrix A is the most influential for this calculation. The third row of
A consists of the partial derivatives of the drift function of Zt . From Eqs.(26) and (30),the drift function of Zt is given by��1f3 (x1,t , x2,t , x3,t ) = µz (st , Vt , Zt ) = (α − 1) rt + αVt e2(γ1 −1)(st −s0 ) Zt .
2

(31)

Since µz is linear in Vt , the second partial derivative with respect to Vt is zero, i.e.,(0,2,0)f3= 0. This is the (3, 5)-element of A, which multiplies the fifth element of Ψ0,u ,(0,2,0)i.e., E0 [(Vu − V0 )2 ], on the right-hand side of Eq.(7). Then, f3

9= 0 indicates thatthe conditional second moment of (an increment of) Vt is not eﬀectively used for thecomputation of E0 [Zt − Z0 ].

After the logarithmic transformation of Vt , on the other hand, this problem is avoided.

After the transformation, µz becomes��1f3 (x1,t , x2,t , x3,t ) = µz (st , vt , Zt ) = (α − 1) rt + αe2(γ1 −1)(st −s0 )+vt Zt ,2(0,2,0)which is nonlinear in vt . Thus, the coeﬃcient of E0 [(vu − v0 )2 ], f3(32), is nonzero. The

numerical results reported below are those after this transformation. For ease of reference
to Table 1, the functional forms of fi and gij (i, j = 1, 2, 3) for the SV-CEV model aresummarized in Table 2.
It is important to note that there are many models proposed in the literature in which
the dynamics of the logarithmic volatility are directly specified: see, e.g., Andersen et al.
(2002), Chernov et al. (2003), and Scott (1987). An advantage of the logarithmic volatility
models is to ensure non-negativity of the volatility process without placing parameter
constraints. A diﬃculty, on the other hand, is to lose analytical tractability of option
prices. The proposed method may then have potential to overcome this diﬃculty.

4.2

Accuracy results for the SV-CEV model

As a base case, the parameter values for the SV-CEV model are set to(γ1 , κ2 , θ2 , σ2 , γ2 , ρ12 ) = (0.75, 1.0, 0.152 , 1.0, 1.0, −0.5) .
In particular, σ2 = γ2 = 1 is determined with reference also to HW (1987, Table II). Someof these values are then changed to further examine the cases in which the approximationis (in)accurate. It is noted that γ1 < 1 and −1 < ρ12 < 0 contribute to capturing theso-called leverage eﬀect: the volatility of the spot price tends to be high when the spot
price falls than rises. Various combinations of (γ1 , ρ12 ) are also tried, the results of which
are similar to those reported below.
Table 3 presents put option prices (multiplied by 100) and relative pricing errors(expressed in %) in various cases. Panel A of Table 3 presents the results of the basecase. First, the second-order approximation, AP2, overvalues the out-of-the-money (OTM,
K/S0 = 0.9) option by 8.4% and undervalues the at-the-money (ATM, K/S0 = 1.0) option by 2.0%. The approximation error for the OTM option appears large, however, theabsolute diﬀerence between the AP2 and MC prices is 0.086×10−2 , or equivalently 0.086%
of the spot price, S0 = 1. It is not surprising that the approximation error for the in-themoney (ITM, K/S0 = 1.1) option is very small, as the option value trivially approaches
e−rt K − S0 , irrespective of the computation methods involved. Second, the third-orderapproximation, AP3, achieves high accuracy. Specifically, the approximation error for the
OTM option is reduced to around 1% and that for the ATM option is negligibly small.

10

Panel B of Table 3 presents the results when the maturity length is extended to one
year, t = 1. In theory, the accuracy of the approximation decreases with increasing
maturity. In reality, however, the results are mixed between the ATM and OTM options.
While the former is not much aﬀected, the accuracy for the latter becomes worse thanthat in the base case. Still, AP3 seems eﬀective.
Next, the initial value of the volatility process is increased to V0 = 0.32 , assuming thecase where the current volatility level is unusually high relative to its long-term mean,
θ2 = 0.152 . Panel C of Table 3 shows that while AP2 becomes worse, especially for the
ATM and ITM options, AP3 continues to work well. The decrease in the accuracy due tothe larger V0 is mitigated when both V0 and θ2 are increased to 0.32 , as shown in Panel Dof Table 3. The results of Panels C and D indicate that the proposed method tends to be
more accurate when the current state values are around the long-term means than when
they are distant.
Increasing the volatility of Vt by doubling the value of σ2 results in the deteriorationof the accuracy, as shown in Panel E of Table 3. In fact, this is the most diﬃcult case forthe proposed method, where even AP3 undervalues the OTM option by nearly 30%. The
increase in σ2 makes the variation in Vt higher, which in turn makes the leptokurtosis of
the return distribution higher. In this case, the diﬃculty of the approximation arises. A
similar problem seems to appear in the HW approximation: HW (1987, p. 294). Note thathigher variation in Vt can also be generated by a larger V0 through the CEV specification.
In this case, however, the deterioration of the accuracy can be avoided by the third-orderapproximation, as shown in Panels C and D of Table 3.
On the other hand, increasing the speed of mean reversion of Vt by doubling thevalue of κ2 improves the accuracy, as shown in Panel F of Table 3. In particular, evenAP2 performs reasonably well. This improvement occurs because faster mean reversion
virtually reduces the variation in Vt and thus the leptokurtosis of the return distribution.
Next is the case where the increase in κ2 does not necessarily reduce the variation in
Vt . Given κ2 = 2, the value of σ2 is raised such that the unconditional variance of Vt ,(θ2 σ2 )2 /(2κ2 − σ22 ), remains the same as that in the base case. This is achieved by setting√
σ2 = 2 with θ2 unchanged. Panel G of Table 3 shows that due to the increase in σ2 ,the error magnitude generally increases from that in Panel F of Table 3. Compared tothe base case, the error magnitude for the OTM option is larger, indicating that the eﬀectof the large σ2 (higher variation) dominates that of the large κ2 (faster mean reversion).
Conversely, the error magnitude for the ATM option is smaller and that for the ITM
option is similar, indicating that the eﬀect of faster mean reversion prevails.
Although these cases are artificial and limited, they seem to reveal when the proposed
method is (in)accurate together with the magnitude of approximation errors. Overall, thegood performance of the approximation can be expected when the volatility process is not
too volatile. In reality, the volatility process does not seem too volatile when it is estimated11

from the time-series of the spot price. When estimated from the cross-section of option
prices, it is reported to be too volatile: see, e.g., Bakshi et al. (1997), and Bates (1996,
2000). This is because high variation in the volatility process is required for generatinghigh leptokurtosis of the return distribution, which in turn is required for explaining theimplied volatility surface, especially at short maturity. This too volatile behavior under
the risk-neutral measure may suggest that the diﬀusion component alone has diﬃculty ingenerating a suﬃcient degree of leptokurtosis. This diﬃculty is mitigated by introducing
jumps. Then, the proposed method may actually be more eﬀective when it is combinedwith the jump component. It is important to note that even after the inclusion of the jumpcomponent, the role of the diﬀusion component is reported to be decisive for capturing
not only time-varying second moment and negative skewness of the return distribution,
but also the term structure of implied volatilities: see Das and Sundaram (1999), Carrand Wu (2003b), and Huang and Wu (2004). Besides, diﬀusion models can be used fordescribing stochastic variations in the arrival rate of jumps: see Carr et al. (2003), and
Huang and Wu (2004). The proposed method may contribute to searching appropriate
models of the diﬀusion component, which is left for future study.

5

Identification of Eﬀective Higher-Order Moments

Based on the previous results that the third-order approximation improves the accuracy,
this section examines which of the conditional third moments added to Ψn=20,t have a greater
contribution to this improvement. It is useful to discriminate such eﬀective moments fromthose of little significance. The latter moments can be removed from Ψn=30,t to reducethe computational burden without much reducing the accuracy. The findings from thisanalysis are particularly useful for high-dimensional models, where the dimension of amoment vector increases rapidly with the number of factors: in general, the length ofΨs,t consisting of up to n-th conditional moments of a d-dimensional process is−1 = (n + d)!/(n!d!) − 1.

�n+d�
n

The following three moment vectors are considered:�
AP3–0** : (Ψn=20,tψ300,0,tψ210,0,tψ201,0,tψ120,0,tψ111,0,tψ102,0,t )� ,

�
AP3–*0* : (Ψn=20,tψ030,0,tψ210,0,tψ120,0,tψ111,0,t

ψ021,0,tψ012,0,t )� ,

�
AP3–**0 : (Ψn=20,tψ003,0,tψ201,0,tψ111,0,tψ102,0,t

ψ021,0,tψ012,0,t )� .

The moment vector labeled AP3–0** contains the conditional third moments except those
not related to st , ψ0∗∗ , leading to the abbreviation as –0**. The moment vectors labeled
AP3–*0* and AP3–**0 are constructed analogously. The dimension of these vectors is
fifteen.
Table 4 presents put option prices (multiplied by 100) and relative pricing errors (expressed in %) for the simplified third-order approximations in the base case: for ease of12

comparison, the previous results in Panel A of Table 3 are also presented. First, in termsof similarity to the performance of AP3, AP3–*0* ranks first. In fact, the error pattern
little changes. The result indicates that the conditional third moments related to vt have afundamental role for the improvement of the accuracy. This makes sense, as the accuracy
is measured in terms of option prices, which are highly sensitive to the volatility process.
Also, the error magnitude for AP3–**0 remains almost the same as that for AP3. On theother hand, while AP3–0** actually improves AP2, the error pattern is more similar to
AP2 than AP3. This result indicates that the conditional third moments related to st do
not contribute much to the improvement of the accuracy.
Further investigation reveals that among the conditional third moments related to vt(and Zt ), ψ021,0,t = E0 [(vt − v0 )2 (Zt − Z0 )] is the most influential. The following moment� ψ�vector of dimension ten, labeled AP2+021, is considered: (Ψn=2021,0,t ) . The relatives,tpricing errors for AP2+021 are also presented in Table 4, showing that the accuracy is
not much decreased from that of AP3 with the computational burden similar to that ofAP2. This result is robust to various cases considered in the previous section.

6

Concluding Remarks

This study proposed an approximation of European option prices that can be eﬃciently
computed under flexibly specified diﬀusion models of the spot price. An analytical approximation to the characteristic function of the log spot price is first obtained as the solution
to ordinary diﬀerential equations, and then Fourier inversion is applied to obtain option
prices. Using a stochastic volatility model in which both the spot price and volatility
processes have the constant elasticity of volatility specification, Monte Carlo simulations
revealed that the third-order approximation generally achieves high accuracy, except whenthe volatility process exhibits high variation, and that the second-order approximation iseﬀective in some cases, especially when mean reversion of the volatility process is relatively
fast.

13

References
Andersen, T. G., L. Benzoni, and J. Lund, 2002, “An Empirical Investigation of ContinuousTime Equity Return Models,h Journal of Finance, 57, 1239-1284.
Andersen, T. G., and J. Lund, 1997, “Estimating continuous-time stochastic volatility
models of the short-term interest rate,h Journal of Econometrics, 77, 343-377.
Bakshi, G., C. Cao, and Z. Chen, 1997, “Empirical performance of alternative option
pricing models,h Journal of Finance, 52, 2003-2049.
Barndorﬀ-Nielsen, O. E., 1998, “Processes of normal inverse Gaussian type,h Finance andStochastics, 2, 41-68.
Bates, D. S., 1996, “Jumps and stochastic volatility: exchange rate processes implicit in
Deutsche mark options,h Review of Financial Studies, 9, 69-107.
Bates, D. S., 2000, “Post-’87 crash fears in the S&P 500 futures option market,h Journalof Econometrics, 94, 181-238.
Brenner, R. J., R. H. Harjes, and K. F. Kroner, 1996, “Another look at models of theshort-term interest rate,h Journal of Financial and Quantitative Analysis, 31, 85-107.
Carr, P., H. Geman, D. B. Madan, and M. Yor, 2002, “The fine structure of asset returns:
an empirical investigation,h Journal of Business, 75, 305-332.
Carr, P., H. Geman, D. B. Madan, and M. Yor, 2003, “Stochastic volatility for Lévy
processes,h Mathematical Finance, 13, 345-382.
Carr, P., and L. Wu, 2003a, “The finite moment log stable process and option pricing,hJournal of Finance, 58, 753-777.
Carr, P., and L. Wu, 2003b, “What type of process underlies options? A simple robust
test,h Journal of Finance, 58, 2581-2610.
Chernov, M., A. R. Gallant, E. Ghysels, and G. Tauchen, 2003, “Alternative models forstock price dynamics,h Journal of Econometrics, 116, 225-257.

14

Das, S. R., and R. K. Sundaram, 1999, “Of smiles and smirks: a term structure perspective,h Journal of Financial and Quantitative Analysis, 34, 211-239.
Duﬃe, D., J. Pan, and K. Singleton, 2000, “Transform analysis and asset pricing for aﬃne
jump-diﬀusions,h Econometrica, 68, 1343-1376.
Gallant, A. R., and G. Tauchen, 1998, “Reprojecting partially observed systems withapplication to interest rate diﬀusions,h Journal of the American Statistical Association,
93, 10-24.
Heston, S. L., 1993, “A closed-form solution for options with stochastic volatility withapplications to bond and currency options,h Review of Financial Studies, 6, 327-343.
Heston, S. L., and S. Nandi, 2000, “A closed-form GARCH option valuation model,hReview of Financial Studies, 13, 585-625.
Huang, J., and L. Wu, 2004, “Specification analysis of option pricing models based ontime-changed Lévy processes,h Journal of Finance, 59, 1405-1440.
Hull, J., and A. White, 1987, “The pricing of options on assets with stochastic volatilities,hJournal of Finance, 42, 281-300.
Jones, C. S., 2003, “The dynamics of stochastic volatility: evidence from underlying andoptions markets,h Journal of Econometrics, 116, 181-224.
Madan, D. B., P. Carr, and E. C. Chang, 1998, “The variance gamma process and option
pricing,h European Finance Review, 2, 79-105.
Madan, D. B., and E. Seneta, 1990, “The variance gamma (V.G.) model for share market
returns,h Journal of Business, 63, 511-524.
Melino, A., and S. M. Turnbull, 1990, “Pricing foreign currency options with stochastic
volatility,h Journal of Econometrics, 45, 239-265.
Pan, J., 2002, “The jump-risk premia implicit in options: evidence from an integrated
time-series study,h Journal of Financial Economics, 63, 3-50.

15

Scott, L. O., 1987, “Option pricing when the variance changes randomly: theory, estimation, and an application,h Journal of Financial and Quantitative Analysis, 22, 419-438.
Scott, L. O., 1997, “Pricing stock options in a jump-diﬀusion model with stochastic volatility and interest rates: applications of Fourier inversion methods,h Mathematical Finance,
7, 413-426.
Shoji, I., 2002, “Approximation of conditional moments of diﬀusion processes,h International Journal of Computational and Numerical Analysis and Applications, 1, 163-190.
Takamizawa, H., and I. Shoji, 2009, “Modeling the Term Structure of Interest Rates with
General Diﬀusion Processes: A Moment Approximation Approach,h Journal of Economic
Dynamics and Control 33, 65-77.

16

[A]i,j

1

2

3

4

5

1

9

(1,0,0)f1(1,0,0)f2(1,0,0)f3(1,0,0)2f1 + g11
(1,0,0)g22(1,0,0)g33(1,0,0)f2 + g12(1,0,0)f3 + g13(1,0,0)g23(0,1,0)f1(0,1,0)f2(0,1,0)f3(0,1,0)g11(0,1,0)2f2 + g22
(0,1,0)g33(0,1,0)f1 + g12(0,1,0)g13(0,1,0)f3 + g23(0,0,1)
f1(0,0,1)
f2(0,0,1)
f3(0,0,1)g11(0,0,1)g22(0,0,1)
2f3 + g33
(0,0,1)g12(0,0,1)
f1 + g13(0,0,1)
f2 + g23(0,2,0)f1/2(0,2,0)f2/2(0,2,0)f3/2(0,2,0)g11/2(0,1,0)(0,2,0)
2f2
+ g22
/2
(0,2,0)g33/2(0,1,0)(0,2,0)f1+ g12/2(0,2,0)g13/2(0,1,0)(0,2,0)
f3
+ g23
/2

10

3g11

0

0

(2,0,0)f1/2(2,0,0)f2/2(2,0,0)f3/2(1,0,0)(2,0,0)
2f1
+ g11
/2
(2,0,0)g22/2(2,0,0)g33/2(1,0,0)(2,0,0)f2+ g12/2(1,0,0)(2,0,0)f3+ g13/2(2,0,0)g23/2(1,0,0)3f1 + 3g11

11

0

3g22

0

0

12

0

0

3g33

0

0
(0,1,0)g11

f1 + 2g12

2
3
4
5
6
7
8

17

Table 1:0(0,1,0)3f2 + 3g22

13

2g12

g11

0

14

2g13

0

g11

15

g22

2g12

0

16

g23

g13

g12

17

g33

0

2g13

(1,0,0)f2 + 2g12(1,0,0)f3 + 2g13(1,0,0)g22(1,0,0)g23(1,0,0)
g33

18

0

2g23

g22

0

f3 + 2g23

19

0

g33

2g23

0

g33

[b]j

f1

f2

f3

g11

0
(0,1,0)(0,1,0)g13

0(0,1,0)(0,1,0)

g22

6

7

(0,0,2)/2

f1(0,0,2)/2

f2(0,0,2)/2

f3(0,0,2)

/2

2f1

g22

(0,0,2)/2

2f2(0,0,1)

+ g33

f1
f2
f3

g11

(0,0,2)2f3(0,0,2)g12(0,0,1)+ g13(0,0,1)

+ g23

f1
f2

(1,1,0)f1(1,1,0)f2(1,1,0)f3(0,1,0)+ g11(1,0,0)+ g22(1,1,0)(1,0,0)(0,0,2)/2

f3(0,0,2)/2

f3(0,1,0)+ f2(1,0,0)+ g23(1,1,0)(0,1,0)(0,0,1)+ g12(1,0,0)+ g33(0,0,1)+ g12(0,1,0)

+ g13

2f3

(1,0,1)

f1(0,0,1)+ f3(1,0,0)f2(0,1,0)(1,0,1)(1,0,1)+ g13(1,0,1)+ g23

f1(0,1,0)f2(0,1,1)(0,1,1)(0,1,1)(0,1,1)(0,0,1)+ f3(0,1,1)

+ g23

18

g13

g23

3g11

(0,1,0)+ 2g12(0,1,0)2g13(0,1,0)+ 2g12(1,0,0)+ g23

f3 + g13(1,0,0)(0,1,0)f1 + 2g13

g33(0,0,1)

2g23

g33

+ g33

f1

+ g22

0

(0,0,1)(0,0,1)(1,0,0)f2(0,0,1)2f2

0

2f2 + g22

f2 + 2g23(0,1,1)g11(1,0,0)3g33(0,0,1)2g12(1,0,0)(0,0,1)
2f1 + g11
+ 2g13(0,0,1)g22(1,0,0)(0,0,1)
f2 + g12+ g23(0,0,1)(1,0,0)2f3 + g33
+ 2g13(1,0,0)g22(1,0,0)
2g23

0

g22

(1,0,1)+ g11

0(1,0,0)(0,0,1)(0,1,1)(0,0,1)3g22(0,1,0)3g33(0,0,1)g11(0,1,0)g11(0,0,1)2g12(0,1,0)(0,0,1)
f1 + g12+ g13(0,1,0)2g13(0,1,0)(0,0,1)
2f2 + g22
+ 2g23(0,0,1)(0,1,0)2f3 + g33
+ 2g23

2f1 + g11

g12

f3

(1,0,0)3g22

0(0,0,1)(1,0,1)(0,1,1)(0,0,1)(0,0,1)
3f3 + 3g33

Table 1 (continued):(1,1,0)(1,1,0)

3g11

g11

f2

(1,0,1)

+ g12

+ g13

0

0

(1,0,1)

2f3(0,1,0)(0,1,1)

g22

g33
f1

f1

(0,0,1)(1,1,0)9(1,0,1)

2f1(1,1,0)/2/2

8(1,0,0)(1,0,0)g33

g12(0,1,0)10

11(3,0,0)/6

f1(3,0,0)

/6

f2

f1
f2

12

(0,3,0)/6

f1(0,3,0)

/6

f2

13

(0,0,3)/6

f1(2,1,0)/2(0,0,3)/6

f2(2,1,0)/2(0,3,0)f3/6(0,3,0)g11/6(0,2,0)(0,3,0)f2+ g22/6(0,3,0)g33/6(0,2,0)(0,3,0)f1/2 + g12/6(0,3,0)g13/6(0,2,0)(0,3,0)f3/2 + g23/6(0,0,3)f3/6(0,0,3)g11/6(0,0,3)g22/6(0,0,2)(0,0,3)f3+ g33/6(0,0,3)g12/6(0,0,2)(0,0,3)f1/2 + g13/6(0,0,2)(0,0,3)f2/2 + g23
/6

0

0

0

(0,1,0)(0,2,0)3f2+ 1.5g22

0(2,1,0)f3/2(1,1,0)(2,1,0)
2f1
+ g11
/2
(2,0,0)(2,1,0)f2+ g22/2(2,1,0)g33/2(2,0,0)(1,1,0)(2,1,0)f1/2 + f2+ g12/2(1,1,0)(2,1,0)f3+ g13/2(2,0,0)(2,1,0)f3/2 + g23/2(0,1,0)(1,1,0)
3f1+ 3g11(2,0,0)1.5g220(0,0,1)(0,0,2)
3f3+ 1.5g33

019(3,0,0)f3/6(2,0,0)(3,0,0)f1+ g11/6(3,0,0)g22/6(3,0,0)g33/6(2,0,0)(3,0,0)f2/2 + g12/6(2,0,0)(3,0,0)f3/2 + g13/6(3,0,0)g23/6(1,0,0)(2,0,0)
3f1+ 1.5g11

0(1,0,0)+ g12(1,0,0)

+ g13

f2
f3

(2,0,0)(2,0,0)(2,0,0)/2(2,0,0)/2(2,0,0)

/2

g22
g23
g33

0
0
0

Table 1 (continued):(0,2,0)g11/2(0,0,2)0(0,1,0)f1

g11(0,2,0)+ g12(0,2,0)g13(0,1,0)f3(0,0,2)+ g13

0(2,0,0)(0,1,0)+ 2g13f3(1,0,0)(0,0,1)(1,0,0)f3(1,1,0)+ g22(2,0,0)

+ g13

g33

g22

/2

g23

(0,0,1)+ g23(0,0,2)g33(0,0,2)f2

0(1,1,0)/2 + 2g12(1,1,0)(0,0,2)

f1
/2

+ g11

2f2
/2

(0,2,0)(0,2,0)(0,1,0)

+ f2

/2

g12

+ g23

g33

2f1

0

/2

0

(1,0,0)0(1,1,0)/2 + g23(1,1,0)(2,0,0)(2,0,0)0(2,0,0)

+ g12

/2

14

15

(2,0,1)/2

f1(2,0,1)

/2

f2

f3

(2,0,1)/2

f3(1,0,1)

+ g11

f1
f2

(2,0,1)

2f1(2,0,1)

g22

+ g33

(1,0,1)f2(2,0,0)f1(2,0,0)

f2

/2

+ g12

(2,0,1)/2(1,0,1)+ g13(2,0,1)(2,0,1)/2 + g23(0,0,1)f1(1,2,0)/2

f2(1,2,0)/2

f3(0,2,0)+ g11(1,1,0)+ g22(1,2,0)/2 + g12(0,2,0)(1,2,0)(1,1,0)(1,2,0)20(1,0,0)(0,0,1)

+ g11

+ f3

(1,0,1)(2,0,0)(0,1,0)2f1(1,0,1)(1,0,1)(1,0,0)f2(1,0,0)2f3(2,0,0)+ g12(0,1,0)(1,0,1)(1,0,1)+ g33(2,0,0)g22/2(2,0,0)(1,0,1)+ f2(1,1,0)+ f3(1,1,0)

+ f3

/2

f1

/2

f2

(0,1,0)f3(0,1,1)+ g13(1,0,1)

+ g23

3g33
(0,2,0)+ g12(0,0,1)+ g11(0,1,0)+ g11(1,1,0)/2 + 2g12+ 2g12(1,1,0)+ 2g13(0,0,1)+ g22(0,1,1)+ 2g12(0,2,0)(0,0,1)(1,1,0)2f2

g33(0,2,0)/2(1,0,0)+ 2g232f2(1,1,0)2f3(1,1,0)0(1,1,1)(1,0,1)

2f1
2f1(0,2,0)0(1,1,1)(1,1,0)

+ g22

g33

(1,1,1)(1,0,1)

+ g23

f3

+ g12

(0,1,1)(1,1,0)(2,0,0)(0,1,1)

0

+ g13

+ g13

(1,1,1)3g22(1,1,0)(1,0,0)(1,1,1)3g11(1,1,0)

+ g11

+ f1

(1,1,1)+ 3g22

g23

Table 1 (continued):

+ g33

f1

g13
2f2

/2 + g23/2(0,2,0)/2 + 2g13

g22(1,1,0)(0,2,0)(2,0,0)(1,0,0)+ g22

1.5g11

1.5g332f1(1,0,1)2f2(1,2,0)

+ g23

3f2

+ 2g12+ g11

2f3/2 + g13

f3(0,1,1)2f1/2(0,2,0)f3

0(0,0,1)/2+ f2(1,1,0)+ 3g11f2(1,1,1)

/2

f1

/2

(1,1,1)(1,2,0)(1,2,0)/2(1,1,1)g33(1,0,1)

3f1

/2

2f2

(2,0,1)/2 + f3(1,2,0)f1/2(2,0,0)

f3

/2

16

/2

(1,0,0)f1(0,1,0)

+ f2

+ f3

(0,1,1)(1,0,1)+ g12(0,1,0)+ g33(1,0,0)+ g22(1,0,0)

+ g33

2f3

(0,1,1)(1,0,1)+ g13(0,1,1)+ 2g13(1,1,0)+ 2g23(1,0,1)+ 2g230(1,1,0)(1,0,1)(1,1,0)(0,1,1)

+ g23

17

18

(1,0,2)/2

f1(1,0,2)

/2

f2

f1
f2

19

(0,2,1)/2

f1(0,1,2)/2(0,2,1)/2

f2(0,1,2)/2(0,2,1)f3/2(0,2,1)g11/2(0,1,1)(0,2,1)
2f2
+ g22
/2
(0,2,0)(0,2,1)f3+ g33/2(0,1,1)(0,2,1)f1+ g12/2(0,2,0)(0,2,1)
f1/2 + g13/2(0,2,0)(0,1,1)(0,2,1)f2/2 + f3+ g23/2(0,1,2)f3/2(0,1,2)g11/2(0,0,2)(0,1,2)f2+ g22/2(0,1,1)(0,1,2)
2f3
+ g33
/2
(0,0,2)(0,1,2)f1/2 + g12/2(0,1,1)(0,1,2)f1+ g13/2(0,1,1)(0,0,2)(0,1,2)f2+ f3/2 + g23
/2

0

0

0

(0,0,1)(1,0,1)(0,0,2)
2f1
+ g11
+ g13
(0,0,2)g22/2(0,0,1)(1,0,1)(0,0,2)
f2
+ g12
+ g23
/2
(0,0,1)(1,0,0)(0,0,2)(1,0,1)
2f3
+ f1
+ g33
/2 + 2g13(1,0,1)g22(1,0,0)(1,0,1)
f2+ 2g23(0,0,1)(0,1,1)
3f2+ 3g22(0,2,0)1.5g33(0,1,1)g11(0,2,0)g11/2(0,0,1)(0,1,1)
f1+ 2g12(0,1,0)(0,2,0)(0,1,1)f1+ g12/2 + g13(0,2,0)g13(0,1,0)(0,0,1)(0,2,0)(0,1,1)
2f2
+ f3
+ g22
/2 + 2g23(0,1,0)(0,1,1)(0,2,0)
2f3
+ g33
+ g23

(0,1,1)g11(0,0,2)g12(0,0,1)(0,1,1)(0,0,2)
f1
+ g12
+ g13
/2
(0,1,0)(0,1,1)
f1+ 2g13(0,0,1)(0,1,1)(0,0,2)
2f2
+ g22
+ g23
(0,0,1)(0,1,0)(0,0,2)(0,1,1)
2f3
+ f2
+ g33
/2 + 2g23

0

0

0

21

(1,0,2)f3/2(0,0,2)(1,0,2)f1+ g11/2(1,0,2)g22/2(1,0,1)(1,0,2)
2f3
+ g33
/2
(0,0,2)(1,0,2)f2/2 + g12/2(1,0,1)(0,0,2)(1,0,2)f1+ f3/2 + g13/2(1,0,1)(1,0,2)f2+ g23/2(0,0,2)
1.5g11(1,0,0)3f3(1,0,1)+ 3g33(0,0,2)

g12

Table 1 (continued):(0,0,2)

1.5g22(0,1,0)3f3(0,1,1)+ 3g33(0,0,2)g11/2

Table 1 (continued): The elements of A(Xs ) and b(Xs ) for the third-order approximation!
The moment vector Ψn=3s,t consists of up to the conditional third moments of a three-dimensional process, Xt = (x1,t x2,t x3,t ) , in the following
order: define ψpqr,s,t = Es [(x1,t − x1,s )p (x2,t − x2,s )q (x3,t − x3,s )r ], andΨn=2= (ψ100,s,ts,tψ010,s,t

ψ001,s,t

ψ200,s,t

ψ020,s,t

ψ002,s,t

ψ110,s,tψ101,s,t

ψ011,s,t )! ,!
Ψn=3
= (Ψn=2s,ts,tψ300,s,t

ψ030,s,t

ψ003,s,t

ψ210,s,tψ201,s,t

ψ120,s,t

ψ111,s,tψ102,s,t

ψ021,s,t

ψ012,s,t )! .

n=3
A(Xs ) and b(Xs ) for Ψn=2
s,t are the 9 × 9 upper left submatrix of A(Xs ) and the 9 × 1 subvector of b(Xs ) for Ψs,t , respectively. The SDEs
for xi,t (i = 1, 2, 3) are dxi,t = fi (Xt )dt + ξi (Xt )! dWt , and the notations in the table are as follows: gij = ξi! ξj ,22(k,l,m)fi.

=∂fik+l+m,∂xk1 ∂xl2 ∂xm
3and(k,l,m)

gij

=

k+l+m
∂gij

∂xk1 ∂xl2 ∂xm
3f1 (st , vt , Zt )=Et [dst ]/dt=

rt − g11 /2

f2 (st , vt , Zt )=Et [dvt ]/dt=

κ2 (θ2 exp (−vt ) − 1) − g22 /2

f3 (st , vt , Zt )=Et [dZt ]/dt=(−rt + αf1 + α2 g11 /2)Zt

g11 (st , vt , Zt )=Et [(dst )2 ]/dt=

g22 (st , vt , Zt )=Et [(dvt )2 ]/dt

exp (2(γ1 − 1)(st − s0 ) + vt )=g33 (st , vt , Zt )=Et [(dZt )2 ]/dt=

(αZt )2 g11

g12 (st , vt , Zt )=Et [dst dvt ]/dt=

σ2 ρ12 exp ((γ1 − 1)(st − s0 ) + (γ2 − 0.5)vt )g13 (st , vt , Zt )=Et [dst dZt ]/dt=

αZt g11g23 (st , vt , Zt )=Et [dvt dZt ]/dt=

αZt g12σ22 exp (2(γ2 − 1)vt )Table 2: The functions based on which the elements of A and b are computed
The SV-CEV model has the following SDEs for the log spot price process, st , and the log
volatility process, vt :

dvt

�

�

1
rt − e2(γ1 −1)(st −s0 )+vt dt + e(γ1 −1)(st −s0 )+vt /2 dW1,t ,
2
�
�
1
=
κ2 (θ2 e−vt − 1) − σ22 e2(γ2 −1)vt dt + σ2 e(γ2 −1)vt dW2,t .
2

dst =

The SDE for Zt is given by Eqs.(25)–(27). The order of the listed functions corresponds
to the moment vector, Ψn=2
s,t , given by Eq.(15).

23

K/S0

Put Option Prices (×100)
0.9
1.0
1.1

Relative Pricing Errors (%)
0.9
1.0
1.1

Panel A: Base case
AP2
1.115
AP3
1.017
MC
1.029

4.057
4.136
4.141

10.89810.84310.865

8.43
−1.17

−2.02
−0.11

0.31
−0.20

Panel B: t = 1
AP2
2.524
AP3
2.172
MC
2.229

5.703
5.785
5.782

12.26412.00012.039

13.28
−2.54

−1.36
0.06

1.87
−0.32

Panel C: V0 = 0.32
AP2
3.998
AP3
3.558
MC
3.605

7.118
7.468
7.533

13.83913.511
13.556

10.91
−1.29

−5.51
−0.86

2.09
−0.33

Panel D: V0 = 0.32 and θ2 = 0.32
AP2
4.394
8.060
AP3
4.173
8.244
MC
4.175
8.238

14.39814.21514.218

5.23
−0.05

−2.16
0.07

1.26
−0.02

Penel E: σ2 = 2
AP2
1.535
AP3
0.769
MC
1.092

3.621
3.813
3.920

10.84010.237
10.66540.65−29.55

−7.63
−2.71

1.65
−4.01

Panel F: κ2 = 2
AP2
1.071
AP3
1.001
MC
1.011

4.128
4.161
4.157

10.916
10.88510.895

5.90
−1.02

−0.70
0.10

0.20
−0.09

√
Panel G: κ2 = 2 and σ2 = 2
AP2
1.177
4.046
AP3
1.000
4.098
MC
1.049
4.099

10.86210.767
10.820

12.26
−4.63

−1.30
−0.02

0.39
−0.49

Table 3: Accuracy results in various cases
The relative pricing error is defined as pap /pmc − 1, where pap and pmc denote the option
prices computed by the approximation and MC methods, respectively. The SV-CEV
model is originally specified as
dSt /St = rt dt +�

Vt (St /S0 )γ1 −1 dW1,t ,dVt = κ2 (θ2 − Vt )dt + σ2 Vtγ2 dW2,t ,with Et [dW1,t dW2,t ] = ρ12 dt. The parameter values are set to (γ1 , κ2 , θ2 , σ2 , γ2 , ρ12 ) =(0.75, 1.0, 0.152 , 1.0, 1.0, −0.5). The input values are as follows: S0 = 1; K = 0.9, 1.0, 1.1;
t = 0.5 (half a year); V0 = 0.152 ; and rt = 0 (constant). These parameter and input values
form a base case. Some of them are then changed as indicated in each panel.

24

K/S0

Put Option Prices (×100)
0.9
1.0
1.1

Relative Pricing Errors (%)
0.9
1.0
1.1

AP2
AP3

1.115
1.017

4.057
4.136

10.89810.843

8.43
−1.17

−2.02
−0.11

0.31
−0.20

AP3–0**
AP3–*0*AP3–**0
AP2+021

1.091
1.014
1.041
1.047

4.066
4.137
4.156
4.147

10.87810.84510.882
10.886

6.10
−1.38
1.20
1.83

−1.82
−0.11
0.35
0.14

0.12
−0.19
0.16
0.20

Table 4: Accuracy results for the simplified third-order approximations in thebase case
The simplified third-order approximations correspond to the following moment vectors:
�
AP3–0** : (Ψn=20,tψ300,0,tψ210,0,tψ201,0,tψ120,0,tψ111,0,tψ102,0,t )� ,

�
AP3–*0* : (Ψn=20,tψ030,0,tψ210,0,tψ120,0,tψ111,0,t

ψ021,0,tψ012,0,t )� ,

�
AP3–**0 : (Ψn=20,tψ003,0,tψ201,0,tψ111,0,tψ102,0,t

ψ021,0,tψ012,0,t )� ,where Ψn=20,t is a nine dimensional vector consisting of the conditional second moments
of increments of (st , vt , Zt ) and ψpqr,0,t = E0 [(st − s0 )p (vt − v0 )q (Zt − Z0 )r ]. AP2+021
� ψ�
corresponds to the moment vector as (Ψn=2021,0,t ) .
0,t25