Monte-Carlo-Simulation zur Optionsbewertung - uni-marburg.de · Monte-Carlo-Simulation zur...

Monte-Carlo-Simulation zurOptionsbewertung

Fachbereich Mathematik und Informatikder Philipps-Universität Marburg

Bachelorthesiszur Erlangung des akademischen Grades

Bachelor of Science

vorgelegt von

Nils Benedikt Kollmargeboren am 13.08.1993 in Adenau

im August 2016

Erstprüfer: Prof. Dr. Dr. Marcus Porembski

Inhaltsverzeichnis1 Einleitung und Motivation 13

1.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.2 Approximation von π . . . . . . . . . . . . . . . . . . . . . . 14

2 Stochastische und Finanzmathematische Grundlagen 172.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Optionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3 Stochastische Grundlagen . . . . . . . . . . . . . . . . . . . 192.4 Wiener Prozess . . . . . . . . . . . . . . . . . . . . . . . . . 212.5 Stochastische Differentialgleichungen und Ito-Lemma . . . . 26

3 Monte-Carlo-Simulation 313.1 Grundidee . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2 Optionspreisbewertung . . . . . . . . . . . . . . . . . . . . . 343.3 Pseudozufallszahlen . . . . . . . . . . . . . . . . . . . . . . . 38

3.3.1 Linearer Kongruenzgenerator . . . . . . . . . . . . . 393.3.2 Weitere Generatoren . . . . . . . . . . . . . . . . . . 42

3.4 Standardnormalverteilte Zufallszahlen . . . . . . . . . . . . . 433.4.1 Invertierung . . . . . . . . . . . . . . . . . . . . . . . 443.4.2 Box-Muller-Algorithmus . . . . . . . . . . . . . . . . 463.4.3 Polar-Algorithmus . . . . . . . . . . . . . . . . . . . 46

4 Varianzreduktion 494.1 Antithetic Variates . . . . . . . . . . . . . . . . . . . . . . . 494.2 Control Variates . . . . . . . . . . . . . . . . . . . . . . . . . 51

5 Diskretisierungsverfahren 555.1 Euler-Maruyama-Verfahren . . . . . . . . . . . . . . . . . . . 555.2 Milsteinverfahren . . . . . . . . . . . . . . . . . . . . . . . . 55

6 Bewertung Asiatischer Optionen mit Monte-Carlo-Simulation 616.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.2 Bewertung einer arithmetischen Asiatischen Option . . . . . 626.3 Varianzreduktion . . . . . . . . . . . . . . . . . . . . . . . . 63

7 Fazit 677.1 Zusammenfassende Bewertung und Fazit . . . . . . . . . . . 67

3

7.2 Persönliches Schlusswort . . . . . . . . . . . . . . . . . . . . 68

8 Appendix: Matlabcodes 698.1 Simulation von stochastischen Prozessen . . . . . . . . . . . 698.2 Zufallszahlen . . . . . . . . . . . . . . . . . . . . . . . . . . 728.3 Invertierung . . . . . . . . . . . . . . . . . . . . . . . . . . . 738.4 Polar-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . 748.5 Vergleich der Konfidenzintervalle . . . . . . . . . . . . . . . 758.6 Vergleich der approximativen Optionswerte der einfachen Monte-

Carlo-Simulation mit varianzreduzierten Methoden . . . . . 788.7 Vergleich der Kurssimulationen mit dem Euler-Maruyama-

Verfahren und mit dem Milstein-Verfahren . . . . . . . . . . 818.8 Bewertung Asiatischer Optionen . . . . . . . . . . . . . . . . 83

Abbildungsverzeichnis2.1 Simulation von drei Wiener Prozessen . . . . . . . . . . . . . 242.2 Simulation von drei Brownschen Bewegungen mit Drift und

Volatilität . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.3 Zerlegung einer Geometrischen Brownschen Bewegung . . . . 292.4 Simulation von drei Geometrischen Brownschen Bewegungen 30

3.1 3-dimensionale Zufallsvektoren mit RANDU I . . . . . . . . 413.2 3-dimensionale Zufallsvektoren mit RANDU II . . . . . . . . 423.3 Zufallszahlen mit dem Beasley-Springer-Moro-Algorithmus . 453.4 Zufallszahlen mit dem Polar-Algorithmus . . . . . . . . . . . 47

4.1 Konfidenzintervalle des Standard Monte-Carlo-Schätzers unddes Schätzers mit Control Variates . . . . . . . . . . . . . . 53

4.2 Approximative Optionsbewertung mit einfacher Monte-Carlo-Methode und varianzreduzierten Monte-Carlo-Methoden . . 54

5.1 Simulation eines Kurses mit dem Milsteinverfahren und mitdem Euler-Maruyama-Verfahren zu jeweils 3 Pfaden . . . . . 59

6.1 Simulation des Werts einer Asiatischen Option mit varianzre-duzierten Schätzern und dem Standard Schätzer . . . . . . . 66

5

Verzeichnis der Matlabcodes1.1 Monte-Carlo-Simulation zur Berechnung von π . . . . . . . . 15

8.1 Simulation der in Kapitel 2 vorgestellten Algorithmen . . . . 698.2 Zerlegung der Geometrischen Brownschen Bewegung . . . . 718.3 Beispiel: RANDU . . . . . . . . . . . . . . . . . . . . . . . . 728.4 Multiplikativer Kongruenzgenerator . . . . . . . . . . . . . . 728.5 Beasley-Springer-Moro-Algorithmus . . . . . . . . . . . . . . 738.6 Polar-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . 748.7 Vergleich der Konfidenzintervalle mit Control Variate . . . . 758.8 Varianz des Schätzers mit Control Variate . . . . . . . . . . 768.9 Varianz des einfachen Monte-Carlo-Schätzers . . . . . . . . . 778.10 Approximativ optimales b . . . . . . . . . . . . . . . . . . . 788.11 Monte-Carlo-Simulation zur Optionsbewertung einfach, mit

Control Variates und mit Antithetic Variates . . . . . . . . . 788.12 Monte-Carlo-Simulation zur Optionsbewertung mit Antithe-

tic Variates . . . . . . . . . . . . . . . . . . . . . . . . . . . 798.13 Monte-Carlo-Simulation zur Optionsbewertung mit dem ein-

fachen Monte-Carlo-Schätzer . . . . . . . . . . . . . . . . . . 808.14 Monte-Carlo-Simulation zur Optionsbewertung mit Control

Variates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 808.15 Simulation eines Kurses mit dem Euler-Maruyama-Verfahren

und mit dem Milsteinverfahren . . . . . . . . . . . . . . . . 818.16 Simulation des Werts einer Asiatischen Option . . . . . . . . 83

7

Liste der Algorithmen1 Simulation eines Wiener Prozesses . . . . . . . . . . . . . . . 232 Simulation einer Brownschen Bewegung mit Drift µ und Vola-

tilität σ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253 Simulation eines Kurses auf Grundlage der Geometrischen Brown-

schen Bewegung . . . . . . . . . . . . . . . . . . . . . . . . . 28

4 Grundalgorithmus zur Approximation des Optionspreises . . 375 Approximation eines pfadabhängigen Optionspreises . . . . . 386 Linearer Kongruenzgenerator . . . . . . . . . . . . . . . . . . 39

7 Simulation eines Kurses auf Grundlage der Geometrischen Brown-schen Bewegung mittels Milstein-Verfahren . . . . . . . . . . 57

8 Approximation des Optionspreises einer Asiatischen Option mitEuler-Maruyama-Verfahren . . . . . . . . . . . . . . . . . . . 62

9

Tabellenverzeichnis1.1 Monte-Carlo-Simulation zur Bestimmung von π . . . . . . . 14

5.1 Vergleich der Endkurse . . . . . . . . . . . . . . . . . . . . 59

6.1 Optionswerte und Varianzen mit drei Schätzern . . . . . . . 64

11

Wahrlich es ist nicht das Wissen,sondern das Lernen, nicht das Besitzensondern das Erwerben, nicht dasDa-Seyn, sondern das Hinkommen, wasden grössten Genuss gewährt.

Carl Friedrich Gauß 1

1 Einleitung und Motivation

1.1 EinleitungWer sich für Optionen auf dem Finanzmarkt interessiert, dem ist immerauch an dem Wert der Option gelegen. Denn, wenn der Marktpreis derOption höher ist als der Wert, den man selbst der Option zuschreibt, sosollte man diese Option nicht kaufen. Also versuchen wir, Optionen ganzallgemein zu bewerten. Bei sogenannten Plain Vanilla Optionen, also Op-tionen mit den einfachsten Auszahlungsarten, gelingt dies sehr einfach mitder Black-Scholes-Formel 2. Betrachten wir aber Optionen mit nur geringfü-gig komplizierteren Auszahlungsfunktionen, so müssen wir meist auf andereMethoden zurückgreifen. Hier kommt die Monte-Carlo-Simulation ins Spiel,die eine relativ gute Näherung des exakten Wertes liefert und dennoch einesehr einfache Struktur bietet. Weil der Grundalgorithmus bestehen bleibt,können wir mit verhältnismäßig geringem Aufwand für die jeweilige Auszah-lungsfunktion der Option die Simulation anpassen.Damit ist auch die Zielstellung für diese Arbeit klar: Wir wollen möglichsteinfach und mit geringem Rechenaufwand komplizierten, exotischen Optio-nen einen Wert zuordnen. Explizit heißt das hier, dass wir die GeometrischeBrownsche Bewegung wie Glasserman [3] als Modell für die Kursentwick-lung des Underlyings herleiten und darauf aufbauend den Grundalgorith-mus der Monte-Carlo-Simulation zur Optionsbewertung nach Glasserman[3] formulieren. Auf dem Weg dorthin - und hier kommen wir auf das Zi-tat Gauss’ zurück - werden uns einige Möglichkeiten zur Verfeinerung derStandardmethode auffallen. Diese werden wir vor allem in Form der Va-rianzreduktionstechniken diskutieren. Außerdem werden wir die Kurse mitstochastischen sowie numerischen Methoden simulieren und diese Methodeninsbesondere nach ihrer Umsetzbarkeit beurteilen. Letztendlich soll uns dieseHerangehensweise dazu führen, eine Asiatische Option bewerten zu können.Gleichzeitig werden wir die Genauigkeit anhand von Konfidenzintervalleneinschätzen und eine Aussage über die Effizienz der Methode machen, dasheißt wie schnell unser Monte-Carlo-Algorithmus den Optionspreis zu einem

1Carl Friedrich Gauss hat diesen Satz in einem Brief vom 2.September 1808 an WolfgangBolyai verwendet, als er nach Bolyais wissenschaftlicher Arbeit fragte.Quelle: Briefwechsel zwischen Carl Friedrich Gauß und Wolfgang Bolyai, herausgege-ben von F. Schmidt, P. Stäckel, Leipzig 1899, S. 94

2siehe Günther und Jüngel [4, Satz 4.12]

13

KAPITEL 1. EINLEITUNG UND MOTIVATION

bestimmten Konfidenzintervall berechnet.Um die Idee der Monte-Carlo-Simulation zu veranschaulichen, beginnen wirmit einem mathematisch eher einfaches Beispiel: der numerischen Berech-nung von π.

1.2 Approximation von πDie populärste Herangehensweise die Monte-Carlo-Simulation zu illustrie-ren ist wohl die näherungsweise Bestimmung von π am Einheitskreis. Dabeiwerden m Paare von Zufallszahlen simuliert, die auf dem Einheitsintervallgleichverteilt sind. Nennen wir diese Zufallsvariablen X und Y . Dement-sprechend betrachten wir nur den positiven Viertelkreis. Mit der Flächen-formel des Kreises A = π · r2 folgt für den Einheitskreis A = π. Sei B =(X(ω), Y (ω)) | X(ω)2 + Y (ω)2 ≤ 1 die Menge, die jene Realisationen re-präsentiert, die in dem Viertelkreis liegen, dann gilt für die Kardinalität vonB:

|B| · 4m

= π

Anhand dieses Beispiels erkennt man die Grundidee der Monte-Carlo-Simulation: Der Erwartungswert eines Zufallsexperiments wird durch Simu-lation der Realisierungen der Zufallsvariablen geschätzt. Außerdem ist be-reits hier eine wichtige Eigenschaft von Monte-Carlo-Verfahren offensichtlich,nämlich dass die Güte der Näherung mit wachsender Anzahl an Simulationensteigt.

π 3,14159265359Anzahl der Simulationen m Näherung von π absoluter Fehler

10 3,6 0,45840734641100 3,28 0,138407346411000 3,2320 0,0904073464110000 3,1404 0,00119265359100000 3,1440 0,002407346411000000 3,1433 0,0017073464110000000 3,1414 0,00019265359

Tabelle 1.1: Monte-Carlo-Simulation zur Bestimmung von π

Um dieses Verhalten zu illustrieren ist in Matlabcode 1.1 die Simulationdes Beispiels implementiert und wurde für die obige Tabelle mit unterschied-licher Anzahl an Simulationen ausgeführt. Diese sehr einfache Simulationkann uns bereits eine auf drei Nachkommastellen genaue Approximation fürπ berechnen.

14

KAPITEL 1. EINLEITUNG UND MOTIVATION

1 rand ( ’ s t a t e ’ , 3 ) ;2 % Berechnung von Pi m i t t e l s Monte−Carlo−Simulat ion3 % n : S imu l a t i on sha eu f i g k e i t ( −> Unendlich )4 n=10000000;5 % auf (0 , 1 ) g l e i c h v e r t e i l t e Zu f a l l s z ah l e n6 x=rand (n , 1 ) ; % Zu f a l l s v e k t o r der Laenge n7 y=rand (n , 1 ) ; % Zu f a l l s v e k t o r der Laenge n8

9 punkte_im_vierte lkre i s = sum(x.^2+y.^2<1) ;10 pi = ( punkte_im_vierte lkre i s /n) ∗4

Matlabcode 1.1: Monte-Carlo-Simulation zur Berechnung von π

15

2 Stochastische undFinanzmathematischeGrundlagen

2.1 EinleitungIn diesem Kapitel wollen wir erklären, was eine Option als Finanzderivat ist.Außerdem ist es vonnöten, die stochastischen Grundlagen für den Monte-Carlo-Schätzer zu legen und, kurz gesagt, die Voraussetzungen zur Simu-lation von Kursverläufen zu schaffen. Eine zentrale Rolle dabei spielt dieGeometrische Brownsche Bewegung, worauf dieses Kapitel hinarbeitet. Da-bei beruht der erste Teil dieses Kapitels, der Optionen behandelt, auf demVorlesungsskript [20] zur Finanzmathematik von Prof. Dr. Dr. Marcus Po-rembski, der stochastische Teil vor allem auf der Literatur von Glasserman[3], Günther und Jüngel [4], sowie Müller-Gronbach et al. [18]. Der Abschnittüber stochastische Differentialgleichungen benötigte Erkenntnisse der spezi-fischen Literatur von Øksendal [19].

2.2 OptionenEine Option ist ein Vertrag, der das Recht beinhaltet einen Basiswert (engl.Underlying oder Underlying Asset) zu einem im Vorhinein vereinbarten PreisX, dem Strike, zu kaufen oder zu verkaufen. Eine Option hat (meist) einenFälligkeitszeitpunkt T > 0, bis zu diesem die Optionen ausgeübt werdenkann. Handelt es sich um ein Kaufrecht, so sprechen wir von einer Callopti-on, bei einem Verkaufsrecht von einer Putoption. Optionen werden außerdemnoch darin unterschieden, zu welchen Zeitpunkten sie ausgeübt werden kön-nen. Sogenannte Europäische Optionen können nur zum FälligkeitszeitpunktT ausgeübt werden, das heißt der Kursverlauf ist für diese Art von Optionenunwichtig. Wichtig ist nur, welchen Wert die Option am Fälligkeitszeitpunkthat. Demgegenüber stehen Amerikanische Optionen, die zu jedem beliebigenZeitpunkt t mit t ∈ [0, T ] ausgeübt werden können. Eine Hybridform sindBermudaoptionen, die zu mehreren vereinbarten Zeitpunkten ausgeübt wer-den können. Bermudaoptionen gehören zu den sogenannten exotischen Op-tionen. Auf eine andere exotische Option, die Asiatische, kommen wir später

17

KAPITEL 2. STOCHASTISCHE UND FINANZMATHEMATISCHEGRUNDLAGEN

zurück. Bei den Amerikanischen Optionen und Bermudaoptionen, verfälltdie Option mit der Ausübung, wenn diese vor dem Fälligkeitszeitpunkt er-folgt ist.

Wenn von einem Underlying die Rede ist, so meint man sehr häufig eineAktie. Es gibt aber auch Optionen, die zum Beispiel auf Indizes, Rohstoffeoder Währungen abgeschlossen werden. In dieser Arbeit beziehen wir unsaber zumeist auf beliebige Assets und wenn von einer Aktie die Rede ist, solässt sich dies auch auf allgemeine Assets anwenden, sofern diese am Marktgehandelt werden.

PreisDer Holder der Option, also derjenige, der die Option erworben hat, hatdie Möglichkeit sein Recht auszuüben, aber nicht die Pflicht auszuüben - imGegensatz zu der Partei, die die Option verkauft hat. Diese Partei muss, so-fern die Option durch den Holder ausgeübt wird, den Kauf oder Verkauf desUnderlying durchführen. Hier erkennen wir, dass der Holder einen gewissenVorteil hat, da er wählen kann, ob das Ausübungsrecht genutzt wird odernicht. Das legt den Schluss nahe, dass eine Option einen Preis haben muss,den der Holder zahlt. Ansonsten könnte man durch den Erwerb einer OptionArbitrage, also risikolosen Gewinn, generieren. Auf dem Finanzmarkt wer-den Optionen natürlich mit Preisen gehandelt. Diese Preise kommen auf demMarkt über Angebot und Nachfrage zustande. Aber wie sollte der Options-preis aussehen, wenn er „fair“ für beide Parteien sein soll? Dieser Frage gehenwir in Abschnitt 3.2 genauer nach. An dieser Stelle ist es wichtig zu diffe-renzieren: Der Preis einer Option wird durch Marktmechanismen bestimmt.Den Wert, also den „fairen“ Preis, können wir unter Umständen berechnen.Auf einem vollständigen Markt fielen diese Preise zusammen.

MarktannahmenUm die stochastischen Zusammenhänge auf dem Finanzmarkt darzustellen,müssen bestimmte Annahmen getroffen werden. Die Wichtigste ist die Ar-bitragefreiheit, die aussagt, dass auf dem Markt, den wir betrachten, keinrisikoloser Gewinn generiert werden kann. Außerdem unterstellen wir in derganzen Arbeit eine stetige Verzinsung und einen zeitunabhängigen Zinssatz.Wir gehen grundsätzlich davon aus, dass keine Dividenden ausgeschüttetwerden und verlangen beliebige Geldanlage und Kreditaufnahme.Zudem wenden wir ein sehr wichtiges Konzept der Finanzmathematik an,die risikoneutrale Bewertung. Das bedeutet, wir diskontieren die erwartetenAuszahlungen mit dem risikolosen Zins r ab. Dem Käufer käme somit eine

18


erwartete Rendite µ in Höhe des risikolosen Zinssatzes zu, wir setzen also

µ = r

. Weil der Holder demzufolge weder risikoavers noch risikoaffin handelt, spre-chen wir von risikoneutraler Bewertung. Die tiefergehende Theorie der risiko-neutralen Bewertung und ihr Sinn in der Anwendung wird von Glasserman[3, s.25f] genauer diskutiert. Für uns soll an dieser Stelle genügen, dass wirdurch Diskontierung mit dem risikolosen Zins einen „fairen“ Preis berechnenkönnen. Im vorherigen Abschnitt haben wir gesehen, dass ein vollständigerMarkt diese Vorgehensweise rechtfertigt.Außerdem unterstellen wir, dass die Kurse der Assets, die wir betrachten,einer geometrischen Brownschen Bewegung folgen, wie wir sie im Abschnitt2.5 herleiten werden.

2.3 Stochastische GrundlagenUm die Monte-Carlo-Simulation und die Optionsbewertung mit ebenjenerSimulation zu verstehen, machen wir uns in diesem Abschnitt einige grundle-gende, stochastische Zusammenhänge klar. Dabei wird unser Ziel sein, einenAktienkursverlauf mathematisch zu beschreiben und einen solchen Verlaufzu simulieren. Wir vergegenwärtigen uns zunächst, was der Erwartungswertund die Varianz sowie die Kovarianz einer Zufallsgröße sind. Diese Definitio-nen tauchen in Kapitel 3 an prominenter Stelle wieder auf.

Definition 2.1 (Erwartungswert, Varianz, Kovarianz). Sei X eine reelleZufallsvariable mit integrierbarer Wahrscheinlichkeitsdichte f(x), dann heißt

E[X] :=∫ ∞−∞

xf(x)dx

der Erwartungswert von X und

V ar(X) := E[(X − E[X])2]

die Varianz von X.Sei darüber hinaus Y auch eine reelle Zufallsvariable und die Erwartungs-werte E[Y ] und E[XY ] existieren, dann heißt

Cov(X, Y ) := E[(X − E[X])(Y − E[Y ])]

die Kovarianz von X und Y .

Grundlegende Eigenschaften und Rechenregeln, wie die Linearität bezie-hungsweise lineare Transformation des Erwarungswertes und der Varianz,setzen wir als bekannt voraus. Auf die Beziehung zwischen Varianz und Ko-varianz gehen wir im folgenen Satz ein, da wir diesen für die Varianzreduktionnoch benötigen werden.

19


Satz 2.2. Seien X,Y zwei Zufallsvariablen, dann gilt für die Varianzen

V ar(X + Y ) = V ar(X) + V ar(Y ) + 2Cov(X, Y ) (2.1)

Für einen Beweis verweisen wir auf Henze [5, Satz 21.2].

Definition 2.3 (Normalverteilung). Eine Zufallsvariable X heißt normal-verteilt mit Erwartungswert µ ∈ R und Varianz σ2 > 0, wenn sie die Wahr-scheinlichkeitsdichte

f(x) = 1σ√

2πexp

(−1

2

(x− µσ

)2)

mit x ∈ R besitzt .Für eine normalverteilte Zufallsvariable mit Erwartungs-wert µ und Varianz σ2 schreiben wir X ∼ N (µ, σ2). Falls für X gilt, dassµ = 0 und σ2 = 1 ist, so sprechen wir von der Standardnormalverteilung undschreiben X ∼ N (0, 1).

Nachdem wir nun die Definition der allgemeinen Normalverteilung undder Standardnormalverteilung vorliegen haben, wollen wir noch herausfin-den, wie wir eine standardnormalverteilte Zufallsvariable aus einer normal-verteilten Zufallszahl gewinnen können. Dazu führen wir folgenden Satz an.

Satz 2.4. Sei X eine normalverteilte Zufallsvariable mit Erwartungswert µund Varianz σ2, dann ist die Zufallsvariable

Y = X − µσ

(2.2)

standardnormalverteilt, d.h. Y ∼ N (0, 1).

Beweis. Der Beweis beruht auf einer einfachen Substitution der Integralgren-zen der Verteilungsfunktionen, denn für die Wahrscheinlichkeit von Y = X−µ

σ

gilt:

P(X − µσ

≤ u)

= P (X ≤ σu+ µ)

=∫ σu+µ

−∞

1√2πσ

exp

(−1

2

(x− µσ

)2)dx

=∫ u

−∞

1√2πexp− 1

2y2dy (2.3)

= P (Y ≤ u)

mit Y ∼ N (0, 1). In der Umformung zur Gleichung 2.3 wurde x = µ + σysowie dx = d(µ+ σy) = σdy substituiert.

20


2.4 Wiener ProzessJedenfalls bin ich überzeugt, daß der[der Alte] nicht würfelt.

Albert Einstein 1

Mit diesem Ausspruch, der oftmals mit „Gott würfelt nicht!“ wiederge-geben wird, steigen wir in die Theorie des Wiener Prozesses ein. Auch wennsich Einstein auf die wahrscheinlichkeitstheoretische Erklärung der Quan-tenmechanik bezogen hat, so ist dieser Ausspruch für uns doch interessant,da wir im nachfolgenden Abschnitt eine zufällige Bewegung mathematischbeschreiben und numerisch simulieren werden. Albert Einstein selbst hat aufdiesem Feld Pionierarbeit geleistet und, wenn auch nicht ganz im Sinne desobigen Ausspruches, eine mathematische Erklärung für den Wiener Prozessgeliefert. Im späteren Verlauf sehen wir nämlich, dass wir zur Simulationtatsächlich in einem gewissen Sinne würfeln, wir arbeiten nämlich mit Zu-fallszahlen. Die zentralen Fragen dieses Abschnittes sollen sein: Was ist einWiener Prozess im mathematischen Sinne und wie können wir einen solchensimulieren?Um diese Fragen zu beantworten, müssen wir zunächst einen stochastischenProzess definieren, um im Folgenden weitere spezielle Prozesse zu betrachten.Auf deren Grundlage werden wir dann die Simulationsvorschriften herleiten.

Definition 2.5 (Stochastischer Prozess). Ein stetiger, stochastischer Pro-zess X(·, t), t ∈ [0,∞), ist eine Familie von Zufallsvariablen X : Ω ×[0,∞)→ R, bei der die Abbildung t 7→ X(ω, t) stetig für alle ω ∈ Ω ist.

Im Folgenden wird die kürzere Schreibweise X(·, t) = Xt genutzt. Wirsollten uns aber bewusst sein, dass wir einen stochastischen Prozess bezüg-lich eines Ereignisses ω betrachten. Wenn im Folgenden von einem Pfad dieRede ist, so ist die Realisation eines stochastischen Prozesses mit einem fes-ten ω, also X(ω, t) gemeint. In diesem Fall ist der stochastische Prozess einegewöhnliche Funktion, abhängig von der Zeit t. Die Pfade zu bestimmtenEreignissen eines Prozesses werden für die Simulation sehr wichtig sein.Betrachten wir nun einen speziellen stochastischen Prozess, der in verschie-denen Wissenschaften teilweise an sehr prominenter Stelle auftritt, wie inder Physik oder der Biologie.

Definition 2.6 (Wiener Prozess). Sei (Wt)t≥0 ein stetiger, stochastischerProzess der folgende Eigenschaften erfüllt:

(i) W0 = 0 mit Wahrscheinlichkeit 11 Albert Einstein in einem Brief an Max Born. Mit dem „Alten“ ist Gott gemeint.Quelle: Albert Einstein, Hedwig und Max Born: Briefwechsel 1916 - 1955, RowohltTaschenbuchverlag, 1972, S. 97f

21


(ii) Wt −Ws ist normalverteilt mit Erwartungswert 0 und Varianz t − s,für beliebige t und s, die die Ungleichung 0 ≤ s < t erfüllen.

(iii) Die Inkremente Wt2 −Wt1, Wt4 −Wt3 sind unabhängige Zufallsgrößenfür alle 0 ≤ t1 < t2 ≤ t3 < t4

Dann heißt Wt Wiener Prozess.

Nachdem Albert Einstein bereits 1905 (vgl. [2] sowie [22]) die BrownscheBewegung nutzte, um die Bewegung von Teilchen in Flüssigkeit zu erklären,bewies Norbert Wiener 1923 (siehe Müller-Gronbach et al. [18, Kapitel 4.5.2])die wahrscheinlichkeitstheoretische Existenz dieses stochastischen Prozesses.Der Name Brownsche Bewegung geht auf Robert Brown zurück, der 1828 dieBewegung von in Flüssigkeit gelösten Pollen beobachtete. Der Wiener Pro-zess ist als mathematisches Modell der Brownschen Bewegung zu verstehenund wird uns im weiteren Verlauf dieser Arbeit helfen, Aktienkurse (oder all-gemein Kurse von Assets) zu simulieren. Machen wir uns vorher noch einengrundlegenden Zusammenhang für den Wiener Prozess klar.

Folgerung 2.7. Seien t ≥ 0 und ∆t > 0, dann gilt, dass der Erwartungswertvon Wt gleich 0 ist, E[Wt] = E[Wt − W0] = 0 und dass die InkrementeWt+∆t −Wt normalverteilt mit Erwartungswert 0 und Varianz ∆t sind.

Mit diesen Erkenntnissen, die direkt aus der Definition folgen, können wireinen Versuch unternehmen, einen Wiener Prozess numerisch zu simulieren.Dafür sei t0 = 0 der Anfangszeitpunkt und ∆t die Schrittweite. Es ergebensich also die diskreten Zeitpunkte ti = i∆t für einen nichtnegativen Index iund die Diskretisierung des Wiener Prozesses:

Wti ≈i∑

k=1(Wtk −Wtk−1) (2.4)

Dabei ist zu bemerken, dass die Wtk −Wtk−1 jeweils nach Definition 2.6unabhängig voneinander sind und normalverteilt mit Erwartungswert 0 undVarianz ∆t. Das machen wir uns bei der Simulation von W (t) zunutze undfolgern mit Gleichung 2.2, dass

Zk =Wtk −Wtk−1√

∆t(2.5)

standardnormalverteilt ist. Stellen wir Gleichung 2.5 nachWtk um, haben wireine iterative Formel für den Wiener Prozess, abhängig von einer standard-normalverteilten Zufallszahl Zk. Nun können wir Algorithmus 1 zur Erzeu-gung von Realisationen eines Wiener Prozesses angeben. Dieser Algorithmusist, wie die zwei folgenden Algorithmen, in Matlabcode 8.1 realisiert. Die Si-mulation von drei unterschiedlichen Wiener Prozessen in Abhängigkeit von

22


Input : ∆t, Zi ∼ N (0, 1) für i = 1, . . . , nset W0 = 0for i = 1, . . . , n do

Wti = Wti−1 +√

∆tZiendOutput : Vektor W ∈ Rn+1, der den Pfad eines Wiener Prozesses

beinhaltetAlgorithmus 1 : Simulation eines Wiener Prozesses

der Anzahl der Zeitschritte können wir in Abbildung 2.1 sehen, später werdenwir aber deutlich mehr als nur drei Pfade zur Optionsbewertung nutzen. Manbeachte hierbei, dass wir mit einem etwaigen Zeitintervall T die Schrittweite∆t := T/n setzen sollten. Man kann aber auch eine bestimmte Schrittweitevorgeben und einen Pfad der Länge n∆t simulieren. Nun ist Algorithmus 1äquivalent zu der Aussage, die aus Gleichung 2.4 und 2.5 folgt:

Wti ≈i∑

k=1Zk√

∆t

Von entscheidender Bedeutung für die Genauigkeit der Diskretisierung isthierbei, wie groß die Schrittlänge ∆t ist. Wir werden später noch einmalsehen, wie eine genauere Diskretisierung aussehen könnte. Nachdem wir nunden Wiener Prozess simulieren können, betrachten wir darauf aufbauendeinen etwas allgemeineren Fall, der einen Parameter eines grundsätzlichenWachstums des Kurses und einen Parameter für eine gewisse Zufälligkeit umden Trendparameter beinhaltet. Wir werden diese Parameter im Folgendenµ und σ nennen.

23


Abbildung 2.1: Simulation von drei Wiener Prozessen

Definition 2.8. Für gegebene Paramter µ und σ > 0 nennen wir einenstochastischen Prozess Xt Brownsche Bewegung mit Drift µ und Volatilitätσ, wenn

Xt − µtσ

(2.6)

ein Wiener Prozess ist.

Stellen wir die Gleichung 2.6 nach Xt um, so definiert zu gegebenem Wie-ner Prozess Wt

Xt = µt+ σWt (2.7)

die Brownsche Bewegung mit Drift µ und Volatilität σ. Wegen der Eigen-schaften des Wiener Prozesses folgt, dass Xt nach N (µt, σ2t) verteilt ist. Mitder Tatsache, dass wir aus einer normalverteilten Zufallsgröße sehr einfach ei-ne standardnormalverteilte Zufallsgrößen erzeugen können, wie in Gleichung2.2, folgt, dass wir Xt wie den Wiener Prozess in Gleichung 2.5 simulierenkönnen. Es gilt, dassXtk−Xtk−1 normalverteilt mit Erwartungswert µ∆t undVarianz σ2∆t ist. Daraus folgt wieder mit der Transformation aus Gleichung2.2, dass die Zufallsgröße Zi standardnormalverteilt ist:

Zti :=(Xtk −Xtk−1)− µ∆t

σ√

∆t(2.8)

24


Input : ∆t, µ, σ, Zi ∼ N (0, 1) für i = 1, . . . , nset Xt0 = 0for i = 1, . . . , n do

Xti = Xti−1 + µ∆t+ σ√

∆tZiendOutput : Vektor X ∈ Rn+1, der den Pfad einer Brownschen Bewegung

enthältAlgorithmus 2 : Simulation einer Brownschen Bewegung mit Drift µ undVolatilität σ

Damit gilt bei gegebenem Startwert Xt0 Algorithmus 2.

Abbildung 2.2: Simulation von drei Brownschen Bewegungen mit Drift undVolatilität

Der stochastische Prozess Xt, der in Algorithmus 2 berechnet wird, löstdie stochastische Differentialgleichung

dXt = µdt+ σdWt (2.9)

näherungsweise. Für eine Implementation des Algorithmus steht der Pro-grammcode 8.1 zur Verfügung, mit welchem die Abbildung 2.2 erzeugt wur-de. Sie zeigt uns drei Pfade einer Brownschen Bewegung mit Drift 0,1 und

25


Volatilität 0,35. Diese Grafik stellt den wesentlichen Nachteil der BrownschenBewegung mit Drift und Volatilität heraus, der Kurs wird negativ. Diese Ei-genschaft ist bei einem Kursmodell zum Beispiel für Aktien oder Rohstoffe,Assets im Allgemeinen, nicht wünschenswert und wird im nächsten Schrittausgebessert. Im Vergleich zur Abbildung 2.1 kann man erkennen, dass dieKursverläufe einen geringeren Trendwachstum haben und nicht so weit aus-einander liegen.

2.5 Stochastische Differentialgleichungen undIto-Lemma

Nun wollen wir aus der Brownschen Bewegung die geometrische BrownscheBewegung herleiten, welche wir dann tatsächlich nutzen, um Aktienkurse zusimulieren. Dazu müssen wir uns aber etwas genauer mit stochastischen Dif-ferentialgleichungen und dem Ito-Lemma beschäftigen. Im Folgenden wirdfür stochastische Differentialgleichungen die Schreibweise SDE genutzt, wasfür stochastic differential equation steht. In dieser Arbeit können und sollenstochastische Differentialgleichungen nicht im Mittelpunkt stehen, deshalbsei für den Beweis des Ito-Lemmas und weitergehende Diskussion der The-matik auf Øksendal [19] sowie Günther und Jüngel [4] verwiesen. Klären wiralso zuerst, was eine SDE nach Ito ist. Zu gegebenem Wiener Prozess Wt,stochastischem Prozess Xt und zwei Funktionen a, b : R × R+ → R nennenwir

dXt = a(Xt, t)dt+ b(Xt, t)dWt (2.10)

mit der Integralschreibweise

Xt = Xt0 +∫ t

t0a(Xs, s)ds+

∫ t

t0b(Xs, s)dWs (2.11)

Ito-SDE oder stochastische Differentialgleichung nach Ito.Hierbei heißt a(Xt, t) Drift Term und b(Xt, t) Diffusion. Eine Lösung vonGleichung 2.11 heißt Ito-Prozess und ist eine Verallgemeinerung des Wiener-Prozesses. Dies ist direkt erkenntlich, wenn wir die triviale SDE dXt = dWt

betrachten, also ist Drift 0 und Diffusion 1. Zu bemerken ist auch, dass wirbeim zweiten Integral in Gleichung 2.11 kein Riemannsches Integral oderLebesgue-Integral nutzen können, da die Integrationsvariable ein WienerProzess ist. Zu Lösung des Problems hat Ito (wie auch Stratonovich, sie-he Günther und Jüngel [4]) ein anderes Integral eingeführt, das man überdWt definieren kann. Da das Integral in dieser Arbeit nicht von allzu großerBedeutung ist, sei für eine genaue Definition des Integrals auf Øksendal [19]verwiesen.

26


Euler-Maruyama-Verfahren

Um eine Ito-SDE zu berechnen, wird häufig auf eine Diskretisierung der Dif-ferentialgleichung zurück gegriffen. Die wohl gängigste und einfachste hierbeiist das Euler-Maruyama-Verfahren. Sei also Xt ein stochastischer Prozess,der die Gleichung 2.10 erfüllt, dann kann man die Lösung durch X(t) wiefolgt approximieren: Sei der Startwert bekannt, also X0 = X(0), und sei-ne Zi standardnormalverteilte Zufallszahlen, dann berechnet man die X(ti)iterativ mit der Vorschrift

X(ti+1) = X(ti) + a(X(ti))(ti+1 − ti) + b(X(ti))√ti+1 − tiZi+1. (2.12)

In Algorithmus 2 wird ein solches Verfahren bereits genutzt, da wir in Glei-chung 2.4 eine einfache Euler-Maruyama-Diskretisierung angewandt haben.Betrachtet man die triviale Differentialgleichung dS(t) = dWt, mit S(0) = 0und äquidistanter Schrittweite ∆t, folgt mit a = 0 und b = 1 in obigerSchreibweise: X(ti+1) = X(ti) +

√∆tZi+1. Es ergibt sich die Gleichung 2.5.

Um die Geometrische Brownsche Bewegung herzuleiten, benötigen wir dasIto-Lemma:Lemma 2.9 (Ito ). Sei Xt ein Ito-Prozess und h : R×R→ R zweimal stetigdifferenzierbar. Dann ist Yt = h(t,X(t)) = h(t, x) ein Ito-Prozess und es giltfolgende Formel (Ito-Formel):

dY (t) =(∂h

∂t(t,X(t)) + ∂h

∂x(t,X(t))a(t,X(t)) + ∂2h

2∂x2 (t,X(t))b(t,X(t))2)dt

+ ∂h

∂x(t,X(t))b(t,X(t))dWt (2.13)

Beweis. Ein Beweis des Lemmas findet sich in den Werken Øksendals [19]und Korns [10].Wenden wir nun das Ito-Lemma auf die Ito-SDE

dS(t)S(t) = µdt+ σdWt (2.14)

an. Diese Differentialgleichung hat Drift a(t, S(t)) = µS(t) und Diffusionb(t, S(t)) = σS(t). Ein Problem, das bei der „normalen“ Brownschen Bewe-gung noch besteht, ist die Negativität. Deshalb logarithmieren wir die SDE,also sei h(t,X(t)) = h(x) = ln(x) mit X(t) = S(t). Mit dem Ito-Lemmafolgt

dYt =(

1xa(t, S(t))− 1

2b(t, S(t))2

S(t)2

)dt+

(1xb(t, S(t)

)dWt

=(µS(t)S(t) −

12σ2S(t)2

S(t)2

)dt+

(1

S(t)σS(t))dWt

=(µ− 1

2σ2)dt+ σdWt

27


Das heißt, es gilt

d ln(S(t)) =(µ− 1

2σ2)dt+ σdWt (2.15)

und wir können S(0) als Startwert definieren, die Integralschreibweise wählenund exponieren. Daraus ergibt sich folgende Definition.

Definition 2.10 (geometrische Brownsche Bewegung). Sei Wt ein WienerProzess, dann ist die geometrische Brownsche Bewegung durch

St = S0 exp((µ− 1

2σ2)t+ σWt

)(2.16)

gegeben.

Input : ∆t, µ, σ,S0, Zi ∼ N (0, 1) für i = 1, . . . , nset S(0) = S0for i = 1, . . . , n do

S(ti) = S(ti−1) + µS(ti−1)∆t+ σS(ti−1)√

∆tZiendOutput : Vektor S ∈ Rn+1, der den Kursverlauf beinhaltetAlgorithmus 3 : Simulation eines Kurses auf Grundlage der Geometri-schen Brownschen Bewegung

Wir haben bereits festgestellt, dass in Gleichung 2.14 µS(t) den Drift dar-stellt, dass heißt für ein positives µ ist ein Anstieg des Kurses zu erwarten,zumindest langfristig im Mittel gesehen, für ein negatives µ gerade ein Fallendes Kurses. Im Gegensatz dazu steht logischerweise das σ für den stochas-tischen Anteil, für die Zufälligkeit. Wenn wir den speziellen Fall betrachten,bei dem σ = 0 ist, fällt in Gleichung 2.14 der stochastische Anteil weg undes bleibt eine deterministische Differentialgleichung. Abbildung 2.3 zeigt eineZerlegung der Geometrischen Brownschen Bewegung in Drift und Diffusion.Die Grafik ist so zu verstehen, dass wenn wir Drift- und Diffusionskurveaddieren und den Startwert subtrahieren, so ergibt es den Kursverlauf. Diegenaue Implementation ist in Matlabcode 8.1 dokumentiert. Zu bemerkenist noch, dass hier 10 Jahre simuliert werden, damit der Einfluss des Driftsdeutlich wird.Der Pfad einer Geometrischen Brownschen Bewegung wird in Algorithmus 3simuliert, dieser ist auch in Matlabcode 8.1 realisiert und erzeugt die Abbil-dung 2.4. In der Grafik sind drei Pfade eines Kurses mit Startwert S0 = 90simuliert. Man erkennt, dass die (mögliche) Negativität der Brownschen Be-wegung verschwunden ist, diese tritt auch nicht bei Startwerten nahe nullauf. Mit Startwert S0 = 0 ergibt sich logischerweise der triviale Kursverlauf,der keinen Wert außer 0 annimmt.

28


Abbildung 2.3: Zerlegung einer Geometrischen Brownschen Bewegung

Den Algorithmus 3 werden wir in Kapitel 3 verwenden, um mehrere Kurseeines Underlying zu simulieren. Es ist aber zu beachten, dass für die SDEaus Gleichung 2.14 die Euler-Maruyama-Diskretisierung aus Gleichung 2.12vorgenommen wurde. Das ist eine recht einfache Approximation, eine weitereMöglichkeit zur Diskretisierung werden wir im späteren Verlauf der Arbeitnoch kennenlernen.

29


Abbildung 2.4: Simulation von drei Geometrischen Brownschen Bewegungen

30

3 Monte-Carlo-Simulation

3.1 GrundideeZiel der Monte-Carlo-Simulation ist es, zu einer gegebenen Zufallsvariable Xeinen Erwartungswert E[X] zu schätzen. Dazu werden m Realisationen desZufallsexperimentes mit computergenerierten Zufallszahlen simuliert. Grund-lage der Monte-Carlo-Simulation ist der folgende Satz.

Theorem 3.1 (Zentraler Grenzwertsatz der Statistik). Seien X1, X2, . . . , Xm

unabhängige und identisch verteilte Zufallszahlen mit E[Xi] = µ und 0 <var(Xi) = σ2 < ∞ für alle Indizes i = 1, · · · ,m , dann gilt für alle reellenZahlen x :

P

(m∑i=1

Xi − µσ√m≤ x

)m→∞−−−→ Φ(x) (3.1)

, wobei Φ(x) die Verteilungsfunktion der Standard-Normalverteilung bezeich-net (vgl. Definition 2.3).

Beweis. Für eine Beweis sei auf Hesse [6, Kapitel 7, Theorem 7.2.1] verwie-sen.

Bemerkung 3.2. Aus Theorem 3.1 folgt direkt, dass Xm := ∑mi=1Xi für

große m näherungsweise normalverteilt ist, Xm ∼ N (mµ,mσ2) mit E[Xm] =nµ und var(Xm) = nσ2. Außerdem stellt Theorem 3.1 keinerlei Anforderungan die Verteilung der unabhängigen Zufallsvariablen Xi, insbesondere müs-sen diese Xi nicht normalverteilt sein.

Stellen wir nun, wie in dem einleitenden Beispiel dieser Arbeit, den Zu-sammenhang zwischen dem Erwartungswert eines Zufallsexperiments undder Flächenberechnung her. Sei dazu f die Dichte der Zufallsvariable X,die auf dem Einheitsintervall definiert sei, und der Erwartungswert gegebendurch

E[X] =∫ 1

0f(x)dx.

Wenn wir das Integral als Flächeninhalt auffassen, können wir die Fläche un-ter der Funktion f(x) durch Summen mit gegebenen Stützstellen U1, . . . Umapproximieren. ∫ 1

0f(x)dx ≈ 1

m

m∑i=1

f(Ui)

31

KAPITEL 3. MONTE-CARLO-SIMULATION

Nutzen wir dazu als Stützstellen Ui gleichverteilte Zufallszahlen auf demEinheitsintervall, so folgt mit dem starken Gesetz der großen Zahlen, dassfür m→∞ die Konvergenz

1m

m∑i=1

f(Ui)→∫ 1

0f(x)dx

mit Wahrscheinlichkeit 1 gilt. Darauf aufbauend wird im Folgenden derMonte-Carlo-Schätzer definiert.

Definition 3.3 (Monte-Carlo-Schätzer). Sei X eine Zufallsvariable und sei-en X1, . . . , Xm Stichproben von X, so ist der (einfache) Monte-Carlo-Schätzergegeben durch

Xm := 1m

m∑i=1

Xi (3.2)

Es folgt sehr schnell aus der Definition und Kolmogorows starkem Gesetzder großen Zahlen, dass der Monte-Carlo-Schätzer erwartungstreu ist, dasheißt es gilt E[Xm] = E[X]. Wichtig für eine Bewertung des Schätzers isteine Aussage der folgenden Form: Mit welcher Wahrscheinlichkeit liegt dasErgebnis des Schätzers in der „Nähe“ des exakten Werts? Solche Aussagenkönnen durch Konfidenzintervalle mathematisch dargestellt werden werden.

Definition 3.4 (p-Konfidenzintervall). Sei ε eine reelle Zahl, p eine reelleZahl aus dem Einheitsintervall, m eine natürliche Zahl und µ = E[Xi], dannnennen wir ein Intervall I der Form I = [µ−ε, µ+ε] ein p-Konfidenzintervallfür die Monte-Carlo-Simulation, wenn für die Wahrscheinlichkeit

P

(1m

m∑i=1

Xi ∈ I)

= p (3.3)

gilt.

Die Konfidenzintervalle hängen, wie wir im nächsten Satz zeigen werden,von den Parametern des Schätzers ab: von der Simulationshäufigkeit m, derStandardabweichung σ und dem Erwartungswert µ der Zufallsvariable.

Satz 3.5. Sei ein p ∈ (0, 1) gegeben. Dann gibt es ein k > 0, sodass dasp-Konfidenzintervall Im folgender Form ist:

Im =[µ− kσ√

m,µ+ kσ√

m

](3.4)

Beweis. In diesem Beweis nehmen wir an, dassm hinreichend groß ist, sodassX approximativ normalverteilt ist. Zuerst berechnen wir die Varianz des

32


Monte-Carlo-Schätzers:

var(Xm) = var

(1m

m∑i=1

Xi

)

= 1m2

m∑i=1

var(Xi)

= σ2

m

Für die Standardabweichung des Monte-Carlo-Schätzers gilt also:√var(Xm) = σ√

m

Transformieren wir nun Xm zu einer standardnormalverteilten Zufallszahl,wie in Gleichung 2.2, das heißt:

Y = Xm − µσ√m

∼ N (0, 1) (3.5)

Geben wir nun den Wert der Wahrscheinlichkeit p an und formen die Glei-chung um:

p = P (−k ≤ Y ≤ k)

= P

−k ≤ Xm − µσ√m

≤ k

= P

(µ− kσ√

m≤ Xm ≤ µ+ kσ√

m

)

Für k gilt dann mit den Eigenschaften der Standardnormalverteilung: MitP (−k ≤ Y ≤ k) = p folgt, dass P (Y ≤ k) = p

2 + 12 , was genau dann der Fall

ist, wenn für k gilt:k = Φ−1

(p

2 + 12

)(3.6)

Häufig ist uns aber die Varianz einer Zufallsvariable unbekannt, was dieNotwendigkeit mit sich bringt, die Varianz zu schätzen. Die geschieht in derRegel mit der sogenannten korrigierten Stichprobenvarianz.

Definition 3.6 (korrigierte Stichprobenvarianz). Sei X eine Zufallsvariablemit E[X] = µ und σ2 = var(X). Weiter seien X1, . . . , Xm unabhängige und

33


identisch verteilte Zufallszahlen mit Xi ∼ X sowie dazu Xm der Monte-Carlo-Schätzer wie in Definition 3.3.Die korrigierte Stichprobenvarianz s2 ist

s2 := 1m− 1

m∑i=1

(Xi − Xm)2

Bemerkung 3.7 (korrigierte Stichprobenvarianz und t-Verteilung). Die kor-rigierte Stichprobenvarianz ist ein erwartungstreuer Schätzer von σ2 = var(X).Da, wie oben bereits angemerkt, häufig auf die korrigierte Stichprobenvarianzzurück gegriffen wird und nicht die tatsächliche Varianz genutzt wird, ist essinnvoll ein Konfidenzintervall für den Fall mit s2 anzugeben. Wichtig ist,dass die standardisierte Zufallsvariable wie in den Gleichungen 2.2 bzw. 3.5

Z = X − µs

√m

nicht standard-normalverteilt ist, sondern nach der sogenannten student-schen t-Verteilung mit m− 1 Freiheitsgraden verteilt ist. Die Herleitung derKonfidenzintervalle erfolgt analog, bis auf die Tatsache, dass man das k inGleichung 3.6 durch die Inverse der t-Verteilung berechnet. Grundsätzlichist die Standardnormalverteilung aber auch an dieser Stelle anwendbar, dadie studentsche t-Verteilung gegen die Standardnormalverteilung für m→∞konvergiert, wie Mosler und Schmid [17, 4.3.2] gezeigt haben. Resultat dar-aus ist letztendlich, dass man für geeignet große m die Konfidenzintervallenach Satz 3.5 berechnen kann, auch wenn man „nur“ die korrigierte Stich-probenvarianz zur Verfügung hat.

An der Form der Konfidenzintervalle erkennen wir nun aber zwei wich-tige Eigenschaften der Monte-Carlo-Simulation. Erstens sehen wir, dass imZähler der Konfidenzintervallgrenzen die Standardabweichung der Zufalls-variable X steht und somit die Konvergenz der Simulation beeinflusst. Jekleiner die Varianz und damit die Standardabweichung, desto schneller kon-vergiert die Methode. Deshalb erscheint es nur folgerichtig zur Verbesserungder Monte-Carlo-Methode bei der Varianz anzusetzen und zu versuchen diesezu reduzieren. Genaueres sehen wir in Kapitel 4, in welchem wir mehrere Me-thoden zur Varianzreduktion kennenlernen und anwenden werden.Beispielemit 95%-Konfidenzintervallen sehen wir in Kapitel 4 über Varianzreduktion.Zweitens sehen wir aber auch, dass die Konvergenz von der Ordnung O( 1√

m)

ist, was die sehr langsame Konvergenz der Standardmethode erklärt.

3.2 OptionspreisbewertungAusgehend von der Beschreibung einer Option aus Kapitel 2 und den Er-kenntnissen, die wir in Abschnitt 3.1 bereits erworben haben, können wir

34


nun Optionspreise mit dem Monte-Carlo-Schätzer approximieren. Dazu ver-gegenwärtigen wir uns die Notationen aus Kapitel 2. Sei S(t) also der Kursdes Underlying Assets zum Zeitpunkt t mit 0 ≤ t ≤ T , X der Strike derOption und T der Fälligkeitszeitpunkt. Betrachten wir hier eine europäischeCalloption. Wir können zwei Fälle unterscheiden:

(i) S(T ) > X, das bedeutet, dass eine Ausübung der Option sinnvoll wäre,denn das Underlying kann zu einem geringeren Preis gekauft werdenals der Preis, zu dem es auf dem Markt gehandelt wird.

(ii) S(T ) ≤ X. In diesem Fall wird der (rational handelnde) Holder dieOption verfallen lassen.

Für Putoptionen gelten diese Aussagen für die negierten Ungleichungenund wenn wir mit europäischen Optionen rechnen, genügt es den ZeitpunktT zu betrachten, da die Option nur in T ausgeübt werden kann. ObigerLogik folgend, können wir für den Payoff der Option zum Zeitpunkt T diegeschlossene Formel

(S(T )−X)+ = max0, S(T )−X (3.7)

angeben. Damit wir Aussagen über den jetzigen Wert der Option machenkönnen, ist es notwendig den Present Value des Payoffs, also den abgezins-ten Wert, zu betrachten. Und da wir keine sichere Aussage darüber machenkönnen, wie der Kurs des Underlying in T ist, müssen wir den erwartetenBarwert, den Erwartungswert des Present Value

E[e−rT (S(T )−X)+] (3.8)

berechnen. Hierbei bezeichnet r den risikolosen Zins wie in Kapitel 2, welcheim Folgenden die erwartete Rendite µ des Underlying ersetzt. Das beruht aufdem Prinzip der risikoneutralen Bewertung (vgl. Glasserman [3, Abschnitt1.2.2]), welches wir in dem Abschnitt 2.2 über die Marktannahmen voraus-gesetzt haben.In Abschnitt 2.5 haben wir den Ansatz zur Beschreibung eines Kursverlaufseines Assets kennengelernt. Gehen wir also davon aus, dass der Kurs S(T )des Underlayings der stochastischen Differentialgleichung

dS(t)S(T ) = rdt+ σdWt (3.9)

genügt (wie in Gleichung 2.14).Wt bezeichnet hier einen Wiener Prozess mitErwartungswert 0 und Varianz t, r den risikolosen Zins und σ die Volatili-tät des Underlyings. Die Lösung der Stochastischen Differentialgleichung istgegeben durch

S(T ) = S(0)exp(

[r − 12σ

2]T + σWT

)

35


(vgl. Glasserman [3, Kapitel 1]), wobei WT normalverteilt mit Erwartungs-wert 0 und Varianz T ist. Mit einer standardnormalverteilten Zufallsvaria-blen Z ∼ N (0, 1) gilt, wie wir in Gleichung 3.5 gezeigt haben, dass

√TZ

verteilt ist wie WT . Vergleiche hierzu auch Seydel [21, S.32]. Damit ergibtsich für den Kurs des Assets zum Zeitpunkt T

S(T ) = S(0)exp(

[r − 12σ

2]T + σ√TZ

)(3.10)

. Dieser Gleichung liegt zugrunde, dass wir zu gegebener Simulation desWiener Prozesses eine Lösung der SDE aus Gleichung 3.9 genau bestimmenkönnen. In manchen Fällen ist es aber so, dass wir nicht nur an dem Endwertdes Kurses interessiert sind, sondern dass wir auch zu gegebenen Zeitschrit-ten den Kurs des Assets benötigen, um den Optionswert angeben zu können.Das kommt zum Beispiel bei Asiatischen Optionen vor, die auf dem geometri-schen oder arithmetischen Mittel des Kurses beruhen. In diesem Fall müssenwir mit einem Kursverlauf rechnen, der zu jedem diskreten Zeitpunkt einenWert angibt. Darauf kommen wir am Ende dieses Kapitels zurück, indemwir einen Algorithmus mit der Euler-Maruyama-Diskretisierung angeben.

Was haben wir in diesem Abschnitt bereits herausgefunden?

1. Wir haben mit Gleichung 3.8 eine Darstellung abhängig vom Kurs desUnderlying für die erwartete Auszahlung.

2. Wir haben mit Gleichung 3.10 eine Formel für den (zufälligen) Kurs desUnderlying abhängig von einer standardnormalverteilten Zufallszahl.

Bringen wir Gleichung 3.8 und 3.10 zunächst zusammen. Dann haben wir,unter der Annahme, wir könnten den Erwartungswert berechnen und diestandardnormalverteilte Zufallszahl erzeugen, alle Werkzeuge die erwarteteAuszahlung C einer Option zu berechnen.

C = E[e−rT

((S(0)exp

([r − 1

2σ2]T + σ

√TZ

))−X

)+](3.11)

Tatsächlich haben wir gerade für die Erwartungswertberechnung den Monte-Carlo-Schätzer eingeführt. Sei nun der Monte-Carlo-Schätzer für die erwar-tete, abgezinste Auszahlung der Option Cm, Zi seien unabhängige, standard-normalverteilte Zufallszahlen und m die Anzahl der Wiener Prozesse. Danngilt

Cm := e−rT[

1m

m∑i=1

((S(0)exp

([r − 1

2σ2]T + σ

√TZi

))−X

)+](3.12)

36


Mithilfe der kompakteren Schreibweise

Cm = e−rT[

1m

m∑i=1

Si(T )−X)+]

(3.13)

mit Si(T ) = S(0)exp([r − 1

2σ2]T + σ

√TZi

)sehen wir direkt, dass man

(Si(T )−X)+ als Stichprobe einer Zufallsvariablen interpretieren kann. Dasrechtfertigt, dass wir in Gleichung 3.12 von einem Monte-Carlo-Schätzersprechen. Daraus folgen nun die in Abschnitt 3.1 hergeleiteten Eigenschaf-ten eines Monte-Carlo-Schätzers, insbesondere, dass Cm → C fast sicherkonvergiert, für m→∞. Zur numerischen Umsetzung dieses Verfahrens giltAlgorithmus 4, wie ihn Glasserman [3, Kapitel 1.1.2] hergeleitet hat.

Input : T , X, S(0), r, σfor i=1, · · · ,m do

generate Zi ∼ N (0, 1)set Si(T ) = S(0)exp

([r − 1

2σ2]T + σ

√TZi

)set Ci = e−rT (Si(T )−X)+

endset Cm = 1

m

∑mi=1 Ci

Result : approximativer Optionspreis Cm

Algorithmus 4 : Grundalgorithmus zur Approximation des Optionspreises

Durch Algorithmus 4 haben wir nun einen Grundalgorithmus hergeleitet,mit dem wir den Wert einer Europäischen Calloption berechnen können, fürein Put müsste Si(T )−X in X−Si(T ) geändert werden. An dieser Stelle istanzumerken, dass wir in Abschnitt 2.5 mit dem Euler-Maruyama-Verfahreneine Approximation des Kursverlaufs S(ti) eingeführt haben. Diesen könnenwir nutzen, um pfadabhängige Optionen zu bewerten. Es liegt also nahe,den Kurs, der in dem dazugehörigen Algorithmus 3 berechnet wird, in Algo-rithmus 4 für S(ti) einzusetzen. Sei dazu Λ(S(t)) die Auszahlungsfunktioneiner beliebigen, pfadabhängigen Option, dann ist durch Algorithmus 5 derAlgorithmus zur Optionsbewertung mit dem Euler-Maruyama-Verfahren ge-geben. Einen Algorithmus für Asiatische Optionen geben wir in Kapitel 6an, in welchem wir den Wert einer Asiatischen Option berechnen.Noch zu erläutern ist, wie die Erzeugung der standardnormalverteilten

Zufallszahlen Zi effizient zu realisieren ist. Ein kurzer Überblick dazu wirduns im nächsten Abschnitt erwarten.Natürlich gelten die in Abschnitt 3.1 angesprochenen Probleme auch fürunseren Algorithmus 4, sodass wir nicht umhin kommen, die Konvergenzdes Verfahrens zu verbessern. Ansatz dazu wird die Varianz des Schätzerssein. Dazu in Kapitel 4 mehr.

37


Input : T , X, S(0), r, σfor i=1, · · · ,m do

for j = 1, · · · ,n dogenerate Zi,j ∼ N (0, 1)set Si(tj) = Si(tj−1) + µSi(tj−1)∆t+ σSi(tj−1)

√∆tZi,j

endset Ci = e−rTΛ(Si(t))

endset Cm = 1

m

∑mi=1Ci

Result : approximativer Optionspreis Cm

Algorithmus 5 : Approximation eines pfadabhängigen Optionspreises

3.3 PseudozufallszahlenRandom numbers should not begenerated with a method chosen atrandom.

Donald E. Knuth 1

Im vorherigen Abschnitt haben wir gesehen, dass Zufallszahlen eine zen-trale Rolle in der Monte-Carlo-Simulation einnehmen. Deswegen erscheint essinnvoll, sich Gedanken darüber zu machen, wie solche Zufallszahlen erzeugtwerden und ob wir tatsächlich identisch verteilte und unabhängige Zufalls-zahlen erzeugen können. Zufallszahlen haben in der Mathematik und in derInformatik sehr viele Anwendungen, weswegen die Generierung von „guten“Zufallszahlen durchaus noch Gegenstand der aktuellen Forschung ist, wie derMersenne-Twister-Algorithmus 2 aus dem Jahr 1997 eindrucksvoll zeigt.

Dieser Abschnitt beruht vor allem auf dem Standardwerk „Art of ComputerProgramming“ von Knuth [9, Kapitel 3], sowie „Tools for ComputationalFinance“ [21, Kapitel 2] von Seydel, das sich speziell auf die finanzmathe-matische Anwendung bezieht.

Für den Monte-Carlo-Algorithmus zur Optionsbewertung benötigen wir stan-dardnormalverteilte Zufallszahlen X ∼ N (0, 1), um die Pfade der WienerProzesse zu simulieren. Dabei ist die Vorgehensweise, um standardnormal-verteilte Zufallszahlen zu generieren, wie folgt:

1siehe [9, Vol. II, Seminumerical Algorithms, Kapitel 3.1]2Die offizielle Seite des Mersenne-Twister-Algorithmus: http://www.math.sci.hiroshima-u.ac.jp/∼m-mat/MT/emt.html

38


(1) Erzeuge gleichverteilte (Pseudo-)Zufallszahlen Y ∼ U [0, 1](2) Transformiere die gleichverteilten Zufallszahlen durch eine geeignete Funk-

tion f(Y ) so, dass die Zufallsgröße X := f(Y ) standardnormalverteiltist, also X ∼ N(0, 1)

Wir sprechen an dieser Stelle von Pseudozufallszahlen, da Zufallszahlen, diemit einem deterministischen Algorithmus erzeugt worden sind, streng ge-nommen nicht unabhängig sind. Sie weisen nur eine sehr geringe Korrelationauf und sind berechenbar. Ziel der folgenden Pseudozufallszahlengeneratorenzur Erzeugung gleichverteilter Zufallszahlen ist demnach die näherungsweiseUnabhängigkeit der Zufallszahlen.

3.3.1 Linearer KongruenzgeneratorDer lineare Kongruenzgenerator nach Lehmer [12, s.141-146] wurde bereits1949 vorgestellt und ist einer der ersten Zufallszahlengeneratoren. Deswegenund weil er die Methodik vieler Zufallszahlengeneratoren verdeutlicht, stellenwir diesen Zufallszahlengenerator hier vor.

Definition 3.8 (Linearer Kongruenzgenerator). Sei X eine endliche Menge,f : X → X eine Funktion, der Modulus m eine natürliche Zahl und reelleZahlen a, b, x0 ∈ 0, . . . ,m− 1, für die gilt:

f(xn) = xn+1 = (axn + b) mod m (3.14)

Dann heißt f linearer Kongruenzgenerator. Für den Fall, dass b = 0 ist,sprechen wir von einem multiplikativen Kongruenzgenerator.

Input : Parameter a, b,m, x0 wie in der Definition, Iterationshäufigkeitk m

for i = 0, . . . , k doset xi+1 = (axi + b) mod mset yi = xi

m

endResult : näherungsweise gleichverteilte Pseudozufallszahlen yi

Algorithmus 6 : Linearer Kongruenzgenerator

Aus der Iterationsvorschrift ergibt sich der Algorithmus 6 und die yi sindauf dem Einheitsintervall gleichverteilte Pseudozufallszahlen. Entscheidendfür die Güte der Zufallszahlen ist hier vor allem die Wahl der Parameter a undm, denn es können maximal m verschiedene Zufallszahlen auftreten. UnterUmständen kann eine Zufallszahl xj zum zweiten Mal in der Iteration erzeugtwerden, dann wiederholt sich die Periode und es wurden nur j verschiedeneZufallszahlen generiert, mit j < k. Die Problematik dessen wird an dem

39


folgenden Zahlenbeispiel deutlich: Seien a = 2, m = 10, x0 = 3, f(xn) wiein Gleichung 3.14 und yi = xi

m, dann folgt

y1 = 0, 6 y2 = 0, 2 y3 = 0, 4 y4 = 0, 8

und für alle Indizes i ∈ 0, 1, . . . , k gilt yi+4 = yi. Das bedeutet, dass nurfünf voneinander verschiedene Zufallszahlen erzeugt wurden. Zur Lösung die-ses Problems bringen wir folgenden Satz von Knuth [9, Theorem A, Kapitel3.2] an, der mit Erkenntnissen der Zahlentheorie bewiesen wurde. Deshalbverzichten wir auf diesen Beweis.

Satz 3.9 (Knuth). Ein Linearer Kongruenzgenerator wie in Gleichung 3.14mit m > 1 liefert für beliebigen Startwert X0 ∈ 0, 1, . . . ,m − 1 genaudann maximale Periodenlänge von m, wenn die folgenden Bedingungen er-füllt sind:(1) b und m sind relativ prim(2) jede Primzahl, die Teiler von m ist, teilt auch a− 1(3) wenn 4 Teiler von m ist, dann ist auch a− 1 durch 4 teilbar

Dieser Satz hilft uns nun dabei, bessere Parameter für einen Kongruenzge-nerator zu finden. Für eine detailliertere Darstellung sei an dieser Stelle aufKnuth [9, Kapitel 3.2] verwiesen. Der Modulusm kann z.B. durchm = 231−1festgelegt sein, da dies die größte darstellbare 32-Bit Zahl ist und gleichtzeitigdie Bedingung einer Mersenne-Primzahl erfüllt. Dieser Logik folgend steht inMatlab ein multiplikativer Kongruenzgenerator mit den Parametern a = 75

und m = 231 − 1 zur Verfügung, wie Lewis, Goodman und Miller [13] 1969empfohlen haben. Diese Parameterwahl gewährleistet zwar nach obigem Satzmaximale Periodenlänge, jedoch ist diese Länge für viele Anwendungen nochzu gering. Ein weiteres Problem bei Linearen Kongruenzgeneratoren ist, dassdie erzeugten Pseudozufallszahlen beziehungsweise die d-Tupel von d aufein-anderfolgenden Zufallszahlen (y1, y2 . . . yd), (yd+1, yd+2 . . . y2d) , . . . ∈ [0, 1]dauf Hyperebenen liegen. Marsaglia spricht von „’crystalline’ nature of multi-plicative generators“ 3. Das entspricht natürlich nicht der Anforderung voneiner Gleichverteilung und Unabhängigkeit der Zufallszahlen. Hierbei ist zubeachten, dass eine höhere Anzahl an Hyperebenen, eine „bessere“ Vertei-lung der Zufallszahlen gewährleistet. Marsaglia [14] zeigte, dass bei gegebe-nem Modulus m und gegebener Dimension d maximal (d!m) 1

d Hyperebenenim d-dimensionalem Einheitswürfel existieren, auf denen die d-Tupel liegen.Diese Berechnungen beziehen sich auf multiplikative Kongruenzgeneratoren,für Generatoren mit b 6= 0 gelten aber analoge Aussagen. Für eine weiter-gehende Analyse des Hyperebenenverhaltens von Kongruenzgeneratoren seiauf Seydel [21, Kapitel 2.1.3] verwiesen. Wir sehen im folgenden Beispiel,dass die d-dimensionalen Vektoren durchaus auf sehr wenigen Hyperebenen

3 Marsaglia [14, Seite 1, Zeile 15]

40


liegen können und die obere Schranke, die Marsaglia angibt, deutlich unter-schritten wird.

Abbildung 3.1: 3-dimensionale Zufallsvektoren mit RANDU I

Beispiel 3.10 (RANDU). Betrachten wir für einen multiplikativen Kon-gruenzgenerator (wie in Gleichung 3.14 defniert) den Fall a = 216 + 3 undm = 231. Dieser Generator heißt RANDU und wurde in den 1960er Jahrenvon IBM verbreitet. Die daraus resultierenden dreidimensionalen Zufallsvek-toren liegen aber, wie auch Seydel [21] gezeigt hat, auf nur 15 Hyperebenen,was diesen Generator für die meisten Anwendungen unbrauchbar macht. InAbbildung 3.1 sind 20000 Zufalllszahlen zu dreidimensionalen Vektoren zu-sammengefasst. Auf den ersten Blick erscheint die Verteilung der Punkte imdreidimensionalen Einheitsintervall annähernd gleichverteilt. Rotieren wiraber die Grafik, was in Abbildung 3.2 dokumentiert ist, lassen sich die (nur)15 Ebenen deutlich erkennen.

41


Abbildung 3.2: 3-dimensionale Zufallsvektoren mit RANDU II

Die Abbildungen sind mit den Matlabcodes 8.3 und 8.4 erzeugt worden.Dabei ist zu beachten, dass die Abbildungen, die Matlab erzeugt, von Handrotiert werden müssen, um das Hyperebenenverhalten erkennbar werden zulassen.Auch wenn die Periodizität der linearen Kongruenzgeneratoren Probleme

bereitet, sind sie dennoch Grundlage einiger weiterentwickelter Generato-ren, die wir aber nicht näher behandeln werden. So haben unter anderemWichmann und Hill [23] kombinierte Generatoren zur Verbesserung des vor-gestellten Algorithmus nach Lehmer [12] vorgestellt, die auf MultiplikativeKongruenzgeneratoren zurückgreifen. Ziel dabei war die Hyperflächenproble-matik der Zufallszahlen zu verringern, sowie die Periodenlänge zu erhöhen,ohne dabei jedoch den Eingabewert des Modulus m zu erhöhen. Letzteresführt bei der Implementation unter anderem hinsichtlich der Darstellbarkeitder Zahlen zu Problemen. Für eine tiefergehende Diskussion dieser Genera-toren kann auf Wichmann und Hill [23] oder Glasserman [3] zurückgegriffenwerden.

3.3.2 Weitere GeneratorenNeben den Linearen Kongruenzgeneratoren gibt es einige weitere, sehr einfa-che Methoden, die relativ gute Ergebnisse hervorbringen, so zum Beispiel dersogenannte verzögerte Fibonacci-Generator. Dieser ist durch die Vorschrift

Ni+1 = (Ni−α −Ni−β) mod m

definiert. Hierbei muss die natürliche Zahl i größer als α und β sein. Das be-deutet, dass die Werte N1, . . . , Nmaxα,β durch einen anderen Zufallszahlen-generator, zum Beispiel durch einen Linearen Kongruenzgenerator, erzeugt

42


werden müssen. Dann sind die Pseudozufallszahlen Ni+1m

annähernd gleich-verteilt auf dem Einheitsintervall.

Mersenne-TwisterDer aktuell wichtigste Generator ist der von Matsumoto und Nishimura[16] entwickelte Mersenne-Twister-Algorithmus, der Elemente der Fibonacci-Generatoren und der Linearen Kongruenzgeneratoren aufnimmt. Ziel war es,die Nachteile der vorhergehenden Generatoren auszumerzen, wie die zu ge-ringe Periodizität und die Gitterstruktur der Zufallsvektoren. Das ist derGrund, dass in vielen mathematischen Programmen der Mersenne-Twister-Algorithmus als Standardgenerator implementiert ist. Matsumoto und Nishi-mura [16] haben eine Periodenlänge von 219937−1 (eine Mersenne-Primzahl)erreicht, konnten aber trotzdem eine sehr schnelle Berechnungszeit und dieGleichverteilung bis in Dimension 623 gewährleisten. Letzteres ist in demSinne, wie wir es auch in Beispiel 3.10 (RANDU) diskutiert haben, zu ver-stehen.

3.4 Standardnormalverteilte ZufallszahlenNachdem wir gesehen haben, auf welche Art und Weise (näherungsweise)gleichverteilte Zufallszahlen generiert werden können, ist die Frage, wie wirvon gleichverteilten zu normalverteilten Zufallszahlen kommen. Dazu bietensich grundsätzlich zwei Vorgehensweisen an:

(i) Invertierung der Verteilungsfunktion

(ii) Transformation der Zufallszahlen

Grundlage der ersten Idee ist folgender Satz.

Satz 3.11. Sei U eine auf dem Einheitsintervall gleichverteilte Zufallsgrößeund F eine stetige, streng monoton wachsende Verteilungsfunktion. Dann istdie, durch

Z := F−1(U)

gegebene Zufallsgröße Z, nach F verteilt.

Beweis. Es sei zuerst angemerkt, dass die Inverse F−1 aufgrund der Mono-tonie existiert. Mit den Eigenschaften der Gleichverteilung folgt für alle zaus dem Einheitsintervall

P (U ≤ z) = z

43


, was man schreiben kann als

P (U ≤ F (z)) = F (z)

und das liefert uns direkt die Behauptung:

P (F−1(Z) ≤ z) = F (z)

An dieser Stelle fällt auf, dass für eine bestimmte Wahrscheinlichkeits-verteilung die Inverse der Verteilungsfunktion bekannt sein muss. Da wirin diesem Kapitel nicht an allgemeinen Verteilungen interessiert sind, son-dern explizit standardnormalverteilte Zufallszahlen erzeugen wollen, ist dasein Problem für uns. Für die Normalverteilung existiert nämlich keine ge-schlossene Formel der Verteilungsfunktion und ihrer Inversen. Aus diesemGrund sind wir hier auf numerische, das heißt approximative Verfahren zurInvertierung der Verteilungsfunktion angewiesen.

3.4.1 InvertierungZur Invertierung der Standardnormalverteilung kann unter anderem auf denBeasley-Springer-Moro Algorithmus zurückgegriffen werden. Dieser ist inMatlabcode 8.5 realisiert und erzeut sehr genaue Werte für die Inverse derStandardnormalverteilung. Hierbei wurde Glassermans Formulierung des Al-gorithmus [3, s.67/68] genutzt. Für Abbildung 3.3 sind 200000 gleichverteilteZufallszahlen dem Algorithmus von Beasley-Springer-Moro übergeben unddie Häufigkeit des Auftretens der Werte als Balkendiagramm dargestellt. Wirerkennen sofort, dass die erzeugten Zahlen approximativ standardnormalver-teilt sind.

44


Abbildung 3.3: Zufallszahlen mit dem Beasley-Springer-Moro-Algorithmus

Um den zweiten Ansatz zu beleuchten, müssen wir uns der Wahrschein-lichkeitstheorie bedienen und einen einfachen Fall des Transformationssatzesbetrachten.

Satz 3.12. Sei U eine Zufallsvariable mit Dichtefunktion f auf der MengeA = x ∈ Rn : f(x) > 0 und sei g eine Transformation mit g : A → B =g(A). Ferner sei die Transformation g umkehrbar und ihre Inverse g−1 seistetig differenzierbar. Dann gibt die Abbildung

z 7→ f(g−1(z))|dg−1

dz(z)|

die Wahrscheinlichkeitsdichte der Zufallsvariable Z := g(U) für z ∈ B an.

Beweis. Für einen Beweis sei auf Devroye [1, Theorem 4.2] verwiesen.

Um diesen Satz auf unsere Problemstellung anwenden zu können, müssenwir für alle reellen Zahlen x aus dem Einheitsintervall f(x) = 1 setzen. DasEinheitsintervall definiert dann die Menge A und Satz 3.12 ergibt für den Fallim Eindimensionalen eine Differentialgleichung erster Ordnung. Wir sucheneine Funktion g(x), die die Gleichung

|dg−1

dz(z)| = 1√

2πexp(−z

2

2 )

erfüllt. Analytisch können wir diesen Fall nicht berechnen. Dafür aber imzweidimensionalen Fall. Das führt uns zum sogenannten Box-Muller-Verfahren.

45


3.4.2 Box-Muller-AlgorithmusBeziehen wir uns in diesem Abschnitt auf die Notation des Satzes 3.12 unddefinieren in diesem Sinne A als den zweidimensionalen Einheitswürfel [0, 1]2und setzen f(x) = 1 für alle x aus A. Man kann dann zeigen, worauf wir andieser Stelle verzichten und auf Seydel [21, Kapitel 2.3.1] verweisen wollen,dass die Transformation

g(x) = √−2ln(x1) cos(2πx2)√−2ln(x1) sin(2πx2)

, x =(x1x2

)∈ A (3.15)

den Anforderungen des Satzes 3.12 genügt und die Determinante der JacobiMatrix genau die Dichtefunktion der zweidimensionalen Standardnormalver-teilung ist. Wenn also U auf [0, 1]2 gleichverteilt ist, so ist g(U) standardnor-malverteilt.

3.4.3 Polar-AlgorithmusUm den Box-Muller-Algorithmus zu verbessern, insbesondere um die tri-gonometrischen Funktionen zu umgehen, entwickelten Marsaglia und Bray[15] den sogenannten Polar-Algorithmus. Dieser ist weniger aufwendig, da erdurch Polartransformation die trigonometrischen Funktionen vermeidet. Da-zu seien U1, U2 auf dem Einheitsintervall gleichverteilte Zufallszahlen, dannsind V1 := 2U1 − 1 und V2 = 2U2 − 1 auf dem Intervall [-1,1] gleichver-teilt. Für den Polar-Algorithmus werden nun alle Zufallszahlen verworfen,die, jeweils zu zweidimensionalen Vektoren zusammengefasst, nicht auf derEinheitskreisscheibe liegen. In Formeln ausgedrückt verwerfen wir diese Zu-fallszahlen V1, V2, für die die Ungleichung V := V 2

1 + V 22 ≥ 1 gilt. Die ak-

zeptierten Zahlen transformieren wir wieder auf das Einheitsquadrat mit derPolartransformation

T

(V1V2

)=(

V 21 + V 2

21

2πarctan(V2V1

)

)

. Nun können wir festhalten, dass, wenn T1 und T2 die erste und die zweiteKomponente der Transformation definiert, Z1 :=

√−2ln(T1) sowie Z2 :=√

−2ln(T2) standardnormalverteilt sind. Durch Umformungen können wirfolgende Form erreichen

Z1 = V1

√−2 ln(V )

V

und

Z2 = V2

√−2 ln(V )

V

46


Damit ergibt sich ein Algorithmus auf Grundlage der Idee des Box-Muller-Verfahrens, der aber deutlich weniger Aufwand benötigt. Auch wenn manmiteinbezieht, dass ein Teil der erzeugten Zufallszahlen verworfen wird unddeshalb weniger Zufallszahlen entstehen als bei dem Box-Muller-Verfahren.Eben jener Grund ist es aber auch, der dazu führt, dass häufig mit numeri-scher Invertierung der Standardnormalverteilung gearbeitet wird und nichtmit der Transformation, zum Beispiel in Form des Polar-Algorithmus. Beider numerischen Invertierung werden keine der bereits erzeugten, gleichver-teilten Zufallszahlen verworfen, was einen Effizienzvorteil darstellen kann.Das ist auch in der Abbildung 3.4 zu erkennen, hierzu wurden mit demPolaralgorithmus aus 200000 gleichverteilten Zufallszahlen standardnormal-verteilte Zufallszahlen erzeugt. Die Grafik ähnelt Abbildung 3.3 sehr stark,was durchaus wünschenswert ist. Aber in Abbildung 3.4 sind zum selben In-put deutlich weniger Zufallszahlen entstanden. Es wurden insgesamt 157000standardnormalverteilte Zufallszahlen erzeugt, im Gegensatz zum Beasley-Springer-Moro-Algorithmus, in dem es 200000 waren. Das liegt an der an-gesprochenen Verwerfung der „unerwünschten“ Zahlen. Die Implementationzur Abbildung 3.4 ist in Matlabcode 8.6 wiederzufinden.

Abbildung 3.4: Zufallszahlen mit dem Polar-Algorithmus

47

4 VarianzreduktionMonte-Carlo-Methoden weisen eine Konvergenz auf, die zum einen nicht all-zu schnell erscheint und zum anderen von der Varianz des Monte-Carlo-Schätzers abhängt. Wir haben darauf schon im Kapitel 3 in Satz 3.5 hin-gewiesen und wollen nun in diesem Kapitel zwei Methoden einführen, diezu einer geringeren Varianz des Schätzers und somit zu einer schnellerenKonvergenz führen. Techniken, die zu einer Reduktion der Varianz führenkönnen, haben Glasserman [3] sowie Müller-Gronbach et al. [18] genaueruntersucht. Das Kapitel beruht zwar in Grundzügen auf den Erkenntnissendieser Lehrwerke, dennoch vermag es diese Arbeit nicht die Thematik derVarianzreduktion in dem Umfang, wie es vor allem bei Glasserman [3] derFall ist, zu beschreiben. Deshalb werden wir nur zwei verbreitete Verfahrenbehandeln, die Methode der Antithetic Variates und die der Control Variates.Weitere Verfahren können in Glassermans Buch [3, Kapitel 4] nachgelesenwerden. Diese Methoden tragen hier wie auch in der (deutschen) Literaturder Einheitlichkeit halber englische Bezeichnungen.

4.1 Antithetic VariatesBetrachten wir eine standardnormalverteilte Zufallsvariable Z, da wir damiteine geometrische Brownsche Bewegung simulieren wollen. Wegen der Sym-metrie der Standardnormalverteilung um 0 betrachten wir zu Z auch −Z.Für eine zweimal differenzierbare Funktion f definieren wir dann den geradenAnteil der Funktion fg(Z) = 1

2(f(Z) + f(−Z)) sowie den ungeraden Anteilfug(Z) = 1

2(f(Z) − f(−Z)). Es gilt die Gleichheit f(Z) = fg(Z) + fug(Z).Wenden wir den Monte-Carlo-Schätzer z.B. zur numerischen Integration oderzur Erwartungswertberechnung auf den geraden Anteil der Funktion an,so nähert diese Approximation denselben Wert an, wie der Monte-Carlos-Schätzer auf die Funktion angewendet. Das liegt an der Symmetrie der Zu-fallsvariablen Z und dem Transformationssatz. Es gilt nämlich die IdentitätE[f(Z)] = E[fg(Z)]. Somit folgt der Monte-Carlo-Schätzer mit AntitheticVariates

Y := 1m

m∑i=1

fg(Zi) = 12m

m∑i=1

(f(Zi) + f(−Zi)) (4.1)

Aber ist durch die Einführung der Antithetic Variates tatsächlich eine Va-rianzreduktion erfolgt? Wir rechnen nach. Dazu machen wir uns klar, dass

49

KAPITEL 4. VARIANZREDUKTION

wir, um ein vergleichbares Ergebnis zu erzielen, doppelt so viele „Standard“Monte-Carlo-Simulationen machen können wie Simulationen mit AntitheticVariates. Das liegt an den in etwa doppelt so hohen Kosten bei der Berech-nung mit Antithetic Variates.Eine geringere Varianz haben wir also genau dann, wenn

V ar(Y ) < V ar

(1

2m

2m∑i=1

f(Zi))

(4.2)

gilt. Die linke Seite berechnet sich wie folgt

V ar(Y ) = V ar

(1

2m

m∑i=1

f(Zi) + f(−Zi))

= 14m2V ar

(m∑i=1

f(Zi) + f(−Zi))

= 14mV ar (f(Zi) + f(−Zi))

, wobei wir die lineare Transformationseigenschaft der Varianz benutzt ha-ben und, dass die Zufallsvariablen f(Zi) + f(−Zi) identisch verteilt undunabhängig sind. Für die rechte Seite ergibt sich:

V ar

(1

2m

2m∑i=1

f(Zi))

= 2m4m2V ar(f(Zi))

= 12mV ar(f(Zi))

Also ist zur Ungleichung 4.2

14mV ar(f(Zi) + f(−Zi)) <

12mV ar(f(Zi))

äquivalent. Mit einigen Umformungen und Satz 2.2 aus Kapitel 2 ergibt sichdann

V ar(f(Zi)) + V ar(f(−Zi)) + 2Cov(f(Zi), f(−Zi)) < 2V ar(f(Zi))

, was genauCov(f(Zi), f(−Zi)) < 0

ergibt. Das bedeutet, wir haben eine Varianzreduktion mittels AntitheticVariates, genau dann, wenn die Kovarianz der Zufallszahlen negativ ist.

50


4.2 Control VariatesFür die Varianzreduktion mittels Control Variates betrachten wir zuersteinen Zufallsvektor (X, Y ), wobei X und Y zweimal integrierbar sind. Au-ßerdem soll Y derart sein, dass wir den Erwartungswert E[Y ] leicht berech-nen können. Sei also E[X] der Erwartungswert, den wir berechnen wollenund sei b eine beliebige reelle Zahl, dann ist zur Control Variate Xb :=X − b(Y − E[Y ]) der Monte-Carlo-Schätzer wie folgt gegeben

Cm,b = 1m

m∑i=1

Xi,b

= 1m

m∑i=1

(Xi − b(Yi − E[Y ]))

= bE[Y ] + 1m

m∑i=1

(Xi − bYi) (4.3)

, wobei (Xi, Zi)) unabhängig und wie (X, Y ) verteilt sind. Dass der SchätzerCm,b gegen den einfachen Monte-Carlo-Schätzer 1

m

∑mi=1Xi konvergiert, ist

offensichtlich. Nun stellt sich die Frage, wie das b aussehen muss, damit wireine möglichst gute Varianzreduktion erhalten. Setzen wir dazu voraus, dassV ar(X) > 0 sowie V ar(Y ) > 0 und ρX,Y den Korrelationskoeffizienten vonX und Y angibt.

Satz 4.1. Für die minimale Varianz des modifizierten Schätzers aus Glei-chung 4.3 gilt:

minb∈R

V ar(Cm,b) = V ar(X(1− ρ2X,Y ))

= V ar(Cm,b∗)

und b∗ = Cov(X, Y )V ar(Y ) . (4.4)

Beweis. Für einen Beweis sei auf Müller-Gronbach et al. [18, Satz 5.8] ver-wiesen.

Nun sagt uns der Satz 4.1 neben dem optimalen Wert für b auch, dass wirmöglichst ein Y wählen sollten, das stark mit X korreliert. Das Vorzeichendes Korrelationskoeffizienten ist hier egal. Ein Problem bei der Anwendungdes optimalen b aus Gleichung 4.4 ist, dass wir meist, wenn der Erwartungs-wert von X unbekannt ist, die Kovarianz von X und Y nicht berechnenkönnen. Deshalb nutzt man häufig anstelle der Kovarianz die korrigierteStichprobenkovarianz, wie auch die korrigierte Stichprobenvarianz für die

51


Varianz von Y . Das heißt unser optimales b∗ schreiben wir näherungsweiseals

bm =∑mi=1(Yi − Y )(Xi − X)∑m

i=1(Yi − Y )2(4.5)

und wir ersetzen b in Gleichung 4.3 durch bm, wobei hier X sowie Y jeweilsdas arithmetische Mittel der Stichprobe bezeichnet. Glasserman [3, s.187]sowie Müller-Gronbach et al. [18, Satz 5.11] haben gezeigt, dass bm gegen b∗mit Wahrscheinlichkeit 1 für m→∞ konvergiert.

Beispiel 4.2. Wenn wir Optionsbewertung mit einer varianzreduzierten Monte-Carlo-Simulation durchführen wollen, so können wir, wie es Glasserman [3,s.188/189] vorgeschlagen hat, neben dem Payoff der Option den Kurs desUnderlying als Y in dem modifizierten Monte-Carlo-Schätzer nutzen. Vor-teil daran ist, dass wir die Kurse zu allen Zeitpunkten und zu allen Pfa-den bereits berechnet haben. Außerdem kennen wir den Erwartungswert desKurses S(t), nämlich exp(rT )S(0), wobei hier wieder r der risikolose Zinsund T der Fälligkeitszeitpunkt ist. Für das b berechnen wir die Varianz vonS(t) sowie die Kovarianz von S(t) und des Payoffs mittels der korrigiertenStichprobenvarianz und der korrigierten Stichprobenkovarianz wie in Glei-chung 4.5 angegeben. Die Implementation des optimalen b ist in Matlabcode8.10 realisiert. Darauf aufbauend können wir die Monte-Carlo-Simulationmit Control Variates berechnen, wie beispielhaft in Matlabcode 8.14. In Ab-bildung 4.1 ist ebenjenes berechnet und die Konfidenzintervalle wie in Satz3.5 des einfachen Monte-Carlo-Schätzers dem Schätzer mit Control Variatesgegenüber gestellt. Es wurde eine europäische Calloption mit Strike X = 100,risikolosem Zins von 10%, Volatilität von 35%, Fälligkeitszeitpunkt T = 1,Startwert S0 = 90 und jeweils 100 Zufallszahlen pro simuliertem Pfad be-schrieben. Wir betrachten das Verhalten der Konfidenzintervalle in Abhän-gigkeit von der Anzahl der simulierten Wiener Prozesse. Die Berechnungen,die zu der Abbildung 4.1 geführt haben, sind mit dem Matlabcode 8.7 erzeugtworden. Zu beachten ist hier, dass der Monte-Carlo-Schätzer mit ControlVariate einen größeren Aufwand hat, wenn man den Algorithmus mit demeinfachen Monte-Carlo-Schätzer vergleicht. Der Mehraufwand ist in etwa derGrößenordnung von m, der Anzahl der Simulationen. Deshalb stellen wir denSchätzer mit Control Variates und Simulationshäufigkeit m dem einfachenMonte-Carlo-Schätzer mit 2m Simulationen gegenüber. In der Abbildung 4.1ist auf der X-Achse dementsprechend die Anzahl der Simulationen mit Con-trol Variates abgetragen, der Standard Schätzer wurde jeweils doppelt so oftsimuliert.

52


Abbildung 4.1: Konfidenzintervalle des Standard Monte-Carlo-Schätzers unddes Schätzers mit Control Variates

Beispiel 4.3. Betrachten wir nun ein nahezu identisches Setting wie in Bei-spiel 4.2, in dem die Parameter bis auf die Anzahl der Zufallszahlen pro Pfadgleich sind. Hier nutzen wir 1500 Zufallszahlen pro Wiener-Prozess. Mit demProgrammcode 8.11 kann der näherungsweise Optionswert der EuropäischenOption aus Beispiel 4.2 in Abhängigkeit von der Anzahl der Pfade berechnetwerden. An dieser Stelle nutzen wir Control Variates und Antithetic Variatesseparat, um die Genauigkeit dieser varianzreduzierten Methoden zu zeigen.Das Ergebnis dieser Implementation liegt mit Abbildung 4.2 vor. Es ist zuerkennen, dass bei den varianzreduzierten Methoden eine deutlich geringereAnzahl an Pfadsimulationen nötig ist, um eine gewisse Genauigkeit zu erhal-ten. Außerdem sehen wir in unserem Beispiel, dass sowohl die Methode derControl Variates als auch die Methode der Antithetic Variates eine Verbes-serung der Konvergenz herbeiführt. Die modifizierten Schätzer erzeugen mitweniger Aufwand als der Standard Schätzer eine vergleichbare Genauigkeitdes Optionswertes.

53


Abbildung 4.2: Approximative Optionsbewertung mit einfacher Monte-Carlo-Methode und varianzreduzierten Monte-Carlo-Methoden

54

5 Diskretisierungsverfahren

5.1 Euler-Maruyama-VerfahrenDieses Kapitel befasst sich damit, Fehler, die bei der Diskretisierung einerSDE auftreten, zu verringern. Wir werden den Fokus auf die Herleitung undUmsetzbarkeit des Milsteinverfahrens legen, nicht jedoch auf die genauenBedingungen und Beweise der Konvergenzordnungen. Für ebenjenes sei derinteressierte Leser auf Kloeden und Platen [8] verwiesen.

Halten wir noch einmal fest, dass wir eine Lösung einer Differentialgleichungder Form

dX(t) = a(X(t))dt+ b(X(t))dW (t) (5.1)suchen. Wir haben in Kapitel 2 bereits eine einfache Approximation dieserLösung mit Gleichung 2.5 gesehen.Kloeden und Platen [8] haben dazu die Konvergenz unter bestimmten Be-dingungen bewiesen. Ziel wird es sein, das Euler-Maruyama-Verfahren

X(ti+1) = X(ti) + a(X(ti))(ti+1 − ti) + b(X(ti))√ti+1 − tiZi+1 (5.2)

zu verbessern.

5.2 MilsteinverfahrenZuerst vergegenwärtigen wir uns, dass die Approximation∫ t+h

tb(X(z))dWz ≈ b(X(t))(Wt+h −Wt) (5.3)

gilt und die Grundidee des Eulerverfahrens darstellt. Wenden wir die Ito-Formel aus Gleichung 2.13 in Integralschreibweise auf Gleichung 5.1 an:

f(X(t)) = f(X(t0) +∫ t

t0

(f ′(X(s))a(X(s)) + 1

2f′′(X(s))b(X(s))2

)ds

+∫ t

t0f ′(X(s))b(X(s))dWs

Für den Spezialfall f(x) = x ergibt sich dann

X(t) = X(t0) +∫ t

t0a(X(s))ds+

∫ t

t0b(X(s))dWs (5.4)

55

KAPITEL 5. DISKRETISIERUNGSVERFAHREN

. Wenden wir nun die Ito-Formel aus Gleichung 2.13 auf a(X(t)) und b(X(t))an:

a(X(s)) = a(X(t0)) +∫ s

t0

(a′(X(z))a(X(z)) + 1

2a′′(X(z))b(X(z))2

)dz

+∫ s

t0a′(X(z))b(X(z))dWz

b(X(s)) = b(X(t0)) +∫ s

t0

(b′(X(z))a(X(z)) + 1

2b′′(X(z))b(X(z))2

)dz

+∫ s

t0b′(X(z))b(X(z))dWz

Wobei b′(X(z)) und a′(X(z)) die ersten Ableitungen nach X(z) bezeichnen.Nun setzen wir diese Formeln in Gleichung 5.4 ein und fassen alle Doppelin-tegrale bis auf die, die bezüglich dWzdWs integriert werden, zu einem RestR1 zusammen. Es ergibt sich folgende Gleichung:

X(t) = X(t0) + a(X(t0))∫ t

t0ds+ b(X(t0))

∫ t

t0dWs (5.5)

+∫ t

t0

∫ s

t0b′(X(z))b(X(z))dWzdWs +R1

Außerdem wenden wir die Ito-Formel auf b′(X(z))b(X(z)) an und erhaltenb′(X(t0))b(X(t0)) + R2, wobei R2 wieder ein Restterm mit Integralen ist.Wie oben ziehen wir b′(X(t0))b(X(t0)) wieder aus dem Doppelintegral alsVorfaktor heraus und schreiben das Übrige als additiven Restterm. Somitfolgt aus Gleichung 5.5 mit Berechnung des ersten (Riemannschen) Integrals

X(t) = X(t0) + a(X(t0))(t− t0) + b(X(t0))∫ t

t0dWs

+ b′(X(t0))b(X(t0))∫ t

t0

∫ s

t0dWzdWs + R

mit einem Restterm R höherer Ordnung. Es bleibt die Frage, wie das Dop-pelintegral ∫ t

t0

∫ s

t0dWzdWs

berechnet wird. Für den Fall, dass t0 = 0 ist, hilft wieder das Ito-Lemma.Vereinfacht man dieses Doppelintegral zu∫ t

0

∫ s

0dWzdWs =

∫ t

0W (s)dWs

und berechnet das einfache Integral mit dem Ito-Lemma wie folgt: Sei f(x) =x2 und X(t) = W (t) also a(X(t)) = 0, b(X(t)) = 1 in der Schreibweise derIto-Formel aus Gleichung 2.13, dann ergibt sich

W (t)2 = W (0)2 + (t− 0) +∫ t

02W (s)dWs

56


, was äquivalent ist zu:

12W (t)2 − 1

2t =∫ t

0W (s)dWs =

∫ t

0

∫ s

0dWtdWs

Unter anderem Seydel [21, s.117] hat gezeigt, dass auch die allgemeine Form∫ t

t0

∫ s

t0dWtdWs = 1

2(W (t)−W (t0))2 − 12(t− t0) (5.6)

gilt. Seydel nutzt dabei Eigenschaften des Ito-Integrals, die wir jedoch nichtbesprochen haben. Es ergibt sich mit der Approximation aus Gleichung 5.3die iterative Formel

X(ti+1) = X(ti) + a(X(ti))(ti+1 − ti) + b(X(ti))√ti+1 − tiZi+1

+ 12b′(X(ti))b(X(ti))((ti+1 − ti)Z2

i+1 − (ti+1 − ti))

. Mit der Berechnung des Doppelintegrals aus Gleichung 5.6 können wirnun Algorithmus 7 für das sogenannte Milsteinverfahren angeben, um einenKursverlauf zu simulieren. Hierbei ist zu beachten, dass der Algorithmus 7sich nur geringfügig von dem Euler-Maryuama-Verfahren unterscheidet, derin Algorithmus 3 realisiert ist. Außerdem gelten in unserem Kursmodell dieGleichungen X(t) = S(t) und b′(S(t))b(S(t)) = σ2S(t) sowie alle weiterenBezeichnungen aus Kapitel 2.

Input : ∆t, µ, σ,S0, Zi ∼ N (0, 1) für i = 1, . . . , nset S(0) = S0for i = 1, . . . , n do

set dWi =√

∆tZiendfor i = 1, . . . , n do

S(ti) = S(ti−1) + µS(ti−1)∆t+ σS(ti−1)dWi + 12σ

2S(ti−1)(dW 2i −∆t)

endOutput : Vektor S ∈ Rn+1, der den Kursverlauf beinhaltetAlgorithmus 7 : Simulation eines Kurses auf Grundlage der Geometri-schen Brownschen Bewegung mittels Milstein-Verfahren

Aber was haben wir nun durch das Milsteinverfahren gewonnen? Zualler-erst sehen wir, dass unsere Diskretisierung genauer geworden ist. Wenn wireine naive Herangehensweise wählen, so könnten wir sagen, die Approximati-on wird besser, weil zu dem Term des Euler-Maruyama-Verfahrens ein Termwie in der Taylorentwicklung hinzukommt. Mathematischer sind Kloedenund Platen [8, Theorem 10.3.5] an die Diskretisierung herangegangen undhaben bewiesen, dass die starke Konvergenz des Milsteinverfahrens gegen-über dem Euler-Maruyama-Verfahren durch den additiven Term doppelt so

57


gut geworden ist. Die starke Konvergenz liegt hier mit Konvergenzordnung1 vor, das Euler-Maruyama-Verfahren hat nur Ordnung 0,5.Dennoch bereitet uns das Milsteinverfahren Probleme bei der Imlementati-on, vor allem im Höherdimensionalen. Da auch in der Praxis zumeist aufdie Verbesserung des Euler-Maruyama-Verfahrens durch das Milsteinverfah-ren verzichtet wird, verfahren wir mit der Euler-Maruyama-Diskretisierungweiter. Das hat den Vorteil, dass der Schätzer immer noch eine einfacheStruktur hat. Insbesondere falls man den Schätzer modifizieren muss, ist dasnützlich. Im mehrdimensionalen Fall ist die Anwendung des Milsteinverfah-rens deutlich komplizierter als das Euler-Maruyama-Verfahren. Ein weitererGrund, der gegen das Milsteinverfahren spricht, ist, dass wir im folgendenBeispiel 5.1 keine signifikanten Veränderung feststellen. Somit erscheint eseinfacher, die Genauigkeit des Euler-Maruyama-Verfahrens durch Erhöhungder Zeitschritte zu verbessern. Wir erzeugen mehr Zufallszahlen.

Beispiel 5.1. Betrachten wir auf Grundlage einer geometrischen Brown-schen Bewegung einen Kursverlauf mit risikolosem Zins r = 0, 1 und Volati-lität σ = 0, 34. Der Startwert des Kurses ist S0 = 75, 54 und wir simulieren1000 Zufallszahlen für eine Zeitraum von zwei Jahren. Das Ergebnis dieserSimulation ist durch Abbildung 5.1 gegeben. Hier sind die Pfade, die mit demEuler-Maruyama-Verfahren erzeugt wurden blau und die Pfade, die mit demMilsteinverfahren erzeugt wurden, sind rot. Die roten Pfade sind aber na-hezu nicht zu erkennen. Das bedeutet, dass die Werte der erzeugten Pfadesich kaum unterscheiden und stützt den Aspekt, dass wir keine signifikanteVerbesserung mit dem Milsteinverfahren erzielt haben. Auch die Endwertevon zehn Kursverläufen in Tabelle 5.1 weisen keine allzu großen Unterschie-de auf. Die größte relative Abweichung des Euler-Maruyama-Verfahrens zumMilsteinverfahrens ist in etwa 1, 1%. Sowohl die Abbildung 5.1 als auch dieErgebnisse der Tabelle 5.1 wurden mit Matlabcode 8.15 realisiert, zur Erzeu-gung der Grafik wurden drei Pfade simuliert, für die Tabelle haben wir zehnsimuliert.

58


Abbildung 5.1: Simulation eines Kurses mit dem Milsteinverfahren und mitdem Euler-Maruyama-Verfahren zu jeweils 3 Pfaden

Kurs mit Milstein Kurs mit Euler-Maruyama Abweichung relative Abweichung62,1452 61,9181 0,2271 0,36543%44,1631 44,2499 0,0868 0,19654%161,3536 161,4898 0,1362 0,08441%138,3841 139,1101 0,7260 0,52464%81,7261 82,6317 0,9056 1,10809%91,5651 92,1687 0,6036 0,65920%52,3274 52,3891 0,0617 0,11787%98,9810 100,0142 1,0332 1,04382%64,2027 64,2644 0,0617 0,09618%94,8874 94,9583 0,0710 0,07480%

Tabelle 5.1: Vergleich der Endkurse

59

6 Bewertung AsiatischerOptionen mitMonte-Carlo-Simulation

6.1 EinleitungDieses Kapitel befasst sich mit der Bewertung von Optionen einer Klasse,deren Wert sich nicht generell durch eine geschlossene Formel angeben lässt,den Asiatischen Optionen. Asiatische Optionen nutzen als Grundlage für dieAuszahlungsfunktion nicht den Schlusskurs des Underlying, sondern den Mit-telwert über die gesamte Periode bis zum Fälligkeitszeitpunkt T oder auchnur über einen Teil dieser Periode. Der Mittelwert kann arithmetisch odergeometrisch berechnet werden. Wir betrachten an dieser Stelle, wie in dergesamten Arbeit, nur europäische Auszahlungsfunktionen. Eine vorzeitigeAusübung zu betrachten erscheint bei einer Option, deren Auszahlungsfunk-tion den Mittelwert über die ganze Periode beinhaltet, nicht allzu sinnvoll.

Geometrische Asiatische OptionenGeometrische Asiatische Calloptionen haben eine Auszahlungsfunktion derForm

Λ(S(t)) = n

√√√√ n∏i=1

S(ti)−X+

. Wegen des Produkts in der Auszahlungsfunktion entstehen interessante Ef-fekte bei dem Versuch eine geschlossene Bewertungsformel herzuleiten. WieKwok [11, Kapitel 4.3.5] gezeigt hat, existiert eine geschlossene Bewertungs-formel für geometrische Asiatische Optionen. Diese Formel soll nicht Themadieser Arbeit sein, zeigt uns aber, wieso wir uns mit numerischer Simulationbeschäftigen:

1. Es gibt Optionen, die geschlossene Bewertungsformel haben, wie dieEuropäischen Plain Vanilla Optionen oder die geometrischen Asiati-schen Optionen. Für solche Optionen macht eine numerische Simulati-on wenig Sinn, da wir den Wert der Option direkt berechnen können.

61

KAPITEL 6. BEWERTUNG ASIATISCHER OPTIONEN MITMONTE-CARLO-SIMULATION

2. Es gibt Optionen, für die es keine geschlossene Bewertungsformeln gibtund für die eine numerische Simulation zur Bewertung notwendig ist.

Eine Option vom zweiten Typ stellt die arithmetische Variante der Asia-tischen Option dar, welche wir im Folgenden bewerten werden.

6.2 Bewertung einer arithmetischen AsiatischenOption

Ziehen wir also für die numerische Simulation zur Bewertung der arithme-tischen Asiatischen Option die Monte-Carlo-Simulation heran. Dazu musszunächst die Auszahlungsfunktion

Λ(S(t)) = 1n

n∑j=1

S(tj)−K+

(6.1)

bekannt sein. Mit den Erkenntnissen aus Kapitel 3 schlussfolgern wir, dass inAlgorithmus 5 nur die Auszahlungsfunktion Λ(S(t)) eingesetzt werden muss.Also ergibt sich Algorithmus 8, hier für eine Calloption.

for i=1, · · · ,m dogenerate Zi,j ∼ N (0, 1)for j = 1, · · · ,n do

set Si(tj) = Si(tj−1) + µSi(tj−1)∆t+ σSi(tj−1)√

∆tZi,jendset Λ(Si(t)) =

(1n

∑nj=1 Si(tj)−K

)+

set Ci = e−rTΛ(Si(t))endset Cm = 1

m

∑mi=1Ci

Result : näherungsweiser Optionspreis Cm

Algorithmus 8 : Approximation des Optionspreises einer Asiatischen Op-tion mit Euler-Maruyama-Verfahren

Diese Standard Methode hat eine relativ schlechte Konvergenz, was derhohen Varianz des Schätzers geschuldet ist. Diesen Effekt haben wir bereitsin Kapitel 2 sowie 4 gesehen. Wenden wir also die Varianzreduktionstechni-ken auf diesen Fall an. Dabei fällt auf, welche positiven Eigenschaften dieMonte-Carlo-Simulation besitzt: Wir können mit recht geringem (mathe-matischen) Aufwand die Schätzer modifizieren und einen Algorithmus zurBewertung einer asiatischen Option recht schnell angeben. Dazu können wirjetzt auch noch den neuen Schätzer modifizieren und Konvergenzverbesse-rung erzielen. Dazu haben wir zwei verschiedene Methoden kennengelernt.,Control Variatesn und Antithetic Variates.

62


6.3 VarianzreduktionControl VariatesZur Varianzreduktion mittels Control Variates müssen wir uns zuerst über-legen, was die Control Variate sein soll. Die Anforderung, die wir an diezusätzliche Zufallsvariable Y in der Control Variate stellen, ist, dass derErwartungswert leicht berechenbar ist. Deshalb könnten wir auf die Ideekommen, die Summen der Kurswerte Y := ∑n

i=0 Sj(ti) zum j-ten Pfad zunutzen. Wir gehen hierbei vor, wie Höllbacher [7, s.81]. Den Erwartungswertkann man recht einfach berechnen, wenn vorher klar ist, dass für Indizes iund j

E[Sj(ti)] = ertS0 (6.2)sowie ti = ∆t · i (6.3)

gilt.

E[Y ] = E[n∑i=0

Sj(ti)]

6.2=n∑i=0

ertiS0

6.3= S0

n∑i=0

(er∆t

)i= S0

er∆t·(n+1) − 1er∆t − 1

In der letzten Zeile haben wir die n-te Partialsumme der geometrischen Rei-he verwendet. Um nun den Monte-Carlo-Schätzer mit Control Variates wiein Gleichung 4.3 angeben zu können, fehlt uns noch das optimale b aus Glei-chung 4.4, welches wir direkt berechnen können:

b∗ = Cov(Λ(S(t)), Y )V ar(Y ) (6.4)

, wobei Λ(S(t)) die Auszahlungsfunktion des Standard Monte-Carlo-Schätzerbezeichnet, wie sie im vorherigem Abschnitt bereits berechnet wurde. Andieser Stelle betrachten wir die Funktion aber als Zufallsvariable zu den mPfaden.Jetzt haben wir alle Bestandteile zusammen getragen, um den Schätzer mitControl Variates analog zur Gleichung 4.3 anwenden zu können.

Cm,b∗ := b∗ · S0er∆t·(n+1) − 1er∆t − 1 + Cm − b∗ · Y (6.5)

63


, wobei Y den Monte-Carlo-Schätzer von Y bezeichnet:

Y = 1m

m∑j=1

n∑i=0

Sj(ti)

Antithetic VariatesDen Monte-Carlo-Schätzer mit Antithetic Variates Y haben wir in folgenderForm kennengelernt:

Y = 12m

m∑i=1

f(Zi) + f(−Zi)

, wobei Zi eine standardnormalverteilte Zufallszahl ist und wir die Funktionf durch die abgezinste Auszahlungsfunktion ersetzen müssen. Es ergibt sichfür die Asiatische Option also der Schätzer

Y = 12me−rT

m∑i=1

Λ(S−i (t)) + Λ(S+i (t)) (6.6)

. Hier bezeichnen Λ(S−i (t)) bzw. Λ(S+i (t)) die Auszahlungsfunktionen der

arithmetischen Asiatischen Calloption aus Gleichung 6.1, wobei die ersteFunktion mit den Zufallszahlen −Zi und die zweite Funktion mit Zi simuliertworden sind.

Abschließendes BeispielMit zwei Varianten der Varianzreduktion, der Antithetic Variates und derControl Variates, haben wir zusammen mit dem Algorithmus 8 die Möglich-keit eine arithmetische Asiatische Option zu bewerten.

Beispiel 6.1. Betrachten wir eine arithmetische Asiatische Calloption mitden Parameter Strike X = 60, risikolosem Zins r = 0, 1, Volatilität σ =0, 34, Fälligkeit in zwei Jahren und Startwert des Kurses S0 = 75, 54. Wirsimulieren den Optionspreis mit 600 Pfaden und jeweils 1500 Zufallszahlen.Dieses Vorgehen ist in Matlabcode 8.16 implementiert, der die Abbildung 6.1erzeugt und die Ergebnisse der Tabelle 6.1 berechnet.

Schätzer Wert der Option VarianzStandard Monte-Carlo Schätzer 20,8473 554,3156

Antithetic Variates 19,9933 59,2964Control Variates 19,9162 7,2071

Tabelle 6.1: Optionswerte und Varianzen mit drei Schätzern

64


Da die Varianz mit Control Variates am geringsten ist, sollten wir diesenSchätzer als besten erachten, wie wir anhand der Konfidenzintervallen sehenwerden. Um eine Aussage über die Genauigkeit der Simulation machen zukönnen, berechnen wir von dem Schätzer mit Control Variates die Ausdeh-nung des 95%-Konfidenzintervalls, wie in Satz 3.5 hergeleitet. Für das k ausGleichung 3.6 berechnen wir mit p = 0, 95%

k = Φ−1(p

2 −12

)= 1, 96

und mit der Varianz 7, 2071 beziehungsweise der Standardabweichung 2, 6846des Schätzers folgt die Ausdehnung des Konfidenzintervalls um den exaktenWert:

k · σ√m

= 1, 96 · 2, 684624, 4949 = 0, 2148

. Das bedeutet, dass der Wert 19, 9162 mit einer Wahrscheinlichkeit von 95%maximal um 0, 2148 von dem exakten Wert abweicht.

Der Wert des Schätzers mit Control Variates für 600 Simulationen, 19, 9162,ist dementsprechend der Wert, den wir als Endergebnis zur Bewertung derarithmetischen Asiatischen Option angeben wollen und das Ergebnis dieserArbeit darstellt. Abbildung 6.1 illustriert, was wir durch die deutlich schlech-tere Varianz des Standard Schätzers bereits gesehen haben, nämlich dass dieDifferenz der Optionswerte mittels varianzreduzierten Methoden gegenüberden Werten ohne Varianzreduktion stark abweichen.

65


Abbildung 6.1: Simulation des Werts einer Asiatischen Option mit varianzre-duzierten Schätzern und dem Standard Schätzer

66

7 Fazit

7.1 Zusammenfassende Bewertung und FazitWir haben begonnen die Monte-Carlo-Simulation zu untersuchen und zuverstehen, motiviert durch das Wissen um Optionen, deren Wert nicht ana-lytisch berechnet werden kann. Mit dem letzten Kapitel haben wir gesehen,dass die Überlegungen aus den vorherigen Kapiteln dazu befähigt haben,eine arithmetische Asiatische Option bewerten zu können. Wir haben dieMethoden, die wir im Laufe dieser Arbeit hergeleitet und diskutiert haben,auf die Asiatische Option angewendet.Dieses Vorgehensweise zeigt uns, dass die Euler-Maruyama-Diskretisierungund der damit einhergehende Algorithmus zur Bewertung pfadabhängigerOptionen relativ schnell für andere Optionen adaptiert werden kann. Als Bei-spiel für eine exotische, pfadabhängige Option konnten wir für die arithme-tische Asiatische Option eine Varianzreduktion mittels Control Variates undAntithetic Variates realisieren und damit die Konfidenzintervalle schrumpfenlassen. Das ist ein grundsätzliches Ziel für Schätzfunktionen und somit auchfür die pfadabhängige Optionsbewertung.Somit stellt die von uns erarbeitete Form der Monte-Carlo-Simulation zurBewertung pfadabhängiger Optionen eine Ergänzung zu den Black-Scholes-Formeln dar. Wir können den Wert einer Option berechnen, für die es keinegeschlossene Bewertungsformel gibt.Dafür war einiges an theoretischer Grundlage zu leisten. Zuallererst habenwir ein zuverlässiges Kursmodell auf Basis von Wiener Prozessen kennenge-lernt.Dazu benötigten wir standardnormalverteilte Zufallszahlen. In Abschnitt 3.3lernten wir Methoden kennen, wie Pseudozufallszahlen als gleichverteilte Zu-fallszahlen erzeugt werden und wie wir von diesen Zufallszahlen zu standard-normalverteilten Zufallszahlen kommen, mit dem Polar-Algorithmus odermit numerischer Invertierung.Darauf aufbauend konnten wir die eigentliche Monte-Carlo-Simulation for-mulieren, die wir in ihrer Varianz und damit in ihrer Konvergenz verbessernkonnten. Dabei fällt uns wieder auf: Wir erreichen relativ gute Ergebnisseund Verbesserungen durch recht einfache Strukturen. Wir haben uns zur Va-rianzreduktion nur auf Control Variates und Antithetic Variates beschränktund konnten damit sehr deutliche Verbesserungen erzielen. Somit ist auchunser anfangs gestecktes Ziel erfüllt: Wir haben einen Algorithmus hergelei-

67

KAPITEL 7. FAZIT

tet, um eine arithmetische Asiatische Option zu bewerten.Gleichzeitig fällt uns aber im Laufe dieser Arbeit auf, dass nicht alles ZurAnwendung kam, was wir behandelt haben. An dieser Stelle ist vor allemdas Milsteinverfahren zu nennen. In Kapitel 5, konnten wir keine deutlicheVerbeserung durch das Milsteinverfahren erkennen und haben diese deshalbfür die abschließende Bewertung der Asiatischen Option verworfen.

7.2 Persönliches SchlusswortIm Kontext dieser Arbeit können wir auf einen weiteren, allgemeinen Aspekthinweisen, der dieses Thema interessant macht: Die große Spannweite anThemengebieten, die zur Optionsbewertung notwendig sind.

Wir benötigten die stochastischen und finanzmathematischen Grundlagenfür die Kursentwicklung und den Monte-Carlo-Schätzer. Die Zahlentheorietrug die Zufallszahlengeneratoren bei und machte eine Bewertung dieser Zu-fallszahlen überhaupt erst möglich. Gleichzeitig griffen wir mit der Simu-lation auf numerische Methoden zurück und modifizierten die Simulation.Diese Verfahren haben wir zu Algorithmen zusammengefasst und konntensie somit in Matlab implementieren. Somit war der Weg zur Bewertung ei-ner Asiatischen Option, ein Weg durch viele verschiedene Forschungsgebieteder Mathematik, ganz im Sinne des einleitenden Zitats von Gauss aus demersten Kapitel.

68

8 Appendix: Matlabcodes

Matlabcodes

8.1 Simulation von stochastischen ProzessenDieser Programmcode realisiert die Algorithmen 1, 2 sowie 3 und gibt jeweilsdazu die Kursentwicklung als Grafik aus. Diese sind in den Abbildungen 2.1,2.2 und 2.4 festgehalten.

1 % Parameter der Brownschen Bewegungen2 r = 0 . 1 ; % r i s i k o l o s e r Zins3 sigma = 0 . 3 5 ; % Vo l a t i l i t a e t4 T = 1 ; % Fa e l l i g k e i t s z e i t p un k t5 S0 = 90 ; % Startwert der Geometrischen

Brownschen Bewegung6

7 n = 100 ;8 delta_t = T/n ;9 f o r k = 1 :3

10 randn ( ’ s t a t e ’ , k∗4) ; % j e 3 ve r s ch i edeneZu f a l l s z ah l e n

11 % erzeugen 3 ve r s ch i edene Pfade12 %% Algorithmus 1 zur Simulat ion e i n e s Wiener Proze s s e s13 W = zero s (1 , n+1) ; % Wiener Prozess14 Z = randn (n , 1 ) ; % n s tanda rdno rma lve r t e i l t e

Zu f a l l s z ah l e n15 W(1 ,1 ) = 0 ; % Startwert 016

17 f o r i = 1 : n18 W(1 , i +1) = W(1 , i ) + sq r t ( de lta_t ) ∗ Z( i , 1 ) ;19 end20 hold on ; g r i d on ;21 f i g u r e (1 )22 x l ab e l ( ’ Anzahl der Z e i t s c h r i t t e ’ ) ;23 y l ab e l ( ’Wert des Wiener Proze s s e s ’ ) ;24 p lo t (W) ; % Plot ausgeben25

69

KAPITEL 8. APPENDIX: MATLABCODES

26 %% Algorithmus 2 zur Simulat ion e i n e r BrownschenBewegung mit Dr i f t r und

27 %Vo l a t i l i t a e t sigma28 X = ze ro s (1 , n+1) ; % Brownsche Bewegung mit Dr i f t r29 % und Vo l a t i l i t a e t sigma30 X(1 ,1 ) = 0 ; % Startwert 031

32 f o r j = 1 : n33 X(1 , j +1) = X(1 , j ) + r ∗ delta_t + sigma ∗ s q r t (

de lta_t ) ∗ Z( j , 1 )34 end35 hold on ; g r i d on ;36 f i g u r e (2 )37 x l ab e l ( ’ Anzahl der Z e i t s c h r i t t e ’ ) ;38 y l ab e l ( ’Wert der Brownschen Bewegung ’ ) ;39 p lo t (X) ; % Plot ausgeben40

41 %% Algorithmus 3 zur Simulat ion e i n e r GeometrischenBrownschen Bewegung

42 % mit Dr i f t r∗S( t ) und D i f f u s i on sigma∗S( t )43 S = ze ro s (1 , n+1) ; % Geometrische Brownsche Bewegung44 S (1 , 1 ) = S0 ; % Startwert S0=9045

46 f o r k = 1 : n47 S (1 , k+1) = S (1 , k ) + S (1 , k )∗ r∗delta_t + S (1 , k )∗

sigma∗ s q r t ( de lta_t ) ∗ Z(k , 1 ) ;48 end49 hold on ; g r i d on ;50 f i g u r e (3 )51 x l ab e l ( ’ Anzahl der Z e i t s c h r i t t e ’ ) ;52 y l ab e l ( ’Wert der Geometrischen Brownschen Bewegung ’ )

;53 p lo t (S) ; % Plot ausgeben54

55 end

Matlabcode 8.1: Simulation der in Kapitel 2 vorgestellten Algorithmen

Zerlegung einer Geometrischen Brownschen BewegungIn Abbildung 2.3 haben wir eine Zerlegung eines Kurses in Trend und dieZufälligkeit um den Trend gesehen. Vorliegender Matlabcode erzeugt dieseGrafik und nutzt dazu eine Euler-Maruyama-Diskretisierung.

70


1 %% Implementation der Euler−Maruyama−Di s k r e t i s i e r ungmit Zer legung

2 % der Geometrischen Brownschen Bewegung in Dr i f t undD i f f u s i on

3 randn ( ’ s t a t e ’ ,857) ;4 r = 0 . 1 ; sigma = 0 . 3 5 ; T = 10 ; S0 = 90 ;5 n = 1000 ; de lta_t = T/n ;6 % Simulat ion e i n e s Pfades e i n e s Wiener Proze s s e s mit n

Zu f a l l s z ah l e n7 dW = sqr t ( de lta_t )∗ randn (n , 1 ) ;8 S = ze ro s (n+1 ,1) ;9 S ( 1 , : ) = S0 ; % Anfangswerte

10 % Berechnung e i n e r Geometrischen Brownschen Bewegungmit

11 % Euler−Maruyama−Di s k r e t i s i e r ung12 f o r i = 1 : n13 S( i +1 , : ) = S( i , : ) .∗ ( 1 + r∗delta_t + sigma∗dW( i , : ) ) ;14 end15 Aktienkurs = S(n+1 , : ) ;16 % Berechnung des D r i f t e i n f l u s s e s17 P = ze ro s (n+1 ,1) ;18 P(1 , 1 ) = S0 ;19 f o r i = 1 : n20 P( i +1 ,1)=P( i , 1 )+r∗delta_t ∗S( i , 1 ) ;21 end22 % Berechnung des D i f f u s i o n s e i n f l u s s e s23 Q = ze ro s (n+1 ,1) ;24 Q(1 ,1 ) = 0 ;25 f o r i = 1 : n26 Q( i +1 ,1)=Q( i , 1 )+ S( i , 1 ) ∗ sigma∗dW( i , : ) ;27 end28 % Es g i l t P+Q = S29 % gr a f i s c h e Ausgabe30 g r id on ;31 hold on ;32 p lo t (S)33 p lo t (P)34 p lo t (Q+S0 ) % Addit ion des Sta r twer t e s wegen der

U eb e r s i c h t l i c h k e i t35 % in der Graf ik36 l egend ( ’Kurs ’ , ’ E i n f l u s s des D r i f t s ’ , ’ E i n f l u s s der

D i f f u s i o n ’ ) ;

71


37 x l ab e l ( ’ Z e i t s c h r i t t e ’ ) ;38 y l ab e l ( ’Kurs ’ ) ;

Matlabcode 8.2: Zerlegung der Geometrischen Brownschen Bewegung

8.2 ZufallszahlenDer vorliegende Programmcode implementiert Beispiel 3.10 mit 20000 Zu-fallszahlen, die mit dem darauffolgendem Matlabcode 8.4 für einen Multi-plikativen Kongruenzgenerator und den Parametern von RANDU erzeugtworden sind. Es werden daraus dreidimensionale Vektoren erzeugt und alsGrafik ausgegeben (Abbildungen 3.1 und 3.2).

1 %% Diese Methode implement ie r t e in e d r e id imens i ona l eV i s u a l i s i e r ung der Zu f a l l s z ah l en ,

2 % die durch RANDU erzeugt s ind3

4 % Aufruf des mu l t i p l i k a t i v en Kongruenzgenerator RANDU5 x = MCG(65539 ,2^31 ,20000 ,1) ;6 g r id on ;7 hold on ;8 ro tate3d on ;9 t i t l e ( ’ Pseudo−Zu f a l l s z ah l e n mit RANDU ’ )

10 f o r i =1:3:1909811 s c a t t e r 3 (x ( i ) , x ( i +1) , x ( i +2) , ’ f i l l e d ’ , ’ b ’ ) ;12 end

Matlabcode 8.3: Beispiel: RANDU

1 %% Diese Funktion e rzeugt zu gegebenen Parametern e inenVektor mit

2 % Zu f a l l s z ah l e n durch e inen mu l t i p l i k a t i v enKongruenzgenerator

3 f unc t i on r e s u l t=MCG(a ,m, k ,X0)4 X = ze ro s (k , 1 ) ;5 X(1 ,1 )=X0 ; % X0 Startwert6 f o r i =1:1 : k−17 X( i +1 ,1)=(mod( a∗X( i , 1 ) ,m) ) ; % I t e r a t i o n s v o r s c h r i f t8 end9 f o r i =1:1 : k

10 X( i , 1 )=X( i , 1 ) /m;11 end12 r e s u l t = X;

Matlabcode 8.4: Multiplikativer Kongruenzgenerator

72


8.3 InvertierungDer Beasley-Springer-Moro-Algorithmus ist in Programmcode 8.5 implemen-tiert und orientiert sich an Glassermans Formulierung [3, s.68]. Hier werdenzur Illustration 200000 gleichverteilte Zufallszahlen übergeben und als Gra-fik, abhängig von der Häufigkeit, dass ein Wert angenommen wird, ausgege-ben.

1 %% Diese Methode implement ie r t den Beasley−Spr inger−Moro−Algorithmus zur

2 % Erzeugung von s tandardnorma lve r t e i l t en Zu f a l l s z ah l e nnach Glasserman

3 Zu f a l l s z ah l e n = 200000; % Anzahl der g l e i c h v e r t e i l t e nZu f a l l s z ah l e n

4 u = rand ( Zu f a l l s z ah l en , 1 ) ;5

6 %% Eingabe Parameter7 a0 = 2.50662823884 ;8 a1 = −18.61500062529;9 a2 = 41.39119773534 ;

10 a3 = −25.44106049637;11

12 b0 = −8.47351093090;13 b1 = 23.08336743743 ;14 b2 = −21.06224101826;15 b3 = 3.13082909833 ;16

17 c0 = 0.3374754822726147 ;18 c1 = 0.9761690190917186 ;19 c2 = 0.1607979714918209 ;20 c3 = 0.0276438810333863 ;21 c4 = 0.0038405729373609 ;22 c5 = 0.0003951896511919 ;23 c6 = 0.0000321767881768 ;24 c7 = 0.0000002888167364 ;25 c8 = 0.0000003960315187 ;26

27 y=u−0.5;28

29 X = ze ro s ( Zu f a l l s z ah l en , 1 ) ;30 i f ( abs ( y ) <0.42)31 r=y . ^ 2 ;32 X=y .∗ po lyva l ( [ a3 , a2 , a1 , a0 ] , r ) . / po lyva l ( [ b3 , b2 , b1 , b0 , 1 ] ,

r ) ;

73


33

34 e l s e35 r=min (u,1−u) ;36 r=log (− l og ( r ) ) ;37 X=po lyva l ( [ c8 , c7 , c6 , c5 , c4 , c3 , c2 , c1 , c0 ] , r ) ;38 X=X.∗ s i gn (y ) ; % −, wenn r <0.5 oben , a l s o min = r39 % +, wenn r >0.5 oben , a l s o min = 1−r40 end41 histogram (X) % g r a f i s c h e Ausgabe42 g r id on ;43 y l ab e l ( ’ Anzahl Zu f a l l s z ah l e n ’ ) ;44 x l ab e l ( ’Wert der Zu f a l l s z ah l e n ’ ) ;

Matlabcode 8.5: Beasley-Springer-Moro-Algorithmus

8.4 Polar-AlgorithmusDer Polar-Algorithmus in Matlabcode 8.6 erzeugt mit der Verwerfungsme-thode aus gleichtverteilten Zufallszahlen standardnormalverteilte Zufallszah-len. Wir haben hier einen Input von 200000 gleichverteilten Zufallszahlen,wie im vorhergehenden Beasley-Springer-Moro-Agorithmus.

1 %% Diese r Code implement ie r t den Polara lgor i thmus zurErzeugung von

2 % standardnorma lve r t e i l t en Zu f a l l s z ah l e n3

4 n=100000;5 U1=rand (n , 1 ) ; % Input : g l e i c h v e r t e i l t e Zu f a l l s z ah l e n6 U2=rand (n , 1 ) ; % Input : g l e i c h v e r t e i l t e Zu f a l l s z ah l e n7 V1=ze ro s (n , 1 ) ;8 V2=ze ro s (n , 1 ) ;9 V=ze ro s (n , 1 ) ;

10 j = 1 ;11 f o r i =1:n12

13 V1( i ) = 2 ∗ U1( i ) − 1 ; % g l e i c h v e r t e i l t e Zu f a l l s z ah l e nauf [−1 ,1]

14 V2( i ) = 2 ∗ U2( i ) − 1 ; % g l e i c h v e r t e i l t e Zu f a l l s z ah l e nauf [−1 ,1]

15 V( i ) = V1( i )∗V1( i ) + V2( i )∗V2( i ) ;16

17 i f (0<V( i ) && V( i )<=1) % Verwerfung oder Akzept ierender Zahlen

74


18 Z1( j ) = V1( i ) ∗ s q r t (−2∗ l og (V( i ) ) /V( i ) ) ;19 Z2( j ) = V2( i ) ∗ s q r t (−2∗ l og (V( i ) ) /V( i ) ) ;20 j = j +1; % j za eh l t d i e ak z ep t i e r t en

Zahlen21 end22 end23

24 T = ze ro s (2∗n , 1 ) ;25 T = [ Z1 , Z2 ] ;26 histogram (T) ; % g r a f i s c h e Ausgabe27 g r id on ;28 y l ab e l ( ’ Anzahl Zu f a l l s z ah l e n ’ ) ;29 x l ab e l ( ’Wert der Zu f a l l s z ah l e n ’ ) ;

Matlabcode 8.6: Polar-Algorithmus

8.5 Vergleich der KonfidenzintervalleDieser Programmcode berechnet zu verschieden vielen Pfaden die Konfidenz-intervalle der Standard Monte-Carlo-Methode und der Methode mit ControlVariates.

1 X = 100 ; r = 0 . 1 ; sigma = 0 . 3 5 ; T = 1 ; S0 = 90 ; n =100 ; t = 0 ;

2 KonfiInterval l_CC1 = ze ro s (5 , 1 ) ;3 KonfiInterval l_CC2 = ze ro s (5 , 1 ) ;4 KonfiInterval l_MC_einfach1 = ze ro s (5 , 1 ) ;5 KonfiInterval l_MC_einfach2 = ze ro s (5 , 1 ) ;6 P = 0 . 9 5 ; % 95%−Kon f i d en z i n t e r va l l e7 % Quanti l der Standardnormalverte i lung8 k = norminv (1/2 + P/2) ;9 f o r i =1:5

10 % Berechnung der Kon f i d en z i n t e r va l l e f u e rve r s ch i edene Anzahlen

11 % an s imu l i e r t e n Pfaden . 100 b i s 500 Pfade12 M = i ∗100 ;13 v1 = VarianzCC (X, r , sigma ,T, S0 , n ,M) ;14 v2 = VarianzMC_einfach (X, r , sigma ,T, S0 , n , 2∗M) ;15 vexakt = c a l l ( S0 , t ,X, r , sigma ,T)16 KonfiInterval l_CC1 ( i , 1 ) = k∗ s q r t ( v1/M ) ;17 KonfiInterval l_CC2 ( i , 1 ) = − k∗ s q r t ( v1/M ) ;18 KonfiInterval l_MC_einfach1 ( i , 1 ) = k∗ s q r t ( v2/M ) ;

75


19 KonfiInterval l_MC_einfach2 ( i , 1 ) = −k∗ s q r t ( v2/M );

20 end21 hold on ;22 x l ab e l ( ’ Anzahl der Simulat ionen ’ ) ; % Achsenbeschr i f tung23 y l ab e l ( ’Ausdehnung der Kon f i d en z i n t e r va l l e ’ ) ; %

Achsenbeschr i f tung24 g r id on ;25

26 ax = gca ;27 bar ( KonfiIntervall_MC_einfach1 , . 5 , ’b ’ ) ;28 bar ( KonfiIntervall_CC1 , . 2 5 , ’ c ’ ) ;29 bar ( KonfiIntervall_MC_einfach2 , . 5 , ’b ’ ) ;30 bar ( KonfiIntervall_CC2 , . 2 5 , ’ c ’ ) ;31 l egend ( ’ Kon f i d en z i n t e r va l l des e in fachen Monte−Carlo−

S c h tzers ’ , ’ Kon f i d en z i n t e r va l l mit Contro le Var ia t e s’ ) ;

32 ax . XTickLabel = ’ ’ , ’ 100 ’ , ’ ’ , ’ 200 ’ , ’ ’ , ’ 300 ’ , ’ ’ , ’ 400 ’ , ’ ’, ’ 500 ’ , ’ ’ , ’ 600 ’ ;

33 ax . YTickLabel = ’ ’ , ’C−4 ’ , ’C−3 ’ , ’C−2 ’ , ’C−1 ’ , ’C = 12 ,296’ , ’C+1 ’ , ’C+2 ’ , ’C+3 ’ , ’C+4 ’ , ’ ’ ;

Matlabcode 8.7: Vergleich der Konfidenzintervalle mit Control Variate

Berechnung der VarianzenDie folgenden beiden Programmcodes sind zur Berechnung der Konfidenz-intervalle in Matlabcode 8.7 nötig. Grundlage dafür bietet wieder die Euler-Maruyama-Diskretisierung. In Matlabcode 8.8 wird das optimale b (wie inGleichung 4.4), das mittels Matlabcode 8.10 berechnet wird, benutzt.

1 f unc t i on Varianz_des_CC_Schaetzers = VarianzCC (X, r ,sigma ,T, S0 , n ,M)

2

3 randn ( ’ s t a t e ’ , 3 ) ;4 h = T/n ;5 % Erzeugung der Wiener−Prozes se f u e r M Pfaden6 dW = sqr t (h)∗ randn (n ,M) ;7 % Berechnung der Aktienkurse f u e r M Pfade8 S = ze ro s (n+1,M) ;9 S ( 1 , : ) = S0 ; % Anfangswerte

10 f o r i = 1 : n11 S( i +1 , : ) = S( i , : ) .∗ ( 1 + r∗h + sigma∗dW( i , : ) ) ;

76


12 end13 Aktienkurse = S(n+1 , : ) ; % Aktienkurs zum Zeitpunkt T

fue r M v i e l e Pfade14 Auszahlung = max(0 , S(n+1 , : )−X) ; % Berechnung der

Payo f f s15

16 E = S(1) .∗ exp ( r∗T) ;17 bn = optimales_b_approximativ ( Auszahlung , S(n+1 , : ) ) ;18 % Berechnung der Control Var iate19 Contro lVar iate = Auszahlung − bn∗(S(n+1 , :) − E) ;20

21 VCC = exp(−r∗T)∗ (mean( Contro lVar iate ) ) ;22 % ko r r i g i e r t e St i chprobenvar ianz b e r

Auszahlungsfunkt ion mit Control Var ia t e s23 Varianz_des_CC_Schaetzers = var ( Contro lVar iate ) ;24 end

Matlabcode 8.8: Varianz des Schätzers mit Control Variate

1 f unc t i on Varianz_des_MC_Schaetzers = VarianzMC_einfach (X, r , sigma ,T, S0 , n ,M)

2

3 randn ( ’ s t a t e ’ , 3 ) ;4 h = T/n ;5 % Erzeugung der Wiener−Prozes se zu M Pfaden6 dW = sqr t (h)∗ randn (n ,M) ;7 %Berechnung der Aktienkurse f u e r a l l e Pfade8 S = ze ro s (n+1,M) ;9 S ( 1 , : ) = S0 ; % Anfangswerte

10 f o r i = 1 : n11 S( i +1 , : ) = S( i , : ) .∗ ( 1 + r∗h + sigma∗dW( i , : ) ) ;12 end13 Aktienkurse = S(n+1 , : ) ; % Aktienkurs zum Zeitpunkt T

fue r M Pfade14 Auszahlung = max(0 , S(n+1 , : )−X) ; % Berechnung der

Payo f f s15

16 % abge z i n s t e r Wert des Monte−Carlo−Schae tze r s17 V = exp(−r∗T)∗ (mean( Auszahlung ) ) ;18

19 % ko r r i g i e r t e St i chprobenvar ianz b e r d i e Auszahlungen20 Varianz_des_MC_Schaetzers = var ( Auszahlung ) ;21 end

Matlabcode 8.9: Varianz des einfachen Monte-Carlo-Schätzers

77


1 f unc t i on r e s u l t= optimales_b_approximativ (X,Y)2 % Berechnung der k o r r i g i e r t e n St i chprobenvar ianz3 % ( ohne Faktor 1/( l ength (Y)−1) )4 Zaehler1 = ze ro s ( l ength (Y) ) ;5 f o r i =1: l ength (Y)6 Zaehler1 ( i ) = (Y( i )−mean(Y) ) ^2;7 end8 Var = sum( Zaeh ler1 ) /( l ength (Y)−1) ;9 % ( wieder ohne Vorfaktor )

10 Zaehler2 = ze ro s ( l ength (Y) ) ;11 f o r i =1: l ength (Y)12 Zaehler2 ( i ) = (X( i )−mean(X) ) ∗ (Y( i )−mean(Y) ) ;13 end14 Cov = sum( Zaeh ler2 ) /( l ength (Y)−1) ;15 r e s u l t = Cov/Var ;16 end

Matlabcode 8.10: Approximativ optimales b

8.6 Vergleich der approximativen Optionswerteder einfachen Monte-Carlo-Simulation mitvarianzreduzierten Methoden

Um die Optionswerte, die approximativ mit einfacher Monte-Carlo-Simulation,mit Antithetic Variates und mit Control Variates berechnet werden, abhän-gig von verschieden vielen Pfaden zu vergleichen, greift der vorliegende Codeauf die drei folgenden Funktionen zurück. Diese berechnen jeweils den Opti-onswert eines europäischen Calls zu gegebener Anzahl an Pfadsimulationendes Wiener Prozesses. Matlabcode 8.11 ruft diese Funktionen mit den Para-metern aus Beispiel 4.2 und Beispiel 4.3 und Werten für M von 100 bis 2000und erstellt dazu eine Grafik, Abilldung 4.2.

1 X = 100 ; r = 0 . 1 ; sigma = 0 . 3 5 ; T = 1 ; S0 = 90 ; n =1500 ; t = 0 ;

2 M = 100 ;3 S ch r i t t e = 20 ;4 Q = ze ro s (4 , S ch r i t t e ) ;5

6 f o r i =1: S ch r i t t e7 Q( : , i ) = Optionswerte (X, r , sigma ,T, S0 , n , i ∗M) ;8 end9

78


10 f i g u r e (1 ) ;11 ax = gca ;12 x = 1 : S ch r i t t e ;13 p lo t (x ,Q) ;14 x l ab e l ( ’ Anzahl der Simulat ionen ’ ) ; % Achsenbeschr i f tung15 y l ab e l ( ’Wert in ’ ) ;16 g r id on17 l egend ( ’ e i n f a ch e r Monte−Carlo− S c h tzers ’ , ’ Ant i th e t i c

Var ia te s ’ , ’ Control Var ia t e s ’ , ’ exakt nach Black−Scho l e s ’ ) ;

18 ax . XTickLabel=0 ,200 ,400 ,600 ,800 ,1000 ,1200 ,1400 ,1600 ,1800 ,2000 ;

Matlabcode 8.11: Monte-Carlo-Simulation zur Optionsbewertung einfach, mitControl Variates und mit Antithetic Variates

Berechnung einer europäischen Calloption mittelsMonte-Carlo mit Antithetic Variates

1 f unc t i on Optionswerte = Optionswerte (X, r , sigma ,T, S0 , n ,M)

2 randn ( ’ s t a t e ’ , 3 )3 h = T/n ;4 % Berechnung der Wiener−Prozes se f u e r M Pfaden5 dW = sqr t (h)∗ randn (n ,M) ;6 % Berechnung der Aktienkurse f u e r M Pfade7 S1 = ze ro s (n+1,M) ;8 S1 ( 1 , : ) = S0 ; % Anfangswerte9 S2 = S1 ;

10 f o r i = 1 : n11 S1 ( i +1 , : ) = S1 ( i , : ) .∗ ( 1 + r∗h + sigma∗dW( i , : ) ) ;12 S2 ( i +1 , : ) = S2 ( i , : ) .∗ ( 1 + r∗h − sigma∗dW( i , : ) ) ;13 end14 % Berechnung der Auszahlungsfunkt ion15 payo f f = 0 . 5∗ (max(0 , S1 (n+1 , : )− X ) + max(0 , S2 (n+1 , :) −

X) ) ;16 % Berechnung des S c h tzers und der Opt ionspre i s e17 % Standard S c h tzer18 V_einfach = MCeuler (X, r , sigma ,T, S0 , n ,2∗M) ;19 % Ant i th e t i c Var iate20 V = exp(−r∗T)∗ mean( payo f f ) ;21 % Control Var iate

79


22 V_cc = MCeuler_cal l_ControlVariates (X, r , sigma ,T, S0 , n ,M);

23 % exakter Wert , mit Black−Scho l e s24 Vexakt = c a l l ( S0 , 0 ,X, r , sigma ,T) ;25 Optionswerte = [ V_einfach V V_cc Vexakt ] ;26

27 end

Matlabcode 8.12: Monte-Carlo-Simulation zur Optionsbewertung mitAntithetic Variates

Berechnung einer europäischen Calloption mittelseinfacher Monte-Carlo-Simulation

1 f unc t i on r e s u l t = MCeuler (X, r , sigma ,T, S0 , n ,M) %Berechnung e i n e s Ca l l s

2 randn ( ’ s t a t e ’ , 3 ) ;3 h = T/n ;4 % Berechnung der Wiener−Prozes se zu M Pfaden5 dW = sqr t (h)∗ randn (n ,M) ;6 % Berechnung der Aktienkurse f u e r a l l eP f ade7 S = ze ro s (n+1,M) ;8 S ( 1 , : ) = S0 ; % Anfangswerte9 f o r i = 1 : n

10 S( i +1 , : ) = S( i , : ) .∗ ( 1 + r∗h + sigma∗dW( i , : ) ) ;11 end12 Aktienkurs = S(n+1 , : ) ;13 % Berechnung der Auszahlungsfunkt ion14 Payof f = max(0 , S(n+1 , : )−X) ;15 % Berechnung des Schae tze r s und der Opt ionspre i s e16 V = exp(−r∗T) ∗(cumsum( Payof f ) . / ( 1 :M) ) ;17 var (V) ;18 r e s u l t = sum(V) /M;

Matlabcode 8.13: Monte-Carlo-Simulation zur Optionsbewertung mit demeinfachen Monte-Carlo-Schätzer

Berechnung einer europäischen Calloption mittelsMonte-Carlo mit Control Variates

1 f unc t i on r e s u l t = MCeuler_cal l_ControlVariates (X, r ,sigma ,T, S0 , n ,M)

80


2 randn ( ’ s t a t e ’ , 3 ) ;3 h = T/n ;4 % Erzeugung der Wiener−Prozes se zu M Pfaden5 dW = sqr t (h)∗ randn (n ,M) ;6 % Berechnung der Aktienkurse f u e r a l l e Pfade7 S = ze ro s (n+1,M) ;8 S ( 1 , : ) = S0 ; % Anfangswerte9 f o r i = 1 : n

10 S( i +1 , : ) = S( i , : ) .∗ ( 1 + r∗h + sigma∗dW( i , : ) ) ;11 end12 Aktienkurs = S(n+1 , : ) ; % Aktienkurs zum Zeitpunkt T f r

M−v i e l e Pfade13

14 % Berechnung der Auszahlungsfunkt ion15 Auszahlung = max(0 , S(n+1 , : )−X) ; % Payof f16 E = S ( 1 , : ) .∗ exp ( r∗T) ;17

18 bn = optimales_b_approximativ ( Auszahlung , S(n+1 , : ) ) ;19

20 % Berechnung der Control Var iate21 Contro lVar iate = Auszahlung − bn∗(S(n+1 , :) − E) ;22 V = exp(−r∗T)∗ (cumsum( Contro lVar iate ) . / ( 1 :M) ) ;23 var (V) ;24 %25 r e s u l t = sum(V) /M;

Matlabcode 8.14: Monte-Carlo-Simulation zur Optionsbewertung mit ControlVariates

8.7 Vergleich der Kurssimulationen mit demEuler-Maruyama-Verfahren und mit demMilstein-Verfahren

1 %% Diese Methode implement ie r t das Mi l s t e i nv e r f ah r en2 % und das Euler−Maruyama−Verfahren3 % zur D i s k r e t i s i e r ung s t o c h a s t i s c h e r D i f f e r e n t i a l−4 % gle ichungen , angwandt auf e inen Kursver lau f m i t t e l s5 % geometr i s che r Brownschen Bewegung6

7 randn ( ’ s t a t e ’ , 9 ) ;8 sigma=0.34; r =0.1 ;T = 2 ; S0 = 75 . 5 4 ;

81


9 n = 1000 ; h = T/n ;10

11 f o r j =1:3 % Zur Erzeugung von d r e i Pfaden12

13 Z = randn (n , 1 ) ;14 dW = ze ro s (n , 1 ) ;15 f o r i =1:n16 dW( i ) = sq r t (h)∗Z( i , 1 ) ;17 end18

19 Mi l s t e i n = ze ro s (n+1 ,1) ;20 Euler = ze ro s (n+1 ,1) ;21 Euler (1 , 1 ) = S0 ; % Anfangswerte22 Mi l s t e i n (1 , 1 ) = S0 ; % Anfangswerte23 F = ze ro s (n , 1 ) ;24 f o r i = 1 : n % i t e r a t i v e Berechnungsvor schr i f t en25 F( i ) = 0 .5 ∗ sigma∗ sigma ∗ ( dW( i )∗dW( i ) − h) ;26 Mi l s t e i n ( i +1) = Mi l s t e i n ( i ) ∗(1 + r∗h + sigma∗dW( i )

+ F( i ) ) ;27 Euler ( i +1) = Euler ( i ) ∗(1 + r∗h + sigma∗dW( i ) ) ;28 end29 % Endkurse mit dem Mi l s t e i nv e r f ah r en30 Aktienkurs_Milste in = Mi l s t e i n (n+1)31 % Endkurse mit dem Euler−Maruyama−Verfahren32 Aktienkurs_Euler = Euler (n+1)33 g r id on ;34 hold on ;35 % Gra f i s che Ausgabe36 p lo t ( Euler , ’ red ’ ) ;37 p lo t ( Mi l s t e in , ’ b lue ’ ) ;38

39 end40 l egend ( ’ Kursver lau f mit Euler−Maruyama−Verfahren ’ , ’

Kursver lau f mit dem Mi l s t e i nv e r f ah r en ’ ) ;41 x l ab e l ( ’ Z e i t s c h r i t t e ’ ) ;42 y l ab e l ( ’ Kurswert ’ ) ;

Matlabcode 8.15: Simulation eines Kurses mit dem Euler-Maruyama-Verfahren und mit dem Milsteinverfahren

Dieser Matlabcode implementiert das Milsteinverfahren und das Euler-Maruyama-Verfahren zu Simulation eines Kurses über zwei Jahre auf Grundlage einerGeometrischen Brownschen Bewegung mit Startwert S0 = 75, 54, risikolosemZins r = 0, 1 und Volatilität σ = 0, 34 Wir erzeugen jeweils drei Pfade und

82


geben den Kurs des Euler-Maruyama-Verfahrens in rot in der Grafik aus undden Verlauf des Milsteinverfahrens in blau.

8.8 Bewertung Asiatischer Optionen

1 %% Berechnung des Wertes e i n e r As i a t i s chen Ca l l opt i onm i t t e l s

2 % Control Var iates , Ant i th e t i c Var ia t e s und StandardMethode

3 % sowie Varianzberechnung der Schaetzer4 randn ( ’ s t a t e ’ , 3 ) ;5 X = 60 ; r = 0 . 1 ; sigma = 0 . 3 4 ; T = 2 ; S0 = 75 . 5 4 ;6 n = 1500 ; h = T/n ;7 Pfade = 600 ;8 % Berechnung der Wiener−Prozes se zu M Pfaden9 dW = sqr t (h)∗ randn (n , Pfade ) ;

10 % Berechnung der Aktienkurse f u e r a l l e Pfade11 S = ze ro s (n+1,Pfade ) ;12 S ( 1 , : ) = S0 ; % Anfangswerte13 S1 = S ; % S1 i s t d i e Ant i th e t i c Var iate14 f o r i = 1 : n15 S( i +1 , : ) = S( i , : ) .∗ ( 1 + r∗h + sigma∗dW( i , : ) ) ;16 S1 ( i +1 , : ) = S1 ( i , : ) .∗ ( 1 + r∗h − sigma∗dW( i , : ) ) ;17 end18

19 F = ze ro s ( Pfade , 1 ) ;20 f o r M = 1 : Pfade21

22 mit te l_ar i thmet i s ch = ze ro s (1 ,M) ;23 mitte l_ar i thmet i sch_1 = ze ro s (1 ,M) ;24

25 f o r j = 1 :M26 mit te l_ar i thmet i s ch (1 , j ) = 1/(n+1) ∗ sum(S ( : , j ) ) ;27 mitte l_ar i thmet i sch_1 (1 , j ) = 1/(n+1) ∗ sum(S1 ( : , j ) ) ;28 end29

30 % Berechnung des Payof f der As i a t i s chen Option31 Auszahlung_arithmetisch = max(0 , mit te l_ar i thmet i sch−X) ;32 Auszahlung_ar i thmet i sch_ant i thet ic = 0 .5 ∗ ( max(0 ,

mit te l_ar i thmet i sch−X) + max(0 , mitte l_ar ithmetisch_1−X) ) ;

33 % Berechnung des Schae tze r s und der Opt ionspre i s e

83


34 V_arithmetisch = exp(−r∗T) ∗(cumsum(Auszahlung_arithmetisch ) . / ( 1 :M) ) ;

35 V_ar i thmet i sch_ant i thet ic = exp(−r∗T) ∗(cumsum(Auszahlung_ar i thmet i sch_ant i thet ic ) . / ( 1 :M) ) ;

36

37

38 E = S0∗( exp ( r∗h∗(n+1) )−1)/( exp ( r∗h)−1) ; %Erwartungswert von Y

39 Y = ze ro s (M, 1 ) ; % zu s a e t z l i c h e Z u f a l l s v a r i a b l e derControl Var iate

40

41 f o r j =1:M42 Y( j , 1 ) = sum(S ( : , j ) ) ;% Summe der Kurswerte zum j−

ten Pfad43 end44 Y_dach = 1/M ∗ sum(Y) ; % Monte−Carlo−Schaetzer der

Kurswerte45

46 % opt imales b47 b = optimales_b_approximativ ( Auszahlung_arithmetisch ,Y)

;48

49 % ge s cha e t z t e r Optionswert mit Standard Monte−Carlo−Methode

50 C_m_dach = V_arithmetisch (M)51

52 % ge s cha e t z t e r Optionswert m i t t e l s Monte−Carlo−53 % Methode mit Ant i th e t i c Var iate54 C_m_dach_antithetic = V_ar i thmet i sch_ant i thet ic (M)55

56 % approximat iver Optionswert mit Control Var iate57

58 F(M)= b∗E + C_m_dach − b∗Y_dach59

60 Varianz_des_Standard_Schaetzer (M) = var (Auszahlung_arithmetisch ) ;

61 Varianz_des_Control_Variate_Schaetzer (M) = var (Auszahlung_arithmetisch . ’ − b∗ Y + b∗E) ;

62 Varianz_des_Antithetic_Variate_Schaetzer (M) = var (Auszahlung_ar i thmet i sch_ant i thet ic ) ;

63

64

65 g r id on ;

84


66 hold on ;67 f i g u r e (3 ) ;68

69 p lo t ( V_arithmetisch ) ; % Standard Schaetzer70 p lo t ( V_ar i thmet i sch_ant i thet ic ) ; % Schaetzer mit

Ant i th e t i c Var iate71 end72 hold on ;73 p lo t (F) ; % Schaetzer mit Control Var ia t e s74 l egend ( ’ Standard S c h tzer ’ , ’ S c h tzer mit Ant i th e t i c

Var ia te s ’ , ’ S c h tzer mit Control Var ia t e s ’ ) ;75 x l ab e l ( ’ S imulat ionen ’ ) ;76 y l ab e l ( ’Wert der Option ’ ) ;77

78 Varianz_des_Standard_Schaetzer ( Pfade )79 Varianz_des_Control_Variate_Schaetzer ( Pfade )80 Varianz_des_Antithetic_Variate_Schaetzer ( Pfade )

Matlabcode 8.16: Simulation des Werts einer Asiatischen Option

Dieser Matlabcode implementiert das Beispiel 6.1 und erzeugt die Abbildung6.1. Hierbei werden 2000 Simulationen durchgeführt und die errechneten Op-tionswerte mittels Standard Schätzer, den Werten des Schätzers mit ControlVariates und des Schätzers mit Antithetic Variates gegenüber gestellt.

85

Literaturverzeichnis[1] Devroye, Luc: Non-uniform random variate generation. Springer, 1986.

[2] Einstein, Albert: Über die von der molekularkinetischen Theorie derWärme geforderte Bewegung von in ruhender Flüssigkeit suspendiertenTeilchen. Annalen der Physik, 1905.

[3] Glasserman, Paul: Monte Carlo Methods in Financial Engineering.Springer, 2003.

[4] Günther, Michael und Ansgar Jüngel: Finanzderivate mit MATLAB R©.Vieweg + Teubner, 2. Auflage, 2010.

[5] Henze, Norbert: Stochastik für Einsteiger. Springer, 10. Auflage, 2013.

[6] Hesse, Christian: Angewandte Wahrscheinlichkeitstheorie. vieweg,1. Auflage, 2003.

[7] Höllbacher, Thomas: Hedging mit Monte Carlo Algorithmen. Diplom-arbeit, Universität Bayreuth, 2011.

[8] Kloeden, Peter und Eckhard Platen: Numerical Solution of StochasticDifferential Equations. Springer, 3. Auflage, 1999.

[9] Knuth, Donald E.: The Art of Computer Programming. Addison-Wesley,3. Auflage, 1998.

[10] Korn, Ralf: Moderne Finanzmathematik, Band 1. Springer Spectrum,2014.

[11] Kwok, Yue Kuen: Mathematical Models of Financial Derivatives. Sprin-ger Finance, 2. Auflage, 2008.

[12] Lehmer, Derrick E.: Mathematical methods in large-scale computingunits. Proceedings of a Second Symposium on Large-Scale Digital Cal-culating Machinery, Seiten 141–146, 1949.

[13] Lewis, Peter, Allen Goodman und James Miller: A peudo-random num-ber generator for the System/360. IBM Systems Journal, 1969.

[14] Marsaglia, George: Random numbers fall mainly in in the planes. Pro-ceedings of the National Academy of Sciences, 1968.

87

Literaturverzeichnis

[15] Marsaglia, George und Thomas A. Bray: A Convenient Method for Ge-nerating Normal Variables. SIAM Review, Seiten 260–264, Juli 1964.

[16] Matsumoto, Makoto und Takuji Nishimura. Letzter Zugriff: 31.08.2016,1997. http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html.

[17] Mosler, Karl und Friedrich Schmid: Wahrscheinlichkeitsrechnung undschließende Statistik. Springer, 2. Auflage, 2006.

[18] Müller-Gronbach, Thomas, Erich Novak und Klaus Ritter:Monte Carlo-Algorithmen. Springer, 2012.

[19] Øksendal, Bernt: Stochastic differential equations : an introduction withapplications. Springer, 6. Auflage, 2007.

[20] Porembski, Marcus: Finanzmathematik. Philipps-Universität Marburg,2015. Vorlesungsskript.

[21] Seydel, Rüdiger: Tools for Computational Finance. Springer, 5. Auflage,2011.

[22] Smulochowski, Marian: Zur kinetischen Theorie der Brownschen Mole-kularbewegung und der Suspension. Annalen der Physik, 1906.

[23] Wichmann, Brian und David A. Hill: An Efficient and Portable Pseudo-Random Number Generator. Applied Statistics, Seiten 188–190, 1982.

88

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

Eidesstattliche ErklärungIch erkläre, dass ich die vorliegende Bachelorarbeit „Monte-Carlo-Simulationzur Optionsbewertung“ selbst verfasst und dazu keine anderen als die ange-führten Hilfsmittel verwendet habe. Ich erkläre weiterhin, dass die vorliegen-de Arbeit weder in gleicher noch in ähnlicher Form einer Prüfungsbehördevorgelegt wurde.

Ellershausen, den 31.08.2016

Nils Benedikt Kollmar

Monte-Carlo-Simulation zur Optionsbewertung - uni-marburg.de · Monte-Carlo-Simulation zur...

Documents

Transcript of Monte-Carlo-Simulation zur Optionsbewertung - uni-marburg.de · Monte-Carlo-Simulation zur...