School of Education, Culture and Communication
Division of Mathematics and Physics
BACHELOR’S DEGREE PROJECT IN MATHEMATICS
Multilevel Monte Carlo Simulation for American Option Pricing
by
Sabina Colakovic and Viktor Ågren
MAA322 — Examensarbete i matematik för kandidatexamen
DIVISION OF MATHEMATICS AND PHYSICS
MÄLARDALEN UNIVERSITYSchool of Education, Culture and Communication
Division of Mathematics and Physics
MAA322 — Bachelor’s Degree Project in Mathematics
Date of presentation:
1 June 2021
Project name:
Multilevel Monte Carlo Simulation for American Option Pricing
Authors:
Sabina Colakovic and Viktor Ågren
Version: 2nd June 2021 Supervisor: Doghonay Arjmand Reviewer: Milica Rancic Examiner: Daniel Andrén Comprising: 15 ECTS credits
Abstract
In this thesis, we center our research around the analytical approximation of American put options with the Multilevel Monte Carlo simulation approach. The focus lies on reducing the computational complexity of estimating an expected value arising from a stochastic differential equation. Numerical results showcase that the simulations are consistent with the theoretical order of convergence of Monte Carlo simulations. The approximations are accurate and considerately more computationally efficient than the standard Monte Carlo simulation method.
Keywords: Multilevel Monte Carlo simulation, Stochastic Differential Equations, Option pricing.
Acknowledgements
We would like to take the opportunity to show our gratitude to our supervisor, Doghonay Arjmand, for all his guidance during the completion of this thesis. We would also like to thank Milica Rancic for all the comments and suggestions.
Contents
Acknowledgements 2
1 Introduction 7
2 Preliminaries 10
2.1 Time-stepping Methods for Ordinary Differential Equations . . . 10
2.1.1 Forward Euler Method . . . 10
2.1.2 Backward Euler Method . . . 11
2.1.3 Runge-Kutta Method . . . 11
2.2 Discretization Methods for Stochastic Differential Equations . . . 12
2.2.1 The Euler Scheme . . . 12
2.2.2 Euler–Maruyama Scheme . . . 14
2.3 Generating Sample Paths . . . 14
2.3.1 Standard Brownian Motion . . . 15
2.3.2 Geometric Brownian Motion . . . 17
2.3.3 Random Walk Construction . . . 18
2.4 The ITÔ Lemma . . . 19
3 Standard Monte Carlo Method 20 3.1 Error Estimate for the Monte Carlo Method . . . 22
3.2 Variance Reduction . . . 24
4 Multilevel Monte Carlo for pricing American options 27 4.1 Black-Scholes Model . . . 27
4.2 Pricing American Options . . . 29
4.3 Multilevel Monte Carlo . . . 31
4.3.1 Complexity Theorem . . . 33
5 Numerical Results 35 6 Conclusion 40 A Algorithm 44 A.1 Basic Monte Carlo for Calculating Simulated Prices . . . 44
A.3 Multilevel Monte Carlo Simulation . . . 46
B Bachelor Degree Objectives 48
B.1 Objective 1: Knowledge and understanding in the major field of the research . 48 B.2 Objective 2: Ability to search, collect and evaluate information . . . 48 B.3 Objective 3: Identify, formulate and solve problems . . . 49 B.4 Objective 4: Ability to present orally and writing problems and solutions in
dialogue to different groups . . . 49 B.5 Objective 5: Ability in the major field to make judgements in respect to
List of Figures
2.1 Simulated standard Brownian motion paths . . . 18 4.1 Hold and exercise region for a sample path which hits the optimal exercise
boundary . . . 30 5.1 Evolution of (a) the absolute error as 𝑀 increases, and (b) computational time
as estimated error converges . . . 37 5.2 Evolution of (a) the absolute error as 𝑀 increases, and (b) the computational
time as estimated error decreases compared with estimated price . . . 38 5.3 Payoff evolution for every initial time point 𝑡 . . . 38 5.4 Payoff evolution for changes in initial stock price 𝑠 . . . 39
List of Tables
5.1 Approximation of put option price as 𝑀 increases and compassion of compu-tational time. . . 36 5.2 Approximation of put option price as 𝑁 increases and compassion of
Chapter 1
Introduction
The pricing of American options is still an active field of research within the quantitative finance world. Dissimilar to European options, there is no analytical solution for the value of the option price [5]. American options can be exercised at any time up until the expiration date. The most important evaluation of an American option is to estimate the optimal exercising time correctly to get the highest possible payoff [19]. Hertz [20] introduced the first Monte Carlo simulation in finance in his article, where he advocated for the need of a new concept as the conventional methods were not satisfactory to give accurate estimates. This statement is a well-discussed topic regarding the future of simulation in pricing options. Tilley [31] preformed the first Monte Carlo simulation on American options and in his article he attempted to weaken the arguments against the use of simulations of pricing American options. His work included an eight-step description of the algorithm and it includes bundling parameters alongside with a backwards approach of processing each time step. Carriere [16] attempted to improve the algorithm by the theory of optimal stopping with a sequential regression algorithm. His result shows that there is no difference in the formula used for American and European options. The well documented algorithm used to approximate the price for American put options is showcased in the paper.
Barraquand and Martineau [6] focused on the numerical valuation in higher dimensions. This is something previous algorithms could not calculate because of the exponentially in-creasing memory demand. The approach is successful by combining the Monte Carlo method simulation with state aggregation. The aggregation works by introducing partitioning the state space into cells. Appropriate partings results in the approximation that would be close to optimal exercising. Stratifying these states makes it possible to determine the strategies that result in the optimal payoff. Raymar and Zwecher [26] continued to work on the model calcu-lated by Barraquand and Martineau [6]. The improvements were done with a maximum price constraint and partitioning based on the second statistic. Their analysis of the applications showed that the two-factor approach is closer to the benchmark values except for cases where the option is too far out of the money. The use of the second statistic proved to be potent at making approximations for early exercising.
Broadie and Glasserman [9] found an issue when working with a single estimator. In order to work around the problem they created a model that calculates two estimates, where one estimate is biased high and the other is biased low. Subsequently, they should converge
between the estimates resulting in a better true price approximation. They found that the relative errors were quite small without a large effort on computation. Broadie, Glasserman and Jain [10] followed up on their previous work [9] with advancements to the algorithm. They complemented it with a confidence interval bounded by the estimates to reduce the variance combined with the use of European options to reduce the branching process, in conjunction with other strategies to get faster convergence rates and shorter simulation times. Carr [15] developed an algorithm to calculate the American option value and the boundaries based on randomization. This included a three-step approach for calculating expected value for a fixed randomly generated parameter. Their results proved to be computationally fast and the option value was approached for a close approximation of the true value. Longstaff and Schwartz [25] in 2000 had a similar approach to [16] and [32] on valuing American options. The unique aspect of their work is the strategy of only calculating the paths where the options are in-the-money. As a results, their LSM algorithm produces an efficient computational process and a true price approximation.
The most naive way to compute the optimal American option value is to exploit Monte Carlo methods. Caflisch and Chaudharry [13] reviewed the basic properties of American options and the difficulties of applying Monte Carlo valuation to American options. This included Branching processes to obtain the upper and lower bounds on the American option price, a Martingale optimization and the Least Squares Monte Carlo (LSM), that provided a direct method for pricing American option. Quasi-random sequences was then used to improve the performance of the LSM method. Their results proved the difficulties involved in applying Monte Carlo evaluation to American options, as well as valuable results in improving Monte Carlo simulations. Birge [7] presented some empirical evidence that Quasi-random streams can produce more accurate results with fewer iterations, compared to the standard Monte Carlo. Quasi-Monte Carlo also offers less fluctuation than standard Monte Carlo. Adding to that, it also makes on-line error approximations possible. Birge’s results showed that Quasi-Monte Carlo has an advantage in a simple model where analytical results are available. The advantage may be even greater in more complex models, such as Multilevel Monte Carlo. As Monte Carlo comes together with high computational burden, there have been other attempts in decreasing the computational cost by exploiting more interesting approaches based on Quasi-Monte Carlo methods. Longstaff and Schwartz [25] suggested that their method might be improved by the use of quasi-random points. However, Caflisch and Chaudharry [13] considered that there are two difficulties with this extension of their method. With this extension of the method, there are two potential difficulties. The issue is high dimensional and the prices along the various paths in LSM method, are correlated. Both of which can be troublesome for Quasi-Monte Carlo. Multilevel Monte Carlo, deals with a more advanced variant of Monte Carlo method. Giles [18] developed the Multilevel Monte Carlo for simulations of European option valuation. The results have shown that a Multilevel approach, using a geometric sequence of time steps, can reduce Monte Carlo’s path simulations order of complexity, thus reducing the computational burden.
Options are a form of financial contract that grants the investor the right to purchase or sell the underlying assets on or before a certain future date, known as maturity, at a mutually negotiated exercise price, called the strike price. There are several kinds of options. The most common ones are European-, Asian- and American options. In this paper, we will focus on
American options, and the key feature of the American options is that they are continuous time instruments because of their availability to be exercised at any time. Early exercising is the decision to exercise the option prior to the expiration date. Thanks to its exercising time flexibility, the most common options that are traded on exchanges are American.
In the world of quantitative finance, the pricing of American options is still an active field of research. There are no general solutions for the valuation of American options, unlike European options. One often utilizes Monte Carlo simulation in option pricing. However, Multilevel Monte Carlo, as proposed by Giles in [18], was used as an alternative method to reduce complexity and computational cost.
Assume that the stock value 𝑆(𝑡) evolves by the following stochastic differential equation, called the Itˆo geometric Brownian motion [13],
𝑑𝑆= 𝜇𝑆𝑑𝑡 + 𝜎𝑆𝑑𝑊 ,
where 𝜇 and 𝜎 are the constant average growth rate and volatility, respectively, and 𝑊 = 𝑊 (𝑡) is a standard Brownian motion. The aim of the project is to implement the Multilevel Monte Carlo method for valuing American options. The research objective is to approximate,
𝑓𝐴(𝑡, 𝑠) ≡ max
𝑡≤𝜏≤𝑇
𝐸[𝑒−𝑟 (𝜏−𝑡)max(𝐾 − 𝑆(𝜏), 0)|𝑆(𝑡) = 𝑠],
by implementing the Multilevel Monte Carlo method, where 𝐾 is the strike price, 𝑆(𝑡) is the spot price of the underlying asset, 𝜏 is the class of admissible stopping times with values in [𝑡, 𝑇 ], 𝑟 is the risk-free interest rate, and 𝑠 is the initial stock price. In the present thesis, we aim to extrapolate the ideas presented in [18] for American option valuation.
This report will start with an introduction to the theoretical preliminaries, stochastic differ-ential equations, as well as examples of numerical methods for solving those. Followed by an explanation of the standard Monte Carlo method with error estimates and variance reduction techniques, and the Multilevel Monte Carlo method to price American options. We present and analyze the numerical results in the following chapter, where we investigate the option prices depending on the values of different variables, as well as discuss the convergence rates and complexity for the different methods. Finally, we conclude and review the results, as well as suggest additional methods and further research to improve the algorithm.
Chapter 2
Preliminaries
2.1
Time-stepping Methods for Ordinary Differential
Equa-tions
The most commonly utilized time-stepping methods are the Forward Euler- and Backward
Euler methods. Forward Euler has the advantage of including an explicit update equation, making it simpler to incorporate in practice. Backward Euler, on the other hand, necessitates the solution of an implicit equation, so it is more expensive, but it has better stability properties in general.
There are two standard forms for expressing the solution of initial value problems in systems of ordinary differential equations, the first one being
𝑦0(𝑥) = 𝑓 (𝑥, 𝑦(𝑥)), 𝑦(𝑥0) = 𝑦0, (2.1) where 𝑦 is the solution assumed to be a differentiable function on the interval [𝑥0,𝑥¯] to a finite dimensional Euclidean space R𝑁
.
2.1.1
Forward Euler Method
We begin our discussion of numerical methods for initial value ordinary differential equations (ODEs) with a review of the fundamental concepts. The discretization Forward Euler [2] is used to demonstrate these principles. The problem to be solved in the general form is,
y0= f(𝑡, y), 0 ≤ 𝑡 ≤ 𝑏, (2.2) where the initial y(0) = c is given. To approximate the equation (2.2), the interval of integration is first discretized by a mesh
0 = 𝑡0 < 𝑡1 <· · · < 𝑡𝑁−1 < 𝑡𝑁 = 𝑏,
and let ℎ𝑛= 𝑡𝑛− 𝑡𝑛−1 be the 𝑛th step size. The approximations can be constructed as
y0(= c), y1, . . . ,y
with y𝑛being an intended approximation of y(𝑡𝑛). Since y0is given and proceed to integrate the
Ordinary Differential Equation (ODE) in distinct steps. Each 𝑛 ∈ [1, 𝑁] is an approximation
y𝑛−1at 𝑡𝑛−1, and the variable of interest is y𝑛at 𝑡𝑛.
When constructing a discretization method, consider the Taylor’s expansion
y(𝑡𝑛) = y(𝑡𝑛−1) + ℎ𝑛y 0(𝑡 𝑛−1) + 1 2ℎ 2 𝑛y 00(𝑡 𝑛−1) + . . . . (2.3)
Using the big O notation, equation (2.3) can also be written as
y(𝑡𝑛) = y(𝑡𝑛−1) + ℎ𝑛y 0
(𝑡𝑛−1) + 𝑂 (ℎ 2 𝑛),
where the rightmost term of the Taylor’s expansion in equation (2.3), the forward Euler method can be derived by replacing y0by f, producing the scheme, [2]
y𝑛 = y𝑛−1+ ℎ𝑛f 𝑡𝑛−1,y𝑛−1
.
2.1.2
Backward Euler Method
The backward Euler method [2] is derived from equation (2.2), just like the forward Euler method. The difference between Forward- and Backward Euler method is that using the Backward Euler method, it centers at 𝑡𝑛, rather than at 𝑡𝑛−1, which produces
y𝑛= y𝑛−1+ ℎ𝑛f(𝑡𝑛,y𝑛). (2.4)
In terms of geometry, the backward Euler method uses the tangent at the future point (𝑡𝑛,y𝑛)
rather than the tangent at (𝑡𝑛−1,y𝑛−1) as in the forward Euler method, which improves stability.
The forward Euler method is explicit, while the backward Euler method is implicit. This implies that the unknown vector y𝑛 at each step appears on both sides of the equation (2.4).
Hence, at each step, a nonlinear system of algebraic equation has to be solved. One advantage of the backward Euler method is the method’s stability. When obtaining equation (2.4) from equation (2.2), it produces
𝑦𝑛= 𝑦𝑛−1+ ℎ𝑐𝑦𝑛, i.e.,
𝑦𝑛 = (1 − ℎ𝑐)−1𝑦𝑛−1,
where 𝑐 is a constant. The backward Euler method requires fewer steps than the forward Euler method to solve a problem. However, each backward Euler step can be more costly in terms of computation time. There are several applications where the implicit method’s total computational cost is much lower than the explicit Euler method.
2.1.3
Runge-Kutta Method
The Runge-Kutta method, in numerical applications, is a widely used method for its effective use of solving problems that include differential equations. The second-order Runge-Kutta methods employs two function evaluations and yields precision proportional to the order of
ℎ2. Runge-Kutta methods are “one step” in the sense that the outcome at the end of one step is only functionally dependent on the result at the end of the previous step. If 𝑦𝑛 denotes a
computed approximation to 𝑦(𝑥𝑛), then 𝑦𝑛can be solved by
𝑦𝑛 = 𝑦𝑛−1+ ℎ
𝑠
Õ
𝑖=1
𝑏𝑖𝐹𝑖,
where 𝐹1, 𝐹2, . . . , 𝐹𝑠 are computed derivatives from approximations 𝑌1, 𝑌2, . . . , 𝑌𝑠to the solu-tion at 𝑥𝑛−1+ ℎ𝑐1, 𝑥𝑛−1+ ℎ𝑐2, . . . , 𝑥𝑛−1+ ℎ𝑐𝑠. This leads to 𝐹𝑖 = 𝑓 (𝑥𝑛−1+ ℎ𝑐𝑖, 𝑌𝑖), 𝑖 = 1, 2, . . . , 𝑠
for the differential equation system in (2.1), 𝐹𝑖 = 𝑓 (𝑌𝑖), 𝑖 = 1, 2, . . . , 𝑠. The values of the
approximations 𝑌𝑖, 𝑖 =1, 2, . . . , 𝑠 can be found from
𝑦𝑖 = 𝑦𝑛−1+ ℎ 𝑠
Õ
𝑗=1
𝑎𝑖 𝑗𝐹𝑗 𝑖 =1, 2, . . . , 𝑠.
It turns out that the 𝑐 vector components are related to the 𝐴 matrix elements by
𝑐𝑖=
𝑠
Õ
𝑗=1
𝑎𝑖 𝑗, 𝑖 =1, 2, . . . , 𝑠.
The number of stages in a method of this type is the number of 𝑌 vectors needed to compute the solution, and it is a measure of the particular method complexity [12].
Higher order methods results in larger stability region. This leads to that the higher order Runge-Kutta methods have larger stability than for example the Forward Euler method. A result in substantially larger stability leads to that larger step sizes can be used. There are also methods based on Taylor’s decomposition on several points, which results in higher order methods for ordinary differential equations. For further research on higher-order methods, see [3, 4, 1].
2.2
Discretization Methods for Stochastic Differential
Equa-tions
The most simplistic method for approximating the simulation of stochastic differential equations is the Euler scheme. One can begin by discussing the properties of the Euler scheme, then undertake an expansion to refine the Euler scheme, and finally present criteria for comparing methods. An extension of the Euler method for ordinary differential equation to stochastic differential equation is the Euler-Maruyama method. The Euler–Maruyama [19], is a method to approximate a numerical solution of a stochastic differential equation.
2.2.1
The Euler Scheme
The processes 𝑋 is satisfying the stochastic differential equation
Consider that 𝑋 take on values in R𝑑
and 𝑊 is an 𝑚-dimensional standard Brownian motion, were 𝑎 takes values in R𝑑
and 𝑏 takes values in R𝑑×𝑚.
Denote ˆ𝑋 as a time-discretized approximation of 𝑋, then the Euler approximation on a time grid 0 = 𝑡0 < 𝑡1 < · · · < 𝑡𝑚 is defined by ˆ𝑋(0) = 𝑋 (0), for 𝑖 = 0, . . . , 𝑚 − 1, and the time-discretized approximation
ˆ
𝑋(𝑡𝑖+1) = ˆ𝑋(𝑡𝑖) + 𝑎 𝑋ˆ(𝑡𝑖) [𝑡𝑖+1− 𝑡𝑖] + 𝑏 𝑋ˆ(𝑡𝑖)
√
𝑡𝑖+1− 𝑡𝑖𝑍𝑖+1,
with 𝑍1, 𝑍2, . . . m-dimensional independent standard normal random vectors, and a fixed spacing ℎ implies that 𝑡𝑖 = 𝑖ℎ.This provides that the largest of the increments 𝑡𝑖+1− 𝑡𝑖decreases
to zero. For any fixed time-step ℎ > 0, ˆ𝑋(𝑖ℎ) may be written as ˆ𝑋(𝑖) and the corresponding Euler scheme as ˆ 𝑋(𝑖 + 1) = ˆ𝑋(𝑖) + 𝑎 ( ˆ𝑋(𝑖)) ℎ + 𝑏 ( ˆ𝑋(𝑖)) √ ℎ 𝑍𝑖+1, [19].
Theorem 1 (Theorem 3.1 in [14]). Let ¯𝑋 and ¯¯𝑋 be forward Euler approximations of the
stochastic process 𝑋 : [0, 𝑇 ] × Ω → R, satisfying the stochastic differential equation
𝑑 𝑋(𝑡) = 𝑎 (𝑡, 𝑋 (𝑡)) 𝑑𝑡 + 𝑏 (𝑡, 𝑋 (𝑡)) 𝑑𝑊 (𝑡), 0 ≤ 𝑡 < 𝑇 ,
with time steps
{¯𝑡𝑛} ¯ 𝑁 𝑛=0, ¯𝑡0 =0, ¯𝑡𝑁¯ = 𝑇 , {¯¯𝑡𝑚} ¯¯ 𝑁 𝑚=0, ¯¯𝑡0=0, ¯¯𝑡𝑁¯¯ = 𝑇 , respectively, and Δ𝑡𝑚 𝑎𝑥 =max max 0≤𝑛≤ ¯𝑁−1 ¯𝑡𝑛+1− ¯𝑡𝑛, max 0≤𝑚≤ ¯¯𝑁−1 ¯¯𝑡𝑚+1− ¯¯𝑡𝑚 .
Suppose that there exists a positive constant 𝐶 such that the initial data and the given functions
𝑎, 𝑏: [0, 𝑇 ] × R → R satisfy 𝐸[| ¯𝑋(0)|2+ | ¯¯𝑋(0)|2] ≤ 𝐶, 𝐸[ ¯ 𝑋(0) − ¯¯𝑋(0) 2 ] ≤ 𝐶Δ𝑡max, and |𝑎 (𝑡, 𝑥) − 𝑎 (𝑡, 𝑦) | < 𝐶 |𝑥 − 𝑦|, |𝑏 (𝑡, 𝑥) − 𝑏 (𝑡, 𝑦) | < 𝐶 |𝑥 − 𝑦|, |𝑎 (𝑡, 𝑥) − 𝑎 (𝑠, 𝑥) | + |𝑏 (𝑡, 𝑥) − 𝑏 (𝑠, 𝑥) | ≤ 𝐶 (1 + |𝑥|)p|𝑡 − 𝑠|.
Then there is a constant 𝐾 such that
max n 𝐸[ ¯𝑋2(𝑡, ·)], 𝐸 [ ¯¯𝑋2(𝑡, ·)] o ≤ 𝐾 (𝑇 + 1), 𝑡 < 𝑇 , and 𝐸[ ¯ 𝑋(𝑡, ·) − ¯¯𝑋(𝑡, ·) 2 ] ≤ 𝐾Δ𝑡max, 𝑡 < 𝑇 .
2.2.2
Euler–Maruyama Scheme
For the scalar Stochastic Differential Equation (SDE)
𝑑𝑆𝑡 = 𝑎 (𝑆𝑡, 𝑡)𝑑𝑡 + 𝑏 (𝑆𝑡)𝑑𝑊𝑡, (2.6) divide a time-discretization of the interval [0, 𝑇 ] into time-steps 0 = 𝑡0 < 𝑡1 < · · · < 𝑡𝑁 = 𝑇. The Euler-Maruyama scheme is the most intuitive way to discretize the scalar SDE in (2.6)
ˆ 𝑋(𝑡𝑖+1) = ˆ𝑋(𝑡𝑖) + 𝑎 ( ˆ𝑋(𝑡𝑖)) (𝑡𝑖+1− 𝑡𝑖) + 𝑏 ( ˆ𝑋(𝑡𝑖)) √ 𝑡𝑖+1− 𝑡𝑖 𝑍𝑖+1, ˆ 𝑋(𝑡0) = 𝑋 (0), (2.7)
where 𝑍𝑖, . . . , 𝑍𝑛 are standard normal random variables that are independent and identically
distributed. Let 𝑡𝑖 := 𝑖ℎ, ℎ > 0, and ˆ𝑋(𝑖) := ˆ𝑋(𝑖ℎ). Then,
ˆ
𝑋(𝑖 + 1) = ˆ𝑋(𝑖) + 𝑎 ( ˆ𝑋(𝑖)) ℎ + 𝑏 ( ˆ𝑋(𝑖)) √
ℎ 𝑍𝑖+1.
Definition 1. A discretization ˆ𝑋 has strong order of convergence 𝛽 > 0 if there exist positive constants 𝑐 and ℎ0such that for all h ∈ (0, ℎ0) we have
𝐸[| ˆ𝑋(𝑛ℎ) − 𝑋 (𝑇 ) |] ≤ 𝑐ℎ𝛽.
Definition 2. A function 𝑔 : R → R is called polynomially bounded if there exist positive
constants 𝑘 and 𝑞 such that for all 𝑥 ∈ R we have |𝑔(𝑥) | ≤ 𝑘 (1 + |𝑥|𝑞
).
Consider 𝐶𝑃2𝛽+2 as the set of functions whose derivatives of orders up to 2𝛽 + 2 are polynomially bounded.
Definition 3. A discretization ˆ𝑋 has weak order of convergence 𝛽 > 0 if there exist positive constant ℎ0 and for any 𝑓 ∈ 𝐶𝑃2𝛽+2 there exists a constant 𝑐 = 𝑐( 𝑓 ) > 0 such that for all
ℎ ∈ (0, ℎ0) we have
|𝐸 [ 𝑓 ( ˆ𝑋(𝑛ℎ))] − 𝐸 [ 𝑓 ( 𝑋 (𝑇 ))] | ≤ 𝑐ℎ𝛽.
The approximation of the SDE in equation (2.7) can in certain instances be used to achieved a more refined discretization, called the Milstein scheme
ˆ 𝑋(𝑖 + 1) ≈ ˆ𝑋(𝑖) + 𝑎 ( ˆ𝑋(𝑖)) ℎ + 𝑏 ( 𝑋 (𝑖)) √ ℎ 𝑍𝑖+1+ 1 2𝑏 0( ˆ 𝑋(𝑖))𝑏 ( ˆ𝑋(𝑖)) ℎ(𝑍2 𝑖+1− 1).
The Milstein scheme has a weak convergence with order of 1 and a strong convergence of order 0.5 [19].
2.3
Generating Sample Paths
The most well-known method to simulate paths for a variety of stochastic process is the
Brownian motion. We dedicate this section to methods of precisely simulating continuous time processes over a discrete set of dates. Additionally, we present the simulation of Brownian motion paths of one dimension and the Geometric Brownian motion.
2.3.1
Standard Brownian Motion
One of the most frequently asked questions is whether the random process reaches a limit when the phase size is finer and finer. We now move on to continuous random processes to see if the model makes sense to the point that the step size reaches zero when we use lattices to model asset prices. The models are based on the SDE:
𝑑𝑆= 𝑎 (𝑡, 𝑆)𝑑𝑡 + 𝑏 (𝑡, 𝑆)𝑑𝑊
for the asset price 𝑆(𝑡), where the 𝑎(𝑡, 𝑆)𝑑𝑡 term accounts for deterministic motions, and the other term 𝑏(𝑡, 𝑆)𝑑𝑊 accounts for random motions. To develop such stochastic models, the first step is to define the Brownian motion 𝑊 (𝑡). Randomly generated piece-wise process 𝑊 (𝑡) is composed of pieces with identical statistical properties, no matter how it is sub-divided. This is because stock prices appear random on time scales. The following are the characteristics of the Brownian Motion:
(i) 𝑊 (𝑡) have independent increments. For any date 𝜏 and for any Δ𝜏 > 0, the value of Δ𝑊 = 𝑊 (𝜏 + Δ𝜏) − 𝑊 (𝜏) is independent of 𝑊 (𝑡) for all 𝑡 ≤ 𝜏. This means that the increments of Brownian motion are independent of everything that happened prior to or at the current date 𝜏.
(ii) Increments Δ𝑊 ≡ 𝑊 (𝑡2) −𝑊 (𝑡1) is Gaussian random variables with mean 0 and variance
Δ𝑡 ≡ 𝑡2− 𝑡1. Consequently,
Δ𝑊 ≡ 𝑊 (𝑡2) − 𝑊 (𝑡1) =√𝑡2− 𝑡1𝜉 ,
where 𝜉 ∼ N(0, 1) is a Gaussian random variable with mean zero and variance 1. The aim is to represent Δ𝑊 with zero mean, as it follows the random part of the asset price movements. Any non-zero term would represent a deterministic piece, hence that could be put in the drift term 𝑎(𝑡, 𝑆)𝑑𝑡.
(iii) It is easily proven that for any 𝛿 > 0, 𝑊 (𝑡) is a continuous random process by,
prob{|𝑊 (𝑡 + Δ𝑡) − 𝑊 (𝑡)| > 𝛿} = prob |𝜉 | > √𝛿 Δ𝑡 −−−−→ Δ𝑡→0 0.
The above equation is the definition of a continuous stochastic process. (iv) 𝑊 (𝑡) is almost surely nowhere differentiable.
For any 𝐾 > 0, we argue that
prob 𝑊(𝑡 − Δ𝑡) − 𝑊 (𝑡) Δ𝑡 < 𝐾 =prob |𝜉 | > 𝐾 √ Δ𝑡 −−−−→ Δ𝑡→0 0.
Hence, the probability that the slope is bounded equals zero as Δ𝑡 → 0. The above property demonstrates that 𝑊 (𝑡) is almost certainly not differentiable anywhere.
(v) The continuity and non-differentiability followed directly from the scaling of Δ𝑊. Since Δ𝑊 =pΔ𝑡𝜉, where 𝜉 is N(0, 1), and can be written as Δ𝑊 ∼ 𝑂(
√
Δ𝑡),or more concisely 𝑑𝑊 ∼ 𝑂 (
√ Δ𝑡).
To show that Δ𝑊 is Gaussian with variance Δ𝑡, which follows directly from our desire to have 𝑊(𝑡) to be sub-dividable into finer and finer intervals, each with identical properties. Consider
Δ𝑊 = 𝑊 (𝑡1) − 𝑊 (𝑡0) ≡ 𝑛−1 Õ 𝑘=1 [𝑊 (𝜏𝑘+1) − 𝑊 (𝜏𝑘)], 𝜏𝑘 = 𝑡0+ 𝑘 𝑛 (𝑡1− 𝑡0),
where each 𝛿𝑊𝑘 ≡ 𝑊 (𝜏𝑘+1) − 𝑊 (𝜏𝑘) are independent random variables with identical
distri-butions. Since the variables are independent, the variances are
Var[Δ𝑊] =
𝑛−1
Õ
𝑘=0
Var[𝑊 (𝜏𝑘+1) − 𝑊 (𝜏𝑘)] = 𝑛 · Var[𝑊 (𝑡1) − 𝑊 (𝑡0)].
Let 𝑣 (𝑦) = Var[𝑊 (𝑡 + 𝑦) − 𝑊 (𝑡)]. Thus,
𝑣(𝑡1− 𝑡0) = 𝑛 · 𝑣
𝑡1− 𝑡0 𝑛
for any 𝑡1, 𝑡0and 𝑛. Therefore, Δ𝑊 is the sum of 𝑛 independent, identically distributed variables with mean 0 and variance (𝑡1− 𝑡0)/𝑛. As 𝑛 approaches ∞, the central limit theorem guarantees
that 𝛿𝑊 is Gaussian with mean zero and variance 𝑡1− 𝑡0[19, 28].
Definition 4. (The Wiener process). [Definition 2.8 in [14]] The one-dimensional Wiener
process 𝑊 : [0, ∞) ×Ω → R, also known as the Brownian motion, has the following properties:
1. with probability 1, the mapping 𝑡 ↦→ 𝑊 (𝑡) is continuous and 𝑊 (0) = 0; 2. if 0 = 𝑡0 < 𝑡1 < · · · < 𝑡𝑁 = 𝑇, then the increments
𝑊(𝑡𝑁) − 𝑊 (𝑡𝑁−1), . . . , 𝑊 (𝑡1) − 𝑊 (𝑡0)
are independent; and
3. for all 𝑡 > 𝑠 the increment 𝑊 (𝑡) −𝑊 (𝑠) is normally distributed, with 𝐸 [𝑊 (𝑡) −𝑊 (𝑠)] = 0 and 𝐸 [(𝑊 (𝑡) − 𝑊 (𝑠))2] = Var[𝑊 (𝑡) − 𝑊 (𝑠)] = 𝑡 − 𝑠.
To generate Brownian motion, determine 𝑊 at finitely many time steps {𝑡𝑛 : 𝑛 = 0, . . . , 𝑁 } of
the form 𝑡0 < 𝑡1 < · · · < 𝑡𝑁 = 𝑇. Definition 4 shows how to generate 𝑊 (𝑡𝑛) by a sum of inde-pendent normally distributed random variables for computational method to generate independ-ent normal distributed random variables. Denote these incremindepend-ents as Δ𝑊𝑛 = 𝑊 (𝑡𝑛+1) − 𝑊 (𝑡𝑛).
By Definition 4, the Brownian motion 𝑊 (𝑡) is itself a normally distributed random variable for fixed time 𝑡. To generate 𝑊 for all 𝑡 ∈ R is computationally infeasible, since it requires infinite computational work [14].
2.3.2
Geometric Brownian Motion
A geometric Brownian motion is an exponentiation of the Brownian motion. This means that a stochastic process 𝑆(𝑡) is a geometric Brownian motion if the term log 𝑆(𝑡) is a Brownian motion with an initial value log 𝑆(0). Since the exponential function only takes positive values, geometric Brownian motion is always positive, compared to ordinary Brownian motion that can take on non-positive values. The percentage changes for geometric Brownian motion is
𝑆(𝑡2) − 𝑆(𝑡1) 𝑆(𝑡1) , 𝑆(𝑡3) − 𝑆(𝑡2) 𝑆(𝑡2) , . . . , 𝑆(𝑡𝑛) − 𝑆(𝑡𝑛−1) 𝑆(𝑡𝑛−1) ,
where the changes are independent for any time-step 𝑡1 < 𝑡2 < · · · < 𝑡𝑛, rather than the absolute changes 𝑆(𝑡𝑖+1) − 𝑆(𝑡𝑖). These characteristics explain why geometric Brownian motion, rather
than ordinary Brownian motion, is used to model asset prices. Suppose 𝑊 is a standard Brownian motion and 𝑋 satisfies
𝑑 𝑋(𝑡) = 𝜇𝑑𝑡 + 𝜎𝑑𝑊 (𝑡),
so that 𝑋 ∼ BM(𝜇, 𝜎2). Set 𝑆(𝑡) = 𝑆(0) exp( 𝑋 (𝑡)) ≡ 𝑓 ( 𝑋 (𝑡)), with an application of the Itô formula it shows that
𝑑𝑆(𝑡) = 𝑓0( 𝑋 (𝑡))𝑑𝑋 (𝑡) +1 2𝜎 2 𝑓00( 𝑋 (𝑡))𝑑𝑡 = 𝑆(0) exp( 𝑋 (𝑡)) [𝜇𝑑𝑡 + 𝜎𝑑𝑊 (𝑡)] + 12𝜎2𝑆(0) exp +( 𝑋 (𝑡))𝑑𝑡 = 𝑆(𝑡) ( 𝜇 +1 2𝜎 2)𝑑𝑡 + 𝑆(𝑡)𝜎𝑑𝑊 (𝑡). (2.8)
A stochastic differential equation of the form 𝑑𝑆(𝑡)
𝑆(𝑡)
= 𝜇𝑑𝑡 + 𝜎𝑑𝑊 (𝑡), (2.9)
is often used to describe a geometric Brownian motion. In equation (2.9), 𝑆(𝑡) has the drift 𝜇 𝑆(𝑡), which implies the following equation that can be verified through Itô formula:
𝑑log 𝑆(𝑡) = (𝜇 − 1
2𝜎
2)𝑑𝑡 + 𝜎𝑑𝑊 (𝑡). (2.10)
To show that 𝑆 is a process of the type in equation (2.9), the notation 𝑆 ∼ GBM(𝜇, 𝜎2) will be used. If 𝑆 has the initial value 𝑆(0), and from equation (2.10) one can see that 𝑆 ∼ GBM(𝜇, 𝜎2) then 𝑆(𝑡) = 𝑆(0) exp 𝜇− 1 2𝜎 2 𝑡+ 𝜎𝑊 (𝑡) . If 𝑢 < 𝑡 then 𝑆(𝑡) = 𝑆(𝑢) exp 𝜇− 1 2𝜎 2 (𝑡 − 𝑢) + 𝜎(𝑊 (𝑡) − 𝑊 (𝑢)) ,
since the increments of 𝑊 are independent and normally distributed, a simple recursive method for simulating 𝑆 values at 0 = 𝑡0 < 𝑡1 < · · · < 𝑡𝑛:
𝑆(𝑡𝑖+1) = 𝑆(𝑡𝑖) exp 𝜇− 1 2𝜎 2 (𝑡 𝑖+1− 𝑡𝑖) + 𝜎 √ 𝑡𝑖+1− 𝑡𝑖𝑍𝑖+1 ,
2.3.3
Random Walk Construction
When simulating the Brownian motion (see Fig. 2.1), one primarily concentrates on simulating values (𝑊 (𝑡1), . . . , 𝑊 (𝑡𝑛)) at a fixed set of points 0 < 𝑡1 < · · · < 𝑡𝑛. Simulating the 𝑊 (𝑡𝑖)
is straightforward since Brownian motion has independent normally distributed increments. Let 𝑍1, . . . , 𝑍𝑛 be independent standard normal random variable, for a Brownian motion with 𝑡0 =0 and 𝑊 (0) = 0. As a result, the following values can be generated
𝑊(𝑡𝑖+1) = 𝑊 (𝑡𝑖) + √
𝑡𝑖+1− 𝑡𝑖 𝑍𝑖+1, 𝑖 =0, . . . , 𝑛 − 1. When 𝑋 ∼ 𝐵𝑀 (𝜇, 𝜎2) with given 𝑋 (0) and constant 𝜇 and 𝜎, set
𝑋(𝑡𝑖+1) = 𝑋 (𝑡𝑖) + 𝜇(𝑡𝑖+1− 𝑡𝑖) + 𝜎 √
𝑡𝑖+1− 𝑡𝑖 𝑍𝑖+1, 𝑖 =0, . . . , 𝑛 − 1. Let the coefficients be time dependent. Hence, the recursion becomes
𝑋(𝑡𝑖+1) = 𝑋 (𝑡𝑖) + ∫ 𝑡𝑖+1 𝑡𝑖 𝜇(𝑠)𝑑𝑠 + s ∫ 𝑡𝑖+1 𝑡𝑖 𝜎2(𝑢)𝑑𝑢 𝑍𝑖+1, 𝑖=0, . . . , 𝑛 − 1. (2.11)
Replacing equation (2.11) with the Euler approximation one gets
𝑋(𝑡𝑖+1) = 𝑋 (𝑡𝑖) + 𝜇(𝑡𝑖) (𝑡𝑖+1− 𝑡𝑖) + 𝜎 (𝑡𝑖)
√
𝑡𝑖+1− 𝑡𝑖 𝑍𝑖+1, 𝑖=0, . . . , 𝑛 − 1.
Since the increments 𝑡1, . . . , 𝑡𝑛mean and variance are no longer exact it will have a discretization error even at the increments [19].
2.4
The ITÔ Lemma
The Itô Lemma is one of the most important formulas in financial analysis and it explains how to distinguish stochastic process functions. Begin with a Taylor expansion to the lowest orders for a function of two variables 𝐹 (𝑡, 𝑋) to understand the Ito formula in its most simple form
𝑑𝐹 = 𝜕 𝐹 𝜕 𝑡 𝑑 𝑡+ 𝜕 𝐹 𝜕 𝑋 𝑑 𝑋 + 1 2 𝜕2𝐹 𝜕 𝑡2 (𝑑𝑡)21 2 𝜕2𝐹 𝜕 𝑋2 (𝑑𝑋)2+ 𝜕 2𝐹 𝜕 𝑡 𝜕 𝑋 𝑑 𝑡 𝑑 𝑋+ . . . ,
for the stochastic process 𝑋
𝑑 𝑋 = 𝜇𝑑𝑡 + 𝜎𝑑𝑊 ,
where 𝜇 is the deterministic drift and 𝜎 the volatility, and 𝑊 is a Wiener process, also called a Brownian motion, with the property 𝑑𝑊2= 𝑑𝑡. For the lowest order it is simplified to
(𝑑𝑋)2= 𝜇2𝑑 𝑡2+ 𝜎2𝑑𝑊2+ 2𝜇𝜎𝑑𝑡𝑑𝑊 → 𝜎2𝑑 𝑡 ,
𝑑 𝑡 𝑑𝑊 and 𝑑𝑡𝑑𝑡 can be ignored in the lowest order. Hence, to the lowest order of 𝑑𝐹:
𝑑𝐹 = 𝜕 𝐹 𝜕 𝑡 + 𝜇𝜕 𝐹 𝜕 𝑋 + 1 2𝜎 2𝜕2𝐹 𝜕 𝑋2 𝑑 𝑡+ 𝜎 𝜕 𝐹 𝜕 𝑋 𝑑𝑊 , (2.12)
and, (2.12) is called the Ito’s formula [28].
Theorem 2. (Theorem 3.9 in [14]). Suppose that the assumption in Theorem(1) holds, and
that 𝑋 satisfies the stochastic differential equation
𝑑 𝑋(𝑠) = 𝑎 (𝑠, 𝑋 (𝑠)) 𝑑𝑠 + 𝑏 (𝑠, 𝑋 (𝑠)) 𝑑𝑊 (𝑠), 𝑠 >0 𝑋(0) = 𝑋0,
and let 𝑔 : (0, +∞) × R → R be a given bounded function in 𝐶2( (0, ∞) × R). Then 𝑦(𝑡) ≡
𝑔(𝑡, 𝑋 (𝑡)) satisfies the stochastic differential equation
𝑑𝑦(𝑡) = 𝜕𝑡𝑔(𝑡, 𝑋 (𝑡)) + 𝑎 (𝑡, 𝑋 (𝑡))𝜕𝑥𝑔(𝑡, 𝑋 (𝑡)) + 𝑏2(𝑡, 𝑋 (𝑡)) 2 𝜕𝑥 𝑥𝑔(𝑡, 𝑋 (𝑡)) 𝑑 𝑡 + 𝑏 (𝑡, 𝑋 (𝑡))𝜕𝑥𝑔(𝑡, 𝑋 (𝑡))𝑑𝑊 (𝑡).
Chapter 3
Standard Monte Carlo Method
Monte Carlo(MC) simulations are named after Monaco’s well-known gambling destination, and it was first developed by Stanislaw Ulam. Ulam was a mathematician who worked on the Manhattan Project in 1939. In collaboration with John von Neumann they developed the Monte Carlo simulation in 1946 [23]. A Monte Carlo simulation works by constantly repeating random samples to receive a result, where it takes the uncertain variable and assigns it a random value, and it repeats this process while assigning the variable in question to different values. When the simulation is complete, the average results provide an estimate. Consider the integral
𝛼 = ∫ 1
0
𝑓(𝑥)𝑑𝑥,
if 𝑈 is a random variable that is uniformly distributed on [0, 1]. Then independent points can be drawn of 𝑈 to evaluate the function 𝑓 at 𝑈1, 𝑈2, ..., 𝑈𝑛, and take the mean of those values to get a Monte Carlo estimate
ˆ 𝛼𝑛= 1 𝑛 𝑛 Õ 𝑖=1 𝑓(𝑈𝑖). If∫1
0 | 𝑓 (𝑥) |𝑑𝑥 < ∞, then by the strong law of large numbers,
𝑃{ lim 𝑛→∞ ˆ 𝛼𝑛 = 𝛼} =1. Assume that∫ 1 0 [ 𝑓 (𝑥)] 2𝑑𝑥 <∞, where, 𝜎2 𝑓 = ∫ 1 0 [ 𝑓 (𝑥) − 𝛼]2𝑑𝑥 ,
then by the Central Limit Theorem [14] (see Theorem 3), the distribution of the error 𝛼𝑛− 𝛼
is approximately normally distributed with mean 0 and standard deviation 𝜎𝑓/
√
𝑛, and the accuracy of this approximation is improving with an increase of 𝑛. The unknown parameter 𝜎𝑓 can be estimated using the sample standard deviation
𝑠𝑓 = v t 1 𝑛− 1 𝑛 Õ 𝑖=1 ( 𝑓 (𝑈𝑖) − ˆ𝛼𝑛)2.
In particular, from the values of the function 𝑓 (𝑈1), 𝑓 (𝑈2), . . . , 𝑓 (𝑈𝑛) an estimate is obtained
of the integral 𝛼 and a measure of the error. The standard error is proportional to 1/√𝑛, and decreasing it twice requires increase of 𝑛 by a factor of 4 [19].
Monte Carlo simulation is used in finance to value and analyse financial instruments, portfolios, investments, and the method can evaluate and price options [18]. It works by simulating the sources of uncertainty that affect their value, The price of the underlying asset is first simulated by random generation for several paths. The value of the option can be approximated by calculating the average of discounted returns over all paths. The option is priced under risk-neutral measure, since the discount rate is equal to the risk-free interest rate. The variance of the estimator should approach zero, in order to get a good estimate from this simulation [30].
In Monte Carlo path simulation, one is interested in the expected value of a quantity, and it is a functional of the solution to a stochastic differential equation. Consider the SDE
𝑑𝑆(𝑡) = 𝑎(𝑆, 𝑡)𝑑𝑡 + 𝑏(𝑆, 𝑡)𝑑𝑊 (𝑡), 0 < 𝑡 < 𝑇 , (3.1) given the initial data 𝑆0 to compute the expected value of 𝑓 (𝑆(𝑇 )), where 𝑓 (𝑆) is a scalar
function with a uniform Lipschitz bound. Hence, there exists a constant 𝑐 such that
| 𝑓 (𝑈) − 𝑓 (𝑉 ) | ≤ 𝑐 k 𝑈 − 𝑉 k, ∀𝑈, 𝑉 . (3.2) Euler discretization of this SDE with time-step ℎ is
ˆ
𝑆𝑛+1 = ˆ𝑆𝑛+ 𝑎 ( ˆ𝑆𝑛, 𝑡𝑛) ℎ + 𝑏 ( ˆ𝑆𝑛, 𝑡𝑛)Δ𝑊𝑛.
The estimate for 𝐸 [ 𝑓 (𝑆𝑇)] is the mean of the payoff values 𝑓 ( ˆ𝑆𝑇/ℎ) from 𝑁 independent path
simulations, ˆ 𝑌 = 𝑁−1 𝑁 Õ 𝑖=1 𝑓 ˆ 𝑆(𝑖) 𝑇/ℎ .
Provided that 𝑎(𝑆, 𝑡) and 𝑏(𝑆, 𝑡) satisfy certain conditions, the expected mean-square-error (MSE) in the estimate ˆ𝑌 is asymptotic to the form
𝑀 𝑆 𝐸 ≈ 𝑐1𝑁−1+ 𝑐2ℎ2, (3.3) where 𝑐1 and 𝑐2 are positive constants. In Monte Carlo sampling, the first term in (3.3)
corresponds to the variance in ˆ𝑌, and the second term introduced by the Euler discretization is the square of the 𝑂 (ℎ) bias [18].
Consider the SDE
𝑑 𝑋(𝑡) = 𝑎 (𝑡, 𝑋 (𝑡)) 𝑑𝑡 + 𝑏 (𝑡, 𝑋 (𝑡)) 𝑑𝑊 (𝑡)
on the interval 𝑡0 ≤ 𝑡 ≤ 𝑇 . To compute the value 𝐸 [𝑔( 𝑋 (𝑇 ))], Monte Carlo method can be
used based on the following approximation,
𝐸[𝑔( 𝑋 (𝑇 ))] ' 𝑁 Õ 𝑗=1 𝑔( ¯𝑋(𝑇 ; 𝜔𝑗)) 𝑁 ,
where ¯𝑋 is an approximation of 𝑋. The error in the Monte Carlo method is given by 𝐸[𝑔( 𝑋 (𝑇 ))] − 𝑁 Õ 𝑗=1 𝑔( ¯𝑋(𝑇 ; 𝜔𝑗)) 𝑁 =𝐸 [𝑔( 𝑋 (𝑇 )) − 𝑔( ¯𝑋(𝑇 ))] − 𝑁 Õ 𝑗=1 𝑔( ¯𝑋(𝑇 ; 𝜔𝑗)) − 𝐸 [𝑔( ¯𝑋(𝑇 ))] 𝑁 . (3.4)
The first term of the error representation in (3.4) is the time discretization error and the second term is the statistical error [14]. To understand the statistical error of Monte Carlo methods, The Central Limit Theorem is the fundamental result.
Theorem 3. (The Central Limit Theorem).[Theorem 5.2 in [14]] Assume 𝜉𝑛, 𝑛 = 1, 2, 3, . . . , are independent, identically distributed and 𝐸 [𝜉𝑛] = 0, 𝐸 [𝜉
2 𝑛] = 1. Then 𝑁 Õ 𝑛=1 𝜉𝑛 √ 𝑁 ⇀ 𝜈, (3.5)
where 𝜈 ∼ N(0, 1) and ⇀ denotes convergence of the distribution. The convergence of (3.5) means 𝐸 [𝑔(Í𝑁
𝑛=1𝜉𝑛/
√
𝑁)] → 𝐸 [𝑔(𝜈)] for all bounded and continuous functions 𝑔.
Theorem 4(Theorem 5.8 in [14]). Assume that 𝑎, 𝑏 and 𝑔 are differentiable to any order and
these derivatives are bounded, then there holds
𝐸[𝑔( 𝑋 (𝑇 )) − 𝑔( ¯𝑋(𝑇 ))] = 𝑂 (max Δ𝑡).
3.1
Error Estimate for the Monte Carlo Method
Let 𝑦 be the output and 𝑁 be the sample size for the Monte Carlo simulation. Let the individual values from each trail be labelled 𝑦𝑖 where 𝑖 = 1, . . . , 𝑁. Consider the error made
when estimating the expected value of 𝑦 when using a sample that has 𝑁 trials. If the probability density function of 𝑦 is 𝑓 (𝑦), then the expected value of 𝑦 is,
𝜇𝑦 = 𝐸 [𝑦] =
∫ ∞
−∞
𝑦 𝑓(𝑦)𝑑𝑦.
For 𝑁 trials, the sample mean, would be a reasonable estimator for 𝜇𝑦,
¯𝑦 = 1 𝑁 𝑁 Õ 𝑖=1 𝑦𝑖.
Monte Carlo simulation involves pseudo-random draws of the inputs. Hence, one will obtain different results each time the simulation is performed, and slightly different results for ¯𝑦. The variability in the results of ¯𝑦 depends on 𝑁, the number of trials in each Monte Carlo
simulation, also known as the moments of the error. The first question to ask is how accurate is ¯𝑦 as an estimate of 𝜇𝑦. To answer this, calculate the expectation of ¯𝑦 − 𝜇𝑦,
𝐸[ ¯𝑦 − 𝜇𝑦] = 𝐸 [ ¯𝑦] − 𝜇𝑦 = 𝐸 [1 𝑁 𝑁 Õ 𝑖=1 𝑦𝑖] − 𝜇𝑦 = 1 𝑁 𝑁 Õ 𝑖=1 𝐸[𝑦𝑖] − 𝜇𝑦.
In the Monte Carlo method, the 𝑦𝑖term occur from a random sampling of the inputs, 𝐸 [𝑦𝑖] = 𝜇𝑦.
Thus,
𝐸[ ¯𝑦 − 𝜇𝑦] = 1 𝑁
𝑁 𝜇𝑦− 𝜇𝑦 =0.
This results shows that on average, the error in using ¯𝑦 to approximate 𝜇𝑦 is zero, also known
as an unbiased estimator. To quantify the variability in ¯𝑦, use the variance of ¯𝑦 − 𝜇𝑦. Since the
term 𝜇𝑦 is a constant, the variance is zero. Therefore,
Var[ ¯𝑦 − 𝜇𝑦] = Var[ ¯𝑦] − Var[𝜇𝑦]
=Var[ ¯𝑦] =Var "Í𝑁 𝑖=1𝑦𝑖 𝑁 # = 1 𝑁2 Var " 𝑁 Õ 𝑖=1 𝑦𝑖 # .
The Monte Carlo method draws independent random samples; Therefore, the variance of the sum from samples 𝑦𝑖is the sum of their variance. Hence,
Var[ ¯𝑦 − 𝜇𝑦] = 1 𝑁2 𝑁 Õ 𝑖=1 Var[𝑦𝑖] = 1 𝑁2 𝑁 𝜎2 𝑦 = 𝜎2 𝑦 𝑁 .
Therefore, the following equations holds:
𝜇𝑦 ≡ 𝐸 [ ¯𝑦] = 𝜇𝑦, 𝜎2 ¯𝑦 ≡ 𝐸 [( ¯𝑦 − 𝜇𝑦) 2] = 𝜎 2 𝑦 𝑁 .
The quantity 𝜎¯𝑦 is known as the standard error of the estimator, and as the square root of the
sample size√𝑁increases, the standard error decreases. To reduce the uncertainty in the mean estimate by a factor of ten, the number of Monte Carlo trials must be increased by a factor of 100.
To approximate the distribution of ¯𝑦, the central limit theorem can be applied for large sample size 𝑁. The central limit theorem says for large 𝑁, the distribution of ¯𝑦 − 𝜇𝑦 will
approach a normal distribution with mean 0 and variance 𝜎
2 𝑦 𝑁, 𝑓( ¯𝑦 − 𝜇𝑦) → N (0, 𝜎¯𝑦) = N 0, 𝜎2 𝑦 𝑁 .
In a practical situation one cannot calculate the above error estimate of the confidence intervals because they depend on 𝜎𝑦, an unknown value. Therefore, an unbiased estimate of 𝜎
2 𝑦 is utilized 𝑠2 𝑦 ≡ 1 𝑁 − 1 𝑁 Õ 𝑖=1 (𝑦𝑖− ¯𝑦) 2 .
This introduces additional uncertainty in the quality of the estimate, and for small sample sizes it could be significant [14, 17].
3.2
Variance Reduction
Monte Carlo methods are simulation algorithms that are used to approximate a numerical quantity in real-world mathematical applications, often programmed in languages such as Matlab, R, and Python. Although computer speeds and computational times has increased dramatically over the years, one still needs variance reduction techniques for increasing the efficiency of Monte Carlo simulation by reducing the variance of simulation estimates. Simple variance reduction techniques are often surprisingly successful and simple to implement. It provides possibilities to what would otherwise be difficult in certain applications, such as rare event simulation and quantum chemistry. Different variance reduction techniques are used in the most advanced Monte Carlo simulations, such as monitor variates, partial integration, systematic sampling, and significance sampling. If we can simplify a version of the problem, it can sometimes be solved directly with control variates approach. This is frequently the case in simple problems, such as quantitative finance pricing problems.
The error variance of Monte Carlo integration is usually in the form 𝜎2/𝑛, where 𝑛 denotes the number of samples. By increasing the value of 𝑛, the approximation of the result is better, but the computational time increases. Instead, one will often find a way to minimize the variance 𝜎2. To accomplish this, create a new Monte Carlo problem with the same solution as our previous one but a lower variance.
The control variates method is one of the most efficient and widely used techniques for increasing the Monte Carlo simulations performance, by reducing the variance. It uses knowledge about the errors in known-quantity estimates to reduce the error in an unknown-quantity estimate. Let 𝑌1, . . . , 𝑌𝑛be the outputs from 𝑛 replications of a simulation. Suppose that 𝑌𝑖, 𝑖 = 1, . . . , 𝑛 are independent and identically distributed, where the aim is to estimate
𝐸[𝑌𝑖]. The sample mean ¯𝑌 = (𝑌1+ · · · + 𝑇𝑛)/𝑛 is the most common estimator for practical use. The estimator is unbiased and converges as 𝑛 → ∞ with a probability of 1. Suppose that another output 𝑋𝑖along with 𝑌𝑖could be calculated on each replication. Assume that the pairs
(𝑋𝑖, 𝑌𝑖), 𝑖 = 1, . . . , 𝑛 are independent and identically distributed the expected value, 𝐸 [𝑋] of
the output 𝑋𝑖are known. For any fixed 𝑏
𝑌𝑖(𝑏) = 𝑌𝑖− 𝑏 ( 𝑋𝑖− 𝐸 [𝑋]),
would be calculated using the 𝑖th replication and thereafter compute the sample mean
¯ 𝑌(𝑏) = ¯𝑌 − 𝑏 ( ¯𝑋− 𝐸 [𝑋]) = 1 𝑛 𝑛 Õ 𝑖=1 (𝑌𝑖− 𝑏 ( 𝑋𝑖− 𝐸 [𝑋])). (3.6)
The observed error ¯𝑋 − 𝐸 [𝑋] is provided as a control in estimating 𝐸 [𝑌 ], and known as the control variate estimator [24].
Lemma 1. The control variate estimator(3.6), is unbiased and consistent, [19].
Proof.
𝐸[ ¯𝑌(𝑏)] = 𝐸 [ ¯𝑌 − 𝑏 ( ¯𝑋 − 𝐸 [𝑋])] = 𝐸 [ ¯𝑌] = 𝐸 [𝑌 ] This is consistent with probability of 1 because
lim 𝑛→∞ 𝑛 Õ 𝑖=1 𝑌𝑖(𝑏) = lim 𝑛→∞ 1 𝑛 𝑛 Õ 𝑖=1 (𝑌𝑖− 𝑏 ( 𝑋𝑖− 𝐸 [𝑋])) = 𝐸 [𝑌 − 𝑏 ( 𝑋 − 𝐸 [𝑋])] = 𝐸 [𝑌 ].
The variance or each 𝑌𝑖(𝑏) is
Var[𝑌𝑖(𝑏)] = Var[𝑌𝑖− 𝑏 ( 𝑋𝑖− 𝐸 [𝑋])] = 𝜎𝑌2− 2𝑏𝜎𝑋𝜎𝑌𝜌𝑋𝑌 + 𝑏 2 𝜎2 𝑋 ≡ 𝜎 2(𝑏), where 𝜎𝑋2 =Var[𝑋], 𝜎 2
𝑌 =Var[𝑌 ] and 𝜌𝑋𝑌 is the correlation coefficient between 𝑋 and 𝑌 .
Since the control variate estimator ¯𝑌(𝑏) has a variance of 𝜎2(𝑏)/𝑛 and the sample mean ¯𝑌 has a variance of 𝜎𝑌2/𝑛. This implies that the control variate estimator has a smaller variance than the standard estimator as long as the constraint 𝑏2𝜎𝑋 < 2𝑏𝜎𝑌𝜌𝑋𝑌 holds. To minimize the variance (3.2) the optimal coefficient 𝑏∗ can be used,
𝑏∗ = 𝜎𝑌 𝜎𝑋 𝜌𝑋𝑌 = Cov[𝑋, 𝑌 ] Var[𝑋] .
Substitute the value of 𝑏∗ from (3.2) into the equation (3.6) and simplify, which produce the ratio of the variance from the optimally controlled estimator to that of the uncontrolled estimator
Var[ ¯𝑌 − 𝑏∗( ¯𝑋− 𝐸 [𝑋])] Var[ ¯𝑌]
=1 − 𝜌2𝑋𝑌.
Several potential control variates can often be combined into a single control variate, although this is not always possible. Since there are several different control variates, the issue of
multiple controls is intriguing. The previous findings can easily be generalized to several control variables. Let ¯Xdenote the vector of sample means of the controls. For a fixed 𝑏 ∈ R𝑑
the control variate estimator ¯𝑌(𝑏) is ¯
𝑌(𝑏) = ¯𝑌 − 𝑏>( ¯X − 𝐸 [𝑋]),
is an unbiased and consistent estimator of 𝐸 [𝑌 ]; see Lemma (1). The optimal parameter values in the vector b are given by
b∗ = Σ−1𝑋 Σ𝑋𝑌,
where 𝚺−1𝑋 and 𝚺𝑋𝑌 corresponds to the covariance matrix of 𝑋 and the vector of covariances
between (𝑌 , 𝑋). The estimation of b∗is given by
ˆ b∗𝑛 = S
−1 𝑋 S𝑋𝑌,
where 𝑆𝑋 is the 𝑑 × 𝑑 matrix with 𝑗 , 𝑘 entry and 𝑆𝑋𝑌 is the 𝑑-vector with 𝑗 th entry.
𝑆 𝑋( 𝑗)𝑋( 𝑘) = 1 𝑛− 1 𝑛 Õ 𝑖=1 𝑋( 𝑗 ) 𝑖 𝑋¯ (𝑘) 𝑖 − 𝑛 ¯𝑋 ( 𝑗 ) ¯ 𝑋(𝑘) ! , 𝑗 , 𝑘 =1, . . . , 𝑑, 𝑆 𝑋𝑌( 𝑗) = 1 𝑛− 1 𝑛 Õ 𝑖=1 𝑋( 𝑗 ) 𝑖 𝑌𝑖− 𝑛 ¯𝑋 ( 𝑗 )¯ 𝑌 ! , 𝑗 =1, . . . , 𝑑,
Chapter 4
Multilevel Monte Carlo for pricing
American options
As seen in previous parts of this study, standard Monte Carlo provides a reasonably straightfor-ward method for computing random variable expectations. In the context of financial analysis, this entails combining traditional Monte Carlo with a numerical method for SDEs to calculate the estimated value of the payout that corresponds to the option value. The only disadvantage is that standard Monte Carlo has a slow convergence rate, which results in the high computational cost. As Giles shows in [18] it is possible to build on this convergence rate and significantly reduce simulation costs. This method, or approach to simulation, is known as Multilevel Monte
Carlo (MLMC). In the following sections the most basic version of MLMC will be explained and using that knowledge we price the American option.
4.1
Black-Scholes Model
In the year 1973, Myron Scholes and Fisher Black published their paper The Pricing of
Options and Corporate Liabilities [8]. Robert C. Merton later expanded the mathematical understanding of the option pricing model, later named Black-Scholes option pricing model. The Black-Scholes partial differential equation is a Nobel Prize winning result in economics (1997) and one of the most famous models in finance.
Let 𝑓 (𝑡, 𝑆(𝑡)) be the price of an European put option where 𝑆(𝑡) is the price of a stock. The price of a stock satisfies the stochastic differential equation 𝑑𝑆 = 𝜇𝑆𝑑𝑡 + 𝜎𝑆𝑑𝑊, where the drift 𝜇 and the volatility 𝜎 are constants. Assume that there exist a risk-free paper, 𝐵, which follows 𝑑𝐵 = 𝑟 𝐵𝑑𝑡, where 𝑟 is the constant risk-free rent. First, consider the portfolio 𝐼 = − 𝑓 + 𝛼𝑆 + 𝛽𝐵 for 𝛼(𝑡), 𝛽(𝑡) ∈ R to find the partial differential equation of the price, and
𝑓(𝑡, 𝑆(𝑡)) of an option. Then by Itô’s formula
𝑑𝐼 = −𝑑𝑓 + 𝛼𝑑𝑆 + 𝛽𝑑𝐵 = −( 𝑓𝑡+ 𝜇𝑆 𝑓𝑆+ 1 2𝜎 2𝑆2𝑓 𝑆 𝑆)𝑑𝑡 − 𝑓𝑆 𝜎𝑆 𝑑𝑊 + 𝛼( 𝜇𝑆𝑑𝑡 + 𝜎𝑆𝑑𝑊 ) + 𝛽𝑟 𝐵𝑑𝑡 = −( 𝑓𝑡+ 𝜇𝑆 𝑓𝑆+ 1 2𝜎 2𝑆2𝑓 𝑆 𝑆) + (𝛼𝜇𝑆 + 𝛽𝑟 𝐵) 𝑑 𝑡+ (− 𝑓𝑆+ 𝛼)𝜎𝑆𝑑𝑊 .
Now choose 𝛼 such that the portfolio 𝐼 becomes risk-less. Hence, 𝛼 = 𝑓𝑆, 𝑑𝐼 = −( 𝑓𝑡+ 𝜇𝑆 𝑓𝑆+ 1 2𝜎 2 𝑆2𝑓𝑆 𝑆) + ( 𝑓𝑆𝜇 𝑆+ 𝛽𝑟 𝐵) 𝑑 𝑡 = −( 𝑓𝑡+ 1 2𝜎 2 𝑆2𝑓𝑆 𝑆) + 𝛽𝑟 𝐵 𝑑 𝑡 . (4.1)
Assume the existence of an arbitrage opportunity is disallowed. Implying that 𝑑𝐼 = 𝑟 𝐼 𝑑𝑡, where 𝑟 is the interest rate for risk-less investment,
𝑑𝐼 = 𝑟 (− 𝑓 + 𝛼𝑆 + 𝛽𝐵)𝑑𝑡 = 𝑟 (− 𝑓 + 𝑓𝑆𝑆+ 𝛽𝐵)𝑑𝑡.
(4.2)
Equations (4.1) and (4.2) show that
𝑓𝑡+ 𝑟 𝑠 𝑓𝑠+ 1
2𝜎 2𝑠2𝑓
𝑠 𝑠 = 𝑟 𝑓 , 𝑡 < 𝑇 ,
at maturity time 𝑇 the contract value is given by the definition,
𝑓(𝑇 , 𝑠) = max(𝐾 − 𝑠, 0).
Equation (4.1) is a deterministic partial differential equation and is called the Black-Scholes equation.
Theorem 5. (Theorem 4.1 in [14]). Suppose that 𝑎, 𝑏, and 𝑔 are differentiable to any order and these derivatives are bounded. Let 𝑋 be the solution of the stochastic differential equation,
𝑑 𝑋(𝑡) = 𝑎 (𝑡, 𝑋 (𝑡)) 𝑑𝑡 + 𝑏 (𝑡, 𝑋 (𝑡)) 𝑑𝑊 (𝑡),
and let 𝑢(𝑥, 𝑡) = 𝐸 [𝑔( 𝑋 (𝑇 )) | 𝑋 (𝑡) = 𝑥]. Then 𝑢 is the solution of the Komogorov backward equation 𝐿∗𝑢 ≡ 𝑢𝑡+ 𝑎𝑢𝑥+ 1 2𝑏 2 𝑢𝑥 𝑥 =0, 𝑡 < 𝑇 𝑢(𝑥, 𝑇 ) = 𝑔(𝑥).
Theorem 6. (Feynman-Ka𝑐ˇ)(Theorem 4.7 in [14]). Suppose that 𝑎, 𝑏, 𝑔, ℎ and 𝑉 are bounded
smooth functions. Let X be the solution of the stochastic differential equation 𝑑𝑋 (𝑡) =
𝑎(𝑡, 𝑋 (𝑡))𝑑𝑡 + 𝑏 (𝑡, 𝑋 (𝑡))𝑑𝑊 (𝑡) and let 𝑢(𝑥, 𝑡) =𝐸 [𝑔( 𝑋 (𝑇 )) exp( ∫ 𝑇 𝑡 𝑉(𝑠, 𝑋 (𝑠))𝑑𝑠) | 𝑋 (𝑡) = 𝑥] + 𝐸 [− ∫ 𝑇 𝑡 ℎ(𝑠, 𝑋 (𝑠))𝑒 ∫𝑠 𝑡 𝑉(𝜏,𝑋 (𝜏))𝑑𝜏 , 𝑑𝑠| 𝑋 (𝑡) = 𝑥].
Then 𝑢 is the solution of the partial differential equation
𝐿∗ 𝑉𝑢 ≡ 𝑢𝑡+ 𝑎𝑢𝑥+ 1 2𝑏 2𝑢 𝑥 𝑥+ 𝑉𝑢 = ℎ, 𝑡 < 𝑇 (4.3) 𝑢(𝑥, 𝑇 ) = 𝑔(𝑥).
Comparing equation (4.1) with equation (4.3), one can reasonably conclude that 𝑢 cor-responds to 𝑓 , 𝑋 to ˜𝑆, 𝑎(𝑡, 𝑥) = 𝑟𝑥, 𝑏(𝑡, 𝑥) = 𝜎𝑥, 𝑉 = −𝑟 and ℎ = 0. Using Feynman-Ka ˇ𝑐 (Theorem 6),
𝑓(𝑡, ˜𝑆(𝑡)) = 𝐸 [𝑒−𝑟 (𝑇 −𝑡)max(𝐾 − ˜𝑆(𝑇 ), 0)] with
𝑑 ˜𝑆 = 𝑟 ˜𝑆 𝑑 𝑡+ 𝜎 ˜𝑆 𝑑𝑊 .
This creates a crucial connection between the approximation based on the Monte Carlo method and the partial differential equations [14].
4.2
Pricing American Options
The price of American options can be obtained by maximizing over all the stopping times 𝜏∈ [0, 𝑇 ] strategies. This depends on the stock price up to the stopping time [14, 29],
𝑓𝐴(𝑡, 𝑠) ≡ max 𝑡≤𝜏≤𝑇
𝐸[𝑒−𝑟 (𝜏−𝑡)max(𝐾 − 𝑆(𝜏), 0)|𝑆(𝑡) = 𝑠]. (4.4)
For calculating the optimal selling strategy for an American option, assume that the option is only allowed to be exercised at the discrete time levels 0, Δ𝑡, 2Δ𝑡, . . . , 𝑇 . Consider the small time step (𝑇 − Δ𝑡, 𝑇 ) and by assumption the option is not exercised in that step. Therefore, the value 𝑓 (𝑡, 𝑠) holds, where 𝑓 (𝑇 , 𝑠) = max(𝐾 − 𝑠, 0) and for 𝑇 − Δ < 𝑡 < 𝑇
𝑓𝑡 + 𝑟 𝑆 𝑓𝑆+ 1 2𝜎
2
𝑆2𝑓𝑆 𝑆 = 𝑟 𝑓 . (4.5)
If a fixed stock price 𝑠 = 𝑆(𝑇 − Δ𝑡), holds the expression 𝑓 (𝑇 − Δ𝑡, 𝑠) < max(𝐾 − 𝑠, 0). When keeping the option, the option gives the expected value 𝑓 (𝑇 − Δ𝑡, 𝑠). Observe that the expected value when keeping the options is less than the value max(𝐾 − 𝑠, 0), which is obtained by selling the option at time 𝑇 − Δ𝑡. Therefore, one can conclude that it is optimal to sell if 𝑓 (𝑇 − Δ𝑡, 𝑠) < max(𝐾 − 𝑠, 0) ≡ 𝑓𝐹. By modifying the initial data at 𝑡 = 𝑇 − Δ𝑡 to
max( 𝑓 (𝑇 − Δ𝑡, 𝑠), 𝑓𝐹) and repeat the step (4.5) for (𝑇 − 2Δ𝑡, 𝑇 − Δ𝑡). The price of an American
option can be obtained as the limit of this solution as Δ𝑡 → 0 [14].
Consider an American put option, which gives the option holder the right to sell the underlying asset at the exercise price at any time before the option expires. Let 𝐾 be the exercise price, 𝑇 the maturity time, 𝑟 the risk free interest rate, and 𝜎 to be the volatility. The price 𝑆(𝑡) of the option can be obtained from the solution of the linear complementary problem
𝜆= 𝜕 𝑆(𝑡) 𝜕 𝑡 + 1 2𝜎 2𝑥2𝜕2𝑆(𝑡) 𝜕 𝑥2 + 𝑟𝑥𝜕 𝑆(𝑡) 𝜕 𝑥 − 𝑟 𝑆, 𝑥 >0 and 𝑡 ∈ [0, 𝑇 ], [𝑆(𝑡) − (𝐾 − 𝑥)]𝜆 = 0, 𝑥 >0 and 𝑡 ∈ [0, 𝑇 ], 𝑆(𝑡) − (𝐾 − 𝑥) ≥ 0, 𝜆 ≥ 0, 𝑥 >0 and 𝑡 ∈ [0, 𝑇 ], 𝑆(𝑡) = max(𝐾 − 𝑥, 0), 𝑥 >0 and 𝑡 = 𝑇 , 𝑆(𝑡) = 𝐾, 𝑥 =0 and 𝑡 ∈ [0, 𝑇 ], 𝑆(𝑡) → 0, 𝑥 → ∞ and 𝑡 ∈ [0, 𝑇 ],
where 𝑡 denote the time and 𝑥 the value of the underlying asset [22].
Consider a general class of continuous-time American option pricing by specifying a process 𝑈 (𝑡), 0 ≤ 𝑡 ≤ 𝑇 that represents the discounted payoff from exercise at time 𝑡, and a class of admissible stopping time 𝜏 with values in the interval [0,𝑇 ]. The optimal expected discounted payoff is
sup
𝜏∈T
𝐸[𝑈 (𝜏)].
The payoff to the holder of the option from exercise at time 𝑡 is ˜ℎ( 𝑋 (𝑡)), where ˜ℎ denotes a non negative payoff function. Suppose the existence of an instantaneous short rate process {𝑟 (𝑡), 0 ≤ 𝑡 ≤ 𝑇 }, the pricing problem are denoted by
sup 𝜏∈T 𝐸 h 𝑒− ∫𝜏 0 𝑟(𝑢)𝑑𝑢ℎ˜( 𝑋 (𝜏)) i .
The assumption is taken with respect to the risk-neutral measure, which is implied in the form of discounting in this expression. By making the discount factor a component of 𝑋 (𝑡) it can absorb the discounting into the function ˜ℎ. It is normal to consider the admissible stopping rules 𝜏 to be functions of the current state in this Markovian context, augmenting the state vector if necessary. Suppose the option expires at 𝑇 . The value of the option at time 0 is
sup
𝜏∈T
𝐸 𝑒−𝑟 𝜏(𝐾 − 𝑆(𝜏))+,
where 𝜏 are stopping times with respect to 𝑆, with values in the interval [0,𝑇 ]. To achieve this supremum, denote 𝜏∗ as the optimal stopping time with the form
𝜏∗ =inf{𝑡 ≥ 0 : 𝑆(𝑡) ≤ 𝑏∗(𝑡)}, with an optimal exercise boundary 𝑏∗(see Fig.4.1).
Using simulation methods for pricing American options, limit the options so that it can only be exercised at a defined collection of exercise opportunities 𝑡1 < 𝑡2 < · · · < 𝑡𝑚. The discrete-time process 𝑋0 = 𝑋 (0), 𝑋1, . . . , 𝑋𝑚 is a Markov chain on R
𝑑
. Consider American option pricing to be based on simulation of this Markov chain, and simulating 𝑋𝑖+1from 𝑋𝑖may
imply simulating 𝑋 (𝑡) values with 𝑡𝑖 < 𝑡 < 𝑡𝑖+1time steps. To apply a time discretization using
a time step smaller than the intervals 𝑡𝑖+1− 𝑡𝑖 between the exercise dates, might be necessary
[19].
4.3
Multilevel Monte Carlo
Consider Monte Carlo simulation with different time-steps ℎℓ = 𝑀 −ℓ
𝑇, with ℓ = 0, 1, . . . , 𝐿, number of simulations 𝑀, and time to maturity 𝑇 , for a given the Brownian motion 𝑊 (𝑡). Let 𝑃 denote the payoff 𝑓 (𝑆(𝑇 )), ˆ𝑆
ℓ, 𝑀ℓ denote the approximations of 𝑆(𝑇 ) and ˆ𝑃ℓ denote
the approximations of 𝑃. This is done by using a numerical discretization with time-step ℎℓ.
The idea is to write the expected payoff with a fine discretization using 2ℓ
uniform time-steps as a telescopic sum. Let ˆ𝑃ℓ be the simulated payoff with a discretization using 2ℓ uniform time-steps, 𝐸[ ˆ𝑃𝐿] = 𝐸 [ ˆ𝑃0] + 𝐿 Õ ℓ=1 𝐸[ ˆ𝑃ℓ − ˆ𝑃ℓ−1].
The Multilevel Monte Carlo method estimates each of the expectations on the right-hand side independently, in a way to minimize the computational complexity. For 𝑁0number of samples,
let ˆ𝑌0be an estimator for 𝐸 [ ˆ𝑃0], and ˆ𝑌ℓ for ℓ > 0 be an estimator for 𝐸 [ ˆ𝑃ℓ − ˆ𝑃ℓ−1] using 𝑁ℓ paths. The estimator used, is a mean of 𝑁ℓ independent samples. For ℓ > 0,
𝐸[ ˆ𝑃ℓ− ˆ𝑃ℓ−1] ≈ ˆ𝑌ℓ = 𝑁 −1 ℓ 𝑁ℓ Õ 𝑖=1 ˆ 𝑃(𝑖) ℓ − ˆ𝑃 (𝑖) ℓ−1 . (4.6) The quantity ˆ𝑃(𝑖) ℓ − ˆ 𝑃(𝑖)
ℓ−1 comes from two discrete approximations with different time-steps, a
fine and a coarse, but the same Brownian motion paths. This is done with various selections of simulations and time-steps for each level. The magnitude depended on the strong convergence properties of the scheme. Let 𝑉ℓ denote the variance of a single sample ˆ𝑃
(𝑖) ℓ − ˆ
𝑃(𝑖)
ℓ−1. The
Complexity Theorem (see Theorem 7) shows that the convergence rate of 𝑉𝑙 as 𝑙 → ∞
determines the efficiency of the Multilevel approach. To ensure a better efficiency, (4.6) may be modified, and different estimators of ˆ𝑃could be used on the fine and coarse discretization levels of ˆ𝑌ℓ, 𝐸[ ˆ𝑃ℓ] = 𝐸 [ ˆ𝑃0] + 𝐿 Õ ℓ=1 𝐸 h ˆ 𝑃 𝑓 ℓ − ˆ𝑃 𝑐 ℓ−1 i , where ˆ𝑃 𝑓 ℓ,𝑃ˆ 𝑐
ℓ−1 are the estimators using 2 ℓ
and 2ℓ−1
steps respectively in the computation of ˆ 𝑌ℓ. Provided that 𝐸[ ˆ𝑃 𝑓 ℓ] = 𝐸 [ ˆ 𝑃𝑐 ℓ], (4.7)
the telescoping sum property is maintained. To implement this, the Brownian motion incre-ments are first constructed for the simulation of the discrete path leading to the evaluation of the term ˆ𝑃(𝑖)
ℓ , continuing to sum them in groups of size 𝑀 to give the discrete Brownian motion
increments for the evaluation of ˆ𝑃(𝑖)
ℓ−1. The variance of this estimator is 𝑉 [ ˆ
𝑌ℓ]=𝑁−1
ℓ 𝑉ℓ, where
𝑉ℓ is the variance of a single sample. Using stratified sampling or a zero-mean control variate
to reduce the variance, the same inverse dependence on 𝑁ℓ would apply. The Variance of the
combined estimator ˆ𝑌 =Í𝐿 ℓ=0𝑌ˆℓ is, 𝑉[ ˆ𝑌] = 𝐿 Õ ℓ=0 𝑁−1 ℓ 𝑉ℓ.
The cost of computation is equal to
𝐿 Õ ℓ=0 𝑁ℓℎ −1 ℓ .
Assume that 𝐿 1 and consider the behaviour of 𝑉ℓ as 𝑙 → ∞. In the particular case of the
Euler discretization and the Lipschitz payoff function, provided that the term 𝑎(𝑆, 𝑡) and 𝑏(𝑆, 𝑡) satisfy certain conditions, there is a 𝑂 (ℎ) weak convergence and 𝑂 (ℎ1/2) strong convergence. Hence, as ℓ → ∞, 𝐸 𝑃ˆℓ − 𝑃 = 𝑂 (ℎℓ), (4.8) and 𝐸 || ˆ𝑆 ℓ, 𝑀ℓ − 𝑆(𝑇 ) || 2 = 𝑂 (ℎ ℓ). (4.9)
From the Lipschitz property (3.2), it follows that
V[ ˆ𝑃ℓ− 𝑃] ≤ 𝐸 [( ˆ𝑃ℓ − 𝑃)2] ≤ 𝑐2𝐸[|| ˆ𝑆
ℓ, 𝑀ℓ− 𝑆(𝑇 ) ||2]. (4.10)
Combining (4.10) with (4.9) gives 𝑉 [ ˆ𝑃ℓ− 𝑃] = 𝑂 (ℎℓ). Further,
( ˆ𝑃ℓ − ˆ𝑃ℓ−1) = ( ˆ𝑃ℓ − 𝑃) − ( ˆ𝑃ℓ−1− 𝑃) −→ V[ ˆ𝑃ℓ − ˆ𝑃ℓ−1] ≤ (V[ ˆ𝑃ℓ− 𝑃])1/2+ (V[ ˆ𝑃ℓ−1− 𝑃])1/2 2 .
Hence, for the simple estimator (4.6), the single sample variance 𝑉ℓ is 𝑂 (ℎℓ) and the optimal
choice for 𝑁ℓ is asymptotically proportional to ℎℓ. Setting 𝑁ℓ = 𝑂 (𝜖 −2𝐿 ℎ
ℓ), the variance of
the combined estimator ¯𝑌 is 𝑂 (𝜖2). If 𝐿 is now chosen such that
𝐿 =
log𝜖−1
log𝑀 + 𝑂 (1), (4.11)
as 𝜖 → 0, then ℎ𝐿 = 𝑀
−𝐿 = 𝑂 (𝜖 ). Therefore the bias error 𝐸 [ ˆ
𝑃𝐿 − 𝑃] is 𝑂 (𝜖), due to (4.6). Consequently, a MSE which is 𝑂 (𝜖2) are obtained, with a computational complexity which is 𝑂(𝜖−2𝐿2) = 𝑂 (𝜖−2(log𝜖)2) [18].
4.3.1
Complexity Theorem
A way to compare the computational cost of the standard Monte Carlo with the Multilevel Monte Carlo is by computational complexity theory. This is the cumulative number of times the random number generator is called for a certain method. The cost of regular Monte Carlo is proportional to the number of paths multiplied by the cost of simulating a single path. To measure the computational complexity of the standard Monte Carlo process, an exact analytical solution to an SDE is rarely feasible in most applications. This means one has to use a combination of regular Monte Carlo and a numerical scheme for solving SDEs, such as the Euler-Maruyama scheme or Milsteins schemes. Let ˆ𝑌 be the numerical scheme’s approximation of 𝑌 . The mean square error gives us the equation:
MSE ≡ 𝐸 [( ˆ𝑌 − 𝑌 )2]. (4.12)
Rewriting equation (4.12) yields
𝐸[( ˆ𝑌 − [𝑌 ])2] = 𝐸 [( ˆ𝑌+ 𝐸 [ ˆ𝑌] − 𝐸 [ ˆ𝑌] − 𝑌 )2)
= 𝐸 [𝐸 [( ˆ𝑌 − 𝐸 [ ˆ𝑌])2) + 𝐸 [(𝐸 [ ˆ𝑌] − 𝑌 )2],
given the first term to be the statistical error of the Monte Carlo variance. Since when approximating a time continuous process with a discrete process, the second term is the bias of the numerical discretization scheme [21].
It is often difficult to precisely simulate the random variable of interest, in this case, the payoff 𝑃 in certain applications. This is often the case when an application that involves the simulation of SDEs is performed. Consider 𝑌 to be an approximation to 𝐸 [𝑃], the mean square error is given by
MSE(𝑌 ) = 𝐸 [(𝑌 − 𝐸 [𝑃])2] = Var[𝑌 ] + (𝐸 [𝑌 ] − 𝐸 [𝑃])2. (4.13)
Let ˆ𝑌 be the Multilevel estimator
ˆ 𝑌 = 𝐿 Õ 𝑙=0 ˆ 𝑌𝑙, where 𝐸[ ˆ𝑌] = 𝐸 [𝑃𝐿], Var[ ˆ𝑌] = 𝐿 Õ 𝑙=0 𝑁−1 𝑙 𝑉𝑙.
From this, the mean square error in (4.13) can be written as
𝐸[( ˆ𝑌 − 𝐸 [𝑃])2] = Var[ ˆ𝑌] + (𝐸 [ ˆ𝑌] − 𝐸 [𝑃])2=Var[ ˆ𝑌] + (𝐸 [𝑃𝐿 − 𝑃])2. This leads to the following complexity theorem [11].
