School of Education, Culture and Communication
Division of Mathematics and Physics
BACHELOR THESIS IN MATHEMATICS / APPLIED MATHEMATICS
Deterministic Quadrature Formulae for the Black–Scholes Model
by
Timo Kudljakov and Sajedeh Saadat
Kandidatarbete i matematik / tillämpad matematik
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
Bachelor thesis in mathematics / applied mathematics
Date:
2021-06-02
Project name:
Deterministic Quadrature Formulae for the Black–Scholes Model
Author(s):
Timo Kudljakov and Sajedeh Saadat
Version: 9th June 2021 Supervisor: Anatoliy Malyarenko Reviewer: Jean-Paul Murara Examiner: Christopher Engström Comprising: 15 ECTS credits
Abstract
There exist many numerical methods for numerical solutions of the systems of stochastic differential equations. We choose the method of deterministic quadrature formulae proposed by Müller–Gronbach, and Yaroslavtseva in 2016. The idea is to apply a simplified version of the cubature in Wiener space. We explain the method and check how good it works in the simplest case of the classical Black–Scholes model.
Keywords: Deterministic quadrature formulae; Stochastic differential equation; Black– Scholes model; Discretization method.
Acknowledgements
First and foremost, we have to thank our research supervisor Prof. Anatoliy Malyarenko. Without his assistance and dedicated involvement in every step throughout the process, this thesis would have never been accomplished. The door to Prof. Anatoliy Malyarenko office was always open whenever we ran into a trouble spot or had a question about our research or writing. He consistently allowed this thesis to be our own work, but steered us in the right direction whenever he thought we needed it. We would like to thank you very much for your support and understanding over these past six months.
Our great gratitude to our examiner Christopher Engström for his insightful suggestions for further improvements for this thesis.
Contents
1 Introduction 6
1.1 The Formulation of the Task . . . 6
1.2 Review of Literature . . . 6
2 Theoretical Considerations 9 2.1 The Idea . . . 9
2.2 The Black–Scholes Model . . . 9
2.2.1 The Concept of an Option . . . 10
2.2.2 The Black–Scholes Model . . . 10
2.2.3 The Black–Scholes formulas . . . 11
2.3 Quadratures . . . 12
2.3.1 First Quadrature Method . . . 12
2.3.2 Second Quadrature Method . . . 13
2.4 The Description of the Algorithm . . . 13
3 Application to the Black-Scholes model 17 3.1 The Non-uniform Time Discretization . . . 17
3.2 The Itô–Taylor Transition . . . 17
3.3 The Support Diameter Reduction Method . . . 20
3.4 The Support Cardinality Reduction Method . . . 20
3.5 Realization . . . 21
4 Numerical Results 24 4.1 Comparison of the Initial Price of the Underlying Stock . . . 25
4.2 Comparison of the Exercise Price . . . 25
4.3 Comparison of the Risk-Free Interest Rate . . . 26
4.4 Comparison of the Maturity of the Option . . . 27
4.5 Comparison of the Volatility of the Stock per Annum . . . 27
4.6 Summary . . . 27
5 Conclusions 29 5.1 Conclusions . . . 29
5.2 Further Research . . . 29
A MATLAB Simulation 32
B Required Objectives for the Thesis 37
B.1 Objective 1: Knowledge and understanding of the main field . . . 37
B.2 Objective 2: Ability to search, collect and interpret . . . 37
B.3 Objective 3: Identify, formulate and solve problems . . . 37
B.4 Objective 4: Ability to present the report by writing and orally, so others
understand it . . . 38
List of Figures
4.1 Comparison of the initial price of the underlying stock . . . 25
4.2 Comparison of the exercise price . . . 26
4.3 Comparison of the risk-free interest rate . . . 26
4.4 Comparison of the maturity of the option . . . 27
List of Tables
2.1 Notation . . . 11
4.1 The parameters of the option . . . 24
Chapter 1
Introduction
1.1
The Formulation of the Task
Let {𝑋𝑡: 0 ≤ 𝑡 ≤ 𝑇 } be the solution of the stochastic integral equation:
𝑋𝑡 = 𝑥 + ∫ 𝑡 0 𝜇( 𝑋𝑠) 𝑑𝑠 + ∫ 𝑡 0 𝜎( 𝑋𝑠) 𝑑𝑊𝑠,
where the first integral is a Riemann one and the second one is called stochastic or Itô integral. This equation describes the dynamics of the risky security price. Let 𝑓 : ℝ → ℝ be the discounted payoff of a financial derivative.
In financial engineering, a principal problem is to calculateE[ 𝑓 ( 𝑋𝑇)], the no-arbitrage price
of the given derivative. In this matter we can use a numerical method based on a discretisation scheme, which is a set of random variables {𝑋𝑡(𝑚) : 𝑡 = 𝑡𝑖} with 0 = 𝑡0 < 𝑡1 < ... < 𝑡𝑛 = 𝑇,
where 𝑚 is the number of random variables.
We would like to use a sophisticated deterministic algorithm for calculation the expectation
of the random variable 𝑋𝑡 and apply it in the Black–Scholes model. More explanation about
this algorithm can be found in Müller-Gronbach and Yaroslavtseva (2016).
1.2
Review of Literature
In recent decades, stochastic differential equations (SDEs) have been appeared in many ap-plications related to engineering, see Henderson and Plaschko (2006), economics and finance, see Pagès (2018), where noise or random perturbations are important. Typically, compared to ordinary differential equations, SDEs contains a stochastic term (random variable) which is often the derivative of the Brownian motion or the Wiener process.
To solve SDEs, several methods have been developed both analytically (Karatzas and Shreve, 1998; Nualart, 2006; Øksendal, 2003) and numerically (Kloeden et al., 1994; Särkkä and Solin, 2019). However, in most cases with more complex coefficients and stochastic terms, explicit solutions are not available and numerical solution methods are needed. In favor of computers, the most efficient and applicable approach to solve SDEs is the simulation of sample paths of time discrete approximations. In this approach, the time interval [0, 𝑇 ] is discretized
and the approximate value of the sample paths is computed recursively at the discretization times. The statistical tools can then be employed to determine how good the approximation is and whether it is enough close to the exact solution.
Itô–Taylor series (Kloeden et al., 1994) are an extension of the Taylor series of ODEs to SDEs. The derivation is basically identical to the Taylor series solution except that the time derivative computations is replaced with application of the Itô formula. The simplest time discrete approximation of SDEs is the Euler–Maruyama approximation which utilizes only the first two terms in the simple Taylor expansion.
With time discretization 0 = 𝜏0 < 𝜏1 < ... < 𝜏𝑁 = 𝑇 of interval [0, 𝑇 ] and step size
𝛿= 𝑇 /𝑁
𝑋𝑛+1 = 𝑋𝑛+ 𝑎 ( 𝑋𝑛)Δ𝑛+ 𝑏 ( 𝑋𝑛)Δ𝑊𝑛. (1.1)
For 𝑛 = 0, 1, . . . , 𝑁 − 1 and initial value of 𝑋0= 𝑥0where
Δ𝑛= 𝜏𝑛+1− 𝜏𝑛 = 𝛿,
Δ𝑊𝑛 = 𝑊𝜏𝑛+1 − 𝑊𝜏𝑛.
(1.2)
The random variable Δ𝑊𝑛 is normally distributed with mean E(Δ𝑊𝑛) = 0 and variance
E( (Δ𝑊𝑛) 2) = Δ
𝑛.
A natural way of classifying time discrete numerical methods for SDEs is to compare them with strong and weak Taylor approximations (Kloeden et al., 1994).
Definition 1. A discrete time approximation 𝑌 converges strongly with order 𝛾 > 0 to 𝑋 at time 𝑇 if there exists a positive constant 𝐶, which does not depend on 𝛿, and a 𝛿0 > 0 such that
𝜖(𝛿) =E(| 𝑋𝑇 − 𝑌 (𝑇 ) |) ≤ 𝐶𝛿𝛾. (1.3)
for each 𝛿 ∈ (0, 𝛿0).
Definition 2. Moreover, a time discrete approximation 𝑌 converges weakly with order 𝛽 > 0 to 𝑋 at time 𝑇 as 𝛿 ↓ 0 if for each polynomial 𝑔 there exists a positive constant 𝐶, which does not depend on 𝛿, and a finite 𝛿0> 0 such that
|E(𝑔( 𝑋𝑇)) −E(𝑔(𝑌 (𝑇 ))) | ≤ 𝐶𝛿 𝛽
. (1.4)
for each 𝛿 ∈ (0, 𝛿0).
For instance, the Euler approximation usually converges with weak order 𝛽 = 1, while it converges with the strong order of 𝛾 = 0.5.
It is possible to derive higher order schemes by including even more terms from the stochastic Taylor expansions. These are usually more of theoretical than practical interest because they involve higher order derivatives which may be complicated.
The desired order of convergence determines where the Taylor expansion must be truncated. However, the weak convergence criterion only concerns probabilistic aspects of the sample path and not the sample path itself. Therefore, for a certain degree of convergence, the required number of terms of the expansion is less for the case of weak convergence than for the case of strong convergence if a certain degree of convergence is desired.
The order 2.0 weak Taylor scheme is given by 𝑌𝑛+1= 𝑌𝑛+ 𝜇Δ + 𝜎Δ𝑊 + 1 2𝜎 𝜎 0 [(Δ𝑊 )2− Δ]+ 𝜇0𝜎Δ𝑍 + 1 2( 𝜇𝜇 0+ 1 2𝜎 2𝜇00 )Δ2+ ( 𝜇𝜎0+ 1 2𝜎 2 𝜎00) [Δ𝑊 Δ − Δ𝑍], (1.5)
where 𝜇 and 𝜎 are drif and difussion respectively and evaluated at 𝑌𝑛 and Δ𝑍 is a random
variable representing the double stochastic integral.
Weak Taylor scheme is simpler, even though the degree of convergence is higher, see Lindström et al. (2015).
Chapter 2
Theoretical Considerations
2.1
The Idea
Consider a 𝑑-dimensional system of stochastic differential equations.
d𝑋 (𝑡) = 𝑎( 𝑋 (𝑡)) d𝑡 + 𝑏( 𝑋 (𝑡)) d𝑊 (𝑡), 𝑡 ∈ [0, 1],
𝑋(0) = 𝑥0
(2.1)
with the initial value 𝑥0 ∈ ℝ𝑑, drift coefficient 𝑎 : ℝ𝑑 → ℝ𝑑, diffusion coefficient 𝑏 : ℝ𝑑 →
ℝ𝑑×𝑚and an 𝑚-dimensional driving Brownian motion 𝑊 = (𝑊 (𝑡))𝑡∈[0,1]as well as the function
𝑓: ℝ𝑑 → ℝ.
Our computational task is to approximate the expected value
𝐼(𝑥0, 𝑎, 𝑏, 𝑓) = E[ 𝑓 ( 𝑋 (1))] (2.2)
by means of a deterministic method that uses 𝑥0and finitely many evaluations of 𝑎, 𝑏 and 𝑓
and their derivatives in ℝ𝑑
. The problem’s formulation is the same as in Müller-Gronbach and Yaroslavtseva (2016, Section 1). In equation (2.2) 𝑋 (1) is a random variable. We can consider 𝑋 (1) to be a price of the risky security. In that case, function 𝑓 becomes the payoff of the financial instrument. Our problem is to calculate the expected value. Müller-Gronbach and Yaroslavtseva (2016) propose a new method, which is suitable for any market models of this type. However we will check this method in this thesis using one of the simplest market models, the Black–Scholes model.
2.2
The Black–Scholes Model
In this section, the idea of basic financial derivatives and the Black–Scholes model origin and equations will be introduced. All definitions and formulas defined in this section for the Black–Scholes model are taken from Hull (2018).
2.2.1
The Concept of an Option
Options are financial derivatives that can be traded either on exchange or in the over-the-counter market. There exist two types of options. The first type is called the call option, and the second type is called the put option. A call option gives the holder the right to buy the underlying asset by specific date/time for a certain price. A put option gives the holder the right to sell the underlying asset by a specific date/time for a certain price. The contract price is known as the
exercise priceor strike price. The date/time in the contract is also known as the expiration date or maturity. Options are divided into distinctive categories: American options and European
options. The only difference between European options and American options is that American options can be exercised at any time before maturity, and European options can be exercised only at the time of maturity. In this paper we will only focus on European options.
2.2.2
The Black–Scholes Model
The Black–Scholes model has been developed at the beginning of the 1970s by the math-ematicians Fischer Black, Myron Scholes, and Robert Merton. The Black–Scholes model simulates the dynamics of a financial market containing derivative financial instruments such as options, futures, forwards, and swaps. The main idea behind the model is to hedge the options in an investment portfolio by buying and selling the underlying asset properly and eliminating risk. In 1997 for the development of the model, Robert Merton and Myron Scholes were awarded the Nobel prize for Economics.
The Black-Scholes stochastic differential equation looks like this
d𝑆(𝑡) = 𝑟 𝑆(𝑡) d𝑡 + 𝜎𝑆(𝑡) d𝑊 (𝑡), 𝑡 ∈ [0, 𝑇 ],
𝑆(0) = 𝑆0.
(2.3) As we can see equation (2.3) is similar to equation (2.1). In more detailed notation, the system (2.1) looks like this:
𝑑 𝑋𝑖(𝑡) = 𝑎𝑖( 𝑋1(𝑡), ..., 𝑋𝑑(𝑡))𝑑𝑡 +
𝑚
∑︁
𝑗=1
𝑏𝑖 𝑗( 𝑋1(𝑡), ..., 𝑋𝑑(𝑡))𝑑𝑊𝑗(𝑡)
where 1 ≤ 𝑖 ≤ 𝑑 and 1 ≤ 𝑗 ≤ 𝑚. The symbol 𝑑 means a positive integer, the number of variables in the market model. In the Black–Scholes model, 𝑑 = 1, and the variable is called
𝑆, the price of risky security. The symbol 𝑚 means another positive integer, the number of
Brownian motionsin the model. In the Black-Scholes model we have 𝑚 = 1. In Black-Scholes function a is constant which is equal to the risk-free rate that is notated by 𝑟. Next constant in the function is b which represents volatility and is notated by 𝜎. The 𝑓 denotes a function of
𝑑 real variables, the discounted payoff of a financial derivative. For example in our case the
European call option in the Black-Scholes model, we have
𝑓(𝑥) = 𝑒−𝑟𝑇max{𝑥 − 𝐾, 0} (2.4)
The random variable 𝑓 ( 𝑋 (1)) is the discounted payoff of the financial derivative of European
𝑐 The European call option price
𝑝 The European put option price
𝑆0 Stock price at time zero
𝐾 Strike price
𝑟 Continuously compounded risk-free rate
𝜎 Stock price volatility
𝑇 Time to maturity of the option
Table 2.1: Notation
there exist two ways of solving this problem. One way is general, this is exactly how Müller-Gronbach and Yaroslavtseva (2016) formulate their general formulation. Our formulation is to take proposed algorithm by Müller-Gronbach and Yaroslavtseva (2016) and to study this algorithm how it works and later apply this algorithm to the Black–Scholes model European call option.
2.2.3
The Black–Scholes formulas
The Black–Scholes–Merton formulas used for the pricing of European call and European put options are: 𝑐 = 𝑆0𝑁(𝑑1) − 𝐾 𝑒−𝑟𝑇𝑁(𝑑2) and 𝑝= 𝐾 𝑒−𝑟𝑇𝑁(−𝑑2) − 𝑆0𝑁(−𝑑1) where 𝑑1= ln(𝑆0 𝐾 + ( 𝑟+𝜎2 2 )𝑇 ) 𝜎 √ 𝑇 and 𝑑2 = ln(𝑆0 𝐾 + ( 𝑟−𝜎2 2 )𝑇 ) 𝜎 √ 𝑇 = 𝑑1− 𝜎 √ 𝑇 ,
Where the notation is given in Table 2.1.
𝑁(𝑥) stands for the cumulative probability function for a variable with standard normal
distribution. That is to say that, it is the probability that a variable with standard normal distribution will be smaller than 𝑥.
2.3
Quadratures
Our task is to calculate the expected value of the European call option. We know that the expected value is integral. Our task in this paper is to calculate that integral. There exist a group of mathematical methods, which are called quadrature methods so that approximate numerical values of integrals. We explain the method proposed by Müller-Gronbach and Yaroslavtseva (2016). The best way to do so is by explaining the simplest quadrature methods in more detail. Then we will describe how this more advanced and complicated method that Müller-Gronbach and Yaroslavtseva (2016), proposed works. We will compare the two methods and demonstrate their similarities and differences between them.
2.3.1
First Quadrature Method
Let us begin with continuous function 𝑓 : [𝑎, 𝑏] → ℝ. Now, let’s consider the points 𝑎 = 𝑥0 < 𝑥1< ... < 𝑥𝑛 = 𝑏.
Time to choose an arbitrary point 𝜉𝑖in the interval [𝑥𝑖, 𝑥𝑖+1] and put
𝐼𝑛(𝑎, 𝑏, 𝑓 ) =
𝑛−1
∑︁
𝑖=0
𝑓(𝜉𝑖) (𝑥𝑖+1− 𝑥𝑖). (2.5)
We know from single variable calculus that when 𝑛 increases without bound and at the same time the length of the longest interval
max
0≤𝑖≤𝑛−1
(𝑥𝑖+1− 𝑥𝑖)
goes to 0, then the sequence { 𝐼𝑛(𝑎, 𝑏, 𝑓 ) : 𝑛 ≥ 1 } converges to a finite limit. A finite limit is
an interval that squeezes closer and closer to. This limit is denoted by∫ 𝑏
𝑎
𝑓(𝑥) 𝑑𝑥 and is called the Riemann integral. The Riemann integral can be evaluated by the fundamental theorem of calculus or approximated by numerical integration. Moreover, the value of the limit does not
depend on the choice of the above sequence and the points 𝜉𝑖.
Now let’s consider the finite set 𝑆 = { 𝜉𝑖: 0 ≤ 𝑖 ≤ 𝑛 − 1 }, and let 𝔉 be the set of all subsets
of 𝑆. The function 𝑄𝑛: 𝔉 → ℝ given by
𝑄𝑛( 𝐴) = ∑︁
𝜉𝑖∈ 𝐴
(𝑥𝑖+1− 𝑥𝑖), 𝐴 ∈ 𝔉
is an example of a measure with finite support 𝑆. In case if 𝑏 − 𝑎 = 1, then the measure 𝑄𝑛 is
probabilistic, that is, 𝑄𝑛 ∈ M1(ℝ1).
Equation (2.5) takes the new form of
𝐼𝑛(𝑎, 𝑏, 𝑓 ) =
𝑛−1
∑︁
𝑖=0
Define ∫ 𝑏 𝑎 𝑓(𝑥) 𝑑𝑄𝑛(𝑥) = 𝑛−1 ∑︁ 𝑖=0 𝑓(𝜉𝑖)𝑄𝑛({𝜉𝑖}). Then we have ∫ 𝑏 𝑎 𝑓(𝑥) 𝑑𝑄𝑛(𝑥) → ∫ 𝑏 𝑎 𝑓(𝑥) 𝑑𝑥, 𝑛 → ∞.
That is, if we cannot calculate the Riemann integral exactly, we can calculate it approximately:
∫ 𝑏 𝑎 𝑓(𝑥) 𝑑𝑥 ≈ ∫ 𝑏 𝑎 𝑓(𝑥) 𝑑𝑄𝑛(𝑥).
2.3.2
Second Quadrature Method
Similarly to the first quadrature method, letP( 𝑋 (1)) be the measure only this time with infinite
support, given by
P( 𝑋 (1))( 𝐴) =P{𝑋 (1) ∈ 𝐴}. Let 𝑓 : ℝ𝑑
→ ℝ be the discounted payoff of a financial derivative. For example, 𝑑 = 1, and
𝑓(𝑥) = 𝑒−𝑟𝑇max{𝑥 − 𝐾, 0},
the discounted payoff of a European call option. The price of the option is
E[ 𝑓 ( 𝑋 (1))] =
∫ ∞
−∞
𝑓(𝑥) 𝑑P( 𝑋 (1))(𝑥). (2.6)
How to find a measure 𝑄𝑛in this case? In other words, we seek for such a 𝑄𝑛that
E[ 𝑓 ( 𝑋 (1))] ≈ ∫ ∞ −∞ 𝑓(𝑥) 𝑑𝑄𝑛(𝑥) = 𝑛−1 ∑︁ 𝑖=0 𝑓(𝑦𝑖)𝑄𝑛({𝑦𝑖}).
In our case price of the European call option is modified equation (2.6). We substitute 𝑓 in
equation (2.6) with equation (2.4) 𝑓 (𝑥) = 𝑒−𝑟𝑇max{𝑥 − 𝐾, 0} and now we have derived exact
formula for calculating the expected value of the European call option.
E[ 𝑓 ( 𝑋 (1))] =
∫ ∞
−∞
𝑒−𝑟𝑇max{𝑥 − 𝐾, 0} 𝑑P( 𝑋 (1))(𝑥).
2.4
The Description of the Algorithm
The following section are stated in Müller-Gronbach and Yaroslavtseva (2016) in the case of description of the algorithm and of course Kloeden and Platen (1992) book for the weak Taylor scheme.
Our goal is to find the probability measureQ𝑛such that:
Q𝑛( 𝑋0, 𝑎, 𝑏) = 𝑅𝜏𝑛◦ 𝐷𝜏𝑛 ◦ 𝑇
𝑎,𝑏
𝜏𝑛 ◦ ... ◦ 𝑅𝜏1◦ 𝐷𝜏1 ◦ 𝑇
𝑎,𝑏 𝜏1 (𝛿𝑋0)
Let M1(ℝ𝑑) be the set of all probabilistic measures on ℝ𝑑 with finite support. We start from
the measure 𝛿𝑥0 ∈ M1(ℝ
𝑑
) such that 𝑆 = {𝑥0}, 𝛿𝑥0(∅) = 0, 𝛿𝑥0(𝑆) = 1, and implement a
non-uniform discretization 0 = 𝑡0 < ... < 𝑡𝑛=1 with a step-size 𝜏𝑖 = 𝑡𝑖− 𝑡𝑖−1 where 𝑡𝑖 =1 − (1 − 𝑖/𝑛) 𝜂 , 0 ≤ 𝑖 ≤ 𝑛 (2.7) for 𝜂 > 0.
A single step of length 𝜏 can be explained by applying the Ito Taylor transition by the following map 𝑇𝑎,𝑏 𝜏1 : M1(ℝ 𝑑 ) → M1(ℝ𝑑 )
To the measure 𝛿𝑥0, and then we apply the map
𝐷𝜏
1: M1(ℝ
𝑑
) → M1(ℝ𝑑
)
which is the support diameter reduction to the result of the previous map and is a probabilistic measure with finite support. Finally, we apply the map
𝑅𝜏
1: M1(ℝ
𝑑
) → M1(ℝ𝑑
)
to the previous result to reduce the size of the support of a discrete measure.
Define a 1-dimensional case of the weak order 𝛾 Itô-Taylor scheme which obtained by adding some parts of the so-called Itô–Taylor expansion to the Euler scheme
ˆ 𝑋𝑦,𝑎,𝑏(𝜏) = 𝑦 + ∑︁ 𝛼∈Γ𝛾 𝜓𝛼(𝑦)𝐽𝛼(𝜏) (2.8) where 𝑦 ∈ ℝ𝑑
is the initial value. 𝐽𝛼(𝜏) define as a ℓ-dimensional Itô-integral
𝐽𝛼(𝜏) = ∫ 𝜏 0 ∫ 𝑢ℓ 0 · · · ∫ 𝑢3 0 ∫ 𝑢2 0 1 d𝑊𝛼1(𝑢1) d𝑊𝛼2(𝑢2) . . . d𝑊𝛼ℓ−1(𝑢ℓ−1) d𝑊𝛼ℓ(𝑢ℓ), (2.9) where 𝜏 ∈ [0, 1], 𝑊0(𝑢) = 𝑢 for 𝑢 ≥ 0 and for ℓ ∈ ℕ, 𝛼 = (𝛼1, ..., 𝛼ℓ) ∈ ℕℓ
0.
The symbol Γ𝛾 denotes the set of multi-indices with 𝛼 = (𝑚 + 1)
𝛾
+ 𝑚 + 1 elements (0), (1), . . . , (𝑚), (0, 0), (0, 1), . . . , (0, 𝑚), (1, 0), . . . , (1, 𝑚), . . . , (𝑚, 0), . . . , (𝑚, 𝑚), such that
𝑚is the number of the Wiener process.
Example 1. How many elements is contained in Γ2when 𝑚 = 1? Write down all of them.
the set Γ2contains 6 elements:
Define the Kloeden–Platen family { 𝜉𝛼: 𝛼 ∈ Γ2} of random variables which is called a
simplification of order 𝛾 if ℙ𝜉𝛼 ∈ M1(ℝ) for all ,𝛼 ∈ Γ𝛾 and
𝔼 𝐿 Ö ℓ=1 𝜉𝛼 (ℓ) = 𝔼 𝐿 Ö ℓ=1 𝐽𝛼 (ℓ)
for all 𝛼(1), ..., 𝛼( 𝐿) ∈ Γ𝛾 with deg(𝛼(1) + ... + 𝛼( 𝐿)) ≤ 2𝛾 + 1 and 𝐿 = 1, ..., 2𝛾 + 1.
Definition 3. The degree of an index 𝛼, (deg(𝛼)), is the number of elements in the index 𝛼 plus the number of zeroes among them.
So we have,
deg(𝛼) = (
2, if 𝛼 = 0,
1, otherwise
and deg(𝛼1, 𝛼2) is the number of zeros in the pair (𝛼1, 𝛼2) plus 2.
Then we have a simplified step by defining a discrete variable 𝜉𝛼(𝜏) = 𝜏
deg(𝛼)/2𝜉 𝛼 and
replacing it by the iterated Itô-integrals 𝐽𝛼(𝜏) in equation (2.8). So we have,
𝑌𝑦,𝑎,𝑏(𝜏) = 𝑦 + ∑︁
𝛼∈Γ𝛾
𝜓𝛼(𝑦)𝜉𝛼(𝜏). (2.10)
According to the Kloeden and Platen (1992) when 𝛾 = 2 and 𝑗 , 𝑗1, 𝑗2∈ {1, ..., 𝑚} we can set
𝜉𝑗 = 𝑁𝑗, 𝜉( 𝑗
1, 𝑗2) =
1
2(𝑁𝑗1𝑁𝑗2+ 𝑉𝑗1, 𝑗2)
we can take independent 𝑁 (0, Δ) Gaussian random variables 𝑁𝑗 for 1 ≤ 𝑗 ≤ 𝑚 with three
point distribution such that:
𝑃(𝑁𝑗 = ± √ 3Δ) = 1 6, 𝑃(𝑁𝑗 =0) = 2 3
as well as 𝜉( 𝑗 ) = 𝜉(0, 𝑗 ) = 𝜉( 𝑗 ,0) = 𝑁𝑗 also we can take another independent random variable
(𝑉𝑗1, 𝑗2)1≤ 𝑗1< 𝑗2≤𝑚 with two-point distribution 𝑃((𝑉𝑗1, 𝑗2) = ±Δ) = 1
2. If 𝑗1 > 𝑗2 then we have
𝑉𝑗1, 𝑗2 = −𝑉𝑗1, 𝑗2 and when 𝑗1 = 𝑗2 and 𝑚 = 1 the 𝑉𝑗1, 𝑗2 are not required. It is summarized as
follow: 𝜉( 𝑗 1, 𝑗2) = 1 2(𝑁𝑗1𝑁𝑗2+ 𝑉𝑗1, 𝑗2), if 𝑗1 < 𝑗2, 1 2(𝑁𝑗1𝑁𝑗2− 𝑉𝑗1, 𝑗2), if 𝑗1 > 𝑗2, 1 2(𝑁2𝑗1− 1), if 𝑗1 = 𝑗2,
The Itô–coefficients function is stated as 𝜓(𝛼)(𝑦) = 𝐿𝛼
1(𝑦1) ◦ ... ◦ 𝐿𝛼ℓ(𝑦𝑑). (2.11)
In our case when 𝛾 = 2 and 𝑚 = 𝑑 = 1 we have 𝜓( 𝑗 )(𝑦) = 𝐿𝑗(𝑦), 𝜓( 𝑗
where 𝐿𝑗 is the operator 𝐿𝑗 = 𝑑 ∑︁ 𝑘=1 𝑏𝑘 , 𝑗 𝜕 𝜕 𝑦𝑘 for 𝑗 = 1, ..., 𝑚 and 𝐿𝑗 = 𝜕 𝜕 𝑡 + 𝑑 ∑︁ 𝑘=1 𝑎𝑘 𝜕 𝜕 𝑦𝑘 +1 2 𝑑 ∑︁ 𝑘 ,𝑙=1 𝑚 ∑︁ 𝑗=1 𝑏𝑘 , 𝑗𝑏𝑙 , 𝑗 𝜕2 𝜕 𝑦𝑘𝜕 𝑦𝑙 for 𝑗 = 0.
The map 𝐷𝜏 is given by
𝐷𝜏( 𝜇) =
∑︁
𝑦∈supp(𝜇)
𝜇(𝑦)𝛿b𝑦c𝐻𝜏
for 𝜇 ∈ M1(ℝ𝑑). We project points 𝑦 ∈ ℝ𝑑on to 𝐻𝜏 = [−𝜅𝜏
−𝜀 , 𝜅 𝜏−𝜀]𝑑for 𝜀, 𝑘 > 0 as follows: b𝑦c𝐻𝜏 = (𝑦𝑖+ (−𝑦𝑖− 𝜅𝜏 −𝜀) +− (𝑦𝑖− 𝜅𝜏 −𝜀) +)𝑖=1,...,𝑑 Put 𝐴𝑗 ,𝜏= [ 𝑗 √ 𝜏,( 𝑗 + 1) √
𝜏) for 𝑗 ∈ ℤ to define finite set of such 𝑗 denoted by 𝐽𝜇, 𝜏which
intersect with interval 𝐻𝜏 or in other word with the support of 𝐷𝜏◦ 𝑇
𝑎,𝑏 𝜏
𝐽𝜇,𝜏 = { 𝐴𝑗 ,𝜏∩ supp(𝜇) ≠ ∅}.
We implement the support cardinality reduction strategy by using the target cardinality bound 𝑘 = 2𝛾+𝑑+1 𝑑 . Choose 𝑢 ∈ ℝ𝑘+1 \ {0} that satisfies 𝑘+1 ∑︁ 𝑖=1 𝑎𝑗(𝑥𝑖)𝑢𝑖 =0, 1 ≤ 𝑗 ≤ 𝑘 . where 𝑎𝑗(𝑥𝑖) = 𝑥 𝑗−1
𝑖 are the entries of a (2𝛾 + 2) × (2𝛾 + 3) Vandermonde matrix 𝐴 and then
solve the system of linear equations 𝐴𝑢 = 0.
Assume that the measure to which we apply 𝑅𝜏, has the form
𝜇= 𝑁 ∑︁ 𝑖=1 𝑤𝑖𝛿𝑥 𝑖. Set 𝛼 = min𝑖:𝑢𝑖<0 −𝑤𝑖 𝑢𝑖
where {𝑖 : 𝑢𝑖 < 0 ≠ ∅}. Keep all weights 𝑤𝑖, if 𝑖 > 𝑘 + 1 and define
a new weight𝑤e𝑖 = 𝑤𝑖+ 𝛼.𝑢𝑖 for 1 < 𝑖 < 𝑘 + 1 which leads to a new measure,
e𝑢= 𝑁 ∑︁ 𝑖=1 e 𝑤𝑖.𝛿𝑥 𝑖
as a result𝑤e𝑖 > 0 and we have 𝑤𝑖 =0 for some 1 < 𝑖 < 𝑘 + 1 in which the minimal possible
value of 𝛼 attains and leads to the elimination of the corresponding point 𝑥𝑖 from the support.
We repeat this procedure 𝑁 − 𝑘 times.
Finally, we describe some theoretical results about the convergence of the above algorithm, see Müller-Gronbach and Yaroslavtseva (2016).
Chapter 3
Application to the Black-Scholes model
In this chapter the Black–Scholes model will be implemented to the deterministic algorithm that explained in chapter 2.
3.1
The Non-uniform Time Discretization
The non-uniform discretization is implemented by using equation 2.7 which described in previous chapter where we set 𝜂 = 2 and 𝑛 = 128.
3.2
The Itô–Taylor Transition
Consider the Black–Scholes stochastic differential equation previously stated in (2.3), compare it with the system of stochastic differential equation (2.1) . The equivalent equation to (2.3) will turn out to the following equation with 𝑇 = 1.
d𝑆(𝑡) = 𝑟𝑇 𝑆(𝑡) d𝑡 + 𝜎 √ 𝑇 𝑆(𝑡) d𝑊 (𝑡), 𝑡 ∈ [0, 1], 𝑆(0) = 𝑆0. we obtain 𝑎(𝑦) = 𝑟𝑇 𝑦, 𝑏(𝑦) = 𝜎 √ 𝑇 𝑦, where 𝑦 = 𝑆(𝑡).
For calculations, we intend to use a particular case of Equation (2.8) with 𝑑 = 𝑚 = 1 and
Γ = 2. The set Γ2 has been explicitly written in Example 1. Using the results of the above
example, we obtain: ˆ
𝑋𝑦,𝑎,𝑏(𝜏) = 𝑦 + 𝜓0(𝑦)𝐽0(𝜏) + 𝜓1(𝑦)𝐽1(𝜏) + 𝜓(0,0)(𝑦)𝐽(0,0)(𝜏)
+ 𝜓(0,1)(𝑦)𝐽(0,1)(𝜏) + 𝜓(1,0)(𝑦)𝐽(1,0)(𝜏) + 𝜓(1,1)(𝑦)𝐽(1,1)(𝜏). (3.1)
The Itô–coefficient functions 𝜓𝛼(𝑦) can be calculated by Equation (2.11), which involves a lot
of work. Instead, we use the fact mentioned in Müller-Gronbach and Yaroslavtseva (2016): Equation (3.1)
. . . is a single step of size 𝜏 with initial value 𝑦 of the weak order 𝛾 Itô–Taylor scheme. . .
The above scheme is given in (Kloeden and Platen, 1992, Section 14.2, Equation (2.1)). In our notation, this equation takes the form
𝑌𝑛+1= 𝑌𝑛+ 𝑎𝜏 + 𝑏𝑊 (𝜏) + 1 2𝑏 𝑏 0 (𝑊2(𝜏) − 𝜏) + 𝑎0𝑏 ∫ 𝜏 0 𝑊(𝑢2) d𝑢2+ 1 2 𝑎 𝑎0+ 1 2𝑎 00 𝑏2 𝜏2 + 𝑎 𝑏0+ 1 2𝑏 00 𝑏2 𝜏𝑊(𝜏) − ∫ 𝜏 0 𝑊(𝑢2) d𝑢2 . (3.2)
We compare Equations (3.1) and (3.2) as follows. For each multi-index 𝛼 ∈ Γ2we calculate
the multiple Itô integral 𝐽𝛼(𝜏), using Equation (2.9). The resulting integral is a part of exactly
one term in the right hand side of Equation (3.2). The remaining part of this term is equal to 𝜓𝛼(𝑦).
Put 𝛾 = 0. Equation (2.9) gives
𝐽0(𝜏) =
∫ 𝜏
0
d𝑢1= 𝜏.
The unique term in the right hand side of Equation (3.2) that contains 𝜏, is 𝑎𝜏. It follows that 𝜓0(𝑦) = 𝑎 = 𝑟𝑇 𝑦. Next, 𝐽1(𝜏) = ∫ 𝜏 0 d𝑊 (𝑢1) = 𝑊 (𝜏).
The corresponding term is 𝑏𝑊 (𝜏), and we obtain 𝜓1(𝑦) = 𝑏 = 𝜎
√ 𝑇 𝑦.
Next term is,
𝐽(0,0)(𝜏) = ∫ 𝜏 0 ∫ 𝑢2 0 d𝑢1d𝑢2 = ∫ 𝜏 0 𝑢2d𝑢2 = 1 2𝜏 2.
The corresponding term of equation (3.2) is 12 𝑎 𝑎0+ 1
2𝑎 00 𝑏2 𝜏2, and we obtain 𝜓(0,0)(𝑦) = 𝑎𝑎0+ 1 2𝑎 00 𝑏2= 𝑟𝑇 𝑦 · 𝑟𝑇 = 𝑟2𝑇2𝑦. Then , 𝐽(1,0)(𝜏) = ∫ 𝜏 0 ∫ 𝑢2 0 d𝑊 (𝑢1) d𝑢2= ∫ 𝜏 0 𝑊(𝑢2) d𝑢2. By equation (3.2), we obtain 𝜓(1,0)(𝑦) = 𝑎0𝑏 = 𝑟𝑇 𝜎 √ 𝑇 𝑦 = 𝑟 𝜎𝑇3/2𝑦.
Next, 𝐽(1,1)(𝜏) = ∫ 𝜏 0 ∫ 𝑢2 0 d𝑊 (𝑢1) d𝑊 (𝑢2) = ∫ 𝜏 0 𝑊(𝑢2) d𝑊 (𝑢2) = 1 2(𝑊 2(𝜏) − 𝜏), By equation (3.2), 𝜓(1,1)(𝑦) = 𝑏𝑏0= 𝜎 √ 𝑇 𝑦𝜎 √ 𝑇 = 𝜎2𝑇 𝑦. Finally 𝐽(0,1)(𝜏) = ∫ 𝜏 0 ∫ 𝑢2 0 d𝑢1d𝑊 (𝑢2) = ∫ 𝜏 0 𝑢2d𝑊 (𝑢2) = 𝜏𝑊 (𝜏) − ∫ 𝜏 0 𝑊(𝑢2) d𝑢2. We recognize Δ𝑊Δ − Δ𝑍 and obtain
𝜓(0,1)(𝑦) = 𝑎𝑏0+ 1 2𝑏 00 𝑏2 = 𝑟𝑇 𝑦𝜎 √ 𝑇 = 𝑟𝑇3/2𝜎𝑦.
Now Equation (3.1) becomes ˆ
𝑋𝑦,𝑎,𝑏(𝜏) = 𝑦 + 𝑟𝑇 𝑦𝐽0(𝜏) + 𝜎 √
𝑇 𝑦 𝐽1(𝜏) + 𝑟2𝑇2𝑦 𝐽(0,0)(𝜏)
+ 𝑟 𝜎𝑇3/2𝑦 𝐽(0,1)(𝜏) + 𝑟 𝜎𝑇3/2𝑦 𝐽(1,0)(𝜏) + 𝜎2𝑇 𝑦 𝐽(1,1)(𝜏).
We replace 𝐽𝛼(𝜏) with 𝜉𝛼(𝜏), where
𝜉0(𝜏) = 𝜏, 𝜉(0,0)(𝜏) = 1 2𝜏 2, 𝜉 1(𝜏) = 𝑁1 √ 𝜏, 𝜉(1,0)(𝜏) = 𝑁1𝜏3/2, 𝜉(0,1)(𝜏) = 𝑁1𝜏3/2, 𝜉(1,1)(𝜏) = 1 2(𝑁 2 1 − 1)𝜏,
Example 2. Write down the Kloeden–Platen family for 𝑑 = 𝑚 = 1 in the case 𝛾 = 2 .
𝜉 = (𝜉(0), 𝜉(1), 𝜉(0,0), 𝜉(0,1), 𝜉(1,0), 𝜉(1,1)) a simplification 𝜉 of order 2 by 𝜉( 𝑗 ) = 𝜉(0) = 𝜉(1) = 𝜉(2) = 𝜉(3) =1.
Similarly,
𝜉( 𝑗 ) = 𝜉(0, 𝑗 ) = 𝜉( 𝑗 ,0) which yields 𝜉(0,1) = 𝜉(1,0) = 1. and 𝜉( 𝑗1, 𝑗2) =
1 2(𝑁𝑗1𝑁𝑗2+ 𝑉𝑗1, 𝑗2) such that, (𝑉𝑗1, 𝑗2)1≤ 𝑗1< 𝑗2≤𝑚, so 𝜉(0,0) =1/2. Furthermore, 𝜉( 𝑗1, 𝑗2) = (𝑁 2 𝑗1− 1)/2, If 𝑗1= 𝑗2So we have 𝜉(1,1) =1/2.
Then we have simplified weak Itô–Taylor
𝑌𝑛+1 = 𝑌𝑛+ 𝑟𝑇 𝑦𝜏 + 𝜎 √ 𝑇 𝑦 𝑁1 √ 𝜏+ 𝑟2𝑇2𝑦1 2𝜏 2 + 𝑟 𝜎𝑇3/2𝑦 𝑁1𝜏3/2+ 𝑟 𝜎𝑇3/2𝑦 𝑁1𝜏3/2+ 𝜎2𝑇 𝑦 1 2(𝑁 2 1 − 1)𝜏
andP{𝑁1 = − √
3} = P{𝑁1 =
√
3} = 16, P{𝑁1 =0} = 23. So it becomes three–point set with elements 1 6[𝑦 + 𝑟𝑇 (𝑦)𝜏 + 𝑟 2 𝑇2(𝑦) 1 2𝜏 2+ (√ 3𝜎 √ 𝑇(𝑦) √ 𝜏) + (2 √ 3𝑟𝜎𝑇32(𝑦)𝜏 3 2) + 𝜎2𝑇(𝑦)𝜏] 1 6[𝑦 + 𝑟𝑇 (𝑦)𝜏 + 𝑟 2𝑇2(𝑦)1 2𝜏 2− (√3𝜎√𝑇(𝑦)√𝜏) − (2√3𝑟𝜎𝑇32(𝑦)𝜏32) + 𝜎2𝑇(𝑦)𝜏] 2 3[𝑦 + 𝑟𝑇 (𝑦)𝜏 + 𝑟 2𝑇2(𝑦)1 2𝜏 2− 1 2𝜎 2𝑇(𝑦)𝜏]
Clearly the number of points in the support of the measure would grow very fast after N steps. By this reason, after each simplified weak Itô–Taylor step, we perform additional steps to reduce the above number.
3.3
The Support Diameter Reduction Method
Based on the algorithm all points of the support move to the fixed box called 𝐻𝜏 = [−𝜅𝜏
−𝜀
, 𝜅 𝜏−𝜀]𝑑 and since in our case 𝑑 = 1 it is just an interval.
The points that belong to 𝐻𝜏 are kept as they are and the points that lies to the left of 𝐻𝜏
are moving to its left end and their weights are added and the same for the points that lie to the
right of 𝐻𝜏 which are moving to its right end and their weights are added. This leads to the
reduction of the number of points in support as the first additional step.
3.4
The Support Cardinality Reduction Method
At the second additional step, we cover the space ℝ1 with the intervals [ 𝑗√𝜏,( 𝑗 + 1)
√ 𝜏) of
length √𝜏. Only finitely many of them intersect with the support of 𝐷𝜏 ◦ 𝑇𝑎,𝑏
𝜏 . Denote by
𝐽𝜇,𝜏 the finite set of such 𝑗 , that the intervals [ 𝑗 √
𝜏,( 𝑗 + 1) √
𝜏) and 𝐻𝜏 intersect and then the support partitioned into the number of disjoint subsets. Put the target cardinality bound
𝑘 = 2𝛾 + 𝑑 + 1 𝑑 = 6 1 =6
in order to apply the sequential support point elimination method. If the number of points in the subset is less than 𝑘 = 6, we do nothing. Otherwise, let 𝑁 > 𝑘 be the number of points in the subset, let 𝑥𝑖be the points in the support, and let 𝑤𝑖be their weights, that is,P{𝜉 = 𝑖} = 𝑤𝑖.
A system of linear equations is given by
𝑘+1 ∑︁ 𝑖=1 𝑎𝑗(𝑥𝑖)𝑢𝑖 =0, 1 ≤ 𝑗 ≤ 𝑘 . where 𝑎𝑗(𝑥𝑖) = 𝑥 𝑗−1
𝑖 .This system contains 𝑘 = 6 equations and 𝑘 + 1 = 7 unknowns 𝑢1, . . . ,
𝑢7. The solutions constitute a line in the space ℝ7. Choose a solution which has at least one
negative component. Define
where we take into account only those ratios −𝑤𝑖/𝑢𝑖 which are positive. Keep all weights 𝑤𝑖,
if 𝑖 > 𝑘 + 1, and replace
𝑤𝑖 ↦→ 𝑤𝑖+ 𝛼𝑢𝑖, 1 ≤ 𝑖 ≤ 𝑘 + 1.
Exactly one new weight will become equal to 0. Its number 𝑖 is exactly the number in which the
minimal possible value of 𝛼 attains. The corresponding point 𝑥𝑖 disappears from the support.
We repeat this procedure 𝑁 − 𝑘 times.
3.5
Realization
In this section, the MATLAB simulation script is explained more in detail. Full MATLAB simulation script can be found at the end of this paper in appendix A. Bold numbers refer to code lines in the main simulation script.
In prior, parameters of the Black–Scholes have been defined: 1. S — The initial price of the underlying stock;
2. K — The exercise price; 3. R — The risk-free interest rate; 4. T — The maturity of the option;
5. Sigma — The volatility of the stock per annum.
Then the parameters 𝜂, 𝜅, 𝜖 used in the deterministic algorithm are defined and N is the number of intervals, we divide interval from 0 to 1 into 𝑁 sub-intervals.
Black-Scholes call option price using MATLAB build in command blsprice is used and the price is displayed in order to compare it with our final result.
Then the parameters of the Black–Scholes SDE recalculated to make 𝑇 = 1 with new values for 𝑅 and 𝜎. The non-uniform time discretization is calculated based on the formula in the description of the Algorithm part of the thesis and the step size 𝜏 is defined.
Within lines 51–53 two working matrices are assigned with two rows and many columns. First row is for coordinates of the points. Second row contains the weights of each point. S is initial price of the security by using it directly in the matrices we obtain much better result. We assign sum of all weights to be equal to 1.
We introduce some necessary constants from lines 55–59 that will be used throughout the loop. It is better to calculate them once before loop starts and then just recall them when necessary from the loop. The loop is started with the number of intervals 𝑁 then some constants are calculated depend on number of 𝜏. Between lines 70–80 the first step of our algorithm is realized. We have several points and now instead of each individual point, we create three new points and create three new weights to this new points.
New matrix measurenew is created, we increase number of points by three times. So, now our matrix has three times more columns. New coordinates and weights are calculated and are assigned into new matrix measurenew.
In the loop between lines 84–91 the Itô-Taylor transition is implemented by using the previous constants defined before. Since we would like to have our data of the points in ascending order.these following lines 93-97 sort the array into the right order. Now we copy new results from matrix measurenew into matrix measure. This way we continue using measure as our main matrix of work.
Now we enter into second part of our algorithm called The support diameter reduction.In line 102 one half of our interval 𝐻𝜏 or the left side of it is calculated and called edge. Since it
is symmetric with respect to the point S it is sufficient to calculate it once.
Within following lines 103–148 move the points which lie outside the interval or the edge and the weights of the points that lie outside left or right side of the interval are added. Then we again transfer matrix measurenew results into matrix measure for further calculations.
In the lines 152–221 we realize third part of the algorithm which called the support
cardinality reduction. In this part, we have all points within interval. This interval is devided
into sub-intervals of equal lengths. The length of each sub interval is√𝜏.
In the line 154 the number of intervals is calculated into which we divide our interval. Then we reserve place for twice as many intervals we have in line 155 because we have left and right side from point S.
Recall from the target cardinality bound which calculated and is equal to 6. If the number of points in interval is less than seven then we just copy these points into measurenew. Otherwise, we have to decrease the number of points in the interval.
In lines 181–216 we work on intervals that contain more than six points and eliminate points until we achieve six points per interval, then the six points are copied into the array
measurenew.
The big loop (62–225) ends after being repeated 𝑁 times.Lines 227–228 illustrate our quadrature formula for calculating European call option price. We take the value of the option. We take measure compare with strike price 𝐾 and choose maximal. Then we multiply by weight and add everything together using the dot product function in MATLAB.
Finally, the discounted payoff of the financial derivative is obtained by using the determin-istic algorithm and compared to MATLAB Black–Scholes function.
Consider the system (2.1) and assume the following assumptions.
Assumption 1. There exists a real number 𝐾 ∈ [1, ∞) such that 𝑥0 ∈ [−𝐾, 𝐾]𝑑, the
compon-ents of the functions 𝑎 and 𝑏 are 6 times continuously differentiable, and the absolute values of all their partial derivatives up to order 6 are bounded by 𝐾 from above.
Assumption 2. There exist a positive integer 𝑟, a nonnegative real number 𝛽, and a positive constant 𝐶 such that the absolute values of all partial derivatives of the function 𝑓 (𝑥) up to
order 𝑟 are bounded by 𝐶 (1 + |𝑥|𝛽
) from above.
Assumption 3. There exists a real number 𝜆 ∈ (0, 1] such that the minimal eigenvalue of the
symmetric positive-definite matrix 𝑏(𝑥)𝑏(𝑥)> is bounded by 𝜆 from below.
Under these assumptions, Müller-Gronbach and Yaroslavtseva (2016) proved the following result. Denote 𝑞 = min{6, 𝑟 }, let 𝑛 be the number of subintervals in the above algorithm, and let 𝑒 be the absolute error of the approximation.
Theorem 1. Let the above assumptions hold true. Then, there exists a positive constant 𝑐 that depends on 𝑓 , such that
• if 𝜂 < 4 𝑞, then 𝑒 ≤ 𝑐 𝑛𝜂; • if 𝜂 = 4 𝑞, then 𝑒 ≤ 𝑐ln(𝑛) 𝑛𝜂 ; • if 𝜂 > 4 𝑞, then 𝑒 ≤ 𝑐 𝑛2.
Note that the payoff function 𝑓 (𝑥) = max{𝑥 − 𝐾, 0} of a European call option does not satisfy Assumption 2. However, the algorithm still can give a good enough approximation for some initial values. We investigate this question in the next chapter.
Chapter 4
Numerical Results
In this chapter, we examine the accuracy of the deterministic quadrature formulae originally proposed by Müller-Gronbach and Yaroslavtseva (2016). By using the deterministic quadrature formulae and the MATLAB simulation. We calculate the European call option price. Next, we compare the deterministic quadrature formulas call option price with the price derived from using the Black—Scholes model. We aim to see how close is the deterministic quadrature formulae result to the absolute value of the European call option price. All the obtained numerical results from MATLAB simulation are presented in graphs and tables. The chapter aims to see how close to real value the deterministic quadrature formulae can yield with the initial and parameter values given in Table 4.1.
The main constants parameters of the Black–Scholes shown in Table 4.2.
Using the MATLAB build in command for the Black-Scholes model blsprice we get the European call option price to be 𝑐 = 9.339017532400462. After running MATLAB simulation a few times using the deterministic quadrature formulae yields the European call option price to be 𝑐 = 9.338303563702935. As we can see from the numerical results, the deterministic quadrature formulae yield a very accurate price for the European call option. In the following sections, we will experiment with different parameters of the option to see if any changes to the accuracy will accrue. The main result of this chapter is numerical simulation under different conditions. 𝑆 100.0 𝐾 95.0 𝑅 0.1 𝑇 0.25 𝜎 0.25
𝜂 2.0
𝜅 33.0
𝜖 0.1
𝑁 128
Table 4.2: The parameters of the method
4.1
Comparison of the Initial Price of the Underlying Stock
In this section, alteration of the initial price of the underlying stock will be considered to illustrate its effects on the deterministic quadrature formulae. We have changed our MATLAB simulation script a little bit to accommodate our desire to perform changes in the initial price of the underlying stock. New MATLAB simulation performs the range of the initial price of the underlying stock from 50 to 120 with an increment of 1 each time. The Graph illustrates the accuracy of the employed algorithm with respect to the MATLAB build in command for the Black–Scholes model.
Figure 4.1: Comparison of the initial price of the underlying stock
From the graph 4.1, we observe that change in the initial price of the underlying stock does not influence on the accuracy of the deterministic quadrature formulae.
4.2
Comparison of the Exercise Price
In this section, variation of the exercise price of the underlying stock will be considered to illustrate its effects on the accuracy of the deterministic quadrature formula. New MATLAB simulation performs the range of the exercise price from 50 to 120 with an increment of 1 each time. The Graph illustrates the accuracy of the employed algorithm with respect to the MATLAB build in command for the Black–Scholes model.
From graph 4.2, we observe that change in the exercise price does not influence the accuracy of the deterministic quadrature formulae.
Figure 4.2: Comparison of the exercise price
4.3
Comparison of the Risk-Free Interest Rate
In this section, we will look closer at how a change in the risk-free interest rate the affects accuracy of the deterministic quadrature formulae. New MATLAB simulation performs for the range of the risk-free interest rate from 0.1 to 1 with an increment of 0.05 each time. The Graph illustrates the accuracy of the employed algorithm with respect to the MATLAB build in command for the Black–Scholes model.
Figure 4.3: Comparison of the risk-free interest rate
From graph 4.3, we observe that change in the risk-free interest rate starts affecting the deterministic quadrature formulae numerical results. Deviation starts to occur when the risk-free interest rate rises are above 60%.
4.4
Comparison of the Maturity of the Option
In this section, we will look closer at how a change in the option’s maturity affects the accuracy of the deterministic quadrature formulae. The new MATLAB simulation performs for the range of the option’s maturity from 0.1 to 1 with an increment of 0.05 each time. The Graph illustrates the accuracy of the employed algorithm with respect to the MATLAB build in command for the Black–Scholes model.
Figure 4.4: Comparison of the maturity of the option
From graph 4.4, we observe that the change in the option’s maturity starts affecting the
deterministic quadrature formulae numerical results. Deviation starts to occur when the
option’s maturity is between 0.4–0.45, and variation only increases with the additional time to maturity.
4.5
Comparison of the Volatility of the Stock per Annum
In this section, we will look closer at how a change in the volatility of the stock per annum affects the accuracy of the deterministic quadrature formulae. New MATLAB simulation performs a range of volatility from 0.1 to 0.5 with an increment of 0.01 each time. The Graph illustrates the accuracy of the employed algorithm with respect to the MATLAB build in command for the Black-Scholes model.
From graph 4.5, we observe that the change in the volatility of the stock per annum starts affecting the deterministic quadrature formulae numerical results. Deviation starts to occur when the volatility of the option has reached 0.3, and variation only increases with the increment of volatility.
4.6
Summary
To summarise the numerical results, we can conclude that the deterministic quadrature formulae accuracy is good with the initial and parameter values. We have studied five different parameters
Figure 4.5: Comparison of the volatility of the stock per annum
to check the accuracy. Change in the initial price and the exercise price of the underlying stock does not affect the accuracy of the deterministic quadrature formulae. Change in the risk-free interest rate has a minimal effect on the accuracy of the deterministic quadrature formulae. The first significant variation in the accuracy of the deterministic quadrature formulae accrues when the option’s maturity surpasses 0.4-0.45, and variation only increases with the additional time to maturity. The second significant variation in the accuracy of the deterministic quadrature formulae accrues when the volatility of the stock per annum exceeds 0.3, and variation only increases with the increment of volatility.
Chapter 5
Conclusions
5.1
Conclusions
In this project, the main idea is estimate the expected value of a function 𝑓 of a numerical solution of system of stochastic differential equation by employing deterministic quadrature formula proposed by Müller-Gronbach and Yaroslavtseva (2016) in Black-Scholes model.
The estimated value which obtain by the deterministic algorithm has compared with the analytical European call option which yield a value close to the exact call option price. However, in numerical result section of the thesis we examine the accuracy of the method by changing initial values of the option. According to our result, we observe that the algorithm is not sensitive to changes in initial price, exercise price and risk-free interest rate, whereas it deviate in case when the maturity is more than approximately 0.45 as well as volatility when it reached 0.3.
5.2
Further Research
The deterministic method can implement in different market model. In this thesis we apply it for one-dimensional case ,while ones can apply it for multi-dimensional case for example Gatheral model where 𝑑 = 3.
Bibliography
D. Henderson and P. Plaschko. Stochastic differential equations in science and engineering. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2006. ISBN 981-256-296-6. doi: 10.1142/9789812774798. URL https://doi.org/10.1142/9789812774798. With 1 CD-ROM (Windows, Macintosh).
J. C. Hull. Options, futures and other derivatives. Pearson, 2018.
I. Karatzas and S. Shreve. Brownian Motion and Stochastic Calculus, volume 113 of Graduate
Texts in Mathematics. Springer, New York, second edition, 1998. ISBN 978-0-387-97655-6.
P. E. Kloeden and E. Platen. Numerical solution of stochastic differential
equa-tions, volume 23 of Applications of Mathematics (New York). Springer-Verlag,
Berlin, 1992. ISBN 3-540-54062-8. doi: 10.1007/978-3-662-12616-5. URL
https://doi.org/10.1007/978-3-662-12616-5.
P. E. Kloeden, E. Platen, and H. Schurz. Numerical solution of SDE through computer
experiments. Universitext. Springer-Verlag, Berlin, 1994. ISBN 3-540-57074-8. doi: 10.1007/978-3-642-57913-4. URL https://doi.org/10.1007/978-3-642-57913-4. With 1 IBM-PC floppy disk (3.5 inch; HD).
E. Lindström, H. Madsen, and J. N. Nielsen. Statistics for finance. Chapman & Hall/CRC Texts in Statistical Science Series. CRC Press, Boca Raton, FL, 2015. ISBN 978-1-4822-2899-1.
T. Müller-Gronbach and L. Yaroslavtseva. Deterministic quadrature formulas for
SDEs based on simplified weak Itô–Taylor steps. Found. Comput. Math., 16(5):
1325–1366, 2016. ISSN 1615-3375. doi: 10.1007/s10208-015-9277-5. URL
https://doi.org/10.1007/s10208-015-9277-5.
D. Nualart. The Malliavin calculus and related topics. Probability and its Applications (New York). Springer-Verlag, Berlin, second edition, 2006. ISBN 978-28328-7; 3-540-28328-5.
B. Øksendal. Stochastic differential equations. An introduction with applications. Universitext. Springer-Verlag, Berlin, sixth edition, 2003. ISBN 3-540-04758-1. doi: 10.1007/978-3-642-14394-6. URL https://doi.org/10.1007/978-3-10.1007/978-3-642-14394-6.
G. Pagès. Numerical probability. An introduction with applications to finance. Universitext. Springer, Cham, 2018. ISBN 978-3-319-90274-6; 978-3-319-90276-0. doi: 10.1007/978-3-319-90276-0. URL https://doi.org/10.1007/978-10.1007/978-3-319-90276-0.
S. Särkkä and A. Solin. Applied stochastic differential equations, volume 10 of Institute of
Mathematical Statistics Textbooks. Cambridge University Press, Cambridge, 2019. ISBN 978-1-316-64946-6; 978-1-316-51008-7.
Appendix A
MATLAB Simulation
1 % Script: Deterministic.m
2 % Authors: Timo Kudljakov, Sajedeh Saadat
3 % Aim: Calculation of an approximate price of the European call option 4 % in the Black-Scholes model by a deterministic quadrature formula 5 % Algorithm: Thomas Muller-Gronbach, Larisa Yaroslavtseva. ...
Deterministic
6 % Quadrature Formulas for SDEs Based on Simpified Weak Ito-Taylor ...
Steps.
7 % Found. Comput. Math 16 (2016), 1325-1366. 8 % Version: April 5, 2021
9 %
10 % Constants 11 %
12 % The parameters of the option 13 %
14 % The initial price of the underlying stock 15 S=100.0;
16 % The exercise price 17 K=95.0;
18 % The risk-free interest rate 19 R=0.1;
20 % The maturity of the option 21 T=0.25;
22 % The volatility of the stock per annum 23 SIGMA=0.25;
24 %
25 % The parameters of the method 26 % 27 ETA=2.0; 28 KAPPA=33.0; 29 EPSILON=0.10; 30 N=128; 31 %
32 % Calculate and show the Black-Scholes price 33 format long;
34 [Call,Put]=blsprice(S,K,R,T,SIGMA); 35 disp(Call);
36 % Recalculate the parameters of the option to make T=1 37 r=R*T;
38 sigma=SIGMA*sqrt(T);
39 % Calculate the non-uniform time discretization 40 t=zeros(size(1:N+1)); 41 for i=1:N+1 42 t(i)=i-1; 43 end 44 t=1-(1-t/N).^ETA; 45 tau=t(2:N+1)-t(1:N);
46 % The number of columns in the working matrices 47 M=1;
48 % the working matrix 49 measure=zeros(2,M);
50 % The first row contains the points 51 measure(1,1)=S;
52 % The second row contains the weights 53 measure(2,1)=1.0;
54 % The necessary constants 55 r2=r*r;
56 sigma2=sigma*sigma; 57 sqrt3=sqrt(3); 58 work=[0 0 0]; 59 work2=0;
60 % Create the waitbar
61 w = waitbar(0,'Please wait...'); 62 for i=1:N
63 % the first step
64 % calculate some local constants 65 tau1=tau(i); 66 tau12=sqrt(tau1); 67 tau2=tau1*tau1; 68 tau32=tau1*tau12; 69 work=[1 1 1]; 70 work=work*(1+r*tau1+r2*tau2/2); 71 work3=sqrt3*tau12*sigma; 72 work(1)=work(1)-work3; 73 work(3)=work(3)+work3; 74 work3=2*sqrt3*r*sigma*tau32; 75 work(1)=work(1)-work3; 76 work(3)=work(3)+work3; 77 work3=sigma2*tau1; 78 work(1)=work(1)+work3; 79 work(2)=work(2)-work3/2; 80 work(3)=work(3)+work3;
81 % create the matrix measurenew 82 measurenew=zeros(2,3*M);
83 % The Ito-Taylor transition 84 for kk=1:M
86 p=measure(2,kk); 87 measurenew(1,3*kk-2:3*kk)=y*work; 88 measurenew(2,3*kk-2)=p/6; 89 measurenew(2,3*kk-1)=2*p/3; 90 measurenew(2,3*kk)=measurenew(2,3*kk-2); 91 end
92 % sort the result 93 work1=measurenew(1,:); 94 work2=measurenew(2,:); 95 [B,Index]=sort(work1); 96 measurenew(1,:)=B;
97 measurenew(2,:)=work2(Index); 98 % copy the result back
99 measure=measurenew; 100 M=3*M;
101 % The support diameter reduction 102 edge=KAPPA*tau1^(-EPSILON);
103 % Are there points outside the left edge? 104 % If yes, how many elements are outside?
105 kk=0;
106 jj=1;
107 while(measure(1,jj)<S-edge)
108 kk=kk+1;
109 jj=jj+1;
110 if (jj>M) % all elements are outside
111 break;
112 end
113 end
114 % Are there points outside the right edge? 115 % If yes, how many elements are outside?
116 ll=0;
117 jj=M;
118 while(measure(1,jj)>edge+S)
119 ll=ll+1;
120 jj=jj-1;
121 if (jj<1) % all elements are outside
122 break;
123 end
124 end
125 % there are kk points outside the left edge, ll points outside the 126 % right edge, and M-kk-ll points inside the support
127 % copy points 128 Mnew=0; 129 if (kk>0) 130 Mnew=1; 131 end 132 Mnew=Mnew+M-kk-ll; 133 if (ll>0) 134 Mnew=Mnew+1; 135 end 136 measurenew=zeros(2,Mnew); 137 if (kk>0)
138 measurenew(1,1)=S-edge; 139 measurenew(2,1)=sum(measure(2,1:kk)); 140 measurenew(:,2:M-ll-kk+1)=measure(:,kk+1:M-ll); 141 else 142 measurenew(:,1:M-ll)=measure(:,1:M-ll); 143 end
144 % copy eventual points outside the left edge 145 if (ll>0)
146 measurenew(1,Mnew)=S+edge;
147 measurenew(2,Mnew)=sum(measure(2,M-ll+1:M)); 148 end
149 % copy the result back 150 M=Mnew;
151 measure=measurenew;
152 % The support cardinality reduction
153 % How many intervals intersect with the support? 154 j_max = ceil(edge/tau12);
155 measurenew=zeros(2,12*j_max); 156 % the index in measure
157 jj=1;
158 % the index in measurenew
159 kk=0;
160 % we work inside each interval 161 for j=-j_max:j_max-1
162 left=j*tau12+S;
163 right=left+tau12;
164 % How many points are inside the interval [left,right)?
165 number = 0; 166 if (jj>M) 167 break; 168 end 169 while (measure(1,jj)<right) 170 number=number+1; 171 jj=jj+1; 172 if (jj>M) 173 break; 174 end 175 end 176 if (number<7)
177 %just copy points
178 measurenew(:,kk+1:kk+number)=measure(:,jj-number:jj-1);
179 kk=kk+number;
180 else
181 % Decrease the number of points
182 newnumber = number;
183 work=measure(:,jj-number:jj-1);
184 while (newnumber>6)
185 % calculate the vector in ker(A) with z=1
186 u=ones(size(1:7));
187 for ll=1:7
188 for mm=1:7
190 u(ll)=u(ll)/(work(1,mm)-work(1,ll)); 191 end 192 end 193 end 194 if (u≥0) 195 u=-u; 196 end 197 % calculate alpha 198 alpha=inf; 199 mm=0; 200 for ll=1:7 201 if (u(ll)<0) 202 beta=-work(2,ll)/u(ll); 203 if (beta<alpha) 204 alpha=beta; 205 mm=ll; 206 end 207 end 208 end
209 % assign new weights
210 for ll=1:7
211 work(2,ll)=work(2,ll)+alpha*u(ll);
212 end
213 % shift the numbers
214 work(:,mm:newnumber-1)=work(:,mm+1:newnumber); 215 newnumber=newnumber-1; 216 end 217 % copy 6 points 218 measurenew(:,kk+1:kk+6)=work(:,1:6); 219 kk=kk+6; 220 end 221 end 222 M=kk; 223 measure=measurenew(:,1:M); 224 waitbar(i/N,w); 225 end 226 close(w); 227 price=dot(max(measure(1,:)-K,0),measure(2,:)); 228 price=price*exp(-r); 229 disp(price);
Appendix B
Required Objectives for the Thesis
In Appendix B, we will prove how our thesis is written under the Swedish Higher Education Agency’s requirements for the Bachelor’s degree.
B.1
Objective 1: Knowledge and understanding of the main
field
Throughout the writing process of this thesis, we have seen many intersections of common definitions, equations, and theories with already passed courses that are within the scope of our program. Which was a great benefit for us as we were able to show our knowledge and understanding in certain areas of mathematics and finance.
B.2
Objective 2: Ability to search, collect and interpret
To complete our thesis, we had to conduct an immense amount of research. We were incredibly fortunate with our main research paper Müller-Gronbach and Yaroslavtseva (2016). Which provided us with all the necessary research literature and gave us a good base and great hints on where to find further needed sources.
B.3
Objective 3: Identify, formulate and solve problems
Durning the writing process of this thesis we have come across many new formulas and equations previously unknown to us. We had to conduct deep research into them to understand their theory and application. Prof. Anatoliy Malyarenko played immense role in explaining and elaborating some heavy concepts to us. Without his help we would not have been able to understand some of the key aspects of the work. Our task was to calculate how close is the deterministic quadrature formulae result to the absolute value of the European call option price calculated by Black–Scholes model. We were able to accomplish our task by using MATLAB
script. Our numerical results prove that the deterministic quadrature formulae results are accurate.
B.4
Objective 4: Ability to present the report by writing and
orally, so others understand it
Objective 4, is considered satisfied by a well-structured thesis outlay where everything is matched with the right order of chapters and sections. To achieve the best results we have followed a preapproved thesis writing layout. By strictly following it we were able to keep the thesis structure in correct place.
B.5
Objective 5: Ability to put the report to social context
This thesis confirms that the European call option price can be accurately calculated using the deterministic quadrature formulae rather than using Black–Scholes. In the numerical results section, we show the limitations of this method and demonstrate which parameters play a key role in the accuracy of the method. We consider our thesis project as a worthy rival to Black–Scholes and in the end, we propose possible expansions of our topic for future research.
