Explicit cubature method in Heston model

School of Education, Culture and Communication

Division of Applied Mathematics

BACHELOR THESIS IN MATHEMATICS / APPLIED MATHEMATICS

Explicit cubature method in Heston model

by

Michaela Törnqvist and Abdukayum Sulaymanov

Kandidatarbete i matematik / tillämpad matematik

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

Date:

2017-06-02

Project name:

Explicit cubature method in Heston model

Authors:

Michaela Törnqvist and Abdukayum Sulaimanov

Supervisor: Anatoliy Malyarenko Reviewer: Milica Ranˇci´c Examiner: Ying Ni Comprising: 15 ECTS credits

Abstract

In this paper we demonstrate one of the most recent methods in financial engineering called Cubature on Wiener space developed by Lyons and Victoir and its applications for Heston model. This model is theoretically a contender to classical models that relies on Monte Carlo simulation method, and has accuracy up to degree 5 of iterated Stratonovich integrals of Wiener process. Contribution of this paper could be seen in the solving of ODEs in finance in order to get the price of contingent claim.

Contents

Acknowledgements . . . 4

1 Introduction 5 1.1 The formulation of the problem . . . 5

1.2 Review of literature . . . 6

2 Theoretical considerations 10 2.1 Important definitions in finance . . . 10

2.2 The Riemann–Stieltjes Integral . . . 11

2.2.1 The sufficient condition of Riemann–Stieltjes integrability . . . 13

2.3 Banach Space . . . 13

3 The solution of the problem 15 3.1 Stochastic Cubature Formulae . . . 15

3.2 Heston model . . . 19

3.3 Conclusion . . . 24 A Criteria of Bachelor Thesis 26

List of Figures

1.1 Ω space . . . 7 3.1 ω₁[2](t) path . . . 24

List of notations

Here we give some explanatory descriptions of the mathematical symbols used throughout this paper.

SDE : stochastic differential equation ODE : ordinary differential equation

r : interest rate T : time of maturity t : current time S(t) : price of stock Ω : sample space F : filtration

P∗ : risk neutral probability measure W(t) : Wiener process for t ≥ 0

Acknowledgements

First and foremost, we would like to show our sincere gratitude to our research supervisor, Pro-fessor Anatoliy Malyarenko for his dedicated support, assistance and continuous encourage-ment throughout the process of writing this thesis. The time our professor took by reschedul-ing parts of his day in order to provide us with unfailreschedul-ing support and immense knowledge is something we appreciated dearly.

We would also like to express our thankfulness to our reviewer, Dr. Milica Ranˇci´c and our examiner, Dr. Ying Ni for their valuable comments on this thesis. By providing us with constructive criticism, they allowed us to stay confident and feel enthusiastic about the im-provements of our thesis.

Throughout our studies at the Mälardalen University, we have been privileged with pas-sionate and encouraging professors and teachers. We would especially like to express our gratitude to Dr. Lars-Göran Larsson, for his valuable advice and kindness.

Finally, we cannot begin to express our genuine gratitude to our student counselor, Anna-Clara Meurling Karlsson, whom has throughout the years at the university been an invaluable academic support and showed devoted involvement in our studies.

Chapter 1

Introduction

1.1

The formulation of the problem

We will throughout this paper refer to Kijima (2013) for necessary mathematical definitions which are mainly presented in Chapter 2.

Financial Engineering is an area of applied mathematics that studies market models. The main objective of this thesis is to price a contingent claim, more precisely, we would like to price a European call option using a cubature method in Heston model.

Let (Ω, F, P) be a probability space, B(t) = ertbe the risk-free bank account and let S(t) be an Rm-valued price process. Assume that the market is complete, that is, all contingent claims are attainable at maturity. Now, since the market is complete, there exists a unique probability measure P∗equivalent to P such that the discounted price processes

S_i∗(t) = e−rtSi(t), 1 ≤ i ≤ m, 0 ≤ t ≤ T

are martingales under P∗. Finally, let X = X (ω) be a contingent claim with maturity T . In words, the fundamental theorem of financial engineering says: the no-arbitrage price of X is equal to its discounted expected payoff under the risk-neutral measureP∗. Mathematically,

C= e−rTE∗[X ]. By definition of the mathematical expectation,

E∗[X ] = Z Ω X(ω) dP∗(ω), therefore C= e−rT Z Ω X(ω) dP∗(ω).

The main objective of this thesis is as follows: instead of calculating the integral over a complicated space Ω with respect to a complicated measure P∗, we calculate an integral over a “simple” subset Ω0⊂ Ω with respect to a “simple” measure Q. In particular, the simplest possible subset is a finite one:

Denote λi= Q({ωi}), 1 ≤ i ≤ N. (1.1) Then we have C≈ e−rT N

∑

i=1 λiX(ωi). (1.2)

Two natural questions arise:

1. If we have two formulae of type (1.2), which is better? 2. What is the “best” approximation, and how to find it?

1.2

Review of literature

To explain the main idea, consider the simplest case when Ω = R, and for an event A P∗(A) = √1

2π

Z

A

e−x2/2dx, that is, the standard normal distribution. To calculate the integral

Z ∞

−∞

f(x) dP∗(x),

assume that the function f is smooth. By Taylor expansion we get f(x) ≈ f (0) + f0(0)x + f 00₍₀₎ 2! x 2_{+ · · · +} f(n)(0) n! x n_.

Assume we found some points x1, . . . , xN and weights in equation (1.1) such that for any

polynomial P(x) of degree 0, 1, . . . , n we have

Z _∞ −∞ P(x) dP∗(x) = N

∑

i=1 λiP(xi). (1.3)

Then we may expect that the absolute error      Z ∞ −∞ f(x) dP∗(x) − N

∑

i=1 λif(xi)      is small at least for some “good” functions f (x).

This is indeed the case. A formula

Z ∞ −∞ f(x) dP∗(x) ≈ N

∑

i=1 λif(xi)

Figure 1.1: Ω space

is called a quadrature formula of degree n, if Equation (1.3) holds true for all polynomials of degrees up to n. Similarly, when d > 1, a formula

Z Rd f(x) dP∗(x) ≈ N

∑

i=1 λif(xi)

is called a cubature formula of degree m, if for all polynomials P(x) of degree 0, 1, . . . , n we have Z Rd P(x) dP∗(x) = N

∑

i=1 λiP(xi).

One of standard references in the theory of cubature formulae is Stroud (1971). The main result of this theory is the so called Tchakaloff theorem.

Denote by Rm[X1, . . . , Xd] the set of all polynomials of degree m in d variables with real

coefficients.

Definition 1. Let µ be a positive measure on Rd, and m be a natural number. We define the points x1, . . . , xnto be in the support of µ, and let the positive weights λ1, . . . , λndefine a

cubature formula of degree m with respect to µ if and only if for all P ∈ Rm[X1, . . . , Xd],

Z Rd P(x)µ(dx) = n

∑

i=1 λiP(xi). (1.4)

Let us now refer to an important theorem (Tchakaloff (1957)) by restating Lyons and Victoir (2004):

Theorem 1. Let m be a positive integer and let µ be a positive measure on Rdwith the property that R

|P(x)|µ(dx) < ∞ and ∀P ∈ Rm[X1, . . . , Xd]. One can then find n points x1, . . . , xn∈ Rd

and n positive real numbers λ1, . . . , λn, with

such that the cubature relation holds in equation_{(1.4) for all P ∈ Rm}[X1, . . . , Xd].

Remark 1. The Tchakaloff theorem is an existence theorem. This means that its proof does not provide any method to construct a cubature formula. In general, constructing cubature formulae is a very complicated task (Stroud (1971)).

Now, we would like to realise the above described ideas in the “complicated” space Ω = C_{0([0, T ]; R}d) of all continuous functions ωω_{ω : [0, T ] → R}d with ωωω (0) = 0. The first question is as follows: How does the Taylor expansion look like in this space? The answer was given by Platen and Wagner (1982), see also Kloeden and Platen (1991) and Kloeden and Platen (1992). In fact, as we will explain below, there are two variants of the so called stochastic Taylor expansion: the first includes iterated Itô integrals, while the second includes iterated Stratonovich integrals. The latter variant is simpler than the former one. The reason is as follows: the stochastic Taylor expansion is obtained by iteration of the Itô formula. In the case of Itô integrals, the Itô formula contains the term with the second derivative, while in the case of Stratonovich integration, the above formula is completely similar to the classical chain rule and contains only the first derivative. In fact, the iterated Statonovich integrals play the same rôle in the space C0([0, T ]; Rd), as polynomials play in the space Rd.

A condition similar to (1.4) was first found by Kusuoka (2001). Later, Lyons and Vic-toir (2004) realised that Kusuoka’s condition is related to cubature formulae in Wiener space C0([0, T ]; Rd). They proved the following extension of the Tchakaloff theorem:

Let k be an nonnegative integer, and let ααα be either the empty set if k = 0 or a multi-index α

αα = (α1, . . . , αk) with 0 ≤ αi≤ d. Define the number kααα k as k plus the number of zeros among

αi’s and call this the degree of ααα . Let dm be the number of multi-indices ααα with kααα k ≤ m, including the empty set. Denote W0(t) = t. The Tchakaloff theorem then takes the following

form:

Theorem 2 (Lyons and Victoir (2004)). For any positive integer m, there exists N ≤ dmpaths

ω

ωω1, . . . , ωωωn of bounded variation, and N positive weights λ1, . . . , λn that define a cubature

formula on Wiener space of degree m at time T . That is, for all multi-indices ααα = (i1, . . . , ik)

withkααα k ≤ m we have E Z 0<t1<···<tk<T ◦dWi1(t1) ◦ · · · ◦ dWik(tk)  = N

∑

j=1 λj Z 0<t1<···<tk<T dωi1_j (t1) · · · dωi_jk(tk), (1.5)

where the integral in the left hand side is an the iterated Stratonovich integral, and the integral in the right hand side is the iterated Riemann–Stiltjes integral.

Remark 2. Again, this theorem is an existence theorem. To find examples of cubature for-mulae, Lyons and Victoir (2004) use complicated mathematical tools which will not be under consideration in this thesis.

Theorem 3. Let’s assume that µ satisfies the moment matching condition (1.5) and f is infin-itely differentiable with all derivatives bounded. Then there exists constant C that depends on m and f such that

     E∗[ f (X (t))] − N

∑

j=1 λjf(ωωωj)      ≤ Ct(m+1)2 .

This Theorem is one of cornerstones of our paper because we rely on this in our future calculation as target of fulfillment.

In Chapter 2 we will explain more in detail the mathematical definitions necessary for understanding Theorem 2. In Chapter 3 we will explain how to apply Theorem 2 to the solution of real-life problems, and consider the pricing of a European option under Heston model as an example.

Chapter 2

Theoretical considerations

2.1

Important definitions in finance

Definition 2 (σ -field). A collection of subsets of a set S is called a σ -field, denoted by B, if and only if 1. ∅ ∈ B, 2. if A ∈ B then S \ A ∈ B, 3. if A1, A2, . . . , ∈ B then ∞ S i=1 A_i∈ B.

Definition 3 (The Borel σ -field). The Borel σ -field of subsets of R, written B, is the smallest σ -field generated by the open sets. That is, ifO denotes the collection of all open subsets of R, then B = σ (O).

Definition 4 (Filtration). A filtration {Ft} is a set of σ -fields satisfying:

Ft ⊂F , f or m ≥ 0,

Fs⊂Ft, f or s≤ t.

Definition 5 (Martingale). Let {X (t)}, 0 ≤ t ≤ T be an integrable stochastic process on a probability space (Ω, F, P) with a filtration {Ft}, that is, E[|X(t)|] < ∞ for all t ∈ [0, T ]. If the

process satisfies

E[X (t) | Fs] = X (s), s≤ t,

it is then called a martingale.

We will now explain what we mean by “complicated space” and “complicated measure”. Definition 6 (Wiener process). Let the probability space be (Ω, F, P) and consider a stochastic process {W (t),t ≥ 0} that is defined on the probability space. Then W (t) is called standard Brownian Motion or Wiener process if and only if

2. the increments are normally distributed with parameters µ = 0 and σ2= s, for 0 ≤ t ≤ t+ s < ∞,

3. W (0) = 0, with continuous sample paths.

Definition 6 is highly imperative and utile for modern finance and is because of that used in many calculations and algorithms.

Definition 7 (d-dimensional Wiener process). This is the process W(t) = (W1(t),W2(t), . . . ,Wd(t))>,t ≥ 0,

where W1(t), . . . , Wd(t) are independent Wiener processes.

We can also refer to the following uniqueness property of Wiener measure (Carkovs et al. (2011)).

Theorem 4. There exist a unique probability measure P on the Borel σ -field of the space Ω = C0([0, T ]; Rd) such that the coordinate mapping process

W(t, ωωω ) = ωωω (t), 0 ≤ t ≤ T is a d- dimensional Wiener process.

Now, if we again consider the price of the contingent claim C= e−rTE∗[X ] = e−rT

Z

Ω

X(ωωω )dP∗(ωωω )

we see that we now can rely on Theorem 4 and write the price of our contingent claim X in integral form.

However, is it possible for us to calculate this integral precisely or should we attempt to simplify our equation further? To answer this we try to apprehend the essence of Wiener space thoroughly (see figure (1.1), where ω- function and λ -weight) and quickly realise that it is indeed complicated to calculate but possible.

2.2

The Riemann–Stieltjes Integral

The Riemann–Stieltjes Integral is one of many widely used integrals in probability and was named after Dutch mathematician Stieltjes and German mathematician Riemann (Carkovs et al. (2011), Rudin (1976)).

Definition 8. A partition P in the given interval [c, d] is defined as a finite set of points x₀, x1. . . , xnwith

Definition 9. A set D ⊂ R is bounded from above if there exists T ∈ R, called an upper bound of D, {x : x6 T } for every x ∈ D. By similar logic, D is bounded from below if there exists l_{∈ R, called a lower bound of D, {x : x > l} for every x ∈ D. So, a bounded set is bounded} from below and from above.

Let us assume that we have a bounded real function f (x) that is smooth on the closed interval [c, d], and then let us introduce one more definition.

T_k= sup xk−16x6xk f(x), l_k= inf xk−16x6xkf(x), U(P, f ) = n

∑

k=1 T_k∆xk, L(P, f ) = n

∑

k=1 l_k∆xk,

where ∆x_k = x − k − x_k−1. We now have a formal definition of upper and lower Riemann Integrals in integral form, which we usually study in the course of calculus,

Z d c f(x) dx = inf P U(P, f ), Z d c f(x) dx = sup P L(P, f ).

It is called Riemann-integrable on [c, d], iff the lower and upper integrals are equal. We write f ∈R, where R denotes the set of Riemann integrable functions. We can visualize it by looking at figure 2.

Z d

c

f(x)dx.

Now, when we understand the trivial Riemann integral, let us consider the general case with the new definition.

Definition 10. Let a function β be defined on [c, d]. It is called a function of bounded variation (BV) if and only if there exists M∈ R such that for all partitions of [c, d] we have

n−1

∑

j=0 |β (xj+1) − β (xj)| ≤ M. Denote ∆_βx_k= β (xk+1) − β (xk) and U(P, f , β ) = n

∑

k=1 T_k∆_βx_k, L(P, f , β ) = n

∑

k=1 l_k∆_βxk.

Now, for all real functions f that are bounded on the interval[c, d], define the upper and lower Riemann–Stiltjes integrals by

Z d c f(x)dβ (x) = inf P U(P, f , β ) Z d c f(x)dβ (x) = sup P L(P, f , β ).

A function f will be called Riemann–Stieltjes integrable when the upper and lower integ-rals are equal. Their common value is called the Riemann–Stiltjes integral of f with respect to β .

2.2.1

The sufficient condition of Riemann–Stieltjes integrability

Theorem 5 (Jordan). A function f has bounded variation on the interval [c, d] if and only if it can be written as the difference of two monotonic non-decreasing functions.

Theorem 6. If f is continuous on the closed interval [c, d], then it is integrable w.r.t β on the interval[c, d]. If f is monotonic and β is continuous on the same closed interval [c, d], then f is integrable w.r.t β on the interval [c, d].

2.3

Banach Space

Banach space will be considered in the Section 3.1 and denoted C0,BV([0, T ]; Rd+1) and here

we will attempt to give some nontrivial concept of this space by referring to (Rudin (1991)). Definition 11. A vector space over R is a special set V that contains 0 and satisfies the scalar multiplication- and addition conditions below:

1. u + v = v + u{∀u, v ∈ V } 2. ∃0 : 0 + u = 0, ∀u ∈ V

3. (v + u) + w = v + (u + w), {∀u, v, w ∈ V } 4. For all u ∈ V,∃ (-u) ∈ V: u + (-u) = 0 5. λ (ku) = (λ k)u, {k, λ ∈ R, u ∈ V } 6. λ (u + v) = λ u + λ v

Definition 12. Let S∗ be vector space over field of complex scalars, then the norm on the vector space S∗is the function of ||·|| such that S∗_{→ R, satisfy below properties:}

1. ||x|| ≥ 0 for all x ∈ S∗, and ||0|| = 0 if and only if x = 0. 2. ||ix|| = |i|||x||, for all x ∈ S∗and i ∈ C.

3. ||x + y|| ≤ ||x|| + ||y||, for all x, y ∈ S∗.

Definition 13. Let us assume that {xn}n∈Nis a series of elements of S∗. Then,

limx→∞kx − xnk = 0 if for every ε > 0 we can find N ∈ N: kx − xnk < ε, for all n ≤ N.

Definition 14. If for every ε > 0 there is an N ∈ N : kxm− xnk < ε for all m, n ≥ N, it is a

Cauchy sequence.

Definition 15. A normed vector space S with property that all Cauchy sequences are conver-gent is complete. A complete normed vector space is called a Banach space.

Chapter 3

The solution of the problem

In this Chapter, we formulate and partially solve the problem of pricing European call option in Heston model by the advanced theory explained above. First, we planned to chose Black– Scholes model instead. However, we found that the stochastic Taylor formula for this model is trivial.

On the other hand, when we applied our theory to the Heston model, we found that about 50 systems of differential equations had to be solved. We solved explicitly 4 of them to demonstrate the method. The remaining systems can be solved similarly. Another important remark here, we used the works of Lyons and Victoir (2004) where they used trick and replaced complicated dW∗ by ωi- piecewise linear continuous path and they also calculated the slope

of this function θi, j, all values could be found on the Appendix of this paper.

3.1

Stochastic Cubature Formulae

A stochastic integral (Itô process) is a stochastic process on the probability space (Ω, F, P). Taylor expansions on stochastic integrals is a tool for better understanding how stochastic integration schemes are constructed, See: Canhanga et al. (2017).

Let’s consider as market model:

dX_t= µ(t, X(t))dt | {z } deterministic part +

_∑

(t, X(t))dW∗(t) | {z } stochastic part , X(0) = X0 |{z}

initial price at time zero

(3.1) where,

X(t) : [0, T ] → Rnis a stochastic process µ : [0, T ] × Rn→ Rnis the drift

∑ : [0, T ] ×Rn→ Rn×d is a Wiener process under P∗∈ (Ω, , F).

an integral form X(t) = X (0) + Z t 0 µ (X (s))ds + Z t 0 σ (X (s))dW∗(s), (3.2) where we assume µ, σ ∈ C∞ _{and satisfy the linear growth bound. By generating the Taylor}

expansion for any function f (x) ∈ C2we can use the Itô lemma. Thus f(X (t)) = f (X (0)) + Z t 0 L 0_{f(X (s))ds +}Z t 0 L 1_{f(X (s))dW}∗_{(s), (3.3)} defining L0_{= µ} ∂ ∂ x+ 1 2σ 2 ∂2 ∂ x2, (3.4) L1_{= σ} ∂ ∂ x.

By applying the Itô to the functions f = µ and f = σ in (2) we obtain the simplest non-trivial Itô–Taylor expansion

X(t) = X (0) + µ(X (0)) Z t 0 ds+ σ (X (0)) Z t 0 dW(s)∗ + Z t 0 Z s 0 L 0 µ (X (u))du · ds + Z t 0 Z s 0 L 1 µ (X (u))dW∗(u)ds + Z t 0 Z s 0 L 0 σ (X (u))du · dW∗(s) + Z t 0 Z s 0 L 1 σ (X (u))dW∗(u) · dW∗(s) (3.5)

If we continue to expand the equations above to higher orders, we will instantly see that they will become very complicated to solve because of second derivative ofL0. So instead we will try to mitigate this obstacle by changing the equation (3.2) for one-dimensional Stra-tonovich SDE in the integral form.

X(t) = X (0) + Z t 0 ˜ µ (X (s))ds + Z t 0 σ (X (s)) ◦ dW∗(s), where the one-dimensional Stratonovich correction ˜µ has the form

˜

µ = µ −1 2σ

2_.

Then we will see that our equation (3.3) takes the form f(X (t)) = f (X (0)) + Z t 0 ˜ L0_{f(X (s))ds +} Z t 0 L 1_{f(X (s)) ◦ dW}∗_(s), where ˜ L0_{= ˜µ} ∂ ∂ x.

The new non-trivial Stratonovich–Taylor expansion will be X(t) = X (0) + ˜µ (X (0)) Z t 0 ds+ σ (X (0)) Z t 0 ◦dW (s)∗ + Z t 0 Z s 0 ˜ L0_{µ (X (u))duds +˜} Z t 0 Z s 0 L 1_˜ µ (X (u)) ◦ dW∗(u)ds + Z t 0 Z s 0 ˜ L0_{σ (X (u))du ◦ dW}∗_{(s) +} Z t 0 Z s 0 L 1 σ (X (u)) ◦ dW∗(u) ◦ dW∗(s). Before proceeding further let’s define important difference between Itô and Stratonovich in-tegral, Kijima (2013) Let’s consider a stochastic integral

I(t) =

Z t

0

ψ (u)dW (u), 0 ≤ t ≤ T, where ψ is a stochastic process with continues paths.

I(t) could be calculated by uniformly dividing the interval into m subintervals , where 0 = t0< t1· · · < tm= t,tj≡

t mj. Now let’s approximate I(t) as Riemann–Stieltjes sum

Jm≡ n−1

∑

j=0 ψ (tj) + ψ(tj+1) 2 {W (tj+1) −W (tj)} (3.6) or Riemann–Stieltjes sum could also be defined as:

I_m≡

n−1

∑

j=0

ψ (tj){W (tj+1) −W (tj)}. (3.7)

However, when n → ∞, limits of equations (3.6) and (3.7) differ. Equation (3.6) called Stra-tonovichintegral and (3.7) called Itô integral are directly connected to martingales and therefor appropriate approximations in finance.

Our market model in the Stratonovich form will be X(t) = d

∑

i=0 Z t 0 Vi(X(s)) ◦ dWi∗(s) (3.8)

where W₀∗(s) = s and where

Vi₀(y) = ˜Vi₀(y) −1 2 d

∑

j=1 N

∑

k=1 V_jk(y)∂V i j ∂ yk (y) is the multidimensional Stratonovich correction. Then

f(Xt(x)) = f (x) + d

∑

i=0 Z t 0 (Vif)(Xs(x)) ◦ dWs∗i,

where (Vif)(y) := N

∑

j=1 V_ij(y)∂ f ∂ yj (y). For k ≥ 1, define the multiple Stratonovich integral by

I(t, ααα , ◦dW∗) := Z t 0 · · · Z t_k−2 0 Z t_k−1 0 ◦dW_tk∗αk◦ dW∗αk−1 tk−1 ◦ · · · ◦ dW ∗α1 t1 . Then f(Xt(x)) =

∑

kααα k≤m I(t, ααα , ◦dW∗)(Vαk· · · Vα1f)(x) + R X m, (3.9)

where x = x0and the remainder RXmcontains multiple integrals of degrees greater than m.

Let C0,BV([0, T ]; Rd+1) be the space of functions ωωω that are Rd+1-valued with ωωω (0) = 0 whose components have bounded variation on [0, 1]. Let ν be a probabilistic measure on this space. Along with the system (3.8), consider the system

Xt(x) = x + d

∑

i=0 Z t 0 Vi(Xs(x)) dωsi,

that is, the system of random ordinary differential equations. Let Xωω_tω(x) be its solution. Define the time-scaled path ωωωs[t] : Ω → C0,BV([0,t]; Rd+1) by

ω_si[t] = (

tω_s/t0 , if i = 0, √

tω_s/ti , otherwise. We would like to estimate the approximation error

|E∗[ f (Xt(x))] − Eν[ f (Xω ω ω [t] t (x))]|.

For f (Xωω_tω [t](x)) we have an equation similar to (3.9): f(Xωω_tω [t](x)) =

_∑

kααα k≤m I(t, ααα , dωωω [t])(Vαk· · · Vα1f)(x) + R X m, where I(t, ααα , dωωω [t]) := Z t 0 · · · Z t_k−2 0 Z t_k−1 0 dω_tαk_k dω_tαk−1_k−1 · · · dωα1 t1 .

The approximation error becomes |

_∑

kααα k≤m E∗[I(t, ααα , ◦dW∗)](Vαk· · · Vα1f)(x) + E ∗_[RX m] −

_∑

kααα k≤m Eν_{[I(t, αα} α , dωωω [t])](Vαk· · · Vα1f)(x) − E∗[R X m]|. (3.10) Kusuoka (2001) proposed the following idea. If the first and third terms in the above formula cancel each other, then the remaining approximation error

|E∗[RX_m] − E∗[RX_m]| is small.

Definition 16 (The moment matching condition, Kusuoka (2001)). The measure ν satisfies the moment matching condition of order m if for all kααα k ≤ m we have

E∗[I(t, ααα , ◦dW∗)](Vαk· · · Vα1f)(x) = Eν[I(t, ααα , dωωω [t])](Vαk· · · Vα1f)(x).

It is easy to check that the sums cancel each other. To discuss the remaining error, consider the simplest case, when the measure ν is supported by finitely many points ωωω1, . . . , ωωωN with

ν (ωωωi) = λi. Then Eν[RXm] can be calculated exactly as a finite sum. This case is called the

cubature on Wiener space, see Lyons and Victoir (2004) and our discussion above. The error is given by Theorem 3 and is indeed small for small values of t.

As we mentioned above we would like to show our justification that application of the stochastic Taylor formula on Black–Scholes formula gives trivial case. So, let’s rewrite our Black–Scholes formula in integral form first:

S(t) = S(0) + Z t 0 rds+ Z t 0 σ dWs (3.11)

Now, we refer to equation on See: Kloeden and Platen (1992) and expand our Black–Scholes equation (3.11) by applying (3.5) we get

S(t) = S(0) + r Z t 0 ds+ σ Z t 0 dW_s+ R (3.12) Now, one can see that as equations (3.11) and (3.12) are equal, R = 0.

3.2

Heston model

Heston model consist of two stochastic differential equations, the first reflects changes in price. See: Carkovs et al. (2011)

dS(t) = µS(t)dt +pV(t)S(t)dW₁∗(t), (3.13) and the second stochastic differential equation that reflects changes in volatility of the stock

dV(t) = [κ(θ −V (t))]dt + σpV(t)dW₂∗(t),

where, (3.14) κ - reversion rate,

θ - long-run mean, V(t)- variance

Now we assign our vector to be

X = S(t) V(t) 

Σ =pV (t)S(t) 0 0 σpV(t)

 and so, our drift in the model could be denoted

µ µµ =  µ S(t) κ (θ − V (t))  . Now, we can write our Heston model in integral form

 S(t) V(t)  = S(0) V(0)  + 2

∑

i=0 Z t 0 V_i(S(t),V (t))dW_i∗(s) where ˜ V0=  µ S(t) κ (θ − V (t))  , V1= pV (t)S(t) 0  and V2=  0 σpV(t)  . (3.15) As we see above, equation in the form of Itô are very inconvenient for us to work with, so by Stratonovich correction, where Vi, 1 ≤ i ≤ d our Stratonovich correction will be in the

form Vi₀(y) = µj(y) − 1 2 2

∑

i=1 2

∑

k=1 Vk_i∂ V j i ∂ yk , 1 ≤ j ≤ 2.

here let us first define our V1₀and apply the Stratonovich correction formula to (3.15), V1₀=µS(t) −1 2 2

∑

i=1 2

∑

k=1 V_ik∂V 1 i ∂ yk =µS(t) −1 2  V₁1∂V 1 1 ∂ S +V 2 1 ∂V₁1 ∂V +V 1 2 ∂V₂1 ∂ S +V 2 2 ∂V₂1 ∂V  =µS(t) −1 2 p V(t)S(t)pV(t) + 0 + 0 + σpV(t)0 =µS(t) −1 2S(t)V (t) (3.16)

and V2₀in the Stratonovich form will be V2₀=κ(θ −V (t)) −1 2 2

∑

i=1 2

∑

k=1 V_ik∂V 2 i ∂ yk =κ(θ −V (t)) −1 2  V₁2∂V 2 1 ∂ S +V 2 1 ∂V₁2 ∂V +V 2 2 ∂V₂2 ∂ S +V 2 2 ∂V₂2 ∂V  =κ(θ −V (t)) −1 2  0 + 0 + 0 +1 2σ p V(t)σp1 V(t)  =κ(θ −V (t)) −1 4σ 2_. (3.17)

Now, let us consider Heston model in Itô form:  S(t) V(t)  = S(0) V(0)  + Z t 0  µ S(s) κ (θ − V (s))  ds+ Z t 0 pV (s)S(s) 0  dW₁∗(s) + Z t 0  0 σpV(s)  dW₂∗(s). (3.18) However, as we know so far, Stratonovich form is more convenient for us to work with and thus we rewrite the above model in Stratonovich form by the result which we got from (3.16) and (3.17)  S(t) V(t)  = S(0) V(0)  + Z t 0  µ S(s) −1₂S(s)V (s) κ (θ − V (s)) −1₄σ2  ◦ ds+ Z t 0 pV (s)S(s) 0  ◦ dW₁∗(s) + Z t 0  0 σpV(s)  ◦ dW₂∗(s). (3.19) We will now use a trick by which we will bring it to the form of an usual Riemann integral and consider the work of Lyons and Victoir (2004) where they simplified the two stochastic part in our equation(dW₁∗(s), dW₂∗(s)) by applying discrete path

 S(t) V(t)  = S(0) V(0)  + Z t 0  µ S(s) −1₂S(s)V (s) κ (θ − V (s)) −1₄σ2  ◦ ds + Z t 0 pV (s)S(s) 0  ◦dω [0] 1 (s) ds ds+ Z t 0  0 σpV(s)  ◦ √ 3dω₂[2](s) ds ds. (3.20) Now , let’s solve this equation. We will try to calculate the volatility by taking the derivative from equation (3.20) with respect to time in order to get rid off all integrals.

V0= κ(θ −V (t)) −1 4σ 2_{+ 4θ} i, jσ p V(t), κ > 0, θ > 0. Let’s denote _√ V = u V0= 2uu0 Insertion yields 2uu0 κ (θ − u2) −1₄σ2+ 4θi, jσ u = 1 2uu0 −κu2₋4θi, jσ κ u− (θ + −1 4 σ2 κ )  = 1 2uu0 u−2θi, jσ κ | {z } α∗>0
2 − θ + −1 4 σ 2 κ + 4σ2θ_{i, j}2 κ2
| {z } β2>0 = −κ (3.21)

We want to clarify here that the value of β2> 0. Then 2uu0 (u − α∗)2_{− β}2 = −κ α∗+β β u− (α∗_{+ β )}+ α∗+β −β u− (α∗_{− β )}
u 0_{= −κ}

Now we can integrate both sides α∗+ β β ln |u − (α ∗_{+ β )| −}α∗− β β ln |u − (α ∗_{− β )| = −κt + eC,} ln |u − (α ∗_{+ β |}αβ∗+1 |u − (α∗_{− β )|}α ∗β −1  = −κt + eC, where Ce= constant we exponentiate exp  ln |u − (α ∗_{+ β |}αβ∗+1 |u − (α∗_{− β )|}α ∗β −1  = exp( eC) exp(−κt),  |u − (α∗_{+ β |}αβ∗+1 |u − (α∗_{− β )|}α ∗β −1  = τ exp(−κt), where τ = constant we then exponentiate by β on both side

|√V_t− (α∗+ β |α∗+β

|√V_t− (α∗_{− β )|}α∗−β = C2exp(−κtβ ),

where C₂= constant Now, in order to find C2let’s denote

√

V0= γ∗and let t = 0

|γ∗− (α∗+ β |α∗+β

|γ∗_{− (α}∗_{− β )|}α∗−β = C2.

Now, we can express our implicit form of solution in the following way,     √ V_t− (α∗+ β √ V₀− (α∗_{+ β )}     α∗+β ·     √ V₀− (α∗− β √ V_t− (α∗_{− β )}     α∗−β = exp(−κtβ ) (3.22) So, let’s consider the second equation on (3.20) and at least try to find implicit solution for S_t and we take derivative w.r.t time

S0= µSt− 1 2VtSt+ 4 √ V_tS_tθi, j, S0= St µ − 1 2Vt+ 4 √ Vtθi, j
, 1 S_tdS= µ − 1 2Vt+ 4 √ V_tθi, j
dt, ln |St| +C3= µt − Z ₁ 2Vt+ Z 4√V_tθi, j, StC1= exp(µt − Z ₁ 2Vt+ Z 4√Vtθi, j), S_t= exp(µt − R ₁ 2Vt+ R 4√V_tθi, j) C1 , where C1> 0,C3= constant (3.23)

Now, let’s return to (1.2)

C= e−rT

n

∑

i=1

λiX(ωi)

where, we will define

X(ω) = max{S(T, ω) − K, 0}. (3.24) Stage 1: Now let’s demonstrate our scheme of calculations. The first we solve i = 1, path (t, 0, 0) with weight 1/2. Our equation on (3.20) will get form:

 S(t) V(t)  = S(0) V(0)  + Z t 0  µ S(s) −1₂S(s)V (s) κ (θ − V (s)) −1₄σ2  ◦ ds if we take derivative w.r.t time we get

 S0(t) V0(t)  =  µ S(t) −1₂S(t)V (t) κ (θ − V (t)) −1₄σ2. 

It could easily be shown that with initial value S(0) and V (0) we would attain some con-stant by solving these two equations, after we find directly our X (ω1) from (3.24).

Stage 2: Now, we move to second stage and calculate i = 2, path (t, ω₁[2](t), ω₂[2](t)) and with weight 1/24. By doing the same procedure we get explicit equations that we showed above (3.23) and (3.22). However, here we divide our path for (0, 1/4), (1/4, 1/2), (1/2, 3/4), (3/4, 1) intervals t ∈ [0, 1/4] : for t = 0 with initial values S(0) and V (0) and 4θ1,1, 4θ1,2 we

solve (3.23) and (3.22) and get S0and V0, for t = 1/4 we calculate again (3.23) and (3.22) and

get values S1and V1. As displayed, we continue using the same manipulation until we get the

value of S4and V4that we are interested in, and after we put our S4and get the value of X (ω2)

in (3.24).We continue in the same matter until we find the values in j = 13 path. The final step in finding our contingent claim would be putting all 13 values of X (ωj) into (1.2).

We will demonstrate the calculation only for path ω₁[2](t), because the path√3ω₂[2](t) will be calculated at same way, see Figure 3.1, where β0= 4θ1,1, β1= 4θ1,2, β2= 4θ1,3, β3= 4θ1,4,

Figure 3.1: ω₁[2](t) path

values of θi, j could be found on the Table A.2 and paths and weights on Table A.2. As you

see as we mentioned above we are interested on the final value S4.

Theorem 3 states existence of constant that has small error for T < 1 ,as order of our Cubature method is m = 5 for small time interval our method is more preferable than usual Monte Carlo method that has standard error of _√σ

n, here n number of simulation and σ

-standard deviation. Our justification could be described as follow: in order to decrease error two time in Monte Carlo method one should do simulation 4 times, however if you divide our time in Theorem 3 by two as result you get 8 times less error compare to Monte Carlo method. Our justification for best method and which type is better could be answered by simple answer; the best method is that method which gives less error and requires less resources to compute contingent claim in our case it is Explicit cubature method in Heston model. So, we conclude that we at least theoretically calculated contingent claim, however, if out time interval is large we can not use this method. Fortunately, there is other methods for T > 1 ,where pricing of european option in a market model with two stochastic volatilities described, see:Canhanga et al. (2017)

3.3

Conclusion

At the beginning of our journey, we thought we had enough theoretical background in order to handle our problem. However, as we proceed we got introduced to the interesting and fascin-ating world of mathematical problems that we never before had encountered. Our task was to demonstrate Explicit cubature method in Heston model in financial engineering, that is much more quicker than usual Monte Carlo method commonly used in the industry. Usual Monte Carlo method handles random variables and even with modern technology it takes time to get results with small error because of the need of increasing number of simulations. Our method relates to the discrete case where our desired information lies on piecewise linear continuous path that is well defined without any randomness. Unfortunately, it was a very tedious task to calculate all ODEs which will be around 50 by hand. We could only demonstrate the general case of calculating scheme up to i = 2. Another important aspect of this paper is to

demon-strate that we could by hand solve contingent claim for small time intervals. Furthermore, it follows that if we were to write our equations with algorithm using any scientific tool such as MATLAB, the time for reaching results would be less than a seconds. Another important remark here is that our accuracy of calculating contingent claim is exact up to degree 5 of iterated Stratonovich integrals, after that order, it will generate some error.

Finally, the answer for question one would be Theorem 3 and for question two we can answer that If we have two formulae of type (1.2), under condition that both have equal m ,the best approximation would be that one which has less N.

Appendix A

Criteria of Bachelor Thesis

The Swedish National Agency for Higher Education emphases certain requirements for bach-elor theses in mathematics, mathematical statistics, financial mathematics and actuarial sci-ence to be fulfilled. Obtained from the report “How to Write a Thesis” written by Sergei Silvestrov, Anatoliy Malyarenko and Dmitrii Silvestrov these requirements will here be con-sidered and explained how they where satisfied.

In this thesis, various areas in mathematics have been applied, investigated and adapted. The main areas have been probability theory, statistics, stochastic calculus, linear algebra and functional analysis. The major field of study has been the Heston model by cubature formula. To form a basis of knowledge to the subject, the authors chose to review the fundament-als of financial engineering mainly by refering to Carkovs et al. (2011) and Kijima (2013). The introduction is written in a inductive way meaning for it to be a fundamental and solid basis for the remaining theory of the thesis. All definitions and method used are presented and explained in comprehensive manner for the readers to easily follow and understand the subject. The theory and information found in this thesis has been taken from books, articles and scientific papers distributed by Professor Anatoliy Malyarenko and also collected by the authors themselves using licenese from the library at Mälardalens University Västerås. Chapter 2 summarize the essence of the formulation of the problem and introduce the readers to the major theory parts of the thesis. In Theorem 3 the conclusion is drawn that if our time interval is large than the mentioned method can not be applied, the authors then encourage and propose other methods by giving the references to other research. Overall this is a high level mathematical thesis where the authors have shown their ability to investigate and critically discuss information, problems and solutions.

An oral presentation of the thesis will be on June 2nd 2017. The authors will then clarify and explain the methods and results of the research and allow the attending audience to ask questions regarding the outcomes of the research.

Paths and weights for d = 2, degree 5 paths weights (t, 0, 0) 1/2 (t, ω₁[2](t),√3ω₂[2](t)) 1/24 (t, −ω₁[2](t),√3ω₂[2](t)) 1/24 (t, ω₁[2](t), −√3ω₂[2](t)) 1/24 (t, −ω₁[2](t), −√3ω₂[2](t)) 1/24 (t, ω₂[2](t),√3ω₁[2](t)) 1/24 (t, −ω₂[2](t),√3ω₁[2](t)) 1/24 (t, ω₂[2](t), −√3ω₁[2](t)) 1/24 (t, −ω₂[2](t), −√3ω₁[2](t)) 1/24 (t, 2ω₁[1](t), 0) 1/24 (t, −2ω₁[1](t), 0) 1/24 (t, 0, 2ω₁[1](t)) 1/24 (t, 0, −2ω₁[1](t)) 1/24

Table A.1: Lyons and Victoir (2004)

Numerical values for d = 2, L = 4

slope values θ1,1 0.048077867969717 θ1,2 0.089335677121310 θ1,3 1.177095041847874 θ1,4 -0.314508586938901 θ2,1 1.012368507495347 θ2,2 -1.61960476379500 θ2,3 1.702104005103968 θ2,4 -0.094867748804311

Bibliography

B. Canhanga, A. Malyarenko, Y. Ni, and S. Silvestrov. Fast Monte Carlo pricing of European options in a market model with two stochastic volatilities. Working paper, 2017.

J. Carkovs, A. Malyarenko, and K. Pärna. Exploring the world of financial engineering. Mälardalen University, 2011.

M. Kijima. Stochastic processes with applications to finance. CRC Press, 2013.

P. E. Kloeden and E. Platen. Stratonovich and Itô stochastic Taylor expansions. Math. Nachr., 151:33–50, 1991. ISSN 0025-584X.

P. E. Kloeden and E. Platen. Numerical solution of stochastic differential equations, volume 23 of Applications of Mathematics (New York). Springer-Verlag, Berlin, 1992.

S. Kusuoka. Approximation of expectation of diffusion process and mathematical finance. In Taniguchi Conference on Mathematics Nara ’98, volume 31 of Adv. Stud. Pure Math., pages 147–165. Math. Soc. Japan, Tokyo, 2001.

T. Lyons and N. Victoir. Cubature on Wiener space. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 460(2041):169–198, 2004.

E. Platen and W. Wagner. On a Taylor formula for a class of Itô processes. Probab. Math. Statist., 3(1):37–51 (1983), 1982. ISSN 0208-4147.

W. Rudin. Principles of mathematical analysis. McGraw-Hill Book Co., New York-Auckland-Düsseldorf, third edition, 1976.

W. Rudin. Functional analysis. International Series in Pure and Applied Mathematics. McGraw-Hill, Inc., New York, second edition, 1991.

A. H. Stroud. Approximate calculation of multiple integrals. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1971.

V. Tchakaloff. Formules de cubatures mécaniques à coefficients non négatifs. Bull. Sci. Math. (2), 81:123–134, 1957. ISSN 0007-4497.