Implied volatility expansion under the generalized Heston model

School of Education, Culture and Communication

Division of Applied Mathematics

BACHELOR THESIS IN MATHEMATICS / APPLIED MATHEMATICS

Implied volatility expansions under the generalized Heston model

by

Hanna Andersson and Ying Wang

School of Education, Culture and Communication

Division of Applied Mathematics

Bachelor thesis in mathematics / applied mathematics

Date:

2020-05-29

Project name:

Implied volatility expansions under the generalized Heston model

Author(s):

Hanna Andersson and Ying Wang

Supervisor(s): Ying Ni Reviewer: Anatoliy Malyarenko Examiner: Linus Carlsson Comprising: 15 ECTS credits

Abstract

In this thesis, we derive a closed-form approximation to the implied volatility for a European option, assuming that the underlying asset follows the generalized Heston model. A new para-meter is added to the Heston model which constructed the generalized Heston model. Based on the results in Lorig, Pagliarani and Pascucci [11], we obtain implied volatility expansions up to third-order. We conduct numerical studies to check the accuracy of our expansions. More specifically we compare the implied volatilities computed using our expansions to the results by Monte Carlo simulation method. Our numerical results show that the third-order implied volatility expansion provides a very good approximation to the true value.

Keyword: Heston model; Generalized Heston model; implied volatility; implied volatility expansion; Black–Scholes; Monte Carlo method; European options

Acknowledgements

First of all we would like to express our deepest gratitude to our supervisor Dr. Ying Ni for in-troducing to, and guiding us through this project. We are thankful for her patient, the valuable advice and she is always available with her consistent support when we come across obstacles along the project.

Our special appreciation to our reviewer Professor Anatoliy Malyarenko for taking his time, and giving his helpful and professional comments to improve the quality of our thesis.

Our great gratitude to our examiner Dr. Linus Carlsson for his insightful suggestions for fur-ther improvements for this thesis.

Last but not least we would like give our sincerely thanks to our fellow students Omar Mo-hammad and Rafi Khaliqi for sharing great advice concerning the numerical analysis. We would also like to show gratitude to each other for always keeping up the good mood and supporting each other during this special circumstances.

Contents

1 Introduction 4

1.1 Motivation . . . 4

1.2 Problem Formulation and Contribution . . . 6

1.3 Methodology . . . 6

1.4 Disposition . . . 7

2 Theoretical Consideration 8 2.1 European Options . . . 8

2.2 General Local Stochastic Volatility Model . . . 9

3 Analytical Results 12 3.1 Generalized Heston Model . . . 12

3.2 The Implied Volatility Expansions . . . 15

4 Numerical Results 22 4.1 Procedure . . . 22

4.2 Third-order Implied Volatility Expansion vs Monte Carlo Approximation . . 24

4.3 Comparison Between Different Orders of The Implied Volatility Expansions . 25 4.4 Summary . . . 27 5 Conclusion 28 5.1 Conclusion . . . 28 5.2 Further Research . . . 28 A Proofs 29 B Matematical Derivations 33 C MATLAB Codes 35

List of Figures

4.1 A comparison among the different orders of the expansion and the Monte Carlo approximation with different maturities. . . 26

List of Tables

4.1 A comparison of the implied volatility between our third-order expansion and Monte Carlo approximation with respect to different strike prices. . . 24 4.2 A comparison of the implied volatility between our third-order expansion and

Monte Carlo approximation with respect to different maturities. . . 24 4.3 A comparison of the implied volatility between our third-order expansion and

Chapter 1

Introduction

1.1

Motivation

In financial market, an option is a contract which gives the holders of the contract the right, but not the obligation, to buy or sell an underlying asset by a certain date (or earlier, depending on option style) and price [2, 6]. There are two types of options, namely the call option and the put option[6]. These two types of options give the holders the opposite positions and rights on the contract. The call option gives the holder the right, but not the obligation to buy the underlying asset by a certain date and price, and on the other hand, the put option gives the holder the right, but not the obligation to sell the underlying asset by a certain date and price [6]. The date and the price in the contract are known as the maturity and the strike price. There are two major styles of the options, which are the American options and the European options. American options can be exercised at any time up to the maturity, and European options can be only exercised at the maturity [6].

The revolution of trading and pricing the derivative securities started in the beginning of 1970’s [17]. Fischer Black, Myron Scholes, and Robert Merton achieved to calculate the price of European stock options in 1973 [3, 15]. At same year, the Chicago Board of Option Exchange started the trading of options in exchanges [17]. The options have become one of the most important instruments in the financial market [9]. The model which Black et al. (1973) developed is well-known as the Black–Scholes (or Black–Scholes–Merton) model. It has had a huge influence in the pricing and hedging derivatives [6]. All parameters which are used in the Black–Scholes pricing formulas can be directly observed in the market except the volatilityof the stock price [6, 17]. One can only observe the historical volatilities but not the implied volatilitiessince the volatility can be volatile at any time moment. It is unfortunately not possible to invert the Black–Scholes formulas to derive the expression of the volatility as a function of the known parameters which are observed in the market [6]. Moreover, there are

some market stylized facts such as volatility smiles1 and skews2 which can not be explained under the Black–Scholes model. Based on these factors, a lot of extensions have been derived to model the implied volatilities. One reasonable direction is making the volatility stochastic [1, 4].

The Heston model is a popular stochastic volatility (SV) model in mathematical finance [8, 13]. It was presented in 1993 by Steven Heston, who derived a closed-solution for the price of a European call option on an asset with stochastic volatilities [16]. The arbitrary correlation between volatility and spot-asset returns3are allowed in this model [16]. One popular exten-sion to the Heston model is a generalized Heston model which uses a mean-reverting constant elasticity of variance (CEV) process as volatility [8], i.e.

dSt = bStdt + √ v_tS_t(ρdWt2+ q 1 − ρ2_dW1 t ) dV_t = α(θ −Vt)dt + βVtξdWt2 with ξ ∈ (1/2, 1).

Many SV models have the undesirable property that moments of order higher than one can explode in finite time [7]. The authors of [7] also examined the moment explosion property across a spectrum of the SV models which includes the moment stability in the Heston model. They demonstrated that the negative correlation between the processes can provide a super-linear stabilizing effect in the variance process, which can prevent a moment explosion. Thus, the correlation coefficient ρ in this project is set to be negative.

In the research of Altmayer and Neuenkirch (2015), they established an integration by parts formula for the quadrature of discontinuous payoffs in a multidimensional Heston model un-der the assumption of a strictly positive volatility [8]. A variance process for the generalized Heston model is presented with the constant elasticity parameter ξ ∈ (1/2, 1) [8, 13] as men-tioned in previous section. The model of interest in this project is constructed partially based on their variance process.

Lorig et al. (2014) derived a family of asymptotic expansions for European-style option prices and implied volatilities for an asset whose risk-neutral dynamics are described by a general class of local stochastic volatility (LSV) models [11]. Their expansions are explicit, hence there are no requirements for any special functions nor numerical integration. The explicit im-plied volatility expansions have been obtained for the imim-plied volatility under following spe-cific LSV models: constant elasticity of variance (CEV), stochastic, alpha, beta, rho (SABR),

model and SV model. In LV models, the volatility is a function of time, while in the SV mod-els the variance process is stochastic. This project is an extension of the paper by Lorig et al. (2014). Our study is based on some of their results to derive the closed-form approximation to the implied volatility for a European option.

1.2

Problem Formulation and Contribution

Our task is to derive a third-order expansion for the implied volatility, with consideration of an asset whose risk-neutral dynamics are described by a generalized Heston model. Based on this main purpose, we apply the general results in [11] to derive explicit implied volatility expansions for our model. Hence, the primarily focus is on the following two problems,

1. Using the formula which given in [11] to derive the third-order expansion of the implied volatility under the generalized Heston model.

2. Checking the accuracy of our expansions using Monte-Carlo Simulation.

As mentioned in the motivation, Lorig et al. (2014) presented implied volatility expansions for five different LSV models. We provide another additional model where the implied volatility expansion for an asset whose risk-neutral dynamics are described by a generalized Heston model. As calibration of a model for the underlying process often involves calibrating to the market quoted/induced implied volatilities. Such an explicit and easy-to-use approximation formula for implied volatility is very useful for the calibration purpose. Our generalized model has one extra parameter therefore it is more flexible to fit the real market data than the standard Heston model. Our expansion is also more efficient since it has closed-form solution which gives a value directly, not only convenient for evaluating implied volatility but also suitable for determining option price. Because of its efficiency it is much more competitive than other approximation methods e.g. Monte Carlo simulation which is relative more time consuming.

1.3

Methodology

Since this thesis is an extension of the paper of Lorig et al. (2014), we applied some results which are provided in their work to obtain our implied volatility expansions. Firstly, a new ad-ditional variable was added to the Heston model to obtain the generalized model [8, 13]. Then the expressions for the implied volatility were achieved by deriving the expansions according to procedure process.

To check the accuracy of the expansion we used Monte Carlo methods to simulate the price for a European call option, then we found the implied volatility. Furthermore, we compared our expansions with this simulated implied volatility. The software which was used for the numerical analyses was MATLAB.

1.4

Disposition

The thesis is structured as follows:

In Chapter 1 the topic is introduced together with the motivation and previous research. Chapter 2 introduces the theoretical work together with the general LSV model. In Chapter 3 we show how we modified the generalized Heston model and the derivations of the implied volatility expansions. In Chapter 4 we provide numerical analysis to show that our model is appropriate. This thesis ends with Chapter 5 where a conclusion takes place, together with suggestions to further research. Appendix A provides the proofs and the formulas, Appendix B provides detailed mathematical derivations and Appendix C contains the MATLAB codes.

Chapter 2

Theoretical Consideration

2.1

European Options

In this thesis we consider only one type of options, namely the European call options. As mentioned in the motivation part, the European options can be only exercised at the maturity. The payoff function from a long position (also called call) of a European option is given by

max(S − K, 0),

where S is the price of the underlying assets at the maturity, and K is the strike price [6]. In this project, we mainly follow [11] that assume a general Heston model for the logarithm of the underlying asset X = log S and the logarithm of the strike price k = log K. For simplicity, we also assume a market with no arbitrage opportunities, no dividends, no transactions costs and zero interest rate. We fix a time to maturity t > 0, and an initial value X0= x, thus the

European call option payoff becomes H(Xt) = max(eXt− ek, 0) where Xt is the final price of

the underlying assets at the maturity.

We also consider moneyness (K/S) for a European call option in this project. The moneyness is used when one consider an option, whether it is profitable to exercise at maturity or not. There are three different states for the moneyness, namely at the money (ATM), in the money (ITM), and out of the money (OTM). ATM appears when stock price equals to the strike price, ITM is when stock price is greater than the strike price, and OTM is the opposite of ITM. Since we are applying the Heston model and it includes change of variables, mentioned above, we are considering the log-moneyness (k − X ). In simple moneyness, ATM is when the value of moneyness is 1, ITM corresponds to moneyness greater than 1 while OTM is when it is less than 1. In the log-moneyness contraiwise, ATM corresponds to moneyness equal to 0, ITM to greater than 0, and OTM to less than 0.

2.2

General Local Stochastic Volatility Model

As mentioned in previous section, we consider European options with strike price K and ma-turity time t on a non-dividend-paying stock with no transaction costs, zero interest rates and an arbitrage free market. We consider an equivalent martingale measure P on a complete filtered probability space (Ω, {Ft,t ≥ 0}, P) which is chosen by the market [11]. The filtration

{Ft,t ≥ 0} is given as the history of the market. All stochastic processes which are defined

below live on this probability space and all expectations are taken with respect to P. Consider the following general form of a LSV model for the strictly positive underlying process,

St= eXt, dXt= − 1 2σ 2_(X t,Yt)dt + σ (Xt,Yt)dWt, X0= x ∈ R, dYt= α(Xt,Yt)dt + β (Xt,Yt)dBt, Y0= y ∈ R, dhW, Bi_t= ρdt, |ρ| < 1. (2.1) The results in [11] can be applied to the instantaneous coefficient correlation function ρ(X ,Y ). Here for simplicity we consider only constant ρ. In [11], the LSV model above is studied. They provided a general procedure for obtaining option prices expansions and the corres-ponding implied volatility expansions.

Definition 1 (The infinitesimal generator of the process (X ,Y )). A (t) = a(t,x,y)(∂2

x − ∂x) + f (t, x, y)∂y+ b(t, x, y)∂y2+ c(t, x, y)∂x∂y, (x, y) ∈ R2, (2.2)

where, a:= σ 2 1 2 , b:= σ₂2 2 , c:= ρσ1σ2, f := µ2. (2.3) and σ1 and σ2 are the diffusion rates of each process i.e. σ1= σ (Xt,Yt) and σ2= β (Xt,Yt)

and µ2is the drift rate of the second process i.e. µ2= α(Xt,Yt)

As in [11], we assume that the stochastic differential equation (SDE) (2.1) has a classical unique solution, the coefficients of the diffusion terms are strictly positive functions and all coefficients in the system of SDE’s are smooth. The equation (2.1) includes basically all one-factor SV models, all LSV models and all one-factor LSV models. Based on those as-sumptions we investigate the generalized Heston model, which is a popular type of one-factor SV model.

Let u(t, x, y) be the European option price at time t with initial values (x, y) = (X0,Y0) := (x, y).

Definition 2 (Black–Scholes price uBS). For fixed initial stock price S0, strike price K and

time to maturity t, the Black-Scholes price uBSis given by

uBS(σ ) := S0N (d1(σ )) − KN (d2(σ )), (2.4)

whereN is the cumulative distribution function of a standard normal variable, and d1(σ ) := 1 σ √ t(log( S0 K) + σ2t 2 ), d2(σ ) := d1− σ √ t. (2.5) Remark1. We consider the log initial stock price x = log S0and the log strike price k = log K

in this project, hence uBS is given by

uBS(σ ) := exN (d1(σ )) − ekN (d2(σ )), (2.6) where d₁(σ ) := 1 σ √ t(x − k + σ2t 2 ), d2(σ ) := d1− σ √ t. (2.7) Definition 3 (Implied volatility [11]). For fixed initial stock price S0, strike price K and time

to maturity t, the implied volatility corresponding to a call price u on its domain

(max(S0− K), S0) is defined as the unique strictly positive real solution σ to the equation

uBS(σ ) = u. (2.8) Definition 4 (η–function). We define the η–function from [11] as below

ηi, j=

∂_xi∂yjη (x, y)

i! j! , η ∈ {a, b, c, f }. (2.9) where a, b, c and f refers to the functions in (2.3). For example

ai, j=

∂_xi∂yja(x, y)

i! j! , where i, j ∈ N

We summarize in Proposition 1 the most relevant results from [11] that are implemented in our study. Note that Proposition 1 holds under some technical conditions on the infinitesimal generatorA (2.2). These conditions are given in rigorous form in [11] and further released in [12]. At this stage of study, we assume that our generalized Heston model satisfies the required conditions and we derive the expansions using Proposition 11. These expansions are then used as approximations to the implied volatility. We conduct numerical experimental studies to investigate the quality of these approximations.

1_{We are confident that the generalized Heston model indeed satisfy the technical conditions. A full discussion}

Proposition 1. In [11], the implied volatility σ := [uBS]−1(u) for a European call option in general LSV model is given by

σ = σ0+ ∞

∑

σn, (2.10)

where the first three orders are given by σ0= σ0,0, σ1= σ1,0+ σ0,1, σ2= σ2,0+ σ1,1+ σ0,2, σ3= σ3,0+ σ2,1+ σ1,2+ σ0,3, (2.11) and σ0,0=p2a0,0, σ1,0= a1,0 2σ0 (k − x), σ0,1= ta0,1(c0,0+ 2 f0,0) 4σ0 +a0,1c0,0 2σ₀3 (k − x), (2.12) where ai, j, ci, j, fi, j(i, j ∈ [0, 1]) are obtained in definition 4 by (2.9).

Remark2. Note that [uBS]−1 is an analytical function on its domain. The lengthy expression on σ2,0, σ1,1, σ0,2et al. are not provided in this project, one can find in [11] and the webpage

http://explicitsolutions.wordpress.com which is provided in [11]. In the proofs we use combined expressions for σ0, σ1etc.

Chapter 3

Analytical Results

3.1

Generalized Heston Model

The model of our interest is the generalized Heston model which is introduced in [7] see also [8], which is given by dSt = √ Z_tS_tdWt, S0= s > 0, dZt = κ(θ − Zt)dt + δ ZtξdBt, Z0= z > 0, dhW, Bit = ρdt. (3.1) The parameters that appear in the model are:

• S – Stock price • Z – Variance

• κ – Mean-reversion rate

• θ – Long-run average of the variance • δ – Volatility of volatility, i.e. vol-vol • ρ – Instantaneous correlation coefficient

• ξ – Constant elasticity for the variance process (abbreviated as constant elasticity) • t – Time to maturity

• W, B – Brownian motions

In [11], the third-order implied volatility expansion is obtained for the standard Heston model, i.e. model (3.1) with ξ = 0.5. Our task is to apply Proposition 1 to derive the third-order

implied volatility expansion for the generalized Heston model. Firstly, change of variables are performed in the system of SDE (3.1) and then we determine the infinitesimal generator to the resulting system.

Let (X ,Y ) := (log S, log Z), by applying Itô’s lemma1 [6], we obtain the following system of equations dXt= − 1 2e Yt_{dt + e}12Yt_dW t, X0= x := logs, dYt= (κθ e−Yt− 1 2δ 2_e2Yt(ξ −1)_{− κ)dt + δ e}Yt(ξ −1)_dB t, Y0= y := logz, dhW, Bit= ρdt. (3.2) From the system of equations (3.2) we can identify each drift- and diffusion rate, namely

µ1= − 1 2e y_, σ1= e 1 2y, µ2= κθ e−y− 1 2δ 2_e2y(ξ −1)_{− κ,} σ2= δ ey(ξ −1), (3.3) where µ denotes the drift rate and σ denotes the diffusion rate.

Proposition 2. The infinitesimal generator for the model (3.2) is A = 1 2e y_(∂2 x − ∂x) + (κθ e−y− 1 2δ 2 e2y(ξ −1))∂y+ 1 2δ 2 e2y(ξ −1)∂_y2+ ρδ ey(ξ − 1 2)_∂ x∂y. (3.4)

Indeed, the coefficients in(2.2) take the following form: a(x, y) =1 2e y_{, b(x, y) =} 1 2δ 2_e2y(ξ −1)_, c(x, y) = ρδ ey(ξ −12), _f(x, y) = κθ e−y−1 2δ 2_e2y(ξ −1)_{− κ.} (3.5) Proof. By applying the formula (2.3) in Definition 1

a(x, y) =σ 2 1 2 = 1 2(e 1 2y₎2 ₌1 2e y_, b(x, y) =σ 2 2 2 = 1 2(δ e y(ξ −1)₎2 ₌1 2δ 2_e2y(ξ −1)_,

The authors of [11] has considered the Taylor series expansion, about the initial values (x, y), of the generator (3.4) to present the formulas that give one type of implied volatility expansion. To derive the implied volatility expansions, we now apply formulas in Proposition 1 (and its completed form in [11] as described in Remark 2). In Proposition 3 below, we calculate the derivatives of the coefficients in the generator which are ingredients for the expansions (2.11) and (2.12).

Proposition 3. The derivatives of the coefficients in the generator (3.4) are a_0,0= 1 2e y_{, a} 0,1= 1 2e y_, a0,2= 1 4e y_{, a} 0,3= 1 12e y_, a_1,0= 0, a_1,1= 0, a_1,2= 0, a_2,0= 0, a2,1= 0, a3,0= 0, b_0,0= 1 2δ 2_ey(2ξ −2)_{, b} 0,1= (ξ − 1)δ2ey(2ξ −2), b0,2= (ξ − 1)2δ2ey(2ξ −2), b0,3= 1 3(ξ − 1) 3 δ2ey(2ξ −2), b_1,0= 0, b_1,1= 0, b1,2= 0, b2,0= 0, b2,1= 0, b3,0= 0, c0,0= ρδ ey(ξ − 1 2), _c0,1= ρδ (ξ −1 2)e y(ξ −1₂)_, c_0,2= 1 2ρ δ (ξ − 1 2) 2_ey(ξ −1₂)_{, c} 0,3= 1 12ρ δ (ξ − 1 2) 3_ey(ξ −1₂)_, c_1,0= 0, c_1,1= 0, c1,2= 0, c2,0= 0, c2,1= 0, c3,0= 0, f_0,0= κθ e−y−1 2δ 2_ey(2ξ −2)_{− κ, f} 0,1= −κθ e−y− (ξ − 1)δ2ey(2ξ −2), f_0,2= 1 2κ θ e −y_{− (ξ − 1)}2 δ2ey(2ξ −2), f0,3= − 1 12κ θ e −y₋1 3(ξ − 1) 3 δ2ey(2ξ −2), f1,0= 0, f1,1= 0, f1,2= 0, f2,0= 0, f_2,1= 0, f_3,0= 0.

Proof. By applying the η–function (2.9), a_0,0= ∂ 0 x∂y0a(x, y) 0!0! = a(x, y) = 1 2e y_, a_0,1= ∂ 0 x∂y1a(x, y) 0!1! = ∂ 1 ya(x, y) = 1 2e y_, a0,2= ∂_x0∂_y2a(x, y) 0!2! = 1 2∂ 2 ya(x, y) = 1 4e y_, a_0,3= ∂ 0 x∂y3a(x, y) 0!3! = 1 6∂ 3 ya(x, y) = 1 12e y_, a_1,0= ∂ 1 x∂y0a(x, y) 1!0! = ∂ 1 xa(x, y) =0, a_1,1= ∂ 1 x∂y1a(x, y) 1!1! = ∂ 1 x∂y1a(x, y) =0, a_1,2= ∂ 1 x∂y2a(x, y) 1!2! = 1 2∂ 1 x∂y2a(x, y) =0, a_2,0= ∂ 2 x∂y0a(x, y) 2!0! = 1 2∂ 2 xa(x, y) =0, a2,1= ∂_x2∂_y1a(x, y) 2!1! = 1 2∂ 2 x∂y1a(x, y) =0, a_3,0= ∂ 3 x∂y0a(x, y) 3!0! = 1 6∂ 3 xa(x, y) =0.

The same procedure is used for calculating the respective equations for b, c and f .

3.2

The Implied Volatility Expansions

In this section we present the analytical results for deriving the implied volatility expansion. We follow the results in [11] to evaluate a third-order expansion under the generalized Heston model. The closed-form solution for the expansions are presented. We truncate the expansion (2.10) by n terms and define the n-th order (approximating) expansion as

σ = σ0+ n

∑

Theorem 1 (Zeroth-order expansion). The zeroth-order implied volatility expansion under the generalized Heston model(3.1) is given by:

σ0= e

1 2y.

Proof. This can be proved by following Proposition 1 σ0= σ0,0=p2a0,0=

r 2 ·1

2e

y_{= e}12y_.

Theorem 2 (First-order expansion). The first-order implied volatility expansion under the generalized Heston model(3.1) is given by:

σ = σ0+ σ1, where σ1is σ1= 1 4(k − x)ρδ e y(ξ −1)₋1 8te 1 2y(δ2e2y(ξ −1)− ρδ ey(ξ − 1 2)+ 2κ(1 − θ e−y)).

Proof. This can be proved by following Proposition 1 σ1= σ1,0+ σ0,1 = (k − x)(2a0,0a1,0+ a0,1c0,0) 4√2a3/2_0,0 +ta0,1(c0,0+ 2 f0,0) 4√2√a_0,0 = (k − x)(2 · 1 2ey· 0 + 1 2eyρ δ ey(ξ − 1 2)) 4√2(1₂ey₎3/2 + t1₂ey(ρδ ey(ξ −12)+ 2κθ e−y− 21 2δ2e2y(ξ −1)− 2κ) 4√2(1₂ey₎1/2 = (k − x)ρδ e yξ +1₂y−3 2y 4 + t(ρδ ey+yξ −12y− 1 2y+ 2κθ ey−y− 1 2y− δ2ey+2yξ −2y− 1 2y− 2κey− 1 2) 8 = (k − x)ρδ e y(ξ −1) 4 + t(ρδ eyξ+ 2κθ e−12y− δ2ey(2ξ − 3 2)− 2κe 1 2y) 8 = 1 4(k − x)ρδ e y(ξ −1)₋1 8te 1 2y(δ2e2y(ξ −1)− ρδ ey(ξ − 1 2)+ 2κ(1 − θ e−y)).

Theorem 3 (Second-order expansion). The second-order implied volatility expansion under the generalized Heston model(3.1) is given by:

where σ2is σ2= 1 12(k − x) 2_{((ξ −}7 4)ρ 2₊1 2)δ 2_ey(2ξ −5₂) + t 24e −5₂y₍₍5 4− 2ξ )ρ 2 δ2ey(2ξ +1)+ 2δ2ey(2ξ +1)+ (3 4(κθ e −y₋1 2δ 2_e2y(ξ −1)_{− κ)}

− κθ e−y− (ξ − 1)δ2ey(2ξ −2)+ ξ (κθ e−y−1 2δ 2_e2y(ξ −1)_{− κ))tρδ e}y(ξ +5₂)_{+ (}5 4(κθ e −y −1 2δ

2_e2y(ξ −1)_{− κ) + 2(−κθ e}−y_{− (ξ − 1)δ}2_ey(2ξ −2)_))te3y_{(κθ e}−y₋1

2δ 2_e2y(ξ −1)_{− κ)} + (1 4ρ 2 ξ +1 4ρ 2₋1 4)tδ 2 e2y(ξ +1)) + 1 24t(k − x)e −5 2y_{((2ξ −}5 4)ρ 2 δ2ey(2ξ +1) + ((2ξ −3 2)(κθ e −y₋1 2δ

2_ey(2ξ −2)_{− κ) + 2(−κθ e}−y_{− (ξ − 1)δ}2_ey(2ξ −2)_{))ρδ e}y(ξ +3₂)_).

Proof. We start with the formula provided for σ2

σ2=σ2,0+ σ1,1+ σ0,2

= (k − x)

2

96√2a7/2_0,0

(32a3_0,0a_2,0− 9a2_0,1c2_0,0+ 4a_0,0(2a2_0,1b_0,0+ 2a_0,2c2_0,0+ a_0,1c_0,0(−5a_1,0+ c_0,1)) + 4a2_0,0(−3a2_1,0+ 4a1,1c0,0+ 2a0,1c1,0))

+ t 96√2a5/2_0,0

(a3_0,0(−2ta2_1,0+ 32a2,0) + 9a20,1c20,0− 2a0,0(8a0,2c20,0− 2a0,1c0,0(a1,0− 2c0,1)

+ a2_0,1(8b0,0+ 3t f0,0(c0,0+ f0,0))) + 2a20,0(−6a1,02 − ta0,1a1,0c0,0+ 2(−2ta20,1b0,0+ 4a1,1c0,0

+ 2a0,2(12b0,0+ t(c0,0+ 2 f0,0)2) + a0,1(−4c1,0+ tc0,0(c0,1+ 2 f0,1) + 2t f0,0(c0,1+ 2 f0,1)))))

+ t(k − x) 96√2a5/2_0,0

(−9a2_0,1c0,0(c0,0+ 2 f0,0) + a0,0(16a0,2c0,0(c0,0+ 2 f0,0) + a0,1(−10a1,0(c0,0

+ 2 f0,0) + 8(c0,1f0,0− c0,0(c0,1+ f0,1)))) + 8a0,02 (2a1,1(c0,0+ 2 f0,0) + a0,1(c1,0+ 2 f1,0))).

Then we divided the whole expression into three parts: σ₂(i)= (k − x)

2

96√2a7/2_0,0

(32a3_0,0a_2,0− 9a2_0,1c2_0,0+ 4a0,0(2a20,1b0,0+ 2a0,2c20,0+ a0,1c0,0(−5a1,0+ c0,1))

σ₂(iii)= t(k − x) 96√2a5/2_0,0

(−9a2_0,1c_0,0(c_0,0+ 2 f_0,0) + a_0,0(16a_0,2c_0,0(c_0,0+ 2 f_0,0) + a_0,1(−10a_1,0(c_0,0 + 2 f0,0) + 8(c0,1f0,0− c0,0(c0,1+ f0,1)))) + 8a0,02 (2a1,1(c0,0+ 2 f0,0) + a0,1(c1,0+ 2 f1,0))).

Then we substitute and simplify them one by one σ₂(i)= (k − x)

2

96√2a7/2_0,0

(32a3_0,0a_2,0− 9a2_0,1c2_0,0+ 4a0,0(2a20,1b0,0+ 2a0,2c20,0+ a0,1c0,0(−5a1,0+ c0,1))

+ 4a2_0,0(−3a2_1,0+ 4a1,1c0,0+ 2a0,1c1,0)) = (k − x) 2 96√2(1₂ey₎7/2(0 − 9( 1 2e y₎2 ρ2δ2e2y(ξ − 1 2)+ 4 ·1 2e y₍₂₍1 2e y₎21 2δ 2_e2y(ξ −1) + 2 ·1 4e y ρ2δ2e2y(ξ − 1 2)+1 2e y ρ δ ey(ξ − 1 2)(0 + ρδ (ξ −1 2)e y(ξ −1₂)_{)) + 4(}1 2e 2₎2_{(0 + 0 + 0))} =(k − x) 2 12e72y (−9 4ρ 2 δ2ey(2ξ +1)+ 2ey(1 4δ 2_e2yξ₊1 2ρ 2 δ2e2yξ+1 2(ξ − 1 2)ρ 2 δ2e2yξ) + 0) =(k − x) 2 12e72y (−9 4ρ 2 δ2ey(2ξ +1)+ (1 2+ ρ 2_{+ (ξ −}1 2)ρ 2_)δ2_ey(2ξ +1)₎ =(k − x) 2 12e72y ((ξ −7 4)ρ 2₊1 2)δ 2_ey(2ξ +1) = 1 12(k − x) 2_{((ξ −}7 4)ρ 2₊1 2)δ 2_ey(2ξ −5₂)_.

The same procedure for proving σ₂(ii), since it is not sensible to provide the complete proof here, thus it is provided in Appendix A.

σ₂(ii)= t 96√2a5/2_0,0

(a3_0,0(−2ta2_1,0+ 32a_2,0) + 9a2_0,1c2_0,0− 2a_0,0(8a_0,2c2_0,0− 2a_0,1c_0,0(a_1,0− 2c_0,1) + a2_0,1(8b0,0+ 3t f0,0(c0,0+ f0,0))) + 2a20,0(−6a21,0− ta0,1a1,0c0,0+ 2(−2ta20,1b0,0 + 4a1,1c0,0+ 2a0,2(12b0,0+ t(c0,0+ 2 f0,0)2) + a0,1(−4c1,0+ tc0,0(c0,1+ 2 f0,1) + 2t f0,0(c0,1+ 2 f0,1))))) = t 24e −5 2y((5 4− 2ξ )ρ 2 δ2ey(2ξ +1)+ 2δ2ey(2ξ +1)+ (3 4(κθ e −y₋1 2δ 2_e2y(ξ −1)_{− κ) − κθ e}−y − (ξ − 1)δ2ey(2ξ −2)+ ξ (κθ e−y−1 2δ 2_e2y(ξ −1)_{− κ))tρδ e}y(ξ +52)_{+ (}5 4(κθ e −y₋1 2δ 2_e2y(ξ −1)

− κ) + 2(−κθ e−y− (ξ − 1)δ2ey(2ξ −2)))te3y(κθ e−y−1 2δ 2_e2y(ξ −1)_{− κ) + (}1 2ρ 2 ξ +1 4ρ 2 −1 4)tδ 2_e2y(ξ +1)_),

σ₂(iii)= + t(k − x) 96√2a5/2_0,0

(−9a2_0,1c_0,0(c0,0+ 2 f0,0) + a0,0(16a0,2c0,0(c0,0+ 2 f0,0) + a0,1(−10a1,0(c0,0

+ 2 f0,0) + 8(c0,1f0,0− c0,0(c0,1+ f0,1)))) + 8a0,02 (2a1,1(c0,0+ 2 f0,0) + a0,1(c1,0+ 2 f1,0)))

= t(k − x) 96√2a5/2_0,0

(−9a2_0,1c2_0,0− 18a2_0,1c_0,0f_0,0+ a0,0(16a0,2c20,0+ 32a0,2c0,0f0,0

+ a0,1(8c0,1f0,0+ 8c0,0c0,1+ 8c0,0f0,1)))

= t(k − x) 96√2a5/2_0,0

(−9a2_0,1c2_0,0− 18a2_0,1c_0,0f_0,0+ 16a0,0a0,2c20,0+ 32a0,0a0,2c0,0f0,0

+ 8a0,0a0,1c0,1f0,0+ 8a0,0a0,1c0,0c0,1+ 8a0,0a0,1c0,0f0,1) = 1 24t(k − x)e −5 2y(−9(1 2e y₎2 ρ2δ2e2y(ξ − 1 2)− 18(1 2e y₎2 ρ δ ey(ξ − 1 2)_f 0,0 + 16 ·1 2e y1 4e y ρ2δ2e2y(ξ − 1 2)_{+ 32 ·}1 2e y1 4e y ρ δ ey(ξ − 1 2)_f 0,0+ 8 · 1 2e y1 2e y ρ δ (ξ −1 2)e y(ξ −1₂)_f 0,0 + 8 ·1 2e y1 2e y ρ δ ey(ξ − 1 2)_{ρ δ (ξ −}1 2)e y(ξ −1₂)_{+ 8}1 2e y1 2e y ρ δ ey(ξ − 1 2)_f0,1) = 1 24t(k − x)e −5₂y₍₋9 4ρ 2 δ2ey(2ξ +1)−9 2ρ δ e y(ξ +3₂)_f 0,0+ 2ρ2δ2ey(2ξ +1)+ 4ρδ ey(ξ + 3 2)_f 0,0 + 2ρδ (ξ −1 2)e y(ξ +3₂)_f 0,0+ 2ρ2δ2(ξ − 1 2)e y(2ξ +1)_{+ 2ρδ e}y(ξ +3₂)_f 0,1) = 1 24t(k − x)e −5₂y_{((2ξ −}5 4)ρ 2 δ2ey(2ξ +1)+ ((2ξ −3 2) f0,0+ 2 f0,1)ρδ e y(ξ +3₂) ) = 1 24t(k − x)e −5 2y((2ξ −5 4)ρ 2 δ2ey(2ξ +1)+ ((2ξ −3 2)(κθ e −y₋1 2δ 2_ey(2ξ −2)_{− κ)}

+ 2(−κθ e−y− (ξ − 1)δ2ey(2ξ −2)))ρδ ey(ξ +32)).

Now we compose the partitions into one expression σ2=σ₂(i)+ σ₂(ii)+ σ₂(iii)

= 1 12(k − x) 2_{((ξ −}7 4)ρ 2₊1 2)δ 2_ey(2ξ −52) + t 24e −5 2y₍₍5 4− 2ξ )ρ 2 δ2ey(2ξ +1)+ 2δ2ey(2ξ +1)+ (3 4(κθ e −y₋1 2δ 2 e2y(ξ −1)− κ) − κθ e−y − (ξ − 1)δ2ey(2ξ −2)+ ξ (κθ e−y−1 2e2y(ξ −1)− κ))tρδ ey(ξ +52)+ (5(κθ e−y−1 2e2y(ξ −1)

Theorem 4 (Third-order expansion). The third-order implied volatility expansion under the generalized Heston model(3.1) is given by:

σ = σ0+ σ1+ σ2+ σ3, where σ3is σ3= 1 192(k − x) 3 δ3ρ ey(3ξ −4)(2ξ − 3)(4ρ2ξ − 10ρ2+ 5) − 1 48δ ρt(k − x)e −9 2y(δ2ey(3ξ + 3 2)(45 8 ρ 2 + 6ρ2ξ2− 12ρ2ξ − 2ξ2+3 2ξ + 3 4) + δ 2_tey(3ξ +5₂)₍1 8ρ 2₊5 8ξ − 3 2ρ 2 ξ2+ ρ2ξ − 7 16) − δ ρtey(2ξ +3)(2κθ e−y(1 − ξ ) + δ2e2y(ξ −1)(−4ξ2+ 7ξ − 3)) + 1

16te

y(ξ +7₂)_(e2y(ξ −1)

δ2 + 2κ − 2κθ e−y)2(−4ξ2+ 4ξ +1

2) − te

y(ξ +72)_(e2y(ξ −1)_{(ξ − 1)δ}2_{+ κθ e}−y₎2₊3

4δ ρte y(2ξ +3) (4ξ2− 4ξ +1 2)( 1 2e 2y(ξ −1) δ2+ κ − κθ e−y) + tey(ξ + 7 2)(δ2e2y(ξ −1)(−6ξ2+ 11ξ − 5) + κθ e−y(3 − 2ξ ))(1 2e 2y(ξ −1) δ2+ κ − κθ e−y)) + 1 48δ 2_{t(k − x)}2_e−9₂y_(δ2_e4ξ y_(14ρ2 ξ −11 2 ρ 2 ξ2 −17 2 ρ 2_{− 2ξ + 2) + δ ρe}y(3ξ +3₂) (3ρ2ξ2− 7ρ2ξ +29 8 ρ 2₊5 4ξ − 7 8) + κe 2y(ξ +1)_{(1 − ξ − 5ρ}2 − 3ρ2ξ2+ 8ρ2ξ ) + κ θ ey(2ξ +1)(−5 2+ 10ρ 2_{+ ξ + 3ρ}2 ξ2− 10ρ2ξ )) + 1 96t 2_e−7₂y_(δ3 ρ ey(3ξ + 3 2) (3 8ρ 2_{+ (2ξ − 1)}2₍1 2− 3 2ρ 2_{) +}11 2 ξ − 11 4 ) − 3 32te 4y_(e2y(ξ −1) δ2+ 2κ − 2κθ e−y)3 + δ2e2y(ξ +1)((δ2e2y(ξ −1)(ξ − 1) + κθ e−y)(4ρ2(ξ −1

2) − 5) + (3 + 3 4ρ 2_{+ 6ρ}2 ξ2− 6ρ2ξ − 6ξ )(1 2e 2y(ξ −1) δ2+ κ − κθ e−y) + (8 − 4ρ2)(1 2κ θ e −y − δ2e2y(ξ −1)(ξ − 1)2)) + δ2tey(2ξ +3)((3 4− ρ 2_{− ρ}2

ξ )(e2y(ξ −1)(ξ − 1)δ2+ κθ e−y) + ρ2(1 2κ θ e −y_{− δ}2_e2y(ξ −1)_(ξ − 1)2) + (1 2+ 1 2ξ − 1 2ρ 2₊1 2ρ 2 ξ − 3ρ2ξ2)(1 2e 2y(ξ −1) δ2+ κ − κθ e−y)) − te4y(δ2e2y(ξ −1) (2ξ − 1)(ξ − 1) +1 2κ θ e

−y_)(e2y(ξ −1)

δ2+ 2κ − 2κθ e−y)2+ δ3ρtey(3ξ +

5 2)₍₋ 1 16− 1 16ρ 2 +1 2ρ 2 ξ −5 8ξ + 1 2ρ 2

ξ2) − 2te4y(e2y(ξ −1)(ξ − 1)δ2+ κθ e−y)2(1 2e

2y(ξ −1)

δ2+ κ − κθ e−y) + δ ρtey(ξ +72)((e2y(ξ −1)(ξ − 1)δ2+ κθ e−y)2+ (1

4ξ 2₊1 4ξ + 3 32)(e 2y(ξ −1) δ2+ 2κ − 2κθ e−y)2 + (e2y(ξ −1)(ξ − 1)δ2(6ξ − 1) + (2ξ + 1)κθ e−y)(1 2e 2y(ξ −1) δ2+ κ − κθ e−y))).

The complete proof to Theorem 4 is not provided in this paper, instead the formula for proving σ3, which is given in [11], is included in Appendix A. The starting point is same as Theorem

3, to have all partial derivatives calculated according to (2.9). Then applying the formula for deriving the expression of σ3. The MATLAB code for the third-order implied volatility

expansion is provided in Appendix C.

Remark3. According to [11, 12], the expansion converges in a certain region in the plane of log-moneyness and the time-to-maturity. We shall not discuss analytically the convergence issues here. Instead we derive these n-th order expansions under our model and use them as n-th order analytical approximations for the implied volatility. In section 4.3, we conduct numerical experimental studies to check the accuracy of these approximations up to the third order.

Chapter 4

Numerical Results

4.1

Procedure

There is no exact explicit formula for the implied volatility under the generalized Heston model (3.1), hence the Monte Carlo methods is applied in this project to implement the numer-ical analysis. We apply the Monto Carlo methods under the generalized Heston model to cal-culate the European call option price, then we use the built-in MATLAB function blsimpv to calculate the implied volatility σ for each European call option price. The implied volat-ility σ is used as the benchmark in this project. We compare the accuracy of our third-order implied volatility approximation ˆσ with the benchmark σ .

The relative error between the third-order implied volatility approximation ˆσ and the bench-mark σ is defined as follows

relative error = | ˆσ − σ | σ .

In this simulation we have a variance reduction technique, namely antithetic variates which improves the accuracy of the simulation [10, 14]. The Cholesky decomposition1is also used as a common method to create the second process so that it is correlated to the first process [10, 14]. We also ensure that the variance takes non-negative values as in the Cox–Ingersoll– Ross(CIR) model [5]. When simulating the option price we simulate one million different paths and the time step is approximately one hour of a business day, we chose these values since we aim for an approximation as close to the true value as possible.

The procedure of applying the Monte Carlo methods to calculate the European call option price is summarized in Algorithm 1 below. The MATLAB code for simulation is provided in Appendix C.

Algorithm 1 The European call option price under the generalized Heston (GHS) model Set: parameters, κ, θ , δ , ξ , ρ,t, K, S, Z in the GHS model (3.1), M = number of paths and N= time steps.

1: set dt = t/N

2: for j = 1 to M do

3: set S1= initial stock price

4: set V1= initial volatility

5: set S11= initial stock price

6: set V 11= initial volatility

7: for i = 1 to N do 8: generate w1_ji 9: generate w2_ji 10: w2_ji= w1_jiρ + w2_jip1 − ρ2 .Cholesky decomposition 11: Si+1= Si+ √ V_iS_i√dtw1_ji 12: Vi+1= Vi+ κ(θ −Vi)dt + δ (Vi)ξ √ dtw2_ji

13: Vi+1= Vi+1(Vi+1> 0) .Ensure non-negative value

14: S1i+1= S1i−

√ V_iS1i

√

dtw1_ji .Generating antithetic paths 15: V1i+1= V 1i+ κ(θ −V 1i)dt − δ (V 1i)ξ

√ dtw2_ji

16: V1i+1= V 1i+1(V 1i+1> 0)

17: end for

18: set calli= max(S(end) − K, 0)

19: set callAi= max(S1(end) − K, 0)

20: end for

4.2

Third-order Implied Volatility Expansion vs Monte Carlo

Approximation

The number of the simulation paths is M = 1 000 000 and the time step N = d2530te, where d.e rounds up to the closest integer of N. The time step N is approximately one hour of a business day as mentioned above. The data which are used in the numerical analyse are: ρ = −0.45, κ = 0.33, θ = 0.30, δ = 0.44, x = 0, and y = logθ , which are similar in the paper of Lorig et al. (2014). In the tables below following notation are used:

• C – European call option price • ε – Relative error

Table 4.1 presents results with t = 0.50, and ξ = 0.73.

k -0.75 -0.50 -0.25 0 0.25 0.5 0.75 C 0.5320 0.4101 0.2758 0.1511 0.0618 0.0175 0.0032 σ 0.5977 0.5763 0.5558 0.5388 0.5232 0.5112 0.5030 ˆ σ 0.5921 0.5740 0.5559 0.5389 0.5238 0.5116 0.5031 ε 0.0093 0.0041 2.1387 · 10−4 1.7045 · 10−4 0.0012 6.5902 · 10−4 9.6983 · 10−5 Table 4.1: A comparison of the implied volatility between our third-order expansion and Monte Carlo approximation with respect to different strike prices.

The results in Table 4.1 shows that when the log-moneyness (k − x) ∈ (−0.75, 0.75) and ξ = 0.73, the relative errors between our third-order implied volatility approximations and the benchmark are less than 1%. The log-moneyness (k − x) has negative correlation with both European call option price and the implied volatility, i.e. the European call option prices and the implied volatilities decrease when the log-moneyness increase. On the other hand, the European call options price and the implied volatility are positively correlated, i.e. the European call option prices decrease when the implied volatilities decrease.

Table 4.2 presents results with k = 0.25 and ξ = 0.73.

t 1/12 3/12 6/12 9/12 1 1.5 2 C 0.0037 0.0272 0.0619 0.0916 0.1168 0.1602 0.1963 σ 0.5308 0.5270 0.5236 0.5208 0.5167 0.5120 0.5073 ˆ σ 0.5303 0.5276 0.5238 0.5204 0.5173 0.5121 0.5080 ε 7.6764 · 10−4 0.0012 4.0672 · 10−4 7.7858 · 10−4 0.0011 2.8763 · 10−4 0.0014 Table 4.2: A comparison of the implied volatility between our third-order expansion and Monte Carlo approximation with respect to different maturities.

The results in Table 4.2 shows that when maturities t ∈ (1/12, 2), k = 0.25 and ξ = 0.73, the relative errors between our third-order implied volatility approximations and the benchmark are less than 0.2%. The European call option price has positive correlation with the maturity while the implied volatility has negative correlation with the maturity.

Table 4.3 presents results with k = 0.25 and t = 0.5. ξ 0.5 0.625 0.75 0.875 1 C 0.0598 0.0610 0.0623 0.0628 0.0638 σ 0.5145 0.5202 0.5251 0.5273 0.5311 ˆ σ 0.5147 0.5202 0.5245 0.5280 0.5309 ε 0.0012 3.0477 · 10−5 0.0011 0.0012 2.9841 · 10−4

Table 4.3: A comparison of the implied volatility between our third-order expansion and Monte Carlo approximation with respect to different ξ .

The results in Table 4.3 shows that when ξ ∈ (0.5, 1), k = 0.25 and t = 0.73, the relative errors between our third-order implied volatility approximations and the benchmark are less than 0.2%. Both European call option price and the implied volatility have positive correlation with ξ .

4.3

Comparison Between Different Orders of The Implied

Volatility Expansions

Figure 4.1 shows the comparison among the first-, the second- and the third-order implied volatility approximation under the generalized Heston model and the benchmark implied volatility which obtained by blsimpv with respect to different maturities while log-moneyness (k − x) ∈ (−0.75, 0.75). The data which are used in the plots are: ρ = −0.45, κ = 0.33, θ = 0.30, δ = 0.44, x = 0, y = logθ , ξ = 0.73, and t = [0.3125 0.625 1.25 2.5].

The results in Figure 4.1 show that when t = 0.3125, both the second- and the third-order ex-pansions are almost identical to the benchmark between log-moneyness (k − x) ∈ (−0.2, 0.4). When the log-moneyness (k − x) ∈ (−0.75, −0.2) ∪ (0.4, 0.75), the third-order expansion is better than the second oder expansions. The same results also show when t = 0.625. When t= 1.25, the third-order expansion is better than the second-order expansion in essentially the whole log-moneyness range except (k − x) ∈ (0.6, 0.75), where the second-order expansion is better than the third-order expansion. While when t = 2.5, there is only the third-order ex-pansion is chose to the benchmark. The first-order exex-pansion is not close to the benchmark whenever t = [0.3125 0.625 1.25 2.5], and the relative error increases when the maturity increases.

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 0.5 0.52 0.54 0.56 0.58 0.6 0.62 t = 0.3125 BS implied volatility First order Second order Third order -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 0.5 0.52 0.54 0.56 0.58 0.6 0.62 t = 0.625 BS implied volatility First order Second order Third order -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 0.5 0.52 0.54 0.56 0.58 0.6 0.62 t = 1.25 BS implied volatility First order Second order Third order -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 0.5 0.52 0.54 0.56 0.58 0.6 0.62 t = 2.5 BS implied volatility First order Second order Third order

Figure 4.1: A comparison among the different orders of the expansion and the Monte Carlo approximation with different maturities.

4.4

Summary

To summarize the numerical results we can conclude that the accuracy of the third-order ex-pansion is good. We studied three different parameters to check the accuracy, we can also con-clude that in all cases the relative error is always less than 1%. From the tables we can see that the strike price and the constant elasticity have a positive correlation with both the European call option price and the implied volatility. On the other hand, as the time to maturity increases we observe that the implied volatility decreases while the call option price increases. When we compared the different orders of the expansion, for log-moneyness between (−0.75, 0.75), it is observed that when time to maturity is longer the third-order expansion is remarkably better than the others. For shorter time to maturity, up to t = 0.625, the second-order is still a relative good approximation for the implied volatility.

Chapter 5

Conclusion

5.1

Conclusion

To conclude this project, by doing extensive calculations we succeeded to derive an explicit formula for the third-order expansion for the implied volatility, which is under the general-ized Heston model. The accuracy of the closed-form approximation is very good compared to the reference method, Monte Carlo simulation. Such expansions can be potentially used for calibrating the generalized Heston model to market implied volatility data. The numerical analysis showed results that confirms our expectations of the accuracy. The relative error is always less than 1% when considering the third-order expansion. The analysis also confirms the result in [11]. In our study the result also showed that both the second- and the third-order approximation are good, especially for short time-to-maturity 0 < t ≤ 0.625. Depending on the level of accuracy that is demanded, the second-order approximation works well for short time-to-maturity. Since the second-order expansion is remarkably shorter than the third-order it is much more convenient to implement. For longer time-to-maturity, t ≥ 1.25, the third-order expansion is required to approach the true value.

5.2

Further Research

Future research can be use this expansion of implied volatility to calibrate the generalized Heston model to market implied volatility surface. One more suggestion for further research is to study the model by varying more parameters in the numerical analysis, to observe how the expansions behave.

Appendix A

Proofs

Below is the proof for σ₂(ii) σ₂(ii)= t

96√2a5/2_0,0

(a3_0,0(−2ta2_1,0+ 32a_2,0) + 9a2_0,1c2_0,0− 2a_0,0(8a_0,2c2_0,0− 2a_0,1c_0,0(a_1,0− 2c_0,1) + a2_0,1(8b0,0+ 3t f0,0(c0,0+ f0,0))) + 2a20,0(−6a1,02 − ta0,1a1,0c0,0+ 2(−2ta20,1b0,0+ 4a1,1c0,0 + 2a0,2(12b0,0+ t(c0,0+ 2 f0,0)2) + a0,1(−4c1,0+ tc0,0(c0,1+ 2 f0,1) + 2t f0,0(c0,1+ 2 f0,1))))) = t 96√2a5/2_0,0 (9a2_0,1c2_0,0− 2a_0,0(8a_0,2c2_0,0+ 4a_0,1c_0,0c_0,1+ a2_0,1(8b_0,0+ 3tc_0,0f_0,0+ 3t f_0,02 )) + 2a2_0,0(2(−2ta2_0,1b_0,0+ 2a0,2(12b0,0+ t(c20,0+ 4c0,0f0,0+ 4 f0,02 )) + a0,1(tc0,0c0,1+ 2tc0,0f0,1 + 2tc0,1f0,0+ 4t f0,0f0,1)))) = t 96√2a5/2_0,0

(9a2_0,1c2_0,0− 2a0,0(8a0,2c20,0+ 4a0,1c0,0c0,1+ 8a20,1b0,0+ 3ta20,1c0,0f0,0+ 3ta20,1f0,02 )

+ 2a2_0,0(2(−2ta2_0,1b_0,0+ 2a0,2(12b0,0+ tc0,02 + 4tc0,0f0,0+ 4t f0,02 ) + ta0,1c0,0c0,1

+ 2ta0,1c0,0f0,1+ 2ta0,1c0,1f0,0+ 4ta0,1f0,0f0,1)))

= t 96√2a5/2_0,0

(9a2_0,1c2_0,0− 16a0,0a0,2c20,0− 8a0,0a0,1c0,0c0,1− 16a0,0a20,1b0,0− 6ta0,0a20,1c0,0f0,0

− 6ta_0,0a2_0,1f_0,02 + 4a2_0,0(−2ta2_0,1b_0,0+ 24a_0,2b_0,0+ 2ta_0,2c2_0,0+ 8ta_0,2c_0,0f_0,0+ 8ta_0,2f_0,02 + ta0,1c0,0c0,1+ 2ta0,1c0,0f0,1+ 2ta0,1c0,1f0,0+ 4ta0,1f0,0f0,1))

= t 96√2a5/2_0,0

= t 24e −5 2y(9(1 2e y₎2 ρ2δ2e2y(ξ − 1 2)− 16 ·1 2e y1 4e y ρ2δ2e2y(ξ − 1 2) − 8 ·1 2e y1 2e y ρ δ ey(ξ − 1 2)_{ρ δ (ξ −}1 2)e y(ξ −1₂)_{− 16 ·}1 2e y₍1 2e y₎21 2δ 2_ey(2ξ −2) − 6t1 2e y₍1 2e y₎2 ρ δ ey(ξ − 1 2)_f0,0− 6t1 2e y₍1 2e y₎2_f2 0,0− 8t( 1 2e y₎2₍1 2e y₎21 2δ 2_ey(2ξ −2) + 96(1 2e y₎21 4e y1 2δ 2 ey(2ξ −2)+ 8t(1 2e y₎21 4e y ρ2δ2e2y(ξ − 1 2)_{+ 32t(}1 2e y₎21 4e y ρ δ ey(ξ − 1 2)_f 0,0 + 32t(1 2e y₎21 4e y_f2 0,0+ 4t( 1 2e y₎21 2e y ρ δ ey(ξ − 1 2)_{ρ δ (ξ −}1 2)e y(ξ −1₂) + 8t(1 2e y₎21 2e y ρ δ ey(ξ − 1 2)_f 0,1+ 8t( 1 2e y₎21 2e y ρ δ (ξ −1 2)e y(ξ −1₂)_f 0,0+ 16t( 1 2e y₎21 2e y_f 0,0f0,1) = t 24e −5₂y₍9 4ρ 2 δ2ey(2ξ +1)− 2ρ2δ2ey(2ξ +1)− 2ρ2δ2(ξ −1 2)e y(2ξ +1)_{− δ}2_ey(2ξ +1) −3 4tρδ e y(ξ +5₂)_f 0,0− 3 4te 3y_f2 0,0− 1 4tδ 2_e2y(ξ +1)_{+ 3δ}2_ey(2ξ +1)₊1 2tρ 2 δ2e2y(ξ +1) + 2tρδ ey(ξ +52)_f 0,0+ 2te3yf0,02 + 1 2tρ 2 δ2(ξ −1 2)e 2y(ξ +1)_{+ tρδ e}y(ξ +5₂)_f 0,1 + t(ξ −1 2)ρδ e y(ξ +5₂)_f 0,0+ 2te3yf0,0f0,1) = t 24e −5₂y₍₍5 4− 2ξ )ρ 2 δ2ey(2ξ +1)+ 2δ2ey(2ξ +1)+ (3 4f0,0+ f0,1+ ξ f0,0)tρδ e y(ξ +5₂) + (5 4f0,0+ 2 f0,1)te 3y_f 0,0+ ( 1 4ρ 2 ξ +1 4ρ 2₋1 4)tδ 2_e2y(ξ +1)₎ = t 24e −5 2y((5 4− 2ξ )ρ 2 δ2ey(2ξ +1)+ 2δ2ey(2ξ +1)+ (3 4(κθ e −y₋1 2δ 2_e2y(ξ −1)_{− κ) − κθ e}−y − (ξ − 1)δ2ey(2ξ −2)+ ξ (κθ e−y−1 2δ 2_e2y(ξ −1)_{− κ))tρδ e}y(ξ +5₂) + (5 4(κθ e −y₋1 2δ 2_e2y(ξ −1)_{− κ)}

+ 2(−κθ e−y− (ξ − 1)δ2ey(2ξ −2)))te3y(κθ e−y−1 2δ 2_e2y(ξ −1)_{− κ) + (}1 2ρ 2 ξ +1 4ρ 2₋1 4)tδ 2_e2y(ξ +1)₎

Below is the formula for calculating σ3

σ3=σ3,0+ σ2,1+ σ1,2+ σ0,3

= (k − x)

3

192√2a11/2_0,0

(48a5_0,0a3,0+ 12a30,1c30,0− 9a0,0a0,1c0,0(2a20,1b0,0+ 2a0,2c20,0

+ a_0,1c_0,0(−4a_1,0+ c_0,1)) + a2_0,0(2c2_0,0(3a_0,3c_0,0+ a_0,2(−11a_1,0+ 3c_0,1)) + a0,1c0,0(36a21,0+ 20a0,2b0,0− 25a1,1c0,0− 11a1,0c0,1+ c20,1+ 2c0,0c0,2)

+ a2_0,1(−24a1,0b0,0+ 2b0,1c0,0+ 6b0,0c0,1− 13c0,0c1,0)) + 2a30,0(6a31,0+ 2a20,1b1,0

− a1,0(16a1,1c0,0+ 9a0,1c1,0) + 3c0,0(2a1,2c0,0+ a1,1c0,1+ 2a0,2c1,0) + a0,1(10a1,1b0,0

− 14a2,0c0,0+ c0,1c1,0+ 2c0,0c1,1)) + 4a40,0(−10a1,0a2,0+ 6a2,1c0,0+ 3a1,1c1,0+ 2a0,1c2,0))

− t(k − x) 768√2a9/2_0,0

(32a5_0,0(ta1,0a2,0− 12a3,0) + 105a30,1c30,0− 6a0,0a0,1c0,0(40a0,2c20,0

− 5a_0,1c_0,0(5a_1,0− 4c_0,1) + a2_0,1(40b_0,0+ t(2c2_0,0+ 15c_0,0f_0,0+ 15 f_0,02 )))

+ 2a2_0,0(−20ta3_0,1b_0,0c_0,0+ 72c2_0,0(a0,3c0,0+ a0,2(−a1,0+ c0,1)) + 2a0,1c0,0(3a21,0

− 18a1,0c0,1+ 2(−9a1,1c0,0+ 3c20,1+ 6c0,0c0,2+ a0,2(96b0,0+ 7tc20,0+ 36tc0,0f0,0+ 36t f0,02 )))

+ a2_0,1(24b0,1c0,0+ 72b0,0c0,1+ 14tc20,0c0,1− 60c0,0c1,0+ 48tc0,0c0,1f0,0+ 24tc0,1f0,02

− a_1,0(96b_0,0+ 7t(c2_0,0+ 6c_0,0f_0,0+ 6 f_0,02 )) + 24tc2_0,0f_0,1+ 48tc_0,0f_0,0f_0,1))

− 2a3_0,0(12a3_1,0+ ta0,1a21,0c0,0+ 2a1,0(8ta20,1b0,0− 2(4a1,1c0,0+ a0,2(40b0,0+ 5tc20,0

+ 16tc0,0f0,0+ 16t f0,02 )) − a0,1(tc0,0(5c0,1+ 8 f0,1) + 4(−3c1,0+ 2t f0,0(c0,1+ 2 f0,1)))) − 2(−2(−2c0,0(2a1,2c0,0+ a1,1c0,1) + 3a0,3c0,0(32b0,0+ 3t(c0,0+ 2 f0,0)2) + a0,2(16b0,1c0,0 + 32b0,0c0,1+ 9tc20,0c0,1− 20c0,0c1,0+ 24tc0,0c0,1f0,0+ 12tc0,1f0,02 + 12tc20,0f0,1 + 24tc_0,0f_0,0f_0,1)) + a_0,1(20ta_0,2b_0,0c_0,0− 3tc_0,0c2_0,1− 32b_0,0c_0,2− 6tc2_0,0c_0,2+ 12c_0,1c_1,0 + 8c0,0c1,1− 4tc20,1f0,0− 16tc0,0c0,2f0,0− 8tc0,2f0,02 + a1,1(56b0,0+ 3t(c0,02 + 8c0,0f0,0+ 8 f0,02 )) − 8tc0,0c0,1f0,1− 8tc0,1f0,0f0,1− 4tc0,0f0,12 − 8tc20,0f0,2− 16tc0,0f0,0f0,2) + a20,1(24b1,0 + t(2b0,1c0,0+ 6b0,0c0,1+ 3c0,0c1,0+ 8c1,0f0,0+ 8c0,0f1,0+ 16 f0,0f1,0)))) + 440,0(ta31,0

+ 2a_1,0(16a_2,0+ 2ta_1,1c_0,0+ ta_0,1c_1,0) + 2(2ta2_0,1b_1,0− 32a_0,2b_1,0− 24a_2,1c_0,0− 3ta_1,1c_0,0c_0,1 + 4a1,1c1,0− 6ta0,2c0,0c1,0− 6ta1,1c0,1f0,0− 12ta0,2c1,0f0,0− 2a1,2(32b0,0+ 3t(c0,0+ 2 f0,0)2)

+ t(k − x)

2

768√2a9/2_0,0

(63a3_0,1c2_0,0(c_0,0+ 2 f_0,0) − 4a_0,0a_0,1(10a2_0,1b_0,0(c_0,0+ 2 f_0,0) + 34a_0,2c2_0,0(c_0,0 + 2 f_0,0) + a_0,1c_0,0(22c_0,1f_0,0− 31a_1,0(c_0,0+ 2 f_0,0) + c_0,0(17c_0,1+ 12 f_0,1))) + 4a2_0,0(6c0,0(3a0,3c0,0(c0,0+ 2 f0,0) + a0,2(4c0,1f0,0− 5a1,0(c0,0+ 2 f0,0) + c0,0(3c0,1+ 2 f0,1))) + a0,1(−27a1,1c20,0+ 3c0,0c20,1+ 6c20,0c0,2− 54a1,1c0,0f0,0+ 2c20,1f0,0+ 8c0,0c0,2f0,0 + 15a2_1,0(c0,0+ 2 f0,0) + 20a0,2b0,0(c0,0+ 2 f0,0) + 4c0,0c0,1f0,1− a1,0(16c0,1f0,0+ c0,0(15c0,1 + 14 f0,1)) + 4c20,0f0,2) + a20,1(−15c0,0c1,0− 16c1,0f0,0+ 2b0,1(c0,0+ 2 f0,0) + 6b0,0(c0,1+ 2 f0,1) − 14c0,0f1,0)) − 8a30,0(a1,0(9a1,1(c0,0+ 2 f0,0) + 5a0,1(c1,0+ 2 f1,0)) + 2(−3(2a1,2c0,0(c0,02 f0,0) + a1,1(c0,1f0,0+ c0,0(c0,1f0,1) + 2a0,2(c1,0f0,0+ c0,0(c1,0f1,0))) − a0,1(2c0,0c1,1+ 2c1,1f0,0− 4a2,0(c0,0+ 2 f0,0) + c1,0f0,1+ c0,1(c1,0+ f1,0+ 2c0,0f1,1))) + 16a4_0,0(6a2,1(c0,0+ 2 f0,0) + 3a1,1(c1,0+ 2 f1,0) + 2a0,1(c2,0+ 2 f2,0))) + t 2 768√2a7/2_0,0 (−45a3_0,1c2_0,0(c0,0+ 2 f0,0+ 3a0,0a0,1(48a0,2c20,0(c0,0+ 2 f0,0+ a20,1(c0,02 f0,0) (16b0,0+ t(−c20,0+ 2c0,0f0,0+ 2 f0,02 )) − 4a0,1c0,0(3a1,0(c0,0+ 2 f0,0) − 2(4c0,1f0,0+ c0,0(3c0,1+ 2 f0,1)))) + 2a20,0(4ta30,1b0,0(c0,0+ 2 f0,0) − 24c0,0(3a0,3c0,0(c0,0+ 2 f0,0 + a0,2(4c0,1f0,0− a1,0(c0,0+ 2 f0,0) + c0,0(3c0,1+ 2 f0,1)) + 2a0,1(3a21,0(c0,0+ 2 f0,0) + a1,0(8c0,1f0,0 + c_0,0(6c_0,1+ 4 f_0,1)) − 2(3c_0,0c2_0,1+ 6c2_0,0c_0,2+ 2c_0,12 f_0,0+ 8c_0,0c_0,2f_0,0− a_1,1c_0,0(c_0,0+ 2 f_0,0) + a0,2(c0,02 f0,0)(32b0,0+ t(−c20,0+ 4c0,0f0,0+ 4 f0,02 )) + 4c0,0c0,1f0,1+ 4c20,0f0,2)) − a2_0,1(8b0,1(c0,0+ 2 f0,0) + ta1,0c0,0(c0,0+ 2 f0,0) + 2(−tc20,0c0,1+ 12b0,0(c0,1+ 2 f0,1) + 4 f_0,0(−4c_1,0+ t f_0,0(c_0,1+ 2 f_0,1)) + 2c_0,0(−5c_1,0+ tc_0,1f_0,0+ 4t f_0,0f_0,1− 2 f_1,0))))

− 2a3_0,0(4(−16a_0,2b_0,1c_0,0− 8a_1,2c2_0,0− 32a_0,2b_0,0c_0,1− 4a_1,1c_0,0c_0,1− 3ta_0,2c2_0,0− 32a_0,2b_0,0c_0,1 − 32a0,2b0,1f0,0− 16a1,2c0,0f0,0− 4a1,1c0,1f0,0− 12ta0,2c0,0c0,1f0,0+ 24a0,2c1,0f0,0

− 12ta0,2c0,1f0,02 − a1,0(6a1,1+ ta0,2c0,0)(c0,0+ 2 f0,0) − 3a0,3(c0,0+ 2 f0,0)(32b0,0+ t(c0,0+ 2 f0,0)2)

− 64a0,2b0,0f0,1− 24ta0,2f0,02 f0,1− 8a0,2c0,0f1,0) + 2ta20,1(2b0,1(c0,0+ 2 f0,0) + 6b0,0(c0,1+ 2 f0,1)

+ c0,0(c1,0+ 2 f1,0)) + a0,1(ta21,0(c0,0+ 2 f0,0) + 2a1,0(tc0,0(c0,1+ 2 f0,1+ 4(c1,0+ 2 f1,0)) − 2(−ta1,1c20,0+ tc0,0c20,1+ 32b0,0c0,2+ 2tc20,0c0,2− 8c0,1c1,0− 2ta1,1c0,0f0,0+ 2tc20,1f0,0 + 8tc_0,0c_0,2f_0,0− 16c_1,1f_0,0+ 8tc_0,2f_0,02 + 8a_2,0(c_0,02 f0,0) − 20ta0,2b0,0(c0,0+ 2 f0,0) + 4tc0,0c0,1f0,1− 8c1,0f0,1+ 8tc0,1f0,0f0,1+ 4tc0,0f0,12 + 8t f0,0f0,12 + 64b0,0f0,2+ 4tc20,0f0,2 + 16c_0,0f_1,1))) − 8a4_0,0(ta_1,0(a_1,1(c_0,0+ 2 f_0,0) + a_0,1(c_1,0+ 2 f_1,0)) − 4(2a_2,1(c_0,0+ 2 f_0,0) + a1,1(c1,0+ 2 f1,0) + 2a0,1(c2,0+ 2 f2,0))))

Appendix B

Matematical Derivations

Itô’s lemma:

Here is Itô’s lemma provided in detail, following [6]. Suppose that x follows the Itô process given by

dx = a(x,t)dt + b(x,t)dW, (B.1) where, dW is a Brownian motion and a and b are functions of x and t.

Here, a is the drift rate of x and b2 is the variance rate. Itô’s lemma shows that a sufficiently smooth function G of x and t follows the process

dG =  ∂ G ∂ xa+ ∂ G ∂ t + 1 2 ∂2G ∂ x2b 2  dt +∂ G ∂ xbdW. (B.2) There exists an extended version of Itô’s lemma which are used in differential calculus where we have

dS = µSdt + σ SdW, (B.3) as a model of stock price movements. From Itô’s lemma, it follows that the process followed by a function G of S and t is dG =  ∂ G ∂ Sµ S + ∂ G ∂ t + 1 2 ∂2G ∂ S2σ 2_S2  dt +∂ G ∂ Sσ SdW. (B.4) Cholesky decomposition

transpose of L. For a 2 × 2 matrix that is real, positive and definite we have a lower triangular matrix as following L=L1,1 0 L_2,1 L_2,2  , (B.6) then we determine A = LLT A=L1,1 0 L2,1 L2,2  L_1,1 L_2,1 0 L2,2  , (B.7) A=  L2_1,1 L1,1L2,1 L_2,1L_1,1 L2_2,1+ L2_2,2  , (B.8) hence, we have  1 ρ ρ 1  =  L2_1,1 L_1,1L_2,1 L_2,1L_1,1 L2_2,1+ L2_2,2  , (B.9)

from this we can see that one solution is

L_1,1= 1, L_2,1= ρ, L_2,2= q

1 − ρ2_.

This gives the following lower triangular matrix L= 1 0

ρ p1 − ρ2 

. (B.10)

Now, we generate two correlated Brownian motions B_t1, B2_t by dW1 t dW_t2  = LdB 1 t dB2_t  , (B.11) dW1 t dW_t2  = 1 0 ρ p1 − ρ2  dB1 t dB2_t  . (B.12) From (B.12) we can conclude the following

dW_t1= dB1_t, dW_t2= ρdB_t1+

q

1 − ρ2dB_t2, (B.13) where, dB1_t and dB2_t are independent Brownian motions. Hence, we have constructed two correlated Brownian motions according to a Cholesky decomposition.

Appendix C

MATLAB Codes

Appendix C contains MATLAB codes from procedures the Monte Carlo methods for evalu-ation the European call option price of the generalized Heston model, functions for third-order implied volatility expansion, the script for generating the model and the scrip for creating the subfigures.

Below is the code for calculating the European call option price under the generalized Heston model by applying the Monte Carlo methods.

function [C] = HEUCall(S0, k, T, V0, kappa, theta, sigma, rho,... N, M, xi)

% European call option pricing by the generalized heston model % S0 : the stock price

% k: the log strike % T: time to maturity % V0: the volatility dt=T/N; %Time increments for j=1:M %Number of paths

%Create zero vectors: S=zeros(N+1,1);

V=zeros(N+1,1); S(1)=S0;

for i=1:N %Number of steps w1=randn;

w2=randn;

w2=w1*rho+w2*sqrt(1-rho*rho); %Cholesky decomposition S(i+1)=S(i)+sqrt(V(i))*S(i)*w1*sqrt(dt);

V(i+1)=V(i)+kappa*(theta-V(i))*dt+sigma*(V(i))^(xi)... *w2*sqrt(dt);

V(i+1)=V(i+1)*(V(i+1)>0); %generating antithetic paths.

S1(i+1)=S1(i)-sqrt(V1(i))*S1(i)*w1*sqrt(dt); V1(i+1)=V1(i)+kappa*(theta-V1(i))*dt-sigma*(V1(i))^(xi)... *w2*sqrt(dt); V1(i+1)=V1(i+1)*(V1(i+1)>0); end call(j)=max(S(end)-exp(k),0); call1(j)=max(S1(end)-exp(k),0); end C=(mean(call)+mean(call1))/2; end

Below is the code for our third-order implied volatility expansion.

function sig = expansionformula(rho, delta, kappa, theta, xi, t,... k, S, Z) x = log(S); y = log(Z); sigma_0 = exp(y/2); sigma_1 = (1/4)*(k-x)*rho*delta*exp(y*(xi-1))-(1/8)*t*exp(y/2)*... (delta^2*exp(2*y*(xi-1))-rho*delta*exp(y*(xi-1/2))+... 2*kappa*(1-theta*exp(-y))); sigma_2 =(1/12)*(k-x)^2*(((xi-(7/4))*rho^2+(1/2))*delta^2*... exp(y*(2*xi-(5/2))))+(t/24)*(exp(-5*y/2))*(((5/4)-2*xi)... *rho^2*delta^2*exp(y*(2*xi+1))+2*delta^2*exp(y*(2*xi+1))... +((3/4)*(kappa*theta*exp(-y)-(1/2)*delta^2*exp(y*(2*xi-2))... -kappa)+(-kappa*theta*exp(-y)-(xi-1)*delta^2*exp(y*(2*xi-2)))... +xi*(kappa*theta*exp(-y)-(1/2)*delta^2*exp(y*... (2*xi-2))-kappa))*t*rho*delta*exp(y*(xi+(5/2)))+...

((5/4)*(kappa*theta*exp(-y)-(1/2)*delta^2*exp(y*(2*xi-2))... -kappa)+2*(-kappa*theta*exp(-y)-(xi-1)*delta^2*exp(y*(2*... xi-2))))*t*exp(3*y)*(kappa*theta*exp(-y)-(1/2)*delta^2*... exp(y*(2*xi-2))-kappa)+((-1/4)+(1/4)*rho^2+(1/2)*... rho^2*xi)*t*delta^2*exp(2*y*(xi+1)))+(1/24)*t*(k-x)*... exp((-5/2)*y)*((2*xi-(5/4))*rho^2*delta^2*exp(y*(2*xi+1))... +((2*xi-(3/2))*(kappa*theta*exp(-y)-(1/2)*delta^2*exp(y*(2*... xi-2))-kappa)+2*f01)*rho*delta*exp(y*(xi+(3/2)))); sigma_3 = (1/192)*(k - x)^3*delta^3*rho*exp(y*(3*xi -4))*... (2*xi - 3)*(4*rho^2*xi- 10*rho^2 + 5)-(t*(k - x))/(48*... exp(y)^(9/2))*((delta^3*rho*exp((3*y*(2*xi +1))/2))*(6*... rho^2*xi^2-12*rho^2*xi+45*rho^2/8-2*xi^2+3*xi/2+3/4)+... (delta^3*rho*t*exp((y*(6*xi + 5))/2))*(5*xi/8-7/16-3*rho^2*... xi^2/2+rho^2*xi+rho^2/8)-delta^2*rho^2*t*exp(y*... (2*xi + 3))*(delta^2*exp(y*(2*xi-2))*(-4*xi^2+7*xi-3)+... 2*kappa*theta*exp(-y)*(1-xi))+ (1/16)*delta*rho*t*exp((y*(2*xi +... 7))/2)*(exp(y*(2*xi - 2))*delta^2 + 2*kappa - 2*kappa*theta...

*exp(-y))^2*(-4*xi^2+4*xi+1/2)+(3/4)*delta^2*rho^2*t*...

exp(y*(2*xi + 3))*((exp(y*(2*xi - 2))*delta^2)/2 + kappa -...

kappa*theta*exp(-y))*(4*xi^2-4*xi+1/2)- delta*rho*t*exp((y*(2*xi... + 7))/2)*(exp(y*(2*xi - 2))*(xi - 1)*delta^2 +...

kappa*theta*exp(-y))^2+t*delta*rho*exp(y*(xi+7/2))...

*((exp(y*(2*xi - 2))*delta^2)/2 + kappa - kappa*theta*exp(-y))... *(delta^2*exp(2*y*(xi-1))*(-6*xi^2+11*xi-5)+kappa*theta*... exp(-y)*(3-2*xi)))+(t*(k - x)^2)/(48*exp(y)^(9/2))*(delta^4*... exp(4*xi*y)*(2-2*xi-(17/2)*rho^2+14*rho^2*xi-(11/2)*rho^2*xi^2)+... delta^3*rho*exp((3*y*(2*xi + 1))/2)*(-(7/8)+(29/8)*... rho^2+3*rho^2*xi^2+(5/4)*xi-7*rho^2*xi)+ delta^2*kappa*... exp(2*y*(xi + 1))*(1-xi-5*rho^2-3*rho^2*xi^2+8*rho^2*xi)+... delta^2*kappa*theta*exp(y*(2*xi +1))*(-5/2+10*rho^2+xi+3*... rho^2*xi^2-10*rho^2*xi))+(t^2/(96*exp(y)^(7/2)))*((delta^3*... rho*exp((3*y*(2*xi+ 1))/2))*((3/8)*rho^2+(1/2-(3/2)*rho^2)*... (2*xi-1)^2+(11/2)*xi-11/4)-(3*t*exp(4*y)*(exp(2*y*(xi - 1))... *delta^2 +2*kappa - 2*kappa*theta*exp(-y))^3)/32+ delta^2*... exp(2*y*(xi + 1))*(-5*(exp(y*(2*xi - 2))...

*delta^2)/2 +kappa - kappa*theta*exp(-y)))/2+ 4*rho^2*... (exp(y*(2*xi - 2))*(xi - 1)*delta^2 + kappa*theta*exp(-y))... *(xi - 1/2))+delta^2*t*exp(y*(2*xi + 3))*((rho^2*xi+rho^2-... (3/4))*f01+ rho^2*(f02)+(3*rho^2*xi^2-(1/2)*rho^2*xi+...

(1/2)*rho^2-(1/2)*xi-1/2)*f00 )+ t*exp(4*y)*((kappa*theta*... exp(-y))/2 -delta^2*exp(y*(2*xi - 2))*(xi - 1)^2- exp(y*... (2*xi - 2))*(xi - 1)*delta^2 - kappa*theta*exp(-y))*(exp(y*(... 2*xi - 2))*delta^2 + 2*kappa - 2*kappa*theta...

*exp(-y))^2+ delta^3*rho*t*exp((y*(6*xi + 5))/2)*((1/2)... *rho^2*xi^2+(1/2)*rho^2*xi-(1/16)*rho^2-(5/8)*xi-1/16) ... - 2*t*exp(4*y)*(exp(y*(2*xi -2))*(xi - 1)*delta^2 +... kappa*theta*exp(-y))^2*((exp(y*(2*xi - 2))*delta^2)... /2 + kappa - kappa*theta*exp(-y))+ delta*rho*t*exp((y*... (2*xi + 7))/2)*((3/32+(1/4)*xi+(1/4)*xi^2)*(exp(2*y*(xi... - 1))*delta^2 + 2*kappa - 2*kappa*theta*exp(-y))^2 +...

(exp(y*(2*xi - 2))*(xi - 1)*delta^2 + kappa*theta*exp(-y))^2 ... +( -kappa*theta*exp(-y)*(2*xi+1)-((6*xi-1)*(xi-1))*delta^2*... exp(2*y*(xi-1)))*(kappa*theta*exp(-y)-(1/2)*delta^2*...

exp(y*(2*xi-2))-kappa)));

sig = sigma_0 + sigma_1 + sigma_2 + sigma_3

Below is one example code for the script which was used to generate the model, calculate the benchmark implied volatility with the MATLAB function blsimpv and the relative errors. rho = -0.45; kappa = 0.33; theta = 0.3; delta = 0.44; xi = 0.5; t = 0.625/2; k = 0.6; S = 1; Z = 0.3; %Volatility N = ceil(t*2530); M = 1000000; r = 0;

sig = expansionformula(rho, delta, kappa, theta, xi, t, k, S, Z) CallPrice = HEUCall(S, k, t, Z, kappa, theta, delta, rho, N, M, xi) value = CallPrice;

blacksig = blsimpv(S, exp(k), r, t, value) err = abs(sig-blacksig)/blacksig

Below is the code for one subfigure in 4.1 when t = 0.3125. rho = -0.45; kappa = 0.33; theta = 0.3; delta = 0.44; xi = 0.73; t = 0.3125; k0 = -0.75; S = 1; Z = 0.3; %Volatility r = 0; n = 11; N = ceil(t*2530); M = 1000000; for i=1:n k(1)=k0;

sig(i) = expansionformula(rho, delta, kappa, theta, xi,... t, k(i), S, Z);

sigfirst(i) = sigmafirst(rho, delta, kappa, theta, xi,... t, k(i), S, Z);

sigsecond(i) = sigmasecond(rho, delta, kappa, theta, xi,... t, k(i), S, Z);

CallPrice(i) = HEUCall(S, k(i), t, Z, kappa, theta, delta,... rho, N, M, xi);

value(i) = CallPrice(i);

blacksig(i) = blsimpv(S, exp(k(i)), r, t, value(i)); k(i+1)=k(i)+abs(2*k0)/n; end k(:,n+1) = []; figure plot(k-log(S),blacksig,’k’) ylim([0.48 0.64]) yticks([0.48:0.02:0.64])

’third-order’},’Location’,’northeast’) % Move axes to origin.

ax = gca; % Get handles to axis ax.YMinorTick = ’on’;

ax.XMinorTick = ’on’;

ax.YAxisLocation = ’origin’; ax.XAxisLocation = ’origin’;

Bibliography

[1] A. Janek, T. Kluge, R. Weron, and U. Wystup, “Fx smile in the heston model,” in Stat-istical tools for finance and insurance, pp. 133–162, Springer, 2011.

[2] A. M. Chisholm, Derivatives demystified: a step-by-step guide to forwards, futures, swaps and options. John Wiley & Sons, 2011.

[3] F. Black and M. Scholes, “The pricing of options and corporate liabilities,” Journal of political economy, vol. 81, no. 3, pp. 637–654, 1973.

[4] F. D. Rouah, The Heston Model and Its Extensions in Matlab and C. John Wiley & Sons, 2013.

[5] J. C. Cox, J. E. Ingersoll Jr, and S. A. Ross, “A theory of the term structure of interest rates,” in Theory of valuation, pp. 129–164, World Scientific, 2005.

[6] J. C. Hull et al., Options, Futures, and Other Derivatives. Pearson Education Limited, 2018.

[7] L. B. Andersen and V. V. Piterbarg, “Moment explosions in stochastic volatility models,” Finance and Stochastics, vol. 11, no. 1, pp. 29–50, 2007.

[8] M. Altmayer and A. Neuenkirch, “Multilevel Monte Carlo quadrature of discontinuous payoffs in the generalized Heston model using Malliavin integration by parts,” SIAM Journal on Financial Mathematics, vol. 6, no. 1, pp. 22–52, 2015.

[9] M. Choudhry, Bond and money markets: strategy, trading, analysis. Butterworth-Heinemann, 2003.

[10] M. Kijima, Stochastic processes with applications to finance. CRC Press, 2013.

local-[14] P. Glasserman, Monte Carlo methods in financial engineering, vol. 53. Springer Science & Business Media, 2013.

[15] R. Merton, “The theory of rational option pricing,” Bell J Econ Manage Sci, vol. 4, pp. 141–183, 03 1973.

[16] S. L. Heston, “A closed-form solution for options with stochastic volatility with ap-plications to bond and currency options,” The review of financial studies, vol. 6, no. 2, pp. 327–343, 1993.