Normal Inverse Gaussian GARCH 模型與選擇權定價

國 立 交 通 大 學

應用數學系

碩 士 論 文

Normal Inverse Gaussian GARCH 模型

與選擇權定價

Normal Inverse Gaussian GARCH Model

And

Option Pricing

研 究 生：莊晉國

指導老師：許元春 教授

中 華 民 國 九 十 五 年 六 月

Normal Inverse Gaussian GARCH 模型

與

選擇權定價

Normal Inverse Gaussian GARCH Model

And

Option Pricing

研究生：莊晉國 Student： Jing-Guo Chuang

指導教授：許元春 Advisor：Dr. Yuan-Chung Sheu

國 立 交 通 大 學

應 用 數 學 系

碩 士 論 文

A thesis

Submitted to Department of Applied Mathematics

College of Science

National Chiao Tung University

In partial Fulfillment of the Requirements

for the Degree of

Master

In

Applied Mathematics

June 2006/6/29

Hsinchu, Taiwan, Republic of China

Normal Inverse Gaussian GARCH 模型

與

選擇權定價

學生：

莊晉國

指導教授

：許元春

國立交通大學

應用數學系﹙研究所﹚

摘要

這篇論文用NIG GARCH 的模型去描述財務市埸資產的log return。在這樣的模

型假設下，我們可以經由Esscher transform的方法來做資產的定價，而這種方法定

出來的價格可以用動態效用函數的架構來說明其合理性。

Normal Inverse Gaussian GARCH Model

And

Option Pricing

Student: Jing-Guo Chuang

Advisor: Yuan-Chung Sheu

Department of Applied Mathematics National Chiao Tung University

Hsinchu, Taiwan, R.O.C.

Abstract

This article uses the NIG GARCH model, the GARCH model with Normal inverse Gaussian innovation, to model the financial asset return. Under this model, we can pricing derivatives via Conditional Esscher transform. The pricing result can be justified by dynamic power utility framework.

誌

謝

首先感謝指導教授許元春老師，這幾年來悉心指導，讓我接觸了很多有趣的

課題。而老師也適時地給我一些建議和方向，讓我能夠抓住問題的重點和決解我

的問題，這讓我在學習研究的過程中受益良多。

其次，感謝研究所同學們，二年來一起讀書、討論問題和努力，有困難時能

夠互想幫助，大家也一起渡過很多快樂的時光。還有修課時的授課老師們，從課

堂上從他們身上學到很多的觀念和想法。

最後要感謝我所以有親人長輩和我的女友，總是掛念我的生活、學業、健康，

感謝他們的支持，陪伴我一起走過這段充實的研究生活。還有我的朋友也給了我

們多的鼓勵，讓我能夠順利的完學業。

Contents

1 Introduction 2 2 Generalized hyperbolic distribution 2

2.1 Generalized hyperbolic distribution . . . 2

2.2 Generalized inverse Gaussian distribution . . . 3

2.3 Alternative parameterization of NIG . . . 3

2.4 The GARCH NIG model . . . 5

3 Pricing Derivative Under NIG GRACH Model 5 3.1 Conditional Esscher Transform . . . 5

3.2 Change of measure for the NIG GARCH(1,1) model . . . 9

4 Estimation 11 5 Numerical Examples 11 6 Conclusion and Further Work 13 7 Appendix 13 7.1 Modified Bessel Functions . . . 13

7.2 Moment structure of Generalized Inverse Gaussian . . . 14

1

Introduction

GARCH is the widely used model to describe the time changed volatility in financial mar-ket. It successfully catch the volatility clustering, but the conditional normal distribution for the innovation still has poor performance to fit the high kurtosis, fat tailed and skewness of the real data. The normal inverse Gaussian , a special subclass of generalized hyper-bolic distribution, has been found the out-performance of fitting the financial return. [2] Barndorff-Nielsen(1997) propose the NIG stochastic volatility model, which is the ARCH type time series for the normal inverse Gaussian innovation. [1] Andersson(2001) extended the idea to GARCH type model with more flexible properties. [9] L. Forsberg using different type of parameterization to model the financial data.

For the purpose of pricing derivatives, we must find a equivalent martingale measure and then take the expectaion of payoff function under such probability measure. In the incom-plete market we have infinitely many equivalent martingale measure, so [11] Gerber (1994) propose a attractive approach to find a reasonable equivalent martingale measure. Using this method, we can find a martingale measure via Esscher transform under the independent increment model. under Garch or other conditional modeling, we can follow the [16] Tak Kuen Siu’s approach via the conditional Esscher transform.

2

Generalized hyperbolic distribution

2.1 Generalized hyperbolic distribution

Distributions that have tails heavier than the normal distribution are ubiquitous in finance. For purposes such as risk management and derivative pricing it is important to use relatively simple models that can capture the heavy tails and other relevant features of financial data. A class of distributions that is very often able to fit the distributions of financial data is the class of generalized hyperbolic distributions. This has been established in numerous investigations. In this paper, the generalized hyperbolic distribution has two different kind of parameter representation. Each representation has good properties in some aspect, and there is a one to one correspondence between two parameterization.

We say X ∼ GH(λ, α, β, δ, μ), if it has the density

fX(x) = ( γ δ)λ √ 2πKλ(δγ) K_λ−1 2(α  δ2+ (x − μ)2) (_δ2_{+ (x − μ)}2_/α)12−λe β(x−μ) _(2.1)

where the parameters satisfies

δ ≥ 0, α ≥ 0, α2> β2 if λ > 0.

δ > 0, α > 0, α2> β2 if λ = 0.

δ > 0, α ≥ 0, α2≥ β2 if λ < 0.

Let γ =α2− β2, γz =α2− (β − z)2

Then we can compute the moment generating function

M (z) = eμzγ

λ_K λ(δγz)

after calculating first and second moments the mean and variance are easily obtained EX = μ +δβKλ+1(δγ) γKλ(δγ) V arX = δKλ+1(δγ) γKλ(δγ) + β2δ2 γ2 ( Kλ+2(δγ) Kλ(δγ) − K_λ+12 (δγ) K_λ2(δγ) )

We focus on the special subclass of GH for λ = 1/2 called Normal Inverse Gaussian distri-bution. A random variable X ∼ N IG(α, β, δ, μ) if the pdf is represented as

fX(x) = α πe (δ√α2_−β2_−βμ)K1(δα  1 + (x−μ_δ )2)  1 + (x−μ_δ )2 eβx (2.2)

The corresponding moment generating function is:

M (z) = eμz+δ( √

α2_−β2₋√_α2_−(β+z)2₎

(2.3) Mean and variance:

EX = μ + βδ

γ (2.4)

V arX = δα2

γ3 (2.5)

IF Xi ∼ NIG(α, β, δi, μi) i = 1, 2, the we have X1+ X2∼ NIG(α, β, δ1+ δ2, μ1+ μ2).

2.2 Generalized inverse Gaussian distribution

Now, we consider the Inverse Gaussian distribution. X ∼ IG(δ, γ), where γ =α2− β2,

if it has the following density:

fX(x) = √ δ

2πx3e

−γ2(x− δγ )2_2x

which has mean and variance as:

EX = δ

γ V arX = δ γ3 2.3 Alternative parameterization of NIG

In order to have scale invariant properties for parameters except the scaling and location parameters, we use another parameterization for NIG distribution suggested by [9] Lars

Forsberg (2002) and [10] Lars Forsberg (2002) with α = _αδ β = _β/α μ = _μ δ = _δ/α ⇔ α = α/δ β = α/δβ μ = μ δ = √αδ

The density of the result parameterization,which we shall denote X ∼ N IG(α, β, μ, δ)

can be written with density

fX(x) = √ α π√δ e α√_1−β2_+β√_α_/δ_(x−μ₎K1  α  1 +(x−μ)_α_δ2   1 +(x−μ)_α_δ2 (2.6)

where 0 ≤ α_{, μ ∈ R, |β}| < 1 and 0 ≤ δ. Under this parameterization the moment

generating function becomes

M (z) = exp  μz + α1− β2−  α2− (αβ+√_α_δ_z)2  (2.7) According to the moment generating function for the new parameterization, we can see that

Lemma 2.1. The scaling properties of the NIG(α_{, β}_{, μ}_{, δ}_{) parameterization are given by}

the following. Let Z1 ∼ NIG(α, β, μ, δ), then cZ1+ d ∼ N IG(α, β, cμ+ d, c2δ), i.e α

and β does not change under scaling and shifting.

The first four central moments can be obtained by the cumulant generating function ln(M (u)) κ1 = EX = μ+ √ αδ  1− β2β  κ2 = V arX = δ  (1− β2)3/2 κ3 = 3δ 3/2_β √ α(1− β2)5/2 κ4 = 3δ 2_(4β2_{+ 1)} α(1− β2)7/2 Here, we have the skewness and kurtosis are given by

skewness = κ3 (κ2)3/2 = 3β √ α(1− β2)1/4 kurtosis = κ4 κ2₂ = 3(4β2+ 1) α1− β2

2.4 The GARCH NIG model

Assume Xt is the conditional return on its variance Zt is normally distributed

Xt|Zt∼ N(μt+



α/δβZt, Zt) (2.8)

where μt = E(Xt|Zt) is the conditonal mean, and the variance Zt is inverse Gaussian

distributed givenF_t−1 Zt|Ft−1∼ IG(  αδt,  α δ_t(1− β2)) (2.9) withE(Z_t|F_t−1) = √δt 1−β2.

Now, the Xt conditionally on Ft−1 are normal inverse Gaussian distributed

Xt|Ft−1∼ NIG(α, β, μt, δt) (2.10)

then the conditional variance of Xtis given by ht= V ar(Xt|Ft−1) = δ

 t

(1−β2₎3/2 which follows

the GARCH type

ht= ρ0+ q  i=1 ρi(Xt−i− EXt−i)2+ p  j=1 πjht−j (2.11)

we have the GARCH-NIG(1,1)

ht= ρ0+ ρ1(Xt−1− EXt−1)2+ π1ht−1 (2.12)

3

Pricing Derivative Under NIG GRACH Model

After constructing the sotck price model. We may be interesting in pricing the derivative. [11] Gerber and Shiu (1994) and [12] Gerber and Shiu(1996) provide a method to pricing options or derivative securities by Esscher transform, which is an efficient tool for pricing many option and contingent claims if the logarithms of the price of the primary securities are certain stochastic process with stationary and independent increment, more generally, the condition of stationarity and independence can be drop, then we can apply this method to GARCH models.

3.1 Conditional Esscher Transform

Our next step is to pricing option under GARCH model, then we need to consider the process without stationarity and independence. Instead of Esscher transform, we need a further tool to obtain the proper equivalent martingale by the so called conditional Esscher transform.

We consider discrete time financial model consisting of one risk-free interest rate r and one risky stock S. For generality, we assume that the innovations process of the underlying stock S is infinitely divisible and that the moment generating function of its distribution exits.

Suppose the filtered probability space (Ω, F, Ft, P) and τ be the time index set 0, 1, 2, . . . , T .

The stock price process follows

where Xt|Ft−1is normal inverse Gaussian distributed ,we denote it by Xt|Ft−1∼ NIG(α, 0, μ, δ)

as our second parameterization of NIG we have{δ_t}_t∈τ a conditional variance process of the

underlying stock. We suppose that δ_t ∈ Ft−1∀t ∈ τ\{0}. Now, assume that {Xt}t∈τ follows

a GARCH(1,1) process. So underP

∀t ∈ τ\{0}

δt = ρ0+ ρ1(Xt−1− μ)2+ π1δt−1 . (3.1)

where ρ0> 0, αi ≥ 0, and π1 ≥ 0 ∀ i = 1, 2

For the covariance stationarity of the GARCH model, we further impose the condition that

ρ1+ π1< 1 (3.2)

[6] Duan (1995) introduced the LRNVR for pricing and assumed that the martingale

measure Q with the LRNVR satisfies some conditions. [16] Tak Kuen Siu (2004) extend

the condition for Normal innovation to infinitely divisible distribution. Normal inverse Gaussian is indeed infinitely divisible, and the distribution is invariant under conditional Esscher transform, we also relax the invariant variance condition, we get

1. Q ∼ P

2. ln St

St−1 is normal inverse Gaussian distributed.

3. E_Q[ St

St−1|Ft−1] = er a.s.

In the sequel, we construct conditional Esscher transforms for the GARCH process{X_t}_t∈τ

associated with a sequence of conditional Esscher parameter{θ_t}_t∈τ.

Suppose{θ_t}_{t∈τ\{0}is a stochastic process with θt}∈ F_{t−1, for all t ∈ τ \{0}. Let MXt|Ft−1(z)}

be the moment generating function of the conditional distribution Xt|Ft−1 underP, where

z ∈ R. That is

MXt|Ft−1(z) := EP[ezXt|Ft−1] (3.3)

For all t ∈ τ \{0}. Let MXt|Ft−1(θ) exists, we define a sequence {Λt}t∈τ with Λ0 = 1 and

Λ_t= t  k=1 eθkXk MXt|Ft−1(θk) , t ∈ τ \{0} (3.4) [16] Tak K. S. (2004) Lemma 3.1. {Λt}t∈τ is a martingale.

Proof. E[Λt|Ft−1] =E t k=1 M_Xk|Fk−1eθkXk(θk)  Ft−1  =t−1_{k=1 M} eθkXk Xk|Fk−1(θk)E eθtYt M_Yt|Ft−1(θt)  Ft−1  = Λ_t−1 (3.5)

Let P_t := P|F_t ∀ t ∈ τ\{0}, and P_T = P. We define a family of probability measure

{Pt,Λt}t∈τ\{0} by the following conditional Esscher transform:

Pt,Λt(A|Ft−1) =EPt  IA e θtXt EPt(eθtXt|Ft−1)  Ft−1  A ∈ Ft (3.6) [16] Tak (2004) Lemma 3.2. Pt,Λt =Pt+1,Λt+1|Ft ∀ t ∈ τ\{0}

Proof. By the martingale property of Λ_{tt∈τ, ifA ∈ Ft}

Pt+1,Λt+1[A|Ft] =EPt+1  IA _E eθt+1Xt+1 Pt+1(eθt+1Xt+1|Ft)  Ft  = IAEPt+1  eθt+1Xt+1 E_Pt+1(eθt+1Xt+1_|Ft)  Ft  = IA (3.7)

takeE_Pt both side, we getP_{t,Λt(A) = Pt+1,Λt+1(A)}

The associated parameter θt is called the conditional Esscher parameter given Ft−1.

Write F (x; θt|Ft−1) for Pt,Λt(Xt≤ x|Ft−1) we have

F (x; θt|Ft−1) =

_y

−∞eθtxdF (x)

MYt|Ft−1(θt)

(3.8)

where F (x) is the cdf of N IG(α, 0, μ, δ).

Let MXt|Ft−1(z; θt) denote the moment generating function of F (x, θt|Ft−1), that is

MXt|Ft−1(z; θt) =

MXt|Ft−1(z + θt)

MXt|Ft−1(θt)

(3.9) For pricing a derivative V , we construct a martingale measure Q equivalent to P by adopt-ing Esscher Transforms.

First, choose a sequence of conditional Esscher parameters {θq_t}_t∈τ\{0} by solving the

fol-lowing equation

then we can define a family of probability measure{P_t,Λq

t}t∈τ\{0}associated with{θ

q

t}t∈τ\{0}.

Again, according to above result,

P_t,Λq

t =Ps,Λqs|Ft s, t ∈ τ with t ≤ s (3.11)

Lemma 3.3. Let Q = P_T,Λq

T the disounted stock price process {e

−rt_S t}t∈τ is a martingale underQ Proof. EQ[e−rtSt|Ft−1] _{= e}−rtE_Q[St|F_t−1] = e−rtEPT St−1eXt _{MXT |FT−1}eθqT XT_(θq T)  Ft−1  = e−rtSt−1EPt eXt eθqtXt M_Xt|Ft−1(θq_t)  Ft−1  = e−rtSt−1MXt|Ft−1(1+θ q t) M_Xt|Ft−1(θtq) = e−rtSt−1er = e−r(t−1)St−1. (3.12)

Then by risk-neutral pricing formula, the price of the derivative V at time t ∈ τ is:

Vt=EQ



e−r(T −t)VT Ft



(3.13)

we callQ a conditional risk neutralized Esscher pricing measure.

We can justify the pricing result by solving a dynamic utility maximization problem.

First, consider the sequence of power utility functions{u_t}_t∈τ with parameter {γ_t}_t∈τ

ut= ⎧ ⎨ ⎩ x1−γt 1−γt if γt = 1 ln x ifγt= 1 (3.14)

We assume that _{Vt is the agent’s price of the derivative V at time t with VT = VT}, such

that it is optimal for the agent not to buy or sell any unit of derivative V at time t.

Proposition 3.1. For all t ∈ τ , Vt= Vt

Proof. consider one stock S and one risk-free interest rate r. The agent owns m shares of

the stock bases his decision on a risk-averse utility function ut(x). Condider a derivative

security V pays Vt at time t. We have

is maximal for η = 0. From φ(0) = 0 ⇒ φ_{(η) = E}_{(Vt+1− e}r_Vt)u t(mSt+1+ η[Vt+1− erVt]) Ft  _(3.16) we obtain  Vt= e−rE[Vt+1u  t(mSt+1)|Ft] E[u t(mSt+1)|Ft] (3.17) note that φ(η) = E[(Vt+1− erVt)2ut(mSt+1+ η[Vt+1− erVt])|Ft] < 0 (3.18) if u_t(x) < 0.

In the particular case of a power utility function with parameter γt > 0, we have

ut(x) = x−γt. Then  Vt= e−rE[Vt+1(mSt+1) −γt|F_t] E[(mSt+1)−γt|F_t] = e −rE[Vt+1(St+1)−γt|Ft] E[(St+1)−γt|F_t] (3.19) since Vt+1 = St+1, we get Vt= St, So St= Vt= e−rE[(St+1) 1−γt_|Ft] E[(St+1)−γt|F_t] = e −r_S tE[e Yt+1(1−γt)_|Ft] E[eYt+1(−γt)|F_t] (3.20) this implies et= MYt+1|Ft(1;−γt) (3.21)

we see that the value of γt is −θt+1q , and by iteration we have γs is−θqs+1 ∀s ≥ t. We have

 Vt= e−r(T −t)EQ  Vt e θq_TYT MYT|FT −1    Ft  (3.22)

3.2 Change of measure for the NIG GARCH(1,1) model

It has been see that from [16] Tak(2004) the conditional Esscher transform perform the same result as [6] Duan (1995).

Now, we apply the conditional Esscher transform to the NIG GARCH(1,1) model.

Sup-pose the log return ln(St/St−1) ={Xt}t∈τ follows a GARCH(1,1) process and Xt|Ft−1 ∼

N IG(α, β, μt, δt) underP.more explicitly,

Xt = r + λ  ht+√αβ(1− β2)1/4  ht+  htεt where _εt∼ NIG(α_{, β}_{, −}√_α_β(1− β2)1/4_{, (1 − β}2)3/2)

where htfollow a GARCH(1,1) process as defined in (2.12), λ is the risk premium and r is

variance E(εt) = −√_α_β(1− β2)1/4+ √ α  1− β2β _{(1− β}2₎3/2_{= 0} V ar(εt) = (1− β 2₎3/2 (1− β2)3/2 = 1

the above specification implies that Xt ∼ NIG(α, β, r + λ√ht, (1 − β2)3/2ht), the

condi-tional mean and variance of Xt are

EXt = r + (λ +

√

αβ(1− β2)1/4)_ht

V arXt = ht

we can see λ +√αβ(1− β2)1/4 as total risk premium. Then we want to find the Esscher

parameter θq_t for this model by solving the following equation

er = M (θt+ 1) M (θt) (3.23) ⇒ r = r + λht+  α2− [αβ+  α(1− β2)3/2_htθt]2 (3.24) −  α2− [αβ+  α(1− β2)3/2_ht(θt+ 1)]2 (3.25) (3.26)

This equation can be solve explicitly by a quadratic form, the solution must satisfies |β+



(1−β2₎3_/2

α htθqt| < 1

Corollary 3.1. Under our assumption if Xt|Ft−1 is N IG(α, β, r + λ√ht, (1 − β2)3/2ht)

distribution under physical measure P, then under the risk neutral measure Q with Ess-cher parameters θq_t, the distribution of Xt|Ft−1 becomes N IG(α, β+



(1−β2₎3/2

α htθq_t, r +

λ√ht, (1 − β2)3/2ht)

Proof. the moment generating function MXt|Ft−1(z; θtq) is given by

MXt|Ft−1(z; θtq) (3.27) = MXt|Ft−1(z + θ q t) MXt|Ft−1(θqt) (3.28) = _{exp(r + λ}_ht+  α2− [αβ+  α(1− β2)3/2_htθt]2) (3.29) × exp(−  α2− [αβ+  α(1− β2)3/2_{ht(θt+ z)]}2) (3.30) = _{exp(r + λ}_{ht+ α}  1− [β+  (1− β2)3/2_ht α θt]2) (3.31) × exp(−  α2− [α(β+  (1− β2)3/2_ht α θt) +  α(1− β2)3/2_htz]2) (3.32)

therefore underQ Yt|Ft−1∼ NIG(α, β+  (1− β)3/2 α htθ q t, r + λ  ht, (1 − β2)3/2ht) (3.33)

Corollary 3.2. We have under original probability measure P the NIG GARCH model is:

Xt = r + λht+√αβ(1− β2)1/4ht+ ηt ηt =  htεt∼ NIG(α, β, −√αβ(1− β2)1/4  ht, ht(1− β2)3/2) ht = ρ0+ ρ1ηt−12 + π1ht−1

After Esscher transform(with Esscher parameters{θ_tq}), under Q the model becomes: Xt = r + λ  ht+√αβ(1− β2)1/4  ht+ ¯ηt ¯ ηt ∼ NIG(α, β+  (1− β2)3/2 α htθ q t, − √ αβ(1− β2)1/4_{ht, ht}(1− β2)3/2) ht = ρ0+ ρ1η¯t−12 + π1ht−1

4

Estimation

The parameters ω = (α, β, ρ0, ρ1, ρ2, π1, λ) ∈ Θ, where α ≥ 0, |β| < 1, ρ0> 0, ρ1, π1 ≥ 0,

ρ1+ π1< 1 and λ ∈ R.

Then we can estimate the parameters by maximizing the log likelihood function

L(α, β, ρ0, ρ1, π1, λ) = T  t=1 Lt(α, β, ρ0, ρ1, π1, λ) = (4.1) n 2ln α _{− n ln π −}1 2 T  t=1 ln δt+ n(α  1− β2_{) + β} T  t=1  α/δ(xt− μt) (4.2) −1 2 T  t=1 ln  1 +(xt− μt) 2 δ_tα  + T  t=1 ln K1  α  1 +(xt− μt)2 δ_tα  (4.3) where μt= r + λ√ht, δt= ht(1− β2)3/2 ht= ρ0+ ρ1(xt−1− mt−1)2+ π1ht−1 (4.4)

where mt= r + (λ +√αβ(1− β2)1/4)√ht All parameters with the following constraint

α ≥ 0, |β| < 1, ρ0 > 0, ρ1 and π1 ≥ 0, ρ0+ π1 < 1, λ ∈ R (4.5)

The gradient is in the appendix.

5

Numerical Examples

We present numerical results of our NIG GARCH pricing model using the colse values of S&P 500 daily index series from Jan 3,2000 to Apr 6,2006, a total of 1578 observation of

log return.

We estimate the parameters of NIG GARCH model via maximum likelihood estimation, the following are the estimated parameters.

α β λ ρ0 ρ1 π1

7.4820 −0.1302 0.3832 2 × 10−6 0.0884 0.8947

After estimating the parameters, we applying the estimated parameters to simulate the op-tion price by Monte Carlo simulaop-tion. In order to find the risk neutral Esscher parameter, solving a nonlinear equation is necessary. So the simulation is time consuming, and we use the basic Monte Carlo simulation without any variate reduction technique. For each option price, we produce 100, 000 paths to calculate it by the soft ware MATLAB 7.1. The fol-lowing tables show the comparison for the option pricing of NIG GARCH model and Black Schole pricing formula.(Assuming the risk free interest rate r = 0 and the initial conditional

volatility h0 is equal to the variance of the whole sample)

Maturity _S/K B-S NIG GARCH

30 _{1300/1000 300.0004 300.0725} 30 _{1300/1100 200.1114 200.1037} 30 _{1300/1200 104.0798 104.7615} 30 _{1300/1300 33.2843 32.0355} 30 _{1300/1400 5.3272 4.9065} 30 _{1300/1500 0.4026 0.5229} 30 _{1300/1600 0.0151 0.0707}

Maturity _S/K B-S NIG GARCH

60 _1300/1000 300.0579 300.0500 60 _1300/1100 201.3983 202.1829 60 _1300/1200 111.7655 111.2655 60 _1300/1300 47.0631 44.3521 60 _1300/1400 14.2924 12.7812 60 _1300/1500 3.1118 3.0055 60 _1300/1600 0.4970 0.7438

Maturity _S/K B-S NIG GARCH

90 _1300/1000 300.0579 300.9501 90 _1300/1100 201.3983 204.6041 90 _1300/1200 111.7655 117.8821 90 _1300/1300 47.0631 53.5759 90 _1300/1400 14.2924 19.5506 90 _1300/1500 7.2563 6.2298 90 _1300/1600 1.9249 2.0444

Comparing with the NIG GARCH option prices, the Black Scholes model always underprices deep out of the money options and it can under price for over price an out of the money options depending on the level of the initial conditional volatility.

6

Conclusion and Further Work

Empirical evidence shows that NIG GARCH model is more flexible to fit the real financial return than the GARCH model with Normal innovation. We can modeling the skewness ,higher Kurtosis for the conditional return. The Esscher transform provide a simple method to find a proper equivalent martingale measure via finding a Esscher parameters θ.

We can find that the Monte Carlo simulation is time consuming, so finding a good method to control variate of the simulation is important task. Since every time step of simulation, we have to solving a nonlinear equation to find the Esscher parameters. Reduce the number of simulation will reduce the simulation time.

7

Appendix

7.1 Modified Bessel Functions

The modified Bessel function of the third kind with order λ, denoted by Kλ(·). Here are

some properties useful.

• Integral representation Kλ(x) = 1₂  _∞ 0 y λ−1_exp(−x 2(y + y −1_{))dy x > 0 (7.1)} • Basic properties Kλ(x) = K−λ(x) (7.2) Kλ+1(x) = 2λ x Kλ(x) + Kλ−1(x) (7.3)

• Relation of Kλ(x) and Iλ(x), asymptotic properties.

Let Iλ(x) be the modified Bessel function of the first kind.

Kλ(x) = π₂ 1

sin (πλ)(I−λ(x) − Iλ(x)) (7.4)

Kλ(x) ∼ Γ(λ)2λ−1x−λ as x ↓ 0 (7.5)

K0(x) ∼ − ln x as x ↓ 0 (7.6)

• Series representation for λ = n + 1

2, n ∈ N K_n+1 2(x) =  π 2x −1/2_e−x_{(1 +}n i=1 (n + i)! (n − i)!i!(2x) −i_{) (7.7)} K−1/2(x) = K1/2(x) =  π 2x −1/2_e−x _(7.8)

• Derivatives w.r.t x K₀(x) = −K1(x) (7.9) Kλ(x) = −1₂(Kλ+1(x) + Kλ−1(x)) = −λ xKλ(x) − Kλ−1(x) (7.10) (ln Kλ(x)) = λ x − Rλ(x) (7.11) (ln Kλ(x)) = Sλ(x) − Rλ(x) x − λ x2 (7.12) where Rλ(x) := Kλ+1(x) Kλ(x) x > 0 (7.13) Sλ(x) := Kλ+2(x)Kλ(x) − K 2 λ+1(x) K_λ2(x) x > 0 (7.14) • Properties of Rλ and §λ R−λ(x) = 1 Kλ−1(x) (7.15) Rλ(x) = 2λ x + R−λ(x) (7.16) R_λ(x) = Rλ(x) x − Sλ(x) (7.17) R−1/2(x) = 1, R1/2(x) = 1 + 1_x, R−3/2(x) = _{x + 1}x (7.18)

7.2 Moment structure of Generalized Inverse Gaussian

If we say Z ∼ GIG(λ, δ, γ). The probability density function of Generalized Inverse Gaussian distribution is denoted by

fZ(z) = _γ δ _{λ z}λ−1 2Kλ(δγ)exp  −1 2( δ2 z + γ 2_z) _(7.19)

The moments of Z are given by

E[Zs_{] =}γ

δ

_sK

λ+s(δγ)

Kλ(δγ) (7.20)

and this formula holds for negative values of s, i.e. for inverse moments, too.

E[ln Z] = ∂E[Zs] ∂s   s=0 (7.21) where ∂E[Zs] ∂s = _γ δ _s ln  δ γ  Kλ+s(δγ) Kλ(δγ) + _γ δ _{s 1} Kλ+s(δγ) ∂ ∂sKλ+s(δγ) (7.22) using ∂ ∂sKλ+s(δγ) = ∂ ∂(λ + s)Kλ+s(δγ) ∂ ∂s(λ + s) (7.23)

and setting s = 0, gives

E[ln Z] = ln  δ γ  + 1 Kλ(δγ) ∂ ∂λKλ(δγ) (7.24)

7.3 Gradients of the GARCH NIG models.

In order to find the maximum likelihood estimator, we must derive the gradients of our likelihood function for the numerical method.

Consider the GARCH-NIG models in section 4, the log likelihood for one observation is given by Lt(α, λ, β, ρ0, ρ1, π1) = (7.25) 1 2ln α _{− ln π −}1 2ln δ  t+ α  1− β2_{+ β}  α/δt(xt− μt) (7.26) − 1 2ln  1 +(xt− μt)2 δtα  + ln K1 ⎛ ⎝α  1 +(xt− μt)2 δtα ⎞ ⎠ (7.27) where μt= r + λ√ht, δt= ht(1− β2)3/2 ht= ρ0+ ρ1(xt−1− mt−1)2+ π1ht−1 (7.28) where mt= r + (λ +√αβ(1− β2)1/4)√ht. ∇Lt(α, β, λ, ρ0, ρ1, π1) = ⎡ ⎢ ⎢ ⎣ ∂Lt(α,β,λ,ρ0,ρ1,ρ2,π1) ∂α .. . ∂Lt(α,β,λ,ρ0,ρ1,ρ2,π1) ∂π 1 ⎤ ⎥ ⎥ ⎦ 6×1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂Lt(α,β,μt,δt) ∂α ∂Lt(α,β,μt,δt) ∂β 0 .. . 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 ∂δt(β_∂β_,ht) ∂μt(α,β,λ,ht) ∂λ 0 0 0 .. . ... 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 6×2 + ⎡ ⎢ ⎢ ⎣ ∂ht(α,β,λ,ρ0,ρ1,ρ2,π1) ∂α .. . ∂ht(α,β,λ,ρ0,ρ1,ρ2,π1) ∂π 1 ⎤ ⎥ ⎥ ⎦ 6×1  ∂μt(α,β,λ,ht) ∂ht ∂δt(β,ht) ∂ht  1×2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ × ⎡ ⎢ ⎣ ∂Lt(α,β,μt,δt) ∂μt ∂Lt(α,β,μt,δt) ∂δt ⎤ ⎥ ⎦ 2×1

where ⎡ ⎢ ⎢ ⎣ ∂ht(α,β,λ,ρ0,ρ1,ρ2,π1) ∂α .. . ∂ht(α,β,λ,ρ0,ρ1,ρ2,π1) ∂π 1 ⎤ ⎥ ⎥ ⎦ 6×1 = ⎡ ⎢ ⎢ ⎣ ∂ht(α,β,λ,ρ0,ρ1,ρ2,π1,ht−1) ∂α .. . ∂ht(α,β,λ,ρ0,ρ1,ρ2,π1,ht−1) ∂π 1 ⎤ ⎥ ⎥ ⎦ 6×1 + ∂ht(α _{, β}_{, λ, ρ0, ρ1, ρ2, π1, ht−1)} ∂ht−1 ⎡ ⎢ ⎢ ⎣ ∂ht−1(α,β,λ,ρ0,ρ1,ρ2,π1) ∂α .. . ∂ht−1(α,β,λ,ρ0,ρ1,ρ2,π1) ∂π 1 ⎤ ⎥ ⎥ ⎦ 6×1

Each components can be written down as

∂Lt(α, β, μt, δt) ∂α = 1 2α +  1− β2+β(x−μt) 2√_α_δt + Jt(  1 +(xt− μt)2 αδt − (xt− μt)2) 2αδt  1 +(xt_α−μ_δt)2 t ) + (xt− μt) 2 2α2δt(1 +(xt_α−μ_δtt)2) ∂Lt(α, β, μt, δt) ∂β = − αβ  1− β2 + (xt− μt)  α δt ∂Lt(α, β, μt, δt) ∂μt = −β   α δt − Jt xt− μt  1 +(xt_α−μ_δt)2 t δt + xt− μt αδt(1 +(xt_α−μ_δtt)2) ∂Lt(α, β, μt, δt) ∂δt = − 1 2δt − β√α(xt− μt) 2δ3/2_t − Jt (xt− μt)2 2  1 +(xt_α−μ_δt)2 t δt2 + (xt− μt) 2 2α_tδ2(1 + (xt_α−μ_δt)2 t ) where Jt=−K0(α _{1 +}(xt−μt)2 α_δt ) K1(α  1 +(xt_α−μ_δtt)2) − 1 α  1 +(xt_α−μ_δtt)2 ∂μt(α, β, λ, ht) ∂λ =  ht ∂μt(α, β, λ, ht) ∂ht = λ 2√_ht ∂δt(β, ht) ∂β = −3ht  1− β2_β ∂δt(β, ht) ∂ht = (1− β 2₎3/2

∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂α = −ρ1 (xt−1− mt−1)β_√(1− β2)1/4ht−1 α ∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂β = −2ρ1(xt−1− mt−1) √ α(1− 2β2)_ht−1 2(1− β2)3/4 ∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂λ = −2ρ1(xt−1− mt−1)  ht−1 ∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂ρ0 = 1 ∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂ρ1 = (xt−1− mt−1) 2 ∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂ρ2 = (xt−1− mt−1) 2_I {xt−1−mt−1<0} ∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂π1 = ht−1 ∂ht(α, β, λ, ρ0, ρ1, ρ2, π1, ht−1) ∂ht−1 = −ρ1 (xt−1− mt−1)(λ +_ √αβ(1 − β2)1/4) ht−1 + π1

References

[1] Andersson J. 2001. ”On The Normal inverse Gaussian stochastic voaltility model”. Journal of Business and Economic Statistics 19: 44-54.

[2] Barndorff-Nielsen, O. E. 1997. ”Normal Inverse Gaussian Distributions and Stochastic Volatility Modelling,” Scandinavian Journal of Statistics, 24, 1-13.

[3] B. M. Bibby and M. Sorensen. ”Hyperbolic processes in finance,” 2001.

[4] Brendan O. Bradley and Murad S. Taqqu. 2001 Financial Risk and Heavy Tails. To Appear in the in 2002 in the volume ”Heavy-tailed distributions in Finance”, Svetlozar T. Rachev, editor, North Holland.

[5] Christian Bauer. 2000. Journal of Economics and Business, 2000, vol. 52, issue 5, pages 455-467 ”Value at risk using hyperbolic distributions,” Journal of Economics and Business, vol. 52, issue 5, pages 455-467.

[6] Duan, Jin-Chuan. 1995. The GARCH Option Pricing Model,” Mathematical Finance 5: 13-32.

[7] Engle, R.F. and K. Sheppard 2001. Theoretical and Empirical Properties of Dynamic Conditional Correlation Multivariate GARCH, UCSD Economics Discussion Papers 2001-15.

[8] Engle, R.F. 2002. Dynamic conditional correlation - a simple class of multivariate GARCH models. Journal of Business and Economic Statistics 20, 339-350.

[9] L. Forsberg. 2002. ”On the Normal Inverse Gaussian Distribution in Modeling Volatility in the Financial Markets,” Studia Statistica Upsaliensia, No. 5 Uppsala University Press.

[10] L. Forsberg and T. Bollerslev. 2002. ”Bridging the gap between the distribution of realized (ecu) volatility and ARCH modeling (of the euro): The GARCH-NIG model,” Journal of Applied Econometrics, 17(5):535V548.

[11] Gerber, Hans U., and Elias S. W. Shiu. 1994. ”Option Pricing by Esscher Transforms (with discussions),” Transactions of the Society of Actuaries 46: 99-191.

[12] Gerber, Hans U., and Elias S. W. Shiu. 1996. ”Actuarial Bridges to Dynamic Hedging and Option Pricing,” Insurance: Mathematics and Economics 18: 183-218.

[13] Oigard T.A., Hanssen A., Hansen R.E., Godtliebsen F. 2005. EM-estimation and mod-eling of heavy-tailed processes with the multivariate normal inverse Gaussian distrib-ution Signal Processing, Volume 85, Issue 8, Pages 1655-1673

[14] Prause, K., 1999, The Generalized Hyperbolic Model: Estimation, Financial Deriva-tives, and Risk Measures, PhD Thesis, Albert-Ludwigs-UniversitAat Freiburg.

[15] Rafael Schmidt ,T. Hrycej, E. Stutzle,2005. Multivariate distribution models with gen-eralized hyperbolic margins. Computational Statistics and Data Analysis, in print. [16] Tak Kuen Siu, Howell Tong and Hailiang Yang 2004. ”On Pricing Derivatives under

GARCH models: A Dynamic Gerber-Shiu’s Approach”, North American Actuarial Journal, 8(3), pp. 17-31.