School of Education, Culture and Communication
Division of Mathematics and Physics
BACHELOR’S DEGREE PROJECT IN MATHEMATICS
Option pricing under Black-Scholes model using stochastic Runge-Kutta
method
by
Ahmad Al-Kadri and Ali Saleh
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:
13 January 2021
Project name:
Option pricing under Black-Scholes model using stochastic Runge-Kutta method
Author(s):
Ahmad Al-Kadri and Ali Saleh
Version: 8th March 2021 Supervisor(s): Anatoliy Malyarenko Reviewer: Achref Bachouch Examiner: Ying Ni Comprising: 15 ECTS credits
Abstract
The purpose of this paper is solving the European option pricing problem under the Black– Scholes model. Our approach is to use the so-called stochastic Runge–Kutta (SRK) numerical scheme to find the corresponding expectation of the functional to the stochastic differential equation under the Black–Scholes model. Several numerical solutions were made to study how quickly the result converges to the theoretical value. Then, we study the order of convergence of the SRK method with the help of MATLAB.
Notation
{𝑋 (𝑡), 𝑡 ≥ 0} stochastic process
{𝑊 (𝑡), 𝑡 ≥ 0} Wiener process also known as the Brownian motion 𝜇 expected return of the underlying asset
𝜎 volatility of the underlying asset
𝑡 time
Í
summation ln natural logarithm IR the set of real numbers
e vector (1; ...; 1)>
ℎ step size of numerical method ODE ordinary differential equation
SDE stochastic ordinary differential equation SRK stochastic Runge–Kutta method
Acknowledgements
We would like to express our sincere gratitude to our supervisor Prof. Anatoliy Malyarenko for his guidance and motivation. His valuable advices and encouragements helped us in the writing of this thesis in many ways.
Contents
Acknowledgements 2 1 Introduction 7 1.1 Literature Review . . . 7 1.2 Research Question . . . 9 1.3 Methodology . . . 92 The Black–Scholes Model 10 2.1 The Concept of an Option. . . 10
2.2 The Black–Scholes Model . . . 10
2.3 The Black–Scholes Formula . . . 11
2.4 Risk-Neutral Valuation . . . 12
3 The Stochastic Runge–Kutta Scheme 15 3.1 The Euler Method . . . 15
3.2 A Runge–Kutta Method of Order 4 . . . 16
3.3 Generalization of Euler Method . . . 17
3.4 Stochastic Runge–Kutta Methods . . . 17
3.4.1 Weak order of convergence . . . 17
3.4.2 A d-dimensional stochastic differential equation system . . . 18
4 Numerical Simulation 22 4.1 The SRK Method . . . 22
4.2 Monte–Carlo Simulation . . . 23
4.3 Order of Convergence . . . 24
5 Numerical Results 25 5.1 The Convergence of the Call Option Price . . . 25
5.2 Order of Convergence of SRK method . . . 26
6 Conclusion 30 A MATLAB Simulation 34 A.1 Script 1 . . . 34
B Required Objectives for the Thesis 42 B.1 Objective 1: Knowledge and understanding of the main field . . . 42
B.2 Objective 2: Ability to search, collect and interpret . . . 42
B.3 Objective 3: Identify, formulate and solve problems . . . 42
B.4 Objective 4: Ability to present the report by writing and orally, so others understand it. . . 43
List of Figures
5.1 Convergence of the SRK method as 𝑚 increases.. . . 27
5.2 Order of convergence of the SRK method as ℎ2decreases. . . 28 5.3 The order of convergence with respect to varying ℎ2. . . 29
List of Tables
3.1 Butcher tableau . . . 19
3.2 Coefficients of the SRK scheme . . . 21
5.1 Approximation the call option value using the Stochastic Runge–Kutta method as 𝑛 increasing and 𝑚 = 105. . . 26
5.2 Approximation the call option value using the Stochastic Runge–Kutta method as 𝑚 increasing and 𝑛 = 75. . . 26
5.3 Order of convergence of the SRK method. . . 27
Chapter 1
Introduction
For many years, several domains such as physics, natural science, and both economics and finance, have been using ordinary differential equations (ODEs) in their studies. Many math-ematicians had ODEs as a point of interest to develop and solve numerical solutions with the development of digital computers. In numerical solution of ODEs, stiffness is a very important concept. An ODE problem is considered stiff, if the solution varies slowly, but there are some nearby solutions that vary rapidly, then the numerical method takes small steps to obtain ideal results.
Scientists and mathematicians faced many hard and stiff problems, that cannot be solved using ODEs, and here came the need for stochastic differential equations (SDEs), seeBurrage and Burrage(1996).
Stochastic differential equations are gaining importance and focus
. . . due to its application for modelling stochastic phenomena in different fields
seeTocino and Ardanuy(2002). Since, in many cases, there is no analytical solutions available for SDEs, mathematicians
. . . are forced to use numerical methods to approximate them
seeTocino and Ardanuy(2002). There are two ways to get these approximation, “strong” and “weak”. In this paper we will only focus on weak approximations.
One of the most famous numerical methods applied to SDEs is the Stochastic Runge–Kutta methods (SRK), which will be introduced in Section3.4later in this thesis.
Many models have been developed, over time, using SDEs for option price, but the most used model in the financial sector stays the Black–Scholes model.
This thesis studies the order of convergence of a numerical method applied for an option pricing under the Black–Scholes model using stochastic Runge–Kutta (SRK) method.
1.1
Literature Review
Mathematicians went through several changes and development of numerical methods begin-ning with ODEs, getting to SDEs, in order to get an approximated results of stiff problems.
These numerical methods were categorized into many groups, such as one-step, multi-steps, implicit, and explicit.
The idea of the approximated solution of the ordinary differential equations depend on the repetition of the initial value. In other words, obtaining the next value of the ODE, depends on the previous or intermediate value. This is the one-step method, since it depends only on one value to get a new approximated one.
The Euler method was the first numerical method in ODE, which was intoduced and derived by Leonard Euler (1707–1783), Fellmann (2007). This method is described as the simplest method used to solve ODEs, but not accurate for a large step-size ,which made many scientists to work on the development of this method, getting the Runge–Kutta methods by generalizing the Euler method. Generalizing the Euler method as an idea is mainly credited toRunge(1895). However, other contributors such as Heun(1900) andKutta (1901) played an important role in carrying on and realizing such generalization, where the latter provided a complete characterization of the set of 4th order Runge–Kutta methods, in addition to the proposition of the first methods of order 5, leaving the 6th order methods to be later introduced byHuťa(1956, 1957). Going back to lower orders, it is worth to mentionNyström(1925); a significant developer of methods of first order equations, who also proposed methods tailored for second order differential equations.
Many researchers have focused on the extension of the Runge–Kutta’s theory and the development of particular methods, since the beginning of technology and especially digital computers. When developing the Runge–Kutta methods, the main interest was devoted to explicit methods. Nowadays, the contributors focus on including implicit methods, which have become suitable for solving stiff differential equations.
Moving to Stochastic differential equations, the numerical methods used to solve ODEs are not suitable for use in solving SDEs, asKloeden and Platen(1989) showed. SDEs have one or many stochastic variables. Differential equations have two ways to be stochastic. One is having stochastic initial value, second is having a random variable, added to the differential equation, called the noise or the diffusion coefficient.
The methods of calculus have been extended to stochastic integrals by two mathematicians, Itô, seeLyons (2010) and Stratonovich, see the original paperStratonovič(1964) or English translationStratonovich(1966), which made stochastic differential equation to be divided into two forms named after them. In this thesis, we will work on the most used form in applied mathematics, the Itô integrals.
In recent years, Rößler Debrabant and Rößler (2008) worked and introduce a new class of second order Runge–Kutta methods for Itô stochastic differential equations, with a multidimensional Wiener process,
seeDebrabant and Rößler(2009).
There are many types of numerical methods for solving SDEs, such as Euler–Maruyama, Runge–Kutta methods and more. In this research, we will take in consideration only the Stochastic Runge–Kutta methods (SRKs).
Stochastic differential equations exist in all the models that price an option, such as Vasicek model Mamon (2004), and Heston model Rouah (2013), but the most used and known the model introduced and developed byBlack and Scholes(1973).
1.2
Research Question
In this thesis, our research task is to implement a SRK scheme to the Black–Scholes (BS) model, then we investigate the properties of convergence and check if the theoretical second-order SRK scheme has indeed a weak of convergence of second-order 2 under our problem setting via experimental studies.
1.3
Methodology
First, we implement the coefficients of the SRK method with weak order 2, given in the Butcher tableau in Debrabant and Rößler(2009), in the general 𝑠-stage Runge–Kutta method. From the obtained Runge–Kutta method, we compute an approximation to the solution. We create an algorithm on MATLAB R2020b, computing the numerical solutions of the SRK model for different numbers of simulations and steps. Finally, we study the convergence of the SRK model, by comparing our numerical results to the theoretical values.
The thesis proceeds as follows: In Chapter 2, we introduce the concepts of options, the Black–Scholes model, and the risk–neutral valuation. In Chapter 3, we review the classical Euler method, Runge–Kutta (RK) method for finding numerical solution of a ODE, and the SRK method for a 𝑑-dimensional system of SDES. In Chapter 4, we describe the implementation of the SRK method. In Chapter 5, we describe our numerical results, where we investigate the European call option prices and the corresponding order of convergence. Finally, in Chapter 6, we conclude our results.
Chapter 2
The Black–Scholes Model
2.1
The Concept of an Option
An option is a financial derivative that depends on an underlying asset’s price. An option gives the holder the right but not obligation to buy or sell at a certain price. The holder of a call option has the right to buy the asset at a stated price within a specific timeframe, while, the holder of a put optionshas the right to sell the asset at a stated price within a specific timeframe. An option can be divided regarding to exercising date into American and European option. The difference between these two options is that a European option can only be exercised at the maturity time, while an American option can be exercised at any time before time of maturity.
2.2
The Black–Scholes Model
The Black–Scholes model is used for pricing European call and put options on a non-dividend paying stock.
This model was developed by Black and Scholes (1973), Merton (1973) and extended in
Garman and Kohlhagen (1983), Orlin Grabbe (1983), and Biger and Hull (1983), and it is considered one of the famous model in the financial market, and was awarded Nobel prize in economics in 1997. The above model contains two securities, the bank account, 𝐵(𝑡), and the stock, 𝑆(𝑡). Their dynamics is given by the following system of equations.
𝑑𝐵(𝑡) = 𝑟 𝐵(𝑡) 𝑑𝑡, 𝐵(0) = 1, 𝑑𝑆(𝑡) = 𝜇𝑆(𝑡) 𝑑𝑡 + 𝜎𝑆(𝑡) 𝑑𝑊 (𝑡), 𝑆(0) = 𝑆0.
(2.1)
The first equation is an ordinary differential equation, where the spot risk–free interest rate 𝑟is constant. The second equation is a stochastic differential equation, just a shortcut which is customary to write instead of a correct stochastic integral equation
𝑆(𝑡) = 𝑆0+ ∫ 𝑡 0 𝜇 𝑆(𝑢) 𝑑𝑢 + ∫ 𝑡 0 𝜎 𝑆(𝑢) 𝑑𝑊 (𝑢),
where the first integral is an ordinary Riemann integral, while the second one is an Itô stochastic integral. The number 𝜇 ∈ R is called the expected rate of return of the stock price, while the number 𝜎 > 0 is called volatility.
It is possible to prove the following result.
Theorem 1. The stochastic process
𝑆(𝑡) = 𝑆0exp(𝜎𝑊 (𝑡) + (𝜇 − 𝜎2/2)𝑡), 𝑡 ∈ [0, 𝑇 ] (2.2)
is the unique solution of the second equation in(2.1).
The stochastic process given by Equation (2.2) is called the Geometric Brownian Motion. The Black–Scholes model, is a model of underlying dynamics. It can price any options, European and American option price. But there is so far, no exact formula for American option under this model.
2.3
The Black–Scholes Formula
The Black–Scholes formula for European call and put option is given respectively by
𝑐= 𝑆0𝑁(𝑑1) − 𝐾 𝑒−𝑟𝑇𝑁(𝑑2), 𝑝= 𝐾 𝑒−𝑟𝑇𝑁(−𝑑2) − 𝑆0𝑁(−𝑑1), 𝑑1= ln(𝑆0 𝐾) + (𝑟 + 𝜎2 2)𝑇 𝜎 √ 𝑇 , 𝑑2 = 𝑑1− 𝜎 √ 𝑇 , (2.3) where
• 𝑐 is the price of a European call option. • 𝑝 is the price of a European put option.
• 𝑁 (𝑑1) and 𝑁 (𝑑2) are the cumulative probability distribution function for the standard
normal distribution.
• 𝑆0in the time 0 price of the stock.
• 𝐾 is the strike price.
• 𝜎 is the volatility of the stock price. • 𝑟 is the risk–free interest rate. • 𝑇 is the time to maturity.
Then we look at the assumption of these five variables and plug them in one of the following equations to know the value of the option.
Black and Scholes (1973) proved the above result using the theory of partial differential equations. Merton(1973) used another approach.
2.4
Risk-Neutral Valuation
The objects of interests in this section are:
• a probability space (Ω, 𝔉, 𝑃) with a physical (historical) probability measure 𝑃;
• the filtration {𝔉𝑡 : 𝑡 ∈ [0, 𝑇 ]} generated by the standard Brownian motion {𝑊 (𝑡), 𝑡 ≤ 0};
• a non-negative 𝔉𝑡–adapted stochastic process 𝑟 (𝑡) — the time 𝑡 instantaneous spot rate;
• a positive 𝔉𝑡–adapted stochastic process 𝜎(𝑡) satisfying
∫𝑇 0 E[𝜎
2(𝑡)]𝑑𝑡 < ∞;
• an 𝔉𝑡–adapted stochastic process 𝛿(𝑡) satisfying
∫𝑇
0 E[|𝛿(𝑡) |]𝑑𝑡 < ∞— the continuous
dividend rate.
The money-market account 𝐵(𝑡) = 𝑆0(𝑡) and the stock price 𝑆(𝑡) = 𝑆1(𝑡) evolve as
𝑑𝐵(𝑡) = 𝑟 (𝑡) 𝐵(𝑡)𝑑𝑡, 𝐵(0) = 1,
𝑑𝑆(𝑡) = [𝜇(𝑡)𝑆 − 𝛿(𝑡)]𝑑𝑡 + 𝜎(𝑡)𝑆𝑑𝑊 , 𝑆(0) = 𝑆0.
(2.4)
Assume that (2.4) has a non-negative, square integrable solution.
Let 𝐶 (𝑡) be the price of a European call option written on the stock with payoff function ℎ(𝑡) and maturity 𝑇 . We would like to calculate 𝐶 (𝑡).
Theorem 2 (Girsanov). Let 𝛽(𝑡) be a stochastic process, adapted to the filtration {𝔉𝑡 : 𝑡 ∈
[0, 𝑇 ]} generated by the standard Brownian motion {𝑊 (𝑡), 𝑡 ≤ 0} on a probability space
(Ω, 𝔉,P), and let E exp 1 2 ∫ 𝑇 0 𝛽2(𝑢)𝑑𝑢 <∞. (2.5) Put 𝑌(𝑇 ) = exp ∫ 𝑇 0 𝛽(𝑡)𝑑𝑊 (𝑡) − 1 2 ∫ 𝑇 0 𝛽2(𝑡)𝑑𝑡 and define ˜ P( 𝐴) = 𝐸 [1𝐴𝑌(𝑇 )], 𝐴 ∈ 𝔉𝑇.
Then, the process
˜
𝑊(𝑡) = 𝑊 (𝑡) − ∫ 𝑡
0
𝛽(𝑢)𝑑𝑢 0 ≤ 𝑡 ≤ 𝑇
is a standard Brownian motion on (Ω, 𝔉𝑇, ˜P).
Replicate the claim by a self-financing dynamic strategy of trading 𝑏(𝑡) units of the money market and 𝜃 (𝑡) units of the underlying stock. That is,
𝐶(𝑡) = 𝑏(𝑡) 𝐵(𝑡) + 𝜃 (𝑡)𝑆(𝑡) = 𝐶 (0) + ∫ 𝑡 0 𝑏(𝑢)𝑑𝐵(𝑢) + ∫ 𝑡 0 𝜃(𝑡)𝑑𝐺 (𝑢), (2.6)
where 𝐺 (𝑡), 𝑡 ∈ [0, 𝑇 ] is the gain process:
𝐺(𝑡) = 𝑆(𝑡) + ∫ 𝑡
0
𝛿(𝑢)𝑑𝑢. (2.7)
Find an equation for the claim price, by calculating the Itô differential of both hand sides of (2.7), using (2.4) and (2.6):
𝑑 𝐺 = 𝑑𝑆 + 𝛿(𝑡)𝑑𝑡 = 𝑆 [𝜇(𝑡)𝑑𝑡 + 𝜎 (𝑡)𝑑𝑊 ]. Define a new process 𝑊∗(𝑡) by
𝜎(𝑡)𝑑𝑊∗ = ( 𝜇(𝑡) − 𝑟 (𝑡))𝑑𝑡 + 𝜎 (𝑡)𝑑𝑊 . (2.8) Then
𝑑 𝐺 = 𝑆 [𝑟 (𝑡)𝑑𝑡 + 𝜎 (𝑡)𝑑𝑊∗]. Substitute this to (2.8) and use (2.4):
𝐶(𝑡) = 𝐶 (0) + ∫ 𝑡 0 𝑏(𝑢)𝑟 (𝑢) 𝐵(𝑢)𝑑𝑢 + ∫ 𝑡 0 𝜃(𝑢)𝑆(𝑢)𝑟 (𝑢)𝑑𝑢 + ∫ 𝑡 0 𝜃(𝑢)𝑆(𝑢)𝜎 (𝑢)𝑑𝑊∗ = 𝐶 (0) + ∫ 𝑡 0 𝑟(𝑢) [𝑏 (𝑢) 𝐵(𝑢) + 𝜎 (𝑢)𝑆(𝑢)]𝑑𝑢 + ∫ 𝑡 0 𝜃(𝑢)𝑆(𝑢)𝜎 (𝑢)𝑑𝑊∗ = 𝐶 (0) + ∫ 𝑡 0 𝑟(𝑢)𝐶 (𝑢)𝑑𝑢 + ∫ 𝑡 0 𝜃(𝑢)𝑆(𝑢)𝜎 (𝑢)𝑑𝑊∗. It follows that 𝑑𝐶 = 𝑟 (𝑡)𝐶 (𝑡)𝑑𝑡 + 𝜃 (𝑡)𝜎 (𝑡)𝑆(𝑡)𝑑𝑊∗. (2.9) Solving (2.9) is unknown, because 𝑊∗(𝑡) is not a Brownian motion. However, one can check that the process
𝛽(𝑡) =
𝑊(𝑡) − 𝜇(𝑡) 𝜎(𝑡)
satisfies condition (2.5). By Girsanov theorem, there exists a probability measure P∗ that makes the process 𝑊∗(𝑡) a standard Brownian motion.
The solution to (2.9) is 𝐶(𝑡) = 𝐶 (0) exp ∫ 𝑡 0 𝑟(𝑢)𝑑𝑢 + ∫ 𝑡 0 𝜃(𝑢)𝜎 (𝑢)𝑆(𝑢)𝑑𝑊∗. Set 𝐵(𝑡) as a numéraire and put 𝐶∗(𝑡) = 𝐶 (𝑡)/𝐵(𝑡). then,
𝐶∗(𝑡) = 𝐶∗(0) + ∫ 𝑡
0
𝜃(𝑢)𝜎 (𝑢)𝑆∗(𝑢)𝑑𝑊∗,
where 𝑆∗(𝑢) = 𝑆(𝑢)/𝐵(𝑢) is the denominated price. The relation of this equation to option pricing is that the denominated claim price is a martingale underP∗. Because 𝐶 (𝑇 ) = ℎ(𝑆(𝑇 )), we obtain by no arbitrage 𝐶(𝑡) = 𝐵(𝑡)E∗ ℎ(𝑆(𝑇 )) 𝐵(𝑇 ) |𝔉𝑡 . (2.10)
Therefore, we have to know two-dimensional distribution of (𝑆(𝑇 ), 𝐵(𝑇 )). Substituting (2.9) into (2.4), we obtain, underP∗,
𝑑𝑆= [𝑟 (𝑡)𝑆 − 𝛿(𝑡)]𝑑𝑡 + 𝜎 (𝑡)𝑆𝑑𝑊∗. (2.11) The economical sense of this equation is as follows: in the world of P∗ the mean rate of return of the stock is equal to that of the risk-free security. In this world, investors would be risk-neutral, since otherwise the risky stock cannot exist. So, the measure P∗ is called the
risk-neutral probability measure.
In particular, for the Black–Scholes model, we have 𝑟 (𝑡) = 𝑟, 𝜇(𝑡) = 𝜇, 𝜎(𝑡) = 𝜎, 𝛿(𝑡) = 0. Equation (2.11) becomes
𝑑𝑆= 𝑟 𝑆𝑑𝑡 + 𝜎𝑆𝑑𝑊∗(𝑡). (2.12) This Black–Scholes stochastic differential equation in the risk-neutral world will be numerically studied below.
Chapter 3
The Stochastic Runge–Kutta Scheme
The Runge–Kutta methods are numerical methods used to approximate numerical solutions of ODEs and SDEs. The most commonly known RK methods for ordinary (non-stochastic) differential equations are the Euler method, Heun’s method (improved Euler’s method) and the classical Runge–Kutta method, known as RK4.
Discretization method
A method is called a discretization method, when it uses a set of discrete points to approx-imate the solution 𝑦(𝑡) of an ODE. Almost all numerical methods are discretization methods. Let 𝐼 = [𝑡0, 𝑇] be an interval, and 𝑡0 < 𝑡1 < · · · < 𝑡𝑁 = 𝑇be an ordered subset of 𝐼, created in 𝑁 intervals, where 𝑡𝑛= 𝑡0+ 𝑛ℎ for 𝑛 = 0, 1, . . . , 𝑁 and ℎ =
(𝑇 −𝑡0)
𝑁 be the step size.
An initial value problem for a first-order ODE has the following form:
𝑦0= 𝑓 (𝑡, 𝑦 (𝑡)), 𝑦(𝑡0) = 𝑦0, (3.1) where the slope of the tangent line of 𝑦(𝑡) is given by the function 𝑓 at each point (𝑡, 𝑦(𝑡)).
A discretization method calculates a sequence 𝑦1, 𝑦2, . . . , 𝑦𝑁 of approximations to the
values 𝑦(𝑡1), 𝑦(𝑡2), . . . , 𝑦(𝑡𝑁). Define a function ˜𝑦 : 𝐼 → R by the formula
˜𝑦(𝑡) = 𝑦𝑘−1+ (𝑡 − 𝑡𝑘−1) 𝑓 (𝑡𝑘−1, 𝑦𝑘−1), 𝑡𝑘−1 < 𝑡 ≤ 𝑡𝑘.
In words, the function ˜𝑦(𝑡) is the piecewise linear interpolation of the sequence of approxima-tions.
Definition 1. Butcher(2000) A discretization method has order of convergence at least 𝑝 for a class of equations, if the global error, max𝑡0≤𝑡 ≤𝑇 |𝑦 (𝑡) − ˜𝑦(𝑡)|, is bounded from above by a
constant times ℎ𝑝
, and exactly 𝑝, if there is an equation for which the global error is bounded from below by another constant times ℎ𝑝
.
3.1
The Euler Method
The Euler method is at least order 1 numerical method, that is known for analyzing an ordinary differential equation (ODE). The Euler method is simple but inaccurate for large step size. This
method was derived by Leonard Euler (1707–1783), seeWilson(2017), and it was published in his bookEuler(1768) in 1768–1770, seeButcher(2016).
The Euler method is described as:
𝑦𝑛+1 = 𝑦𝑛+ ℎ 𝑓 (𝑡𝑛, 𝑦𝑛). (3.2) The method is based on repetitive iteration, when evaluating 𝑓 at (𝑡𝑛, 𝑦𝑛) to approximate the
next value 𝑦𝑛+1. At an interval 𝐼 = [𝑡0, 𝑇], the method applied will initiate an approximation
of discrete values 𝑦𝑖 for 𝑖 = 1, 2, . . . until it reaches the ideal approximation of the numerical
solution at 𝑇 .
To see why the Euler method has the form (3.2), we proceed as follows. By integrating (3.1) from 𝑡𝑛to 𝑡𝑛+1we get the following:
∫ 𝑡𝑛+1 𝑡𝑛 𝑦0(𝑡)𝑑𝑡 = ∫ 𝑡𝑛+1 𝑡𝑛 𝑓(𝑡, 𝑦 (𝑡))𝑑𝑡. From the fundamental theorem of calculus:
𝑦(𝑡𝑛+1) − 𝑦 (𝑡𝑛) = ∫ 𝑡𝑛+1
𝑡𝑛
𝑓(𝑡, 𝑦 (𝑡))𝑑𝑡.
Replacing 𝑦(𝑡𝑛+1) and 𝑦(𝑡𝑛) with 𝑦𝑛+1 and 𝑦𝑛, in the left hand side, and applying the
Rectangular Rule, in the right hand side, gives: ∫ 𝑡𝑛+1
𝑡𝑛
𝑓(𝑡, 𝑦 (𝑡))𝑑𝑡 ≈ ℎ 𝑓 (𝑡𝑛, 𝑦(𝑡𝑛)) ≈ ℎ 𝑓 (𝑡𝑛, 𝑦𝑛), and finally we get the Euler method:
𝑦𝑛+1 ≈ 𝑦𝑛+ ℎ 𝑓 (𝑡𝑛, 𝑦𝑛).
There are many hard and stiff problems that cannot be solved using the Euler method; problems which demand high–order methods. The ideal method, in this case, is instead that of the implicit Runge–Kutta, seeButcher(2000).
3.2
A Runge–Kutta Method of Order 4
The most known and used numerical method is the classical Runge–Kutta method, also known as the Runge–Kutta of order 4 (RK4), for its accuracy and efficiency.
The RK4 method is described as follow:
𝑦𝑛+1 = 𝑦𝑛+ ℎ 6(𝑘1+ 2𝑘2+ 2𝑘3+ 𝑘4) (3.3) where, 𝑘1= 𝑓 (𝑡𝑛, 𝑦𝑛,), 𝑘2= 𝑓 (𝑡𝑛+ ℎ 2, 𝑦𝑛+ 𝑘1 2 ), 𝑘3= 𝑓 (𝑡𝑛+ ℎ 2, 𝑦𝑛+ 𝑘2 2 ), 𝑘4= 𝑓 (𝑡𝑛+ ℎ, 𝑦𝑛+ 𝑘3),
where 𝑘1, . . . , 𝑘4are four points used to evaluate the derivative, and the approximation of 𝑦𝑛+1 is found using the weighted average.
3.3
Generalization of Euler Method
Kutta formulated the general scheme of the RK methods, seeIxaru and Van Den Berghe(2004, Chapter 6), the generalization of the Euler method with 𝑠-stage.
Let the number of stages, 𝑠, be an integer, and 𝑎21, 𝑎31, 𝑎32, . . . , 𝑎𝑠1, 𝑎𝑠2, . . . , 𝑎𝑠,𝑠−1, 𝑏1,
𝑏2, . . . , 𝑏𝑠, 𝑐1, 𝑐2, . . . , 𝑐𝑠 be real coefficients.
The formulated Runge–Kutta method is given by
𝑦𝑛+1 = 𝑦𝑛+ ℎ(𝑏1𝑘1+ . . . + 𝑏𝑠𝑘𝑠) = 𝑦𝑛+ ℎ 𝑠 ∑︁ 𝑖=1 𝑏𝑖𝑘𝑖 (3.4) where 𝑘1= 𝑓 (𝑡𝑛, 𝑦𝑛), 𝑘2= 𝑓 (𝑡𝑛+ 𝑐2ℎ, 𝑦𝑛+ ℎ𝑎21𝑘1), 𝑘3= 𝑓 (𝑡𝑛+ 𝑐3ℎ, 𝑦𝑛+ ℎ(𝑎31𝑘1+ 𝑎32𝑘2)), . . . 𝑘𝑠 = 𝑓 (𝑡𝑛+ 𝑐𝑠ℎ, 𝑦𝑛+ ℎ(𝑎𝑠1𝑘1+ ... + 𝑎𝑠,𝑠−1𝑘𝑠−1)), where, 𝑐𝑖 = Í𝑖−1 𝑗=1𝑎𝑖 𝑗 and Í𝑠 𝑖=1𝑏𝑖 =1 for 𝑖 = 1, 2, . . . , 𝑠.
The scheme for an explicit Runge–Kutta method is defined in the following so–called Butcher tableau; 𝑐1 𝑐2 𝑎21 𝑐3 𝑎31 𝑎32 . . . . . . . . . . . . . . . 𝑐𝑠 𝑎𝑠1 𝑎𝑠2 𝑎𝑠3 . . . 𝑎𝑠,𝑠−1 𝑏1 𝑏2 𝑏3 . . . 𝑏𝑠−1 𝑏𝑠
3.4
Stochastic Runge–Kutta Methods
”A stochastic differential equation is a differential equation whose coefficients are random numbers or random functions of the independent variable (or variables).” van Kampen(1981) Highlighting of the Stochastic Runge–Kutta (SRK) for weak approximation — ” focuses on the expectation of functionals of the solution”Debrabant and Rößler(2009)” — has become a persistent need.
3.4.1
Weak order of convergence
Assuming that the Borel-measurable coefficients a : 𝐼 × R𝑑 → R𝑑
and the diffusion b: 𝐼×R𝑑 →R𝑑×𝑚are measurable functions that satisfy a Lipschitz and a linear growth condition then the Existence and Uniqueness TheoremGao(2012) applies. In the following, let 𝑏𝑗
(𝑡, 𝑥) = (𝑏𝑖 . 𝑗(𝑡, 𝑥))
Let 𝐶𝑙 𝑃(R
𝑑
,R) denote the space of all 𝑔 ∈ 𝐶𝑙(R𝑑,R) that fulfill a polynomial growth condition and let 𝑔 ∈ 𝐶𝑘 ,𝑙
𝑃 (𝐼 ×R 𝑑 ,R) if 𝑔(., 𝑥) ∈ 𝐶𝑘(𝑖,R) and 𝑔(𝑡, .) ∈ 𝐶𝑙 𝑃(R 𝑑 ,R) for all 𝑡 ∈ 𝐼 and 𝑥 ∈R𝑑
Debrabant and Rößler(2009). 𝑌 = (𝑌 (𝑡))𝑡∈𝐼ℎ, 𝑋, 𝐼ℎ, 𝐶
2( 𝑝+1) 𝑝 (R
𝑑
,R).
” A time discrete approximation 𝑌 = (𝑌 (𝑡))𝑡∈𝐼ℎ converges weakly with order
𝑝to 𝑋 as ℎ → 0 at time 𝑡 ∈ 𝐼ℎif for each 𝑓 ∈ 𝐶2( 𝑝+1)𝑝 (R𝑑,R) exist a constant 𝐶𝑓 and a finite 𝛿0 > 0 such that |E( 𝑓 ( 𝑋 (𝑡))) −E( 𝑓 (𝑌 (𝑡))) | ≤ 𝐶𝑓ℎ
𝑝
holds for each ℎ ∈]0, 𝛿0[.”Debrabant and Rößler(2009).
3.4.2
A d-dimensional stochastic differential equation system
Let 𝑑 and 𝑚 be two positive integers. Let 𝑋1(𝑡), . . . , 𝑋𝑑(𝑡) be stochastic processes defined on
the time interval [0, 𝑇 ]. In applications to Financial Engineering, the real number 𝑇 denotes the maturity of a financial instrument. The stochastic processes describe the prices of securities, the volatilities, etc. Let 𝑊1(𝑡), . . . , 𝑊𝑚(𝑡) be independent Brownian motions. Consider a
market model described by the system of stochastic differential equations
𝑑 𝑋𝑖(𝑡) = 𝑎𝑖(𝑡, 𝑋1(𝑡), . . . , 𝑋𝑑(𝑡)) 𝑑𝑡 + 𝑚 ∑︁ 𝑗=1 𝑏𝑖 𝑗(𝑡, 𝑋1(𝑡), . . . , 𝑋𝑑(𝑡)) 𝑑𝑊𝑗(𝑡), 𝑋𝑖(0) = 𝑥𝑖.
To shorten notation, we join the processes 𝑋1(𝑡), . . . , 𝑋𝑑(𝑡) into the singleR 𝑑
-valued stochastic process X(𝑡) = ( 𝑋1(𝑡), . . . , 𝑋𝑑(𝑡))
>
. Similarly, we join the functions 𝑎𝑖 into the R 𝑑
-valued function a, the functions 𝑏𝑖 𝑗 into the function 𝑏 taken values in the set of rectangular matrices
with 𝑑 rows and 𝑚 columns, the Brownian motions 𝑊1(𝑡), . . . , 𝑊𝑚(𝑡) into the R 𝑚
-valued Brownian motion W(𝑡), and the numbers 𝑥𝑖into the vector x0∈R
𝑑
. The system takes the form
𝑑X(𝑡) = a(𝑡, X(𝑡)) 𝑑𝑡 + 𝑏(𝑡, X(𝑡)) 𝑑W(𝑡),
X(𝑡0) = x0.
(3.5)
The notation used by Debrabant and Rößler (2008) differs from that of Butcher (2000) explained in Section3.3. Specifically, let 𝑠 be a positive integer called the number of stages. Choose 5 vectors inR𝑠
and call them 𝜶, 𝜷(1), 𝜷(2), 𝜷(3), and 𝜷(4). Choose 6 square matrices with 𝑠 rows and real-valued entries and call them 𝐴(0), 𝐴(1), 𝐴(2), 𝐵(0), 𝐵(1), and 𝐵(2). Assume that the matrices 𝐴(0), 𝐴(1), 𝐵(0), and 𝐵(1) are strictly lower triangular, that is, all their matrix entries lying on the main diagonal and over it, are equal to 0. Denote c(𝑞) = 𝐴(𝑞)e, where 0 ≤ 𝑞 ≤ 2 and e is the vector with 𝑑 components all equal to 1.
Let the points 0 = 𝑡0 < 𝑡1 < · · · < 𝑡𝑁 = 𝑇 be given, and let ℎ𝑛 = 𝑡𝑛+1− 𝑡𝑛. Afterwards, an explicit 𝑠-stage SRK method is proposed to define a 𝑑-dimensional approximation process
𝑌𝑛=𝑌 (𝑡𝑛) with 𝑌0=𝑥0and 𝑌𝑛+1= 𝑌𝑛+ 𝑠 ∑︁ 𝑖=1 𝛼𝑖𝑎(𝑡𝑛+ 𝑐 (0) 𝑖 ℎ𝑛, 𝐻 (0) 𝑖 ) ℎ𝑛+ 𝑠 ∑︁ 𝑖=1 𝑚 ∑︁ 𝑘=1 𝛽(1) 𝑖 𝑏 𝑘 (𝑡𝑛+ 𝑐 (1) 𝑖 ℎ𝑛, 𝐻 (𝑘) 𝑖 ) ˆ𝐼(𝑘) + 𝑠 ∑︁ 𝑖=1 𝑚 ∑︁ 𝑘=1 𝛽(2) 𝑖 𝑏 𝑘 (𝑡𝑛+ 𝑐 (1) 𝑖 ℎ𝑛, 𝐻 (𝑘) 𝑖 ) ˆ 𝐼(𝑘,𝑘) √ ℎ𝑛 + 𝑠 ∑︁ 𝑖=1 𝑚 ∑︁ 𝑘=1 𝛽(3) 𝑖 𝑏 𝑘 (𝑡𝑛+ 𝑐 (2) 𝑖 ℎ𝑛,𝐻ˆ (𝑘) 𝑖 ) ˆ𝐼(𝑘) + 𝑠 ∑︁ 𝑖=1 𝑚 ∑︁ 𝑘=1 𝛽(4) 𝑖 𝑏 𝑘 (𝑡𝑛+ 𝑐 (2) 𝑖 ℎ𝑛,𝐻ˆ (𝑘) 𝑖 ) √ ℎ𝑛, (3.6)
for 𝑛 = 0, 1, . . . , 𝑁 − 1 with the values
𝐻(0) 𝑖 = 𝑌𝑛+ 𝑖−1 ∑︁ 𝑗=1 𝐴(0) 𝑖 𝑗 𝑎(𝑡𝑛+ 𝑐 (0) 𝑗 ℎ𝑛, 𝐻 (0) 𝑗 ) ℎ𝑛+ 𝑖−1 ∑︁ 𝑗=1 𝑚 ∑︁ 𝑙=1 𝐵(0) 𝑖 𝑗 𝑏 𝑙 (𝑡𝑛+ 𝑐 (1) 𝑗 ℎ𝑛, 𝐻 (𝑙) 𝑗 ) ˆ𝐼(𝑙), 𝐻(𝑘) 𝑖 = 𝑌𝑛+ 𝑖−1 ∑︁ 𝑗=1 𝐴(1) 𝑖 𝑗 𝑎(𝑡𝑛+ 𝑐 (0) 𝑗 ℎ𝑛, 𝐻 (0) 𝑗 ) ℎ𝑛+ 𝑖−1 ∑︁ 𝑗=1 𝐵(1) 𝑖 𝑗 𝑏 𝑘 (𝑡𝑛+ 𝑐 (1) 𝑗 ℎ𝑛, 𝐻 (𝑘) 𝑗 ) √ ℎ(𝑛), ˆ 𝐻(𝑘) 𝑖 = 𝑌𝑛+ 𝑠 ∑︁ 𝑗=1 𝐴(2) 𝑖 𝑗 𝑎(𝑡𝑛+ 𝑐 (0) 𝑗 ℎ𝑛, 𝐻 (0) 𝑗 ) ℎ𝑛+ 𝑠 ∑︁ 𝑗=1 𝑚 ∑︁ 𝑙=1 𝑙≠𝑘 𝐵(2) 𝑖 𝑗 𝑏 𝑙 (𝑡𝑛+ 𝑐 (1) 𝑗 ℎ𝑛, 𝐻 (𝑙) 𝑗 ) ˆ 𝐼(𝑘,𝑙) √ ℎ𝑛 (3.7)
for 𝑖 = 1, . . . , 𝑠 and 𝑘 = 1, . . . , 𝑚, seeDebrabant and Rößler(2009).
A way to simplify the presentation of Runge–Kutta methods is by using the Butcher tableau, which in this case, is used to resolve the coefficients of SRK method.
Table 3.1: Butcher tableau
c(0) 𝐴(0) 𝐵(0)
c(1) 𝐴(1) 𝐵(1)
c(2) 𝐴(2) 𝐵(2)
𝜶> 𝜷(1)> 𝜷(2)> 𝜷(3)> 𝜷(4)>
This SRK method has two random variables, ˆ𝐼(𝑘) and ˆ𝐼(𝑘,𝑙), defined by
P( ˆ𝐼(𝑘) = ± √︁ 3ℎ𝑛) = 1 6, P( ˆ𝐼(𝑘) =0) = 2 3, ˆ 𝐼(𝑘,𝑙) = 1 2( ˆ𝐼(𝑘)𝐼ˆ(𝑙)− √ ℎ𝑛𝐼˜(𝑘)), 𝑘 < 𝑙 , 1 2( ˆ𝐼(𝑘)𝐼ˆ(𝑙)+ √ ℎ𝑛𝐼˜(𝑙)), 𝑙 < 𝑘 , 1 2( ˆ𝐼 2 (𝑘) − ℎ𝑛), 𝑘 = 𝑙,
where the random variable ˜𝐼(𝑘) has probability distribution
P( ˜𝐼(𝑘) = ± √
ℎ𝑛) = 1 2.
Theorem 3(Debrabant and Rößler(2009)). Let 𝑎𝑖
, 𝑏𝑖, 𝑗 ∈ 𝐶2,4
𝑃 (𝐼 ×R 𝑑
,R) for 1 ≤ 𝑖 ≤ 𝑑, 1 ≤ 𝑗 ≤ 𝑚. If the coefficients of the SRK method satisfy the following conditions
1. 𝜶>e =1. 2. 𝜷4>e =0. 3. 𝜷(3)>e =0 4. ( 𝜷(1)>e)2=1. 5. 𝜷(2)>e =0. 6. 𝜷(1)>𝐵(1)e =0. 7. 𝜷(4)>𝐴(2)e =0. 8. 𝜷(3)>𝐵(2)e =0. 9. 𝜷(4)>(𝐵(2)e)2 =0,
then the method has weak order 1. If in addition, 𝑎𝑖
, 𝑏𝑖, 𝑗 ∈ 𝐶3,6
𝑃 (𝐼 ×R 𝑑
,R) for 1 ≤ 𝑖 ≤ 𝑑, 1 ≤ 𝑗 ≤ 𝑚, satisfy 50 additional
conditions given at Debrabant and Rößler (2009, pp. 584–585), then the method has weak order 2.
In this theorem, the notation 𝐶𝑃2,4(𝐼 ×R𝑑
,R) means that the corresponding function is twice continuously differentiable in the first argument and 4 times continuously differentiable in the second one, moreover, the function and all its partial derivatives up to the given order grow at infinity not faster than a polynomial. The sense of the notation 𝐶𝑃3,6(𝐼 ×R𝑑
,R) is similar. The coefficients of weak second order used in this paper, are proposed in Debrabant and Rößler(2009), and presented in the following table,
Note that Table3.1 is just a set of arbitrary matrices and vectors written in a convenient compact form. Table3.2is the set of particular matrices and vectors chosen carefully in such a way that under some conditions the corresponding stochastic RK method has order 2.
Table 3.2: Coefficients of the SRK scheme 0 1 2 1 2 6−√6 10 1 -1 2 3+2 √ 6 5 0 0 342 491 342 491 3 √︃ 38 491 342 491 342 491 0 −3 √︃ 38 491 0 0 0 0 0 0 0 −214 513 √︃ 1105 991 − 491 513 √︃ 221 4955 − 491 513 √︃ 221 4955 0 0 0 214513 √︃ 1105 991 491 513 √︃ 221 4955 491 513 √︃ 221 4955 1 6 2 3 1 6 193 684 491 1368 491 1368 0 1 6 √︃ 491 38 − 1 6 √︃ 491 38 −4955 7072 4955 14144 4955 14144 0 − 1 8 √︃ 4955 221 1 8 √︃ 4955 221
Chapter 4
Numerical Simulation
In the previous chapters, we discussed the stochastic Runge–Kutta schemes and considered a special scheme described inDebrabant and Rößler(2009). The above scheme is applicable to a system of stochastic differential equations (3.5). In the case when the system contains only one equation, the scheme simplifies. We show the simplifications in Section4.1.
4.1
The SRK Method
In system (3.5), we put X(𝑡) = 𝑆(𝑡), the price of a risky security, a(𝑡, X(𝑡)) = 𝑟 𝑆(𝑡), 𝑏(𝑡, X(𝑡)) = 𝜎𝑆(𝑡). The above system takes the form
𝑑𝑆(𝑡) = 𝑟 𝑆(𝑡) 𝑑𝑡 + 𝜎𝑆(𝑡) 𝑑𝑊 (𝑡), 𝑆(0) = 𝑆0.
This system is the celebrated Black–Scholes model, under the risk–neutral probabilityP∗ with zero dividend yield for the risky asset. Here 𝑟 is the constant risk-free spot interest rate and 𝜎 is the constant volatility parameter, see also Equation (2.12) above.
The next step is to apply the numerical method (3.6)–(3.7) to the Black–Scholes model. First, according to (Debrabant and Rößler,2009, Section 4), we can put 𝐴(2) to be equal the zero matrix if 𝑚 = 1, that is, when the system contains only one Wiener process. Second, we introduce a uniform discretization, that is, we solve the Black–Scholes stochastic differential equation numerically on the interval [0, 𝑇 ] and put 𝑡𝑖 = 𝑖ℎ, where 0 ≤ 𝑖 ≤ 𝑁 and ℎ = 𝑇 /𝑁.
Applying this in the SRK scheme for 𝑚 = 1, 𝑠 = 3, 𝑖 = 1, 2, 3, 𝑘 = 𝑗 = 1 with a constant time step ℎ = 𝑇 /𝑁, we get the following
𝑌𝑛+1 = 𝑌𝑛+ 3 ∑︁ 𝑖=1 𝛼𝑖𝑟 𝐻 (0) 𝑖 ℎ+ 3 ∑︁ 𝑖=1 𝛽(1) 𝑖 𝜎 𝐻 (1) 𝑖 ˆ 𝐼(1)+ 3 ∑︁ 𝑖=1 𝛽(2) 𝑖 𝜎 𝐻 (1) 𝑖 ˆ 𝐼(1,1) √ ℎ + 3 ∑︁ 𝑖=1 𝛽(3) 𝑖 𝜎 ˆ 𝐻(1) 𝑖 ˆ 𝐼(1) + 3 ∑︁ 𝑖=1 𝛽(4) 𝑖 𝜎𝐻ˆ (1) 𝑖 √ ℎ
with 𝐻(0) 𝑖 = 𝑌𝑛+ 𝑖−1 ∑︁ 𝑗=1 𝐴(0) 𝑖 𝑗 𝑟 𝐻 (0) 𝑗 ℎ+ 𝑖−1 ∑︁ 𝑗=1 𝐵(0) 𝑖 𝑗 𝜎 𝐻 (1) 𝑗 𝐼ˆ(1), 𝐻(1) 𝑖 = 𝑌𝑛+ 𝑖−1 ∑︁ 𝑗=1 𝐴(1) 𝑖 𝑗 𝑟 𝐻(0) 𝑗 ℎ+ 𝑖−1 ∑︁ 𝑗=1 𝐵(1) 𝑖 𝑗 𝜎 𝐻(1) 𝑗 √ ℎ, ˆ 𝐻(1) 𝑖 = 𝑌𝑛.
Taking into account the values in Table3.2, we obtain
𝐻(0) 1 = 𝐻 (1) 1 = ˆ𝐻 (1) 𝑖 = 𝑌𝑛, 𝐻(0) 2 = 𝑌𝑛+ 𝐴 (0) 21𝑟𝑌𝑛ℎ+ 𝐵 (0) 21 𝜎𝑌𝑛𝐼ˆ(1), 𝐻(0) 3 = 𝑌𝑛+ 𝐴 (0) 31𝑟𝑌𝑛ℎ+ 𝐴 (0) 32𝑟 𝐻 (0) 2 ℎ+ 𝐵 (0) 31𝜎𝑌𝑛𝐼ˆ(1), 𝐻(1) 2 = 𝑌𝑛+ 𝐴 (1) 21𝑟𝑌𝑛ℎ+ 𝐵 (1) 21 𝜎𝑌𝑛 √ ℎ, 𝐻(1) 3 = 𝑌𝑛+ 𝐴 (1) 31𝑟𝑌𝑛ℎ+ 𝐵 (1) 31 𝜎𝑌𝑛 √ ℎ.
4.2
Monte–Carlo Simulation
Let 𝑀 be the number of simulations of the trajectory of the stock price 𝑆(𝑡). Let 𝑌𝑁(1), 𝑌𝑁(2), . . . , 𝑌𝑁( 𝑀) be the values of the variable 𝑌𝑁 obtained by using the equations of Section4.1.
Denote the theoretical price of the call option and its approximation by 𝑐 and ˆ𝑐𝑀respectively and apply the Monte Carlo estimate that corresponds to the risk-neutral valuation formula (2.10): ˆ 𝑐𝑀 = 1 𝑀 𝑒−𝑟𝑇 𝑀 ∑︁ 𝑖=1 max{𝑌𝑁(𝑖) − 𝐾, 0}.
According to the Law of Large Numbers, the estimator ˆ𝑐𝑀 is unbiased. According to the Central Limit Theorem, the distribution of the estimator ˆ𝑐𝑀 converges to normal distribution
function, as 𝑀 increases without bound. We may compute its confidence interval with a confidence level 100(1-𝛼) by ( ˆ𝑐𝑀 − 𝑧𝛼 2 𝑠𝑐 √ 𝑀 , 𝑐ˆ𝑀 + 𝑧𝛼 2 𝑠𝑐 √ 𝑀 ), where 𝑠𝑐 = v u t 1 𝑀− 1 𝑀 ∑︁ 𝑖=1 (𝑌(𝑖) 𝑁 − ˆ 𝑐𝑀)2.
4.3
Order of Convergence
For any ℎ1and ℎ2and for the corresponding discrete approximations 𝑌1(𝑡) and 𝑌2(𝑡),
|E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌1(𝑡))] | ≤ 𝐶𝑓ℎ 𝑝 1, |E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌2(𝑡))] | ≤ 𝐶𝑓ℎ 𝑝 2.
Assume that this order is exact, that is,
|E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌1(𝑡))] | ≈ 𝐶𝑓ℎ 𝑝 1, |E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌2(𝑡))] | ≈ 𝐶𝑓ℎ 𝑝 2. 00
If we divide the first equation by the second one, the unknown constant 𝐶𝑓 cancels:
|E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌1(𝑡))] | |E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌2(𝑡))] | ≈ ℎ 𝑝 1 ℎ 𝑝 2
Take logarithm of both hand sides:
log(|E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌1(𝑡))] |) − log(|E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌2(𝑡))] |) ≈ 𝑝 (log ℎ1− log ℎ2),
and we propose an estimate for 𝑝:
𝑝 ≈ log(|
E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌1(𝑡))] |) − log(|E[ 𝑓 ( 𝑋 (𝑡))] −E[ 𝑓 (𝑌2(𝑡))] |) log ℎ1− log ℎ2
. (4.1)
The global order for numerical methods applied to stochastic differential equations is 𝑂 (ℎ𝑝
). In section 5.3, we investigate the behaviour of this estimator by varying the size of ℎ2.
Chapter 5
Numerical Results
In this chapter, we study the convergence of the stochastic Runge–Kutta (SRK) method. Using the SRK method and the MATLAB simulation, we estimate several prices of a call option by increasing the number of simulations 𝑚 and the number of intervals 𝑁. Then, we study the order of convergence for the method and see how what it is and how close it is for the theoretical order of convergence. All the obtained results are shown in graphs and tables. the purpose of this section is to see how quickly the result converges to the true value with the following initial and parameter values.
• 𝑆0=100 𝑇 =0.25 • 𝐾 = 95 𝑟 =0.1 • 𝜎 = 0.5
The main variables, in our numerical simulation, are the number of simulations, the number of intervals and the number of attempts which remains constant.
By the first equation in (2.3), we get the theoretical European call option price, 𝑐 = 13.6953$. We ran the program several times with different data, and of course, using the Stochastic Runge–Kutta method, with main variables 𝑚 and 𝑛 equal to 106and 75 respectively, we obtained the approximate European call option value, where 𝑐 = 13.6933$.
5.1
The Convergence of the Call Option Price
In this section, we check the convergence of the European call option price by changing some data and fixing others. Table5.1 shows the convergence of 𝑐, as the number of intervals (𝑛) increases, while Table 5.2 shows the convergence of 𝑐, as the number of simulations (𝑚) increases.
We ran the program several times, first by changing 𝑛 and fixing 𝑚, then we fixed 𝑛 with changing 𝑚. We got our results for the convergence of the call option in the tables. Then we ran the program to the see if we will get to the theoretical order of convergence of the Runge–Kutta method.
We notice, in Table5.1, that the European call option price convergence to an approximate value 13.6933 when 𝑛 = 75. When we ran the program at n>75, we got close results as 13.6933, so we stopped at 𝑛 = 75 because it was the optimal approximated value to the theoretical one. In Table5.2, we estimate the call option price using different number of simulations. We notice that the call option also converges to 13.6933. In both cases, when 𝑚 and 𝑛 are changed, we got a result that is very close to the theoretical value of the call option that we mentioned and computed in the beginning of this chapter.
Number of intervals (𝑛) Call option value
15 13.7764 25 13.7345 35 13.7111 45 13.7001 55 13.7087 65 13.7201 75 13.6933
Table 5.1: Approximation the call option value using the Stochastic Runge–Kutta method as 𝑛 increasing and 𝑚 = 105.
Number of simulations (𝑚) Call option value 100000 13.7073 250000 13.7063 500000 13.6494 750000 13.6850 1000000 13.6933
Table 5.2: Approximation the call option value using the Stochastic Runge–Kutta method as 𝑚increasing and 𝑛 = 75.
Table5.1shows the European call option price as the number of simulations 𝑚 increases. We see that the result, from our program, of the call option price converges to the true value as 𝑚increases. We notice that we get close approximations of the call option when 𝑚 ≥ 7.6 · 105. When talking about the 95% upper and lower bound, we see that the estimate value lies between them. these bounds show us that we are 95% sure that the confidence interval contains the true value, which we need to estimate.
5.2
Order of Convergence of SRK method
Recalling the convergence for our method from Section3.4and Section4.3, and since we are talking about weak second order Runge–Kutta in our paper, then it has a global order 𝑂 (ℎ2), where 𝑝 = 2. In this section, we will investigate if the order of convergence of the SRK method for call option pricing is indeed the second-order.
Figure 5.1: Convergence of the SRK method as 𝑚 increases.
To be able to find the order of convergence in our method, we added — in our MATLAB script — the MATLAB realisation of Equation (4.1) to solve for 𝑝.
In this realisation, we fix the value of ℎ1from Equation (4.1) and run the program 6 times
with different values of ℎ2.
Applying the initial values presented in the beginning of this chapter, with fixed 𝑚 and ℎ1,
and varying 𝑛 and ℎ2 = 𝑇 /𝑛, after 6 attempts we get an approximated call option price and
an order of convergence equals to 1.88. We realize that this order of convergence tends to the theoretical value, 2, and we see these results in Table5.3and Figure5.2.
Exact call poption price Appr. price after 6 attempts Order of convergence
13.6933 13.7008 1.88
Table 5.3: Order of convergence of the SRK method.
To make sure that the order of convergence really tends to 2, we ran the program many times with various initial and parameter values. Table5.4 and Figure5.3 show the following results.
After looking to the results, in Table 5.4, we notice different values for the order of convergence, where some values are not close to the exact order of convergence. But by looking at Figure5.3, we see that the order of convergence also tends to the theoretical value, in other words get bigger, when ℎ2, the length of the interval, increases.
Figure 5.2: Order of convergence of the SRK method as ℎ2decreases.
changed data Call option price Order of convergence Exact value Appr. with ℎ2= 𝑇 /(6 × 75)
𝑠=150 57.6473 57.6632 0.922 𝑠 =90 7.8187 7.8187 0.895 𝐾 =105 8.9089 8.9113 1.216 𝐾 =85 20.0705 20.828 1.331 𝜎=0.75 18.2915 18.2945 1.394 𝜎=0.25 9.3390 9.3455 0.666
𝑠=150 𝑠 =90
𝑘 =105 𝑘 =85
𝜎=0.75 𝜎 =0.25
Chapter 6
Conclusion
In this paper, the main idea is to solve the Black–Scholes model with certain initial values in order to obtain a theoretical European call option which will be used to compare it with the obtained value of the SRK method. As we expected, after many changes in data and investigating, we got an approximated value to the true call option price as the number of simulations and number of intervals increase.
As mentioned before, we studied the weak second order Runge–Kutta in this thesis. From our results and according to the global order, we notice that the theoretical order of convergence is 2, and after computing our method we come to an approximated order, 1.88, which is very close to the theoretical one. Even when the data changed, the graphs still showed us a convergence in the order to the real one as ℎ2increases.
Bibliography
Nahum Biger and John Hull. The valuation of currency options. Financial Management, 12 (1):24–28, 1983. ISSN 00463892, 1755053X. URLhttp://www.jstor.org/stable/
3664834.
Fischer Black and Myron Scholes. The pricing of options and corporate liabilities. J. Polit.
Econ., 81(3):637–654, 1973. ISSN 0022-3808. doi: 10.1086/260062. URL https:
//doi.org/10.1086/260062.
Kevin Burrage and Pamela M. Burrage. High strong order explicit Runge–Kutta methods for stochastic ordinary differential equations. Appl. Numer. Math., 22(1-3):81–101, 1996. ISSN 0168-9274. doi: 10.1016/S0168-9274(96)00027-X. URL https://doi.org/10.
1016/S0168-9274(96)00027-X. Special issue celebrating the centenary of Runge–Kutta
methods.
John C. Butcher. Numerical methods for ordinary differential equations in the 20th century. J. Comput. Appl. Math., 125(1-2):1–29, 2000. ISSN 0377-0427. doi: 10.1016/S0377-0427(00)00455-6. URLhttps://doi.org/10.1016/S0377-0427(00)
00455-6. Numerical analysis 2000, Vol. VI, Ordinary differential equations and integral
equations.
John C. Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, Ltd., Chichester, third edition, 2016. ISBN 978-1-119-12150-3. doi: 10.1002/ 9781119121534. URL https://doi.org/10.1002/9781119121534. With a foreword by J. M. Sanz-Serna.
Kristian Debrabant and Andreas Rößler. Classification of stochastic Runge-Kutta methods for the weak approximation of stochastic differential equations. Math. Comput. Simulation, 77 (4):408–420, 2008. ISSN 0378-4754. doi: 10.1016/j.matcom.2007.04.016. URL https:
//doi.org/10.1016/j.matcom.2007.04.016.
Kristian Debrabant and Andreas Rößler. Families of efficient second order Runge-Kutta methods for the weak approximation of Itô stochastic differential equations. Appl. Numer.
Math., 59(3-4):582–594, 2009. ISSN 0168-9274. doi: 10.1016/j.apnum.2008.03.012. URL
Leonhard Euler. Institutionum calculi integralis ... impensis Academiae Imperialis Scien-tiarum, Petropoli, 1768. doi: 10.3931/e-rara-29427. URLhttps://doi.org/10.3931/
e-rara-29427.
Emil A. Fellmann. Leonhard Euler. Birkhäuser Verlag, Basel, 2007. ISBN 978-3-7643-7538-6; 3-7643-7538-8. Translated and with notes by Erika and Walter Gautschi.
Yuan Gao. Existence and uniqueness theorem on uncertain differential equations with local Lipschitz condition. Journal of Uncertain Systems, 6(3):223–232, 2012.
Mark B. Garman and Steven W. Kohlhagen. Foreign currency option values. Journal of
International Money and Finance, 2(3):231–237, 1983. URL https://EconPapers.
repec.org/RePEc:eee:jimfin:v:2:y:1983:i:3:p:231-237.
Karl Heun. Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen. Z. Math. Phys. I, 45:23–38, 1900.
Anton Huťa. Une amélioration de la méthode de Runge-Kutta-Nyström pour la résolution numérique des équations différentielles du premier ordre. Acta Fac. Natur. Univ. Comenian.
Math., 1:201–224, 1956. ISSN 0373-8183.
Anton Huťa. Contribution à la formule de sixième ordre dans la méthode de Runge–Kutta– Nyström. Acta Fac. Natur. Univ. Comenian. Math., 1:21–24, 1957. ISSN 0373-8183.
Liviu Gr. Ixaru and Guido Van Den Berghe. Exponential fitting, volume 568 of
Mathem-atics and its Applications. Kluwer Academic Publishers, Dordrecht, 2004. ISBN 1-4020-2099-6. doi: 10.1007/978-1-4020-2100-8. URL https://doi.org/10.1007/
978-1-4020-2100-8. With 1 CD-ROM (Windows, Macintosh and UNIX).
Peter E. Kloeden and Eckhard Platen. A survey of numerical methods for stochastic differential equations. Stochastic Hydrology and Hydraulics, 3(3):155–178, 1989. URL https://
EconPapers.repec.org/RePEc:uts:ppaper:1989-1.
Wilhelm Kutta. Beitrag zur näherungsweisen Integration totaler Differentialgleichungen. Teubner, 1901. URLhttps://books.google.se/books?id=K5e6kQEACAAJ.
Terence J. Lyons. Kiyoshi Itô (1915–2008). Probab. Theory Related Fields, 148(1-2):1–4, 2010. ISSN 0178-8051. doi: 10.1007/s00440-010-0286-7. URLhttps://doi.org/10.
1007/s00440-010-0286-7.
Rogemar S. Mamon. Three ways to solve for bond prices in the Vasicek model. J. Appl. Math.
Decis. Sci., 8(1):1–14, 2004. ISSN 1173-9126. doi: 10.1207/s15327612jamd0801_1. URL
https://doi.org/10.1207/s15327612jamd0801_1.
Robert C. Merton. Theory of rational option pricing. Bell J. Econom. and Management Sci., 4:141–183, 1973. ISSN 0005-8556.
Evert Johannes Nyström. Über die numerische Integration von Differentialgleichungen:
(Mit-geteilt am 23. Sept. 1925 von E. Lindelöf und K.F. Sundman). Societas scientiarum Fennica, 1925.
James Orlin Grabbe. The pricing of call and put options on foreign exchange. Journal of
International Money and Finance, 2(3):239–253, 1983. URL https://EconPapers.
repec.org/RePEc:eee:jimfin:v:2:y:1983:i:3:p:239-253.
Fabrice D. Rouah. The Heston Model and Its Extensions in MATLAB and C#. Wiley, September 2013. ISBN 1118548256.
Carl Runge. Über die numerische Auflösung von Differentialgleichungen. Zenodo, June 1895. doi: 10.1007/bf01446807. URLhttps://doi.org/10.1007/bf01446807.
Ruslan L. Stratonovich. A new representation for stochastic integrals and equations. SIAM J.
Control, 4:362–371, 1966. ISSN 0363-0129.
Ruslan L. Stratonovič. A new form of representing stochastic integrals andequations. Vestnik
Moskov. Univ. Ser. I Mat. Meh., 1964(1):3–12, 1964. ISSN 0201-7385.
Ángel Tocino and Ramón Ardanuy. Runge-Kutta methods for numerical solution of stochastic differential equations. J. Comput. Appl. Math., 138(2):219–241, 2002. ISSN 0377-0427. doi: 10.1016/S0377-0427(01)00380-6. URLhttps://doi.org/10.1016/S0377-0427(01)
00380-6.
Nico G. van Kampen. Stochastic processes in physics and chemistry, volume 888 of Lecture
Notes in Mathematics. North-Holland Publishing Co., Amsterdam-New York, 1981. ISBN 0-444-86200-5.
Robin Wilson. Leonhard Euler. Math. Intelligencer, 39(2):108, 2017. ISSN 0343-6993. doi: 10.1007/s00283-016-9646-1. URLhttps://doi.org/10.1007/s00283-016-9646-1.
Appendix A
MATLAB Simulation
A.1
Script 1
This program is built to calculate an approximate price of the European call option in the Black–Scholes model by a version of a stochastic Runge–Kutta method.
% Script : RungeKutta .m
% Authors : Ahmad Al - Kadri , Ali Saleh
% Aim : Calculation of an approximate price of the European call option
% in the Black - Scholes model by a version of a stochastic Runge - Kutta
% method
% Version : November 5, 2020
% Introduce constants and variables
% Constants , to denote by capital letters
C0 =[0 0.5 1];
C1 =[0 342 342]/491.0;
% C2 is not important , it is equal to 0
A0 =[ 0 0 0 0.5 0 0 -1 2 0]; A1 =[ 0 0 0 342 0 0 342 0 0]/491.0;
% A2 is not important , it is equal to 0
B0 =[ 0 0 0
6- sqrt (6) 0 0
B1 =[ 0 0 0 sqrt (38/491) 0 0 - sqrt (38/491) 0 0]*3; B2 =[0 0 0 -214* sqrt (1105/991) -491* sqrt (221/4955) -491* sqrt (221/4955) 214* sqrt (1105/991) 491* sqrt (221/4955) 491* sqrt (221/4955) ]/513; alphaT =[1 4 1]/6.0; bta1T =[386 491 491]/1368.0; bta2T =[0 sqrt (491/38) - sqrt (491/38) ]/6.0; bta3T =[ -9910 4955 4955]/14144.0; bta4T =[0 - sqrt (4955/221) sqrt (4955/221) ]/8.0; STAGES = 3;
% Variables , to denote by small letters % The initial stock price
s =100.0;
% the strike price
k =95.0;
% the interest rate
r =0.1;
% the volatility
sig =0.5;
% the maturity time
t =0.25;
% the number of intervals
n =155;
% the number of simulation
numSimulations =[100000:10000:1000000];
% Working variables
lowEstimate = zeros ( size ( numSimulations )); upEstimate = zeros ( size ( numSimulations )); estimate = zeros ( size ( numSimulations )); BSValue = ones ( size ( numSimulations )); [ Call , Put ]= blsprice (s ,k ,r ,t , sig ); disp ( Call );
BSValue = BSValue * Call ; control =1;
cumulativeValue =0; cumulativeSquare =0; h=t/n;
sqrth = sqrt (h);
value = sqrt (3) * sqrth ; discount = exp (-r*t);
w = waitbar (0 ,'Please wait...');
% A loop with respect to the number of simulations
for i =1: numSimulations ( end )
% the price of the security
y=s;
% a loop with respect to the number of intervals
for j =1: n
h0 = ones ( size (1: STAGES ))*y; h1 = ones ( size (1: STAGES ))*y; h1hat = ones ( size (1: STAGES ))*y;
% simulate hat I hatI =0; random = rand (1) ; if ( random <=1/6) hatI = value ; end if ( random >=5/6) hatI =- value ; end
% calculate stage values
for ii =1: STAGES for jj =1: ii -1 h0 ( ii )= h0 ( ii )+ A0 (ii , jj )*r* h0 ( jj )*h+ B0 (ii , jj )* sig * h1 ( jj )* hatI ; h1 ( ii )= h1 ( ii )+ A1 (ii , jj )*r* h0 ( jj )*h+ B1 (ii , jj )* sig * h1 ( jj )* sqrth ; end end
% Simulate the stock price
for ii =1: STAGES
y=y+ alphaT ( ii )*r* h0 ( ii )*h; y=y+ bta1T ( ii )* sig * h1 ( ii )* hatI ;
y=y+ bta2T ( ii )* sig * h1 ( ii ) *( hatI * hatI -h) /(2* sqrth ); y=y+ bta3T ( ii )* sig * h1hat ( ii )* hatI ;
y=y+ bta4T ( ii )* sig * h1hat ( ii )* sqrth ;
end end
optionValue = discount * max (y -k ,0) ;
cumulativeValue = cumulativeValue + optionValue ;
if (i == numSimulations ( control ))
% calculate intermediate results
estimate ( control )= cumulativeValue / numSimulations ( control );
variance = cumulativeSquare - numSimulations ( control )* estimate ( control )* estimate ( control );
variance = sqrt ( variance /( numSimulations ( control ) -1)); lowEstimate ( control )= estimate ( control ) -1.96* variance /
sqrt ( numSimulations ( control ));
upEstimate ( control )= estimate ( control ) +1.96* variance / sqrt ( numSimulations ( control ));
waitbar ( numSimulations ( control )/ numSimulations ( end ) ,w) ; control = control +1; end end close (w); % Show results figure ;
plot ( numSimulations , estimate , numSimulations , lowEstimate , numSimulations , upEstimate , numSimulations , BSValue ); title ('Monte Carlo simulations');
grid ;
legend ('Estimate','95% lower bound','95% upper bound','Exact
value');
xlabel ('Number of simulations'); ylabel ('Results');
disp ( estimate ( length ( numSimulations )));
A.2
Script 2
This program is built to calculate the order of convergence of our method.
% Script : RungeKutta .m
% Authors : Ahmad Al - Kadri , Ali Saleh
% Aim : Calculation of an approximate price of the European call option
% in the Black - Scholes model by a version of a stochastic Runge - Kutta
% method
% Introduce constants and variables
% Constants , to denote by capital letters
C0 =[0 0.5 1];
C1 =[0 342 342]/491.0;
% C2 is not important , it is equal to 0
A0 =[ 0 0 0 0.5 0 0 -1 2 0]; A1 =[ 0 0 0 342 0 0 342 0 0]/491.0;
% A2 is not important , it is equal to 0
B0 =[ 0 0 0 6- sqrt (6) 0 0 6+4* sqrt (6) 0 0]/10.0; B1 =[ 0 0 0 sqrt (38/491) 0 0 - sqrt (38/491) 0 0]*3; B2 =[0 0 0 -214* sqrt (1105/991) -491* sqrt (221/4955) -491* sqrt (221/4955) 214* sqrt (1105/991) 491* sqrt (221/4955) 491* sqrt (221/4955) ]/513; alphaT =[1 4 1]/6.0; bta1T =[386 491 491]/1368.0; bta2T =[0 sqrt (491/38) - sqrt (491/38) ]/6.0; bta3T =[ -9910 4955 4955]/14144.0; bta4T =[0 - sqrt (4955/221) sqrt (4955/221) ]/8.0; STAGES = 3;
% Variables , to denote by small letters % The initial stock price
s =100.0;
% the strike price
k =95.0;
% the interest rate
r =0.1;
% the volatility
sig =0.5;
% the maturity time
t =0.25;
at =7;
% order of convergence
order = zeros ( size (1: at ));
% absolute error of the approximation
error = zeros ( size (1: at ));
% the number of intervals
n =75:75:75* at ;
% The length of the interval
h= zeros ( size (1: at ));
% the number of simulation
numSimulations =[100000:10000:1000000];
% Working variables
lowEstimate = zeros ( size ( numSimulations )); upEstimate = zeros ( size ( numSimulations )); estimate = zeros ( size ( numSimulations )); BSValue = ones ( size ( numSimulations )); [ Call , Put ]= blsprice (s ,k ,r ,t , sig ); disp ( Call );
BSValue = BSValue * Call ; discount = exp (-r*t); for kk =1: at control =1; cumulativeValue =0; cumulativeSquare =0; h( kk )=t/n( kk ); sqrth = sqrt (h( kk )); value = sqrt (3) * sqrth ;
w = waitbar (0 ,'Please wait...');
% A loop with respect to the number of simulations
for i =1: numSimulations ( end )
% the price of the security
y=s;
% a loop with respect to the number of intervals
for j =1: n( kk )
h0 = ones ( size (1: STAGES ))*y; h1 = ones ( size (1: STAGES ))*y; h1hat = ones ( size (1: STAGES ))*y;
% simulate hat I hatI =0; random = rand (1) ; if ( random <=1/6) hatI = value ; end if ( random >=5/6)
hatI =- value ;
end
% calculate stage values
for ii =1: STAGES for jj =1: ii -1 h0 ( ii )= h0 ( ii )+ A0 (ii , jj )*r* h0 ( jj )*h( kk )+ B0 ( ii , jj )* sig * h1 ( jj )* hatI ; h1 ( ii )= h1 ( ii )+ A1 (ii , jj )*r* h0 ( jj )*h( kk )+ B1 ( ii , jj )* sig * h1 ( jj )* sqrth ; end end
% Simulate the stock price
for ii =1: STAGES
y=y+ alphaT ( ii )*r* h0 ( ii )*h( kk ); y=y+ bta1T ( ii )* sig * h1 ( ii )* hatI ;
y=y+ bta2T ( ii )* sig * h1 ( ii ) *( hatI * hatI -h( kk )) /(2* sqrth );
y=y+ bta3T ( ii )* sig * h1hat ( ii )* hatI ; y=y+ bta4T ( ii )* sig * h1hat ( ii )* sqrth ;
end end
optionValue = discount * max (y -k ,0) ;
cumulativeValue = cumulativeValue + optionValue ; cumulativeSquare = cumulativeSquare + optionValue *
optionValue ;
if (i == numSimulations ( control ))
% calculate intermediate results
estimate ( control )= cumulativeValue / numSimulations ( control );
variance = cumulativeSquare - numSimulations ( control )* estimate ( control )* estimate ( control );
variance = sqrt ( variance /( numSimulations ( control ) -1) );
lowEstimate ( control )= estimate ( control ) -1.96* variance / sqrt ( numSimulations ( control )); upEstimate ( control )= estimate ( control ) +1.96*
variance / sqrt ( numSimulations ( control ));
waitbar ( numSimulations ( control )/ numSimulations ( end ) ,w);
control = control +1;
end end
close (w);
% Show results
figure ;
plot ( numSimulations , estimate , numSimulations , lowEstimate , numSimulations , upEstimate , numSimulations , BSValue ); title ('Monte Carlo simulations');
grid ;
legend ('Estimate','95% lower bound','95% upper bound','
Exact value');
xlabel ('Number of simulations'); ylabel ('Results');
disp ( estimate ( length ( numSimulations )));
end
for kk =2: at
order ( kk ) =( log ( error ( kk )) -log ( error (1) )) /( log (h( kk ) -log (h (1) )));
end
figure ;
plot (h (2: at ) , order (2: at )); title ('Order of convergence'); grid ;
xlabel ('Length of h_2');
Appendix B
Required Objectives for the Thesis
This appendix shows that our work meets the Swedish Higher Education Agency’s requirements form getting a Bachelor degree.
B.1
Objective 1: Knowledge and understanding of the main
field
Since we read many articles and have studied some courses related to what we wrote about, we have shown that we understood and are familiar with the definitions and the given equations and formulas.
B.2
Objective 2: Ability to search, collect and interpret
In our thesis, we had to collect a lot of information about the stochastic Runge–Kutta, so we had to read many articles about numerical methods. the most article that helped us with our work and research wasDebrabant and Rößler(2009). We had almost all the information needed to understand the model. We analysed and interpreted our results in chapter 5, where we focused on the order of convergence and the convergence of the European call option price.
B.3
Objective 3: Identify, formulate and solve problems
We came across some mentionned equations in previous courses, in which we did not have enough information and knowledge about. Writing this thesis helped us to understand these equations and have all the knowledge we need to solve and compute them.
Our problem was to find an approximate a stochastic Runge–Kutta method to solve a Black– Scholes equation. We did what was required from us by solving and simulate our model, and then implementing it in MATLAB to get an approximated price of the call option.
B.4
Objective 4: Ability to present the report by writing and
orally, so others understand it
This objective is satisfied through the different sections and chapters presented in this paper. We showed the purpose of this thesis and its methodology followed by a chapter where we explained about numerical methods and the gave numerical solutions.
B.5
Objective 5: Ability to put the report to social context
Our thesis shows that the European call option price can also be calculated using other model, than the Black–Scholes, such as the stochastic Runge–Kutta method. This work is new, interesting and can be useful in the financial market.
