**School of Education, Culture and Communication**

**Division of Applied Mathematics**

### Bachelor Thesis in Mathematics / Applied Mathematics

### Option Pricing using the Fast Fourier Transform Method

### by

### Abaynesh Berta

### Kandidatarbete i matematik / tillämpad matematik

### DIVISION OF APPLIED MATHEMATICS

MÄLARDALEN UNIVERSITY**School of Education, Culture and Communication**

**Division of Applied Mathematics**

Bachelor thesis in mathematics / applied mathematics Date:

3rd October 2020 Project name:

Option Pricing using the Fast Fourier Transform Method Author: Abaynesh Berta Supervisor: Milica Ranˇci´c Reviewer: Siyang Wang Examiner: Doghonay Arjmand Course name:

Degree Project in Mathematics Course code:

MMA390 Comprising: 15 ECTS credits

Abstract

The fast Fourier transform (FFT), even though it has been widely applicable in Physics and Engineering, it has become attractive in Finance as well for it’s enhancement of computa-tional speed. Carr and Madan succeeded in implementing the FFT for pricing of an option. This project, inspired by Carr and Madan’s paper, attempts to elaborate and connect the vari-ous mathematical and theoretical concepts that are helpful in understanding of the derivation. Further, we derive the characteristic function of the risk neutral probability for the logarithmic terminal stock price. The Black-Scholes-Merton (BSM) model is also revised including de-rivation of the partial deferential equation and the formula. Finally, comparison of the BSM numerical implementation with and without the FFT method is done using MATLAB.

Keywords: Black–Scholes–Merton model, Characteristic function, Fast Fourier transform, Fourier/Inverse Fourier transform, Option pricing

Acknowledgments

I would like to take the opportunity to sincerely acknowledge the people who contributed to the project. Firstly, I would like to express my special thanks to Milica Ranˇci´c who introduced me to the topic and dedicated her valuable time; thank you for the continuous inspiration, patient guidance and support throughout this project. Secondly, I would like to thank Siyang Wang for the valuable remarks and suggestions. I would also like to thank Doghonay Arjmand and other instructors in the UKK department.

## Contents

1 Introduction 7

1.1 Background and Literature Review . . . 7

1.2 Project Description . . . 9

2 The BSM Model and Fourier Transform 10 2.1 Basic Concepts . . . 10

2.1.1 Probability Theory- and some Mathematical Concepts . . . 10

2.1.2 Option Valuation . . . 12

2.1.3 Characteristic Functions . . . 13

2.2 The Black–Scholes–Merton Model . . . 16

2.3 Fourier Transformation of the BSM Formula . . . 21

2.4 The Fourier Inversion . . . 23

3 Option Pricing Using the FFT 24 3.1 DFT and FFT of the BSM Formula . . . 24

3.1.1 Discretization using the Trapezoidal Rule . . . 24

3.1.2 Discretization using Simpson’s Rule . . . 26

3.2 Numerical Implementation . . . 26

4 Conclusion and Future Work 35 4.1 Conclusion . . . 35

4.2 Future Work . . . 35

References 36 Appendix 38 A MATLAB Code 39 A.1 MATLAB Script of the Implementation . . . 39

A.2 MATLAB Function Files used in the Implementation . . . 42

A.2.1 Function "optByBSM_FFT" . . . 42

A.2.2 Function "optSensByBSM_FFT" . . . 45

A.2.3 Function "MyOptSensByFFT" . . . 48

A.2.5 Function "optByBsmNI" . . . 57 A.2.6 Function "optSensByBsmNI" . . . 59 A.2.7 Function "MyOptSensByNI" . . . 62

## List of Tables

3.1 Time elapsed when N = 213. . . 27 3.2 Time elapsed when N = 216. . . 27 3.3 Comparison of FFT using Simpson and FFT using Trapezoidal verses call

prices obtained using a numerical method without the FT at N = 213. . . 28 3.4 Comparison of FFT using Simpson and FFT using Trapezoidal verses call

prices obtained using a numerical method without the FT at N = 216. . . 29 3.5 Comparison of FFT call prices with the call prices obtained using blsprice

at N = 213. . . 29 3.6 Comparison of FFT call prices with the call prices obtained using blsprice

## List of Figures

3.1 Plot of strike price vs call price at different values of α. . . 30 3.2 Plot of log strike price vs call price at different values of α. . . 31 3.3 Option price surface: (X, Y, Z)=(Time to maturity, Strike, Call) when α = 3. . 31 3.4 Option price surface. (X, Y, Z)=(Time to maturity, Log strike, Call) when α = 3. 32 3.5 Option price surface: (X, Y, Z)=(Time to maturity, Strike, Call) when α = 8.5. 32 3.6 Option price surface. (X, Y, Z)=(Time to maturity, Log strike, Call) when

α = 8.5. . . 33 3.7 Option price surface. (X, Y, Z)=(Time to maturity, Strike, Call) when η =

0.0002. . . 34 3.8 Option price surface. (X, Y, Z)=(Time to maturity, Strike, Call) when η = 1.25. 34

## Abbreviations and Symbols

ATM At–The–Money ITM In–The–Money OTM Out of–The–Money BSM Black–Scholes–Merton

CDF Cumulative Distribution Function CF Characteristic Function

FT Fourier Transform

DFT Discrete Fourier Transform FFT Fast Fourier Transform GBM Geometric Brownian Motion PDF Probability Density Function RV Random Variable

SDE Stochastic Deferential Equation PDE Partial Deferential Equation Φ Normal CDF

φ Characteristic function

## Chapter 1

## Introduction

While dealing with option pricing, accuracy and efficiency are crucial matters because the analysis involves large data manipulation. There are plenty of models available, each with their advantages, disadvantages and different levels of complexity. There are also different approaches to the derivation and implementation of the models. We intend to discuss a number of alternative models. However, we particularly focus on the fast Fourier transform approach applied to the Black–Scholes–Merton (BSM) formula [4].

### 1.1

### Background and Literature Review

An option is a derivative that gives the buyer of the contract the right to carrying out a trans-action within a pre-specified date, called maturity T at agreed exercise price, called a strike price K, where a derivative is a financial instrument whose value is derived from the values of other underlying variables [12]. Options are of two types; an American option which can be exercised at any time before the maturity and a European option which can be exercised only at the maturity. Both the American and the European options are divided into a call- and a put option. A call option gives the holder of the contract the right to buy the underlying asset while a put option gives the right to sell it. This paper considers the European call option only. The buyer of an option is said to take a long position while the seller takes a short position.

Estimating the price of an option helps the actors in financial market to adjust their strategy accordingly. Pricing models are financial tools that enable estimating the option price given existing inputs. Among a number of models available, we discuss some of them in this section. To compute an option price, we need to approximate the movement of the stock price, denoted by S = {St: 0 ≤ t ≤ T }. However, one can not tell before hand whether S will go up

or down with confidence. We call the part of the return which can be determined before hand deterministicwhile the unanticipated part non–deterministic. Assuming that the instantaneous rate of return of a security is expressed by a sum of a deterministic trend and unexpected vibrations, Kijima [13] claims that modeling the impact of such noise by a standard Brownian motion (see Definition 7, Sec. 2.1.1) would be plausible because the intention is to know the marginal change in the price.

model for a non-dividend paying European call option assuming that the stock follows a con-tinuous random walk with constant variance rate [4]. Regarding the market, the model as-sumes that transactions are made continuously in time, no transaction costs are involved, short selling is allowed and it is possible to buy/sell any fraction of a stock. Further, it assumes a constant short-term interest rate. A related derivation executed by Merton, [20] is applied on warrants. In a later publication, Merton showed that the BSM formula applies even in cases when the interest rate is stochastic, when there is dividend payment, and for the American option [21]. Because the BSM formula is a pioneer in deriving the option price in a closed form, it serves as a reference and most of the models developed later focus on relaxation of the assumptions of the BSM and on broader applications. Even though the assumptions in the BSM help simplifying the analysis, they are restrictive and do not characterize the stock price dynamics realistically.

The assumption of a stock return following continuous sample path is relaxed in the later publications. In [22], for example, Merton derived a more general option pricing formula taking into consideration both continuous and jump processes. In this model, the continuous sample path is modeled by the geometric Brownian motion while the jumps are modeled by a Poisson process. David Bates [3] came up with a more sophisticated model for pricing American options with stochastic volatility and jump-diffusion processes under systematic jump and volatility risk.

Among models applying stochastic volatility, Heston [9] derived the European call price in terms a characteristic function. The stochastic volatility in this model follows Ornstein Uh-lenbeckprocess [9]. In addition to the volatility, Heston assumes a stochastic interest rate and having correlation between the asset return and the volatility. The assumption of correlation between the volatility and the asset return can result in a different kurtosis and skewness than the log normal, hence the need for a characteristic function. Similarly, Backshi [2] developed a closed-form pricing formula that he claimed leads to useful analytical hedge ratios, allow-ing stochastic volatility, stochastic interest rates and random jumps. Schoutens, Simonsy and Tistaertz [24] consider as well a non-normal returns and stochastic volatility. They specific-ally apply a stochastic time on Heston stochastic volatility both with and without jumps, the Barndor Nielsen-Shephard model and Lévy models.

In search for a better accuracy and efficiency, transformed versions of already existing models have been developed, e.g. in [9], [5], [14], [15]. Computational time wise, Carr and Madan [5] used the FFT method which they claim is considerably faster than many of the methods thus far. They specifically recommend the FFT method whenever the characteristic function of the underlying uncertainty is available in the closed form. Carr and Madan analyt-ically derived FFT of the European call price as a function of logarithmic strike in terms of the characteristic function of the log of the terminal stock price. The process involves modifying the European call by multiplying it with a damping factor to ensure the integrability of the Fourier transform (FT). Then, they apply an inverse FT and divide by the damping factor to retrieve the original call price back. This project is also inspired by [5] and will attempt to elaborate and connect the various theoretical and mathematical aspects related to the method. Lee [14] generalizes Carr and Madan’s work and incorporates it with the later publications on the topic and extends to account stochastic interest rates. Further, they develop error bounds and error minimization strategies associated with the discretization of the transformed option

pricing to numerous strikes at ones. Lord and Kahl [see 15, and references therein] have proved an algorithm using the principal branch of the complex logarithm of jump–diffusion models, unlike other practitioners that use branch switching to eliminate the discontinuities in the characteristics function which can lead to wrong computational results. Despite some modifications that have been done after Carr and Madan, the FFT method still appears to be a significantly efficient method.

### 1.2

### Project Description

Primary purpose of the project is pricing an option. Among a number of pricing models available, we want to focus on the BSM partly because of it’s relative simplicity and on the other hand it lets us explore fundamental concepts in financial engineering. We want also to experiment the computational speed of FFT method in pricing the European call price in compassion to other methods. Although the use of FFT for option pricing is done previously by several researchers, e.g. [5], [15] and [14], we have added more details to the derivations and have attempted to bring all relevant aspects together so that it becomes easier for readers of any level to grasp the concepts.

In the next Chapter, in basic concepts section, definitions to various probability, mathem-atical and option valuation concepts that are necessary to the topic are presented. Moreover, the proof of characteristic function to risk neutral probability of the logarithmic terminal stock price is done. The second section revises the BSM model and attempts to prove the BSM PDE and the BSM formula. In the third section, we apply a Fourier transform to a modi-fied European call price by a damping factor as recommended by Carr and Madan [5]. In the last section of Chapter 2, we retrieve the original call price using the inverse FT. Further, conditions for square integrablity are discussed in this last section.

Discretization of the inverted FT, using both trapezoidal and Simpson rules, is done in the first section of Chapter 3. Moreover, the FFT algorithm is examined and FFT of the call price is derived in this section. In the second section of this chapter we experiment numerical implementation of the FFT algorithm for the European call prices. Then we compare the result with a MATLAB built-in function and call prices obtained using numerical integration without the transform method. Finally, the fourth Chapter concludes our work and gives indications of possible future studies.

## Chapter 2

## The BSM Model and Fourier Transform

### 2.1

### Basic Concepts

Before jumping to the main topic of the paper, here we want to define some terminology related to probability and option pricing. The definitions are entirely taken from text books, mostly from the book “Stochastic Processes with Applications to Finance” [13].

### 2.1.1

### Probability Theory- and some Mathematical Concepts

Definition 1. Given a sample space Ω and a σ -field F, if a set function P defined on F satisfies i) P(Ω) = 1,

ii) 0 ≤ P(A) ≤ 1 for any event A ∈ F, and

iii) for mutually exclusive events An∈ F, n = 1, 2, ..., that is, Ai∩ Aj= /0 for i 6= j,

P(A1∪ A2∪ · · · ) = ∑∞n=1P(An),

we call P a probability measure (or probability for short). The triplet (Ω, F, P) is called prob-ability space[13].

σ -field is family of events that incorporates all subset events in Ω and the compliments of the events. See [13] for the formal definition.

Definition 2. [13] Given a probability space (Ω, F, P), let X be a mapping from Ω to an interval I of the real line R. The mapping X is said to be a random variable if, for any a< b, {ω : a < X (ω) ≤ b} ∈ F. In this case, X is said to be F − measurable.

When X gets a value, we call it realization and the probability of realization of X is in the interval (a, b], P{a < X ≤ b} is given by:

P{a < X ≤ b} = Z b

a

where f (x) is called probability density function (PDF) of X . While the cumulative distribu-tion funcdistribu-tion (CDF) F(x) defined as P(X ≤ x), x ∈ I, is given by:

F(x) = Z x

−∞

f(x)dx.

Definition 3. A family of random variables {X (t);t ∈ T } (or X (t) for short) parameterized by time t ∈ T is called a stochastic process [13].

Remark1. In most of the equations in this paper, stochastic process is denoted with a subscript as in the notation Xt for clarity.

Definition 4. An increasing family of σ -fields (Ft)t∈[0,T ]: ∀t ≥ s ≥ 0, Fs⊆ Ft⊆ F, is called a

filtration or information flow_{on (Ω, F, P). F}_{t} can be interpreted as the information known at
time t and it increases as time progresses. [19]

Definition 5. [13] For a stochastic process X (t), its quadratic variation is defined, if it exists, as

[X , X ](T ) = lim

N→∞|X(t

N

i ) − X (ti−1N )|2,

where the limit is taken over all partitions such that for each N,
0 = t_{0}N < t_{1}N < · · · < t_{i−1}N < t_{N}N= T and δN= maxit_{i−1}N − t_{i}N as N → ∞.

Definition 6. The conditional expectation of X under the event {Y = y} is given by: E[X |Y = y] =

Z ∞ −∞

xdF(x|y), y_{∈ R,}
provided that R∞

−∞|x|dF(x|y) < ∞, where F(x|y) = P{X ≤ x|Y = y} denotes the conditional

distribution function of X under {Y = y} [11].

Definition 7. [13] A stochastic process z(t),t ≥ 0, defined on the probability space (Ω, F, P) is called a standard Brownian motion if:

• it has independent increments,

• the increment z(t + s) − z(t) is normally distributed with mean 0 and variance s, inde-pendently of time t, and

• it has continuous sample paths and z(0) = 0.

z(t) ∼ Φ(0, T ), √z(t)

T ∼ Φ(0, 1),

where Φ(∗) represents the normal distribution, and the normal PDF has the following form: f(x) = 1 σ √ 2πe −1 2( x−µ σ ) 2 . (2.1)

Additionally, the Brownian motion has the following properties [13]: it is a continuous function of time t, it is not differentiable at any point in time, it has infinite variation on any interval no matter how small the interval is; and it has quadratic variation on [0,t] equal to t, for any t ∈ [0, T ], where T is the time domain. Because of these properties, the integration of z(t) is beyond the usual methods of calculus. Hence, we need a stochastic calculus, namely Itô’s lemma [10] which is the stochastic counter part of a chain rule in deterministic calculus. Definition 8 (Itô’s lemma). Let G(St,t) be a twice differentiable function of t and of the

random process St (a stock price in our case) given by the stochastic defferential equation

(SDE):

dS_{t}= µStdt+ σ Stdzt, (2.2)

where µ and σ are rates of drift (the stock’s expected rate of return) and diffusion (the volatility
of the stock price) parameters, respectively. Then the Itô’s lemma (formula), e.g. in [12], is
defined as:
dF_{t}= ∂ G
∂ St
dS_{t}+∂ G
∂ t dt+
1
2
∂ G
∂ S2_{t} σ
2_{S}2
tdt.

Substituting equation (2.2) and collecting the dt coefficients together, we get:
dG=
∂ G
∂ St
µ St+
∂ G
∂ t +
1
2
∂2G
∂ S_{t}2σ
2_{S}2
t
dt+∂ G
∂ St
σ Stdzt. (2.3)

More explanations regarding Itô’s calculus and proofs can be found e.g. in [10].

### 2.1.2

### Option Valuation

Because an option gives a right but not an obligation, the holder of a call option is likely to exercise if the stock price ST at maturity T is greater than the strike price K. In this case, the

option is said to be in–the–money (ITM) where the reverse case is called out of–the–money (OTM). Thus, the payoff function w of a European call option is given by:

w(S,t) = {ST− K}+. (2.4)

The expression {∗}+denotes the maximum of (∗) or zero [13].

From (2.4), we can see that the payoff function is always positive. Therefore, the writer of an option contract charges the purchaser a fee called an option premium for taking the risk. Having the information about ST, the price of the option might be obvious. However, because

the stock price varies in time randomly, it is a stochastic process. Therefore, the exact value of ST is not known in advance. Hence, the discounted expected payoff function under a risk

neutral probability is considered, see [13] and [7] for example:

C= e−rTEQ[{ST− K}+], (2.5)

where C is the price of the call option, E denotes the expected value, Q is a risk neutral probability measure and r is a risk neutral interest rate.

The risk neutral rate and probability are needed to avoid arbitrage opportunities and to keep the market in equilibrium. An arbitrage opportunity is a possibility of making profit without assuming any risk. The existence of such opportunity might attract investors to take same strategies and would affect the stock price, which in turn results an imbalance in the supply and demand market equilibrium.

The dynamics of a stock price is modeled by the the following SDE:

dS= µSdt + σ Sdz, (2.6)

where S is the underlying stock price, µ is the expected value of S , σ is the volatility of S and zis a Wiener process.

If we consider a function G(St,t) = ln S, then ∂ G_{∂ S} = 1_{S}, ∂

2_{G}

∂ S2 = − 1

S2, ∂ G_{∂ t} = 0. Using the Itô’s
lemma, we rewrite it as:

dG= µ −σ 2 2 dt+ σ dz. (2.7)

When we integrate equation (2.7) over time t = 0 to t = T , we see that ln ST is normally

distributed with mean ln S0+

µ −σ 2 2 T and variance σ2T: ln ST− ln S0∼ Φ µ −σ 2 2 T dt, σ2T dz ln ST ∼ Φ ln S0+ µ −σ 2 2 T dt, σ2T dz ,

where S0 and ST are stock price at time moments t = 0 and t = T , respectively. Given the

normal PDF in equation (2.1), we write the risk-neutral PDF of ln ST as follows:

Q(ln ST)= 1 σ √ 2πT e − ln ST − ln s0+ µ −σ 2 2 T 2 2σ 2T .

Letting sT = ln ST, we can also writeQ as:

Q(sT) = 1 σ √ 2πT e − sT − s0+ µ −σ 2 2 T 2 2σ 2T . (2.8)

### 2.1.3

### Characteristic Functions

In option valuation and in most problems involving uncertainty, the distribution function of the random variable in question plays a significant roll. The characteristic function (CF) of a random variable depicts the distribution function uniquely.

Definition 9. The characteristic function φX(u) of a real valued random variable X is defined

as: φX(u) =E[eiuX], i =

√ −1.

By the definition of the expected value, φX(u) can be written as: φX(u) = Z ∞ −∞ eiuXdF(x) = Z ∞ −∞ eiuXf(x)dx, where F(x) and f (x) denote CDF and PDF of X respectively.

Let us introduce a lemma which will be used in the prove of a theorem stated right after it. Lemma 1. Let y = ˆµ + iu ˆσ2, then

ius_{T} −(sT− ˆµ )
2
2 ˆσ2 =
−(sT− y)2+ 2 ˆµ ˆσ2iu− u2σˆ4
2 ˆσ2 .
Proof.
RHS= −[s
2
T− 2sT( ˆµ + iu ˆσ2) + ( ˆµ2+ 2 ˆµ iu ˆσ2− u2σˆ4)] + 2 ˆµ ˆσ2iu− u2σˆ4
2 ˆσ2
= (−s
2
T+ 2sTµ + 2sˆ Tiu ˆσ2− ˆµ2− 2 ˆµ iu ˆσ2+ u2σˆ4) + 2 ˆµ ˆσ2iu− u2σˆ4
2 ˆσ2
=−s
2
T + 2sTµ + 2sˆ Tiu ˆσ2− ˆµ2
2 ˆσ2
LHS= iusT −
(sT− ˆµ )2
2 ˆσ2 =
2 ˆσ2iusT− (s2T− 2sTµ + ˆˆ µ2)
2 ˆσ2
= −s
2
T+ 2sTµ + 2sˆ Tiu ˆσ2− ˆµ2
2 ˆσ2 = RHS.

Theorem 1. The CF of the risk neutral probability density of the logarithmic terminal spot price of an underlying asset is expressed as:

φsT(u) = e i{s0+(µ−12σ2)T}u− (σ 2T )u2 2 , where s= ln S.

Proof. We start the proof with the definition of a characteristic function:

φsT(u) =E[e
iu sT_{] =}
Z ∞
−∞
eiu sT Q(s
T)dsT
=
Z ∞
−∞
eiu sT 1
σ
√
2πT e
−
sT −
s0+
µ −σ_{2}2
T
2
2σ 2T ds_{T}
= 1
σ
√
2πT
Z ∞
−∞
eiu sT−
sT −
s0+
µ −σ
2
2
T
2
2σ 2T ds_{T}.

Let ˆµ =
h
s_{0}+µ −σ
2
2
T
i
and ˆσ = σ
√
T, then
φsT(u) =
1
ˆ
σ
√
2π
Z ∞
−∞
eiu sT−{sT − ˆµ}
2
2 ˆσ 2 ds_{T}.

Now we apply Lemma 1:
φsT(u) =
1
ˆ
σ
√
2π
Z ∞
−∞
e−(sT −y)
2+2 ˆµ ˆσ2iu−u2 ˆσ4
2 ˆσ 2 ds_{T}
= e
2 ˆµ ˆσ2iu−u2 ˆσ4
2 ˆσ 2 1
ˆ
σ
√
2π
Z ∞
−∞
e−(sT −y)
2
2 ˆσ 2 ds_{T}
= eµ iu−ˆ 1_{2}u2σˆ2 1
ˆ
σ
√
2π
Z ∞
−∞
e−
sT −y
√
2 ˆσ 2
2
dsT.
Let z = _{√}sT−y

2 ˆσ2, and let us consider a closed upper- and lower bounds of sT, [−γ, γ] for the

integral:
z= sT − ( ˆµ + iu ˆσ
2_{)}
√
2 ˆσ2
= s√T− ˆµ
2 ˆσ2
−iu ˆ√σ
2,
φsT(u) = e
ˆ
µ iu−1_{2}u2σˆ2 1
ˆ
σ
√
2π

### Z

√γ − ˆµ 2 ˆσ 2 −iu ˆ_{√}σ 2 −γ− ˆµ √ 2 ˆσ 2 −iu ˆ

_{√}σ 2 e−z2 √ 2 ˆσ2dz.

Theorem 2 (Cauchy’s integral theorem). [1] If a function f is analytic in a simple connected domain D and Γ is any loop (connected contour) in D, then

Z

Γ

f(z)dz = 0.

We have a connected (rectangle) domain and an analytic integrand e−z2. Thus, we can apply Cauchy’s integral theorem. Let ai, i = 1, 2, 3, 4 denote the corners of the rectangular

contour Γ: a1= γ − ˆµ √ 2 ˆσ2 , a2= −γ − ˆµ √ 2 ˆσ2 , a3= −γ − ˆµ √ 2 ˆσ2 −iu ˆ√σ 2, a4= γ − ˆµ √ 2 ˆσ2 −iu ˆ√σ 2. By the Cauchy’s theorem, we have:

I
z∈D
√
2 ˆσ2e−z
2
dz= 0
⇔
Z a_{2}
a1
√
2 ˆσ2e−z
2
+
Z a_{3}
a2
√
2 ˆσ2e−z
2
+
Z a_{4}
a3
√
2 ˆσ2e−z
2
dz+
Z a_{4}
a1
√
2 ˆσ2e−z
2
= 0.
Now we take each integrals separately and check how they behave as γ → ∞:

lim
γ →∞
Z a_{2}
a1
√
2 ˆσ2e−z
2
dz= − lim
γ →∞
√
2 ˆσ2

### Z

a1 a2 e−z2dz= − √ 2π ˆσ2.Note thatR∞ −∞e−z

2

dz=√π . The proof can be found in e.g. [8, Theorem 2.1.5]. Let us now consider the contour between a4and a1.

lim
γ →∞
Z a_{1}
a4
√
2 ˆσ2e−z
2
dz
≤ lim
γ →∞
√
2 ˆσ2
Z a_{1}
a4
e
−z2
dz
≤ lim
γ →∞
√
2 ˆσ2
Z a_{1}
a4
max
z∈[a4,a1]
e
−z2
dz
≤ lim
γ →∞
√
2 ˆσ2
Z a_{1}
a4
e−
γ − ˆµ
√
2 ˆσ 2
_{
}
dz
≤ lim
γ →∞
e−
γ − ˆµ
√
2 ˆσ 2
_{
}
√
2 ˆσ2
Z a_{1}
a4
dz
= 0.

It can be shown in a similar way for that the contour region between a2and a3is 0. Putting all

together:
lim
γ →∞
Z a_{2}
a1
√
2 ˆσ2e−z2dz+ lim
γ →∞
Z a_{3}
a2
√
2 ˆσ2e−z2dz+ lim
γ →∞
Z a_{4}
a3
√
2 ˆσ2e−z2dz
+ lim
γ →∞
Z a_{4}
a1
√
2 ˆσ2e−z
2
dz= 0
⇔ −√2π ˆσ2+ 0 + lim
γ →∞
Z a_{4}
a3
√
2 ˆσ2e−z
2
dz+ 0 = 0
⇔ lim
γ →∞
Z a_{4}
a3
√
2 ˆσ2e−z
2
dz=
√
2π ˆσ2.

Using this final result, we finalize the proof of the characteristic function for the logarithm of the terminal stock price:

φsT(u) = e
ˆ
µ iu−1_{2}u2σˆ2 1
ˆ
σ
√
2π

### Z

√γ − ˆµ 2 ˆσ 2 −iu ˆ_{√}σ 2 −γ− ˆµ √ 2 ˆσ 2 −iu ˆ

_{√}σ 2 √ 2 ˆσ2e−z 2 dz = eµ iu−ˆ 1

_{2}u2σˆ2 1 ˆ σ √ 2π √ 2π ˆσ2= eµ iu−ˆ 1 2u 2

_{ˆ}σ2 = e h iuns0+ µ −σ 2 2 To−1 2u 2 σ2T i .

### 2.2

### The Black–Scholes–Merton Model

Black, Scholes and Merton claim that the price of the stock is positively related with the value of the option and the option is almost surely to be exercised if the stock price is greater than the strike price. Thus, the current value of the option is approximated as the stock price

minus the price of a pure discount bond having the same maturity date and face value equal to the option’s strike price. Contrarily, if the stock price is less than the strike price the option will not be exercised, that is the option value will be zero. This explanation satisfies the payoff function we stated in equation (2.4). Further, they claim that keeping the stock price unchanged the value of an option decreases as its maturity date approaches. We notice from (2.4) that w(S,t) is bounded by [0, ST]. Graphically, Black and Scholes have shown that the

value of option is more volatile than the stock because w(S,t) lays below the 45◦line, that is where the option price equals the stock price. Additionally, they state that the relative volatility of the option is not constant and it depends on the stock price and the time to maturity.

Under the BSM’s assumptions mentioned in the introduction, the value of the option V is a function of the price of the stock, time and the variables considered to be known constants. Mathematically, it can be written as V = V (St,t, r, σ , K) or simply V = V (St,t) by ignoring the

constant terms [4].

Theorem 3. The BSM partial differential equation (PDE) relates the rate of return in the money market account in infinitesimal period with the change in the option price by eliminat-ing the individual stock price expectation µ. The PDE is given by:

∂V
∂ t +
1
2
∂2V
∂ S2σ
2_{S}2_{+}∂V
∂ SrS− rV = 0,

where V is the value of the option, S is the underlying stock price, σ is the volatility of the stock price and r is a risk free interest rate.

Proof. The stock price dynamic is modeled by the SDE in (2.6), according to the continuous random walk assumption which is proven to be a Brownian motion. Using the Itô’s lemma, dV can be written as:

dV =∂V
∂ SdS+
∂V
∂ t dt+
1
2
∂2V
∂ S2σ
2_{S}2_{dt}
= ∂V
∂ S(µSdt + σ Sdz) +
∂V
∂ t
dt+1
2
∂2V
∂ S2σ
2_{S}2_{dt}_{.}

Collecting the dt terms gives:
dV =
∂V
∂ Sµ S +
∂V
∂ t +
1
2
∂2V
∂ S2σ
2_{S}2
dt+∂V
∂ Sσ Sdz. (2.9)
The coefficient of dt in (2.9) is called a drift and is deterministic while the coefficient of dz is a
stochastic term. The BSM PDE derivation underlies on eliminating the stochastic term which
is a source of uncertainty. For this purpose we use a strategy called delta hedging to counter
balance the risk that would arise due to the stochastic behavior.

Suppose we take a short position1on a call option and we want to hedge it by buying ∆ amounts of stocks. Suppose also that we borrow m units of money-market Mt with a risk free

1_{Taking a short position occurs when an investor sells assets, e.g. shares of a stock, that are not possessed.}

interest rate r to finance the transaction (dMt is defined by dMt = rMtdt, [13]). We drop the

time index in the notations hereafter for simplicity. Then our portfolio will be: Π = ∆S + mM.

Taking the derivative and substituting the dS from (2.6) and dM = rMdt, we get: dΠ = ∆dS + mdM.

= ∆(µSdt + σ Sdz) + mrMdt

= (∆µS + mrM)dt + ∆σ Sdz. (2.10) Now we want to get the value of ∆ that makes stochastic terms (i.e. coefficients of dz) in the total portfolio (i.e. in dV + dΠ) zero:

σ S ∂V

∂ S + ∆σ S = 0 ∆ = −∂V

∂ S. (2.11)

Adding together dV from (2.9) and dΠ from (2.10) and substituting the ∆ from 2.11 gives:
dV+ dΠ =
∂V
∂ Sµ S +
∂V
∂ t +
1
2
∂2V
∂ S2σ
2_{S}2
dt+∂V
∂ Sσ Sdz
+
−∂V
∂ Sµ S + mrM
dt+ ∆σ Sdz
=
∂V
∂ Sµ S +
∂V
∂ t +
1
2
∂2V
∂ S2σ
2_{S}2_{+ −}∂V
∂ Sµ S + mrM
dt+
∂V
∂ Sσ + −
∂V
∂ Sσ
Sdz
=
∂V
∂ t +
1
2
∂2V
∂ S2σ
2_{S}2_{+ mrM}
dt. (2.12)

In this last equation of the total portfolio (2.12), we have a deterministic coefficient only. Thus, the portfolio is risk-less and must instantaneously earn the same rate of return as other short-term risk-free securities [12]. Therefore, the left hand side (LHS) of the total portfolio dV+ dΠ, which can also be written as d(V + Π) will be:

d(V + Π) = (V + Π)rdt = (V −∂V

∂ SS+ mM)rdt. (2.13)

By equating (2.12) and (2.13), we get:
∂V
∂ t +
1
2
∂2V
∂ S2σ
2_{S}2_{+ mrM}
dt = (V −∂V
∂ S
S+ mM)rdt.

Canceling the dt from both sides and rearranging the above expression gives us the Black–
Scholes–Merton PDE:
∂V
∂ t +
1
2
∂2V
∂ S2σ
2_{S}2_{+}∂V
∂ SrS− rV = 0.

In this PDE the expected value of the stock price µ is eliminated, implying that the in-vestors expectation or attitude towards risk do not have an impact in the BSM model.

Theorem 4. Black–Scholes–Merton formula which is a closed form formula for the price of the European call option is given by:

C(S0, T ) = S0Φ(d1) − Ke−rTΦ(d2),
where
d_{1}=
lnS0
K
+ rT
σ
√
T +
σ
√
T
2
d_{2}=
ln
S0
K
+ rT
σ
√
T −
σ
√
T
2 .

Proof. We have equation (2.5) and from the geometric Brownian motion (GBM)2, we have :
S_{T} = S0e(r−

1

2σ2)T +σ z_{.} _{(2.14)}

Equation (2.5) can also be written using an indicator function1_{{∗}} as:
C= e−rTEQ[{S_{T} − K}1_{S}_{T}>K].

Because the expected value is a linear operator, we can split the equation: C= e−rTEQ[ST1ST>K] − e

−rT_{E}Q_{[K}_{1}
ST>K].

Starting from the most left part of the equation, we can factor out K from the expected value since it is a constant. Then, we will leave with the indicator function in the expected value which is the probability of finishing ITM [13]:

Ke−rTEQ[1ST>K] = Ke

−rT_{Q(S}

T > K).

Substituting ST from equation (2.14), taking the natural logarithm of both ST and K and

re-arranging, gives us:

Ke−rTQ(ST > K) = Ke−rTQ
S0e(r−
1
2σ
2_{)T +σ z}
> K
= Ke−rTQ
ln S0+ (r −
1
2σ
2_{)T + σ z > ln K}
= Ke−rTQ z> ln K − ln S0− (r −
1
2σ
2_{)T}
σ
!
.
2_{The GBM of S}

Tcan be proven by applying Itô’s formula on the derivative of ln Stand integrating overt [0, T ]

Then, we divide by√T to get the standard normal:
Ke−rTQ(S_{T} > K) = Ke−rTQ √z
T >
ln K − ln S0− (r −1_{2}σ2)T
σ
√
T
!
= Ke−rTΦ
ln K − ln S0− (r −1_{2}σ2)T
σ
√
T
!
.

Multiplying the expression in Φ(∗) by −1 and rearranging using the logarithm rules, will give us the d2 term of the BSM formula. Note that the multiplication by −1 does not change the

value because the standard normal is symmetric.

Ke−rTΦ
− ln K + ln S0+ (r −1_{2}σ2)T
σ
√
T
!
= Ke−rTΦ
ln
S0
K
+ rT
σ
√
T −
σ
√
T
2
.
Now, let us proceed to the first part:

e−rTEQ[ST1ST>K] = e
−rT_{E}Qh_{S}
0e(r−
1
2σ
2_{)T +σ z}
1ST>K
i
= S0EQ
h
e(12σ
2_{)T +σ z}
1ST>K
i
.
Under a new probability measure ˜Q, by Girsanov’s theorem[13], we get:

S_{0}EQhe−12σ
2_{T}_{+σ z}
1ST>K
i
= S0E
˜
Q_{[}_{1}
ST>K]
= S0Q(S˜ T > K).

Under the new probability measure, the Brownian motion becomes: ˜z = z − σ T . See [13, Section 12.4] for the proof. Now we modify ST to the new Brownian motion then we apply

the same steps as we did for the d2part:

S_{T} = S0e(r−
1
2σ2)T+σ (˜z+σ T )
= S0e(r+
1
2σ
2_{)}_{T}_{+σ ˜z}
S0Q(S˜ T > K) = S0Q˜
S0e(r+
1
2σ
2_{)}_{T}_{+σ ˜z}
> K
= S0Q˜
ln S0+ (r +
1
2σ
2_{)T + σ ˜z > ln K}
= S0Q˜
˜z
√
T >
ln K − ln S0− (r +1_{2}σ2)T
σ
√
T
!
= S0Φ
ln K − ln S0− (r +1_{2}σ2)T
σ
√
T
!

= S0Φ
− ln K + ln S0+ (r +1_{2}σ2)T
σ
√
T
!
= S0Φ
ln
S0
K
+ rT
σ
√
T +
σ
√
T
2
.

The expression in Φ(∗) is the d1term of the BSM formula, which concludes the proof.

### 2.3

### Fourier Transformation of the BSM Formula

Carr and Madan [5] used a Fourier transform to derive a call option price which is defined as a function of the natural logarithm of the strike price. They considered a modified version call price C(k), where k = ln K. The expected discounted payoff give in equation (2.5) is written in integral form and es is substituted in the place of ST, where s = ln ST, and ek is substituted

for K.
C= e−rTEQ[{ST− K}+]
=
Z ∞
K
e−rT(ST− K)Q(lnST)d ln ST
C_{T}(k) =
Z ∞
k
e−rT(es− ek)Q(s)ds.

Definition 10. A sufficient (but not necessary) condition for the existence of the Fourier trans-form of a function g(x) and its inverse is:

Z ∞

−∞|g(x)|dx < ∞,

which is called an integrability condition.

However, this CT(k) given above is not integrable as can be shown in the following lines.

As k → −∞, ek→ 0 and we get:
C_{T}(k) =
Z ∞
−∞
e−rTesQ(s)ds
= e−rTEQ[es]
= e−rTEQ[ST]
= e−rTS0erT = S0.

To tackle this problem, Carr and Madan [5] modified C(k) as follows: cT(k) = eα kCT(k),

Now we check whether ct(k) is integrable:
Z ∞
−∞|cT(k)|ds =
Z ∞
−∞|e
α k_{C}
T(k)|ds
= e−rT
Z ∞
−∞
|eα k_{(e}s_{− e}k_{)}_{Q(s)|ds}
→ 0 as k → −∞.
This last expression is integrable, thus the FT exists.

Definition 11. Given a function g(x), the FT is defined by: F [g(x)](v) = ψ(v) =Z ∞

−∞

g(x)e−ivxdx Then, Fourier transform ψ(v) of cT(k) derived as follows:

F [cT(k)] = ψ(v) =
Z ∞
−∞
eivkc_{T}(k)dk
=
Z ∞
−∞
eivk
Z ∞
k
eα k_{e}−rT_{(e}s_{− e}k_{)}Q(s)dsdk
=
Z ∞
−∞
e−rTQ(s)
Z s
−∞
eivkeα k_{(e}s_{− e}k_{)dk}
ds
=
Z ∞
−∞
e−rTQ(s)
es
Z s
−∞
e(iv+α)kdk−
Z s
−∞
e(iv+α+1)kdk
ds
=
Z ∞
−∞
e−rTQ(s)
"
es e(iv+α)k
iv+ α
#s
−∞
−
"
e(iv+α+1)k
iv+ α + 1
#s
−∞
!
ds
=
Z ∞
−∞
e−rTQ(s) e
(1+iv+α)s
iv+ α −
e(iv+α+1)s
iv+ α + 1
!
ds
=
Z ∞
−∞
e−rTQ(s) e

(1+iv+α)s_{[(iv + α + 1) − (iv + α)]}

(iv + α)(iv + α + 1)
!
ds
= e−rT
Z _{∞}
−∞
e(1+iv+α)s
(iv + α)(iv + α + 1) Q(s)ds.
Note that (1 + iv + α) = i(v − i(1 + α)).

ψ(v) = e−rT
Z _{∞}
−∞
ei(v−i(1+α))s
(iv + α)(iv + α + 1) Q(s)ds
= e
−rT
(iv + α)(iv + α + 1)
Z ∞
−∞
ei(v−i(1+α))s Q(s)ds.

The expression in the integral will become φ (v − i(1 + α)) by the definition of a characteristic function. Then, we have:

ψ(v) = e −rT φ (v − i(1 + α )) (iv + α)(iv + α + 1) = e−rTφ (v − i(1 + α )) α2+ α − v2+ iv(2α + 1) (2.15)

### 2.4

### The Fourier Inversion

We retrieve CT(k) by the method of inverse FT defined as follows.

Definition 12. The inverse FT is given by: F−1 [ψ(v)](k) = 1 2π Z Re −ivk ψ(v)dv.

In our case, since we transformed cT(k) = eα kCT(k), we usee−αk/2π instead of1/2π to get

the original function C(k).
CT(k) =
e−αk
2π
Z
Re
−ivk_{ψ(v)dv =} e−αk
π
Z ∞
0
e−ivkψ(v)dv (2.16)
= e
−αk
π
Z ∞
0
e−ivk e
−rT_{φ (v − i(1 + α ))}
α2+ α − v2+ iv(2α + 1)dv.

The modification CT(k) to cT(k) makes it integrable along the negative k axis and at the

same time it affects the integrability along the positive k. For cT(k) to be integrable on both

ends, the following sufficient condition has to be fulfilled: ψ(0) < ∞.

We call this condition a sufficient condition for square-integrability ψ(0) = e

−rT

φ (−i(1 + α )) α2+ α . The expression e−rT

α2+α is finite because both r and T are non–negative. Therefore, we need

to check whether φ (−i(1 + α)) < ∞ for the sufficient condition of square-integrability. To simplify, we use the definition of a characteristic function [19]:

φsT(u) = E[e

iusT_{] = E[e}iuln ST_{]}
E[(eln ST_{)}iu_{] = E[S}iu

T].

By replacing u by −i(1 + α), we proceed following the same manner. φ (−i(1 + α )) < ∞

m

E[S_{T}{i(−i(1+α)}] < ∞
m

E[S_{T}1+α] < ∞.

For equation (2.16) to be square integrable the (1 + α)th moment has to exist and be finite. The value of α has an impact on both the accuracy and efficiency of our derivation. We will check the preferable value of α in the numerical implementation section.

## Chapter 3

## Option Pricing Using the FFT

### 3.1

### DFT and FFT of the BSM Formula

The implementation of the FT requires discretization of the inverse FT using approximation methods for computing integrals numerically, for example trapezoidal- and Simpson’s rule.

### 3.1.1

### Discretization using the Trapezoidal Rule

Theorem 5. [16] Suppose that we have an integrand f (x) on the interval [a, b] , and the inter-val is subdivided into N sub-interinter-vals [xj, xj+1] of width h =(b−a)/N by using the equally

spaced nodes x_{j} = a + jh, for j = 0, 1, ..., N. The composite trapezoidal rule for N
sub-intervals can be expressed as:

Z b

a

f(x)dx ≈ T ( f , h) = h

2( f (x0) + 2 f (x1) + 2 f (x2) + · · · + 2 f (xN−2) + 2 f (xN−1) + f (xN)). The proof can be found in the referenced book [16].

Truncating the upper limit of the integration for the inverse FT (2.16) to a finite number, say B, and applying the trapezoidal rule, gives the discrete Fourier transform (DFT).

CT(k) =
e−αk
π
Z _{∞}
0
e−ivkψ(v)dv
≈e
−αk
π
Z B
0
e−ivkψ(v)dv
CT(k) ≈
e−αk
π
N

### ∑

j=1 e−ivjk_{ψ}T(vj)η, (3.1) where vj= η( j − 1) and η =B/N.

Remark2. Note that η and B correspond to h and b respectively in Theorem 5. The count of summation starts in equation (3.1) from j = 1, as a result we would expect it to run to N + 1.

However, the last term is discarded because it’s contribution will become insignificant as it becomes very small. Additionally, the lower limit for the integration of C(k) is 0 thus a = 0 . That is why we have vj= η( j − 1) , while xj= a + jh in the theorem.

According to [10] however, the weights are expressed in a slightly different way.:

CT(k) ≈ e−αk π N

### ∑

j=1 e−ivjk_{ψ}T(vj) η 2(2 − δj−1), where δn= ( 1 for n = 0 0 for n 6= 0.

While the DFT requires n multiplication, i.e. O(n) for a single option price, the most effi-cient way would be to compute a number of option prices at once. The fast Fourier transform FFT algorithm, taking advantage of the complex roots of unity (see Definition 13 below), en-ables computation of a set of DFT with an n-point series that are power of two in O(n log n) multiplications that would have required O(n2) multiplications using the regular integration method. Cooley and Tukey introduced the FFT algorithm in 1965 and has gone through im-provements [10].

Definition 13. A complex nth root of unity is a complex number ω such that ω2= 1 . There
are exactly n complex nth roots of unity: e2πik/n _{for k = 0, 1, . . . , n − 1. An exponential of a}
complex number can also be written using Euler’s formula: eiu= cos(u) + i sin(u) [6].

Definition 14. [23] The FFT maps a vector x = (xj)N_{j=1} on to some vector g(xu) and the

algorithm is given by:

g(xu) = N

### ∑

j=1

e−i2πN( j−1)(u−1)x_{j} for u = 1, . . . , N.

The conversion of the DFT call price to FFT is done in a few step:

Step 1. Create a range, say [−b, b], of the logarithm of the strike price k partitioned into N with equal spacing λ :

k_{u}= −b + λ (u − 1), for u = 1, . . . , N.

Step 2. Substitute the vjand kuin the DFT equation (3.1):

C_{T}(k_{u}) ≈ e
−αku
π
N

### ∑

j=1 e−ivjku_{ψ}T(vj)η =e −αku π N

### ∑

j=1 e−ivj(−b+λ (u−1))_{ψ}T(vj)η = e −αku π N

### ∑

j=1eibvj−ivjλ (u−1)_{ψ}
T(vj)η

= e −αku π N

### ∑

j=1e−ivjλ (u−1)_{e}ibvj_{ψ}

T(vj)η = e −αku π N

### ∑

j=1e−iη( j−1)λ (u−1)eibvj_{ψ}

T(vj)η
C_{T}(ku) ≈
e−αku
π
N

### ∑

j=1e−iλ η( j−1)(u−1)eibvj_{ψ}

T(vj)η, u= 1, . . . , N. (3.2)

Step 3. Let λ η =2π_{/}_{N}_{and e}ibvj_{ψ}

T(vj)η = xjto get the FFT that corresponds to Definition 14

for the European call price CT(ku) at N number of strike points.

### 3.1.2

### Discretization using Simpson’s Rule

The fact that η = B_{N} and λ η =2π_{N} induces a trade off between η (which determines the
accur-acy of the integration) and λ (which is spacing of the logarithmic strike price and determines
the chance of getting relevant strike prices) for a fixed computational budget N [10]. To give
more weights to the more relevant points the Simpson rule as in Theorem 6 below can be
applied to the inverse FT in equation (2.16).

Theorem 6. [23] Suppose that [a, b] is subdivided into 2N sub-intervals [xj, xj+1] of equal

width h= (b − a)/(2N) by using xj= a + jh for j = 0, 1, . . . , 2N. The composite Simpson rule

for2N sub-intervals can be expressed as: Z b

a

f(x)dx ≈ S( f , h) = h

3( f (x0) + 4 f (x1) + 2 f (x2) + 4 f (x3) + . . . +2 f (x2N−2) + 4 f (x2N−1) + f (x2N)).

Applying Theorem 6, expression (3.1) becomes: C(ku) = e−αku π N

### ∑

j=1e−i2πN ( j−1)(u−1)eibvjψ_{T}(v_{j})η

3(3 + (−1)
j_{− δ}
j−1), (3.3)
where vj= η( j − 1), η =B/Nand δn=
(
1 for n = 0
0 for n 6= 0.

In this computation, Carr and Madan [5] are particularly interested at the option price of at-the-money (ATM) condition. For the OTM case, they introduce a damping sinh(αk) function as they state that the former becomes unstable as k is close to zero. This part will not be covered in this paper.

### 3.2

### Numerical Implementation

In this section we attempt to compute the call prices obtained by the FFT method. For this purpose, we utilize the MATLAB function optByMertonFFT [18] modified for the

BSM model. Related functions (i.e. functions called from optByMertonFFT) such as optSensByMertonFFT,optSensByFFT, and characteristicFcnMerton76 [18] are also modified to serve our purpose.

Comparisons are done using blsprice[17] and a numerical integration method without the transformation. The numerical integration is done by modifying MATLAB functions optByMertonNI, optSensByMertonNI, and optSensByNI [18] to BSM model.

In Tables 3.1-2, the time elapsed for computation of call prices using FFT methods (derived
using Simpson’s and Trapezoidal rule), blsprice function, and numerical integration are
presented when N = 213and when N = 216, respectively. Other parameters are fixed as: α = 3,
η = 0.1, λ = 2π /Nη for all computations presented in Tables 3.1-6, where N is the number
of grid points of the integration, α is the damping factor, η is step size of the integration, and
λ is the step size of logarithmic strike price. We consider values of the stock price at time zero
S_{0}, the risk free rate, and volatility of the stock price 100, 0.03, and 0.16, respectively.

Table 3.1: Time elapsed when N = 213. Elapsed time is 0.015856 seconds. % FFT using Simpson

Elapsed time is 0.015700 seconds. % FFT using Trapezoidal

Elapsed time is 0.016305 seconds. % blsprice

Elapsed time is 13.021780 seconds.% Numerical integration

Table 3.2: Time elapsed when N = 216. Elapsed time is 0.065797 seconds. % FFT using Simpson

Elapsed time is 0.061066 seconds. % FFT using Trapezoidal

Elapsed time is 0.083717 seconds. % blsprice

Elapsed time is 93.773157 seconds.% Numerical integration

Remark 3. The time might differ depending on the performance of the computer. This ex-periment is done with a 64-bit window 10 operating system, 8 RAM, and 2.5GHz processor personal computer.

Keeping other parameters constant, the value of N affects the computation speed signi-ficantly. Further, we notice that the FFT derived using the Trapezoidal rule followed by the FFT by Simpson’s rule are faster than the other two. The order of speeds are consistent in both cases, i.e. when N = 213 and N = 216. When N is smaller however, say N = 210, the blspriceoutperforms both FFT methods. The time elapsed varies every time we run the code. Thus, we can not induce that the times show how fast the method is. However, we can infer the that the FFT method is faster than the other two because the orders are consistent.

Tables 3.3-4 shows the comparison of the FFT call prices with call prices computed using numerical integration (with out the FT) when N = 213 and N = 216, respectively. Error1 are relative errors between the FFT call by Simpson’s rule and the numerically obtained call prices while Error2 are relative errors between the FFT call prices by Trapezoidal rule and the call prices by numerical integration with respect to the numerical method (in %). The relative errors, e.g. Error1 and Error2 are calculated as follows:

Error1 = abs(SimpFFT-CallNI)./CallNI.∗100; Error2 = abs(TrapFFT-CallNI)./CallNI.∗100;

Table 3.3: Comparison of FFT using Simpson and FFT using Trapezoidal verses call prices obtained using a numerical method without the FT at N = 213.

Strike SimpFFT TrapFFT CallNI Error1 Error2

______ _______ _______ ______ __________ __________ 96.238 7.4943 7.4943 7.4943 6.5489e-09 6.5491e-09 96.979 7.0214 7.0214 7.0214 6.6323e-09 6.6321e-09 97.725 6.5634 6.5634 6.5634 6.865e-09 6.8653e-09 98.478 6.1211 6.1211 6.1211 7.072e-09 7.0718e-09 99.236 5.6951 5.6951 5.6951 7.3921e-09 7.3924e-09 100 5.2858 5.2858 5.2858 7.8501e-09 7.8499e-09 100.77 4.8938 4.8938 4.8938 8.391e-09 8.3913e-09 101.55 4.5193 4.5193 4.5193 9.1197e-09 9.1195e-09 102.33 4.1626 4.1626 4.1626 9.9824e-09 9.9827e-09 103.12 3.8239 3.8239 3.8239 1.1012e-08 1.1011e-08 103.91 3.5033 3.5033 3.5033 1.2343e-08 1.2343e-08 104.71 3.2006 3.2006 3.2006 1.4045e-08 1.4045e-08

Moreover, Tables 3.5-6 show comparison of the FFT call prices with call prices obtained using the blsprice MATLAB built-in function when N = 213 and N = 216, respectively. Error3are relative errors between the FFT call by Simpson’s rule and the call prices by blspricewhile Error4 are relative errors between the FFT call prices by Trapezoidal rule and the call prices by blsprice function, both with respect to the prices by the blsprice function.

All relative errors shown in Tables 3.5-6 are below 1%. Additionally, it can be observed graphically that call prices by blsprice and prices by the FFT method entirely overlaps for the suggested parameter values.This implies that the FFT pricing method has a remarkable ac-curacy. We do not notice significant differences in the accuracy of Simpson’s and Trapezoidal derivation of the FFT.

The strike grids become finer as the vale of N increases. Therefore the values we have in the ranges of the strikes when N = 213 and when N = 216 are not the same.Thus we can not compare the relative errors but we notice that the method preforms well for the ATM condition, i.e. when K < 100.

Table 3.4: Comparison of FFT using Simpson and FFT using Trapezoidal verses call prices obtained using a numerical method without the FT at N = 216.

Strike SimpFFT TrapFFT CallNI Error1 Error2

______ _______ _______ ______ __________ __________ 99.522 5.5396 5.5396 5.5396 7.7898e-09 7.7884e-09 99.617 5.4883 5.4883 5.4883 7.7564e-09 7.7579e-09 99.713 5.4373 5.4373 5.4373 7.5114e-09 7.5098e-09 99.808 5.3865 5.3865 5.3865 7.6359e-09 7.6377e-09 99.904 5.336 5.336 5.336 7.5498e-09 7.548e-09 100 5.2858 5.2858 5.2858 7.84e-09 7.8418e-09 100.1 5.2359 5.2359 5.2359 7.9191e-09 7.9173e-09 100.19 5.1862 5.1862 5.1862 7.7675e-09 7.7694e-09 100.29 5.1368 5.1368 5.1368 8.0136e-09 8.0118e-09 100.38 5.0876 5.0876 5.0876 8.0285e-09 8.0304e-09

Table 3.5: Comparison of FFT call prices with the call prices obtained using blsprice at N= 213.

Strike SimpFFT TrapFFT Callbls Error3 Error4 ______ _______ _______ _______ _______ _______ 96.238 7.4943 7.4943 7.4705 0.31894 0.31894 96.979 7.0214 7.0214 6.9972 0.34495 0.34495 97.725 6.5634 6.5634 6.539 0.3728 0.3728 98.478 6.1211 6.1211 6.0966 0.40253 0.40253 99.236 5.6951 5.6951 5.6704 0.43424 0.43424 100 5.2858 5.2858 5.2612 0.468 0.468 100.77 4.8938 4.8938 4.8692 0.50387 0.50387 101.55 4.5193 4.5193 4.4949 0.54192 0.54192 102.33 4.1626 4.1626 4.1385 0.58223 0.58223 103.12 3.8239 3.8239 3.8001 0.62487 0.62487 103.91 3.5033 3.5033 3.4799 0.6699 0.6699 104.71 3.2006 3.2006 3.1779 0.71739 0.71739

Table 3.6: Comparison of FFT call prices with the call prices obtained using blsprice at N= 216.

Strike SimpFFT TrapFFT Callbls Error3 Error4 ______ _______ _______ _______ _______ _______ 99.522 5.5396 5.5396 5.515 0.44666 0.44666 99.617 5.4883 5.4883 5.4637 0.45086 0.45086 99.713 5.4373 5.4373 5.4127 0.45509 0.45509 99.808 5.3865 5.3865 5.3619 0.45936 0.45936 99.904 5.336 5.336 5.3114 0.46366 0.46366 100 5.2858 5.2858 5.2612 0.468 0.468 100.1 5.2359 5.2359 5.2112 0.47236 0.47236

Although the parameters α, and η are fixed for the computations presented in the Tables 3.1-6, we have also experimented with different values of α, and η for a fixed N. When α = 0 results with call prices equal to infinity. In this simple experiment, the FFT call prices do not show significant deviation (i.e. from what we expect call price to be [0, S0] by the conditions

we discussed in Chapter 2) for values of α between 0.02 and 7 when η = 0.01 and N = 216. When α > 7.5 however, the prices become oscillatory especially when the strike grids approach zero. This deviation can be easily noticed from Figs. 3.1-2, the curves in a red dashed lines. For α values < 0.003 the curve lies above the actual prices whereas the curve lies below the actual call price, as curves in the cyan doted line in Figs. 3.1-2, when values of α between 0.004 and 0.02. These intervals of α do not hold always, they vary significantly as we change the values of η, N or some of the input variables such as S0.

Figures 3.1-2 show plots of strike points verses call prices and logarithmic strike points verses call prices, respectively at arbitrary α values that reflect the property of the curves discussed in the previous paragraph.

Figure 3.1: Plot of strike price vs call price at different values of α.

Surface plots of the call prices against the time to maturity and strike (log strike) points are presented in Figs. 3.3-8. The parameters N = 216, η = 0.01, and λ = 2π/Nη for Figs. 3.3-6. Additionally, the input variables S0= 100, the risk free rate r = 0.03, the volatility of

the stock price σ = 0.16, and the strike price range K = [1 : 200] for all surface plots. Figures 3.3 and 3.5 show 3D plots of the call prices using FFT by Simpson’s rule against the time to maturity and strike points when α = 3 and α = 8.5, respectively.

Figure 3.2: Plot of log strike price vs call price at different values of α.

Figure 3.4: Option price surface. (X, Y, Z)=(Time to maturity, Log strike, Call) when α = 3.

Figure 3.6: Option price surface. (X, Y, Z)=(Time to maturity, Log strike, Call) when α = 8.5. Similarly, Figs. 3.4 and 3.6 show surface plots of the call prices using FFT by Simpson’s rule against the time to maturity and logarithmic strike points when α = 3 and α = 8.5, re-spectively. From Figs. 3.5-6, it is clearly visible that the surfaces do not behave well when α = 8.5, involving oscillations and jumps. When α is small, say α > 0.003 all other variables kept constant, the call price becomes way higher than the actual values at every points of K.

Changing the value of η also impacts both the accuracy and the integrity of the call prices. Keeping the other parameters constant (N = 216, α = 3 and λ = 2π/Nη), η values between 0.01 and 0.1 give lower relative errors for the comparisons presented in the tables above. Graphically however, we see significant deviations when η ≤ 0.001 and η ≥ 0.6 everything else kept constant. This is more visible when the maturity is shorter, see Fig. 3.7-8. However, if we increase the value of N as we decrease η the computation will not be affected because the impact on λ will be counterbalanced.

The oscillations and jumps are evident when the time to maturity is short. Carr and Madan [5] state that

for very short maturities, the call value approaches its non-analytic intrinsic value causing the integrand in the Fourier inversion to be highly oscillatory, and there-fore difficult to numerically integrate.

This situation is tackled by multiplying the original call by sinh(αk) instead of eα k_{, where α}

is the damping factor and k is the logarithmic strike price. Carr and Madan call this method "Fourier transform of the modified time value".

Figure 3.7: Option price surface. (X, Y, Z)=(Time to maturity, Strike, Call) when η = 0.0002.

## Chapter 4

## Conclusion and Future Work

### 4.1

### Conclusion

In this project, we have elaborated the derivation of the European call option pricing using the FFT method introduced by Carr and Medan [5] by connecting various mathematical and financial aspects related to the method. We have derived the CF of logarithmic terminal stock price by implementing Cauchy’s contour integral. Additionally, we have revised the BSM model and derived the BSM PDE using delta hedging while the BSM formula is derived by applying Girsanov’s theorem.

For the implementation, we have modified MATLAB functions that are made for the Mer-ton’s option pricing model when underlying returns are discontinuous [22] to apply to the BSM model. The FFT call prices by Simpson’s and the Trapezoidal rule are compared and from the relative errors we conclude that the differences are insignificant. Comparisons of call prices using the FFT method with call prices obtained using blsprice and a numerical method without the implementation of FT show small relative errors, all below 1%. This veri-fies the accuracy of the FFT method that we have implemented. We can also conclude that the FFT method is significantly faster than the numerical method without a transformation, as can be seen in Tables 3.1-2.

Further, we attempt to show how the the call prices by the FFT method behave graphically as we change the different parameter values. As we change one parameter, e.g. α or η the curves/surfaces behave quite differently. To our scope, we could not implement precision algorithms that would enable us to choose efficient values of the parameters, which we believe is essential for implementation of the FFT pricing method.

Carr and Madan have a different algorithm for short maturities. Thus one is supposed to implement different algorithms for different scenarios. This can be marked as a draw back of the method.

### 4.2

### Future Work

As we have stated in the conclusion, the choice of efficient parameters values is crucial part of the pricing method using the FFT. We urge further studies on the decision of exact parameter

values that conveniently give efficient and effective computational results. The research paper [14] would be worth investigating in this regard. Lee introduces bounds to the numerical pricing error of DFT and FFT that "enable algorithms to select efficient quadrature parameters and to price with guaranteed numerical accuracy" as he claims.

This thesis, apart from reviewing different pricing models, is restricted to the BSM model which is criticised for it’s limiting assumptions. We believe that other models that involve discontinuity, stochastic volatility and interest rate, e.g. Merton’s option pricing model when underlying returns are discontinuous [22], Heston [9], Bate [3], and so on would also be interesting to study. Considering the paper by Lord ang Kahl [15] will be helpful for studies on jump–diffusion models. Unlike other practitioners that use branch switching of the complex logarithm to eliminate the discontinuities that might arise in the characteristics function, Lord and Kahl have proved an algorithm using the principal branch [see 15, and references therein].

## Bibliography

[1] Edward B. Saff and Arthur D Snider. Fundamentals of complex analysis with applica-tions to engineering, science and mathematics, 2002.

[2] Gurdip Bakshi, Charles Cao, and Zhiwu Chen. Empirical performance of alternative option pricing models. The Journal of finance, 52(5):2003–2049, 1997.

[3] David S Bates. Jumps and stochastic volatility: Exchange rate processes implicit in deutsche mark options. The Review of Financial Studies, 9(1):69–107, 1996.

[4] Fischer Black and Myron Scholes. The pricing of options and corporate liabilities. Journal of political economy, 81(3):637–654, 1973.

[5] Peter Carr and Dilip Madan. Option valuation using the fast fourier transform. Journal of computational finance, 2(4):61–73, 1999.

[6] Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, and Clifford Stein. Introduc-tion to algorithms. MIT press, 2009.

[7] John C Cox and Stephen A Ross. The valuation of options for alternative stochastic processes. Journal of financial economics, 3(1-2):145–166, 1976.

[8] Leon M Hall. Special functions. Available in: http://web. mst. edu/˜ lmhall/SPFNS/spfns. pdf, 1995.

[9] Steven L Heston. A closed-form solution for options with stochastic volatility with applications to bond and currency options. The review of financial studies, 6(2):327– 343, 1993.

[10] Ali Hirsa and Salih N Neftci. An introduction to the mathematics of financial derivatives. Academic press, 2013.

[11] Robert V Hogg, Joseph McKean, and Allen T Craig. Introduction to mathematical stat-istics. Pearson Education, 2005.

[12] John Hull et al. Options, futures and other derivatives/John C. Hull. Upper Saddle River, NJ: Prentice Hall„ 2009.

[13] Masaaki Kijima. Stochastic processes with applications to finance. Chapman and Hall/CRC, 2016.

[14] Roger W Lee et al. Option pricing by transform methods: extensions, unification and error control. Journal of Computational Finance, 7(3):51–86, 2004.

[15] Roger Lord and Christian Kahl. Complex logarithms in heston-like models. Math-ematical Finance: An International Journal of Mathematics, Statistics and Financial Economics, 20(4):671–694, 2010.

[16] John H Mathews, Kurtis D Fink, et al. Numerical methods using MATLAB, volume 4. Pearson prentice hall Upper Saddle River, NJ, 2004.

[17] Inc. MathWorks. MATLAB 2016. The MathWorks, Inc. Natick, MA, USA, 1995-2016. [18] Inc. MathWorks. Financial Instruments Toolbox. The MathWorks, Inc. Natick, MA,

USA, 2020. pg. 3-28.

[19] Kazuhisa Matsuda. Introduction to option pricing with fourier transform: Option pricing with exponential lévy models. Department of Economics The Graduate Center, The City University of New York, pages 1–241, 2004.

[20] Robert C Merton. Continuous-time speculative processes. SIAM Review, 15(1):1–42, 1973.

[21] Robert C Merton. Theory of rational option pricing. The Bell Journal of economics and management science, pages 141–183, 1973.

[22] Robert C Merton. Option prices when underlying stock returns are discontinuous. Journal of Financial Economics, 3:125–144, 01 1976.

[23] Martin Schmelzle. Option pricing formulae using fourier transform: Theory and applic-ation. Preprint, http://pfadintegral. com, 2010.

[24] Wim Schoutens, Erwin Simons, and Jurgen Tistaert. A perfect calibration! now what? The best of Wilmott, page 281, 2003.

## Appendix A

## MATLAB Code

### A.1

### MATLAB Script of the Implementation

%Use optByBSM_FFT to compute option prices, and plot an option price surface.

S_0 = 100;%Stoke price

Rate = 0.03;%Risk free rate

OptSpec = ’call’;

Sigma = 0.16; %votatility of the Stoke price

N = 2^16;%number of strike points

alpha=3; %damping factor

eta=0.1; %step size of the Fourier transform

lambda=2*pi/N/eta; %step size of the logarithmic strike points %FRFTParam = (lambda .* eta) ./2 ./pi;

%(abs(FRFTParam - 1/N) < 1e-14) %check that this is true to use fft.

Settle = datenum(’31-Jul-2020’); Maturity = datemnth(Settle, 6);

K = []; % K is not specified (will use the entire FFT K grid)

%% Compute option prices for the entire FFT K grid using Simpsons rule

tic;

[Call, Kout] = optByBSM_FFT(Rate, S_0, Settle, Maturity, OptSpec, K,Sigma, ...

’eta’, eta, ’lambda’, lambda, ’alpha’,alpha,’N’,N); toc

%% Compute option prices for the entire FFT K grid using TRapezoidal rule

tic;

[Callt, Koutt] = optByBSM_FFT(Rate, S_0, Settle, Maturity, OptSpec, K, Sigma, ...

’eta’, eta, ’lambda’, lambda, ’alpha’,alpha,...

’Quadrature’,’trapezoidal’,’N’,N); toc

tic;

Call_bls=blsprice(S_0,Kout, Rate, 0.5, Sigma); toc

%% % tic;

% CallNI = optByBsmNI(Rate, S_0, Settle, Maturity, OptSpec, Kout,Sigma); % toc

%% Compare with numerical integration method (with "BSM_NI" framework)

Range = (32764:32773); % when N=2^13 Range=(4092:4103)

Strike = Kout(Range); % This range is chosen to get an interval around ATM point.

SimpFFT = Call(Range); TrapFFT=Callt(Range);

CallNI = optByBsmNI(Rate, S_0, Settle, Maturity, OptSpec, Strike,Sigma); Error1 = abs(SimpFFT-CallNI)./CallNI*100;

Error2 = abs(TrapFFT-CallNI)./CallNI*100;

table(Strike, SimpFFT, TrapFFT,CallNI, Error1, Error2)

%% Compare with the blsprice built in function

Callbls=Call_bls(Range);

Error3 = abs(SimpFFT-Callbls)./Callbls*100; Error4 = abs(TrapFFT-Callbls)./Callbls*100;

table(Strike, SimpFFT, TrapFFT,Callbls, Error3, Error4)

%% Plot an Option Price Surface

Settle = datenum(’31-Jul-2020’);

Maturity = datemnth(Settle, 12*[1/12 0.25 (0.5:0.5:3)]’); Times = yearfrac(Settle, Maturity);

K = (1:200)’; eta=0.01;

lambda=2*pi/N/eta;

Call = optByBSM_FFT(Rate, S_0, Settle, Maturity, OptSpec, K, ... Sigma, ’N’, N, ...

’eta’, eta, ’lambda’, lambda, ’alpha’,alpha,’N’,N, ’ExpandOutput’, true);

[X,Y] = meshgrid(Times,K); figure;

surf(X,Y,Call);

title(’Call Price Surface Plot when \alpha=3’); xlabel(’Years to maturity’);

ylabel(’Strike price’); zlabel(’Call price’); view(-115,35);

xlim([0 Times(end)]); zlim([0 100]);

%% 3D fig maturity vs log K price vs Call price

[X,Y] = meshgrid(Times,log(K)); figure;

surf(X,Y,Call);

title(’Call Price Surface Plot when \alpha=3’); xlabel(’Years to maturity’);

ylabel(’Log strike price’); zlabel(’Call price’);

view(-115,35);

xlim([0 Times(end)]); zlim([0 100]);

%% Fig K price vs Call price at different alpha values

Call_1 = optByBSM_FFT(Rate, S_0, Settle, Maturity, OptSpec, K, ... Sigma, ’N’, N, ...

’eta’, eta, ’lambda’, lambda, ’alpha’,0.003,’N’,N, ’ExpandOutput’, true);

Call_2 = optByBSM_FFT(Rate, S_0, Settle, Maturity, OptSpec, K, ... Sigma, ’N’, N, ...

’eta’, eta, ’lambda’, lambda, ’alpha’,0.005,’N’,N, ’ExpandOutput’, true);

Call_3 = optByBSM_FFT(Rate, S_0, Settle, Maturity, OptSpec, K, ... Sigma, ’N’, N, ...

’eta’, eta, ’lambda’, lambda, ’alpha’,3,’N’,N, ’ExpandOutput’, true);

Call_4 = optByBSM_FFT(Rate, S_0, Settle, Maturity, OptSpec, K, ... Sigma, ’N’, N, ...

’eta’, eta, ’lambda’, lambda, ’alpha’,8.5,’N’,N, ’ExpandOutput’, true) ; figure; plot(K,Call_1(:,3),’b-.’,’linewidth’,0.5) hold on plot(K,Call_2(:,3),’c:’,’linewidth’,1.5) hold on

plot(K,Call_3(:,3),’y’,’linewidth’,1.5) hold on

plot(K,Call_4(:,3),’r--’,’linewidth’,0.5) hold off

ylim([-10 110]);

legend(’\alpha=0.003’, ’\alpha=0.005’, ’\alpha=3’, ’\alpha=8.5’) title(’Strike price vs Call price’);

xlabel(’Strike price’); ylabel(’Call price’);

%% Fig log K price vs Call price at different alpha values

figure;

plot(log(K),Call_1(:,3),’b-.’,’linewidth’,0.5) hold on

plot(log(K),Call_2(:,3),’c:’,’linewidth’,1.5) hold on

plot(log(K),Call_3(:,3),’y’,’linewidth’,1.5) hold on

plot(log(K),Call_4(:,3),’r--’,’linewidth’,0.5) hold off

ylim([-10 110]);

legend(’\alpha=0.003’, ’\alpha=0.005’, ’\alpha=3’, ’\alpha=8.5’) title(’Log strike price vs Call price’);

xlabel(’Log strike price’); ylabel(’Call price’);

### A.2

### MATLAB Function Files used in the Implementation

### A.2.1

**Function "optByBSM_FFT"**

function [Price, KOut] = optByBSM_FFT(Rate, S_0, Settle, Maturity, ... OptSpec, K, Sigma, varargin)

%optBYBSM_FFT Option price by BSM model, using FFT and FRFT

% Vanilla European option price by BSM model, using Carr-Madan FFT

% and Chourdakis FRFT methods.

%

% Inputs:

% Rate - Continuously compounded risk-free interest rate. (

Scalar)

% S_0 - Current underlying asset price. (NINST x 1, or NColumns x 1)

% Settle - Settlement date in date numbers, date character vectors

,

% or date strings. (NINST x 1, or NColumns x 1)

% Maturity - Option maturity date in date numbers, date character

% vectors, or date strings. (NINST x 1, or NColumns x 1)

% OptSpec - A vector of strings or a cell array of character

vectors

% indicating ’call’ or ’put’. (NINST x 1, or NColumns x

1)

% K - Strikes for the options.

% (NINST x 1, NRows x 1, NRows x NColumns, or [])

% If this input is an empty array ([]), option prices

will

% be computed on the entire FFT (or FRFT) K grid,

% which is determined as exp(log-K grid). Each column

% of the log-strike grid (k_u) has ’N’ points with

% ’lambda’ spacing that are roughly centered around

% each element of log(S_0), where ’k’ is the log strike,

’u’

% is a parametre at which we evaluate the characteristic

% function, and ’lambda’ is the log stike step size.

%

% _{* Note: See name-value pair ’ExpandOuput’ for more information}

about

% the function behavior and proper dimensions for S_0,

%

% Sigma - Volatility of the underlying asset. (Scalar)

% %

% Optional Inputs:

% ’Basis’ - Day-count basis. (Scalar)

% 0 = actual/actual % 1 = 30/360 (SIA) % 2 = actual/360 % 3 = actual/365 % 4 = 30/360 (BMA) % 5 = 30/360 (ISDA) % 6 = 30/360 (European) % 7 = actual/365 (Japanese) % 8 = actual/actual(ICMA) % 9 = actual/360 (ICMA) % 10 = actual/365(ICMA) % 11 = 30/360E (ICMA) % 12 = actual/actual(ISDA) % 13 = BUS/252 % Default: 0 %

% ’DividendYield’ - Continuously compounded underlying asset yield. (

Scalar)

% Default: 0

%

% ’N’ - Number of grid points in the characteristic function

% variable and in each column of the log-K grid. (Scalar)

% Default: 4096 (i.e., 2^12-point FFT or FRFT)

%

% ’eta’ - Characteristic function variable grid spacing. (Scalar)

% Default: 0.01

%

% ’lambda’ - Log-K grid spacing. (Scalar)

% Default: 2*pi/N/eta

%

% _{* Note: If (lambda*eta) is 2*pi/N, FFT is used.}

% Otherwise, FRFT is used.

%

% ’alpha’ - Damping factor for Carr-Madan formulation. (Scalar)

% Default: 1.5

%

% ’Quadrature’ - Single string or character vector indicating the type

% of quadrature, including: "trapezoidal", "simpson"

% Default: "simpson"

%

% ’ExpandOutput’ - Logical indicating whether to expand the outputs. (

Scalar)

% false: (outputs are NINST x 1 vectors)

% If false, the outputs are NINST x 1 vectors. Also,

the