Pricing Put Options with Multilevel Monte Carlo Simulation

School of Education, Culture and Communication

Division of Mathematics and Physics

BACHELOR’S DEGREE PROJECT IN MATHEMATICS

Pricing European Put Options with Multilevel Monte Carlo Simulation

by

Jonathan Schöön

MAA322 — Examensarbete i matematik för kandidatexamen

DIVISION OF MATHEMATICS AND PHYSICS

MÄLARDALEN UNIVERSITY SE-721 23 VÄSTERÅS, SWEDEN

School of Education, Culture and Communication

Division of Mathematics and Physics

MAA322 — Bachelor’s Degree Project in Mathematics

Date of presentation:

3rd June 2021

Project name:

Pricing European Put Options with Multilevel Monte Carlo Simulation

Author: Jonathan Schöön Version: 21st June 2021 Supervisor: Doghonay Arjmand Reviewer: Daniel Andrén Examiner: Fredrik Jansson Comprising: 15 ECTS credits

Abstract

Monte Carlo path simulations are common in mathematical and computational finance as a way of estimating the expected values of a quantity such as a European put option, which is functional to the solution of a stochastic differential equation (SDE). The computational complexity of the standard Monte Carlo (MC) method grows quite large quickly, so in this thesis we focus on the Multilevel Monte Carlo (MLMC) method by Giles, which uses multigrid ideas to reduce the computational complexity. We use a Euler-Maruyama time discretisation for the approximation of the SDE and investigate how the convergence rate of the MLMC method improves the computational times and cost in comparison with the standard MC method. We perform a numerical analysis on the computational times and costs in order to achieve the desired accuracy and present our findings on the performance of the MLMC method on a European put option compared to the standard MC method.

Acknowledgements

I would like to thank my supervisor Doghonay Arjmand for his support throughout the process of writing this thesis. He was always available for meetings when I needed it and his knowledge on the subject helped me whenever I needed guidance.

Contents

Acknowledgements 2 1 Introduction 6 2 Preliminaries 10 2.1 Brownian Motion . . . 10 2.2 Itô’s Lemma . . . 10

2.3 Geometric Brownian Motion . . . 12

2.4 Time discretization for deterministic problems . . . 13

2.4.1 Local and Global Truncation Error . . . 13

2.5 Approximation of SDE . . . 14

3 Standard Monte Carlo 16 3.1 Error estimates . . . 17

4 Multilevel Monte Carlo 19 4.1 Giles Complexity Theorem . . . 22

5 Numerical Results 24 5.1 Verification of strong and weak convergence . . . 24

5.2 Computational Times . . . 24

5.3 Pricing the European put option with MC and MLMC . . . 25

6 Conclusion 31 A MATLAB Codes 34 A.1 Convergence check . . . 34

A.2 MLMC function . . . 35

A.3 Computational times . . . 35

List of Figures

2.1 Brownian Paths . . . 11

3.1 Discretisation Error RMSE . . . 18

5.1 Strong and weak convergence . . . 25

5.2 Computational times with different values of M . . . 26

5.3 Computational times with different values of L . . . 26

5.4 MC behaviour with 𝑀 = 218, 𝐿 = 5 . . . 27

5.5 MLMC behaviour with 𝑀 = 218, 𝐿 = 5 . . . 27

5.6 MC behaviour with 𝑀 = 214, 𝐿 = 10 . . . 28

List of Tables

5.1 Expected values of European put option with MLMC for different levels 𝑀 and 𝐿 . . . 29 5.2 MSE’s of European put option with MLMC for different levels 𝑀 and 𝐿 . . . 29 5.3 Computational Times (in seconds) of European put option with MLMC for

Chapter 1

Introduction

Contracts similar to options have been known to man since ancient times with the earliest known occurence being the Greek mathematician and philosopher Thales of Miletus. Thales speculated in the harvest of olives and decided to reserve the exclusive rights to use the olive presses in his region and then made a fortune when the crop was successful due to olive presses being in high demand.[10] Other examples throughout early history include the D¯ojima rice exchange in Edo period Japan[21], and the trading of "opsies" on the Amsterdam stock ex-change in the 17th century.[9]. Options as we now know them have been around for decades but they were first established and standarized on the Chicago Board Options exchange in April of 1973 on the 125th birthday of the Chicago Board of Trade[19] and the popularity of options continued to grow from then onwards. In the early 1970s the groundbreaking work of Fischer Black, Myron Scholes, and Robert Merton developed the model known as the Black-Scholes-Merton model, or simply the Black-Scholes model (BS), which changed the way we price European stock options by figuring out how to determine a relationship between the market’s required return on the option to the required return on the stock which was very tricky due to the nature of the relationship depending on both the stock price and the time of the option[13]. The BS-model is dependent on geometric Brownian motion (GBM) which is a continous-time stochastic process that follows a Brownian motion with drift. This is also known as the Wiener process. GBM must satisfy a stochastic differential equation (SDE) and is logarithmic in nature. In this thesis we will focus on European stock options[13], more specifically on European put options and if there is a way to implement a working multilevel Monte Carlo simulation on them. Options are a class of financial contracts which gives the owner the right to, within a specified timeframe, buy or sell a good at a predetermined price. These options can be divided into two major classes, call options which gives the owner the right to buy an asset, and put

options which gives the owner the right to sell an asset. The predetermined price is called the

strike price 𝐾 and the date of maturity 𝑇 is the specified end date of the option. It is important to notice that owners have the right to exercise their option, but are not obligated to do so if they do not want to.

It is common practice that options contracts consists of 100 shares of the underlying asset, and the cost of the contract, which is called the premium, is charged per share. Therefore, if the premium is 10 cents per share the total cost of the contract would be $0.10 ∗ 100 = $10 .

The premium is based on the strike and the time of maturity of the option and is therefore different for each contract depending on how it is drawn. [17]

There are a number of different types of options but the some of the most common ones are European options, American options, and Asian options. European and Asian options can only be exercised at the time of maturity, whereas American options may be exercised at any point between the purchase date and the time of maturity. The payoff, where 𝑆(𝑇 ) is the price of the asset at the time of maturity 𝑇 , for a European call option is

𝑆_𝐶(𝑇 ) = 𝑚𝑎𝑥 (𝑆(𝑇 ) − 𝐾, 0) = (

𝑆(𝑇 ) − 𝐾, 𝑆(𝑇 ) ≥ 𝐾 0, 𝑆(𝑇 ) < 𝐾 .

(1.1) Similarly for a European put option the payoff is

𝑆𝑃(𝑇 ) = 𝑚𝑎𝑥 (𝐾 − 𝑆(𝑇 ), 0) =

(

𝐾− 𝑆(𝑇 ), 𝐾 ≥ 𝑆(𝑇 ) 0, 𝐾 < 𝑆(𝑇 ).

(1.2) As it can be argued that option prices can be expressed as the expected value, due to the case that the solution to the SDE used in the pricing model is also expressed as an expected value based on the Feynman-Kaˆc formula, it can be shown that the payoff of a European call option can be modified as such [5]

𝑓𝐶𝐸(𝑡, 𝑠) ≡ E

h

𝑒−𝑟 (𝑇 −𝑡)𝑚 𝑎𝑥(𝑆(𝑇 ) − 𝐾, 0|𝑆(𝑡) = 𝑠 i

(1.3) and similarly for the European put option

𝑓_𝑃 𝐸(𝑡, 𝑠) ≡ E h 𝑒−𝑟 (𝑇 −𝑡)𝑚 𝑎𝑥(𝐾 − 𝑆(𝑇 ), 0|𝑆(𝑡) = 𝑠 i (1.4) Where 𝑟 is the riskless interest rate.

The payoff formula for American options are the same as European options with the only difference being that American options can be exercised at any time 0 ≤ 𝑡 ≤ 𝑇 .

𝑓_𝐶 𝐴𝑚(𝑡, 𝑠) ≡ 𝑚𝑎𝑥 E h 𝑒−𝑟 (𝑇 −𝑡)𝑚 𝑎𝑥(𝑆(𝑇 ) − 𝐾, 0|𝑆(𝑡) = 𝑠 i (1.5) 𝑓_𝑃 𝐴𝑚(𝑡, 𝑠) ≡ 𝑚𝑎𝑥 E h 𝑒−𝑟 (𝑇 −𝑡)𝑚 𝑎𝑥(𝐾 − 𝑆(𝑇 ), 0|𝑆(𝑡) = 𝑠 i (1.6) Asian options differ in that the payoff is based on the average price of the underlying asset during the time from purchase to maturity. The average price for the option can be determined by either the arithmetic average or the geometric average and therefore one can purchase either a arithmetic Asian option or a geometric Asian option. For simplicity’s sake we choose to only present the payoff formula for the arithmetic Asian option.

𝑓_𝐶 𝐴𝑠(𝑡, 𝑠) ≡ E h 𝑒−𝑟 (𝑇 −𝑡)𝑚 𝑎𝑥( 𝐴(0, 𝑇 ) − 𝐾, 0|𝑆(𝑡) = 𝑠 i (1.7) 𝑓_𝑃 𝐴𝑠(𝑡, 𝑠) ≡ E h 𝑒−𝑟 (𝑇 −𝑡)𝑚 𝑎𝑥(𝐾 − 𝐴(0, 𝑇 ), 0|𝑆(𝑡) = 𝑠 i (1.8) Where 𝐴(0, 𝑇 ) = _𝑇1 ∫𝑇

0 𝑆(𝑡)𝑑𝑡. Since we are focusing on the European put option in this thesis

the most important payoff function to remember is (1.4)

In financial mathematics Monte Carlo (MC) simulations are computational algorithms used to simulate expected payoffs in financial derivatives by repeated random sampling and statistical analysis. The name Monte Carlo comes from the casinos of Monte Carlo and was named as such as a code name by nuclear physicists working on the Manhattan Project during World War II who needed a way to simulate the behavior of neutron diffusion in fissionable material[20]. MC simulations are advantageous in that it can accommodate for any stochastic process of a single market variable 𝑆 that provides a payoff at time 𝑇 both when the payoff depends on the path followed by the underlying variable 𝑆 and when it only depends on the final value of 𝑆, but a disadvantage to this method is that it can be very time consuming and complex computationally[13]. The MC method is inherently slow due to its use of numerical approximation of stochastic processes whenever the true solution to the underlying SDE is not apparent. In MC the estimated variance converges as 1/𝑁, where 𝑁 is the number of samples, and thus the inherent estimation error converges as 1/√𝑛[24]. With the addition of approximation another error, the discretization error, is introduced and therefore additional computational complexity is required as one seeks to discount for this error by finding a high level of accuracy of the estimation. The higher the level of accuracy sought, the quicker the increase in computational cost and time grows. Multiple researchers have developed different methods to improve the cost of the MC method by so called variance reduction techniques, see [12] [15]. The problem with these techniques is that while improving computational complexity, they do not necessarily improve computational cost. A method that does im-prove computational complexity however is the so called quasi MC method. The quasi MC method uses quasi-random sequences such as the Halton and Hammersley, Faure, and Sobol sequences, instead of the pseudorandom sequences of the MC method, in order to reduce the estimation error by selecting samples instead of using random sampling.[12] The drawback of this method is that it requires a lot more knowledge to be implemented efficiently. We will try to develop a method to decrease this computational time by implementing a Multilevel Monte Carlo method. Multilevel Monte Carlo (MLMC) methods are similar to the Monte Carlo in that it uses random sampling, but in order to reduce the computational complexity it takes these random samples on different levels of accuracy. The MLMC method is attributed to the work of professor Michael B. Giles and we will use his work as a basis for this thesis[11].

In the next chapter we will introduce the necessary theory needed for later sections by giving introducing the stochastic differential equation used for this thesis, as well as conducting a numerical discretization on the SDE. Chapter 3 will focus on the Monte Carlo method and chapter 4 on Multilevel Monte Carlo simulations based on the European put option. In chapter

5 we will present our findings and conduct a numerical analysis on our results, and finally in chapter 6 we present our conclusion on the subject.

Chapter 2

Preliminaries

2.1

Brownian Motion

A common example of a stochastic process is the Brownian motion, which is also known as the Wiener Process.

Definition 1. (Definition 2.8 from [8]) The one-dimensional Wiener process 𝑊 : [0, ∞) ×𝜔 →

R, has the following properties:

1. with probability 1, the mapping 𝑡 → 𝑊 (𝑡) is continuous and 𝑊 (0) = 0; 2. if 0 = 𝑡0 < 𝑡₁ < ... < 𝑡_𝑁 = 𝑇, then the increments

𝑊(𝑡_𝑁) − 𝑊 (𝑡_𝑁−1), ..., 𝑊 (𝑡₁) − 𝑊 (𝑡₀) are independent; and

3. for all 𝑡 > 𝑠 the increment 𝑊 (𝑡) − 𝑊 (𝑠) has the normal distribution, with E[𝑊 (𝑡) − 𝑊(𝑠)] = 0 and E[(𝑊 (𝑡) − 𝑊 (𝑠))2] = 𝑡 − 𝑠

If we intend to generate a Wiener process we need to determine 𝑊. We do so by generating a finite amount of time steps {𝑡𝑛 : 𝑛 = 0, . . . , 𝑁 } in the form 0 = 𝑡0 < 𝑡1 < . . . < 𝑡𝑁 = 𝑇.

With definition 1 we then generate 𝑊 (𝑡𝑛) as the sum of the independent normally distributed

random variables denoted by 4𝑊𝑛 = 𝑊 (𝑡𝑛+1) −𝑊 (𝑡𝑛) By the definition above we see that when

the time 𝑡 is fixed, the Brownian motion 𝑊 (𝑡) is then a normally distributed random variable as well. It is important to note that it computationally unrealistic to generate 𝑊 for all 𝑡 ∈ R as the process would require infinite amounts of work done in the computation [8]. See Fig(2.1) for an example of 50 simulated Brownian paths with 100 timesteps.

2.2

Itô’s Lemma

Itô’s lemma is a commonly used identity in finance and is of important use in the derivation of the BS model. It is used to find the differential of time-dependent functions of stochastic

processes and its characteristics are similar to the chain rule used in calculus. Itô’s lemma can be derived by expanding a Taylor series and applying the rules of stochastic calculus.

Suppose that the value of a variable 𝑥 follows the Itô process

𝑑𝑥 = 𝑎 (𝑥, 𝑡)𝑑𝑡 + 𝑏 (𝑥, 𝑡)𝑑𝑧 (2.1) where 𝑑𝑧 is a Wiener process and 𝑎 and 𝑏 are functions of 𝑥 and 𝑡. The variable 𝑥 has a drift rate of 𝑎 and a variance rate of 𝑏2. Itô’s lemma shows that a function 𝐺 of 𝑥 and 𝑡 follows the process 𝑑 𝐺 =  𝜕 𝐺 𝜕 𝑥 𝑎+ 𝜕 𝐺 𝜕 𝑡 + 1 2 𝜕2𝐺 𝜕 𝑥2 𝑏2  𝑑 𝑡+ 𝜕 𝐺 𝜕 𝑥 𝑏 𝑑 𝑧 (2.2) where the 𝑑𝑧 is the same Wiener process as in equation (2.1): Thus, 𝐺 also follows an Itô process, with a drift rate of

𝜕 𝐺 𝜕 𝑥 𝑎+ 𝜕 𝐺 𝜕 𝑡 + 1 2 𝜕2𝐺 𝜕 𝑥2 𝑏2 and a variance rate of

 𝜕 𝐺

𝜕 𝑥 2

𝑏2 For proof see [22].

2.3

Geometric Brownian Motion

In order to model stock prices according to the BS-model we use the continuous-time stochastic process known as Geometric Brownian motion (GBM). In order for a stochastic process 𝑆(𝑡) to follow GBM it must satisfy the following SDE

𝑑𝑆(𝑡) = 𝜇𝑆(𝑡)𝑑𝑡 + 𝜎𝑆(𝑡)𝑑𝑊 (𝑡) (2.3) where the constant 𝜇 is the stochastic drift, the constant 𝜎 is the volatility, and 𝑊 (𝑡) is a Wiener process.

Equation (2.3) when applied with Itô’s lemma implies that 𝑑 𝑙 𝑜𝑔 𝑆(𝑡) = (𝜇 −

1 2𝜎

2_{)𝑑𝑡 + 𝜎𝑑𝑊 (𝑡) (2.4)}

From (2.4) we see that if 𝑆 ∼ 𝐺 𝐵𝑀 (𝜇, 𝜎2) and if S has initial value 𝑆(0) then 𝑆(𝑡) = 𝑆(0)𝑒( 𝜇−

𝜎 2

2 )𝑡+𝜎𝑊 (𝑡) (2.5)

Or a bit more general, if 𝑢 < 𝑡 then

𝑆(𝑡) = 𝑆(𝑢)𝑒( [𝜇−

𝜎 2

2 ]) (𝑡−𝑢)+𝜎 (𝑊 (𝑡)−𝑊 (𝑢)) (2.6)

2.4

Time discretization for deterministic problems

The Euler method is a method of solving ordinary differential equations (ODE) with a given initial value, and is named after the mathematician Leonhard Euler whom developed the theory in the 18th century[4]. The idea of the Euler method is to approximate the curve of an ODE by computing the slope of the tangent line at any point of the curve using small, equidistant, time-steps ℎ on an interval (0, 𝑇 ). Consider the first order form ODE

𝑑𝑦 𝑑 𝑡

= 𝑓 (𝑡, 𝑦 (𝑡)), 𝑦(𝑡₀) = 𝑦₀. (2.7) Through Euler discretization we can write this as

𝑦(𝑡 + ℎ) − (𝑦 (𝑡) ℎ

≈ 𝑓 (𝑡, 𝑦 (𝑡)). (2.8) or

𝑦(𝑡 + ℎ) ≈ 𝑦 (𝑡) + ℎ 𝑓 (𝑡, 𝑦 (𝑡)) (2.9) Equation (2.10) is called the Forward Euler’s method and the value of 𝑦(𝑡) is the approximation to the solution of the ODE at time 𝑡 : 𝑦(𝑡) ≈ 𝑦(𝑡 + ℎ) and is explicit as the solution 𝑦(𝑡 + ℎ) is an explicit function of 𝑦(𝑡) for each timestep ℎ. Similarly another way of approximating (2.8) is 𝑦(𝑡 + ℎ) − 𝑦 (𝑡) ℎ ≈ 𝑓 (𝑡 + ℎ, 𝑦 (𝑡 + ℎ)) (2.10) which is equal to 𝑦(𝑡 + ℎ) ≈ 𝑦 (𝑡) + ℎ 𝑓 (𝑡 + ℎ, 𝑦 (𝑡 + ℎ)). (2.11) This is known as the Backward Euler’s method and it is implicit as the approximation 𝑦(𝑡 + ℎ) appears on both sides of the equation. When the right hand side is nonlinear one needs to use some other method such as fixed point iteration or the Newton-Raphson method to solve the algebraic equation of 𝑦(𝑡 + ℎ).

2.4.1

Local and Global Truncation Error

The local truncation error (LTE) is the error made in a single step of the Euler method and it is the difference between the numerical solution and the exact solution that occurs after one step. Given the numerical solution

𝑦₁= 𝑦₀+ ℎ 𝑓 (𝑡₀, 𝑦₀) (2.12) and the exact solution which uses a Taylor expansion assuming 𝑦(𝑡) being a smooth function

𝑦(𝑡₀+ ℎ) = 𝑦(𝑡₀) + ℎ𝑦0(𝑡₀) + ℎ2

2 𝑦

00_(𝑡

0) + O (ℎ3). (2.13)

the local truncation error (LTE) can be written as the difference between these two equations 𝐿𝑇 𝐸 = 𝑦 (𝑡₀+ ℎ) − 𝑦₁ =

1 2ℎ

2

Equation(2.14) shows that when the timestep ℎ is small the LTE is approximately proportional to ℎ2. Finding the LTE of the Backward’s Euler method follows the same principal and leads to the result that the LTE is approximately proportional to O (ℎ2).

The global truncation error (GTE) is the accumulation of LTE commited in each step from the intial time to a fixed time 𝑡 The number of steps is (𝑡 − 𝑡0)/ℎ which is proportional to 1/ℎ.

As the LTE is proportional to ℎ2we can deduce that the GTE is expected to be proportional to ℎ for the Forward Euler method, and similarly the GTE is proportional to O (ℎ) for the Backward Euler method.

2.5

Approximation of SDE

The Forward Euler method can be extended to approximate the solution of an SDE. This extended method is called the Euler-Maruyama method and is named as such after Japanese mathematician Gisiro Maruyama.

If we first consider the SDE

𝑑 𝑋(𝑡) = 𝑎( 𝑋 (𝑡))𝑑𝑡 + 𝑏( 𝑋 (𝑡))𝑑𝑊 (𝑡) (2.15) where 𝑋 is a random variable and 𝑊 is a Wiener process, the Euler-Maruyama time

discretisation of this equation can then be written as ˆ

𝑋(𝑡_𝑖+1) = ˆ𝑋(𝑡_𝑖) + 𝑎 ( ˆ𝑋(𝑡_𝑖)) [𝑡_𝑖+1− 𝑡_𝑖] + 𝑏 ( ˆ𝑋(𝑡_𝑖)) √

𝑡_𝑖+1− 𝑡_𝑖𝑍_𝑖+1 (2.16) For this approximation the time grid we consider is 0 = 𝑡0< 𝑡₁< . . . < 𝑡_𝑚 for 𝑖 = 0, . . . , 𝑚 − 1. The value of the approximation of 𝑋 at time 0 is ˆ𝑋(0) = 𝑋 (0), and 𝑍₁, 𝑍₂, . . . , 𝑍_𝑛 are independent 𝑚-dimensional standard normal random vectors[12]. See also Theorem 1 in [8]. Since we are interested in generating option prices using GBM we can apply the Euler-Maruyama method to equation (2.3) and get the following equation for the approximation where we substitute 𝑆(𝑡) for 𝑋 (𝑡), 𝜇 for 𝑎, and 𝜎 for 𝑏

ˆ

𝑆(𝑡_𝑖+1) = ˆ𝑆(𝑡_𝑖) + 𝜇( ˆ𝑆(𝑡_𝑖)) [𝑡_𝑖+1− 𝑡_𝑖] + 𝜎 ( ˆ𝑆(𝑡_𝑖)) √

𝑡_𝑖+1− 𝑡_𝑖𝑍_𝑖+1 (2.17) Two important concepts to consider concerning time discretisations are strong and weak convergence. These concepts are used to discern both how the numerical method converges and how fast it does so.

Definition 2. Strong Convergence.

The discretisation ˆ𝑋 is said to have a strong order of convergence 𝛽 > 0 if

E || ˆ𝑋(𝑛ℎ) − 𝑋 (𝑇 ) || ≤ 𝑐ℎ𝛽 (2.18) where ℎ is the small enough time step, 𝑐 is some constant, and 𝑛 = 𝑇 /ℎ.[12]

Definition 3. Weak Convergence.

Similarly, ˆ𝑋 has a weak order of convergence 𝛽 if  _{E[ 𝑓 ( ˆ}𝑋(𝑛ℎ))] − E[ 𝑓 (𝑋 (𝑇))]   ≤ 𝑐ℎ 𝛽 (2.19) where 𝑓 represents a function in a set 𝐶𝑃2𝛽+2 which are polynomially bounded functions from

R𝑑to R with derivatives of order 0, 1, . . . , 2𝛽 + 2.[12]

In the case of the Euler Maruyama method its strong order of convergence is usually 1/2 or√4𝑡 whereas the weak order of convergence is usually 1 or 4𝑡 [8]

As the Euler methods are of the first order discretization the convergence rate is significantly slower in convergence than methods of higher order. The Milstein method is another technique for a time discretization approximation of a SDE. While it is also a first order discretization having the same weak order of convergence as the Euler-Maruyama method, its strong order of convergence is shown to be 4𝑡 which is superior to the strong order of convergence for the Euler-Maruyama method√4𝑡. [16]. Other methods can be implemented for higher order discretizations such as the Runge-Kutta method for higher orders[7], as well as [1],[2] ,[3] , where methods based on the Taylor’s decomposition are presented for numerical solutions for initial value problems of higher order.

Chapter 3

Standard Monte Carlo

As mentioned in the introduction, the Monte Carlo method as we know it today was developed by physicists working on the Manhattan project at Los Alamos laboratory during World War II.[20] However, the theory behind it dates back to the Buffon’s needle problem[6] which was formulated by Georges-Louis Leclerc, Comte de Buffon back in the 18th century. Without going into detail, the problem can be adapted to approximate 𝜋 by a crude version of the Monte Carlo method. For this thesis we are interested in generating the expected value of a European option and the Monte Carlo method is a very useful tool in finance to do so. Simply put we are interested in the expected value of the payoff of this option. If we denote the payoff as 𝜃, the crude estimator of the payoff is defined as

ˆ Θ𝑁 = 1 𝑁 𝑁 ∑︁ 𝑖=1 𝑔(𝑆_𝑖(𝑇 )) (3.1) where N is an integer representing the set of stochastic variables, (𝑖|𝑖 = 1, 2, ..., 𝑁) and 𝑔(𝑆_𝑖(𝑇 )) 𝑓 𝑜𝑟 (𝑆_𝑖|𝑖 = 1, ..., 𝑁) is the payoff function of the GBM discussed in the previous section. The expected value of the estimator is then

E[ ˆΘ𝑁] = E " 1 𝑁 𝑁 ∑︁ 𝑖=1 𝑔(𝑆(𝑇 )) # = 1 𝑁 𝑁 ∑︁ 𝑖=1 E[𝑔 (𝑆 (𝑇 ))] = 𝜃 (3.2) It is however important to note that since the GBM includes both individual paths 𝑁 and the step size for each 𝑆(𝑡) we need to modify the crude estimator to fit our specific case

ˆ Θ4𝑡𝑁 = 1 𝑁 𝑁 ∑︁ 𝑖=1 𝑔(𝑆4𝑡 𝑖 ) (3.3)

and so the expected value then becomes E[ ˆΘ4𝑡𝑁] = E " 1 𝑁 𝑁 ∑︁ 𝑖=1 𝑔(𝑆4𝑡 𝑖 ) # = 1 𝑁 𝑁 ∑︁ 𝑖=1 E[𝑔 (𝑆4𝑡𝑖 )] = E[𝑔(𝑆 4𝑡_)] (3.4) Now we note that due to Kolmogorov’s law, the strong law of large numbers[12], our crude estimator, ˆΘ4𝑡_{𝑁, approaches the expected value, E[𝑔(𝑆}4𝑡)], as 𝑁 approaches ∞.

Since we are using the Euler-Maruyama scheme to discretize the GBM, we can therefore deduce that, due to the local and global truncation errors discussed in the previous section, equation (3.4) is not equal to 𝜃. However, we can calculate the bias of the estimator using equation (3.4) 𝐵𝑖 𝑎 𝑠[ ˆΘ4𝑡 𝑁] := E _ˆ Θ4𝑡𝑁 − 𝜃 = E[𝑔(𝑆 4𝑡 )] − E[𝑔(𝑆(𝑇))] (3.5) The variance of our estimator is calculated as follows

V[ ˆΘ4𝑡𝑁] = V " 1 𝑁 𝑁 ∑︁ 𝑖=1 𝑔(𝑆4𝑡 𝑖 ) # = 1 𝑁2 𝑁 ∑︁ 𝑖=1 V[𝑔 (𝑆𝑖4𝑡)] = V[𝑔 (𝑆4𝑡)] 𝑁 (3.6)

3.1

Error estimates

Now that we have established the equations for the expected value of the crude estimator and the variance we can use them to calculate the mean square error (MSE). Since the Monte Carlo method is only an approximation of 𝜃 we need a way of establishing how efficient the estimation is. This is where the MSE is a useful technique to do so since it takes into consideration both the variance and the bias of our estimation. The MSE of our crude estimator is

𝑀 𝑆 𝐸[ ˆΘ4𝑡 𝑁] = V[𝑔 (𝑆4𝑡)] 𝑁 + (E[𝑔(𝑆 4𝑡 )] − E[𝑔(𝑆(𝑇))])2 (3.7) Due to the last term of the MSE being a squared dimension, we can take the square root of the MSE and obtain the root mean square error (RMSE) which further enhances the measurement of accuracy. See Fig(3.1) for an example of the discretisation error that occurs.

𝜖 = 𝑅𝑀 𝑆𝐸 [ ˆΘ4𝑡 𝑁] = √︃ 𝑀 𝑆 𝐸[ ˆΘ4𝑡 𝑁] = √︂ V[𝑔 (𝑆4𝑡)] 𝑁 + (E[𝑔(𝑆 4𝑡_{)] − E[𝑔(𝑆(𝑇))])}2 _(3.8)

For the sake of simplicity we can note that the RMSE (𝜖 ) is the square root of the MSE, and therefore we can denote the MSE as 𝜖2. If we now look at equation (3.7) we can notice that there are seemingly two different ways to improve the accuracy of the estimation. We can either reduce the variance by increasing the 𝑁 in our simulation, or we can reduce the discretization error by making the time steps smaller since we know from equation (2.13) that it is proportional to 4𝑡2since improving the accuracy of the LTE by definition also improves the accuracy of the GTE . However doing so comes at a price since increasing either variable means that the computational cost will also increase. For simulations requiring high accuracy this can be very time consuming and it increases the requirements of a computer’s RAM and CPU. If we want to minimize the cost of computation while still achieving an as accurate estimation as possible we can set up the following inequality statement in order to ensure that neither the bias, nor variance becomes to dominant. We know that the variance is proportional to 1/𝑁 and the bias is proportional to 4𝑡2and so

𝜖2 ≤ O (1 𝑁

Which we rearrange for both terms 𝜖2= O (

1 𝑁

), 𝜖2= O (4𝑡2)

Giving us the following proportional requirements for the most efficient accuracy of 𝜖

𝑁 = O (𝜖−2), 4𝑡 = O (𝜖 ). (3.9) Since the computational cost depends on the level of accuracy we seek, in this case the RMSE of 𝜖 we can see that the Cost, 𝐶 would be

𝐶 = O  𝑁 4𝑡  = O  𝜖−2 𝜖  = O (𝜖−3) (3.10) As can be seen in the equation above the computational complexity to perform a Monte Carlo simulation at this level of accuracy is very large. There is however a way to reduce the computational cost required quite substantially, and this is what we will discuss in the next section.

Chapter 4

Multilevel Monte Carlo

The Multilevel Monte Carlo method was first presented in 2008 by Michael B.Giles in the journal Operations Research. In his article Multilevel Monte Carlo Path Simulation[11] he showed that the MLMC can be used to improve the computational complexity of the Monte Carlo method by simulating paths at different levels of discretization. As mentioned in the previous section the computational complexity needed to achieve a RMSE of 𝜖 was O (𝜖−3) and with Giles new method he was able to bring it down to O (𝜖−2) which is the same complexity of the standard Monte Carlo method when it does not need to approximate SDE’s[11]. The basic idea of the MLMC is to use a multigrid method, which is a geometric sequence of grids that normally doubles in size in each direction compared to its previous iterations. On very fine grids the discretization error is very small, but it comes at the cost of a high computational cost. On the other hand a coarse grid is then less accurate, but at the same time more cost effective. The multigrid method then solves the equation on the fine grids by computing all the grids using corrections and by doing so achieving a higher accuracy at a lower computational cost.

As a first step we introduce the concept of levels in the MLMC method. Each separate level of the MLMC is approximated at varying accuracies where for a level 𝑙 each step size is given by

ℎ_𝑙 = 𝑀−𝑙𝑇 , 𝑙 =0, 1, . . . , 𝐿 (4.1) for the fixed integer 𝑀 ≥ 2. For the purpose of this thesis we set 𝑀 = 2 as it has the effect of decreasing the stepsize by 1/2 for each level. In the equation above 𝑙 = 0 represents the most coarse level where the step size is of the same size as the entire interval, meaning ℎ0 = 𝑇and it

logically follows that when 𝑙 = 𝐿 the simulation is at the finest level which in this case, as we are using the Euler-Maruyama method, means that we obtain a weak order of convergence 1 at the finest level, and therefore our bias at the finest level is ℎ𝐿 = O (𝜖 )We note here that this last step

size has the same requirement as the standard MC method for it to achieve maximum efficiency In the MLMC method the approximation depends on both the parameter 𝑀 and the accur-acy 𝜖 . Higher values of 𝑀 increases the precision faster for each level and as a consequence less levels are needed in the approximation, and similarly for higher sought accuracies 𝜖 more

levels are needed as as the precision required for each path increases.

By following the reasoning of Giles [11] we can set up a simple version of a Multilevel Monte Carlo method with two levels where we want to estimate E[ ˆΘ₁] for the random variable

ˆ

Θ₁, where it is more cost effective to simulate the random variable ˆΘ₀where ˆΘ₀≈ ˆΘ₁

E[ ˆΘ1] = E[ ˆΘ0] + E[ ˆΘ1− ˆΘ0]. (4.2)

We can estimate E[ ˆΘ_{0] and E[ ˆΘ}𝑙− ˆΘ𝑙−1] independently with

𝑁−1 0 𝑁₀ ∑︁ 𝑖=1 [ ˆΘ₀(𝑖)] (4.3) and 𝑁−1 1 𝑁₁ ∑︁ 𝑖=1  ˆ Θ₁(𝑖) − ˆΘ₀(𝑖)  (4.4) Combining (4.2) and (4.3) we get the Monte Carlo estimate of E[ ˆΘ₁] with two levels

ˆ 𝑌_𝑙 = 𝑁−1 0 𝑁₀ ∑︁ 𝑖=1 ˆ Θ(𝑖) 0 + 𝑁 −1 1 𝑁₁ ∑︁ 𝑖=1  ˆ Θ(𝑖) 1 − ˆΘ (𝑖) 0  . (4.5) The difference of the MLMC method and the standard MC method is that in MLMC we calculate the terms on the right hand side of (4.4) independently. We do this by using a different number of samples 𝑁 for the two terms, 𝑁0for E[ ˆΘ0] and 𝑁1for E[ ˆΘ1− ˆΘ0], where

the difference between ˆΘ₁and ˆΘ₀is denoted by ˆΘ(𝑖)

1 and ˆΘ (𝑖)

0 .The lower index is the level and

is typically denoted 𝑙 and 𝑙1where 𝑙 is the fine level and 𝑙1is the coarse level[11]. This means

that ˆΘ(𝑖)

𝑙 and ˆΘ (𝑖)

𝑙−1 comes from the same Wiener path but with different stepsizes ℎ𝑙 = 𝑀 −𝑙

𝑇 and ℎ𝑙−1= 𝑀

−(𝑙−1)_𝑇.

If we now consider the general case of a Multilevel Monte Carlo simulation where we add levels 𝑙 until we have 𝐿 levels, we can generalize equation (4.1) as

E[ ˆΘ𝐿] = E[ ˆΘ0] + 𝐿

∑︁

𝑙=1

E[ ˆΘ𝑙− ˆΘ𝑙−1]. (4.6)

where we independently estimate the terms on the right hand side as 𝑁−1 0 𝑁₀ ∑︁ 𝑖=1 ˆ Θ(𝑖) 0 + 𝐿 ∑︁ 𝑙=1 𝑁−1 𝑙 𝑁₁ ∑︁ 𝑖=1 ( ˆΘ(𝑖) 𝑙 − ˆΘ (𝑖) 𝑙−1) ! . (4.7) where the index 𝑙 represents the levels where 0 is the coarsest and 𝐿 is the finest with increasing amounts of time steps to generate greater accuracy. We now define

ˆ 𝑌₀ = 𝑁−1 0 = 𝑁₀ ∑︁ 𝑖=1 ˆ Θ(𝑖) 0 (4.8)

and ˆ 𝑌_𝑙 = 𝑁−1 𝑙 = 𝑁𝑙 ∑︁ 𝑖=1  ˆ Θ(𝑖) 𝑙 − ˆΘ (𝑖) 𝑙−1  (4.9) for 𝑙 > 0 we can write the MLMC estimator as

ˆ 𝑌 = 𝐿 ∑︁ 𝑙=0 ˆ 𝑌_𝑙. (4.10) Similarly to the bias of the standard MC method (eq 3.5), the MLMC approximation is also a biased estimation of 𝜃 according to the following formula

𝐵𝑖 𝑎 𝑠[ ˆ𝑌] = E[ ˆΘ_𝐿 − Θ] (4.11) Because of the independence of the levels we can then write the variance of the estimation as

V[ ˆ𝑌] = V 𝐿 ∑︁ 𝑙=0 ˆ 𝑌_𝑙 ! = 𝐿 ∑︁ 𝑙=0 𝑁−1 𝑙 V𝑙. (4.12)

We can minimise the variance by choosing 𝑁𝑙, which in this case is the optimal number of

samples, to be proportional to√𝑉_𝑙ℎ_𝑙[11]. As we need to consider the effect the optimal number of samples, 𝑁𝑙, has on the computational cost, combined with each path being simulated using

𝑇/ℎ𝑙timesteps, we can then formulate the total computational cost for the simulation as a sum

over all levels.

𝐿

∑︁

𝑙=0

𝑁_𝑙ℎ−1

𝑙 . (4.13)

Due to the Euler discretization, there is O (ℎ) weak convergence and O (ℎ1/2) strong convergence[23]. With a strong convergence of order 1/2 and under the assumption that the function is Lipschitz, see eq. (2) in [11], then the variance of a single sample is equal to O (ℎ𝑙). As such the optimal

𝑁_𝑙 is therefore proportional to ℎ_𝑙and if we then chose 𝑁_𝑙 = O (𝜖−2𝐿 ℎ_𝑙) we get a total variance of the estimator, ˆ𝑌, that is proportional to 𝜖2[11] as seen below.

V[ ˆ𝑌] = 𝐿 ∑︁ 𝑙=0 𝑁−1 𝑙 V𝑙. = 𝐿 ∑︁ 𝑙=0 O (𝜖2𝐿−1ℎ−1 𝑙 )O (ℎ𝑙) = 𝐿 ∑︁ 𝑙=0 O (𝜖2₎ 𝐿 = O (𝜖2).

With the weak convergence for each level 𝑙 being O (ℎ𝑙) we can deduce that the total bias (4.11)

for this method is of order O (ℎ𝐿) = O ( 𝑀

−𝐿_{). In order to reach the wanted accuracy,𝜖, we}

need the bias to be proportional to 𝜖 , and we can therefore set 𝐿 = −log

𝑀𝜖 =

log 𝜖−1

log 𝑀 = O (log 𝜖

−1_).

By applying these concepts to equation (4.13) we can see that the total computational com-plexity is proportional to 𝐿 ∑︁ 𝑙=0 𝑁_𝑙ℎ−1 𝑙 = 𝐿 ∑︁ 𝑙=0 O (𝜖−2𝐿 ℎ_𝑙) ℎ−1 𝑙 = O (𝜖 −2 𝐿2) = O (𝜖−2(log 𝜖)2) [11]

4.1

Giles Complexity Theorem

In his report [11] Giles proposes the Complexity Theorem as a new technique to calculate the computational complexity of the multilevel method depending on the underlying numerical approximations and multilevel estimators. When we are working with simulations of SDE’s it is important to note that it is infeasible to get an exact simulation of our chosen random variable. As a consequence of this the required paths, and correspondingly the accuracies of these paths, tends to require significant amounts of computational power to simulate accurately. We can quantify the computational complexity required by using the concepts we have introduced in earlier sections. Therefore, since we have the MLMC estimator

ˆ 𝑌 = 𝐿 ∑︁ 𝑙=0 ˆ 𝑌_𝑙 and E[ ˆ𝑌] = E[ ˆΘ𝐿], V[ ˆ𝑌] = 𝐿 ∑︁ 𝑙=0 𝑁−1 𝑙 V𝑙,

we can write the mean square error of the MLMC estimation by summarization of the variance and the squared bias of the estimator

𝑀 𝑆 𝐸_{𝑀 𝐿 𝑀 𝐶} = 𝐿 ∑︁ 𝑙=0 𝑁−1 𝑙 V𝑙+ (E[ ˆΘ𝐿− Θ]) 2 (4.14) Knowing this we can ensure a MSE < 𝜖2if both the variance and the squared bias are ≤ 1₂𝜖2. Giles suggests the following complexity theorem

Theorem 1. (Theorem 3.1 in [11])

Let 𝜃 denote a functional of the solution of SDE (2.3) for a given Brownian path W(t), and let ˆΘ𝑙 denote the corresponding approximation using a numerical discretisation with timestep

ℎ_𝑙 = 𝑀−𝑙𝑇.

If there exists independent estimators ˆ𝑌_𝑙 based on 𝑁𝑙 Monte Carlo samples, and positive

constants 𝛼 ≥ 1₂, 𝛽, 𝑐₁, 𝑐₂, 𝑐₃such that

1. E[ ˆΘ𝑙− 𝜃] ≤ 𝑐1ℎ 𝛼 𝑙, 2. E[ ˆ𝑌_𝑙] = ( E[ ˆΘ0], 𝑙 =0, E[ ˆΘ𝑙− ˆΘ𝑙−1], 𝑙 >0, 3. V[ ˆ𝑌_𝑙] ≤ 𝑐₂𝑁−1 𝑙 ℎ 𝛽 𝑙 ,

4. 𝐶𝑙, the computational complexity of ˆ𝑌𝑙, is bounded by 𝐶𝑙 ≤ 𝑐₃𝑁_𝑙ℎ

𝛽 𝑙,

then there exists a positive constant 𝑐4such that for any 𝜖 < 𝑒−1, there are values 𝐿 and 𝑁𝑙 for which the multilevel estimator

ˆ 𝑌 =Í𝐿

𝑙=0𝑌ˆ𝑙

has an MSE with bound MSE≡ E[( ˆ𝑌 − E[𝜃])2] < 𝜖2,

with a computational complexity 𝐶 with bound

𝐶 ≤          𝑐₄𝜖−2, 𝛽 >1, 𝑐₄𝜖−2(log 𝜖)2, 𝛽=1, 𝑐₄𝜖−2−(1−𝛽)/𝛼, 0 < 𝛽 < 1. See [11] for full proof.

An important takeaway from this theorem is that it is not specific to the case that we have presented, in general any estimator or numerical method can be used for the MLMC as long as it follows the requirements of the theorem. Some other possibilities include discontinous functions or path-dependent functions [11]. The hardest assumption to prove is the third, as it tends to vary depending on the method used for estimation. We see that 𝛽 is the most important parameter in the theorem as it is key in the bounding of both the variances 𝑉𝑙 and

the individual variance of each estimator ˆΘ𝑙. As the optimal number of realizations for each

level 𝑁𝑙 is proportional to

√

Vℎ𝑙 the computational cost per level is then correspondingly

equal to O (ℎ( 𝛽−1)/2_𝑙 ). As such the computational cost needed for the approximation is spent differently depending on the value of 𝛽. When 𝛽 < 1 most of the computational power is spent at the coarser levels and then decreases with the levels, consequently when 𝛽 > 1 most of the power is spent at the finer levels and therefore it increases with the levels, and finally when 𝛽 = 1 it is spent close to even across all levels of the grid[11]. As we are using the MLMC method with a Euler discretization we know that the weak convergence is O (ℎ) and the strong convergence is O (ℎ1/2) which gives a value of 𝛽 = 1 due to the conditions of the complexity theorem. Therefore the computational complexity necessary for this MLMC method is O (𝜖−2.5), which is a significant improvment of the standard MC method which requires a computational complexity of O (𝜖−3). For further information see [11].

Chapter 5

Numerical Results

All results in this section are simulated on an Acer Aspire A315-56 laptop with an Intel Core i5-1035G1 CPU, and a 8GB RAM. The coding and simulations are done in the MATLAB environment.

5.1

Verification of strong and weak convergence

With the explicit formula of the GBM (2.4) we can verify the strong and weak convergence for the Euler-Maruyama method numerically. In Fig(5.1) we compare the weak and strong error of the Euler simulation and the explicit value (2.4) by calculating the values at the final time for several values. This is tested several times with different time steps and we can confirm the resulting weak convergence to be of order 1 and the corresponding strong convergence to be of order 1/2.

5.2

Computational Times

The aim of the MLMC method is to decrease the computational cost, and correspondingly the computational time, of the standard MC method. By generating an algorithm which compares the computational times of the MC and MLMC method we can test this by simulating for different values of the integer 𝑀, the factor by which the timestep is refined at each level with a fixed value of 𝐿, and similarly by simulating for different values of 𝐿, the factor by which the level of accuracy is refined. For the simulation with different integers 𝑀, we fix 𝐿 to be 5, and for the simulation with different values for 𝐿, we fix 𝑀 at 214. For both simulations we set 𝑁 = 40 Fig(5.2) shows the computational times with different values for 𝑀, and Fig(5.3) shows the computational times with different values for 𝐿. In Fig(5.2) we see that, while initially the MC method is faster than the MLMC method, as 𝑀 grows larger the difference in computational times for the MC and MLMC is about 4 times less time consuming for the MLMC method where the computational time at the highest level of 𝑀 = 218, the MC method takes 321.6346 seconds compared to 81.8435 seconds for the MLMC method. In Fig(5.3) we see a similar behaviour where initially the MC method outperforms the MLMC method, but as the sought level of accuracy increases the MC method grows exponentially larger than the

Figure 5.1: Strong and weak convergence

MLMC method. For the highest level of accuracy 𝐿 = 10, the MLMC method only takes 13.8600 seconds, compared to the MC method which takes 644.2551 seconds which is almost 47 times longer than the MLMC method.

5.3

Pricing the European put option with MC and MLMC

We seek to find if the behaviour of the expected values of the European put option, equation (1.4), for the MLMC method is comparable to the more well established MC method. To compare this behaviour we generate an algorithm where we set up the required criteria for the European put option with the following parameters. Initial price 𝑆(0) = 1, Strike 𝐾 = 10, Time to maturity (days) 𝑁 = 40, riskless rate 𝑟 = 0.1 and volatility 𝜎 = 0.5. With these parameters established we can then generate the expected value of the option by setting up the GBM (2.4) and defining our 𝑀 for how many simulations we want to generate. Figures (5.4) and (5.5) shows the behaviour of the MC and MLMC for 𝑀 = 218and 𝐿 = 5, meaning we’re simulating more paths but at a lower level of accuracy, whereas figures (5.6) and (5.7) shows the behaviour when 𝑀 = 214 and 𝐿 = 10 where we simulate fewer paths but at a higher level of accuracy. As can be seen both the MC and MLMC methods behave similar to eachother where the minor difference between the two graphs is due to the statistical error that occurs when generating the Brownian motions as they are random in nature.

Figure 5.2: Computational times with different values of M

Figure 5.4: MC behaviour with 𝑀 = 218, 𝐿 = 5

Figure 5.6: MC behaviour with 𝑀 = 214, 𝐿 = 10

Finally we seek to compare the results of the expected values of the MLMC method with the generally accepted Scholes model for valuating European put options. The Black-Scholes price is calculated with the built-in MATLAB function blsprice which takes the input of the required parameters initial price, strike, riskless rate, time to maturity, and volatility, and the output of the function gives the call and put prices for the European option. For a European put option with these specific parameters, the BS-model prices the option at 8.8910 and we will use this value to compare it with the results of the MLMC method. Since we are interested in the put price, the call price is omitted in the results. As before the same values for the parameters are chosen: Initial price 𝑆(0) = 1, Strike 𝐾 = 10, Time to maturity (days) 𝑁 =40, riskless rate 𝑟 = 0.1 and volatility 𝜎 = 0.5, and we then generate the expected values for the MLMC method with our previously created algorithm. See appendix A for the codes used in this report.

L 2 3 4 5 6 7 8 10 8,7442 8,5150 8,2920 8,1649 8,0185 7,8547 7,5670 11 8,7747 8,7046 8,6733 8,5435 8,5188 8,5383 8,6572 12 9,0120 8,9596 8,8558 8,8848 8,9545 8,8585 8,7116 M = 2𝑥 13 9,0102 8,8558 8,8773 8,8349 8,8558 8,9016 8,9861 14 8,7108 8,6824 8,6758 8,5903 8,6987 8,7314 8,7087 15 8,9382 9,0047 8,8976 8,9718 9,0489 9,0014 8,9987 16 8,8826 8,8387 8,7695 8,8833 8,7604 8,8145 8,7541 17 8,7256 8,7256 8,8835 8,8541 8,6972 8,7773 8,7068 18 8,9461 8,9457 8,9196 8,8819 8,8757 8,9442 8,9744

Table 5.1: Expected values of European put option with MLMC for different levels 𝑀 and 𝐿

L 2 3 4 5 6 7 8 10 0,02155 0,14138 0,35880 0,52722 0,76126 1,07392 1,75298 11 0,01353 0,03474 0,04739 0,12076 0,13853 0,12440 0,05466 12 0,01464 0,00471 0,00124 0,00004 0,00403 0,00106 0,03218 𝑀 =2𝑥 13 0,01421 0,00124 0,00019 0,00315 0,00124 0,00011 0,00904 14 0,03247 0,04351 0,04631 0,09042 0,03698 0,02547 0,03323 15 0,00223 0,01293 0,00004 0,00653 0,02493 0,01219 0,01160 16 0,00007 0,00274 0,01476 0,00006 0,01706 0,00585 0,01874 17 0,02736 0,02736 0,00006 0,00136 0,03756 0,01293 0,03393 18 0,00304 0,00299 0,00082 0,00008 0,00023 0,00283 0,00696 Table 5.2: MSE’s of European put option with MLMC for different levels 𝑀 and 𝐿

L 2 3 4 5 6 7 8 10 0,1410 0,2056 0,2884 0,5208 0,7957 1,1837 1,9781 11 0,2359 0,3622 0,4956 0,7056 0,9501 1,4780 2,3284 12 0,4455 0,5883 0,8801 1,0908 1,4679 2,0862 2,9050 𝑀 =2𝑥 13 0,935 1,2027 1,6072 2,0756 2,4999 3,2384 4,2038 14 1,8127 2,4063 3,0216 3,6731 4,6236 5,7475 7,1224 15 3,6215 5,1922 6,1072 7,3589 8,7815 11,0079 12,4160 16 6,9665 9,9194 12,6521 14,9887 18,1139 20,9563 26,4706 17 14,4984 16,1690 28,2441 33,0084 40,7868 43,1077 52,5120 18 43,9855 98,3635 75,6636 93,337 104,4432 113,7656 140,0623 Table 5.3: Computational Times (in seconds) of European put option with MLMC for different levels 𝑀 and 𝐿

The results of Table(5.2) suggests that for a sufficiently large enough 𝑀, the MLMC method in general performs very well in estimating the price for a European put option. Especially simulations with a level of accuracy 𝐿 = 4 and 𝐿 = 5 seems to generate very accurate results when ran with high levels of 𝑀. Since the computational cost quickly grows with higher levels of 𝑀 and 𝐿. Table(5.3) compares the different computational times for each simulation. We see that at the highest level of 𝑀, the computational time is quite large, especially for higher levels of 𝐿. However, the result of the simulation where 𝑀 = 216and 𝐿 = 5 seems to generate both a very accurate result and a reasonable computational time, especially compared to a standard MC method which would take significantly longer to simulate with such precision.

Chapter 6

Conclusion

The aim of this thesis was to research if it is possible to implement a working MLMC simulation on European put options. We first gave some historical background on options, and introduced relevant theory necessary for developing a MLMC simulation. The theory behind the MLMC method suggests that it can greatly improve the computational cost and complexity of running simulations compared to the standard MC method. It is shown that with an Euler-Maruyama time discretization the MLMC method converges at a faster rate than the MC method and as a result it requires less simulations in order to reach the required accuracy. The results of the numerical analysis suggests that it is indeed possible to implement the MLMC method to European put options. Further research needs to be conducted on ways to improve the computational costs and complexity even further. The Milstein scheme for time discretisation enhances the convergence rate faster than the Euler-Maruyama method and could be of interest to look into. Other methods could include higher order time discretizations or adapting quasi Monte Carlo methods in combination with the MLMC method to further improve computational times but it is beyond the scope of this thesis.

Bibliography

[1] Doghonay Arjmand. Highly Accurate Difference Schemes for the Numerical Solution of

Third-Order Ordinary and Partial Differential Equations.Royal Institute of Technology, School of Computer Science and Communication, Stockholm, ISSN-1653-5715, 2010 [2] Allaberen Ashyralyev & Doghonay Arjmand. A note on the Taylor’s decomposition on

four points for a third-order differential equation.Applied mathematics and computation, vol 188(2). ISSN 0096-3003

[3] Allaberen Ashyralyev, Doghonay Arjmand & Muhammet Koksal. Taylor’s decomposition

on four points for solving third-order linear time-varying systemsJournal of the Franklin Institute, vol. 346, Elsevier, pp. 651-662, 2009

[4] Robert E. Bradley & C. Edward Sandifer. Leonhard Euler: Life, Work and Legacy (Studies

in the History and Philosophy of Mathematics, Volume 5Elsevier Science, Amsterdam, 2007.

[5] Paolo Brandimarte. Numerical Methods in Finance and Economics: A MATLAB-Based

IntroductionJohn Wiley & Sons, New York, pp. 3-21, 2006.

[6] Georges. L. L. Buffon Histoire Naturelle, generale et particuilere, servant de suite a

la Theorie de la Terre et d’introduction a l’Historie des Mineraux... Supplement Tome premier [septieme], vol. 4de l’omprimiere royale, Paris, 1777.

[7] John C. Butcher. Numerical methods for ordinary differential equations in the 20th

centuryJournal of Computational and Applied Mathematics, Elsevier, vol. 125, pp. 1-29, 2000

[8] Jesper Carlsson, Kyong-Sook Moon, et al. Stochastic Differential

Equations: Models and Numerics.Lecture notes used at KTH, 2009. [9] Josef de la Vega. Confusion de Confusiones 1688.

[10] George Crawford & Bidyut Sen. Derivatives For Decision Makers

Strategic Management IssuesJohn Wiley & Sons, New York, pp. 7-9, 1996. [11] Michael B. Giles. Multilevel Monte Carlo path simulation. In Acta

[12] Paul Glasserman. Monte Carlo Methods in Financial Engineering Springer, New York, Vol 53, 2003

[13] John C. Hull. Options, futures, and other derivatives. Prentice Hall, New York, 8. ed., 2012.

[14] Nabil Kahalé. General multilevel Monte Carlo methods for pricing

discretely monitored Asian Options.European Journal of Operational Research, Vol. 287, Elsevier B.V, pp. 739-748. 2020.

[15] Jack P.C. Kleijnen, Ad A.N. Ridder et al. Variance Reduction Techniques in Monte Carlo

MethodsSSRN Electronic Journal, No. 2010-117, 2010.

[16] Peter E. Kloeden & Eckhard Platen. Numerical Solution of Stochastic Differential

Equa-tionsSpringer, Berlin, 1999

[17] Robert W. Kolb & James A. Overdahl. Financial Derivatives John Wiley & Sons, New York, 3. ed., 2003.

[18] Funmilayo E. Kazeem. Multilevel Monte Carlo Simulation in Options Pricing. Master Thesis, University of the Western Cape, South Africa, November, 2014.

[19] Jerry W. Markham. A Financial History of the United States M.E. Sharpe, New York, p 52, 2002.

[20] N. Metropolis. The Beginning of the Monte Carlo Method Los Alamos Science, Special Issue, 1987.

[21] Ulrike Schaede. Forwards and futures in tokugawa-period Japan: A new perspective on

the D ¯ojima rice market.Journal of Banking & Finance, Vol. 13, Issues 4-5, Elsevier, pp. 487-513, 1989.

[22] René L. Schilling & Lothar Partzsch. Brownian Motion: An Introduction to Stochastic

ProcessesDe Gruyter, Berlin, 2nd ed., 2014

[23] Denis Talay & Luciano Tubaro. Expansion of the global error for numerical schemes

solv-ing stochastic differential equationsRapports de Recherche, Vol. 1069, Institut National de Recherche en Informatique et en Automatique, 1989.

[24] Jochen Voss. An Introduction to Statistical Computing John Wiley & Sons, New York, 2014.

Appendix A

MATLAB Codes

A.1

Convergence check

% S t r o n g and weak c o n v e r g e n c e f o r t h e E u l e r method s t e p s = [ 1 : 6 ] ; f o r i = s t e p s N = 2^ i % number o f t i m e s t e p s r a n d n ( ’ s t a t e ’ , 0 ) ; T = 1 ; d t = T / N ; t = 0 : d t : T ; r = 0 . 1 ; s i g m a = 0 . 5 ; S0 = 1 0 0 ; M = 1E6 ; % number o f r e a l i s a t i o n s S = S0∗ ones (M, 1 ) ; % S ( 0 ) f o r a l l r e a l i z a t i o n s W = z e r o s (M, 1 ) ; % W( 0 ) f o r a l l r e a l i z a t i o n s f o r j = 1 :N dW = s q r t ( d t )∗ randn (M, 1 ) ; % Wiener i n c r e m e n t s S = S + S .∗ ( r ∗ d t +sigma ∗dW) ; % p r o c e s s e s a t n e x t time s t e p W = W + dW; % B r o w n i a n p a t h s a t n e x t s t e p end

ST = S0∗ exp ( ( r −sigma ^ 2 / 2 ) ∗ T + sigma ∗W ) ; % e x a c t f i n a l v a l u e w E r r o r ( i ) = mean ( S−ST ) ; % weak e r r o r s E r r o r ( i ) = s q r t ( mean ( ( S−ST ) . ^ 2 ) ) ; % s t r o n g e r r o r end d t = T . / ( 2 . ^ s t e p s ) ; l o g l o g ( d t , a b s ( w E r r o r ) , ’ o− − ’ , dt , dt , ’ − − ’ , dt , abs ( s E r r o r ) , ’ o − ’ , dt , s q r t ( d t ) ) t i t l e ( ’ C o n v e r g e n c e R a t e s ’ ) x l a b e l ( ’ T i m e s t e p s ’ ) y l a b e l ( ’ E x p e c t e d Value ’ ) l e g e n d ( ’ Weak E r r o r ( S i m u l a t e d ) ’ , ’ Weak E r r o r ( E x p l i c i t ) ’ , . . . ’ S t r o n g E r r o r ( S i m u l a t e d ) ’ , ’ S t r o n g E r r o r ( E x p l i c i t ) ’ , ’ l o c a t i o n ’ , ’ b e s t ’ ) g r i d on

A.2

MLMC function

f u n c t i o n o u t p u t = MLMC_Code (M, N , s , t 0 ) r a n d n ( ’ s t a t e ’ , 0 ) ; K = 1 0 ; T = 1 ; d t = ( T− t 0 ) / N; t = t 0 : d t : T ; r = 0 . 1 ; s i g m a = 0 . 5 ; S = s∗ ones (M, 1 ) ; % S ( 0 ) f o r a l l r e a l i z a t i o n s W = z e r o s (M, 1 ) ; % W( 0 ) f o r a l l r e a l i z a t i o n s f o r j = 1 :N dW = s q r t ( d t )∗ randn (M, 1 ) ; % Wiener i n c r e m e n t s S = S + S .∗ ( r ∗ d t +sigma ∗dW) ; % p r o c e s s e s a t n e x t time s t e p W = W + dW; % B r o w n i a n p a t h s a t n e x t s t e p end

o u t p u t = mean ( exp (− r ∗(T− t 0 ) ) ∗ max (K−S ) ) ; %Expected Value European Put Option end

A.3

Computational times

%% f u n c t i o n f o r C o m p u t a t i o n a l t i m e o f MLMC & MC f u n c t i o n Comp_Time = MeshMC (M0) N0 = 4 0 ; %T i m e s t e p s L= 4 ; %L e v e l s o f a c c u r a c y s = 1 ; % I n i t i a l p r i c e t 0 = 0 ; % I n i t i a l t i m e A r r a y _ s = l i n s p a c e ( 0 , 2 , 1 0 ) ; %v e c t o r o f s A r r a y _ t 0 = l i n s p a c e ( 0 , 0 . 9 , 1 0 ) ; %v e c t o r o f t 0 f = z e r o s ( l e n g t h ( A r r a y _ s ) , l e n g t h ( A r r a y _ t 0 ) ) ; %z e r o− m a t r i x f o r Option p r i c e t i c %s t a r t i n g c l o c k f o r MLMC c o u n t e r 1 = 1 ; f o r s = A r r a y _ s %MLMC l o o p c o u n t e r 2 = 1 ; f o r t 0 = A r r a y _ t 0 ,

Ep0 = MLMC_Code (M0, N0 , s , t 0 ) ; % c a l l i n g MLMC f u n c t i o n f o r f i r s t s t e p Sum = Ep0 ;

f o r k = 1 : L , %R e f i n i n g g r i d s f o r d i f f e r e n t l e v e l s L

E p 1 _ f = MLMC_Code (M0 / 2 ^ ( k ) , 2 ^ ( k )∗N0 , s , t 0 ) ; %R e s u l t of f i n e g r i d

Ep0_c = MLMC_Code (M0 / 2 ^ ( k ) , 2 ^ ( k−1)∗N0 , s , t 0 ) ; %R e s u l t of c o a r s e g r i d Sum = Sum + E p 1 _ f − Ep0_c ; %A s s i g n i n g new Sum

end f ( c o u n t e r 1 , c o u n t e r 2 ) = Sum ; %A s s i g n i n g Sum t o f u n c t i o n f o r o p t i o n p r i c e c o u n t e r 2 = c o u n t e r 2 + 1 ; end c o u n t e r 1 = c o u n t e r 1 + 1 ; end Comp_Time ( 1 , 1 ) = t o c ; %E n d i n g c l o c k f o r MLMC and s t o r i n g t h e r e s u l t %% R e p e a t i n g p r e v i o u s p r o c e s s b u t f o r t h e MC method A r r a y _ s _ 2 = l i n s p a c e ( 0 , 2 , 1 0 ) ; A r r a y _ t 0 _ 2 = l i n s p a c e ( 0 , 0 . 9 , 1 0 ) ; f 2 = z e r o s ( l e n g t h ( A r r a y _ s _ 2 ) , l e n g t h ( A r r a y _ t 0 _ 2 ) ) ; t i c c o u n t e r 3 = 1 ; f o r s = A r r a y _ s _ 2 c o u n t e r 4 = 1 ; f o r t 0 = A r r a y _ t 0 _ 2 Sum2 = MLMC_Code (M0, 2 ^ L∗N0 , s , t 0 ) ; f 2 ( c o u n t e r 3 , c o u n t e r 4 ) = Sum2 ; c o u n t e r 4 = c o u n t e r 4 + 1 ; end

c o u n t e r 3 = c o u n t e r 3 + 1 ; end Comp_Time ( 1 , 2 ) = t o c ; %% Comp t i m e w i t h f i x e d L Array_M0 = 2 . ^ [ 6 : 1 8 ] ; %g e n e r a t e d i f f e r e n t l e v e l s o f M o u t p u t = z e r o s ( l e n g t h ( Array_M0 ) , 2 ) ; j = 1 ; f o r M0=Array_M0 , %R u n n i n g l o o p on c o m p u t a t i o n a l t i m e s f o r e a c h l e v e l M o u t p u t ( j , : ) = MeshMC (M0 ) ; j = j + 1 ; end y1 = o u t p u t ( : , 1 ) ; y2 = o u t p u t ( : , 2 ) ; x = Array_M0 ; l o g l o g ( x , y1 , ’ o− − ’ ,x , y2 , ’ d − − ’) x l a b e l ( ’M’ ) y l a b e l ( ’ C o m p u t a t i o n a l Time ’ ) l e g e n d ( ’MLMC’ , ’MC’ , ’ l o c a t i o n ’ , ’ b e s t ’ ) g r i d on

Appendix B

Bachelor Degree Objectives

In the document "How to Write a Thesis?" by Anatoliy Malyarenko, Dmitrii Silvestrov and Sergei Silvestrov, objectives for a successful bachelor thesis are presented. In this section we present our fullfillments of these objectives.

1. Objective 1: The student should demonstrate knowledge and understanding in the major

field of study, including knowledge of the field’s scientific basis, knowledge of applicable methods in the field, specialization in some part of the field and orientation in current research questions.

This thesis satisfies this objective as we have done a thorough litterature review on key financial concepts, and especially on Monte Carlo and Multilevel Monte Carlo methods. We also provide orientation for current research questions on how to improve these methods.

2. Objective 2: The student demonstrates the ability to search, collect, evaluate and

critically interpret relevant information in a problem formation and to critically discuss the phenomena, problem formulation and situations.

Numerous references have been investigated throughout the process of this thesis to find suitable information relevant to the topic. My work reflects this thorough investigation and I’ve successfully discussed the relevant problems and situations.

3. Objective 3: The student should demonstrate the ability to independently identify,

formulate and solve problems and to perform tasks within specified time frames.

I have solved the problems i have met throughout the work of this thesis. I have also performed the tasks given to me by my supervisor Doghonay Arjmand successfully and on time for the deadline of this thesis.

4. Objective 4: The student should demonstrate the ability to present orally, in writing,

and discuss information, problems and solutions in dialogue with different groups.

I have had weekly meetings with my supervisor throughout the work of this thesis where we have had useful discussions about the topics of this report. This thesis serves as proof of my written presentation skills and I intend to prove my ability to present orally during the presentation.

5. Objective 5: The student should demonstrate ability in the major field of study, make

judgments with respect to scientific, societal and ethical aspects.

I have done my absolute best to make sure that any work of another individual is acknowledged and referenced to throughout the entirety of this thesis.