Fourier methods for pricing early-exercise options under levy dynamics

(1)

by

Tolulope Rhoda Fadina

Thesis presented in partial fulfilment of the requirements for the degree of Master of Science in Mathematics at Stellenbosch University

Department of Mathematics, University of Stellenbosch,

Private Bag X1, Matieland 7602, South Africa.

Supervisor: Dr. P.W Ouwehand

(2)

Declaration

By submitting this thesis electronically, I declare that the entirety of the work con-tained therein is my own, original work, that I am the owner of the copyright thereof (unless to the extent explicitly otherwise stated) and that I have not previously in its entirety or in part submitted it for obtaining any qualification.

Signature: . . . . T.R Fadina

2012/01/08

Date: . . . .

(3)

Abstract

The pricing of plain vanilla options, including early exercise options, such as Bermu-dan and American options, forms the basis for the calibration of financial models. As such, it is important to be able to price these options quickly and accurately. Empirical studies suggest that asset dynamics have jump components which can be modelled by exponential Lévy processes. As such models often have characteristic functions available in closed form, it is possible to use Fourier transform methods, and particularly, the Fast Fourier Transform, to price such options efficiently. In this dissertation we investigate and implement four such methods, dubbed the Carr-Madan method, the convolution method, the COS method and the Fourier space-time stepping method. We begin by pricing European options using these Fourier methods in the Black-Scholes, Variance Gamma and Normal Inverse Gaussian mod-els. Thereafter, we investigate the pricing of Bermudan and American options in the Black-Scholes and Variance Gamma models. Throughout, we compare the four Fourier pricing methods for accuracy and computational efficiency.

(4)

Opsomming

Die prysbepaling van gewone vanilla opsies, insluitende opsies wat vroeg uitgeoefen kan word, soos Bermuda-en Amerikaanse opsies, is grondliggend vir die kalibrering van finansiële modelle. Dit is daarom belangrik dat die pryse van sulke opsies vin-nig en akkuraat bepaal kan word. Empiriese studies toon aan dat batebewegings sprongkomponente besit, wat gemodelleer kan word met behulp van exponensiëele Lévyprosesse. Aangesien hierdie modelle dikwels karakteristieke funksies het wat beskikbaar is in geslote vorm, is dit moontlik om Fourier-transform metodes, en in besonders die vinnige Fourier-transform, te gebruik om opsiepryse doeltreffend te bepaal. In hierdie proefskrif ondersoek en implementeer ons vier sulke metodes, genaamd die Carr-Madan metode, die konvolusiemetode, die COS-metode en die Fourier ruimte-tydstap metode. Ons begin deur die pryse van Europese opsies in die Black-Scholes, Gammavariansie (Engels: Variance gamma) en Normaal In-vers Gauss (Engels: Normal InIn-verse Gaussian)-modelle te bepaal met behulp van die vier Fourier-metodes. Daarna ondersoek ons die prysbepaling van Bermuda-en Amerikaanse opsies in die Black-Scholes en Gammavariansiemodelle. Deurlopend vergelyk ons die vier Fourier-metodes vir akkuraatheid en berekeningsdoeltreffend-heid.

(5)

Acknowledgements

To God be the glory great things he has done. I thank God for the gift of life and for his wisdom throughout this research. It has been God all the way.

My profound gratitude goes to my supervisor, Dr Peter Ouwehand. Thank you for the foundation you gave me and for your advice and support. I sincerely appreciate your efforts.

I sincerely appreciate the African Institute for Mathematical Sciences. And to the entire staff of Department of Mathematics, Stellenbosch University, I say THANK YOU.

To my parents, Prince and Mrs Michael Fadina , I say a very big thank you for the best legacy you set for me. To my siblings and in-laws, thanks for always been there. I love you all. To the entire RCCG family in Cape-town, thanks for your prayers and spiritual supports, God bless you all.

This research work was jointly funded by the African Institute for Mathematical Sciences (AIMS) and the Department of Mathematical Sciences at the University of Stellenbosch.

(6)

Dedications

To my loving parents Prince and Mrs Michael Fadina

(7)

List of Figures

2.1 Sample path of a stock driven by geometric Brownian motion. Parameters used: T = 1, r = 0.02, and σ = 0.7. . . 17 2.2 Sample path of Normal inverse Gaussian process: Parameters used: T =

1, σ = 0, θ = 4 and α = 0.05. . . 18 2.3 Sample path of a Variance Gamma process on a fixed grid. Parameters

used: T = 1, σ = 0.3, θ = 4 and κ = 0.05. . . 20 3.1 European calls under BS model. Parameters used: K = [50, 200], S0 =

100, σ = 0.25, r = 0.1, N = 29_{, β = 0.25 and α = 2.00.} _{. . . .} ₄₁

3.2 The absolute differences between the references prices and Carr-Madan method for pricing European calls under the BS model with respect to alpha (α). Parameters used: S0 = 100, K = [80, 100, 120], σ = 0.25,

r = 0.1, N = 210 and β = 0.25. . . 43 3.3 Error convergence of Carr-Madan method for pricing European call under

the BS model with respect to N. Parameters used: S0 = 100, K = 80,

T = 1, σ = 0.25, r = 0.1, α = 2.25 and β = 0.25. . . 44 3.4 European calls under BS model. Parameters used: K = [50, 200], S0 =

100, σ = 0.25, r = 0.1, N = 26 _{and K = 10.} _{. . . .} ₄₉

3.5 The absolute differences between the references prices and CONV method for pricing European calls under the BS model with respect to K. Pa-rameters used: S0 = 100, K = [80, 100, 120], σ = 0.25, r = 0.1 and

N = 26. . . 50 3.6 Error convergence of the CONV method for pricing European calls under

the BS, NIG and VG models with respect to K: Parameters used: T = 1, N = 26, S0 = 100, K = 100, σBS = 0.25, σN IG = 0.12, σV G = 0.12,

θV G= −0.14, κ = 0.2, α = 28.42141, β = −15.086 and δ = 0.31694. . . 50

3.7 Recovered density functions of BS, NIG and VG models. Parameters used: L = 10, N = 27_{, σ}

BS = 0.25, σN IG = 0.12, σV G = 0.12, r = 0.1,

α = 28.42, β = −15.08, δ = 0.31, θV G= −0.14 and κ = 0.2. . . 55

3.8 Effect of T on the COS method for pricing short-dated (T = 0.01) Euro-pean calls under BS model. Parameters used: S0 = 100, K = [50, 200],

σ = 0.25and r = 0.1. . . 59 3.9 European calls under BS model. Parameters used: S0 = 100, K =

[50, 200], σ = 0.25 and r = 0.1. . . 59

(10)

3.10 Effect of L on the COS method when pricing European calls under BS model with the put-call parity and without the put-call parity. Parame-ters used: S0, K = [0, 400], σ = 0.25, r = 0.1 and N = 64. . . 60

3.11 Error convergence of COS method for pricing European calls under BS model with respect to N: Parameters used: S0 = 100, K = [80, 100, 120],

r = 0.1, T = 1, L = 10. . . 61 3.12 Effect of T on the FST method for pricing European calls under BS

model. Parameters used: S0 = 100, K = [50, 200], σ = 0.25, r = 0.1,

N = 64 and x = 1.5. . . 71 3.13 Effect of x on European call prices in the BS model using FST. Parameters

used: S0= 100, r = 0.1, N = 128. . . 72

3.14 Error convergence of the FST method for pricing European calls under the BS, NIG and VG models with respect to (x). Parameters used: T = 1, N = 128, S0 = 100, KBS = 100, KN IG = 100, KV G = 100, σBS = 0.25,

σN IG = 0.12, σV G = 0.12, θV G = −0.14, κ = 0.2, α = 28.42141, β =

−15.086 and δ = 0.31694. . . 73 3.15 Short-dated European calls under BS model with different strikes. . . . 77 3.16 One year European calls under BS model with different strikes. . . 78 3.17 CPU-time for pricing one-year European calls under BS model, when

K = 100. . . 79 3.18 CPU-time for pricing one-year European calls under NIG model, when

K = 100. . . 80 3.19 One year European calls under NIG model with different strikes. . . 81 3.20 European calls under VG model with different maturities. . . 83 3.21 CPU-time for pricing a one-year European call under VG model, when

K = 90. . . 84 4.1 Early-exercise points x∗_{’s for pricing a one-year Bermudan put under the}

BS model and the VG model with M = 10, T = 1, N = 64, r = 0.1, L = 12, σBS= 0.2, σV G= 0.12, κ = 0.2 and θ = −0.14. . . 102

4.2 Error convergence for pricing a one-year Bermudan put under the BS model and VG model with respect to L with M = 10, T = 1, K = 110, S0 = 100, σBS = 0.2, σV G = 0.12, r = 0.1, θV G = −0.14, κ = 0.2, for

different values of N. . . 107 4.3 Error convergence for pricing a one-year Bermudan put using the FST

method under BS and VG models, with respect to xgrid. Parameters used:

S0 = 100, K = 110, r = 0.1, σBS = 0.25, σV G = 0.12, θV G = −0.14,

κ = 0.2, for different values of N. . . 110 4.4 One-year Bermudan put under the BS model. Parameters used: T = 1,

α = 0, L = 10 and x = 1.5 . . . 113 4.5 One-year Bermudan put under the VG model. Parameters used: T = 1,

α = 0, L = 13 and x = 1.5 . . . 115 4.6 One-year American put under the VG model. Parameters used: T = 1,

α = 0, L = 13 and x = 1.5 . . . 117

(11)

A.1 Sample path of a compound Poisson process with 10 jumps. Parameters used: T = 1, λ = 10. Jump sizes are drawn from standard normal distribution.. . . 123

List of Tables

2.1 The characteristic functions of the log-asset price of the BS, NIG and VG models, where %N IG and %V G are the characteristic exponents. . . 21

2.2 The cumulants of the log-asset price of the BS, NIG and VG models, where w is the drift correlation term which satisfies e−wt_{= Φ(−i, t)}_. _. ₂₁

3.1 Error convergence (log10 of the absolute error) of Carr-Madan method

for pricing European calls under the BS model with respect to alpha (α). Parameters used: S0 = 100, K = [80, 100, 120], α = [0.25, 1.5], σ = 0.25,

r = 0.1, N = 210 and β = 0.25. . . 43 3.2 Error convergence of the COS method for pricing European calls under BS

model with respect to L. Parameters used: S0 = 100, K = [80, 100, 120],

σ = 0.25, r = 0.1, N = 64 and T = 1. . . 60 3.3 Error convergence of the FST method for pricing European calls under BS

model with respect to N. Parameters used: S0 = 100, K = [80, 100, 120],

σ = 0.25, r = 0.1 and x = 1.5 . . . 71 3.4 Error convergence of the FST method for pricing European calls under

the BS model with respect to x: Parameters used: T = 1, σ = 0.25, r = 0.1and N = 128. . . 73 3.5 Parameters used for the implementation. Otherwise, stated. . . 76 3.6 Error and CPU(time-milli-seconds) for pricing a one-year European call

under BS model, K=100. . . 79 3.7 Error and CPU(time-milli-seconds) for pricing a one-year European call

under NIG model, K=100. . . 82 3.8 Error and CPU-time (milliseconds) for pricing a one-year European call

under VG model, K=90. . . 84 3.9 CPU-time (milliseconds) for pricing a one-year European call that

con-verges to an error of order O(10−6₎ _{under the BS, NIG and VG models}

using the CONV, COS and FST methods. . . 85 4.1 Parameters used for the implementation. Otherwise, stated. . . 112 4.2 Error and CPU (time-seconds) for pricing a one-year Bermudian put

un-der the BS model, K=110. . . 114 4.3 Error and CPU (time-seconds) for pricing a one-year Bermudian put

(12)

4.4 CPU-time (milliseconds) for pricing a one-year Bermudian put that con-verges to an error of order O(10−6₎ _{under the BS model and VG model}

(13)

Chapter 1

Introduction

In financial markets, the act of determining fast and accurate prices and sensitivities of options is an active research field. An option contract is an agreement between two parties - the holder of the option and the writer of the option. The holder of the option is given the right to make certain decisions in order to receive a certain payoff or she/he loses the premium paid for the option. The payoff is the difference between the underlying asset price St at a prescribed date t in the future, and the

predetermined strike price K. The two basic types of options are call options and put options. A call option gives the holder the right, not obligation, to buy St in

the future t, for K, while a put option gives the holder the right, not obligation, to sell St in the future t, for K. A call option has a positive value when St> K, and

a put option has a positive value when K > St. The payoff function of a call option

is given by

(St− K)+=

(

St− K if St≥ K,

0 otherwise .

And the payoff function of a put option is given by (K − St)+=

(

K − St if St≤ K,

0 otherwise .

Put-call parity establishes the relationship between the value of a call and a put option with similar strike price and maturity time T . The two commonly traded forms of vanilla options are European and early-exercise options.

A European option is a standard type of option contract with no special features except the simple maturity date T and the strike price K. This option can only be exercised at the maturity date, that is, the holder of this option can only make a decision at time T . The Bermudan option can be exercised at predetermined dates prior to, or on the maturity date while an American option can be exercised at any time before, or on the maturity date. In financial markets, practitioners prefer to know the worth of their options at all times in order to manage the risk associated with uncertainty, since the parameters of financial products constantly change over time. The Bermudan and American options are termed early-exercise options.

(14)

Starting from the work ofBlack and Scholes (1973), and Merton (1973), the dy-namics of the underlying asset price is based on the assumption that underlying asset returns are normally distributed with Brownian motion noise but fixed volatility 1_.

The pricing function for the European option in the form of a partial differential equation can be derived by hedging the option payoff function by the risk-neutral delta-hedging strategy. Observed option prices from financial markets show that dif-ferent volatilities should be used for difdif-ferent strikes and maturities. The behaviour of observed prices distinguishes between maturities which contradict option prices from the model driven by Brownian motion with constant volatility, for observed prices move by jumps especially in short-dated options (options traded on in days and months). This reveals that the normality assumption in the Black-Scholes (BS) theory cannot capture heavy tails and discontinuities present in short-term trading in the practical log-returns. The real densities are usually too peaked in short-term trading compared to the normal density Cont and Tankov (2004); interpreting the BS model as unpredictable.

Diverse models have been established in literature to improve the fit of the option prices in the BS model. Hull et al. (1987), Stein et al. (1991) and Heston (1993) assume that the stochastic volatility of the underlying asset is a mean-reverting diffusion process, typically correlated with the underlying process itself. Dupire (1994) and Derman et al. (1994) pioneered a model that describes volatility as the deterministic function of the asset price, known as the local volatility model. The local volatility model retains the pure one-factor diffusion approach of the BS model, but with an extension, as it defines the underlying volatility function as a function of asset price and time. Other models that replace the source of continuity in the behaviour of the underlying asset with some discontinuity, are the Lévy models, the Variance Gamma model (VG) by Madan and Seneta(1990), the Normal Inverse Gaussian model (NIG) byBarndorff-Nielsen(1997), and the Carr-Geman-Madan-Yor model (CGMY) by Carr et al. (2002). All these models have their advantages and disadvantages. In this dissertation, the focus is on models of asset dynamics that do not suffer from the disadvantages of models driven by pure Brownian motion. Empirical studies reveal that underlying asset prices do jump and the phenomenon of implied volatility smile in financial markets shows that the risk-neutral returns are not normally distributedBarndorff-Nielsen(1997),Madan and Seneta(1990),Bates (1996), Cont and Tankov(2004), Tankov(2010). Thus, the focus here is on pricing options where the underlying assets are driven by the exponential Lévy model.

Exponential Lévy models can be divided into three; the continuous exponential Lévy models - processes with no jumps, the only example is the geometric Brownian motion, the BS model; the finite activity exponential Lévy model - a process with continuous sample paths and few jumps, examples are the Merton Jump Diffusion model Merton (1976) and the Kou model Kou (2002); and the infinite activity ex-ponential Lévy model - a process with pure jumps, examples are the VG model, the NIG model and the CGMY model. Under these models, the characteristic functions

1_{The volatility measures the uncertainty associated with underlying asset pricing which is the} same as annualized standard deviation of returns.

(15)

(CFs) of their distributions are available in closed-form. The Lévy-Khinchine repre-sentation is frequently used to obtain the CF’s of Lévy processes,Schoutens(2003), ?,Cont and Tankov(2004).

In the option pricing theory under the BS model, the risk-neutral valuation formula (see Equation 2.38) presents the pricing formula of options as an expectation of the product of the discounted payoff and the probability density function of the underlying process except for options with early exercise features. But in the case of Lévy processes, the risk-neutral densities often comprise special functions and infinite summations. Thus, it is difficult to find their integrals directly. The Feynman-Kač theorem relates the risk-neutral valuation formula to the solution of partial differential equations (PDEs) under the BS model, and partial integro-differential equation (PIDE) under the Lévy model. The PDE equation can be solved based on the terminal conditions. The PDE can be used to price several types of options once the payoff function is known. The PIDE comprises a diffusion term and an integral term. These terms are often treated separately.

Based on the risk-neutral valuation formula, many numerical methods have been employed to value options and to determine their sensitivities e.g. PIDE-based meth-ods, Monte Carlo methodsBoyle(1977), lattice methodsKéllezi and Webber(2004) and many more. Option pricing in Lévy models takes into account empirical facts which often result in a more complicated computation of option prices. In view of this, a more efficient method and fast algorithm have to be developed to cope with these models. For early-exercise options, the solution to the PIDE is not trivial. Diverse finite difference schemes have been proposed in literature Cont and Tankov (2004),Hirsa and Madan(2004),Andersen and Andreasen(2000),Forsyth and Vet-zal (2002), but these methods are vulnerable to errors in convergence due to some factors that will be stated in Chapter 3 below, Lord et al. (2008), Jackson. et al. (2008),Surkov (2009). Thus, their results are not satisfactory.

Instead of applying the direct discounted expectation approach to evaluate the integral of the discounted payoff and risk neutral density function of the underlying process, it is easy to compute the integral of their Fourier transform once the CF of the distribution of the underlying asset is known. CF is the Fourier transform of the probability density function. In mathematical finance, diverse numerical integration methods based on Fourier methods have been established in literature. Some of these methods are the transform-based methods; the Carr Madan method Carr and Madan (1999), the Raible method Raible (2000), the convolution method (CONV method) Lord et al. (2008), the Lewis method Lewis (2001) and the Fourier space time-stepping method (FST method) Jackson. et al. (2008), Surkov (2009), and the Fourier-cosine series expansion (COS method)Fang and Oosterlee(2008), Fang (2010). The idea behind using the transform methods is to take an integral of the discounted payoff over the probability distribution obtained by inverting the corresponding Fourier transform; this originates in the convolution theorem (see Theorem 2.5.2 ) Carr and Madan (1999), Raible (2000), Lee (2004), Lord et al. (2008). The COS method substitutes the probability density function in the

(16)

risk-neutral valuation formula with the Fourier-cosine series expansion. In Lévy processes, the CF is often easier to handle than the density function itself Carr and Madan (1999), Lord et al. (2008), Cherubini et al. (2010). Thus, the Fourier methods are said to be effective approaches for pricing options in Lévy models because the CFs are readily available or can be calculatedRaible(2000),Lord et al.(2008),Fang and Oosterlee (2008) andFang and Oosterlee(2009b).

Fourier methods are computationally very efficient due to the availability of the Fast Fourier Transform (FFT), that is implemented in most of these methods Carr and Madan (1999), Lord et al. (2008), Surkov (2009). These methods are unique in pricing option with series of strikes in a single computation Carr and Madan (1999),Raible (2000),Surkov(2009),Fang and Oosterlee(2008),Lord et al.(2008). They can also be used to determine the hedge parameters at almost no additional computational cost, but in this dissertation, we will only focus on pricing, and not hedging, even though we know that pricing is very closely related to the ability to hedge. The FFT is an efficient algorithm for computing the sum of a finite sequence of complex number. It approximates a continuous Fourier Transform (CFT) by its discrete counterpart, Discrete Fourier Transform (DFT), for a carefully chosen vector Walker(1996), Cherubini et al.(2010).

Carr and Madan(1999) initiated the idea of pricing options using the FFT algorithm. An alternative formulation of transform-based methods for pricing and hedging of options in exponential Lévy models with the risk-neutral valuation formula is the FST method. The FST method involves solving the diffusion and the integral terms of the PIDE separately but in a proportional manner. The Fourier transform of the PIDE featured out the logarithm of the CF. Thus, the Fourier transform of the PIDE results in a system of ordinary differential equations (ODEs) that can be easily solved by breaking the equations up into simpler portions, solving each portion independently, and adding up the solutions. To price exotic options 2_{, the}

FST method is shown to handle prices and sensitivities effectively, Jackson. et al. (2008),Surkov (2009).

Furthermore, in pricing early exercise options, the Carr-Madan method is extended into a method that employs one of the crucial properties of the Fourier methods, the convolution theorem (Theorem2.5.2). As explained earlier, the main idea is to recog-nise the option prices as the convolution of the risk-neutral density and the discounted payoffs. The COS method is also extended to price the early exercise option Fang and Oosterlee(2009b). The FST method allows the prices from one exercise time to be projected back to a second exercise time in one step of the algorithm. Thus, the algorithm for pricing European options remains unchanged for pricing early-exercise options. The price of an American option can be approximated directly by that of a Bermudian option with a larger exercise dates. To price American options effi-ciently using the CONV and COS methods, the Richardson extrapolation method can be applied to the price of the Bermudan options of its counterpart. Lord et al.

2

Exotic options are path dependent options whose payoffs at maturity do not depend only on the price of the underlying asset at maturity but at several times before maturity, e.g. barrier, Bermudan and American options.

(17)

(2008) state that the pricing of American options using the PIDE approach is more favourable than the CONV method and the finite difference methods.

1.1 Motivation for research

Financial models always depend on unknown parameters. In order to provide param-eters matching the real market price of liquid instruments (calibration), the pricing of series of a options is required. In financial markets, early-exercise options are often used as building blocks for more complicated products and for effective calibration. This is quite time-consuming. Thus, we aim to compute the values of these options as fast as possible. As stated earlier, diverse numerical methods have been imple-mented in literature to adequately price these options and derive their hedge. The efficiency of a numerical method for pricing options depends on the following:

• Accuracy - Precise calculation of option prices and their sensitivities are the bedrock for using numerical methods, although numerical methods give ap-proximated results. In a financial market where huge number of assets change hands, a difference of 0.01 in the computed price can lead to a substantial loss or gain in a portfolio. Thus, the importance of good approximation cannot be overestimated.

• Speed - Practitioners operate in a highly competitive market. The ability to re-spond quickly to change in market conditions can give the practitioner an edge over competitors, especially in electronic financial markets where algorithm-trading is used in decision making on when and how options should be exer-cised.

• Convergence - The performance of numerical methods wholly depends on the rate at which the numerical solution converges to the exact solution as the number of points (N) involved in the computation increases. The difference between the approximated value and the exact value is the error. Higher-order convergence is preferred in most numerical methods, for the error decreases faster as N increases. In the case of exponential convergence, the error de-creases exponentially as N inde-creases. The computational speed is related to linear computational complexity, meaning the computational time (CPU-time) grows linearly only with respect to an increase in N. The convergence of the majority of numerical methods for solving financial problems are of first and second orders,Lord et al. (2008),Surkov(2009).

In this dissertation, Fourier’s methods for pricing plain-vanilla and early-exercise options that give exponential convergence and second-order convergence under ex-ponential Lévy models will be presented. The computer used for all programs is TOSHIBA with Pentium (R) Dual-core CPU, RAM 2.00GB (1.87 usable) and the system type is a 64 bits operating system. All programs are written in Python 2.6 on Linux 10.10.

(18)

1.2 Organisation of thesis

The general framework of our research, the exponential Lévy model, is explained in detail in Chapter 2 with applications to the option-pricing theory. Then, the risk-neutral valuation formula is presented. Chapter 3 introduces and explains the Fourier methods for approximating the risk neutral valuation integral. We conclude this chapter by presenting the numerical results obtained by pricing European style options in the BS, NIG and VG models using the Fourier methods. In Chapter 4, we extend the Fourier methods discussed in Chapter 3 to price Bermudan and American-style options. Lastly, we will summarize our findings and give concluding remarks.

(19)

Chapter 2

2. Models of asset prices

2.1 Introduction

In an arbitrage-free market the price of a contingent claim is given by the risk-neutral valuation formula. This formula has three main components: the risk free rate, the future payoff function and the risk-neutral distribution of the underlying asset returns. In this chapter, we discuss the risk-neutral distribution of the underlying asset.

The history of stochastic processes in finance can be traced back to Bachelier (1900), when Brownian motion was originally introduced as a stock-price model. One of the drawbacks of this model was that it accommodated negative stock prices which is not feasible in a real world market. After Bachelier’s work, Samuelson (1965), corrected this flaw by modelling the dynamics of stock prices using geometric Brownian motion instead of arithmetic Brownian motion. This, Samuelson (1965) carried out prior to the derivation of the popularly used BS formula.

Brownian motion is a popularly used Lévy process and is often used to derive ana-lytical formulae for solving option-pricing problems. The independent and stationary increments of Lévy processes motivate their applications to solving these problems. As stated earlier, financial models completely driven by geometric Brownian motion make up the BS model. Empirical studies show that an underlying asset driven by geometric Brownian motion does not work well in practice due to path-continuity and other factors Cont and Tankov(2004); many practitioners maintain that one need not to give up the continuity path if one accepts that a underlying asset is driven by stochastic volatility. Thus, it is important to work in a model where discontinuities cannot be ignored.

In financial modelling, the stationary increments, independent increments, stochas-tic continuity and the infinite divisible properties of asset returns motivate the use of Lévy processes to model asset prices. Lévy processes were founded and pioneered by Paul Lévy in 1930s. The class of infinitely divisible distribution is a crucial class of statistical distribution of Lévy processes. This is as a result of some properties: An infinitely divisible random variable can be expressed as a sum of large number of

(20)

independent and identically distributed (i.i.d) random variables. For every infinitely divisible distribution, there is an associated Lévy process Cont and Tankov(2004). These properties provide a means of representing changes in asset price as a result of a high rate of randomness in the financial market . The price changes in Lévy processes are consistent with the no-arbitrage assumption, and Lévy processes pro-vide a good fit with observed prices from financial markets,Barndorff-Nielsen(1997), Carr et al. (2002),Cont and Tankov (2004). The availability of the CF of the Lévy processes is of great importance in option pricing.

There are two major types of Lévy processes: The jump-diffusion processes and infinite-activity processes. Jump-diffusion processes exhibit much small-scale finite variate of jumps due to diffusion, with a few larger jumps. Infinite-activity processes also exhibit small-scale variation due to jumps with infinitely many small jumps in any non-zero time interval. Nevertheless, in any bounded time interval, there are only a finite number of jumps of larger magnitude than a given ε > 0. A good example of a process with finite variation and infinite activity is the variance-gamma process.

Intuitively, we can say in jump-diffusion models, jumps arrive at distinct times and the interval between the jump times can be modelled as diffusion terms. Carr et al. (2002), speculate in their investigation that the presence of diffusion terms in the dynamics of the stock price is not necessary. In practice, jump-diffusion processes are rarely used because of their weak response to the addition of extra parameters to the activities of the diffusion processes. That is, they still, almost surely generate sample paths which are continuous with respect to time. In view of this, we focus on models driven by infinite activity processes, that is, NIG and VG models. Their empirical performance in fitting equity and asset prices have been remarkably proven in literature.

The majority of theorems and definitions given in this dissertation are fairly gen-eral. Only examples and proofs for theorems which are of interest to us are given, whereas references to proofs of more general results will be given and some proofs will be provided in the appendix. For detailed study, books in the literature of Lévy processes in finance are Schoutens (2003), Cont and Tankov (2004) and Kypriano (2010), and, some general books on Lévy processes areSato(1999) andApplebaum (2004).

Section 2.2 presents the definition and properties of Lévy processes. We also present Lévy-Itô Decomposition that provides basic information on simulation of a Lévy sample path, and the Lévy-Khintchin formula that links Lévy processes to distributions. In section 2.3 , we shall discuss financial models driven by exponen-tial Lévy models, and how to deal with complex discontinuities in the numerical computation of the CFs of their distributions. Section 2.4explains the risk-neutral valuation formula. Lastly, we introduce the fundamental theory of Fourier methods.

(21)

2.2 Definitions and properties of Lévy processes

For a better understanding of the formal definition of Lévy processes, we need to know what càdlàg functions are. These kinds of functions appear throughout in financial mathematics and they play a vital roles in the theory of Lévy processes. Càdlàg functions are natural models for paths of processes with jumps.

Definition 2.2.1 (Càdlàg Function (continue á droite avec des limites á gauche)). A stochastic process (Xt)t≥0 is said to be càdlàg if, for every t ∈ [0, T ], the paths

t → Xt are continuous from the right and limited from the left at every point. That

is,

Left limit Xt−= lim s→t,s<tXt

Right limit Xt+= lim s→t,s>tXt

exist and X_t= Xt+ Sato(1999). Also the jump at time t is written as

∆Xt= Xt+− Xt−.

It is important to note that all continuous functions are càdlàg functions but not all càdlàg functions are continuous.

Definition 2.2.2 (Lévy process). A Lévy process is a càdlàg stochastic process (Xt)t≥0 with X0 = 0 defined on a filtered probability space (Ω, F, F, P) taking

val-ues in R, where F = (Ft)t≥0, and satisfies the following conditions Sato(1999):

1. Independent increments: For every increasing sequence of times t0, . . . , tn,

given j = 1, . . . , n, the random variables X_t_j− X_t_j−1 are independent.

2. Stationary Increments: Xt+u− Xt does not depend on t but Xt+u− Xt∼ Xu.

So that, the distribution of X_t+u− X_t depends only on u, ∀u > 0. 3. Stochastic Continuity: For all ε > 0, limh→0P(|Xt+h− Xt| ≥ ε) = 0.

Condition 1 implies for a given information at time t, a change in the stochastic process Xt+u− Xt is independent of the past information. Processes with stationary

increments imply the increments in Xt+u− Xt, ∀ t, u > 0 have the same distribution

for every time t. i.e, a change in the distribution is independent on time. Condition 3 means that for a given time t, the probability of seeing a jump at t is zero. In other words, a jump appears at random times.

It is important to recall that Brownian motion is an example of a Lévy process, for it satisfies the above conditions. Sato(1999) gives a detailed study on Brownian motion as an example of Lévy processes. The class of Lévy processes include some other processes, for example; Poisson processes, Compound Poisson processes and many more, but the Poisson processes and the Brownian motion are the fundamental examples of Lévy processes. They can be thought of as the building blocks of Lévy processes because every Lévy process is a superposition of Brownian motion and, possibly an infinite number of independent Poisson processes, Cont and Tankov (2004). For more information on the examples of Lévy processes, see AppendixA.1.

(22)

2.2.1 Infinitely divisible distribution

Lévy processes are processes with some strong properties and one of the properties is infinite divisibility.

Definition 2.2.3 (Infinitely divisibility). A random variable Y (or distribution) is said to be infinitely divisible if it can be written as a sum of n independent and identically distributed random variables {Yj}. That is

Y ∼ Y1+ Y2+ · · · + Yn, for all n ≥ 2.

This implies the distribution of Yj depends on n not on j. The distribution function

P is infinitely divisible if and only if the CF Φ(·) is, for every n, the n-th power of some CF Φn(·).

Proposition 2.2.4. Let X_tbe a Lévy process that is sampled at a set of evenly spaced discrete time. For any n ∈ N,

Yj = Xtj n − Xtj−1 n and Xt= n X j=1 Yj. (2.1)

Equation2.1implies Yj is independently and identically distributed because Xt

sat-isfies the independent stationary increments property, then, Xt has an infinitely

divisible law. Conversely, given an infinitely divisible distribution P, there exists a Lévy process Xtwhere the distribution of increments is governed by P Sato(1999).

The normal, gamma, Poisson and α-stable distribution are examples of infinitely divisible distributions. This is known by studying their CFs, see Sato (1999) for details.

2.2.2 Characteristic function

One of the arguments that strengthens the use of CFs is that the CF of a random variable may be known in closed form, particularly for infinitely divisible random variables. As a result of the one-to-one relationship between the CF of a process and the probability density function, once the CF function is known, then the distribution function can be determined by the Fourier inversion (see Section 2.5for detail). The COS method has shown to be a good method for obtaining densities from CF’s. Definition 2.2.5 (Characteristic function). Let X be a random variable, then its characteristic function Φ : R → C is defined as

ΦX(ω) := E[eiω·X], ω ∈ R. (2.2)

If the random variable has a probability density, fX(x), then the CF is the Fourier

transform of fX(x). Thus, equation2.2 becomes

ΦX(ω) := E[eiω·X] =

Z

R

(23)

Proposition 2.2.6 (Characteristic exponent). Let (Xt)t≥0 be a Lévy process in R

and f_X(x) be the probability density of X. Then, there exists a continuous function % : R → C known as the characteristic exponent of X such that

E[eiω·Xt] := et%X(ω), ω ∈ R (2.4) and ΦX(ω) := E[eiω·Xt] = Z R eiω·Xt_f X(x)dx, ω ∈ R (2.5)

is the characteristic function of X_t. By virtue of equation2.5, %Xt(0) = 1.

The relationship between the moment-generating function M(ω) and the CF Φ(ω) is given by

M(ω) = Φ(−iω). (2.6)

Proposition 2.2.7 (Properties of Characteristic function). Let X be a Lévy process in R and ΦX(ω) denote the CF. The following are properties of CF Sato(1999);

• Let X1 and X2 be two Lévy processes and Y = X1+ X2. The CF of ΦY(ω) is

the product of the CF of the two Lévy processes. i.e. Φ_Y(ω) = ΦX1(ω)·ΦX2(ω).

• Let X1 and X2 be two Lévy processes. For any linear function X1 = p + X2q,

ΦX1(ω) = e

ipω_{· Φ}

X2(qω).

• E[eiωX_{] always exists since |e}iωX_{| is continuous and bounded by 1 ∀ ω.}

• ΦX(ω) = ΦX(−ω) .

• ΦX(0) = 1 , for any distribution.

Example 2.2.8 (Characteristic function of a normal distribution). The probability density function of a normal distribution is given by

fX(x) = 1 √ 2πσ2e −(x−µ)2 2σ2 ,

where X is a random variable. Using equation 2.5, the CF is given by ΦX(ω) = Z R eiω√ 1 2πσ2e −(x−µ)2 2σ2 dx = eiµω− 1 2σ 2_ω2 . (2.7)

And the characteristic exponent of a normal distribution is given as %(ω) = log(eiµω−σ2ω22 ) = iµω −σ

2_ω2

(24)

The moments

The derivation of the moments of a distribution from the CF is another crucial tool. The mean, variance, skewness and the kurtosis are the moments we care about in many risk-management applications. Let X be a real-valued random variable with probability distribution P.

• The n-th moment M(Xn₎ _{of the distribution P is}_Sato₍₁₉₉₉₎

M = (i−n)d

n_Φ X(ω)

dωn |ω=0. (2.8)

Recall the characteristic exponent of X is the logarithm of the CF ΦX.

• The n-th cumulant, c_n is defined as cn= (i−n)

dnln ΦX(ω)

dωn |ω=0. (2.9)

• Mean of X: This can be calculated by setting n = 1 in equation (2.9), that is the first cumulant. It is also the expectation of the possible values of X in the distribution Cont and Tankov(2004).

Mean = c1 = (i−1)

d ln ΦX(ω)

dω |ω=0= E[X]. (2.10)

• Variance of X: Is the second cumulant, that is, n = 2 in equation (2.9). The expectation of the square of the difference between the random variable X and the expectation of X, is also known as the variance of X Cont and Tankov (2004). That is,

Var = c2= (i−2)

d2ln ΦX(ω)

dω2 |ω=0= E[(X − E(X))

2_]. _(2.11)

• Skewness of X : This measures the level of asymmetry of probability distribu-tion P of a real-valued random variable X. If the tail on the left side of the probability density function is larger than the tail on the right side and more values lie to the right of the average, this is called negative skewness. If the tail on the right side of the probability density function is larger than the tail on the left side and more values lie to the left of the average, this is called positive skewness. If the size of the right tail is the same as with the left tail, that is, the values are evenly distributed on both sides of the average, this is called undefined skewness, Cont and Tankov (2004). Let n = 3 in equation (2.9),

c3 = (i−3)

d3ln ΦX(ω)

dω3 |ω=0. (2.12)

The skewness can be represented as: Skew = c3 (√c2)3 = E[(X − E(X)) 3_] (pE[(X − E(X))2])3 (2.13)

(25)

• Kurtosis of X : This measures the peakedness of X in P. The fourth cumulant that is, when n = 4 in equation (2.9) divided by the square of the second cumulant (variance) in equation (2.11) minus 3 is known as Kurtosis, Cont and Tankov(2004). Let

c4 = (i−4)

d4_{ln Φ} X(ω)

dω4 |ω=0, (2.14)

the kurtosis is given by Kurt = c4 (c2)2 − 3 = E[(X − E(X)) 4_] (E[(X − E(X))2_])2 − 3. (2.15) 2.2.3 Lévy-Itô Decomposition

As the name implies, Lévy-Itô Decomposition decomposes a Lévy process X as the sum of continuous terms and discontinuous terms. The decomposition is useful since it allows one to compute the CF of any Lévy process with ease. Every Lévy process can be approximated by arbitrary precision of a jump-diffusion process that is, a sum of a compound Poisson process and Brownian motion with drift, possibly more than one, Cont and Tankov(2004). Note that the distribution of every Lévy process is uniquely determined by its characteristic triplet (γ, σ, ν) called Lévy triplet of the process. γ ∈ R is the drift parameter, σ2 _{≥ 0} _{is the diffusion parameter and ν is}

the Lévy measure. The Lévy measure is the expected number of jumps in a given interval per unit time. The pure jump component is characterized by the jumps density, which is called the Lévy density Sato(1999).

Theorem 2.2.9. Every Lévy process (X_t)t≥0 on R can be written as

Xt= γt + σBt+ Xtl+ lim ε↓0 ˆ X_tε (2.16) where • (B_t)t≥0 is a Brownian motion.

• σ is the diffusion parameter.

• γ is the deterministic drift parameter.

• ν is a positive measure (Lévy measure) on R that satisfies: Z

R−{0}

(1 ∧ |x|2)ν(dx) < ∞. (2.17)

Equation 2.17 implies this measure has no mass at the origin, but infinitely many jumps can occur around the origin. Thus, the measure must be a squared integrable around the origin.

• From equation (2.16), γt + σBt is a continuous Gaussian function and every

(26)

• Xl

t+ limε↓0Xˆtε is the discontinuous part of equation (2.16).

•

∆Xl= |∆X| ≥ 1 and X_tl= X

0≤s≤t

∆X_sl,

where ∆X is the jump size, ∆Xl describes large jumps with an absolute size greater than 1, and X_tl is the sum of a finite number of jumps in the time interval of 0 ≤ s ≤ t. • ∆Xε= ε ≤ |∆X| < 1 and Xˆ_tε= X 0≤s≤t ∆X_sε. (2.18) ˆ

X_tε is the sum of a possibly infinite number of small jumps in the limit ε ↓ 0 in the time interval of 0 ≤ s ≤ t. In a situation where ε ↓ 0, the process can have infinitely many small jumps, therefore X_tε may not converge. Thus, the infinite arrival rate of small jumps at zero exists.

The case where the Lévy process Xt is of finite variation is a crucial condition

for simplifying a Lévy-Itô Decomposition. For the proof of Theorem2.2.9, seeSato (1999) and Cont and Tankov(2004).

2.2.4 Lévy-Khinchine Representation

The CF of a Lévy process is easily derived with the aid of the Lévy-Khinchine Representation, once the result of the Lévy-Itô Decomposition is available. This plays a crucial role in option pricing for it forms the basis for using Fourier methods. Recall the expression for Lévy-Khinchine representation:

ΦX(ω) = E[eiω·Xt] = et%X(ω), ω ∈ R (2.19) with %X(ω) = − 1 2ω 2_σ2_{+ iγω +} Z R (eiωx− 1 − iω · xI|x|≤1)ν(dx), (2.20)

where (γ, σ, ν) is the characteristic triplet of a Lévy process Xt with respect to the

truncating function, %X(ω) is the characteristic or Lévy exponent, σ2 ≥ 0, γ ∈ R

and ν is a positive measure on R such that R_R−{0}(1 ∧ |x|2_{)ν(dx) < ∞}_{. For detail,}

see Theorem 8.1 in Sato(1999).

2.3 Models driven by exponential Lévy processes

In this section, we present some models driven by the exponential Lévy model. i.e, BS model, NIG model and VG model. These models are selected due to their empirical behaviour. We pay less attention to the density function and concentrate more on the CF of the distribution because the Fourier pricing formula mostly depend on the CF of the distribution and not the density function. We will compute the moments of the distribution, i.e, mean, variance, skewness and the kurtosis, and also discuss the

(27)

numerical behaviour of CF of the distribution of these exponential Lévy models. We continue by presenting the dynamics of the underlying stock price under exponential Lévy models.

Assume a market consisting of one risk-less asset (the bond) with a price process given by Bt = ert, where r is the compound interest rate, and one risky asset (the

stock), where the stock-price process (St)_t∈R is assumed to have the form:

St:= S0eXt (2.21)

and (Xt)t≥0 is a stochastic process that satisfies some integrability condition (see

Section2.5).

2.3.1 Black-Scholes model as a continuous exponential Lévy model

The BS model is one of the most popular and famous pricing models in finance. This was pioneered by Black and Scholes (1973). The BS model is the only exponential Lévy model with a pure continuous sample path and the model is completely driven by Brownian motion that is normally distributed. Thus, the CF of the log -asset price is based on equation2.7. This model is based on some assumptions: The interest rate and the volatility are functions of time which are fixed; no dividend is paid during the life of the option; trading can be done on any number of underlying assets; the transaction cost of hedging a portfolio is nothing; all risk-free portfolios must have the same return; the market is continuous which is not feasible in real-world markets because there are trading hours of stock exchanges. Let S be a function of t such that

∆S

S (2.22)

is the relative change of the underlying asset price at small time-step ∆t. Due to the instantaneous response of the market to new information about an underlying asset price, the change in the asset price from t → t + ∆t is S → S + ∆S. Recall that a market is made up of risk-free and risky assets. Thus, the relative change in an underlying asset price, equation 2.22, is a combination of two different equations:

• The first part involves the underlying asset that would yield return with cer-tainty (has no risk). For example, money kept in a bank account. The param-eter µ (drift) is the average rate of growth of the underlying asset price. It is always fixed:

µ∆t.

• The second part involves the random change in the price of the underlying asset by external factors. The parameter σ (volatility) often denotes the level of uncertainty in the model and ∆B is a sample from the normal distribution with the mean equals to zero and the variance equal to ∆t:

(28)

Thus, equation 2.22can be written as ∆S

S = µ∆t + σ∆B. (2.23)

∆B can be written as B√∆twhere B is a standard normal distribution with mean zero and variance 1. Thus, ∆B2_{→ ∆t} _{as ∆t → 0.}

We can now use the Itô lemma on the asset’s price (equation 2.23) and some assumptions based on the portfolio-hedging strategy to derive the Black-Scholes PDE equation: ∂V (St, t) ∂t + rS ∂V (St, t) ∂S + 1 2σ 2_S2∂2V (St, t) ∂S2 − rV (St, t) = 0. (2.24)

For detail on the derivation of the Black-Scholes PDE, see Hull (2007). Equation 2.24may have many solutions depending on the terminal condition. V (St, t)can be

a call or put option. For a European call and put option, the terminal conditions are given by the payoff function;

Ψc:= (ST − K)+ Ψp := (K − ST)+.

The Black-Scholes partial differential equation for early exercise options (Bermu-dan and American options) does not satisfy the earlier proposed assumptions because there exists no arbitrage argument for Bermudan options unlike European options. Thus, arbitrage opportunities can occur. The equality condition in equation 2.24is now inequality: ∂V (St, t) ∂t + rS ∂V (St, t) ∂S + 1 2σ 2_S2∂2V (St, t) ∂S2 − rV (St, t) ≤ 0. (2.25)

2.3.2 Normal inverse Gaussian as an infinite exponential Lévy model

The normal inverse Gaussian process was firstly introduced and proposed for mod-elling financial data by Barndorff-Nielsen (1997). Since then, this model has been widely studied in much literature for modelling the returns of asset prices and pricing of options under different frameworks, Rydberg(1997),Eberlein and Prause(2000), Webber and Ribeiro (2003), Kéllezi and Webber (2004), H.Albrecher and Predota (2004). The NIG is an infinitely divisible process with stationary independent incre-ments, and an infinite variation process with stable-like small jumps. This model is a variance-mean mixture representation of the normal inverse Gaussian distribution with the inverse Gaussian as the mixing distribution. As a result of the mixture representation of the NIG distribution, we find that the NIG Lévy process Xt may

be represented via a time change of Brownian motion as Xt:= X(t, α, β, δ) = µt + BLt,

(29)

Figure 2.1: Sample path of a stock driven by geometric Brownian motion. Parame-ters used: T = 1, r = 0.02, and σ = 0.7.

where Lt is the inverse Gaussian Lévy process with a scaling parameter δ and

p

α2_{− β}2, for α is the steepness parameter and β is the drift parameter, and B t

is the Brownian motion with the drift parameter and diffusion coefficient 1, and µ ∈ R Barndorff-Nielsen(1997). This implies NIG is a time-changed Brownian mo-tion and the time change t may be chosen as an inverse Gaussian process independent of the directing Brownian motion Geman(2002). The Lévy measure can be written as

νN IG= exp(βx)

δα

π|x|K1(α|x|)dx (2.26)

where α > 0, −α < β < α, δ > 0 and Kρ(x) is the modified Bessel function of the

third kind with index ρSchoutens(2003). The NIG is described by the characteristic triplet (γ, 0, ν), where γ = 2γα π Z 1 0 sin h(βx)K1(αx)dx,

σ = 0 and ν is given by equation 2.26. σ = 0 implies NIG has no Brownian component. Thus, the process is a pure jump process.

2.3.3 Variance-gamma model as an infinite exponential Lévy model

The variance-gamma process was introduced and firstly proposed to model the dy-namics of an underlying asset price byCarr and Madan.(1998). In their paper, they describe the risk-neutral dynamics of asset prices driven by the VG process, and also derive a closed-form formula for the return density and for pricing European options. Afterwards, this model has been investigated and applied to price different types of

(30)

Figure 2.2: Sample path of Normal inverse Gaussian process: Parameters used: T = 1, σ = 0, θ = 4 and α = 0.05.

options by several authors: Madan and Seneta (1990),Carr et al.(2002),Hirsa and Madan (2004),F.Hubalek et al. (2006),Almendral and Oosterlee(2007).

The VG model is a three-parameter stochastic process that models the dynamics of the logarithm of an underlying asset price. θ is the drift term, σ is the volatility and κ is the variance rate. The VG process is obtained by evaluating Brownian motion with constant drift and volatility at a random time change given by the gamma process, that is, replacing the time in Brownian motion with gamma process Madan and Seneta (1990). Thus, the VG model is characterised with discontinuous sample path. Let

β(t; θ, σ) = θt + σBt

be a Brownian motion with drift θ, volatility σ and a standard Brownian motion Bt.

Also, there exists a process of independent gamma increments over non-overlapping intervals of time (t, t + h), that is, a gamma process γ(t; µ, κ) with mean rate µ and variance rate κ, Carr and Madan. (1998). Besides the volatility and the drift of the time-changed Brownian motion, there exists an additional parameter in this model that controls kurtosis and skewness. The θ controls the skewness while the κ controls the kurtosis. This is the variance of the gamma-distributed time. Under the VG process, the unit period continuously-compound return is normally distributed, conditional on the realisation of a random time that has a gamma density. The VG process X(t; θ, σ, κ) in terms of the Brownian motion, with drift β(t; θ, σ) and gamma process with unit mean rate γ(t; 1, κ, ) can be written as Carr and Madan. (1998);

X(t; θ, σ, κ) = β(γ(t, 1, κ), θ, σ).

The VG process has no Brownian motion and it can be written as the difference of two independent gamma processes. With these properties, the Lévy density can be

(31)

expressed as νV G(dx) =        C exp(G|x|) |x| dx for x < 0, C exp(−M |x|) |x| dx for x > 0 where C = 1 κ, G = r θ2_κ2 4 + σ2_κ 2 − θκ 2 !−1 , M = r θ2_κ2 4 + σ2_κ 2 + θκ 2 !−1 .

The Lévy triplet of the VG model is given by (γ, 0, νV G(dx)) where

γ = −C(G exp(−M ) − 1) − M (exp(−G) − 1)

M G .

The drift parameter θ measures the directional premium since it affects the skewness of the process Carl et al.(2003). The parameter C = 1

κ controls the overall activity

rate of the process while κ = 1

C controls the kurtosis of the model. The parameters

G and M control the rate at which arrival rate reduces with the size of the move. The smaller the value of κ, the more the sample path of the VG model resemble the sample path of a typical stock-price. When θ > 0, then G > M and this results in a positive skewness in the distribution of the VG model. The reverse holds when θ < 0. When θ = 0, then G = M and this makes the distribution of the VG process to be symmetric.

(32)

Figure 2.3: Sample path of a Variance Gamma process on a fixed grid. Parameters used: T = 1, σ = 0.3, θ = 4 and κ = 0.05.

The direct computation of CF of the log-asset price of the NIG model and the VG model is not trivial. The use of a subordinator makes the computation of these CFs possible. This is as a result of the conditional Gaussian structure of the normal in-verse process and the variance-gamma process which also simplifies computations and simulations significantly. A subordinator (St)t≥0 is almost surely a non-decreasing

Lévy process such as time, (a sequence with non-negative increments). Thus, a neg-ative jump does not exist. The Lévy triplet can be expressed as (γs, σs, νs) where

σs= 0. That is, Sthas no diffusion parameter. Sthas a positive drift and jump of

fi-nite variation. The increasing nature of a subordinator qualifies it to be represented as a time indicator. The transformation of an old Lévy process into a new Lévy process with the new characteristic triplet given in terms of the old Lévy process, is termed Subordination, Cont and Tankov (2004) and this is illustrated by Theorem 2.3.1.

Theorem 2.3.1 (Subordination of Lévy processes). Let (Xt)t≥0 be a Lévy process

on R with the Lévy triplet (γ, σ, ν) and its CF ΦX(ω) and characteristic exponent

ψX(ω) as given in equations (2.19) and (2.20). Let (St)t≥0 be a subordinator with

triplet (γs, 0, νs) and Laplace exponent Ls(ω). We assume (Xt)t≥0 and (St)t≥0 are

independent. Let (Y_t)t≥0 be a stochastic process by random time changing the old

Lévy process (X_t)t≥0 such that;

Yt= XSt

is a Lévy process and its CF Φ_Y(ω) is given by;

(33)

Thus, the CF of Yt is derived by substituting the characteristic exponent ψX(ω) of

the old Lévy process into the Laplace exponent of the subordinator L_s(ω). The char-acteristic triplet of Yt, (γy, σy, νy) is given by

γy = γsγ + Z ∞ 0 νs(ds) Z 1|x|≤1xPxt(dx), σ2_y = γs2σ2, νy(dx) = γsν(dx) + Z ∞ 0 Pxt(dx)νs(ds),

where Pxt is the probability distribution of Xt. (Yt)t≥0 is said to be subordinate to the

process (X_t)t≥0 Cont and Tankov (2004).

For more information on the types of subordination, see Appendix A.1.1. Sato (1999) give a detailed proof of Theorem 2.3.1 and Cont and Tankov(2004) give an outline of discussion of the theorem.

Below are the CF’s of the BS, NIG and VG models1 _{Barndorff-Nielsen} ₍₁₉₉₇_),

Madan and Seneta (1990).

BS model Φ(ω, t) = exp(iωµt − 1₂σ2ω2t) NIG model Φ(ω, t) = exp(iωµt − 1

2σ 2_ω2_{t) · %} N IG(ω; t; α, β, δ) %N IG(ω, t, α, β, δ) = exp[δt( p α2_{− β}2_{) − (pα}2_{− (β + iω)}2_)] VG model Φ(ω, t) = exp(iµωt) · %V G(ω; t; ω, κ, θ) %V G(ω; t; θ, κ) = (1 − iωθκ +1₂σ2κω2)− t κ

Table 2.1: The characteristic functions of the log-asset price of the BS, NIG and VG models, where %N IG and %V G are the characteristic exponents.

Knowing the CFs of the BS, NIG and VG models, we can easily compute their cumulants using equation 2.9. Table2.2 presents the cumulants of these models.

BS model NIG model VG model

c1 = (r −1₂σ2)t c1= (r −₂1σ2+ w)t + δtβ(α2− β2)−1/2 c1= (µ + θ)t c2 = σ2t c2= δtα2(α2− β2)−3/2 c2= (σ2+ κθ2)t c3 = 0 c3= 3δtα2β(α2− β2)−5/2 c3= 3(σ2θκ + 2θ3κ2)t c4 = 0 c4= 3δtα2(α2+ 4β2)(α2− β2)−7/2 c4= 3(σ4κ + 2θ4κ3+ 4σ2θ2κ2)t w = δ(pα2_{− β}2_{) − (pα}2_{− (β + iω)}2₎ _{w =} 1 κln(1 − θκ − σ2_κ 2 )

Table 2.2: The cumulants of the log-asset price of the BS, NIG and VG models, where w is the drift correlation term which satisfies e−wt _{= Φ(−i, t)}_.

1_{µ is the drift in BS and NIG models, σ is the volatility, t is the time to maturity, α =} µ σ2 and β =

r µ2+2σ2

µ

(34)

2.3.4 Complex discontinuities in characteristic functions

In financial modelling, CFs of the log-asset price of many option pricing distributions involve multivalued functions. These functions may introduce discontinuities into the CFs of the distributions of financial models, which leads to incorrect option prices, especially if the options are priced by Fourier inversion methods Lord (2008). An alternative to avoid discontinuity is to integrate the system of ODEs that evolves under the risk-neutral measure of the underlying asset of the model by numerical methods, as this will result in a continuous solution Lord and Kahl(2010).

From a computational perspective, Fourier methods are advantageous. Therefore, it is worthwhile to investigate how to do away with discontinuity in the solution. Recall the log of a complex variable z = x+iy = rei(arg(z)+2πn)_{with arg(z) ∈ [−π, π),}

n ∈ Z, and r ∈ R can be expressed as

log z = log r + i(arg(z) + 2πn), where r is the radius. (2.28) If we restrict the log z to its principal branch, that is, setting arg(z) to be the principal argument, arg(z) ∈ [−π, π), and n = 0, then the branch cut of the complex logis (−∞, 0] and the complex log is discontinuous along it,Kahl and Jäckel(2005). Thus restricting log z to its principal branch yield discontinuities in the CF. Hence, we have wrong option prices if we use the discontinuous CF. The restriction on arg(z) ∈ [0, π) implies we cut the complex plane along the negative real axis.

In a popular model of asset-price dynamics known as the Heston model, the log of complex arguments in the CF gives rise to numerical instability, and as a result of this, most implementations in the Heston model are not robust for long dated options, Kahl and Jäckel(2005). To show that exponential Lévy models may have complex discontinuities, we illustrate the effect of the principal branch on the NIG and VG models.

Example 2.3.2 (Variance-gamma model). Let the expression for an underlying asset modelled by the VG process be given by

St= Ptexp(wt + θGt+ σBGt) (2.29)

where Btis a standard Brownian motion, Gtis the gamma process with the variance

rate ν > 0, Pt is the forward price of St, σ > 0 and the parameter w (see Table

2.2) is chosen so that the expectation of the exponential part of equation 2.29 is a unit. For simplicity, let pt = ln Pt and ˆpt = pt+ wt. Recall the conditional CF of

the distribution of log-asset returns under the VG model ΦV G =

exp(iω ˆpt)

(1 − iωθκ +1₂σ2_κω2₎κt

. (2.30)

The ςth _{moment of the underlying asset exists as long as the minimum ς}

− and the maximum ς+ is defined by ς± = − θ σ2 ± r θ2 σ4 + 2 νσ2

(35)

so that the extended CF in equation 2.30is well-defined for ω ∈ ∧x which is defined

by {ω ∈ C| − I(ω) ∈ (ς−, ς+)}. The bounded maximum and minimum moments

does not depend on the maturity time T since the VG model is completely time-homogeneous. The multivalued function is present in

(1 − iωθκ +1 2σ

2_κω2_). _(2.31)

If we restrict the multivalued function to its principal branch, the CF of the dis-tribution of the VG model will only be continuous if (1 − iωθκ + 1

2σ2κω2) does

not cross the negative real line. To satisfy this condition, the imaginary part of (1 − iωθκ + 1₂σ2κω2) must be positive. This can only occur when ω = x + iy, with x = 0 or y = _σθ2. Then 1 − iωθκ +1 2σ 2_κω2 _{= 1 − ixθκ + θκy +}1 2σ 2_κx2₊1 2σ 2_{κ · 2ixy −} 1 2σ 2_κy2_. _(2.32)

For x = 0, equation 2.32 becomes 1 − ixθκ + θκy +1 2σ 2_κx2₊1 2σ 2_{κ · 2ixy −} 1 2σ 2_κy2_{= 1 + θκy −}1 2σ 2_κy2_, _{for θ > 0.} (2.33) For y = θ σ2, 1 − ixθκ + θκy + 1 2σ 2_κx2₊1 2σ 2_{κ · 2ixy −} 1 2σ 2_κy2_{= 1 +} κθ2 2σ + 1 2σ 2_κx2 > 1, for κ, θ > 0 (2.34) so that (1 − iωθκ + 1 2σ

2_κω2₎ _{can never be negative. Now, we can say the principal}

branch of the complex power function is the right one, as this is the only one that gives real values for equation2.31. Thus, the CF of the distribution of the VG model is continuous when the complex log is restricted to the principal branch despite the presence of the multi-valued function in the CF. Thus, the principal branch of a multivalued function works well for the VG model Lord and Kahl (2010).

2.4 Risk-neutral valuation formula

Pricing financial instruments is a subject which lies at the heart of Financial Math-ematics. When valuing financial instruments, a risk-neutral framework must be used to avoid arbitrage opportunities. An arbitrage opportunity is a self-financing portfo-lio which admits a terminal gain, without any probability of intermediate or terminal loss. In order to price options in our models, the model must satisfy the fundamen-tal theorems of asset pricing. The First Fundamenfundamen-tal Theorem of Asset Pricing relates the notion of the existence of equivalent martingale measures (EMM) to an arbitrage-free opportunity, and The Second Fundamental Theorem of Asset Pricing relates the notion of uniqueness of the EMM to market completeness, Harrison and Pliska (1981).

The popularly known option pricing theory of the BS model depends on the fact that the payoff of every contingent claim can be replicated by a self-financing

(36)

portfolio. In other words, the risk of holding or selling an option can be completely replicated against it. Thus, the market is arbitrage-free. In such model, there exists a unique measure Q (risk neutral-measure) which is equivalent to the “ real world measure ” P, and this makes the discounted price process a martingale:

S_s∗ = EQ[S_t∗|Fs],

where

S_t∗= Ste−rt, for 0 ≤ s < t ≤ T.

Recall the expression for the underlying asset price from equation 2.21, St= S0eXt

where Xt is a Lévy process under the real world measure P.

A financial market is risk-neutral whenever the expected return of investment on the underlying asset is the same as the riskless rate of interest. Studying the Black-Scholes PDE, equation 2.24, the function depends on the underlying asset price S and time t, not on the investor’s preferences. And the relevant parameters needed for the valuation of the options are the risk-free rate of interest r and the volatility σ, not the drift term µ. This leads to the work ofCox and Ross(1976) who developed a technique known as the risk-neutral valuation technique . This technique is based on the assumption that under option pricing, investors are indifferent to risk. i.e, the world is risk-neutral.

For a stock-price driven by the NIG model and the VIG model, there are many equivalent measures under which the discounted asset-price process is a martingale because the noise process has jumps of random sizes. A market built on a Lévy process cannot be replicated by a suitable portfolio. Instead, an additional condi-tion must be used to select the appropriate martingale measure from the multiple equivalent martingale measure (EMM). Thus, financial markets driven by Lévy pro-cesses are incomplete, and the EMM is not unique. Only Lévy markets that are built on pure Brownian motion (BS model) or Poisson process are complete, Cont and Tankov (2004). To select the appropriate EMM under exponential Lévy mod-els, six measures have been proposed in literature; the minimal martingale measure, Föllmer and Schweizer (1991), the Esscher martingale measure, Gerber and Shiu (1994), H.Bühlmann et al. (1996), the minimal entropy martingale measure, Miya-hara (1999), the variance optimal martingale measure, Schweizer (1995), the mean correcting martingale measure and the utility martingale measure. Here, we discuss Boyarchenko and Levendorskii (2002) universal approach called the EMM-condition based on the Esscher martingale measure. Since Xt is a Lévy process, using the

characteristic exponent of a Lévy process, %(·), S0 = S0e−t(r+% Q_(−i)) for %Q(0) = 0 and r + %Q(−i) = 0 (2.35)

(37)

is called the EMM-condition. Using the Esscher transform, (see Miyahara (2002)) in terms of the characteristic exponent,

%Q(ω) = %P(ω − iς) − %P(−iς) (2.36) where ς is real, i =√−1and the EMM-condition becomes

r + %P(ω − iς) − %P(−iς) = 0. (2.37)

Thus,

S_t∗ = S0e−t(r+%

P_{(ω−iς)−%}P_(−iς))

, and Q ∼ P.

In Lévy models, we assume under the real-world measure P, the existence of some equivalent martingale measure Q such that the discounted asset price process is a martingale. In the VG model, the three main parameters σ (volatility), ν (variance) and θ (drift), are also chosen to be risk-neutral. The unique pricing formula of a contingent claim is the expectation under the martingale measure of the discounted payoff at the maturity date.

Given the payoff function Ψ(ST) for a European option, the existence of an

equivalent martingale measure Q, not necessarily unique, in an arbitrage-free market, leads to the risk-neutral pricing formula V (t0, S0)at time t0 where ∆t = T − t0 such

that

V (S0, t0) := e−r∆tEQ[Ψ(ST)|F0] (2.38)

where ST is an adapted stochastic process, r is the risk-free rate of return and

F0(F0 = 0) is the filtration generated at time t by the risk-neutral distribution Q.

This expectation makes it possible to value the equation 2.38 by the means of the Monte Carlo methods. Since equation2.38 is an expectation, we represent equation 2.38 as

V (x, t0) := e−r∆t

Z ∞

−∞

Ψ(y)f(y|x; ∆t)dy (2.39)

and f(y|x; ∆t) is the conditional probability density function representing the tran-sition of x at time t0 to y at time T where x and y are any monotonic functions of

the underlying asset’s price. For more information see Harrison and Pliska(1981). By means of the risk-neutral pricing formula, the price of any option except for early exercise options can be written as an expectation of the product of the dis-counted payoff of the option price and the risk-neutral density function. Recall that the probability density function of most Lévy processes is not available in a closed form but can be calculated from the CF by the Fourier inversion method. This leads us to study some basic concepts of Fourier methods.

2.5 Some fundamentals of Fourier methods

Here, we present some mathematical properties of the Fourier methods that are useful to our discussion. For detailed study on Fourier methods in finance, see Cherubini et al. (2010), and for general books in literature on Fourier methods, see Brigham (1974),Boyd (2000),Press et al.(2007).

Fourier methods for pricing early-exercise options under levy dynamics

Tolulope Rhoda Fadina

Declaration

Abstract

Opsomming

Acknowledgements

Dedications

Contents

List of Figures

List of Tables

Chapter 1

Introduction

1.1

Motivation for research

1.2

Organisation of thesis

Chapter 2

2. Models of asset prices

2.1

Introduction

2.2

Definitions and properties of Lévy processes

2.3

Models driven by exponential Lévy processes

2.4

Risk-neutral valuation formula

2.5

Some fundamentals of Fourier methods