School of Education, Culture and Communication
Division of Applied Mathematics
BACHELOR THESIS IN MATHEMATICS / APPLIED MATHEMATICS —
MMA390
A second order Runge–Kutta method for the Gatheral model
by
Jérémy Auffredic
Kandidatarbete i matematik / tillämpad matematik MMA390
DIVISION OF APPLIED MATHEMATICS
MÄLARDALEN UNIVERSITY SE-721 23 VÄSTERÅS, SWEDEN
School of Education, Culture and Communication
Division of Applied Mathematics
Bachelor thesis in mathematics / applied mathematics MMA390
Date:
5 June 2020
Project name:
A second order Runge–Kutta method for the Gatheral model
Author(s): Jérémy Auffredic Version: 26th June 2020 Supervisor(s): Anatoliy Malyarenko Reviewer: Siyang Wang Examiner: Ying Ni Comprising: 15 ECTS credits
Abstract
In this thesis, our research focus on a weak second order stochastic Runge–Kutta method applied to a system of stochastic differential equations known as the Gatheral Model. We approximate numerical solutions to this system and investigate the rate of convergence of our method. Both call and put options are priced using Monte-Carlo simulation to investigate the order of convergence. The numerical results show that our method is consistent with the theoretical order of convergence of the Monte-Carlo simulation. However, in terms of the Runge-Kutta method, we cannot accept the consistency of our method with the theoretical order of convergence without fur-ther research.
Keywords: Stochastic differential equation, Runge–Kutta method, Gatheral model, Monte-Carlo simulation, weak second order.
Acknowledgements
First, I would like to thank my supervisor, Anatoliy Malyarenko for all his dedicated time during this six last months. Secondly, I am thankful for the comments and re-view of Dr. Siyang Wang on my work, and the improvement made from it. Finally, I am also grateful to have been one of the students at Mälardalen University.
Contents
1 Introduction 6
1.1 Background and literature review . . . 6
1.2 Research question . . . 9
1.3 Methodology . . . 9
1.4 Further outline . . . 9
2 The Runge–Kutta scheme 11 2.1 Runge–Kutta methods . . . 11
2.1.1 The case of deterministic differential equations . . . 11
2.1.2 The case of stochastic differential equations . . . 17
2.2 The scheme by Xiao Tang and Aiguo Xiao . . . 22
3 Numerical simulation 23 3.1 Presentation of the Gatheral model . . . 23
3.2 Scheme applied to the Gatheral model . . . 24
3.3 Option pricing . . . 25
3.4 Monte-Carlo simulation . . . 26
3.5 Implementation . . . 27
3.5.1 Preallocation . . . 27
3.5.2 Loop invariant variable . . . 28
3.5.3 Program description . . . 28
3.5.4 Pseudocode: option pricing as the number of simulations varies 29 3.5.5 Pseudocode: option pricing as the number of steps varies . . . 30
4 Numerical Results 31 4.1 Convergence of the results . . . 32
4.2 Convergence of Monte-Carlo . . . 32
4.3 Convergence of Runge–Kutta method . . . 36
5 Conclusion 39 A Algorithm 43 A.1 Simulation algorithm . . . 43
B Bachelor Degree Objectives 55 B.1 Objective 1: Knowledge and understanding of the main field . . . 55 B.2 Objective 2: Ability to search, collect and interpret informations . . . . 55 B.3 Objective 3: Identify, formulate and solve problems . . . 56 B.4 Objective 4: Ability to present orally and by writing, problems and
solutions to different groups . . . 56 B.5 Objective 5: Ability to make judgements in respect to scientific, societal
and ethnical aspects . . . 56
List of Figures
4.1 Convergence of Method (3.4) as m increase . . . 33 4.2 Evolution of the mean absolute error as m increase . . . 35 4.3 Evolution of the mean absolute error as h Ñ 0 . . . 38
List of Tables
2.1 Butcher’s Tableau . . . 17 4.1 Approximation of call and put option prices as m increase . . . 32 4.2 Approximation of call option price using Method (3.4) as m increase . 34 4.3 Approximation of put option price using Method (3.4) as m increase . 34 4.4 Approximation of the call option price using the Runge–Kutta method
on the Gatheral model . . . 36 4.5 Approximation of the put option price using the Runge–Kutta method
Chapter 1
Introduction
For hundreds of years, ordinary differential equations (ODEs) have been used to model phenomena in several domains such as natural science, physics or even econom-ics. Ordinary differential equations have been a topic of interest for many mathem-aticians, who proposed closed-form solutions to these problems and also numerical solutions in the latest years with the improvement of digital computers. Unfortu-nately, ODEs lack of realism when phenomena arise from stochastic effects, such as in finance with the study of stock price.
This is where stochastic differential equations (SDEs) have been further developed to avoid the idealized situation created by ODEs, leading to deterministic solutions, when in reality several millions of solutions could be obtained. In the latest years, stochastic differential equations have been popularized into many domains where stochastic is required such as finance, physics, biology, renewable energy and more. However, SDEs increased the complexity of the new models. Even if some of them have closed-form solutions, most of them do not and require highly complex numer-ical methods and powerful computers. One of the largest family of numernumer-ical meth-ods applied to SDEs, is called Stochastic Runge–Kutta methmeth-ods (SRK). This family of methods has been developed to avoid complex high-order derivatives used in nu-merical methods based on Taylor series.
In recent years, SDEs have been largely developed in the financial sector after the work of Fischer Black and Myron Scholes [2], who found a stochastic model to price options. Since then, several models have been developed using SDEs to price options, but the Black–Scholes model stays the most used one.
This thesis studies the order of convergence of a new numerical method applied on a pricing model. New pricing models are developed each year, continually seeking to improve the comprehension of the financial markets.
1.1
Background and literature review
There exist several numerical methods classified into groups such as explicit or impli-cit, one-step or multi-steps. These numerical methods have been first developed for
ODEs and later on for SDEs, necessary to approximate solutions of complex models. Let us consider an ordinary differential equation with given initial value by,
y1
ptq “ f pt, yptqq, ypt0q “ y0.
The approximation of solutions of ODEs with numerical methods is based on it-eration, starting from the given initial value. This means that the new approximate value will always be computed using only the value obtained from the previous step or also intermediate value within the step. These methods are called one-step meth-ods, because they require only one initial value for each new approximation.
The first numerical method to solve ordinary differential equations, is proposed by Leonhard Euler in 1768 [7], who create a method using the ODE itself to approx-imate solution near the given initial value. The method, although simple, can be sometimes unstable and inaccurate for large step-size between the successive approx-imations. Most numerical methods are characterized by their order of convergence [8]. This method is of order 1, leading to expensive computational effort if accurate results are desired. Later on, Karl Heun [10] improved Euler method to a second order two-stage method. He adds an intermediate approximation of the solution to compute the final approximation, leading to a better accuracy.
The methods construct with several stages based on Euler method are called Runge-Kutta methods. It has been attributed to Runge [25], who was the first one to general-ize the Euler method in 1900 and Kutta [18] for his contribution of methods of order 4 and 5. The generalization of the Euler method give the opportunity to research-ers to find methods with a higher order of convergence up to 8 and to approximate solutions of higher order differential equations. For example, Nyström [21] proposed some methods for first and second order differential equations, while Huta [12] intro-duced a sixth order Runge–Kutta method.
In opposition to ordinary differential equations which are deterministic, stochastic differential equations have one or several stochastic variables. There are two common ways for a differential equation to be stochastic. The first one is to have a stochastic initial value. The second one is to have a stochastic part added to the differential equation also called the diffusion coefficient or the noise, which is a random variable. The solution of a SDE is a stochastic process. We estimated its probability distribu-tion by sampling a large amount of soludistribu-tions and investigate the results, to find the expected value and variance.
Let us consider a stochastic differential equation with given initial value by, dyptq “ f pyptqqdt ` d ÿ j“1 gjpyptqqdWjptq, ypt0q “ y0, t P rt0, T s, (1.1)
where yptq P Rn, and f pyq and gjpyqare Rn-valued functions. The function f is the
deterministic component of the equation, called the drift coefficient, while g is the dif-fusion coefficient, and Wjptq, 1 ď j ď d, are independent Wiener processes.
An SDE can be interpreted into two different forms, Itô or Stratonovich [14], named after the two mathematicians who extended the methods of calculus to stochastic in-tegrals. In this thesis, we will focus on the Itô integral form, the most used one in applied mathematics.
With the rise of SDEs inside predictive models and the complexity to find explicit solutions, many mathematicians developed new numerical methods to approximate solutions of SDEs. Kloeden and Platen [14] showed that numerical methods for ODEs are incompatible and fail when applied to SDEs. So, numerical methods for SDEs have been required to be developed on their own. One of the constraints of all nu-merical methods for SDEs, is to approximate the solution of the diffusion coefficient. The most used method is to sample the Wiener process using pseudo-random num-bers generated by a computer. Komori [15] showed that how these random values are generated, can have an influence on the order of convergence of the method.
Euler–Maruyama method has been the first numerical method to approximate solutions of SDEs. It has been developed by Gisiro
¯ Maruyama [19], adapted from the standard Euler Method. If this numerical method is simple to apply, it is only of weak order 1 and restricted to SDEs of one dimensional Wiener process. In 1974, Mil0šte˘ın [20] developed a method similar to Euler–Maruyama method with a higher accuracy. Since then, several methods have been developed and extended to the class of Runge– Kutta methods for stochastic differential equations, such as the weak second order Runge-Kutta method (SRK) by Kloeden and Platen [14], Komori [16] and Tocino and Vigo-Aguiar [28].
All these methods are good approximation methods for SDEs, but restricted to small dimensional Wiener processes. When the dimension d of Wiener process is large, the methods are becoming too expensive because of the number of evaluations of diffusion coefficients and generation of random variables. Mathematicians such as Rößler [23, 24], Komori and Burrage [16, 17] worked on these issues and developed new efficient SRK methods of weak second order. Xiao Tang and Aiguo Xiao pro-posed recently an improved SRK method of weak order two based on the papers of Rößler, reducing the computational effort of the original method. In this thesis, our work will be focused around this last method.
In the recent years, the financial market has been a place where stochastic pro-cesses had an important role. Stochastic differential equations are present in all recent models used to price options and bonds or to estimate interest rates and volatility. The most famous model has been proposed by Black and Scholes [2] in 1973. However, newer models such as Heston Model [9], Vasicek model [29] or even the CIR model
[4] have also taken an important place. Most of these models have been developed to price financial derivatives. Hull [11] defined a financial derivative as a financial instru-ment whose value depends on the values of other underlying variables. One of the most used derivatives are options. There exist two types of options, call and put. A call option gives to the holder the right, but not the obligation, to buy the underlying asset at a certain date for a certain price. A put option, in opposition to the call, gives the right to sell the underlying asset. These types of financial derivatives are priced in a market that respects some assumptions, such as complete and arbitrage-free mar-ket. These assumptions create a strong link between the price of the asset and the price of the option, which is the existence of a probability measure P˚such that all the
discounted price processes are martingales under P˚. In this case, the no-arbitrage
price of the option is the expectation under P˚ of its discounted final payoff.
1.2
Research question
In this thesis, we wish to estimate the order of convergence of the stochastic Runge– Kutta method of weak second order proposed by Xiao Tang and Aiguo Xiao [27] applied to the Gatheral model [1]. The Gatheral model is a model usually used to price options on the SPX and VIX, respectively a stock index and a volatility index. We wish to analyze if this model can be used to price European stock options under the SRK method. Our main study, is to determine if the SRK method applied to the Gatheral model converges to accurate solutions in the same order of convergence as the theoretical value of the work of Xiao Tang and Aiguo Xiao.
1.3
Methodology
We first implement the coefficients of the explicit method with weak order two given in a Butcher tableau in [27] denoted W2Ito1, in the general s-stage Runge-Kutta method proposed by Xiao Tang and Aiguo Xiao. From the newly constructed Runge-Kutta method, applied to the system of equations of the Gatheral model, we compute an approximation to the solution. We use Monte Carlo simulation to compute the ex-pected value of the stock at maturity. We create two algorithms on MATLAB R2020, computing the numerical solutions to the Gatheral model for different numbers of simulations and steps. Finally, we study the convergence of our method, by compar-ing our numerical results to an estimated true value obtained with a large number of simulations and steps.
1.4
Further outline
The thesis is organized as follows. In Chapter 2, we present the mathematical back-ground of the Runge–Kutta methods in the case of deterministic and stochastic
dif-ferential equations. We also present in this chapter the new scheme by Xiao Tang and Aiguo Xiao improved from Rößler, followed by its scheme under the W2Ito1 explicit method. In Chapter 3, we describe the Gatheral model and our implementation of the SRK method. The two algorithms created in the MATLAB environment to test our method as the parameters vary are also included in Chapter 3. Our numerical results are described in Chapter 4, where we investigate the option prices depending on the values of different variables. Finally, we conclude in Chapter 5, where we re-view our results and give some explanation of unexpected results with some further possible research.
Chapter 2
The Runge–Kutta scheme
2.1
Runge–Kutta methods
The Runge–Kutta (RK) methods are numerical methods based on recursion to approx-imate solutions of Ordinary Differential Equations (ODEs) and Stochastic Differential Equations (SDEs). These methods have the advantages to be free of derivatives and to reduce the complexity of approximation. The theory provides that higher order of convergence leads to better approximation. The most used RK methods are Euler method, Heun’s method and the classical Runge–Kutta method also referred as RK4.
2.1.1
The case of deterministic differential equations
We define as deterministic, a system, an equation or a model which does not depend on random factors. We mean that all specific parameter inputs will always lead to the same outputs.
Definition 1. A first-order ordinary differential equation is given by y1
ptq “ f pt, yptqq, (2.1)
where the function f gives the slope of the tangent line of yptq at each point pt, yptqq. The solution of an ordinary differential equation is usually not unique. Ordinary differential equations can be of higher order leading to higher derivatives y2ptq, . . . ,
ypnqptq. We will focus uniquely on first order (2.1) in this thesis. The equation (2.1)
together with an initial condition such as ypt0q “ y0 is called an Initial Value Problem
(IVP). When considering the study of an initial value problem, two questions arise. First, if there exists a solution to the IVP and second if this solution is unique.
Definition 2. A function f defined on a domain D satisfies the Lipschitz condition, if there exists some constant L such that
|f pt, y2q ´ f pt, y1q| ď L|y2´ y1|
Theorem 1(Existence and Uniqueness for ODEs). Let an initial value problem be defined on the domain D, given by
y1ptq “ f pt, yq and ypt0q “ y0,
where f is continuous in t and satisfy a Lipschitz condition with respect to y. Then for any points of D, there is a neighborhood where there exists a unique solution yptq to the initial value problem.
A proof of Theorem 1 can be found in [3], but is out of the scope in this thesis. Before introducing some numerical methods to approximate solutions of initial value problems, we need to define three important concepts in numerical analysis. We will describe the discretization, the truncation error and the order of convergence of a numerical method.
Discretization
Most of the numerical methods are discretization methods, meaning that they use a set of discrete points to approximate the solution of the function yptq.
Definition 3. Consider an interval I “ rt0, T s. We create an ordered subset of I in N
intervals such as
t0 ă t1 ă ¨ ¨ ¨ ă tN “ T
where tn“ t0` nhfor n “ 0, 1, ¨ ¨ ¨ , N , and where h “ pT ´ t0q{N is the step size.
The subset is then used to approximate the values of yptq at t1, t2, ¨ ¨ ¨ denoted
by y1, y2, . . .. In numerical methods, the notation yptnq refers to the exact solution,
whereas ynrefers to the approximation of yptnq.The advantage of using discretization
is that it reduces the number of evaluations of the derivative to obtain an approx-imation. Unfortunately, discretization leads to approximation errors called local and global truncation errors.
Local and Global Truncation Errors
The discretization error or local truncation error results from the difference between the exact solution yptn`1qand the approximation yn`1at each discrete step and is defined
as
ln`1 “ yn`1´ yptn`1q.
An upper bound of the local truncation error can be approximated using the Taylor’s expansion with the Lagrange remainder [26].
The global truncation error is the accumulation of all local truncation errors of a method and is defined at the final step N by
Order of convergence
Most numerical methods are described by their order of convergence. The order of convergence defines the rate at which the approximation method tends to the exact solution as the step-size h Ó 0.
Definition 4. A numerical method has order p if the Taylor series for the exact solution yptn` hqand the numerical approximation yn`1coincide up to and including the term
hp, taking into account the localizing assumption that y
n “ yptnq. Then there exists a
constant K such that for small h,
|yptn` hq ´ yn`1| ď Khp`1.
In numerical methods, the order of convergence is defined by the notation Big-O. Definition 5. A function f phq is said to be big-O of a function gphq denoted Opgphqq as h Ñ 0, if it satisfies
|f phq| ď C|gphq| for some constant C in a neighbourhood of the point 0.
Numerical methods have a property that, if the local truncation error is Ophp`1
q then the global truncation error is Ophpq. Now that we have defined all necessary
con-cepts to analyze numerical methods, we will look at the simplest numerical method to approximate solution of an IVP, the Euler method.
The Euler Method
Definition 6. The Euler method [7] is described as follows:
yn`1“ yn` hf ptn, ynq, (2.2)
where h is the step-size of the discretization.
The method is based on recursive iteration, evaluating f at a point ptn, ynq to
ap-proximate the next value yn`1. The method applied to an interval rt0, T s, will generate
an approximation of discrete value yi for each ti ě t0, for i “ 1, 2, ¨ ¨ ¨ until reaching
an approximation of the solution at the point T .
The derivation of Euler’s method [30] is made from the fundamental theorem of calculus. Indeed, integrate Equation (2.1) with respect to t from tnto tn`1:
żtn`1 tn y1 ptq dt “ żtn`1 tn f pt, yptqq dt. The fundamental theorem of calculus gives
yptn`1q ´ yptnq “
żtn`1 tn
In the left hand side, replace the exact values yptn`1qand yptnqwith their
approxima-tions yn`1and yn. In the right hand side, apply the rectangle rule
żtn`1 tn
f pt, yptqq dt « hf ptn, yptnqq « hf ptn, ynq.
And finally we obtain
yn`1 « yn` hf ptn, ynq,
i.e. the Euler method.
Example 1. Let us approximate the value yp1.5q from the following IVP, y1
ptq “ 2y, yp1q “ 1,
using Euler method with step size h “ 0.1, knowing that the exact solution is yptq “ e2t´2.
From Euler method, on the interval r1, 1.5s and a step size h “ 0.1, we obtain the following, y1 “ 1 ` p0.1q2p1q “ 1.2, y2 “ 1.2 ` p0.1q2p1.2q “ 1.44, y3 “ 1.44 ` p0.1q2p1.44q “ 1.728, y4 “ 1.728 ` p0.1q2p1.728q “ 2.0736, y5 “ 2.0736 ` p0.1q2p2.0736q “ 2.48832, (2.3)
and the exact solution is yp1.5q “ e2p1.5q´2« 2.72.
Example 1 highlights the simplicity of the Euler method, but also its inaccuracy for large step size h. The method is Ophq, so the step-size need to be halved for the global truncation error to also be halved. The Euler method can become really expensive for small step size h. To improve the accuracy of the Euler method, Karl Heun proposed an improved Euler method [10] with a second order of convergence.
Heun’s Method or Improved Euler Method
Heun’s method is a two-stage method where the trapezoid rule replaces the rectangu-lar one and adds one more derivative evaluation at the point that we want to estimate with the help of a pre-approximation with Euler method.
Definition 7. Heun’s methodis given by
yn`1 “ yn` h 2pf ptn, ynq ` f ptn`1, ˆyn`1qq, ˆ yn`1 “ yn` hf ptn, ynq, hn “ tn`1´ tn. (2.4)
Heun improved Euler method with the evaluation of the derivative at two points ptn, ynqand ptn`1, ˆyn`1q, where ˆyn`1 is first approximated using Euler method. With
this method, Heun reached a local truncation error Oph3q and a global truncation
error Oph2q. A reduction in the step-size of 1
2 will decrease the global truncation error
of 14, allowing the Heun’s method to be twice more accurate than the Euler method.
The classical Runge–Kutta Method
The most used numerical method is the Runge–Kutta method of order 4 (RK4), for its accuracy, efficiency and stability. It is a four-stage method, where the derivative is evaluated at four different points, and then the value of yn`1 is approximated using
the weighted average.
Definition 8. The RK4 method is given by
yn`1“ yn` h 6pk1` 2k2` 2k3` k4q, where k1 “ f ptn, ynq, k2 “ f ptn` h 2, yn` k1 2 q k3 “ f ptn` h 2, yn` k2 2 q k4 “ f ptn` h, yn` k3q (2.5)
The method achieves order 4, so the local truncation error is Oph5
qand the global truncation error is Oph4q. All previous numerical methods are part of the Runge–
Kutta family because they satisfy the generalization of the Euler method developed by Runge.
Generalization of Euler method
After several works to improve the Euler method, Kutta formulated the general scheme of the Runge–Kutta methods [13], the generalization of the Euler method with s-stage. It is basically a weighted average of slopes at different points of the interval rtn, tn`1s.
This generalization has led to methods of order up to 8.
Definition 9. Let s, the number of stages be an integer, and the coefficients a21, a31, a32,
¨ ¨ ¨ , as1, as2, ¨ ¨ ¨ , as,s´1, b1, b2, ¨ ¨ ¨ , c2, ¨ ¨ ¨ , cssatisfying the Taylor polynomial coefficients
of degree s,řs
i“1bi “ 1and
ři´1
j“1aij “ cifor i “ 1, 2, ¨ ¨ ¨ , s.
yn`1 “ yn` hpb1k1` ¨ ¨ ¨ ` bsksq “ yn` h s ÿ i“1 biki, where k1 “ f ptn, ynq, k2 “ f ptn` c2h, yn` ha21k1q, k3 “ f ptn` c3h, yn` hpa31k1` a32k2qq, ¨ ¨ ¨ ks“ f ptn` csh, yn` hpas1k1` ¨ ¨ ¨ ` as,s´1ks´1qq, (2.6)
For a method to be consistent and convergent, the coefficients have to satisfy the order conditions [3] derived from
yn` h s ÿ i“1 biki “ yptnq ` hdyptnq ` h2 2!d 2ypt nq ` h3 3!d 3ypt nq ` ¨ ¨ ¨ ` hs s!d sypt nq,
where the Lagrange remainder has been ignored, considering h small enough.
Remark 1. A s-stage Runge–Kutta method is of order s only for methods up to order 4. Above the fourth order, the method requires more stages that the order of the method. A particular scheme for an explicit Runge–Kutta method satisfying the order con-dition is summarized in a so-called Butcher tableau defined as follows
c1 c2 a21 c3 a31 a32 .. . ... ... ... ... cs as1 as2 as3 ¨ ¨ ¨ as,s´1 b1 b2 b3 ¨ ¨ ¨ bs´1 bs
Example 2. The Butcher tableau of the previously seen methods are given in Table 2.1. It remains to discuss the question of stability of the introduced numerical methods. Following [3], consider a linear differential equation
y1
pxq “ qypxq.
Denote the product hq by z. In a single step of length h, the exact solution behaves as follows:
yptn`1q “ exppzqyptnq.
On the other hand, the approximate solution satisfies the following equation: yn`1“ Rpzqyn,
0 1
(a) Euler method
0 1 1 1{2 1{2 (b) Heun’s method 0 1{2 1{2 1{2 0 1{2 1 0 0 1 1{6 1{3 1{3 1{6 (c) RK4
Table 2.1: Butcher’s Tableau
where the function Rpzq depends on the involved numerical method. For example,
Rpzq “ $ ’ & ’ %
1 ` z, for Euler method
1 ` z ` z2{2, for Heun’s method 1 ` z ` z2{2 ` z3{6 ` z4{24, for RK4,
see [3, p. 101].
Definition 10. The stability region is the set
t z P C : |Rpzq| ď 1 u.
The sense of this definition is as follows. At any point in the stability region, the computed solution remains bounded after any number of time steps.
2.1.2
The case of stochastic differential equations
Stochastic differential equations [26] are defined by the combination of an ordinary differential equation with a stochastic process. This stochastic process also called white noise is usually the derivative of the Wiener process.
Definition 11. An initial value problem for a one-dimensional Itô stochastic differential equation is given by
dyptq “ f pyptqqdt ` gpyptqqdWt, ypt0q “ y0, t P rt0, T s, (2.7)
It can also be write under the Itô integral form yptq “ ypt0q ` żt t0 f pypsqqds ` żt t0 gpypsqqdWs. The integralştt
0f pypsqqdsis a Riemann integral while the integral
şt
t0gpypsqqdWs is
a stochastic integral. As mentioned in the introduction, we will study the stochastic equations in Itô calculus [22].
As seen in Theorem 1, an ODE has an existing unique solution if the Lipschitz condition is satisfied. In the case of SDEs, another condition has to be satisfied to confirm the existence of a unique solution called the linear growth condition.
Definition 12. A function f defined on a domain D P R satisfies the linear growth condition, if there exists a constant c P R such that
|f pt, yq| ď cp1 ` |y|q, for any points pt, yq P D.
Theorem 2(Existence and Uniqueness of solutions for SDEs). Let a stochastic differen-tial equation defined on the domain D be given by
dyptq “ f pyptqqdt ` gpyptqqdWt, ypt0q “ y0,
where f and g are functions defined and continuous for all t P rt0, T sand all pairs of points y1
and y2 P R satisfy both the Lipschitz condition and the linear growth condition. Then, there
exists a constant K such that
|f py1q ´ f py2q| ` |gpy1q ´ gpy2q| ď K|y1´ y2|,
and
|f pyq| ` |gpyq| ď Kp1 ` |y|q.
(2.8)
Then the SDE with initial value admits a unique solution yptq .
In numerical methods applied to SDEs, the convergence is either weak or strong. We focus uniquely on the weak convergence in our thesis.
Weak order of convergence Let CL
PpRd, Rq denote the space of all functions in R being L time continuously
differ-entiable with polynomial growth.
Definition 13. A stochastic numerical method of order p converges with weak order p to the exact solution yptnqif for each f P CP2pp`1qpRd, Rq, there exists a constant Cf and
δ0 ą 0such that
|Epf pynqq ´ Epf pyptnqqq| ď Cfhp, h P p0, δ0q, (2.9)
Here, Ep¨q represents the expected value. Numerical methods applied to SDEs have a property that if the global order is Ophp
q then the local order is Ophp`12q.
The weak order of convergence is beneficial because it allows the replacement of the Wiener Process dWt by other approximation probability distribution less costly and
Euler–Maruyama method
In the case of SDEs, the simplest method is Euler–Maruyama method, using the Euler method adapted to SDEs.
Definition 14. Euler–Maruyama method is given by
yn`1 “ yn` f ptn, ynqh ` gptn, ynq∆Wn,
ypt0q “ y0,
(2.10) where hn “ tn`1 ´ tn, and ∆Wn “ Wtn`1´ Wtn is normally distributed with mean 0
and variance h, denote N p0, hq.
The weak order of convergence of Euler–Maruyama method is 1. Like Euler method, this method is simple but inaccurate for large step-size. Heun’s method seen in Subsection 2.1.1 has also been adapted to the stochastic case.
Heun’s method
As in the case of ODEs, Heun’s method averages the value of f and g at the initial point tn and final point tn`1. An additional stage is computed to estimate ˆyn`1 with
Euler–Maruyama.
Definition 15. Heun’s method is given by yn`1 “ yn` h 2pf ptn, ynq ` f ptn`1, ˆyn`1qq ` 1 2pgptn, ynq ` gptn`1, ˆyn`1qq∆Wn, ˆ yn`1 “ yn` hf pynq ` gpynq∆Wn. (2.11)
Numerical methods for SDEs have the particularity to have two different orders of convergence, one for the deterministic part and another for the stochastic part, respectively denoted by pPd, Psq. Even if Heun’s method has the same Ps order of
convergence as the Euler–Maruyama method, it is more accurate because of higher convergence order of the deterministic part.
In this thesis, we only present these two numerical methods in case of one-di-mensional Wiener process. However, several other methods with higher order can be found in [14]. These numerical methods, even if accurate and efficient in the case of one-dimensional Wiener process, are not suitable for higher one. Rößler [5, 24] is one of the first to propose an efficient new class of SRK methods of weak order 2 suitable for high-dimensional Wiener process.
Definition 16. A d-dimensional stochastic differential equation system is given by
dyptq “ f pyptqqdt ` m ÿ j“1 gjpyptqqdWjptq, ypt0q “ y0, t P rt0, T s, (2.12) where f, g : Rd
Ñ Rd satisfy a Lipschitz condition and a linear growth condition, W ptqis an m-dimensional Wiener process.
General s-stage SRK method of weak approximation of SDEs with m-dimensional Wiener process by Rößler
The method developed by Rößler in [6] is given by
yn`1 “ yn` s ÿ i“1 αif ptn` c p0q i hn, H p0q i qhn` s ÿ i“1 m ÿ k“1 βip1qgkptn` c p0q i hn, H pkq i qˆIpkq ` s ÿ i“1 m ÿ k“1 βip2qgkptn` cp1qi hn, Hpkqi q ˆI pk,kq ? hn ` s ÿ i“1 m ÿ k“1 βip3qgkptn` cp2qi hn, ˆHpkqi qˆIpkq, ` s ÿ i“1 m ÿ k“1 βip4qgkptn` c p2q i hn, ˆH pkq i q a hn, Hp0qi “ Yn` i´1 ÿ j“1 Ap0qij f ptn` c p0q j hn, H p0q j qhn` i´1 ÿ j“1 m ÿ l“1 Bp0qij glptn` c p1q j hn, H plq j q ˆIplq, Hpkqi “ Yn` i´1 ÿ j“1 Ap1qij f ptn` cp0qj hn, Hp0qj qhn` i´1 ÿ j“1 Bp1qij gkptn` cp1qj hn, Hpkqj q a hn, ˆ Hpkqi “ Yn` i´1 ÿ j“1 Ap2qij f ptn` cp0qj hn, Hp0qj qhn` s ÿ j“1 m ÿ l“1, l‰k Bp2qij glptn` cp1qj hn, Hplqj q ˆ Ipk,lq ? hn , (2.13) where i “ 1, ¨ ¨ ¨ , s, k “ 1, ¨ ¨ ¨ , m, ypt0q “ y0 and α, βp1q, ¨ ¨ ¨ , βp4q, cpqq, Apqq, Bpqq, for
q “ 0, 1, 2are the vectors and matrices of coefficients of the SRK method. The method is summarized in the following Butcher’s tableau,
cp0q Ap0q Bp0q
cp1q Ap2q Bp1q
cp2q Ap2q Bp2q
αJ βp1qJ βp2qJ
βp3qJ βp4qJ
In the SRK method, ˆIk, ˆIpk,lqare random variables defined by
P p ˆIpkq “ ˘ a 3hnq “ 1 6, P p ˆIpkq“ 0q “ 2 3, (2.14) ˆ Ipk,lq “ $ ’ ’ ’ ’ ’ & ’ ’ ’ ’ ’ % 1 2p ˆIpkq ˆ Iplq´ a hnI˜pkqq, k ă l, 1 2p ˆIpkqIˆplq` a hnI˜plqq, k ą l, 1 2p ˆI 2 pkq´ hnq, k “ l, k, l P t1, 2, ¨ ¨ ¨ , mu (2.15)
and P p ˜Ipkq“ ˘
?
hnq “ 12.
A weak first order method can be developed from (2.13) if the coefficients satisfy 9 order conditions, and a weak second order if they satisfy 59 order conditions. The order conditions are out of the scope but can be found in [5]. The new method (2.13) if suitable for high dimensional Wiener process, has a large number of order condi-tions to be satisfied and a number of simulated random variables per step of 2m ´ 1. Xiao Tang and Aiguo Xiao [27] based on the method of Rößler proposed a new weak second order method with only 15 order conditions and m ` 2 simulated random variables per step.
New weak second order SRK method by Xiao Tang and Aiguo Xiao The new SRK method is given by
yn`1 “ yn` s ÿ i“1 αif pH p0q i qh ` s ÿ i“1 m ÿ k“1 βip0qgkpHpkqi q ˆIk` s ÿ i“1 m ÿ k“1 βip1qgkpHpkqi q ˆIpk,kq Hp0qi “ yn` s ÿ j“1 ap0qij f pHp0qj qh ` s ÿ i“1 m ÿ k“1 bp0qij gkpHpkqj q ˆIk Hpkqi “ yn` s ÿ j“1 ap1qij f pHp0qj qh ` s ÿ i“1 bp1qij gkpHpkqj qξ ` s ÿ i“1 m ÿ l“1, l‰k bp2qij glpHplqi q ˆIpk,lq, i “ 1, 2, ¨ ¨ ¨ , s, k “ 1, 2, ¨ ¨ ¨ , m, (2.16)
where h is the stepsize, n “ 0, 1, ¨ ¨ ¨ , N ´1, y0 “ ypt0q, ˆIk, ξ, ˆIpk,mqare random variables
defined by P p ˆIk“ ˘ ? 3hq “ 1 6, P p ˆIk “ 0q “ 2 3, ξ “ η1 ? h, (2.17) ˆ Ipk,lq “ $ ’ ’ ’ ’ ’ ’ & ’ ’ ’ ’ ’ ’ % 1 2p ˆIl´ η2 ˆ Ilq, k ă l, 1 2p ˆIl´ η2 ˆ Ilq, k ą l, 1 2p ˆ I2 l ξ ´ ξq, k “ l, k, l P t1, 2, ¨ ¨ ¨ , mu (2.18)
and η1, η2are random variables defined by P pηi “ ˘1q “ 12, i P t1, 2u.
The method is summarized in the following Butcher’s tableau, Ap0q Bp0q
Ap1q Bp1q Bp2q
αJ βp0qJ
Xiao Tang and Aiguo Xiao proposed two suitable explicit methods coefficients of weak second order. The coefficients used in our thesis, denoted by W 2Ito1, are present in the following Butcher’s tableau,
0 0 1 2 0 6´?6 10 0 ´1 2 0 3`2 ? 6 5 0 0 0 0 0 1 4 0 1 2 0 1 0 1 4 0 0 ´ 1 2 0 0 0 0 0 1 6 2 3 1 6 ´1 1 1 2 0 ´2
2.2
The scheme by Xiao Tang and Aiguo Xiao
The explicit method of weak second order W 2Ito1 when implement in (2.16), is given by yn`1 “ yn` 1 6f pH p0q 1 qh ` 2 3f pH p0q 2 qh ` 1 6f pH p0q 3 qh ´ 3 ÿ k“1 gkpHpkq1 q ˆIk` 3 ÿ k“1 gkpHpkq2 q ˆIk ` 3 ÿ k“1 gkpHpkq3 q ˆIk` 2 3 ÿ k“1 gkpHpkq1 q ˆIpk,kq´ 2 3 ÿ k“1 gkpHpkq3 q ˆIpk,kq, Hp0q1 “Hp1q1 “ H p2q 1 “ H p3q 1 “ yn, Hp0q2 “yn` 6 ´?6 10 ´ g1pHp1q1 q ˆI1` g2pH p2q 1 q ˆI2` g3pH p3q 1 q ˆI3 ¯ ` 1 2f pH p0q 1 qh, Hp0q3 “yn´ hf pH p0q 1 q ` 2f pH p0q 2 qh ` 3 ` 2?6 5 ´ g1pHp1q1 q ˆI1` g2pH p2q 1 q ˆI2` g3pH p3q 1 q ˆI3 ¯ , Hp1q2 “yn` 1 4f pH p0q 1 qh ` 1 2g 1 pHp1q1 qξ ` g2pHp2q2 q ˆIp1,2q` g3pH p3q 2 qq ˆIp1,3q, Hp2q2 “yn` 1 4f pH p0q 1 qh ` 1 2g 2 pHp2q1 qξ ` g1pHp1q2 q ˆIp2,1q` g3pH p3q 2 q ˆIp2,3q, Hp3q2 “yn` 1 4f pH p0q 1 qh ` 1 2g 3 pHp3q1 qξ ` g1pHp1q2 q ˆIp3,1q` g2pH p2q 2 q ˆIp3,2q, Hp1q3 “yn` 1 4f pH p0q 1 qh ´ 1 2g 1 pHp1q1 qξ, Hp2q3 “yn` 1 4f pH p0q 1 qh ´ 1 2g 2 pHp2q1 qξ, Hp3q3 “yn` 1 4f pH p0q 1 qh ´ 1 2g 3 pHp3q1 qξ, i “ 1, 2, 3, k “ 1, 2, 3. (2.19)
Chapter 3
Numerical simulation
In this section, we present the Gatheral model as well as the implementation of method (2.19) to this model. The obtained method gives us the possibility to estimate the solu-tion to the Gatheral model on a desired interval, which can be used to price both call and put options.
3.1
Presentation of the Gatheral model
The Gatheral model is a three factor asset pricing model introduced in 2008 by Jim Gatheral in the Fifth world congress of the Bachelier Finance Society. The model, also called Double-Mean-Reverting model (DMR), has been the result of research on a new model to price options on the SPX, VIX and realized variance where the Black-Scholes model was inconsistent. In our project, we use this model to price equity stock options.
The Gatheral model is given by, dSt“ ? νtStdWt1, dνt“ κ1pνt1´ νtqdt ` ξ1νtα1dW 2 t, dν1 t“ κ2pθ ´ νt1qdt ` ξ2νt1α2dWt3, (3.1)
where St is the price of the asset, νt is the volatility and νt1 is the volatility of the
volatility.
A calibration of the different parameters κ1, κ2, θ, α1, α2, ξ1, ξ2 has to be done, to
simulate solutions to the Gatheral model. The derivation of the calibration is slow and complex, so we decide to not include it in our thesis. However, we want to inform the reader of the signification of each parameter, κ1, κ2are the mean reversion
rate, ξ1, ξ2 are the parameters of the variance processes, θ is the long-run volatility
and α1, α2 control how the volatility of volatility changes with the volatility level.
Some good estimates can be found in [1]. In our research, we will use the following estimates for all simulations.
θ “ 0.078, κ1 “ 5.50, κ2 “ 0.10, ξ1 “ 2.6, ξ2 “ 0.45, α1 “ 0.94, α2 “ 0.94,
3.2
Scheme applied to the Gatheral model
We rewrite the system of stochastic differential equation (3.1) with yptq “ pSt, νt, νt1qJ
in term of (2.12), we obtain d » — – St νt ν1 t fi ffi fl “ » — – 0 κ1pνt1´ νtq κ2pθ ´ νt1q fi ffi fldt ` » — – ? νtSt 0 0 fi ffi fldW 1 t ` » — – 0 ξ1νtα1 0 fi ffi fldW 2 t ` » — – 0 0 ξ2νt1α2 fi ffi fldW 3 t. (3.3)
The method (2.19) combined with (3.3) becomes
Sn`1 “ Sn´ ? νnpH p1q 1 q1Iˆ1` ? νnpH p1q 2 q1Iˆ1` ? νnpH p1q 3 q1Iˆ1 ` 2?νnpH p1q 1 q1Iˆp1,1q´ 2 ? νnpH p1q 3 q1Iˆp1,1q, νn`1 “ νn` 1 6 ´ κ1pνn1 ´ pH p0q 1 q2q ¯ h `2 3 ´ κ1pνn1 ´ pH p0q 2 q2q ¯ h ` 1 6 ´ κ1pνn1 ´ pH p0q 3 q2q ¯ h ´ ξ1pH p2q 1 q α1 2 Iˆ2` ξ1pH p2q 2 q α1 2 Iˆ2` ξ1pH p2q 3 q α1 2 Iˆ2` 2ξ1pH p2q 1 q α1 2 Iˆp2,2q´ 2ξ1pH p2q 3 q α1 2 Iˆp2,2q ν1 n`1 “ ν 1 n` 1 6 ´ κ2pθ ´ pH p0q 1 q3q ¯ h ` 2 3 ´ κ2pθ ´ pH p0q 2 q3q ¯ h ` 1 6 ´ κ2pθ ´ pH p0q 3 q3q ¯ h ´ ξ2pH1p3qq α2 3 Iˆ3` ξ2pH2p3qq α2 3 Iˆ3` ξ2pH3p3qq α2 3 Iˆ3` 2ξ2pH1p3qq α2 3 Iˆp3,3q´ 2ξ2pH3p3qq α2 3 Iˆp3,3q Hp0q1 “Hp1q1 “ H p2q 1 “ H p3q 1 “ pSn, νn, νn1q J, Hp0q2 “pSn, νn, νn1qJ` 6 ´?6 10 p ? νnSnIˆ1, ξ1νnα1Iˆ2, ξ2νn1α2Iˆ3qJ ` 1 2p0, κ1pν 1 n´ νnq, κ2pθ ´ νn1qqJh, Hp0q3 “pSn, νn, νn1q J ´ hp0, κ1pνn1 ´ νnq, κ2pθ ´ νn1qq J ` 2p0, κ1pνn1 ´ pH p0q 2 q2qq, κ2pθ ´ pH p0q 2 q3qqJh ` 3 ` 2 ? 6 5 p ? νnSnIˆ1, ξ1νnα1Iˆ2, ξ2νn1α2Iˆ3qJ, Hp1q2 “pSn, νn, νn1q J ` 1 4p0, κ1pν 1 n´ νnq, κ2pθ ´ νn1qq Jh ` p1 2 ? νnSnξ, ξ1pH p2q 2 q α1 2 Iˆp1,2q, ξ2pH p3q 2 q α2 3 Iˆp1,3qqJ,
Hp2q2 “pSn, νn, νn1qJ` 1 4p0, κ1pν 1 n´ νnq, κ2pθ ´ νn1qqJh ` p?νnpH p1q 2 q1Iˆp2,1q, 1 2ξ1ν α1 n ξ, ξ2pH p3q 2 q α2 3 Iˆp2,3qqJ, Hp3q2 “pSn, νn, νn1qJ` 1 4p0, κ1pν 1 n´ νnq, κ2pθ ´ νn1qqJh ` p?νnpH p1q 2 q1Iˆp3,1q, ξ1pH p2q 2 q α1 2 Iˆp3,2q, 1 2ξ2ν 1α2 n ξqJ, Hp1q3 “pSn, νn, νn1qJ` 1 4p0, κ1pν 1 n´ νnq, κ2pθ ´ νn1qqJh ´ p1 2 ? νnSnξ, 0, 0qJ, Hp2q3 “pSn, νn, νn1qJ` 1 4p0, κ1pν 1 n´ νnq, κ2pθ ´ νn1qqJh ´ p0,1 2ξ1ν α1 n ξ, 0qJ, Hp3q3 “pSn, νn, νn1qJ` 1 4p0, κ1pν 1 n´ νnq, κ2pθ ´ νn1qqJh ´ p0, 0,1 2ξ2ν 1α2 n ξq J , (3.4)
This final version of the method allows us to compute from an initial condition, ypt0q “ pS0, ν0, ν01qJ, the solution to the Gatheral model for a specific time t.
3.3
Option pricing
Let us recall that a European option is a financial instrument giving the right, but not the obligation, to the holder to buy or sell an asset for a certain price at a specific time. An option is characterized by the following parameters [11]:
• The price of the asset, denoted by S.
• The agreed future price also called the strike price, denoted by K.
• The time where the contract takes place also called the maturity date, denoted by T .
• The risk-free interest rate, denoted by r. • The volatility of the asset price, denoted by σ.
As mentioned in the introduction, if the financial market is complete and arbitrage-free, then there exists a probability measure P˚such that the value of a European call
option is given by
c “ e´rTE˚
rmaxtST ´ K, 0us
and a put option by
p “ e´rTE˚
3.4
Monte-Carlo simulation
Consider the interval I “ rt0, T swith uniform discretization
t0 ă t1 ă ¨ ¨ ¨ ă tN “ T
where h “ T ´t0
N is the step-size.
We approximate solutions to the system (3.3) with (3.4) on the interval I using the discretization. As seen in Chapter 2, as the step-size h of the discretization decreases, the approximation of solution at time T tends to the exact solution of the SDE.
Let us denote c and p, the exact prices of the call and put option and ˆcand ˆp, the approximations. By the Monte-Carlo estimate,
ˆ cm “ 1 m m ÿ i“1 e´rTErmaxtS Tpwiq ´ K, 0us, and ˆ pm “ 1 m m ÿ i“1 e´rTErmaxtK ´ S Tpwiq, 0us (3.5)
where m is the number of simulated paths wi for 1 ď i ď m of ST. We know that by
the strong law of large numbers, ˆ
cm Ñ c and pˆm Ñ p with probability 1 as m Ñ 8.
In our thesis, we decided to use the same methodology used by Kloeden and Platen [14], by arranging our simulations into L batches, each of them composed of mindependent simulations of ST. By using this methodology, we estimate the mean
option price of each batch by
ˆ µcj “ 1 m m ÿ i“1 e´rTErmaxtS T,i,j ´ K, 0us, ˆ µpj “ 1 m m ÿ i“1 e´rTErmaxtK ´ S T ,i,j, 0us (3.6)
where ST,i,j is the stock price at time T for the ith simulation of the jth batch, and
j “ 1, 2, ¨ ¨ ¨ , L.
From (3.6), we estimate the mean option price for both call and put option by
ˆ µ “ 1 L L ÿ j“1 ˆ µj,
and the variance σ2 µby ˆ σ2µ“ 1 L ´ 1 L ÿ j“1 pˆµj ´ ˆµq2.
By the central limit theorem, for L ˆ m ě 30, ˆµis assumed to be standard normal distributed, giving us the possibility to compute a confidence interval of ˆµ with a confidence level of 100p1 ´ αq by » –µ ´ zˆ 1´α 2 d ˆ σµ2 L ; ˆµ ` z1´ α 2 d ˆ σµ2 L fi fl,
where µ lies inside this interval with probability 1´α and z is the upper 1´α
2 quantile
of the standard normal distribution. Let ∆ “ z1´α2
b
ˆ σµ2
L , we obtained the confidence interval
rˆµ ´ ∆, ˆµ ` ∆s.
It has been proved following the Monte Carlo simulation that |ˆµ ´ µ|decreases at the rate of pLmq´12 as m Ñ 8.
The Gatheral model has not any closed-form explicit solution, so we cannot com-pute the exact error of our method. However, we can run our method for a large number of simulations and batches, and consider the estimated call and put option prices as true values. We compare our results for several numbers of simulations and estimate an order of accuracy of our method. The estimator
b
ˆ σµ2
L can also be
a good indicator of convergence of our method. If the estimator decreases, then the confidence interval will also decrease showing us more accurate results. We want to analyze, how our method results in price accuracy as the number of simulations or steps are changed.
3.5
Implementation
In this section, we describe the implementation of our method using the software MATLAB R2020. We created two similar scripts to compute the price of a call and put European option. The first script investigates approximation of mean, variance and the confidence interval of the call and put option prices as the number of simulations varies and the second script investigates the same estimators as the number of steps increases.
The computational cost of our scripts can be really expensive for a large number of batches, simulations and steps. We construct our algorithms with the use of prealloc-ation space and use of loop invariant variables to optimize the computprealloc-ational time.
3.5.1
Preallocation
The preallocation space is the fact to create a variable memory size before entering a loop where this variable is computed and stored. It will save memory in the system dedicate for this variable during the full time that the program is running. The preal-location will be made one time and then only the storage will be effective each loop. The use of preallocation reduce the computational time.
3.5.2
Loop invariant variable
A loop invariant variable has the property to have the same value before, during and after a loop. This variable, if present inside a loop will be computed each cycle, which increases the computational expense of the program. Such variables can simply be computed, stored before the loop and used as constants in the loop. For small scripts, the computational time without loop invariant variables will not be significant, but for larger scripts like ours, with four and three loops, a variable can be computed in the main inner loop at most,
Number of step value test ˆ LmN.
3.5.3
Program description
Our scripts are both constructed following Method (3.4) to compute stock prices from t0 to T . They are set with the same option parameters and the same Gatheral model
parameters. The main difference between the two scripts is the fact that one runs for various choices of step sizes and the other one for various numbers of simulations. In terms of random variables, we use the built-in function “rand” and “randi” in MATLAB to generate random variables.
3.5.4
Pseudocode: option pricing as the number of simulations
var-ies
Data:
- Define N - number of steps; - Define y - number of batches;
- Define m - number of simulations as a discrete ordered set; - Define initial values of the option:
S0 - Current stock price, σ - Volatility, σσ - Volatility of volatility, T - Maturity
time, r - Interest rate, K - Strike price; - Define Gatheral model parameters ;
- Define invariant variables and preallocation variables; for a “ 1 to y do
p “ 1;
for b “ 1 to lengthpmq do for i “ 1 to N do
Generate random variables η1, η2;
Compute ξ;
Generate random variables ˆIifor i “ 1, 2, 3;
Compute ˆIpk,lqfor k, l “ 1, 2, 3;
Generate Si`1, νi`1, νi`11 ;
end
Compute the call and put option payoff;
Compute the cumulative sum of the call and put option payoff; if b = m(p) then
Compute and store the discounted expected call and put option payoff ;
p “ p ` 1; end
Store the discounted expected payoff in a matrix where the rows are the batches and columns the simulations
end end
Compute the mean, the variance and the confidence interval of each batch obtained from the same number of simulations;
Compute the absolute error and ratio of error of the method for the specific number of simulations
3.5.5
Pseudocode: option pricing as the number of steps varies
Data:
- Define N - number of steps as a discrete set; - Define y - number of batches;
- Define m - number of simulations; - Define the initial values for the option:
S0 - Current stock price, σ - Volatility, σσ - Volatility of volatility, T - Maturity
time, r - Interest rate, K - Strike price; - Define Gatheral model parameters ;
- Define invariant variables and preallocation variables; for step “ 1 to lengthptestq do
- Define N = test(step); - Define h = T/N;
- Define invariants variables dependent of h; for a “ 1 to y do
for b “ 1 to m do for i “ 1 to N do
- Generate random variables η1, η2;
- Compute ξ;
- Generate random variables ˆIifor i “ 1, 2, 3;
- Compute ˆIpk,lq for k, l “ 1, 2, 3;
- Generate Si`1, νi`1, νi`11 ;
end
- Compute the call and put option payoff;
- Compute the cumulative sum of the call and put option payoff; end
- Compute and store the discounted expected call and put option payoff
end
- Store the discounted expected payoff in a matrix where the rows are the batches and columns the number of steps
end
- Compute the mean, the variance and the confidence interval of each batch obtained from the same number of steps;
- Compute the absolute error and the ratio error for the specific number of steps;
Chapter 4
Numerical Results
In this section, we investigate the convergence of our method in terms of Monte-Carlo simulation and Runge-Kutta method. We approximate several prices of a call and put option using our constructed method and algorithms as m and N increase. Then, we compute the absolute errors of these estimations to investigate the convergence order. All the results are presented in the form of tables and the convergence errors are also presented in graphs for both algorithms approximation errors. Our interest is to observe the change in price as the parameters vary. We decide to test our method with the following option parameters and initial values:
• St0 “ 100 • K “ 105 • νt0 “ 0.21 • ν1 t0 “ 0.11 • T “ 1 • r “ 0.05
As described in Section 3.1, the estimated parameters (3.2) of the Gatheral model will remain constant during all the process. In the numerical simulations, the main variables are the number of batches L, the number of simulations m and the number of steps N . However, in this thesis, we do not investigate the influence of the number of batches, which remains constant, L “ 5.
By lack of analytical solution of the Gatheral model, we cannot have the exact call and put option prices to compute the errors of our method. The true prices are ne-cessary to investigate the convergence of our method. We decide to run our program once, with high value parameters m and N and to consider the approximate values as exact prices of the call and put option. From our method, with m “ 106 simulations
and N “ 100 steps, we obtain the following results:
c “ 11.0773$ and p “ 15.8128$. (4.1) All simulations have been done with a CPU Intel I5-9400, with the operating sys-tem Windows 10 on MATLAB R2020a.
4.1
Convergence of the results
Let us recall the theoretical convergence of our method from Chapters 2 and 3. The Runge–Kutta method used in our method is Oph2
qand the Monte-Carlo simulation is OppLmq´12q. From our method, we want to analyze in Section 4.2 and 4.3 if our
method converges with the same order of convergence as m Ñ 8 and N Ñ 8. We estimate the error of our method by the mean absolute error given by,
c “ |ˆµc´ c| and p “ |ˆµp´ p|.
4.2
Convergence of Monte-Carlo
In this section, we investigate the call and put option prices in function of the num-ber of simulations. By the Monte-Carlo theory, our approximations should converge to the true values as m get large. If our method is well constructed each increase of the number of simulations by 4, should halve the mean absolute error. However, our method can be really expensive for large values m. With the available computational performance, we have limited our analysis to L “ 5 batches, each of them composed of m simulations from 100 to 409600 and a fix number of steps, N “ 252. These values could be increased, leading to better approximations if more computational power would be available. In the frame of this thesis, the values considered, are large enough to study the convergence of our method.
Our first results presented in Table 4.1, give us an idea of the convergence of the Monte-Carlo simulation. We approximate call and put option prices for different numbers of simulations, for three distinctive experiments denoted by ex.i, i = 1, 2, 3. We observed in both cases, the call and put option, that the price converges to a common value as m increases. The call option price converges to approximately 11 and the put option price to approximately 15.80. These results are close to the true values that we compute in the introduction of Section 4.
Number of simulation (m)
Call option price Put option price ex. 1 ex. 2 ex. 3 ex. 1 ex. 2 ex. 3 100 11.3805 13.2599 10.0109 14.9313 15.4066 18.7657 1000 11.5745 11.1292 10.2743 17.7272 14.6924 16.1224 10000 11.3176 11.0150 10.7380 15.9536 15.6682 15.9784 100000 10.9749 11.0527 11.0103 15.8806 15.7618 15.8120
Table 4.1: Approximation of call and put option prices as m increase
Figure 4.1 shows the price of the call option as the number of simulations increases for eight experiments. We observe that each experiment converges to a common
value. The same observation is made with the put option price, we omit the figures here. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Number of simulations 104 9 10 11 12 13 14 15 16 Price
Convergence of the call option price by Monte Carlo simulation
Figure 4.1: Convergence of Method (3.4) as m increase
In Tables 4.2 and 4.3, we investigate our method for quadruple numbers of simu-lations starting with m “ 100 to 409600. We present our results with the following es-timates: the mean price of the option, the variance, a 95% confidence interval around the mean price and the mean absolute error.
We observe that the results from the first algorithm converge in terms of call and put option prices to the true values, c and p as m increases. Our method achieves reasonable approximations only for m ě 4 ¨ 105. The standard deviations of the call
and put option prices are reduced respectively to 10´3 and 10´2. The 95% confidence
intervals give us good intervals for the mean prices, of less than 0.02. We also ob-serve in both cases, a convergence of the mean absolute errors to 0. However, we notice one inconsistency in the mean absolute error of the put option price where it increases at 102400 simulations to further decreases again. This can be explained
by the fact that because the process is random, the mean absolute error is also ran-dom. Even if the law of large numbers in Section 3.4 confirmed that the estimates ˆ
cm Ñ c,and ˆpm Ñ pas m Ñ 8, the estimates converge in an erratic rate which can
result in some irregular behavior. We give an idea of the convergence of our method by computing the ratios of the mean absolute errors. Theoretically, the ratios should all be equal to 2, however, we observe that it is not exactly the case of our method. For the call option ratios, we see an oscillation of the estimates close to 2 which could maybe converge closer to 2 if the number of simulations could be increased even more. Unfortunately, we cannot draw any conclusion from the ratios of the put op-tion which are too irregular to see any convergence. Later on, we will investigate the log-error in function of the log-simulation to obtain the exact order of our method.
Number of
simulation (m)
µ
ˆ
cσ
ˆ
µˆ[ˆ
µ
c´ ∆; ˆ
µ
c` ∆]
cRatio
100
10.4462
2.8146
[7.9791;12.9133]
0.6311
-400
10.8351
1.4763
[9.5411;12.1292]
0.2422
2.6060
1600
10.9510
0.3810
[10.6170;11.2850]
0.1263
1.9170
6400
10.9872
0.3702
[10.6627;11.3116]
0.0901
1.4013
25600
11.0263
0.1024
[10.9365;11.1161]
0.0510
1.7675
102400
11.0552
0.0558
[11.0063;11.1041]
0.0221
2.3107
409600
11.0643
0.0069
[11.0582;11.0704]
0.0130
1.6984
Table 4.2: Approximation of call option price using Method (3.4) as m increase
Number of
simulation (m)
µ
ˆ
pσ
ˆ
µˆ[ˆ
µ
p´ ∆; ˆ
µ
p` ∆]
pRatio
100
15.2721
2.3378
[13.2229;17.3213]
0.5407
-400
16.0515
1.7234
[14.5409;17.5621]
0.2387
2.2650
1600
16.0590
0.4571
[15.6583;16.4597]
0.2462
0.9696
6400
15.8477
0.1352
[15.7293;15.9662]
0.0349
7.0446
25600
15.7992
0.1012
[15.7105;15.8879]
0.0136
2.5710
102400
15.7905
0.0574
[15.7401;15.8408]
0.0223
0.6083
409600
15.7962
0.0130
[15.7848;15.8076]
0.0166
1.3457
Table 4.3: Approximation of put option price using Method (3.4) as m increase We investigate the order of convergence of the Monte-Carlo simulation combined with our method using loglog and polyfit build-in functions in MATLAB. The function loglog plots a graphic of the errors and numbers of simulations with a log-arithmic scale, and the function polyfit computes the slope of the linear regression.
Let us define that if
“ k ¨ pLmqp then log “ log k ` p ¨ log m
The parameter p is now the slope of the linear regression that we compute using polyfit function, but also the order of convergence of the Monte-Carlo. In Figure 4.2, the first graphic presents the evolution of the mean absolute errors as the number of simulations increases. We clearly observe that the mean absolute errors converge to 0 which confirms the convergence of our method to the true values. The second graph show us the logarithmic scale of the errors depending on m, the number of simulations. The previous noticed irregularity is clearly visible in the graph. By linear regression with polyfit, we estimated a slope of ´0.44 for the call option errors and ´0.46for the put option. The exact order of convergence of the Monte-Carlo is ´0.5. We observe that Monte-Carlo simulation converges at this theoretical order for our method. 0 1 2 3 4 m 105 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 error error -- m call put 102 104 106 log(m) 10-2 10-1 100 log(error) log(error)--log(m) call put
4.3
Convergence of Runge–Kutta method
We want to analyze the evolution of the call and put option prices as the number of steps in our method increases. Let us recall that the SRK method used in our thesis has the property to be weak second order. We want to know if the SRK method ap-plied to the Gatheral model keeps the same order of convergence. We study the order of convergence by approximating prices of the call and put option with different step-size h. If our method is well constructed, we should observed a convergence of our estimates at a rate of convergence of weak order 2. In a method like ours, where vari-ables are square-root such as the volatility, the results can become complex numbers if the step-size is too large. This is due to the volatility that become negative for large step-size, then the square-root of a negative number is a complex number. We restric-ted the number of steps to be N ě 25, to avoid the volatility to become negative and obtain real values approximation for our results. We compute approximation of the option price with a constant number of batches, L “ 5, the number of simulations m “ 2 ¨ 105, and the number of steps in the interval from 25 to 95 steps.
The results are presented in Table 4.4 and 4.5 in the same way as for the number of simulations. The option price in both call and put converges to the approximated true values c and p as the number of step increases, so our method converges as desired. The standard deviation remains approximately constant or slightly increases as h de-creases. The mean absolute error decreases significantly as h decreases, respectively to an order of 10´4 and 10´3 for the call and put option prices. Some irregularities
are also observed in both call and put option’s mean absolute error, for 35 steps in the call and 45, 65 steps for the put. However, the mean absolute error is mainly de-creasing as the number of steps increases, which is in agreement with the property of convergence of the method that we use.
Number of
steps (N)
µ
ˆ
cσ
ˆ
µˆ[ˆ
µ
c´ ∆; ˆ
µ
c` ∆
]
c25
11.1688
0.0177
[11.1532;11.1843]
0.0915
35
11.1161
0.0409
[11.0802;11.1520]
0.0388
45
11.1231
0.0464
[11.0824;11.1638]
0.0458
55
11.1033
0.0474
[11.0617;11.1448]
0.0260
65
11.0876
0.0228
[11.0676;11.1076]
0.0103
75
11.0724
0.0745
[11.0071;11.1376]
0.0049
85
11.0758
0.0610
11.0223;11.1292]
0.0015
95
11.0778
0.0329
[11.0490;11.1066]
5.0337e-04
Table 4.4: Approximation of the call option price using the Runge–Kutta method on the Gatheral model
Number of
steps (N)
µp
ˆ
σ
ˆ
µˆ[ˆ
µp
´ ∆; ˆ
µp
` ∆]
p25
15.9204
0.0381
[15.8870;15.9537]
0.1076
35
15.8978
0.0310
[15.8707;15.9250]
0.0850
45
15.8219
0.0429
[15.7843;15.8595]
0.0091
55
15.8410
0.0465
[15.8002;15.8818]
0.0282
65
15.8276
0.0230
[15.8074;15.8478]
0.0148
75
15.8326
0.0433
[15.7946;15.8706]
0.0198
85
15.8033
0.0390
[15.7692;15.8375]
0.0095
95
15.8206
0.0287
[15.7954;15.8458]
0.0078
Table 4.5: Approximation of the put option price using the Runge–Kutta method on the Gatheral model
We use the same approach than in the case of number of simulations, to estimate the order of convergence of our method as h Ñ 0. We use the built-in functions loglogto plot our data and polyfit to approximate the slope of the curves. The first graph of figure 4.3 describes the evolution of the mean absolute error of our approximation in function of h. We observe that the mean absolute error tends to 0 as h decreases. In the second graph, we study the logarithmic evolution of the mean absolute error in function of h. We observe a convergence in both call and put option price, however we see that the call option converges slightly faster then the put option price. The function polyfit is used to estimate the slope of the curves in the second graph, which are in the same time the order of convergence. We estimate an order of convergence of 3.61 and 1.83 respectively, for the call and put option. In comparison with the theoretical order of convergence of the SRK method, our method converges faster for the call option price and slower for the put option price. The inconsistencies observed in the mean absolute error and the unexpected order of convergence of our method can be explained by a too small amount of simulations, restricted to m “ 2 ¨ 105, leading to some irregular behavior. Nevertheless, our numerical results indicate an obvious order of convergence higher than one.
0.01 0.015 0.02 0.025 0.03 0.035 0.04 h 0 0.02 0.04 0.06 0.08 0.1 0.12 error error -- h Call Put 0.015 0.02 0.025 0.03 0.035 0.04 log(h) 10-3 10-2 10-1 log(error) log(error) -- log(h) Call Put
Chapter 5
Conclusion
In this thesis, we proposed a second order SRK method applied to the Gatheral model to approximate prices of call and put options. The results in Section 4 confirm us that our method converges to the true option prices as m and N increases. We conclude that our method converges well in terms of Monte-Carlo simulation, with an order of convergence close to the theoretical one. However, our method if suitable to price op-tions in reasonable accuracy, do not converges at the same order of convergence that the SRK method proposed by Xiao Tang and Aiguo Xiao. We observed that the order of convergence is higher in terms of call option price and lower in terms of put option price. A new study based on our method could be done with optimized algorithms, running faster for large numbers of simulations and steps. This study could help us determine if better results could be obtained with larger values of parameters. The main issue of our method, is the amount of time to run an experiment. The average time for an experiment in our thesis is about three hours, for a reasonable accuracy. However, this accuracy could have been increased if more performance was avail-able.
The results are mainly consistent with the expected ones, and lead us to accept our method to price stock options. Some inconsistencies have been observed in the study of the convergence of the Runge-Kutta method, mainly created by a number of simulations too small. These inconsistencies could maybe also disappear if a new study with higher values is made.
Finally, for future research, we may consider optimizing the two algorithms used in this thesis to be less time-consuming with the same degree of accuracy. The gener-ation of random variables is the main consuming part of our algorithms. One optim-ization could be to generate the random variables outside the loops to increase the efficiency of the algorithms, however the memory bound could be reached for large numbers of simulations and time steps.
Bibliography
[1] Christian Bayer, Jim Gatheral, and Morten Karlsmark. Fast Ninomiya– Victoir calibration of the double-mean-reverting model. Quantitative Finance, 13(11):1813–1829, 2013.
[2] Fischer Black and Myron Scholes. The pricing of options and corporate liabilties. Journal of Political Economy, 81:637–654, 1973.
[3] John Charles Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, 2016.
[4] John C Cox, Jonathan E Ingersoll Jr, and Stephen A Ross. An intertemporal general equilibrium model of asset prices. Econometrica: Journal of the Econometric Society, pages 363–384, 1985.
[5] Kristian Debrabant and Andreas Rößler. Classification of stochastic Runge– Kutta methods for the weak approximation of stochastic differential equations. Mathematics and Computers in Simulation, 77(4):408–420, 2008.
[6] Kristian Debrabant and Andreas Rößler. Families of efficient second order Runge–kutta methods for the weak approximation of Itô stochastic differential equations. Applied numerical mathematics, 59(3-4):582–594, 2009.
[7] Leonhard Euler. Institutionum calculi integralis ... impensis Academiae Imperialis Scientiarum, Petropoli, 1768.
[8] Paul Glasserman. Monte Carlo methods in financial engineering, volume 53 of Ap-plications of Mathematics (New York). Springer-Verlag, New York, 2004.
[9] Steven L Heston. A closed-form solution for options with stochastic volatility with applications to bond and currency options. The review of financial studies, 6(2):327–343, 1993.
[10] Karl Heun. Neue methoden zur approximativen integration der differen-tialgleichungen einer unabhängigen veränderlichen. Z. Math. Phys, 45:23–38, 1900.
