Latext Writer

University of Wollongong Research Online

Faculty of Engineering and Information Sciences - Papers: Part A

Faculty of Engineering and Information Sciences

2014

Exact and approximate solutions for options with time-dependent stochastic volatility Joanna Goard University of Wollongong, [email protected]

Research Online is the open access institutional repository for the University of Wollongong. For further information contact the UOW Library:

[email protected]

Publication Details Goard, J. (2014). Exact and approximate solutions for options with time-dependent stochastic volatility. Applied Mathematical Modelling, 38 (11-12), 2771-2780.

Exact and approximate solutions for options with time-dependent stochastic volatility

Abstract In this paper it is shown how symmetry methods can be used to >nd exact solutions for European option pricing under a time-dependent 3/2-stochastic volatility model View the MathML source. =is model with A(t) constant has been proven by many authors to outperform the Heston model in its ability to capture the behaviour of volatility and >t option prices. Further, singular perturbation techniques are used to derive a simple analytic approximation suitable for pricing options with short tenor, a common feature of most options traded in the market.

Keywords solutions, approximate, options, time, volatility, exact, dependent, stochastic

Disciplines Engineering | Science and Technology Studies

Publication Details Goard, J. (2014). Exact and approximate solutions for options with time-dependent stochastic volatility. Applied Mathematical Modelling, 38 (11-12), 2771-2780.

=is journal article is available at Research Online: h?p://ro.uow.edu.au/eispapers/2283

Exact and Approximate Solutions for Options with

Time-dependent Stochastic Volatility

Joanna Goard

University of Wollongong,

Northfields Ave, Wollongong, NSW, 2522, Australia

email: [email protected]

Ph: +61 2 42214188, Fax: +61 2 42214845

Abstract

In this paper it is shown how symmetry methods can be used to find exact solutions

for European option pricing under a time-dependent 3/2-stochastic volatility model dv =

kv(A(t)°v)dt+bv 3 2 dZ. This model with A(t) constant has been proven by many authors

to outperform the Heston model in its ability to capture the behaviour of volatility and

fit option prices. Further, singular perturbation techniques are used to derive a simple

analytic approximation suitable for pricing options with short tenor, a common feature

of most options traded in the market.

Keywords and phrases: stochastic volatility, volatility model, option pricing

Mathematics subject classification: 91G20, 35A22, 22E05

1 Introduction

Stochastic volatility models have become particularly popular for derivative pricing and hedg-

ing over the last 20 years. Empirical evidence on underlying asset prices and on their deriva-

tives strongly suggest that asset volatility is stochastic and that the Black-Scholes (BS) model

[1] is inadequate to describe features of asset returns such as skewness, leptokurtosis and pro-

nounced conditional heteroskedasticity. The literature supports the use of stochastic volatility

in an eÆort to reproduce the implied volatility smile observed in markets and to avoid many

of the shortcomings of the constant variance diÆusion assumed by BS.

Work on stochastic volatility models include that of Hull and White [2], Scott [3] and Wig-

gins [4], whose pricing solutions depended on extensive use of numerical techniques to solve

two-dimensional partial diÆerential equations (PDEs). Assuming that volatility is uncorre-

lated with the spot asset, Stein and Stein [5] manage to analytically value options. However,

volatility shocks are known to be negatively correlated with asset price shocks so that when

volatility increases, stock prices decrease and vice-versa. This is commonly referred to as

leverage and at least partially accounts for a skewed distribution for the asset price. Thus

without correlation, the prices cannot hope to capture such important skewness eÆects arising

from such correlation.

One of the most popular stochastic volatility models to date is the mean-reverting model

of Heston [6]. Using Fourier inversion methods, Heston provides a closed form solution for

the price of European options when the spot asset is correlated with volatility, and variance

v, follows the ‘square-root’ process dv = k(µ ° v)dt + æ p

vdZ, which is commonly referred

to as the ‘Heston’ model. [Note that here and elsewhere in the paper, Z and also Z1 and

Z2 will refer to Wiener processes under a real measure P while Z̃, Z̃1 and Z̃2 will refer

to Wiener processes under a corresponding equivalent risk-neutral measure Q.] However,

numerous empirical studies show that the Heston model is misspecified. Chacko and Viceira

1

[7] performed a comprehensive empirical analysis on variance models of the form

d∫ = (a + b∫)dt + c∫∞dZ (1.1)

where the instantaneous standard deviation of variance was allowed to be proportional to

any power of variance. Using the estimation technique of spectral GMM (generalised method

of moments) they found that the best value of ∞ was between 1 and 2, with the standard

errors indicating that the diÆerences between the values found for ∞ and one half (as in the

Heston model) were statistically significant. This means that compared to the Heston model,

the volatility of the variance process is more highly sensitive to the level of variance. Other

empirical studies have yielded similar results. For example, Jones [8] analyses the more general

CEV (constant elasticity of variance) model using a bivariate series of returns and an at-the-

money short maturity option. He finds that periods of high volatility coincide with periods

of volatile volatility which is in disagreement with the Heston model, in which the volatility

of instantaneous volatility is not level-dependent. In particular a number of his specification

tests favour the nona±ne CEV model over the Heston model with best estimates of gamma

as in (1.1), that are greater than one.

More recently, using various samples of S&P500 index return data, ChrisoÆersen et al [9]

compare the performance of the Heston model with 5 simple alternatives, all of which can be

described by:

dv = kva(µ ° v)dt + ævbdZ a = {0, 1}, b = {1/2, 1, 3/2}. (1.2)

They found that ‘the 3/2N’ model (a = 1, b = 3/2) and ‘the ONE’ model (a = 0, b = 1)

performed best as far as maximising model fit for the samples. The Heston model ranked

4th or 5th. Using out-of-the-money S&P500 index option data for 1996-2004, in minimising

option implied volatility mean square error, the ONE and 3/2N model outperformed the

2

Heston model for both in -sample and out-of-sample experiments. Thus the ‘3/2N’ model

namely

dS = µSdt + p

vSdZ1 (1.3a)

dv = (!v ° µv2)dt + ªv 3 2 dZ2, (1.3b)

has proven empirically to outperform the Heston model. Similarly, Goard and Mazur [10]

show that the 3/2-model outperforms the Heston (and all other models considered in their

test) in its ability to fit the VIX time series. The novel features of (1.3b) are 1) the specification

for the diÆusion having a high power law of 1.5 which can reduce the heteroskedasticity of

volatility1 and 2) a nonlinear drift so that it exhibits substantial nonlinear mean-reverting

behaviour when the volatility is above its long-run mean. Hence after a large volatility spike,

the volatility can potentially quickly decrease while after a low volatility period it can be slow

to increase. It has also been shown that with (1.3b), with µ > 0, v will always remain positive.

Further aspects of the model have been studied by Lewis [12], who using a risk-adjustment via

utility theory provides an analytic solution for option prices under this model, namely for a call

option with strike price K, time to maturity ø and with ∞ ∑ 1 and ∞(1 ° ∞)ª2 ∑ (µ + ª2/2)2

C(S, v, ø ) = S ° Ke°rø

2º

Z

ik

i

°1

ik

i

°1 e°iuy

H(u, v, ø )

u2 ° iu du (1.4)

where 0 < k i

< 1, y = log(S/K) + rø,

H(u, v, ø ) = °(Ø ° Æ)

°(Ø)

h

X( !̄

v , !ø )

i

Æ

M h

Æ, Ø,°X( !̄

v , !ø )

i

,

X(x, t) = x

et ° 1 , µ =

1

2 (1 + µ̂), ± = [µ2 + c̃]

1 2 ,

Æ = °µ + ±, Ø = 1 + 2±, !̃ = 2!

ª2 , c̄ =

(u2 ° iu) ª2

,

µ̂(u) = °1 + 2

ª2 [ p

(µ + ª2/2)2 ° ∞(1 ° ∞)ª2 + (1 ° ∞ + iu)Ωª].

where M (a, b, x) is the Kummer-M confluent hypergeometric function (see e.g. Abramowitz

1See Campbell et al [11] for a discussion on this for the short interest rate model of the same form.

3

and Stegun [13]) and r is the risk-free interest rate. Heston [14] also independently developed

the solution but with 2 integral functions corresponding to 2 cumulative distribution functions.

However, as Jones [8] and many other authors have found from empirical studies, although

the CEV processes greatly improve over the square-root specification, they still fail to match

observed prices of short-dated options and there still remain questions regarding long-term

memory in volatility.

In Section 2 of this paper we outline the Lie symmetries method to solve PDEs and in

Section 3 of this paper show how symmetries can be used to find solutions to European options

on an underlying with price S, under the time-dependent risk-neutral 3/2-model

dS = rSdt + p

vSdZ̃1 (1.5a)

dv = kv(A(t) ° v)dt + bv 3 2 dZ̃2, (1.5b)

with non-zero correlation, Ω, between price and volatility shocks, with |Ω| < 1. This is

a promising candidate for a realistic model for volatility from a financial viewpoint. It not

only includes a diÆusion term in agreement with empirical findings but also includes a realistic

time-dependent drift component in which we are free to choose the moving target. Thus model

(1.5a,b) with its infinitely more degrees of freedom would be able to describe a wide range of

volatility data. For example, using a truncated Fourier series for A(t) would be able to pick up

the main trends of cyclic behaviour from historical records. Further the free function of time

lends itself to calibration so that theoretical and current market prices can be matched for all

maturities (see e.g Wilmott [15] for discussions on calibration). [We note that with F t

denoting

the futures price with maturity T at time t, X t

= ln( Ft F0

), < X t

>= R

t

0 v

s

ds, Carr and Sun [16]

provide the joint conditional Fourier Laplace transform of X T

and < X T

> ° < X t

>. This

relies in part in ‘guessing’ certain forms for transformations.] The form of the exact solution

however still involves the valuation of an integral with complex arguments. This means that

care needs to be taken in numerical calculations and calibration of market and theoretical

4

prices may not be easy. In Section 4 of this paper we obtain a simple, analytic approximation

for prices of put options with short tenor. This approximation provides quick and accurate

values for options with short expiries, such as one or two months. In fact it is these types

of options that dominate the options markets. Then in Section 5, we compare the exact and

approximate solutions using diÆerent maturities and parameters. In Section 6 we present a

short conclusion.

2 Method

A symmetry of a diÆerential equation is a transformation mapping an arbitrary solution to

another solution of the diÆerential equation. The classical Lie groups of point invariance

transformations depend on continuous parameters and act on the system’s graph space that

is co-ordinatised by the independent and dependent variables. As these symmetries can be

determined by an explicit computational algorithm (known as Lie’s algorithm or Lie’s clas-

sical method), many automated computer algebra packages (see e.g Sherring [17]) have been

developed to find them. Thus they are the most extensively used of all symmetries.

If a PDE is an invariant under a point symmetry, one can often find similarity solutions

or invariant solutions which are invariant under some subgroup of the full group admitted

by the PDE. These solutions result from solving a reduced equation in fewer variables.

In essence, the classical method for finding symmetry reductions of a second-order PDE

in one dependent variable P and 3 independent variables (S, v, ø )

¢(S, v, ø, P S

, P v

, P ø

, P SS

, P Sv

, P Sø

, P vv

, P vø

, P ø ø

) = 0, (2.1)

5

is to find a one-parameter Lie group of transformations in infinitesimal form

S§ = S + "£(S, v, ø, P ) + O("2) (2.2a)

v§ = v + "V (S, v, ø, P ) + O("2) (2.2b)

ø§ = ø + "T̄ (S, v, ø, P ) + O("2) (2.2c)

P § = P + "¥(S, v, ø, P ) + O("2) (2.2d)

which leaves (2.1) invariant. The coe±cients £, V, T̄ and ¥ of the infinitesimal symmetry are

often referred to as the ‘infinitesimals’. This invariance requirement is determined by

G(2)¢|¢=0 = 0, (2.3)

where

G = £(S, v, ø, P ) @

@S + V (S, v, ø, P )

@

@v + T̄ (S, v, ø, P )

@

@ø + ¥(S, v, ø, P )

@

@P (2.4)

are vector fields which span the associated Lie algebra, and are called the infinitesimal gen-

erators of the transformation (2.2a-d) and G(2) is the second extension (or second prolonga-

tion) of G, extended to the second jet space, co-ordinatised by S, v, ø, P, P S

, P v

, P ø

, P SS

· · · (see

Chapter 2 in the book of Bluman and Kumei [18]).

Equation (2.3) is a polynomial equation in a set of independent functions of the derivatives

of P . As the equation must be true for arbitrary values of these independent functions, their

coe±cients must vanish, leading to an over-determined linear system of equations, known

as the determining equations for the coe±cients £(S, v, ø, P ), V (S, v, ø, P ), T̄ (S, v, ø, P ) and

¥(S, v, ø, P ). Then for known functions £, V, T̄ , ¥, invariant solutions P corresponding to

6

(2.2a-d) satisfy the invariant surface condition (ISC)

≠ = £(S, v, ø, P ) @P

@S + V (S, v, ø, P )

@P

@v + T̄ (S, v, ø, P )

@P

@ø ° ¥(S, v, ø, P ) = 0, (2.5)

which when solved as a first-order PDE by the method of characteristics, yields the functional

form of the similarity solution in terms of an arbitrary function, i.e

P = q(S, v, ø, ¡(z1, z2)),

where

z1 = z1(S, v, ø ), z2 = z2(S, v, ø ),

and where ¡ is an arbitrary function of invariants z1, z2 for the symmetry. Substituting this

functional form into (2.1) produces a quotient PDE in two independent variables which one

solves for the function ¡(z1, z2).

Further, for an initial-value problem with the initial condition P (S, v, 0) = j(S, v), then

we need a linear combination of generators such that condition (2.5) is satisfied at ø = 0,

P = j(S, v) i.e

£(S, v, 0, j(S, v))P S

(S, v, 0) + T̄ (S, v, 0, j(S, v))P ø

(S, v, 0)

+V (S, v, 0, j(S, v))P v

(S, v, 0) = ¥(S, v, 0, j(S, v)). (2.6)

P S

(S, v, 0) and P v

(S, v, 0) can be found from the final condition. As well, for evolution equa-

tions P ø

(S, v, 0) can be found from the governing PDE (see Goard [19] for details).

7

3 Exact Solution for European Put Option Price under

the Time-Dependent 3/2-Model

We assume that the risk-neutral process for the asset price is given by (1.5a,b) where

corr (dZ̃1, dZ̃2) = Ωdt, r is the constant risk-free rate and v the variance. In (1.5b), A(t) is the

long-run target of the variance process, k is the reversion rate to the long-run target and b is

the volatility of volatility.

Theorem 3.1: The value of a European put option with strike K and expiry T , when the

underlying asset follows the risk-neutral process (1.5a,b), is given by

P (S, v, ø ) = e°rø [K + √(Serø , vT̄ (ø ))] (3.1)

where

T̄ (ø ) =

8

>

>

>

>

<

>

>

>

>

:

1 r

h

ek R

ø

0 a(u)du ° a0

a(ø ) ° a0ek

R

ø

0 a(u)du

R

ø

0

n

a

0(u) a(u)2

e°k R

u

0 a(u1)du1

o

du i

, if a = a(ø )

1 r

(eakø ° 1), if a is constant,

a(ø ) = A(T ° ø ), a0 = a(0),

√(x, y) = °K 2º

Z

iw

i

+1

iw

i

°1

∑

x°iwy°m Kiw

w2 ° iw °(m + 1 ° c1

b

)

°(2m + 1 ° c1 b

)

µ

2ka0 b2r

∂

m

£M µ

m, 2m + 1 ° c1 b

, °2ka0 b2ry

∂∏

dw

(see Note 1 below),

M is the confluent hypergeometric function, (see e.g Abramowitz and Stegun [13]),

and ø = T ° t, 0 < w i

< 1, c1 = °b ° 2Ωiw ° 2k

b , i2 = °1, m =

c1 + p

c21 ° 4wi + 4w2 2b

.

8

Note 1: In practice, it can be convenient to perform the integration in the real plane and

so we compute √(x, y) with w ! [i/2 + w r

] and use

√(x, y) = °K º

Z

w̄

0

R Ω

x°i[i/2+wr]y°m Ki[i/2+wr]

[i/2 + w r

]2 ° i[i/2 + w r

]

°(m + 1 ° c1 b

)

°(2m + 1 ° c1 b

)

µ

2ka0 b2r

∂

m

M

µ

m, 2m + 1 ° c1 b

, °2ka0 b2ry

∂æ

dw r

where w̄ is su±ciently large.

Proof: Given that the stock price follows (1.5a,b), the price of a European put option

P (S, v, ø ) with expiry T and strike K satisfies the following PDE

vS2

2

@2P

@S2 + Ωbv2S

@2P

@S@v +

b2v3

2

@2P

@v2 + rS

@P

@S + kv(a(ø ) ° v)

@P

@v ° rP °

@P

@ø = 0 (3.2a)

where a(ø ) = A(T ° ø ), and subject to

P (S, v, 0) = max(K ° S, 0) (3.2b)

lim S!0P (S, v, ø ) = Ke

°rø (3.2c)

lim S!1P (S, v, ø ) = 0 (3.2d)

P (S, 0, ø ) = max(Ke°rø ° S, 0) (3.2e)

lim v!1P (S, v, ø ) = Ke

°rø (3.2f)

where ø = T ° t.

The boundary conditions in the S direction are the standard conditions for European put

options (see e.g Wilmott [15]). The boundary condition at v = 0 is based on the riskless asset

growing argument and the boundary condition as v ! 1 is explained by Zhu and Chen [20]

and is based on using the BS value for the put in the limit as æ = p

v ! 1 i.e Ke°rø .

Using Lie’s classical symmetries method, we find PDE (3.2a) has a finite-dimensional

9

symmetry with generator

G = F (v, ø )P @

@P + T̄ (ø )

@

@ø ° vT̄ 0(ø )

@

@v ° Sg(ø )(1 +

2Ω

b )

@

@S

where F (v, ø ), T̄ (ø ) and g(ø ) satisfy the determining equations listed in Appendix A. From

these determining equations and with consideration of the given initial and boundary condi-

tions we use the following admitted symmetry generator:

G = °f (ø )P @

@P + T̄ (ø )

@

@ø ° vT̄ 0(ø )

@

@v ° Sf (ø )

@

@S

where T̄ (ø ) = 1 r

f (ø ) and

f (ø ) = ek R

ø

0 a(u)du °

a(0)

a(ø ) ° a(0)ek

R

ø

0 a(u)du

Z

ø

0

µ

a0(u)

a(u)2 e°k

R

u

0 a(u1)du1

∂

du. (3.3)

The corresponding invariant solution is

P = ¡(x, y)e°rø , x = Serø , y = vT̄ (ø ). (3.4)

Substitution of this into (3.2a) we find that ¡(x, y) needs to satisfy the reduced equation

x2¡ xx

+ 2Ωbxy¡ xy

+ b2y2¡ yy

° 2k(qa0 + y)¡y = 0 (3.5)

where a0 = a(0), q = 1 r

and as T̄ (0) = 0, from (3.2b-f) the initial and boundary conditions

for ¡ are given as ¡(x, 0) = max(K°x, 0), ¡(0, y) = K, lim x!1 ¡(x, y) = 0, limy!1 ¡(x, y) =

K.

Letting x = exp(X), y = exp(bY ), (3.5) becomes

¡ XX

+ ¡ Y Y

+ 2Ω¡ XY

° ¡ X

° b¡ Y

° 2k

b ¡

Y

° 2kqa0e

°bY

b ¡

Y

= 0 (3.6)

to be solved on the infinite domain subject to

10

lim Y !°1

¡(X, Y ) = max(K ° eX , 0) (3.7)

and that ¡ ! 0 as X ! 1. Taking the generalised Fourier Transform2 of (3.6) and (3.7) with

respect to X where F{¡(X, Y )} = F (w, Y ) we get that for Im(w) < 0,

@2F

@Y 2 °

@F

@Y

∑

b + 2Ωiw + 2k

b +

2kqa0e °bY

b

∏

+ (iw ° w2)F = 0 (3.8a)

lim Y !°1

F (w, Y ) = °K1+iw

w2 ° iw . (3.8b)

The solution to the above (see e.g Polyanin and Zaitsev [21]) is

F (w, Y ) = (3.9a)

°e°mbY °(m + 1 ° c1

b

)

°(2m + 1 ° c1 b

)

µ

2kqa0 b2

∂

m

K1+iw

w2 ° iw M

µ

m, 2m + 1 ° c1 b

, °2kqa0e°bY

b2

∂

where

c1 = °b ° 2Ωiw ° 2k

b and m =

c1 + p

c21 ° 4wi + 4w2 2b

. (3.9b)

Taking the Fourier inverse of (3.9a) then gives

¡(X, Y ) = °1 2º

Z

iw

i

+1

iw

i

°1

∑

e°iwX e°mbY K1+iw

w2 ° iw °(m + 1 ° c1

b

)

°(2m + 1 ° c1 b

)

µ

2kqa0 b2

∂

m

M

µ

m, 2m + 1 ° c1 b

, °2ka0qe°bY

b2

∂∏

dw

2Given that k1, k2 2 R, k1 < k2 and R 1 °1 e

°k1x|g(x)|dx < 1, R 1 °1 e

°k2x|g(x)|dx < 1, then the generalised Fourier transform Fg(z) =

R 1 °1 e

izx

g(x)dx, z 2 C exists and is analytic for all z in the strip {z 2 C : k1 < Im(z) < k2}. The inversion formula is given by g(x) =

R

iu+1 iu°1 e

°izxFg(z)dz with k1 < u < k2.

11

and hence with X = ln x, Y = 1 b

ln y, q = 1 r

we get

¡(x, y) = °1 2º

Z

iw

i

+1

iw

i

°1

∑

x°iwy°m K1+iw

w2 ° iw °(m + 1 ° c1

b

)

°(2m + 1 ° c1 b

)

µ

2ka0 rb2

∂

m

£M µ

m, 2m + 1 ° c1 b

, °2ka0 b2ry

∂∏

dw. (3.10)

The solution for the European put option price in its natural domain of definition, from

(3.4) is then given by P (S, v, ø ) = ¡(Serø , vT̄ (ø ))e°rø where T̄ (ø ) = 1 r

f (ø ) and f (ø ) is given

in (3.3) and ¡(x, y) is given in (3.10).

This solution agrees with the European put solution P1 in Lewis [12] when a(t) is constant

and ∞ = 1 so that (1.3a,b) (with µ replaced by r), corresponds to a risk-neutral process.

However, as stated and shown by Lewis [12], in practice we do the w°integration in

0 < Im(w) < 1 as within this strip the integrand is well-behaved with Re(w2 °iw) ∏ 0. To do

this, we use the put-call parity, P = Ke°rø ° [S ° C(S, v, ø )], where the expression in square

brackets represents the covered call option with payoÆ min(S T

, K). The payoÆ of the covered

call has a transform which is simply the negative of the payoÆ transform of the put option,

but with the restriction 0 < Im(w) < 1. This then leads to the solution (3.1). We note that

this solution satisfies (3.2b-f). ##

A comparison of put option prices when v = 0.2 using the BS formula and (3.1) with

a(ø ) = 0.3 and a(ø ) = 0.2[cos(4ºø ) + sin(4ºø )] + 0.4 is given in Figure 1. Other parameter

values used were K = 20, b = 2, T = 1, k = 10, Ω = °0.75, w̄ = 1000. Calculations were

performed using the mathematics computation package MAPLE [22] with the integration

performed numerically with a relative error tolerance of 0.5 £ 10°9. From the figure, we see

that in this example with the negative correlation, out-of-the-money put option prices with

the stochastic volatility models are higher than BS prices.

12

Figure 1: Put option values with T = 1, K = 20 using BS formula (bold line) and 3/2 stochastic volatility formula (3.1) with a(ø ) = 0.3 (dotted line) and a(ø ) = 0.2[cos(4ºø ) + sin(4ºø )] + 0.4 (dashed line).

4 Analytic Approximation

To find an approximation to (3.2a-f) for small ø we follow the method outlined by Howison

[23]. We let ø = ≤t0 where 0 < ≤ ø 1 and assume the solution can be expanded as a series

P (S, v, ø ) = 1

X

i=0

≤iP i

(S, v, t0). (4.1)

Substituting (4.1) into (3.2a) we get

° @P

@ø + ≤rS

@P

@S +

≤vS2

2

@2P

@S2 + ≤Ωbv2S

@2P

@S@v +

≤b2v3

2

@2P

@v2 + ≤kv(a(ø )°v)

@P

@v °≤rP = 0 (4.2)

13