Variance reduced Monte-Carlo Simulations of Stochastic Differential Equations

(1)

deposit_hagen

Variance reduced

Monte-Carlo Simulations of Stochastic Differential Equations

Wirtschafts- wissenschaft

Dissertation

Armin Heinrich Müller

(2)

(3)

Inaugural-Dissertation zur Erlangung des Doktorgrades der Fakult¨at f ¨ur Wirtschaftswissenschaft

an der FernUniversit¨at in Hagen

Variance reduced

Monte-Carlo Simulations of Stochastic Differential Equations

Vorgelegt von Dipl.-Phys. Armin Heinrich M¨uller, M.Sc.

Tag der Disputation: 5. September 2017

(4)

(5)

(6)

(7)

(8)

(9)

I Introduction and Theoretical Background 13

1 Introduction 15

2 Theoretical Background 19

2.1 Stochastic Processes . . . 19

2.2 It¯o’s Lemma . . . 21

2.3 Geometric Brownian Motion . . . 22

2.4 Option Contracts . . . 24

2.5 The Black-Scholes Model for Option Pricing . . . 24

2.6 The Feyman-Kac Formula . . . 27

2.7 Risk-neutral Asset Pricing . . . 30

2.8 Monte-Carlo Simulations . . . 30

2.9 Variance Reduction Techniques . . . 31

II Results and Discussion 39

3 An Application of the Put-Call-Parity to Variance reduced Monte- Carlo Option Pricing 41 3.0 Abstract . . . 41

3.1 Introduction . . . 41

3.2 Monte-Carlo Simulation of the Feynman-Kac Formula . . . 42

3.3 Put-Call-Parities for several Options . . . 46

(10)

Contents

3.4 Numerical Results . . . 49

3.5 Discussion . . . 54

3.6 Conclusion . . . 56

3.7 Acknowledgments . . . 56

4 A joint Application of the Put-Call-Parity and Importance Sam- pling to Variance reduced Option Pricing 57 4.0 Abstract . . . 57

4.2 Importance Sampling . . . 58

4.3 Variance Reduction with Put-Call-Parities for European and Asian Options . . . 62

4.4 Joint Application of the Put-Call-Parity and Importance Sampling . 63 4.5 Numerical Results . . . 63

4.6 Discussion . . . 69

5 Improved Variance reduced Monte-Carlo Simulation of in-the-Money Options 71 5.0 Abstract . . . 71

5.3 Variance Reduction with Put-Call-Parities for European and arithmetic Asian Options . . . 75

5.4 Joint Application of the Put-Call-Parity and Importance Sampling . 76 5.5 Numerical Results . . . 77

6 Approximations of Option Price Elasticities for Importance Sam- pling 81 6.0 Abstract . . . 81

6.3 Approximations of Option Price Elasticities . . . 91

6.3.1 General Remarks . . . 91

6.3.2 Approximation by constant Values . . . 93

8

(11)

Contents

6.3.3 Approximation by Step Function . . . 94

6.3.4 Approximation by lower Bound . . . 96

6.3.5 Approximation by Black-Scholes Formula . . . 101

6.3.6 Approximation by Regression . . . 104

6.3.7 Approximation by Integration . . . 109

7 Variance reduced Value at Risk Monte-Carlo Simulations 117 7.0 Abstract . . . 117

7.2 The Risk Measure VaR . . . 118

7.3 Importance Sampling of VaR Estimators . . . 118

7.3.1 Basics of Importance Sampling . . . 119

7.3.2 A Variance reduced VaR Monte-Carlo Estimator . . . 120

7.3.3 Elasticity Approximations for Variance reduced VaR Esti- mates . . . 121

7.4 Numerical Results . . . 125

7.5 Discussion of Results . . . 126

III Conclusion and Annex 129

8 Conclusion 131

List of Figures 134

List of Tables 137

List of Abbreviations 139

List of Symbols 141

Bibliography 145

Curriculum Vitae 153

Acknowledgments 155

(12)

(13)

Foreword

The Ph.D. dissertation of Armin M¨uller, entitled Variance reduced Monte-Carlo Simulations of Stochastic Differential Equations, treats the topic of stochastic rep- resentations of solutions of deterministic partial differential equations, in this case the so-called Black-Scholes equation. It first appeared in the famous paper of Black and Scholes (1973) [1], derived from an arbitrage argument. A clever riskless portfolio combination of stock and option in variable proportions leads to a partial differential equation and the random elements in the stock price movement cancel.

Nevertheless, following Cox and Ross (1976) [2], the solution can be represented by a discounted expectation value of the terminal price, averaged over pseudo price movements, where the drift is the riskless rate.

This is the starting point of the thesis, namely the estimation of the expectation value by statistical simulation. A naive mean value estimate leads to large standard errors, only controlled by large sample size. The well known idea of importance sampling is applied in this context, namely to simulate distorted trajectories which yield a small variance of the weighted average of the payoff function. In order to simulate from the optimal probability measure, one can derive a drift correction to the Itˆo equation for the underlying, which involves both the diffusion matrix and the unknown solution.

In the thesis, various approximations of the true drift correction are discussed, ranging from the inner value to Black-Scholes type fit functions and integration based approximations. In many cases, strong variance reduction can be obtained, often in combination with the put-call parity. Furthermore, other types of functionals, such as the value at risk, are discussed and computed with small error.

(14)

Foreword

I hope, that the text finds an interested audience and will improve the pricing of options and other functionals in economics and physics.

Prof. Dr. Hermann Singer

Chair of Applied Statistics and Empirical Social Research FernUniversit¨at in Hagen

Germany

12

(15)

Part I

Introduction and Theoretical

Background

(16)

(17)

1

Introduction

In 1973, Fischer Black and Myron Scholes published their famous work “The Pric- ing of Options and Corporate Liabilities” [1]. They introduced a deterministic partial differential equation – the so-called Black-Scholes differential equation – which the price of any financial derivative must fulfill in order to avoid arbitrage oppor- tunities. For the special case of a European option on an underlying whose price is supposed to follow a geometric Brownian motion, Black and Scholes gave an explicit solution of their differential equation. Thus, an analytic price formula exists. Subsequently, price formulae for several other options and price models were introduced [3].

However, there is a broad range of financial derivatives for which no analytic formula exists – even in comparably simple underlying price models as the geometric Brownian motion. Furthermore, there are complex price models, e.g., models with stochastic volatility, were the existence of closed-form solutions depends on the particular model [3, 4].

In 1976, John C. Cox and Stephen Ross introduced an alternative approach to determine the prices of financial derivatives. They showed that under certain conditions, option prices can be calculated as the expected value of the discounted pay-off function conditioned on the current underlying price [2].

The no-arbitrage argument employed by Black and Scholes yields a deterministic partial differential equation also for more complicated derivatives or underlying price models. This equation is equivalent to the Cauchy problem, named after the French mathematician Augustin-Louis Cauchy. Ultimately, this partial differential equation is an inhomogeneous Kolmogorov backward equation [5].

The Cauchy problem has a unique solution with a stochastic representation called the Feynman-Kac formula, named after the American physicist Richard Feynman [6, 7] and the Polish-American mathematician Mark Kac [8, 9]. The formula is the expected value of the problem’s boundary condition and further terms conditioned

(18)

1 Introduction

on the current value of the stochastic process underlying the problem [5]. It allows to represent a broad range of financial derivatives as the conditional expected value of the discounted pay-off function [10]. For economists, this introduces a more intuitive way of thinking about a derivative price than the sometimes cumbersomely interpretable analytic formulae.

In many cases, it is not possible to explicitly calculate this conditional expected value. A possible approach is to write the expected value as the integral of the expected value’s argument times the transition density of the underlying stochastic process integrated over all possible final values. This integral can be approximated numerically by a Monte-Carlo simulation. Sample paths simulated from the transition density kernel yield a distribution of possible outcomes of the pay-off function.

The discounted arithmetic average of this distribution is an unbiased estimator of the option price [11].

A fundamental property of a Monte-Carlo simulation is that it yields no unam- biguous solution but rather a probability distribution of the quantity of interest. The standard error of a Monte-Carlo estimator depends on the sample size N of the simulation. Typically, it decreases at a rate proportional to1/√

N. Consequently, to decrease the standard error by one order of magnitude, the sample size must increase by a factor of 100 [11].

Ceteris paribus, this means that also the required computing time increases by a factor of 100. In practical applications, e.g., the pricing of derivatives by financial market participants, long computing times are not acceptable [12]. One approach to cope with this problem is to address its technical side by increasing the computing speed. Also, as a Monte-Carlo simulation consists of many realizations of i.i.d. random variables, in many cases, it is suitable for a parallelization of computations.

Another approach is to address the variance of Monte-Carlo estimators mathematically. In many cases, more than one unbiased estimator of a quantity of interest exists. Then, it is possible to determine (unbiased) estimators with lower variance than the initial (“naive”) estimator. A broad range of techniques exists to come up with so-called variance reduced Monte-Carlo estimators [11].

The focus of this work lies on the analysis of such estimators, treating two different approaches. The first – mathematically rather simple – approach is the application of put-call-parities that allow to represent a random variable as the sum of deterministic elements and another random variable. Under certain conditions, it is possible to employ this property for variance reduction. The second – mathematically more challenging – approach is called importance sampling. The basic idea is to generate sample paths using an alternative measure of probability. As a result,

16

(19)

areas of interest are sampled more accurately.

An intuitive example, where importance sampling is appropriate, is the Monte- Carlo simulation of a deep-out-of-the-money European call. In a conventional Monte-Carlo simulation, many sample paths end below the strike price. Yielding a pay-off function value of 0, they do not contribute to the Monte-Carlo estimator. However, they do increase its variance. Contrastingly, a Monte-Carlo simulation enhanced by importance sampling artificially pushes sample paths above the strike price. This approach can decrease the dispersion of the constituents of the Monte-Carlo estimator, thus reducing its variance. Of course, further mathematical requirements must be fulfilled in order to obtain an unbiased estimator [13].

The presented analysis focuses on the examination of a formula introduced by Singer in 2014 [14]. This theoretically optimal formula for importance sampling entails the difficulty that knowledge of the elasticities of the quantity of interest with respect to the stochastic elements underlying the model is required. Like the quantity itself, its elasticities are generally unknown. To apply the formula, it is necessary to conduct approximations of the elasticities. This approach is applicable to option pricing, but also to the Monte-Carlo simulation of other functionals [14, 15, 16, 17].

The work is structured as follows: Chapter 2 lays the theoretical foundation of this doctoral thesis. It introduces basics of stochastic processes and provides a derivation of It¯o’s lemma. The geometric Brownian motion is introduced as commonly used stochastic process to model stock prices. A description of the Black- Scholes model follows an overview of basic properties of options on financial assets. The Feynman-Kac formula is introduced as the unique solution of the Cauchy problem. It yields a stochastic representation of prices of many different option types. Consequently, the concept of risk-neutral asset pricing is motivated. An introduction to Monte-Carlo simulations is followed by a broad overview on different variance reduction techniques.

The article reprinted in chapter 3 derives put-call-parities for several option types.

These put-call-parities are then applied to conduct variance reduced Monte-Carlo simulations of the considered options. The article was originally published as follows:

[18] Müller, Armin. An application of the put-call-parity to variance reduced Monte-Carlo option pricing. FernUniversität in Hagen – Fakultät für Wirtschaftswissenschaft – Discussion paper No. 495;

(2016).

Chapter 4 contains an article that analyzes the joint application of the put-call-

(20)

1 Introduction

parity approach already discussed in the previous chapter and the importance sampling approach introduced by Singer [13, 14]. It introduces the basic concepts of the applied importance sampling approach. The joint approach generates synergies in the sense that the variance reduction factors achieved surpass the products of the variance reduction factors achieved by the standalone approaches. Both European and arithmetic Asian options are analyzed. Originally, the article was published as follows:

[19] Müller, Armin. A joint application of the put-call-parity and importance sampling to variance reduced option pricing. FernUniver- sität in Hagen – Fakultät für Wirtschaftswissenschaft – Discussion paper No. 496; (2016).

The next chapter contains a summarized version of the two previous chapters. It contains further numerical examples that allow to compare our approach with other approaches discussed in literature. It was originally published as follows:

[20] M¨uller, Armin. Improved variance reduced Monte-Carlo simulation of in-the-money options. Journal of Mathematical Finance;

6(3):361–367; (2016).

The article included in chapter 6 gives a more profound derivation of the importance sampling approach analyzed in this thesis than the two previous chapters. It introduces several approaches suitable to approximate option price elasticities. Nu- merical examples illustrate the different approximations provided. The article was originally published as follows:

[21] M¨uller, Armin. Approximations of option price elasticities for im-

portance sampling. FernUniversität in Hagen – Fakultät für Wirtschaftswis- senschaft – Discussion paper No. 499; (2016).

The last article of this thesis – reprinted in chapter 7 – analyzes the application of Singer’s importance sampling formula to a functional other than options: It demon- strates how the approach can be applied to Monte-Carlo simulations of the risk measure Value at Risk. Originally, the article was published as follows:

[22] Müller, Armin. Variance reduced Value at Risk Monte-Carlo simulations. FernUniversität in Hagen – Fakultät für Wirtschaftswis- senschaft – Discussion paper No. 500; (2016).

The thesis concludes in chapter 8 with a short summary and references to future research areas.

18

(21)

2

Theoretical Background

2.1 Stochastic Processes

Brownian Motion In 1827, the British botanist Robert Brown discovered the stochastic movements of microscopically small particles in liquid solutions. First, he interpreted these dynamics as the motion of small micro-organisms. Later, he realized that he actually observed the random stochastic movements of inanimate particles. This process calledBrownian motion can be observed by analyzing the movement of large molecules under the microscope [23].

Several decades later, in 1900, the French mathematician Louis Bachelier first applied the concept of Brownian motions to describe the statistic movement of security prices on stock markets [24, 25]. Bachelier’s work did not yet fully describe the dynamics of financial markets in a satisfactory manner. However, his work can be considered to be the starting point of a century full of advances in the theory of mathematical finance [26]. Around the same time, in 1905, nobel laureate Albert Einstein applied the equations of Brownian motion to the heat flow in thermody- namic systems [27].

Markov Property An important property of stochastic processes that is assumed to hold true in many models describing the dynamics of financial markets over time is the Markov property, named after the Russian mathematician Andrei Andreje- witsch Markov. It states that only the current value of a stochastic process X(t) influences its future development, while the past of its trajectory has no explicit

(22)

2 Theoretical Background

influence¹. This property can be expressed as follows:

p[X(t+s)|X(τ),0≤τ ≤s] =p[X(t+s)|X(s)] (2.1) The Markov property reflects the weak form of the efficient market hypothesis, which assumes that the current price of a security reflects all the available informa- tion from the past [3, 5].

Wiener Processes The Wiener processW(t), named after the American mathematician Norbert Wiener [28], is an example for a process that shows the Markov property. It has the following properties:

1. The increment∆W(t)over an incremental time interval∆tis given by

∆W(t) :=W(t+ ∆t)−W(t) =ξ_i√

∆t (2.2)

whereξ_iis a standard normally distributed random variable.

2. For any two given intervals (∆t)_i and(∆t)_j, the corresponding increments (∆W(t))_i and(∆W(t))_j are stochastically independent.

The changeW(T)−W(0) of the Wiener processW(t)over a longer period of timeT can be calculated as the sum of the incremental changes. By setting

n= T

∆t, (2.3)

the equation

W(T)−W(0) =

n

X

i=1

ξ_i√

∆t (2.4)

results, where again ξi,(i = 1,2, . . . , n) is a sequence of independent standard normally distributed random variables. From the linearity of expected values, we obtain

IE[W(T)−W(0)] =

n

X

i=1

IE[ξ_i]√

∆t= 0. (2.5)

According to the second property of a Wiener process, increments are stochasti-

1The time argument of stochastic processes is written in parenthesis in this work. Thus,X(t)is used instead of Xtto avoid confusion with the abbreviation used for partial derivatives, e.g., Ct≡∂C/∂t.

20

(23)

2.2 It¯o’s Lemma

cally independent. Therefore, the Bienaym´e formula can be applied and Var [W(T)−W(0)] =

n

X

i=1

Varh ξ_i√

∆ti

=n∆t=T (2.6)

results [3].

Typically, many of the asset pricing models introduced in the second half of the 20^th century use a continuous-time setting instead of a discrete-time setting. By taking the limit∆t → 0, the continuous-time standard Brownian motion with increments dW(t)results [10].

It ¯o Processes An example for a more general stochastic process is the It¯o process, named after the Japanese mathematician It¯o Kiyoshi [29, 30]. Its stochastic properties strongly depend on the properties of the underlying Wiener process. It can be described by the stochastic differential equation

dX(t) =f[X(t)]dt+g[X(t)]dW(t). (2.7) In this equation, f[X(t)] is the drift vector, g[X(t)] is the diffusion coefficient withgg^T = Ω, where Ωis the diffusion matrix, andW(t)is a multivariate Wiener process [3, 13].

2.2 It ¯ o’s Lemma

In section 2.4, we introduce options on underlying financial instruments. Mathe- matically, an option is a non-linear function of a stochastic function. In this section, we introduce a general stochastic differential equation for any functionY [X(t), t]

of a stochastic processX(t) and timet. This result is called It¯o’s lemma or It¯o’s formula.

We consider the functionY =Y [X(t), t]. With∆X(t) :=X(t+ ∆t)−X(t), a Taylor series expansion up to the second order yields the increment

∆Y(t) :=Y [X(t) + ∆X(t), t+ ∆t]−Y [X(t), t]

= ∂Y

∂X(t)[X(t), t] ∆X(t) + ∂Y

∂t [X(t), t] ∆t +1

2

∂²Y

∂X(t)² [X(t), t] ∆X(t)²+ ∂²Y

∂X(t)∂t[X(t), t] ∆X(t)∆t +1

2

∂²Y

∂t² [X(t), t] ∆t²+. . . .

(2.8)

(24)

Discretizing equation (2.7) and using equation (2.2) yields²

∆X(t) =f∆t+gξ_i√

∆t. (2.9)

By squaring this equation, one obtains

∆X(t)² =g²ξ_i²∆t+O

∆t³²

. (2.10)

ξ_i is a random variable that follows a standard normal distribution. Therefore, IE(ξ_i) = 0andVar(ξ_i) = 1result by definition. From

Var(ξ_i) :=IE(ξ_i²)−IE(ξ_i)² = 1 (2.11) it follows that IE(ξ_i²) = 1 and therefore IE(ξ_i²∆t) = ∆t. It can be shown that Var(ξ_i²∆t)is of order∆t². Therefore, for small∆t, equation (2.10) can be considered deterministic, i.e.

∆X(t)² ≈g²∆t. (2.12)

By taking the limits∆t → 0and∆X → 0, using equation (2.12) and ignoring terms of order∆t³² or higher, equation (2.8) simplifies to

dY(t) = ∂Y

∂X(t)[X(t), t]dX(t) + ∂Y

∂t [X(t), t]dt+1 2

∂²Y

∂X(t)² [X(t), t]g²dt

= ∂Y

∂X(t)[X(t), t]f +∂Y

∂t [X(t), t] + 1 2

∂²Y

∂X(t)² [X(t), t]g²

dt + ∂Y

∂X(t)[X(t), t]gdW(t)

= ∂Y

∂Xf +∂Y

∂t + 1 2

∂²Y

∂X²g²

dt+ ∂Y

∂XgdW.

(2.13) In the second line, definition (2.7) has been inserted. In the last line, arguments have been suppressed for simplicity [3, 16, 31].

2.3 Geometric Brownian Motion

As mentioned in section 2.1, the Brownian motion does not satisfactorily describe the price movements of financial products. One major shortcoming is that a stochastic process following a Brownian motion can take on negative values. However, it

2For simplicity, we suppress arguments of the coefficientsf[X(t), t]andg[X(t), t].

22

(25)

2.3 Geometric Brownian Motion would not be appropriate to allow for negative values when modeling the price of a stock. A second disadvantage is that there is no empirical evidence for normally distributed stock prices [32].

A stochastic process commonly used to model stock prices is the geometric Brownian motion. For instance, the Black-Scholes model employs this stochastic process [1]. The stock price S is assumed to follow the stochastic differential equation

dS =µSdt+σSdW (2.14)

where µ is the expected rate of return on the stock, σ its volatility and W again a Wiener process. This equation corresponds to equation (2.7), choosingX = S, f =µS andg =σS.

To obtain a solution forS, due to the diffusion term, it is not possible to integrate equation (2.14) directly. A possible approach is to set

Y = lnS (2.15)

and to calculate dY from It¯o’s formula (2.13). Because of

∂Y

∂S = 1

S, ∂²Y

∂S² =− 1

S² and ∂Y

∂t = 0, (2.16)

dY can be expressed as follows:

dY =

µ− σ² 2

dt+σdW (2.17)

Integration over the interval(t, T)and substitution yields S(T) =S(t) exp

µ− σ²

2

(T −t) +σ[W(T)−W(t)]

. (2.18)

Equation (2.17) implies that if S follows a geometric Brownian motion, its log- return is normally distributed. The additional term −^σ₂² is usually called Wong- Zakai correction[3, 13].

A valid criticism with regards to the application of the geometric Brownian motion in asset pricing is the empirical evidence for heavy-tailed return distributions.

In the case of heavy-tailed returns, extreme events are not appropriately reflected by the geometric Brownian motion. For financial risk management, this is particu- larly challenging and dangerous in the case of extreme losses. However, due to its simplicity, the model of geometric Brownian motion is frequently used [32].

(26)

2.4 Option Contracts

Over the last decades, financial derivatives like future and option contracts have become increasingly important. In 1973, the Chicago Board Options Exchange was founded and became the first exchange market for trading options [33]. At about the same time, Fischer Black and Myron Scholes published their famous work “The Pricing of Options and Corporate Liabilities” [1]. Several works from other authors like Robert Merton, John C. Cox, Stephen Ross and Mark Rubinstein followed, giving a theoretical background to the pricing of derivative assets [2, 34, 35, 36]. A noble price in economics was granted to Robert Merton and Myron Scholes in 1997 [37].

Since 1973, a rich landscape of different option types with a broad range of underlying assets developed. In this introduction, the scope is limited to the description of European call options. A comprehensive overview of different derivative types can be found in [3].

A European call options grants the right to buy an underlying asset at a predefined strike priceK at the options date of maturity T. The option buyer (the long side) pays a premium (the option price) to the option seller (the short side) in order to receive this right. In contrast to other derivative products like futures or forwards, the option buyer adopts no binding position to buy the underlying. He has the right to decide whether or not to exercise the option. The option seller must sell the underlying at the price K if the option buyer decides to exercise the option.

Obviously, the option buyer only exercises the option if the underlying’s price at maturityS(T)exceeds the strike priceK. In this case, he realizes a gain ofS(T)− K. Otherwise, he allows the option to expire without exercising it. In this case, he can buy the underlying asset at a cheaper price on an exchange market. From the option buyer’s perspective, the option’s value at maturity can be defined as follows:

C=C[S(T), T] = max [S(T)−K,0] =: [S(T)−K]⁺ (2.19) The determination of the option price at earlier instancestwitht < T is not trivially obvious and will be the content of the next section [3].

2.5 The Black-Scholes Model for Option Pricing

In 1973, Black and Scholes presented an option pricing model that allows to determine prices of European options based on a no-arbitrage argument and assuming a geometric Brownian motion for the underlying price S(t). Following [13], the

24

(27)

2.5 The Black-Scholes Model for Option Pricing Black-Scholes differential equation can be derived by considering thehedge portfolio

V =C+τ S. (2.20)

It contains one European call with priceC andτ units of the underlying stock with priceS. The option price is considered to be a function of the underlying and time, i.e.C ≡C(S(t), t). The total differential dV result from a Taylor series expansion or from directly applying It¯o’s formula (2.13)³:

dV =dC+τdS

=CSdS+Ctdt+1

2CSSdS²+1

2Cttdt²+CStdSdt+τdS

= (C_S+τ)dS+

C_t+1

2C_SS(σS)²

dt+O dt^3/2

(2.21)

Note that in the first line, no elementSdτappears. This results from the assumption thatτis kept constant over a short interval. It is then adjusted to offset the increment ofS in this interval [1]. In the second line, equations (2.2) and (2.14) were used, suppressing terms of the order dt³² or higher.

In equation (2.21), the first term in parenthesis can be forced to vanish by choosing τ = −C_S. The increments of the hedge portfolio V then do not explicitly depend on changes of the underlying priceSany more. Thus,V becomes riskless.

Consequently, the return on the portfolioV must equal the risk-free rater:

dV =

C_t+ 1

2C_SS(σS)²

dt

=! rVdt

= (rC−rC_SS)dt

(2.22)

The first row represents a repetition of equation (2.21) after setting τ = −C_S. In the third row, the definition ofV from equation (2.20) is inserted. Rearranging elements yields the Black-Scholes differential equation

C_t+rSC_S+1

2σ²S²C_SS −rC = 0. (2.23) For European calls with boundary condition C[S(T), T] = [S(T)−K]⁺, Black and Scholes presented a closed form solution [1]:

C(S(0),0) =S(0)Φ (d₁)−Ke^−rTΦ (d₂) (2.24)

3Going forward, the abbreviation^∂C_∂S ≡CSis frequently used.

(28)

with

d_1,2 =

ln^S(0)_K +

r± ^σ₂² T σ√

T .

Black and Scholes found this solution by substituting

C[S(t), t] =e^−r(T^−t)·u(x, y) (2.25) with

x= 2 σ²

r−σ²

2

d₂ (2.26)

and

y = 2 σ²

r− σ²

2 2

(T −t). (2.27)

By using this substitution, the Black-Scholes partial differential equation (2.23) transforms into the heat-transfer equation

∂u

∂y = ∂²u

∂x² (2.28)

with boundary condition

u(x,0) =







0, x <0

K

exp

x ^σ

2 2

r−^σ₂²

−1

, x≥0. (2.29)

The heat-transfer equation has a well-known solution [38] which after resubsti- tution yields the Black-Scholes formula presented in equation (2.24).

A less instructive, but mathematically less challenging way to verify the accuracy of the Black-Scholes formula (2.24) is to calculate its partial derivatives (the so-called “Greeks”) and insert them into the Black-Scholes equation (2.23). One obtains [39]

∆ :=C_S = Φ(d₁), (2.30)

Θ :=C_t =− σ 2√

T −tSφ(d₁)−rKe^−r(T^−t)Φ(d₂), (2.31) and

Γ :=C_SS = 1 Sσ√

T −tφ(d₁) (2.32)

which together with the Black-Scholes formula from equation (2.24) obviously solve the Black-Scholes differential equation (2.23).

26

(29)

2.6 The Feyman-Kac Formula

An alternative approach to solve the Black-Scholes partial differential equation under given boundary conditions is to express the option price C as a probabilistic term. More specifically,C can be expressed as a conditioned expected value of the discounted boundary condition. First, we will introduce the Cauchy problem and the Feynman-Kac formula as the stochastic representation of its solution. Second, we will specify this solution to the special case of a European call in the Black- Scholes model. Third, we will consider the Feynman-Kac formula of an exotic option.

The Cauchy Problem We consider ad-dimensional multivariate stochastic pro- cessX(t)that follows the It¯o stochastic differential equation (2.7). For t ≥ 0we introduce the Kolmogorov backward operatorLby defining

(Lv) (X, t) :=

d

X

i=1

f_i(X, t)∂v(X, t)

∂X_i + 1 2

d

X

i=1 d

X

k=1

Ω_ik(X, t)∂²v(X, t)

∂X_i∂X_k (2.33) wherev is a twice continuously differentiable function,f_i are the elements of the drift vector and Ω_ik are the elements of the diffusion matrix Ω = gg^T. The drift vector f and the diffusion coefficient g correspond to those from equation (2.7).

Note that the definition ofLdoes not include partial derivatives with respect to time t.

Under further technical assumptions⁴, the Cauchy problem

∂v

∂t +Lv−kv+j = 0 (2.34)

with potentialk, Lagrangianj and with boundary condition

v[X(T), T] =h[X(T)] (2.35) has a unique solutionv(X, t)with the stochastic representation

v(X(t), t) =IE

h[X(T)] exp

− Z T

t

k[X(τ), τ]dτ

+ Z T

t

j[X(s), s] exp

− Z s

t

k[X(τ), τ]dτ

ds

X(t) = X_t

. (2.36)

4See [5], p. 366 for details

(30)

The mathematically challenging proof can be found at [5]. A less technical, less rigorous verification of the result is given by [10].

The Feynman-Kac Representation of the Black-Scholes Formula In the Black-Scholes model, for the underlying’s stochastic processS(t), we choose the geometric Brownian motion

dS =rSdt+σSdW. (2.37)

Note that for reasons that will be discussed in section 2.7, we use the risk-free rate rinstead of the underlying’s returnµ.

Equation (2.34) simplifies to a1-dimensional special case with X(t)≡S(t),

v ≡C[S(t), t], k ≡r,

j ≡0and Lv ≡rSC_S+1

2σ²S²C_SS.

(2.38)

Thus, one recovers the Black-Scholes equation (2.23). For a European call with boundary condition

h[X(T)] =C[S(T), T] = [S(T)−K]⁺, (2.39) one obtains the Feynman-Kac representation of the Black-Scholes formula

C[S(t), t] =e^−r(T^−t)IE

[S(T)−K]⁺ |S(t) =S_t

. (2.40)

It can be shown explicitly [10, 40], that this representation is equivalent to the solution in equation (2.24) presented by Black and Scholes [1].

The Feynman-Kac Representation of exotic Options Also for more complex, so-called “exotic” options, with more complicated boundary conditions, the stochastic representation of the unique solution of the Cauchy problem is a pow- erful approach. It yields a (deterministic) partial differential equation for an exotic option’s price function and a corresponding solution.

As an example, we consider the case of an arithmetic Asian option. The pay- off function does not only depend on the underlying’s final value S(T). It is a

28

(31)

2.6 The Feyman-Kac Formula function of the arithmetic average value of the price, i.e. the pay-off function is path-dependent.

For the underlying S(t), we choose again the stochastic process described in equation (2.37). Additionally, we introduce a stochastic processY(t)which is the integral of the underlying’s price with respect to time:

dY =Sdt, Y(0) = 0 (2.41) The option considered has the following pay-off function:

h[X(T)] =C[S(T), Y(T), T] =

Y(T) T −K

+

, (2.42)

Choosing the same approach as for the European call, for equation (2.34), we obtain a2-dimensional model with

X₁(t)≡S(t), X₂(t)≡Y(t),

v ≡C[S(t), Y(t), t], k ≡r,

j ≡0and

Lv ≡rSC_S+SC_Y +1

2σ²S²C_SS.

(2.43)

An equation similar to the Black-Scholes equation (2.23) with an additional term SC_Y results⁵:

Ct+rSCS+SCY + 1

2σ²S²CSS−rC = 0. (2.44) Consequently, one obtains the following stochastic representation of the option price⁶:

C[S(t), Y(t), t] =e^−r(T^−t)IE

"

Y(T) T −K

+

S(t) =St, Y(t) =Yt

#

. (2.45)

5The same partial differential equation results from a no-arbitrage argument, for details see chapter 3 or [18].

6Equation (2.41) introduces a bivariate stochastic process[dS,dY]^T. On the one hand, this allows to preserve the notation of the previous paragraph on the stochastic representation of the Black- Scholes formula. On the other hand, this comes at the price of a singular diffusion matrix. For details on the singular multinormal distribution see [41].

(32)

2.7 Risk-neutral Asset Pricing

At first glance, it may appear surprising that the Black-Scholes formula (2.24) does only depend on the risk-free raterand not on the underlying’s returnµ. The representation of option prices as expected values presented in the previous section was introduced by Cox and Ross in 1976 [2, 42]. In their work, they showed that it is useful to assume that investors are risk-neutral. The authors argue that the expected rates on the option and its underlying must coincide with the risk-free rate in a risk-neutral world:

IE[S(T)/S(0)|S(0) =S₀] =e^rT

IE[C[S(T), T]/C(S(0),0)|S(0) =S₀] =C[S(0),0]⁻¹IE[C[S(T), T]|S(0) =S₀]

=e^rT

(2.46) Rearranging terms again yields the Feynman-Kac formula for European call options [13]:

C[S(0),0] =e^−rTIE[C[S(T), T]|S(0) =S₀] (2.47) Therefore, the valuation of derivative prices can follow a risk-neutral approach [3]:

1. The underlying’s return is assumed to correspond to the risk-free rate, i.e.µ= r.

2. The expected value of the options pay-off function is calculated in this risk- free world.

3. The expected pay-off is then discounted with the risk-free rate yielding the derivative price.

2.8 Monte-Carlo Simulations

For several options, e.g., for exotic options with complicated pay-off functions, no analytical solution to the deterministic partial differential equation of the option price exists. In these cases, the stochastic representation of the option price can be approximated by a numerical approach conducting a Monte-Carlo simulation.

To obtain a Monte-Carlo estimator, N sample paths are sampled from the It¯o process governing the model considered. For instance, in the Black-Scholes model, one simulates sample paths from a discretized version of equation (2.14) withµ=r

30

(33)

2.9 Variance Reduction Techniques using i.i.d. random numbers. So, for each of the sample paths, one obtains a pay-off valueCˆ_i. The arithmetic average of these valuesCˆ_i yields an unbiased estimator of the derivatives expected value:

Cˆ= 1 N

N

X

i=1

Cˆ_i. (2.48)

As the sample paths are stochastically independent, the estimator’s variance can be calculated from the Bienaym´e formula:

Varh Cˆi

= 1 N²

N

X

i=1

Varh Cˆ_ii

= Varh

Cˆ_ii

N (2.49)

The estimator

Vardh Cˆi

= 1

N(N −1)

N

X

i=1

Cˆ_i²−NCˆ²

!

(2.50) results [43, 11].

2.9 Variance Reduction Techniques

A fundamental property of Monte-Carlo simulation is that the Monte-Carlo estimators’ empirical variances typically decrease at a rate proporational toN⁻¹ whereN is the sample size. This can be seen from equation (2.49).

When conducting a simulation, practitioners demand a certain level of accuracy for Monte-Carlo estimators. Very low sample sizesN may lead to very high variances and therefore to an unacceptably low level of accuracy. One approach to address this problem is to increase the sample sizeN. To avoid more time-consuming simulations, one has to enhance available computing power by increasing the computing speed or by parallelization of calculations. However, as computing power in many cases is a sparse resource, the question arises, whether there are other approaches to reduce the empirical variance.

There is a broad range of so-called “variance reduction techniques” that allow to reduce the empirical variance of Monte-Carlo estimatorsCˆ for a given sample sizeN. These approaches have in common that they do not reduce variance by increasing the sample size. Rather, they aim to construct alternative unbiased Monte- Carlo estimators Cˆ⁰ with lower variance than the naive estimator from equation (2.48). This section introduces several variance reduction techniques summarizing the overview from [11].

(34)

The general aim is to compute an unbiased Monte-Carlo estimator with lowest variance possible within a given computing time. Of course, the variance reduced Monte-Carlo estimator Cˆ⁰ might demand more computing power per replication than the naive estimatorC.ˆ

As the estimators are unbiased by definition, the i.i.d. replications of both estimators have the same expected value, i.e. IE

hCˆ_ii

= IE hCˆ_i⁰i

= C, but σ² = Varh

Cˆ_ii

>Varh Cˆ_i⁰i

=σ⁰².

We assume that the estimatorCˆrequires a computing time of∆tper replication and the variance reduced estimatorCˆ⁰ requires a computing time of∆t⁰ per replication. Under these conditions, it can be shown that the variance reduced estimator performs more efficiently if the inequality

σ⁰²∆t⁰ < σ²∆t (2.51) holds true. In other words, if the product of the variance and the required computing time per replication of the variance reduced estimator exceeds the corresponding product of the naive estimator, it’s not worth to perform the variance reduction technique. In this case, the additional computational effort required per replication to run the variance reduced Monte-Carlo simulation cannibalizes the gains resulting from the lower variance per replication [44].

We now give an overview over several variance reduction techniques.

Control Variates The control variates approach modifies the naive Monte-Carlo estimator by adding an additional term with an expected value of0. Under certain conditions, this modification leads to variance reduction.

Consider a setting whereYˆ₁, . . . ,Yˆ_N areN outputs of a simulation ofY. The aim is to estimate the expected value IE[Y]. The naive approach would be to approximate the expected value by the sample mean

Yˆ =

Yˆ1 +· · ·+ ˆYN

N . (2.52)

One obtains an unbiased estimator, as

IE[ ˆY] =IE[Y]. (2.53)

Now, we assume that with every outputYˆ_i we calculate the outputXˆ_i of another simulated quantity with known expected value IE[ ˆX_i] = IE[X]. For any constantb,

32

(35)

2.9 Variance Reduction Techniques we compute

Yˆ_i(b) = ˆY_i−b( ˆX_i−E[X]) (2.54) and

Yˆ(b) = ˆY −b( ˆX−E[X]) = 1 N

N

X

i=1

( ˆY_i−b( ˆX_i−E[X])) = 1 N

N

X

i=1

Yˆ_i(b). (2.55)

It can be easily verified that

IE[ ˆY(b)] = IE[ ˆY]−bIE[ ˆX−E[X]] =IE[ ˆY] =IE[Y] (2.56) and

Var[ ˆYi(b)] = Var[ ˆYi−b( ˆXi−IE[X])]

=σ_Y² −2bσ_Xσ_Yρ_XY +b²σ²_X =:σ²(b) (2.57) withσ_X² = Var[X], σ_Y² = Var[Y]and the correlation coefficient ρ_XY betweenX andY.

For the variance of the control variate estimatorYˆ(b), one obtainsσ²(b)/N. For the variance of the naive estimator,σ_Y²/N results. Consequently, variance reduction is achieved if the condition

bσ_X <2σ_Yρ_XY (2.58)

holds.

Differentiating equation (2.57) with respect tobyields the optimal coefficient b^∗ = Cov[X, Y]

Var[X] . (2.59)

Using this coefficient b^∗, for the variance reduction ratio, i.e. the ratio between the variance of the naive estimator and the variance of the control variate estimator, one obtains

Var[ ˆY]

Var[ ˆY(b)] = (1−ρ²_XY)⁻¹. (2.60) Thus, a strong correlation betweenXandY can lead to significant variance reduction [11].

Antithetic Variates The antithetic variates approach aims to introduce a negative correlation between the outputs of the simulation resulting in variance reduction.

Again, the objective is to estimate the expected value IE[Y]of a random variable

(36)

Y. To run the simulation, we produce a sequence ofN i.i.d. pairs of observations ( ˆY₁,Yˆ₁⁰),( ˆY₂,Yˆ₂⁰), . . . ,( ˆY_N,Yˆ_N⁰) of the random variable Y. For any given i, both constituents of the corresponding pair have the same distribution, however, generally they are correlated.

The antithetic variates estimator is defined as the sample mean of all2N observations ofY:

Yˆ_AV = 1 2N

N

X

i=1

Yˆ_i+

N

X

i=1

Yˆ_i⁰

!

= 1 N

N

X

i=1

Yˆ_i+ ˆY_i⁰ 2

!

(2.61)

Thus, the estimator can be interpreted as the sample mean ofN i.i.d. observations Yˆ₁+ ˆY₁⁰

2

! ,

Yˆ₂+ ˆY₂⁰ 2

! , . . . ,

Yˆ_N + ˆY_N⁰ 2

!

. (2.62)

Variance reduction is achieved, if for anyi

Var[ ˆY_i+ ˆY_i⁰]<2 Var[ ˆY_i] (2.63) holds, i.e. if the antithetic variates estimator consisting ofN pairs of observations ofY has lower variance than the naive estimator consisting of2N observations of Y.

For arbitrary( ˆY_i,Yˆ_i⁰), we have

Var[ ˆY_i+ ˆY_i⁰] = Var[ ˆY_i] + Var[ ˆY_i⁰] + 2 Cov[ ˆY_i,Yˆ_i⁰]

= 2 Var[ ˆY_i] + 2 Cov[ ˆY_i,Yˆ_i⁰]. (2.64) Thus, variance reduction is achieved if the constituents of the pairs are negatively correlated for any giveni, i.e. if

Cov[ ˆY_i,Yˆ_i⁰]<0. (2.65) The question arises, how the correlation betweenYˆ_iandYˆ_i⁰can be achieved without disturbing the independence between pairs. For observations of a uniformly distributed random variable U over [0,1], it is possible to calculate the sequence 1 − Uˆ₁,1 − Uˆ₂, . . . ,1 − Uˆ_N from the simulated sequence Uˆ₁,Uˆ₂, . . . ,Uˆ_N. For observations of a standard normally distributed random variableZ, “flipping” the sign by calculating the sequence−Zˆ₁,−Zˆ₂, . . . ,−Zˆ_N from the simulated sequence Zˆ₁,Zˆ₂, . . . ,Zˆ_N fulfills this condition [11].

34

(37)

2.9 Variance Reduction Techniques Stratified Sampling The stratified sampling approach partitions the sample space AintoK disjoint subsets A₁, . . . , A_K withS

iA_i = A. In a random Monte-Carlo simulation of Y, a sequence of independent Yˆ₁,Yˆ₂, ...,Yˆ_N with the same distribution asY is generated. In a stratified sampling approach, one decides upfront, what share of observations will be drawn from each of the subsetsAi, the so-calledstrata.

So, each observation drawn from a given stratumAihas the distribution ofY conditioned onY ∈A_i.

Again, the aim is to estimate IE[Y]of a random variableY withY :A→R: IE[Y] =

K

X

i=1

P(Y ∈A_i)IE[Y |Y ∈A_i] =:

K

X

i=1

p_iIE[Y |Y ∈A_i] (2.66)

The simplest approach to stratified sampling is to choose the portion of observations drawn from a stratum A_i to correspond to the theoretical probability p_i = P(Y ∈ A_i), i.e.N_i = p_iN. For each giveni, a sequenceYˆ_i1,Yˆ_i2, . . . ,Yˆ_iN_i of i.i.d. observations is drawn from the conditional distribution ofY onY ∈Ai. The estimator

Yˆ_SS =

K

X

i=1

p_i· 1 N_i

Ni

X

j=1

Yˆ_ij = 1 N

K

X

i=1 Ni

X

j=1

Yˆ_ij (2.67)

is an unbiased estimator for IE[Y].

The concept is most easily demonstrated by an example. Consider a Monte-Carlo simulation drawingN uniformly distributed random numbers from the intervalA= (0,1). Taking the limitN → ∞, the fraction of observations in any intervalAi ⊂A becomes infinitesimally close to the corresponding probability p_i. However, for finite sample sizesN, this is generally not the case. E.g., when choosingN = 1000 andK = 10 commensurate subsets A_i withp_i = 1/K = 0.1, in many cases not each of the subsetsA_i contains exactlyN_i =p_iN = 100observations ofY.

In contrast, in a stratified Monte-Carlo simulation, the strata A_i, i ∈ K are defined in advance andN_i observations are sampled fromY conditional onY ∈ A_i for eachi. Returning to the above example, per definition, each stratum contains 100observations ofY. The stratified approach yields a better approximation of the distribution from which the random numbers are drawn than the naive approach.

Compared to the naive estimator, the stratified sampling estimator reduces the variability of observationsacross strata. In contrast, the variancewithinone given stratum is not altered.

For advanced approaches, a detailed analysis of the variance reduction achieved and optimal stratification directions, refer to [11].

(38)

Latin Hypercube Sampling The Latin hypercube sampling approach extends the stratified sampling approach to high-dimensional Monte-Carlo simulations. In principle, one can apply the stratified sampling approach also to high-dimensional models. However, in such situations the computational effort increases exponen- tially with the number of dimensions d. The Latin hypercube sampling approach mitigates this computational shortcoming by stratifying only the marginal distributions of the joint multivariate distribution.

We introduce the method based on the example of thed-dimensional unit hypercube. For each dimensioni = 1,2, . . . , d, a stratified sample Vˆ_i⁽¹⁾,Vˆ_i⁽²⁾, . . . ,Vˆ_i^(K) is simulated from theK-fold equidistantly stratified unit interval. Arranging thed samples in columns permits to interpret the rows of the resulting table as points in the subcubes along the diagonal of the d-dimensional unit hypercube. In the next step, the entries of each column are permuted randomly. The resulting points are uniformly distributed over the unit hypercube.

The permutation introduces a dependence between theKpoints. Therefore, care must be taken when implementing the Latin hypercube sampling approach: it is permissible to introduce dependenceacrosspaths of a Brownian motion. However, it is not permissible to introduce dependence between the random numbers within one path as this can strongly bias the result [11].

Moment Matching The moment matching approach transforms sample paths to match the moments of the distribution from which sample paths are simulated.

For instance, for a non-dividend paying underlying S(T), the identity IE[S(T)] = e^rTS(0)holds. SimulatingN outputsSˆ₁(T),Sˆ₂(T), . . . ,Sˆ_N(T)fromS permits to calculate the sample mean

S(Tˆ ) = 1 N

N

X

i=1

Sˆ_i(T). (2.68)

The discounted estimator almost surely deviates from the current value, i.e.

e^−rTS(Tˆ )6=S(0). (2.69) However, by conducting the transformation

Sˆ_i⁰(T) =S_i(T)IE[S(T)]

Sˆ , (2.70)

the sample mean of the generated sequenceSˆ₁⁰(T),Sˆ₂⁰(T), . . . ,Sˆ_N⁰ (T)exactly match-

36

(39)

2.9 Variance Reduction Techniques

es the population mean IE[S(T)][11].

Weighted Monte-Carlo The weighted Monte-Carlo approach is a method similar to the matching underlying assets approach. However, instead of changing the simulated valuesSˆ₁(T),Sˆ₂(T), . . . ,Sˆ_N(T), weightsw₁, w₂, . . . , w_N are introduced that ensure [11]

N

X

i=1

wiSˆi(T) = e^rTS(0). (2.71)

Importance Sampling The importance sampling approach introduces a change of the probability measure from which the sample paths are generated, thus reducing the variance of Monte-Carlo estimators. To avoid redundancies within this work, we refer to the broader discussion of importance sampling in section 6.2.

Comparison of Variance Reduction Techniques The different approaches presented in this section differ with regards to their effectiveness and the complexity involved. The trivial antithetic variates approach yields limited variance reduction.

In uncommon situations it can increase variance. Control variates, moment matching and weighted Monte-Carlo generally perform better than antithetic variates with slightly increased complexity. The stratified sampling approach and its generaliza- tion, the Latin hypercube sampling, can be tailored to the specific model simulated, thus achieving significant variance variation but involving increased complexity.

Importance sampling approach is the mathematically most challenging approach.

It bears the risk that the empirical variance of the sample can be increased signif- icantly, when not handled properly. However, if handled correctly, it enables very strong variance reduction results. The focus of this work is on importance sampling.

(40)

(41)

Part II

Results and Discussion

(42)

(43)

3

An Application of the Put-Call-Parity to Variance reduced Monte-Carlo Option Pricing

3.0 Abstract

The standard error of Monte Carlo estimators for derivatives typically decreases at a rate ∝ 1/√

N whereN is the sample size. To reduce empirical variance for estimators of several in-the-money options, an application of the put-call-parity is analyzed. Instead of directly simulating a call option, first the corresponding put option is simulated. By employing the put-call-parity, the desired call price is calculated. Of course, the approach can also be applied vice versa. By employing this approach for in-the-money options, significant variance reductions are observed.

3.1 Introduction

It is a well known-property of certain options that under the assumption of arbitrage free markets a put-call-parity holds, relating the price of a put option to a call option. The put-call-parity can be derived based on the assumption that two financial products with the same pay-off and the same risk-profile must have the same fair market value. Otherwise, arbitrage would be possible [3, 10].

As described by Reider (1994) [45], put-call-parities can be applied to reduce the variance of Monte-Carlo estimators for option prices. E.g., when pricing an in-the- money call, one can instead simulate the price of a put option with same parameters and independent variables. The put-call-parity then yields the price of the in-the- money call, while the variance of the estimator can be dramatically reduced.

This paper is organized as follows: Monte-Carlo simulations of the Feynman-Kac formula as an approach to option pricing are introduced in section 3.2. In section 3.3 put-call-parities for European, arithmetic Asian, digital and basket options are

(44)

3 An Application of the Put-Call-Parity to Variance reduced Monte-Carlo. . . derived. Results for variance reduced Monte-Carlo simulations of these options are presented in section 3.4 and subsequently discussed in section 3.5. Section 3.6 concludes this paper.

3.2 Monte-Carlo Simulation of the Feynman-Kac Formula

Black-Scholes PDE In 1973, Fischer Black and Myron Scholes published an important contribution to option pricing theory [1]. Their option pricing model had tremendous impact on the further development of the trading of financial derivatives. Their work was honored with a noble price in economics later [3]. To derive the central constituent of their model, Black and Scholes derived a partial differential equation resulting from a no-arbitrage argument¹.

Thehedge portfolio

V =C+τ S (3.1)

consists of one call and τ units of the underlying with C ≡ C(S(t), t). The differential dV can be calculated by applying It¯o’s formula²:

dV =dC+τdS

=C_SdS+C_tdt+ 1

2C_SSdS²+ 1

2C_ttdt²+C_StdSdt+τdS

= (CS+τ)dS+

Ct+1 2CSSg²

dt+O dt³²

(3.2)

Here, the differential equation

dS =fdt+gdW (3.3)

and dW² ∝dtwere used. Terms of the order dt³² or higher were suppressed.

By choosingτ =−C_S the hedge portfolioV becomes riskless: The differential dV then depends only on deterministic terms as the first bracket term in the last row

1The representation follows Singer (1999) [13].

2In this paper, the abbreviation^∂C_∂S ≡CS is frequently used.

42

(45)

3.2 Monte-Carlo Simulation of the Feynman-Kac Formula

of equation (3.2) vanishes. In consequence,V must earn the risk-free rater:

dV =

C_t+ 1 2C_SSg²

dt

=! rVdt

= (rC−rC_SS)dt

(3.4)

In the first row, the term (3.2) was repeated and in the third row, definition (3.1) was inserted. Rearranging elements and choosing a geometric Brownian motion (f =rSandg =σS) yields the well known Black-Scholes differential equation

Ct+rSCS+1

2σ²S²CSS −rC = 0. (3.5) Black and Scholes (1973) [1] presented a closed form solution for an European call with boundary conditionC_T = (S_T −K)⁺:

C(S(0),0) =S(0)Φ (d₁)−Ke^−rTΦ (d₂) (3.6) with

d_1,2 =

ln^S(0)_K +

r± ^σ₂² T σ√

T .

Feynman-Kac Formula An alternative approach to option pricing was introduced by Cox and Ross in 1976 [2]. The authors state that in a risk-neutral world, the expected rates on the underlying and on the option must equal the risk-free rate:

IE[S(T)/S(0)|S(0) = S₀] =e^rT

IE[C(S(T), T)/C(S(0),0)|S(0) = S0] =C(S(0),0)⁻¹IE[C(S(T), T)|S(0) =S0]

=e^rT

(3.7)

Setting C(S(T), T) = h(S(T)) and rearranging yields the Feynman-Kac for-

(46)

3 An Application of the Put-Call-Parity to Variance reduced Monte-Carlo. . . mula for option pricing³:

C(S(0),0) = e^−rTIE[h(S(T))|S(0) =S0]

=e^−rT Z

h(S(T))p(S(T), T|S(0),0)dS(T) (3.8) The advantage of this representation is that for more difficult boundary condi- tionsh(S(T)), where no analytical solution to the Black-Scholes PDE (3.5) exists, option prices can be approximated by Monte-Carlo simulations. To obtain a Monte- Carlo estimator, N i.i.d. random numbers are drawn from the transition density p(S(T), T|S(0),0)yielding

Cˆ ≡C(S(0),ˆ 0) = 1 N

N

X

i=1

Cˆ_i (3.9)

with standard error

δCˆ = v u u t

1 N(N −1)

N

X

i=1

Cˆ_i²−NCˆ²

!

. (3.10)

Path-dependent Options A similar approach can be applied to path-dependent options [46]. Here, we consider an arithmetic Asian option depending on the entire path of the underlying:

h(S(T), Y(T), T)≡C(S(T), Y(T), T) =

Y(T) T −K

+

(3.11) with

Y(t) :=

Z t 0

S(u)du, fort∈[0, T]. (3.12) Differentiating equation (3.12) yields

dY =S(t)dt. (3.13)

A partial differential equation for the option price can be derived as above apply-

3This representation follows Singer (1999) [13]. It can be shown that this intuitive representation of the option price as discounted expected value of the pay-off function solves the partial differential equation (3.5). For further details on the Feynman-Kac formula see [5, 10].

44