Pricing options under the non-affine stochastic volatility models: An extension of the high-order compact numerical scheme

Accepted Manuscript

Pricing options under the non-affine stochastic volatility models: An extension of the high-order compact numerical scheme Guangping Shi , Xiaoxing Liu , Pan Tang PII: DOI: Reference:

S1544-6123(15)00136-1 10.1016/j.frl.2015.12.004 FRL 447

To appear in:

Finance Research Letters

Received date: Revised date: Accepted date:

14 September 2015 20 November 2015 8 December 2015

Please cite this article as: Guangping Shi , Xiaoxing Liu , Pan Tang , Pricing options under the nonaffine stochastic volatility models: An extension of the high-order compact numerical scheme, Finance Research Letters (2015), doi: 10.1016/j.frl.2015.12.004

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Title: Pricing options under the non-affine stochastic volatility models: An extension of the high-order compact numerical scheme

Author names and affiliations: Authors: Guangping Shi, Xiaoxing Liu, Pan Tang

CR IP T

Affiliations: School of Economics and Management, Southeast University Postal address: School of Economics and Management, Southeast University, Nanjing, Jiangsu, 211189, china

AN US

E-mail address: Guangping Shi ([email protected]) Xiaoxing Liu ([email protected]) Pan Tang ([email protected])

M

Corresponding author: Corresponding Author: Guangping Shi

AC

CE

PT

ED

Corresponding Author's Email: [email protected]

ACCEPTED MANUSCRIPT

Pricing options under the non-affine stochastic volatility models: An extension of the high-order compact numerical scheme Guangping Shi 1, Xiaoxing Liu, Pan Tang School of Economics and Management, Southeast University, Nanjing, Jiangsu, 211189, china

CR IP T

Highlights

AN US

1. We extend the high-order compact finite difference scheme to price options in non-affine stochastic models. 2. A proper variable substitution is introduced to convert the different nonlinear coefficients of second-order non-cross derivatives into an identical one; then, the high-order compact scheme for option pricing is structured. 3. Based on the local von Neumann stability analysis, a theoretical stability result is obtained under certain restrictions. 4. Numerical tests show that our scheme is convergent and effective on option pricing and that the non-affine parameter and volatility of volatility have important influences on option prices.

M

Abstract: We consider an improvement of a high-order compact finite difference scheme for option pricing in non-affine stochastic volatility models. Upon applying a proper transformation to

ED

equate the different coefficients of second-order non-cross derivatives, a high-order compact finite

PT

difference scheme is developed to solve the partial differential equation with nonlinear coefficients that the option values satisfied. Based on the local von Neumann stability analysis, a theoretical

CE

stability result is obtained under certain restrictions. Numerical experiments are presented showing

AC

the convergence and validity of the expansion methods and the important effects of the non-affine coefficient and volatility of volatility on option values.

Keywords: Non-affine stochastic volatility; Option pricing; High-order compact finite

difference method; Variable mixed derivatives; Nonlinear coefficients

Keywords and JEL codes 1

Corresponding author. E-mail address: [email protected]

ACCEPTED MANUSCRIPT

Keywords: Non-affine stochastic volatility; Option pricing; High-order compact finite difference method; Variable mixed derivatives; Nonlinear coefficients

M

AN US

CR IP T

JEL codes: C02 C52 G13

1. Introduction

ED

The correct specification of volatility dynamics is of crucial importance for option pricing. Most

PT

existing studies have analyzed the dynamics of volatility using an affine stochastic process, among which the Heston (1993) model is the most prominent work. Such models derive an analytically

CE

closed-form solution of the option price and have attracted greater attention for their advantage in

AC

explaining the time-varying, aggregation and leverage effect of the volatility of the asset return. Therefore, these models are also used in other areas, such as econometric estimation (Bates (2006)), term structure modeling (Dai and Singleton (2000)) and dynamic asset allocation (Liu (2007)). However, the affine stochastic volatility is relatively simple because of its square-root assumption, which fails in describing the nonlinear nature (Andersen et al., 2002) of a financial time series. Recently, many studies have focused on non-affine volatility models (Jones,

ACCEPTED MANUSCRIPT

2003; Kaeck and Alexander, 2012), which have been found to capture the dynamics of major equity indices much better than the affine models. Because the exact closed-form solution cannot be obtained in the non-affine stochastic volatility models or in the American option pricing problems, many numerical methods for option pricing

CR IP T

have been developed using the finite difference method (Hout and Foulon, 2010), element-finite volume method (Zhang et al., 2015) or series expansion (Hu and Kanniainen, 2015). Most of the finite difference methods for option pricing addresses the one-dimensional case and use the

AN US

standard, second-order finite difference, whereas Tangman et al. (2008) proposed a new high-order compact finite difference scheme for option pricing, that is only suitable for the constant volatility case. Hout and Foulon (2010) considered the numerical analysis for option

M

pricing in stochastic volatility; however, the final scheme there was of second-order solely due to the low order approximation of the cross diffusion term. Recently, a high-order compact finite

ED

difference scheme for option pricing in stochastic volatility models has been proposed by Düring

PT

and Fournié (2012), in which the stochastic volatility is considered to be affine models, which lead to an identical coefficient of second-order non-cross derivatives of the PDE. However, the

CE

non-affine stochastic volatility models complicate the problem, and we extend the approach

AC

proposed by Düring and Fournié to solve it. An outline of this paper is as follows. In section 2, we introduce the option pricing problems

under non-affine models. In section 3, we demonstrate the discretization of the pricing equation using a high-order compact scheme and analyze its stability. In Section 4, we present numerical experiments and results, and concluding remarks are provided in the last section.

2. Non-affine stochastic volatility models

ACCEPTED MANUSCRIPT

According to the non-affine stochastic volatility models (Chourdakis and Dotsis, 2011), the dynamics of an asset's price and stochastic volatility satisfy the following stochastic diffusions: dSt  rSt dt   t St dWt1 ,

(1)



d t  k (   t )dt  v t2 dWt 2 ,

(2)

CR IP T

where r , k , , v are the risk-free interest rate, the mean reversion speed, the long time mean of volatility, and the volatility of volatility, respectively, Wt1 ,Wt 2 are two standard Brownian motions with correlation   [1,1] , and  reveals the exponent in the diffusion function, which is named

AN US

as a non-affine parameter in our paper. It is obvious that the Heston model is one special case for   1 . For brevity, we only consider the case of European put option pricing problems, however,

the case of call options can be treated analogously. Let V ( S ,  , t ) denote the value of a European

M

put option with strike price K , application of Ito's lemma and standard arbitrage arguments leads to:

ED

 1  V 1  2V  2V 1  2V V V 2   S 2 2   vS 2    v2  rS  [ k (    )   v  ]  rV  0 , t 2 S  2 S  S  2

(3)

PT

for S ,   0 and 0  t  T , where  is the market price of risk. The terminal and boundary

CE

conditions to PDE (3) are provided as in Düring and Fournié (2012): V (S ,  , T )  max( K  S ,0) ,

lim V (S ,  , t )  0,limV (S ,  , t )  Ke r (T t ) ;

AC

S 

S 0

limV (S ,  , t )  0, lim V ( S ,  , t )  0;

 0

 

for

for

(4) T  t  0,   0 , T  t  0, S  0 .

(5) (6)

3. Extending the method of the high-order compact scheme 3.1 Derivation of the high-order compact scheme under non-affine models In this section, we extend the high-order compact method to price options under non-affine volatility models. A proper variable substitution is introduced to convert the different coefficients

ACCEPTED MANUSCRIPT

of second-order non-cross derivatives in Eq. (3) into an identical one, which is the key step to structure the high-order compact scheme for PDE. For this, we apply the following changes in variables x  ln  S K  ,  T  t , u  e r (V K ), y  2

1  (  3) 2

v(3   ) , and we obtain:

u  2u  2u  2u u u  a( y)( 2  2 )  b( y)  c( y )  d ( y )  0,  xy x y x y

(7)

2

2

CR IP T

for x  R, y  0,0    T , where 2

   1 1 1 1 1 a( y)   ( v(3   ) y)  3 , b( y )    ( v(3   ) y)  3 , c( y )  r  ( v(3   ) y)  3 , 2 2 2 2 2

 1

d ( y) 

 1

  v(  1) 1 k 1 k 1 1 ( v(3   ) y)  3  ( v(3   ) y)  3  ( v(3   ) y)   ( v(3   ) y)  3 . 4 2 v 2 v 2 2 1

AN US

Note that the above transformation, which can apply to a more general class of non-affine models, is an extension of the one used in During and Fournié (2012). The approach in our paper contrasts with that in the During and Fournié's work mainly because of the nonlinearity of the

M

coefficients a( y), b( y), c( y) and d ( y) ; it is obvious that by imposing restrictions  =1 , our method

ED

reduces to During and Fournié's. The problem is completed by the new initial and boundary conditions:

PT

u( x, y,0)  max(1  e x ) ; for x  R, y  0 ,

(8)

lim u( x, y, )  0; lim u( x, y, )  1 ; for y  0,0    T ,

(9)

lim u y ( x, y, )  0; lim u y ( x, y, )  0 ; for x  R, 0    T .

(10)

CE

x 

y 0

x 

y 

AC

To solve the numerical problem (7)-(10) at the point ( x, y , ) in the unbounded domain

(, )  (0, )  (0, T ) , we introduce xmax =ln5 and  max =1 according to Ikonen (2008), which

are sufficiently large to eliminate the boundary effect; then, the corresponding max value of y is 2/v(3- ) . With the step size in the x, y, directions being denoted by h1 , h2 ,  , respectively,

here, we set h1  h2  h to conveniently apply a high-order compact scheme; the step numbers are correspondingly N  [ xmax h]  1, M  [ ymax h]  1 , and L  [T  ]  1 . The approximate value

ACCEPTED MANUSCRIPT

of

u

at the grid point is thus denoted by

uin, j  u( xi , y j , n )  u(ih, jh, n ) , where

xi  ih, i  1, 2,..., N , yi  jh, j  1, 2,..., M ,  i  n , n  1, 2,..., L . First, we introduce the high-order

compact finite difference discretization for the two-dimensional equation after the reduction of the

 direction by introducing the notation f ( x, y)  u  : (11)

CR IP T

a( y)(uxx  u yy )  b( y)uxy  c( y)ux  d ( y )u y  f ( x, y ) ,

for {( x, y)  ( xmax , xmax )  (0, ymax )} . Using the first- and second-order central difference operator  x ,  y and  x2 ,  y2 , the discrete approximation to Eq. (11) at point ( xi , y j ) can be written as:

AN US

a( y)( x2ui , j   y2ui , j )  b( y) x y ui , j  c( y) x ui , j  d ( y) y ui , j  ei , j  fi , j ,

(12)

where the associated truncation error is given by: ei , j 

1 1 1 1 a( y)h2 (uxxxx  u yyyy )  b( y )h2 (uxyyy  uxxxy )  c( y )h 2u xxx  d ( y )h 2u yyy  o(h 4 ) . 12 12 6 6

(13)

M

Note that there are four fourth-order derivatives in Eq. (13), but three relevant relations could be derived from auxiliary Eq. (11). However, the derivative becomes achievable by regarding

ED

uxxxx  u yyyy and uxxxy  uxyyy as a whole, respectively. Differentiating Eq. (11), we obtain:

PT

uxxx  uxyy  b( y)uxxy +c( y)uxx +d ( y)uxy  f x  a( y ) , u yyy  uxxy  a( y)uxx + a( y)  d ( y)u yy +b( y)  c( y)uxy +b( y)u xyy +c( y)u x +d ( y )u y  f y  a( y) ,

(14) (15)

CE

uxxxy  uxyyy   a( y)uxxx  a( y)  d ( y) uxyy  [b( y)  c( y)]uxxy  b( y)uxxyy  c( y )u xx  d ( y )u xy  f xy  a( y ) , (16)

AC

a( y)(uxxxx  u yyyy )  b( y)(uxxxy  uxyyy )  c( y)uxxx  2a( y)  d ( y)u yyy  2a( y)uxxyy  d ( y)  2a ( y) u xxy  c( y)  2b( y) uxyy  a( y)uxx  a( y)  2d ( y)uxy  c( y)ux  d ( y)u y  f xx  f yy ,

where

a( y), b( y), c( y), d ( y)

(17)

and a( y), b( y), c( y), d ( y) are the first- and second-order

derivatives of a( y), b( y), c( y), d ( y) with respect to y . For simplicity in notation, we denote the coefficients a( y), b( y), c( y), d ( y) as a, b, c, d in the following. Substituting approximations (14)-(17) into (13), a new expression for the truncation error term can easily be deduced;

ACCEPTED MANUSCRIPT

inserting it into (12) yields the o(h4 ) approximation of PDE (11): h2  b2  h2  b a d    x2 y2 ui , j   2d   b  2c   x2 y ui , j  2a    12  a  y y 12  a a   y  y j

j



h    h  2b bc c  a b  a   2c  2b   d  a     x y2 ui , j    a     c     d  2a     a   x2ui , j  12  a a a a  a  y y   y y 12  2

2

j

j

h   h  a  d d    a   2d    d  2a  a   y2 ui , j   d    d  2a  d   y ui , j a a   y yj   y yj 12  12  2

2

 fi , j 

h2 2 1  a b  1 ( x   y2 ) fi , j   c   x f i , j   d  2a    y fi, j  12 a a  y  yi a y yj

CR IP T

 h2   h2    bd  d  a b  b  c c      b  2c   c   d  2 a  b   u  c    d  2a     c   x ui , j     x y i , j    a a a  a a   y yj   y yj  12  12 

(18)

AN US

The fourth-order compact finite difference scheme considered at the mesh point (i, j) involves the nearest eight neighboring mesh points. Associated with the shape of the computational stencil, the indexes for each node from zero to nine are introduced:

ED

M

 ui 1, j 1 ui , j 1 ui 1, j 1   u6    ui , j ui 1, j    u3  ui 1, j u    i 1, j 1 ui , j 1 ui 1, j 1   u7 8

Hence, scheme (18) is defined by: 2 h2 2 h2

1  b2  2 a    a  y y 6h 2  j

1  2 h

 2 j ,4 j   

lj

fl , l  0,1,...,8, j  1, 2..., M , where

 h2   a  d 1  b2   d  2a   a  + 2  2a   ,   a   2d   a a  y y   y  y j 3h  12  j

AC

1 j ,3 j  

l 0

(19)

 h2  bc c  a b  a     c    a      d  2a     a  a a a  a   y yj 12 

CE 

lj l

u5   u1  , u8 

PT

l 0

0 j  

8

  u  

u2 u0 u4

 h2   bc c  a b  a   1  h2  c     c   d  2 a  a      a    c   d  2a     c    a a a  a a   y yj  12   y  y j 2h 12 

1  b2  2a   2  a  y y 6h  j

1 h2

1  2b  2c  2b   d  a     12h  a  y yj

1  b ad   2d   b  2c   12h  a a   y  y

j

h    a  d 1 h d  d  2a   a    d    d  2a   d  ,   a   2d   12 a 2 h 12 a      y yj    y yj 2

2

,

ACCEPTED MANUSCRIPT

1  b2  1  b a d   1  2b  2a    2d   b  2c   2c  2b   d  a     2  a  y  y 24h  a a   y  y 24h  a 12h   y yj j j



6 j ,8 j 

1  h2  bd  d  a b  b  c   b  2c   c   d  2a     b    2  a a a  a 4h 12    y yj

1  b2  1  b a d   2 a   2d   b  2c      2 a  y  y 24h  a a   y  y 12h  j j



1 4h 2

,

1  2b  2c  2b   d  a     24h  a  y yj

 h2  bd  d  a b  b  c    c   d  2a     b   b  2c   a a a  a   y yj 12 

,

CR IP T

5 j ,7 j 

2 1 h  ab  1 h b 0 j  , 1 j ,3 j   , 2 j ,4 j   (d  2a) , 5 j ,7 j  c   3 12 24a  a  y y 12 24a 48 a y y j

j

where multiple indexes are used for simplicity with  and

y y j

b 48a

,

y y j

signs , and the first index

AN US

corresponds to the upper sign.

, 6 j ,8 j  

The high-order compact scheme of the original problem (7) can be obtained by extending the two-dimensional Eq. (11) considering f ( x, y)  u  . Applying the implicit difference to the 8

   u

n 1 lj l

l 0

+(1   )lj uln   lj  uln 1  uln   ,

M

time derivative u   (uln 1  uln )  yields:

8

l 0

ED

where uln denotes the nine point uin, j in the computational stencil, and   0,1 2,1 represents the forward Euler, Crank-Nicolson and backward Euler schemes, respectively. Thus, the resulting

CE

PT

discrete scheme for node (i, j ) at the nth time can be rewritten as: 8

 u l 0

n 1 lj l

8

  lj uln ,

(20)

l 0

where the coefficients  lj and lj are given by: lj  lj 

lj  , lj  (1   )lj  lj .  

AC

3.2 Stability analysis The proof of the stability is not trivial because our scheme is applied to a problem with

nonlinear variable coefficients, which is much more complex than the linear one in During and Fournié (2012). However, the principle of ''frozen coefficients'' and the von Neumann stability analysis are also used to discuss the stability of our scheme. If each frozen coefficient problem is stable, the variable coefficient problem is also stable. We now rewrite uin, j as

ACCEPTED MANUSCRIPT

uin, j  g n eIiz1  Ijz2 ,

(21)

where I is the imaginary unit, z1 , z2  [0,  ] . Then, the scheme is stable if, for all z1 and z2 , the amplification factor G  g n 1 g n satisfies the relation G 1  0 . 2

(22)

unachievable, however, we can show the following result.

CR IP T

As During and Fournié (2012) pointed out, the complete stability analysis is currently

Theorem 1. For r      0 and   1 2 , scheme (20) satisfies the stability condition (22) under the restrictions that   1 and the step size h is sufficiently small.

( A  B)  (C  D) I , where ( A  B)  (C +D) I

AN US

Proof. Using (21) in (20), an expression for G can be obtained, G =

A  0, j  1, j  3, j  cos z1  2, j  4, j  cos z2  5, j  7, j  cos( z1  z2 )  6, j  8, j  cos( z2  z1 )   , B  0, j  1, j  3, j  cos z1  2, j  4, j  cos z2  5, j  7, j  cos( z1  z2 )  6, j  8, j  cos( z2  z1 )  2 ,

M

C  1, j  3, j  sin z1  2, j  4, j  sin z2  5, j  7, j  sin( z1  z2 )  6, j  8, j  sin( z2  z1 )   , D  1, j  3, j  sin z1  2, j  4, j  sin z2  5, j  7, j  sin( z1 +z2 )  6, j  8, j  sin( z2  z1 )  2 .

ED





Moreover, (22) can be written as 4( AB  CD) ( A  B)2  (C  D)2 

2

( A  B) 2  (C  D) 2   0 ,

PT

which is equivalent to AB  CD  0 . After calculation and simplification, we found 1 h h (4  cos z1  cos z2 ) , C =  sin z1  (d  a)sin z2 , 6 12 12a

B=

z  a z  1  v 2 (  1) a  v 2 (  1) v 2a  (1  cos z )  (1  cos z2 )+2sin 2 1   2  (d  2a )  2 sin 2 2  1   2 2 48  4a 3 2  h 48 2 192a 3h

CE

A=

AC

 2sin 2

z2  a 1  z z z  a  d  d  2   2d   (d  2a )   sin 1 sin 1 cot 2  ,  2  h 12  a 2 2 2  3

h hv h hd  d a  a  d  D=    a  (d  2a)  sin z1    d   (d  2a)  sin z2  cos z2 sin z1 . 3 h 24 96 3 h 24 24 a 3h    

It is obvious that A  0 , A  C . For B , note that the first two items are not less than zero when   1 ; and that the latter two items are nonnegative when the step size h is sufficiently small. Hence, B  0 when   1 and h is sufficiently small. Moreover, we verify that

ACCEPTED MANUSCRIPT

B  v 2 (  1) 96a 2 ,

D  h(d   a) 24  h(d   2a)(av  4d ) 96a 2 ; it is easy to observe that

B  D when   1 and h is sufficiently small. Under the specific conditions in Theorem 1,

we have AB  0 and AB  CD ; therefore, AB  CD  0 , which completes the proof. □ For the more involved situations that the correlation  and the long time mean of volatility 

CR IP T

are non-zero, we reference the lemma 2 in During and Fournié (2012), which is not repeated here. We found that one of the sufficient conditions in lemma 2 that G  1=0 for y =0 is true; 2

consequently, we conjecture that the stability condition (22) is also satisfied for non-zero  and

AN US

 , although it will be difficult to provide an analytical proof. Therefore, we perform additional

numerical tests to validate the stability property of scheme (20) in the next section.

4. Numerical experiments

M

In this section, numerical experiments are performed to test the effectiveness of the scheme (20). According to During and Fournié (2012) and Zhu and Chen (2011), the parameters are set as

ED

follows:   1 2, r  0.1,   0.5, k  2.5,  0.16, v  0.45, T  0.5, K  $10 . We set  to zero to

PT

streamline, although our scheme applies to the case that  is non-zero by extending the coefficient straightforward. Moreover, for the non-affine parameter  , the commonly used values

CE

in previous literature are  =1, 2,3 , and others (such as in Ballestra and Cecere, 2015) that are no

AC

larger than 3. Furthermore, Ignatieva (2009) and Kaeck and Alexander (2012) indicated that moving from the square root diffusion towards 3/2 diffusion allows to capture the behavior of outliners at tails. Therefore, it is reasonable to set   3 , which makes y positive. For comparison, we calculate four sets of European put options with different non-affine parameters

  0.5,1,1.5, 2.5 . For the step size, h needs to be sufficiently small to meet stability conditions without restrictions on the time step size, and h and  are set to 0.1 and 0.01, respectively.

ACCEPTED MANUSCRIPT

Fig. 1 shows that h=0.1 is sufficiently small to guarantee the stability of our scheme. However, the monotonicity of the solution is not granted in Fig. 1. (a) because   1 , whereas the stability appears suitable in Fig. 1. (b), (c) and (d) with   1 , which is exactly consistent with the theoretical result. From Fig. 1, we can also find that the option prices increase with  but not

CR IP T

obviously, the specific option prices are shown in Table 1. This table shows that the European option prices increase with the non-affine parameter, which indicates the importance of the non-affine parameter for option pricing and could not be ignored. (a) European put option value with   0.5

AN US

(c) European put option value with   1.5

(b) European put option value with   1

(d) European put option value with   2.5

Fig. 1. European put options prices with different  values Comparison of the computed option prices with different  and 



Volatility

0.2

9

10

11

12

20

5.0972

2.9480

2.7669

2.5883

2.4128

2.2404

1.4244

5.4609

2.9674

2.7851

2.6061

2.4303

2.2576

1.4400

5.7379

2.9868

2.8041

2.6246

2.4484

2.2753

1.4561

1.5

6.4686

3.0899

2.9041

2.7219

2.5432

2.368

1.5393

2

6.9895

3.3432

3.1468

2.9555

2.7687

2.5863

1.7275

2.5

7.3914

4.2923

4.0737

3.8582

3.6458

3.4364

2.4192

0.5

6.2593

3.0476

2.8623

2.6822

2.5046

2.3303

1.5057

0.8

6.4686

3.0899

2.9041

2.7219

2.5432

2.368

1.5393

1

6.6206

3.1334

2.9460

2.7625

2.5826

2.4062

1.5732

1.5

6.9334

3.2949

3.1007

2.9112

2.7261

2.5452

1.6930

2

7.1871

3.6181

3.4108

3.2901

3.0127

2.8211

1.9195

2.5

7.4476

4.6237

4.4078

4.1926

3.9783

3.7649

2.6985

0.5

6.6206

3.3133

2.9459

2.7623

2.5825

2.4061

1.5731

0.5 0.8

CE

PT

1

AC

0.5

0.7

Asset price

8

ED

2

M

Table 1

ACCEPTED MANUSCRIPT 3.1782

2.9889

2.8039

2.6227

2.4451

1.6072

1

6.8662

3.2474

3.0552

2.8676

2.6842

2.5047

1.6585

1.5

7.0586

3.4169

3.2174

3.0231

2.8338

2.6489

1.7795

2

7.2523

3.7694

3.5575

3.3511

3.1498

2.9534

2.2058

2.5

7.4637

4.7366

4.5226

4.3086

4.0946

3.8807

2.7997

0.5

6.9119

3.2642

3.0714

2.8832

2.6992

2.5193

1.6713

0.8

6.9872

3.3306

3.1347

2.9439

2.7576

2.5756

1.7187

1

7.0353

3.3800

3.182

2.9892

2.8011

2.6175

1.7535

1.5

7.1651

3.5640

3.3584

3.1586

2.9640

2.7742

1.8815

2

7.3084

3.9335

3.7187

3.5071

3.3013

3.0999

2.1433

2.5

7.4809

4.8608

4.6494

4.4371

4.2239

4.0098

2.9139

CR IP T

6.7362

AN US

1

0.8

Furthermore, denote error n  uin, j (2h, )  u2ni ,2 j (h, ), E (h, )  max error n , order  log 2 

1 n T / 

E (2h, ) ; E (h, )

we should exam the ratios of the convergence errors with grid steps along spatial directions being

M

successively decreased, whereas the grid step along the time direction is fixed. For this target, we fix the time grid step to be   0.01 and vary the sizes of the spatial step from 0.4 to 0.05. Table

ED

2 shows that the convergence orders of the scheme in the spatial direction are approaching to 4,

PT

hence approximately o(h4 ) , and errors approach 0 as h  0 . Table 2

h

CE

 =0.5

Convergence orders in the space variable with different 

 =1.5

 =1

values

 =2

E (h,  )

Order

E (h,  )

Order

E (h,  )

Order

E (h,  )

Order

0.000853

-

0.000926

-

0.001006

-

0.001192

-

0.1

0.01343

3.9471

0.01482

3.8799

0.01481

3.8127

0.01675

3.9768

0.2

0.2111

3.8836

0.2115

3.8367

0.2116

3.6591

0.2116

3.9744

0.4

2.2328

3.158

1.9236

3.1899

1.931

3.3792

2.2017

3.4029

AC

0.05

Numerical analysis is also performed to obtain insight into the influences of the non-affine parameter and the volatility of volatility on European put option prices. Fig. 2 shows that the European option price is an increasing function of the volatility  , and it increases faster at the

ACCEPTED MANUSCRIPT

low volatility level than at the high volatility level. Furthermore, Fig. 2 displays that the non-affine parameter influences the option price by impacting the volatility process, when  takes a larger value (for example,   2 ), the influence of volatility on the option price is much more obvious than when  has a smaller value (for example,   1 ). Fig. 3 illustrates the European option

CR IP T

prices under different volatility of volatility v and different volatility level  . Every subgraph of Fig. 3 displays that the European option price is a decreasing function of v . One more result from the comparison of the two subgraphs in Fig. 3 is that the option prices changes more with as

AN US

regard to the volatility of volatility difference when the volatility takes a higher value,

observed in Fig. 3. (b); however, the changes is less when the value of volatility is lower, as observed in Fig. 3. (a).

M

(a) option values with different S at   1

(b) option values with different S at   2

Figure. 2. The impact of volatility and the non-affine coefficient on the European put options prices (b) option values with different v at   0.7056

ED

(a) option values with different v at   0.2304

PT

Figure. 3. European put options prices vs. underlying asset price with different volatilities of volatility

CE

5. Conclusion

This paper presents an improvement of a high-order compact finite difference scheme for the

AC

European option pricing problem under the non-affine stochastic volatility model. The key features of the current scheme are its high efficiency with fourth-order accuracy in space and second-order accuracy in time for the PDE with nonlinear coefficients and inequality coefficients of the second-order non-cross derivatives. Numerical examples demonstrate the efficiency of the proposed scheme. This scheme could be extended to other stochastic volatility models and the American option pricing problem. For the first extension, an appropriate transformation would be

ACCEPTED MANUSCRIPT

introduced to equal the different coefficients of second-order non-cross derivatives. In the second case, one must take a front-fixing transformation and combine the high-order discretization with a high-order resolution of the optimal exercise price. We leave extensions for future research.

Acknowledgments

CR IP T

The authors are highly grateful to any anonymous referee for their careful reading and insightful comments, and the views and opinions expressed are those of the authors. The work is supported by the National Nature Science Foundation of China (Nos.71273048, 71473036,71401031).

AN US

References

Andersen, T.G., Benzoni, L., Lund, J., 2002. Estimating jump-diffusion for equity returns. J. Finan. 57, 1239-1284.

M

Bates, D., 2006. Maximum likelihood estimation of latent affine processes. Rev. Financ. Stud. 19,

ED

909-965.

Ballestra, L.V., Cecere, L., 2015. Pricing American options under the constant elasticity of

CE

45-55.

PT

variance model: An extension of the method by Barone-Adesi and Whaley. Finance. Res. Lett. 14,

Chourdakis, K., Dotsis, G., 2011. Maximum likelihood estimation of non-affine volatility

AC

processes. J. Empirical Finance. 18, 533-545. Dai, Q., Singleton, K.J., 2000. Specification analysis of affine term structure models. J. Finan. 55, 1943-1978. Düring, B., Fournié, M., 2012. High-order compact finite difference scheme for option pricing in stochastic volatility models. J.Comput. Appl. Math. 236, 4462-4473. Heston, S.L., 1993. A closed-form solution for options with stochastic volatility with application

ACCEPTED MANUSCRIPT

to bound and currency options. Rev. Financ. Stud. 6, 327-343. Hu, J., Kanniainen, J., 2015. Asymptotic expansion of European options with mean-reverting stochastic volatility dynamics. Finance. Res. Lett. 14, 1-10. In't Hout, K.J., Foulon, S., 2010. ADI finite difference schemes for option pricing in the Heston

CR IP T

model with correlation. Int. J. Numer. Anal. Model. 7 (2), 303-320. Ikonen, S., Toivanen, J., 2008. Effcient numerical methods for pricing American options under stochastic volatility. Numer. Methods Partial Differential Equations. 24, 104-126.

AN US

Ignatieva, K., Rodrigues, P., Seeger, N., 2009. Stochastic volatility and jumps: exponentially affine yes or no? An empirical analysis of S&P 500 dynamics. Working Paper, Goethe University Frankfurt.

markets. J. Econ. 116, 181-224.

M

Jones, C.S., 2003. The dynamics of stochastic volatility: evidence from underlying and options

ED

Kaeck, A., Alexander, C., 2012. Volatility dynamics for the S&P 500: Further evidence from

PT

non-affine,multi-factor jump diffusions. J. Bank. Finan. 36, 3110-3121. Liu J., 2007. Portfolio selection in stochastic environments. Rev. Finance. Stud. 20, 1-39.

CE

Zhu, S.P., Chen, W.T., 2011. A predictor-corrector scheme based on the ADI method for pricing

AC

American puts with stochastic volatility. Comput. Math. Appl. 62, 1-26. Tangman, D.Y., Gopaul, A., Bhuruth, M., 2008. Numerical pricing of options using high-order compact finite difference schemes. J. Comput. Appl. Math. 218, 270-280. Zhang, R., Zhang, Q., Son, H., 2015. An efficient finite element method for pricing American multi-asset put options. Commun. Nonlinear. Sci. Numer. Simulat. 19, 25-36.

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT