Quasi-Monte-Carlo-Verfahren zur Bewertung von · PDF fileQuasi-Monte-Carlo-Verfahren zur...

Quasi-Monte-Carlo-Verfahren zur

Bewertung von Finanzderivaten

Bachelorarbeit

eingereicht bei

Prof. Dr. Thomas Gerstner

Fachbereich Informatik und Mathematik

Institut fur Mathematik

Goethe Universitat

Frankfurt am Main

von:

Benny Xiang Li

An den Hohlgarten. 23

61138 Niederdorfelden

Matrikelnummer: 3516274

Studienrichtung: Mathematik

Niederdorfelden, 14. Juni 2011

Inhaltsverzeichnis

Abbildungsverzeichnis III

Danksagung V

Ehrenwortliche Erklarung V

1 Einleitung 1

1.1 Aufbau der Arbeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Grundlagen 5

2.1 Marktannahmen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Optionstypen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2.1 Pfadunabhangige Optionen . . . . . . . . . . . . . . . . . . . . 6

2.2.2 Pfadabhangige Optionen . . . . . . . . . . . . . . . . . . . . . . 6

2.3 Options-Bewertungsmethodik . . . . . . . . . . . . . . . . . . . . . . . 8

2.3.1 Risikoneutraler Ansatz . . . . . . . . . . . . . . . . . . . . . . . 8

2.3.2 Black-Scholes-Ansatz . . . . . . . . . . . . . . . . . . . . . . . . 9

2.4 Stochastik und ihre numerische Behandlung . . . . . . . . . . . . . . . 10

3 Quasi-Monte-Carlo-Verfahren 15

3.1 Niederdiskrepanz-Folgen . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.1.1 Halton-Folge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.1.2 Faure-Folge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.1.3 Sobol-Folge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2 Theoretische Eigenschaft . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.2.1 Diskrepanz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2.2 Fehlerschranken . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.3 Praktische Eigenschaft . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

I

II INHALTSVERZEICHNIS

3.3.1 Integral-Approximation . . . . . . . . . . . . . . . . . . . . . . . 37

3.3.2 Vergleich Konvergenzverhalten . . . . . . . . . . . . . . . . . . . 38

3.4 Randomisierte Quasi-Monte-Carlo-Verfahren . . . . . . . . . . . . . . . 40

3.4.1 Randomisierungsansatz von Tuffin . . . . . . . . . . . . . . . . . 40

3.4.2 Erwartungswert und Varianz . . . . . . . . . . . . . . . . . . . . 40

3.4.3 Fehler des Schatzers . . . . . . . . . . . . . . . . . . . . . . . . 41

4 Numerische Ergebnisse 43

4.1 Implementierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.2 Bewertung von Optionen . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.2.1 Europaische Optionen . . . . . . . . . . . . . . . . . . . . . . . 45

4.2.2 Asiatische Optionen diskret geometrisch . . . . . . . . . . . . . 48

4.2.3 Asiatische Optionen kontinuierlich geometrisch . . . . . . . . . . 52

5 Zusammenfassung und Ausblick 57

A Anhang: Programmcodes 59

A.1 Niederdiskrepanz-Folgen-Generatoren . . . . . . . . . . . . . . . . . . . 59

A.2 Black-Scholes-Formel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

A.3 Integral Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

A.4 Europaische Optionen . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

A.5 Asiatische Optionen diskret geometrisch . . . . . . . . . . . . . . . . . 71

A.6 Asiatische Optionen kontinuierlich geometrisch . . . . . . . . . . . . . . 75

Literaturverzeichnis 79

Abbildungsverzeichnis

3.1 Halton-2d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.2 Halton-100d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.3 Faure-2d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.4 Faure-100d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.5 Sobol-2d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.6 Sobol-100d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.7 Integral Konvergenzverhalten . . . . . . . . . . . . . . . . . . . . . . . 38

3.8 Integral Lineare Regressionsanalyse . . . . . . . . . . . . . . . . . . . . 39

4.1 Euro-Call Euler-Maruyama-Ansatz Konvergenzverhalten . . . . . . . . 45

4.2 Euro-Call Euler-Maruyama-Ansatz Lineare Regressionsanalyse . . . . . 46

4.3 Euro-Call Ito-Ansatz Konvergenzverhalten . . . . . . . . . . . . . . . . 47

4.4 Euro-Call Ito-Ansatz Lineare Regressionsanalyse . . . . . . . . . . . . . 47

4.5 Asia-Call diskret Euler-Maruyama-Ansatz Konvergenzverhalten . . . . 49

4.6 Asia-Call diskret Euler-Maruyama-Ansatz Lineare Regressionsanalyse . 49

4.7 Asia-Call diskret Ito-Ansatz Konvergenzverhalten . . . . . . . . . . . . 50

4.8 Asia-Call diskret Ito-Ansatz Lineare Regressionsanalyse . . . . . . . . . 51

4.9 Asia-Call kontinuierlich Euler-Maruyama-Ansatz Konvergenzverhalten . 52

4.10 Asia-Call kontinuierlich Euler-Maruyama-Ansatz Lineare Regressions-

analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.11 Asia-Call kontinuierlich Ito-Ansatz Konvergenzverhalten . . . . . . . . 54

4.12 Asia-Call kontinuierlich Ito-Ansatz Lineare Regressionsanalyse . . . . . 54

III

Danksagung

Sehr bedanken mochte ich mich an dieser Stelle bei Prof. Dr. Thomas Gerstner, der

mich beispielhaft betreute und mir auf der Suche nach Losungswegen immer Rede und

Antwort stand. Auch mochte ich mich bei Stefan Heinz bedanken, der bezuglich dieses

Themas ein ausgezeichneter Diskussionspartner war und ebenso wie Matthias Gartner

beim Korrekturlesen eine große Hilfe darstellte.

Ehrenwortliche Erklarung

Hiermit erklare ich, dass ich die vorliegende Arbeit selbstandig abgefaßt und keine an-

deren Hilfsmittel als die angegebenen benutzt habe. Ich erklare ferner, dass diejenigen

Stellen der Arbeit, die anderen Werken wortlich oder dem Sinn nach entnommen sind,

in jedem einzelnen Fall unter Angabe der Quellen kenntlich gemacht sind.

Niederdorfelden, 14. Juni 2011

Kapitel 1

Einleitung

Finanzderivate haben sich seit der Entwicklung des arbitragefreien Bewertungsmodells,

und im Jahre 1973 durch Fischer Black und Myron Scholes sowie unabhangig davon

Robert Merton zu einem wichtigen Bestandteil der modernen Finanzmarkte entwickelt.

Bereits im Jahre 1630 wurden in den Niederlanden Optionen auf Tulpen ausgegeben,

mit denen sich die Tulpenhandler gegen schwankende Preise absichern wollten. Aller-

dings eroffnete erst im Jahre 1973 der erste amtlich geregelte Derivatehandel Board

Option Exchange in Chicago, 1990 entstand in Frankfurt die Deutsche Terminborse.

Im Wesentlichen gibt es drei Klassen von Finanzderivaten

• Optionen: Eine Option ist ein Vertrag, der dem Halter das Recht, aber nicht die

Pflicht gibt, eine bestimmte Transaktion am oder bis zum Verfallstag zu einem

bestimmten Preis, dem Ausubungspreis, zu tatigen.

• Forwards und Futures : Ein Forward ist ein Vertrag, eine Vereinbarung zwischen

zwei Institutionen oder Personen, einen Basiswert untereinander am Verfallstag zu

einem bestimmten Preis zu kaufen oder zu verkaufen. Der Unterschied zur Option

ist, dass der Basiswert bei einem Forward geliefert und bezahlt werden muss.

Ein Future ist im Wesentlichen ein standardisierter Forward, mit dem gehandelt

werden kann. Der Unterschied zum Forward ist, dass der Future-Wert taglich

berechnet und von den Vertragsparteien ausgeglichen wird.

• Swaps : Ein Swap ist ein Vertrag, zu festgelegten Zeitpunkten gewisse finanzielle

Transaktionen zu tatigen, die durch eine vorgegebene Formel bestimmt werden.

Beispiele sind Zinsswaps, die Vereinbarungen zwischen zwei Parteien darstellen,

1

2 KAPITEL 1. EINLEITUNG

Zinszahlungen fur einen bestimmten Betrag innerhalb eines bestimmten Zeitrau-

mes zu leisten.

Bisher gibt es nur fur wenige Derivatetypen analytische Losungen oder Approximations-

Formeln, aus diesem Grund muss auf numerische Verfahren zur Approximation des

Preises zuruckgegriffen werden. An solche Approximationsverfahren ergeben sich aus

der Praxis folgende Anforderungen

1. Konvergenz: Der Berechnungsaufwand sollte moglichst gering sein und bei einem

Einsatz relativ genaue Ergebnisse innerhalb kurzer Zeit liefern.

2. Approximationsfehler: Eine Abschatzung des Approximationsfehlers ist wichtig,

damit das Risiko der Verwendung des Approximationswerts in weiteren Berech-

nungen quantifiziert werden kann.

3. Flexibilitat: In der Praxis werden standig neue Derivatetypen entwickelt. Ein

flexibles Bewertungsverfahren reduziert den Aufwand und ermoglicht schneller

auf Anderungen zu reagieren.

Zur Approximation der Preise von Finanzderivaten existieren verschiedene Ansatze,

z.B. gitterbasierte Verfahren die Binomialbaum-Methode, Finitedifferenzen-Verfahren

usw. In der vorliegenden Arbeit beschaftige ich mich mit einem simulationsbasierten

Verfahren, das Monte-Carlo- bzw Quasi-Monte-Carlo-Verfahren. Der Grund fur die

Wahl von simulationsbasierten Verfahren ist deren Flexibilitat bezuglich immer kom-

plexer werdender Finanzinstrumente, und deren einfache Implementierung.

Ein Problem des Monte-Carlo-Verfahrens ist die geringe Konvergenzgeschwindigkeit

von O(1/√N). Mit Varianzreduktionsmaßnahmen wird die Konvergenzgeschwindigkeit

des Monte-Carlo-Verfahrens um einen konstanten Faktor verbessert. Zwei gute Metho-

den sind Antithetic Variates und Control Variates, auf die ich in dieser Arbeit nicht

eingehen werde. Fur diese Methoden verweise ich den Leser auf P. Glassermann Monte

Carlo Methods in Financial Engineering.

Erweiterung zum Monte-Carlo-Verfahren ist die Quasi-Monte-Carlo-Verfahren. Das Quasi-

Monte-Carlo Verfahren arbeitet nach dem Schema des Monte-Carlo Verfahrens, verwen-

det dabei aber deterministische Niederdiskrepanz-Folgen anstelle von Pseudo-Zufallsfolgen

zur Simulation der Basiswerte. Durch die gleichmaßigere Verteilung der Folgenglieder

im Wertebereich erreicht das Quasi-Monte-Carlo Verfahren ein asymptotisch besseres

Konvergenzverhalten als das Monte-Carlo-Verfahren, die in dieser Arbeit genauer be-

trachtet wird.

1.1. AUFBAU DER ARBEIT 3

1.1 Aufbau der Arbeit

Die vorliegende Arbeit folgt in ihrem Aufbau den bereits dargestellten Punkten

• Im Kapitel 2 werden die allgemeinen Grundlagen und Prinzipien, wie Marktan-

nahmen, Optionstypen und die arbitragefreie Options Bewertungsmethodik des

verbreiteten Black-Scholes-Modells vorgestellt. Das Black-Scholes-Modell ist fur

deterministische Portfolios eine modellspezifische stochastische Differentialglei-

chung und ihre Losung ist ein stochastisches Optimierungsproblem. Zur losung

des stochastischen Optimierungsproblems wird die Stochastik und ihre numeri-

sche Behandlung zum Schluss des Kapitels eingefuhrt.

• Kapitel 3 beschaftigt sich ausschließlich mit dem Quasi-Monte-Carlo-Verfahren.

In einem ersten Schritt werden die Niederdiskrepanz-Folgen vorgestellt, darauf

aufbauend werden die theoretische Grundlagen und die Eigenschaften des Quasi-

Monte-Carlo-Verfahrens diskutiert. Zum Abschluss des Kapitels werden das Quasi-

Monte-Carlo-Verfahren fur die Approximation des Werts eines Integrals als prak-

tische Anwendung vorgestellt, und die Ergebnisse mit dem klassischen Monte-

Carlo-Verfahren hinsichtlich Fehler und Konvergenzverhalten verglichen.

• In Kapitel 4 werden die numerischen Ergebnisse vorgestellt. In einem ersten

Schritt werden die Implementierung des Quasi-Monte-Carlo-Verfahrens vorge-

stellt, daraufhin werden mit dem Verfahren Pfadunabhangige und Pfadabhangige

Optionen bewertet und Fehler und Konvergenzverhalten des Verfahrens disku-

tiert.

• Die Arbeit schließt in Kapitel 5 mit einer Zusammenfassung der dargestellten

Verfahren und einem Ausblick auf mogliche Erweiterungen.

Kapitel 2

Grundlagen

In diesem Kapitel werden die allgemeinen Grundlagen und Prinzipien eingefuhrt, die

spater in dieser Arbeit von Bedeutung sind.

2.1 Marktannahmen

• Auf dem Markt gibt es keine Reibungsverluste, Transaktionskosten, Steuern usw.

• Auf dem Markt gibt es keine Arbitrage1-Moglichkeiten.

• Marktteilnehmer bevorzugen bei gleichen Kosten und Risiken immer die Strategie

mit dem großten Gewinn.

• Marktteilnehmer konnen beliebige Stuckzahlen von Finanzinstrumenten kaufen

und verkaufen, der Preis des Instruments wird dadurch nicht beeinflusst.

• Der Preis eines Instruments ist unabhangig von der Bonitat der Handelspartner.

• Aktien sind in beliebigen Mengen handelbar, kontinuierliches Handeln ist moglich.

• Der Optionshandel ist zu jedem Zeitpunkt t ∈ [0, T ] moglich .

2.2 Optionstypen

Im wesentlichen gibt es zwei Klassen von Optionstypen, Pfadunabhangige Optionen

und Pfadabhangige Optionen.

1Arbitrage: Ein Handelsgeschaft bei dem Preisdifferenzen fur ertrags- und risikoidentische Anlagenausgenutzt werden, um einen risikolosen Gewinn ohne Kapitaleinsatz zu erwirtschaften.

5

6 KAPITEL 2. GRUNDLAGEN

2.2.1 Pfadunabhangige Optionen

Definition 2.1 (Europaische Optionen) Eine Europaische Option ist ein Vertrag,

der einen Basiswert ST am Verfallstag T zum Ausubungspreis K getatigt werden kann.

Die Auszahlungsfunktion eines Europaischen Calls zum Zeitpunkt T lautet

CT = (ST −K)+.

Die Auszahlungsfunktion eines Europaischen Puts zum Zeitpunkt T lautet

PT = (K − ST )+.

Definition 2.2 (Binare Optionen) Eine Binare Option ist ein Vertrag, der wertlos

wird, wenn der Kurs des Basiswerts ST zum Verfallstag T eine festgelegte Schranke K

uber- oder unterschreitet. Im Gegensatz zur Europaischen Optionen ist die Hohe des

Auszahlungsbetrags B unabhangig vom Kurs des Basiswerts. Die Auszahlungsfunktion

eines Binaren Calls zum Zeitpunkt T lautet

CT =

B falls ST > K

0 sonst.

Die Auszahlungsfunktion eines binaren Puts zum Zeitpunkt T lautet

PT =

B falls ST < K

0 sonst.

Definition 2.3 (Chooser Optionen) Eine Chooser Option ist ein Vertrag, der den

Kaufer zum Verfallstag T wahlt lasst, ob er einen Europaischen Call CT oder einen

Europaischen Put PT erhalten mochte. Die Auszahlungsfunktion einer Choose Option

zum Zeitpunkt T lautet

VT = maxCT , PT.

2.2.2 Pfadabhangige Optionen

Definition 2.4 (Asiatische Optionen mit diskretem geometrischen Mittel) Eine

Asiatische Option ist ein Vertrag, dessen Wert von dem Durchschnitts-Kurs des Basis-

werts abhangt. Die Auszahlungsfunktion eines Asiatischen Calls fur M Mittelungszeit-

2.2. OPTIONSTYPEN 7

punkte mit diskretem geometrischen Mittel zum Zeitpunkt T lautet

CT =

( M∏i=1

Sti

) 1M

−K

+

.

Die Auszahlungsfunktion eines Asiatischen Puts fur M Mittelungszeitpunkte mit dis-

kretem geometrischen Mittel zum Zeitpunkt T lautet

PT =

K −( M∏i=1

Sti

) 1M

+

.

Definition 2.5 (Asiatische Optionen mit kontinuierlichem geometrischen Mittel)

Eine Asiatische Option ist ein Vertrag, dessen Wert von dem Durchschnitts Kurs des

Basiswerts abhangt. Wenn die Zahl der Mittelungszeitpunkte sehr groß wird, konnen

stattdessen auch entsprechende kontinuierliche Mittel betrachtet werden. Die Auszah-

lungsfunktion eines Asiatischen Calls fur das kontinuierliche geometrische Mittel zum

Zeitpunkt T lautet

CT =

(exp

(1

T

∫ T

0

ln(Sτ ) dτ

)−K

)+

.

Die Auszahlungsfunktion eines Asiatischen Puts fur das kontinuierliche geometrische

Mittel zum Zeitpunkt T lautet

PT =

(K − exp

(1

T

∫ T

0

ln(Sτ ) dτ

))+

.

Definition 2.6 (Amerikanische Optionen) Eine Amerikanische Option ist ein Ver-

trag, der fur einen Zeitraum 0 ≤ t ≤ T gilt einen Basiswert St bis spatestens zum

Verfallstag T zum Ausubungspreis K getatigt werden kann. Die Auszahlungsfunktion

eines Amerikanischen Calls zum Zeitpunkt t lautet

Ct = max0≤t≤T

(St −K)+.

Die Auszahlungsfunktion eines Amerikanischen Puts zum Zeitpunkt t lautet

Pt = max0≤t≤T

(K − St)+.


2.3 Options-Bewertungsmethodik

2.3.1 Risikoneutraler Ansatz

Den Risikoneutralen Ansatz haben Cox, Ross und Rubinstein entwickelt. Sie beobachten

den Optionspreis nicht direkt von den Risikopraferenzen der einzelnen Marktteilnehmer,

sondern von einem Modell, der von Aktie und risikolosen Bonds abhangt.

Definition 2.7 (Aquivalentes Martingalmaß) Das aquivalente Martingalmaß P ∗

zu der Wahrscheinlichkeitsverteilung P von St ist dasjenige Wahrscheinlichkeitsmaß

unter dem der diskontierte Prozess e−rt · St ein Martingal ist

e−r∆t · E∗[St+∆t] = St fur alle t, ∆t > 0,

wobei E∗ den Erwartungswert und r die risikolose Zinsrate.

Das heißt der kunftige Gewinn der Aktie ist der risikoneutrale Zins.

Satz 2.1 (Martingal-Ansatz fur Europaische Optionen) Der faire Wert einer Op-

tion ohne vorzeitiges Ausubungsrecht ist der diskontierte Erwartungswert der Auszah-

lung unter dem aquivalenten Martingalmaß gegeben durch

V (S0) = e−rT · E∗[V (ST )],

wobei V (ST ) den Wert der Option zum Verfallstag T bezeichnet.

Beweis: Siehe Harrison, Pliska Selbstfinanzierung, Girsanov-Theorem.

Satz 2.2 (Martingal-Ansatz fur Amerikanische Optionen) Der faire Wert einer

Option mit vorzeitigem Ausubungsrecht ist ein optimale-Stoppzeit-Problem der diskon-

tierte Erwartungswert der Auszahlung unter dem aquivalenten Martingalmaß gegeben

durch

V (S0) = max0≤t≤T

e−rt · E∗[V (St)],

wobei V (St) den Wert der Option zum Zeitpunkt t ∈ [0, T ] bezeichnet.

Beweis: Siehe Harrison, Pliska Selbstfinanzierung, Girsanov-Theorem.

2.3. OPTIONS-BEWERTUNGSMETHODIK 9

2.3.2 Black-Scholes-Ansatz

Mit Hilfe der Black-Scholes-Formel lassen sich einige Optionspreise berechnen, wie von

Europaischen Optionen, Binaren Optionen, Asiatischen Optionen mit diskretem geome-

trischen Mittel, Asiatischen Optionen mit kontinuierlichem geometrischen Mittel usw.

Satz 2.3 (Black-Scholes-Formel fur Europaische Call)

C0 = S0Φ(d1)−Ke−rTΦ(d2)

wobei

Φ(x) =1√2π

∫ x

−∞e−s

2/2 ds, x ∈ R,

d1 =ln(S0/K) + (r + σ2/2)T

σ√T

,

d2 = d1 − σ√T .

Satz 2.4 (Black-Scholes-Formel Asiatische Call diskret geometrisch)

VS,0 = S0AΦ(d+ σ

√T1

)−Ke−rTΦ(d)

wobei

Φ(x) =1√2π

∫ x

−∞e−s

2/2 ds, x ∈ R,

A = exp

(−r(T − T2)− σ2(T2 − T1)

2

),

d =ln(S0/K) + (r − σ2/2)T2

σ√T1

,

T1 = T − M(M − 1)(4M + 1)

6M2∆t,

T2 = T − (M − 1)

2∆t.

Satz 2.5 (Black-Scholes-Formel Asiatische Call kontinuierlich geometrisch)

VS,0 = S0e−12

(r+σ2/6)TΦ(d+ σ

√T/3

)−Ke−rTΦ(d)


wobei

Φ(x) =1√2π

∫ x

−∞e−s

2/2ds, x ∈ R,

d =ln(S0/K) + 1

2(r − σ2/2)T

σ√T/3

.

Definition 2.8 (Black-Scholes-Modell) Black und Scholes modellieren die Kurs-

entwicklung St einer Aktie durch eine stochastische Differentialgleichung

dSt = µdt+ σdWt

wobei

St : Kurs der Aktie S zum Zeitpunkt t

µ : Drift

dt : Zeitdifferential

σ : Volatilitat

dWt : Standard-Wiener-Prozess.

2.4 Stochastik und ihre numerische Behandlung

Um die Black-Scholes stochastische Differentialgleichung zu losen, wird die Stochastik

und ihre numerische Behandlung benotigt.

Definition 2.9 (Stochastischer Prozess) Sei (Ω,Σ, P ) ein Wahscheilichkeitsraum,

mit Ereignismenge Ω, Ereignisalgebra Σ und Wahrscheinlichkeitsmaß P auf Σ. Sei

(Z,Z) ein weiterer mit einer Sigma-Algebra versehner Raum und T eine Indexmenge.

Ein stochastischer Prozess ist eine Familie von Zufallsvariablen

Xt : Ω→ Z

fur alle t ∈ T , wobei die Abbildung Xt Z-messbar sein muss.

Definition 2.10 (Wiener-Prozess) Ein Wiener-Prozess ist ein zeitstetiger stochas-

tischer Prozess, der normalverteilte unabhangige Zuwachse hat und die folgenden Ei-

genschaften

• W0 = 0 fast sicher

2.4. STOCHASTIK UND IHRE NUMERISCHE BEHANDLUNG 11

• Wt ist stetig fast sicher

• Wt −Ws ∼√t− s · Φ(0, 1) = Φ(0, t− s)

wobei Φ(0, 1) eine normalverteilte Zufallszahl mit Mittelwert 0 und Varianz 1 ist.

Definition 2.11 (Euler-Maruyama-Verfahren) Sei (Wt)t≥0 ein Wiener-Prozess und

a, b : R×[0, T ]→ R zwei Funktionen mit folgendem stochastischem Anfangswertproblem

dSt = a(t, St)dt+ b(t, St)dWt, St0 = A.

Diese wird als Integralgleichung

St = St0 +

∫ t

t0

a(u, Su) du+

∫ t

t0

b(u, Su) dWu, t ≥ tt0 ,

interpretiert, wobei das zweite Integral ein Ito-Integral ist. Die stochastische Differenti-

algleichung wird diskretisiert,

Sti+1= Sti + a(ti, Sti)h+ b(ti, Sti)∆Wi, i=0,...,N−1,

wobei h = T/N die Schrittweite zu N ∈ N auf dem Gitter ti = i · h, i=0,...,N−1, ist und

∆Wi+1 = Wti+1−Wti, i=0,...,N−1.

Das wichtigste theoretische Resultat der Naherung St, bezuglich des exakten Wertes St

des Euler-Maruyama-Verfahrens, wird durch die Konvergenz des Verfahrens erfasst.

Definition 2.12 (Starke und schwache Konvergenz) Ein numerisches Verfahren

zur Losung einer stochastischen Differentialgleichung konvergiert stark mit der Ordnung

α > 0, falls fur alle p ∈ N

maxn

E[|Stn − Stn|p

] 1p ≤ K1 · hα

gilt mit einer von der Schrittweite h unabhangigen Konstante K1 > 0.

Das Verfahren konvergiert schwach mit der Ordnung β > 0, falls fur alle Polynome P

maxn

∣∣∣E[P (Stn)]− E[P (Stn)]∣∣∣ ≤ K2 · hβ

gilt mit einer von der Schrittweite h unabhangigen Konstanten K2 > 0.


Satz 2.6 (Konvergenz des Euler-Maruyama-Verfahrens) Gilt fur die Funktio-

nen a(t, S) und b(t, S) mit Konstanten K1, K2, K3, K4 > 0

1. |a(t, x) + a(t, y)|+ |b(t, x) + b(t, y)| ≤ K1|x− y|,

2. |a(t, x)|+ |b(t, x)| ≤ K2(1 + |x|),

3. |a(s, x)− a(t, x)|+ |b(s, x)− b(t, x)| ≤ K3(1 + |x|)|s− t|1/2,

4. a, b ∈ C(2+ε) fur ein ε > 0 (mehr als 2-mal differenzierbar),

dann existiert eine eindeutige Losung der stochastischen Differentialgleichung und das

Euler-Maruyama-Verfahren konvergiert stark mit der Ordnung 1/2 und schwach mit

der Ordnung 1.

Beweis: Siehe P. Kloeden Numerical Solution of Stochastic Differential Equations.

Satz 2.7 (Ito-Lemma) Sei St ein Ito-Prozess, a und b zwei Funktionen mit

dSt = a(t, St)dt+ b(t, St)dWt

und V (t, St) eine Funktion mit stetigen Ableitungen ∂V∂t

, ∂V∂St

und ∂2V∂S2

t. Dann ist V (t, St)

ein Ito-Prozess mit

dVt =

(∂V

∂t+ a(t, St)

∂V

∂St+

1

2(b(t, St))

2∂2V

∂S2t

)dt+ b(t, St)

∂V

∂StdWt

Beweis: Taylor-Entwicklung von V (t, St) in t und St

dV =∂V

∂tdt+

∂V

∂StdSt +

1

2

∂2V

∂S2t

dS2t + . . . .

ersetzen dSt = adt+ bdWt und erhalten

dV =∂V

∂tdt+

∂V

∂St(adt+ bdWt) +

1

2

∂2V

∂S2t

(a2dt2 + 2abdtdWt + b2dW 2t ) + . . . .

Fur dt → 0 konnen dt2 und dtdWt vernachlassigt werden und dW 2t strebt nach dt,

denn dW 2t → E[dW 2

t ] = dt. Nach Zusammenfassen der ubrigen Terme, ist das Lemma

bewiesen.

2.4. STOCHASTIK UND IHRE NUMERISCHE BEHANDLUNG 13

Bemerkung: Im Black-Scholes-Modell ist a(S, t) = µSt und b(S, t) = σSt.

Satz 2.8 (Ito-Formel) Die Losung der Black-Scholes-SDE ist gegeben durch

St = S0 · e(µ− 12σ2)t+σWt .

Diskretisierung

Sti+1= Sti · e(µ− 1

2σ2)h+σ

√h∆Wi ,

wobei h = T/N die Schrittweite zu N ∈ N auf dem Gitter ti = i · h, i=0,...,N−1, ist und

∆Wi+1 = Wti+1−Wti, i=0,...,N−1.

Beweis: Setze V (t, St) = ln(StS0

), dann ist

d ln

(StS0

)= 0 +

(Stµ

1

St+

1

2S2t σ

2

(− 1

S2t

)dt+ Stσ

1

StdWt

)

=

(µ− σ2

2

)dt+ σdWt.

Daraus folgt

ln

(StS0

)=

(µ− σ2

2

)t+ σWt

⇔ St = S0 · e(µ− 12σ2)t+σWt ,

wegen limt→0

St = S0 folgt die Behauptung.

Bemerkung: Der logarithmische Zuwachs ln(StS0

)ist normalverteilt mit Mittelwert(

µ− σ2

2

)t und Varianz σ2t, und somit ist St lognormalverteilt.

Es gilt E[eσWt ] = eσ2t2 , folglich ist der Erwartungswert von St

E[St] = S0 · eµt

und die Varianz von St

V ar[St] = E[S2t ]− (E[St])

2


= S20 · e(2µ+σ2)t − (S0eµt)2

= S0 · e2µt(eσ2t − 1).

Kapitel 3

Quasi-Monte-Carlo-Verfahren

Das Monte-Carlo-Verfahren kann fur Approximation des Wertes eines Integrals verwen-

det werden. Sei Ω ein zusammenhangendes Gebiet aus Rd. Die Mote-Carlo-Approximation

ist definiert als

I(f) =

∫Ω

f(u) du ≈ V ol(Ω) · 1

N

N∑i=1

f(xi),

also erfolgt zunachst die Auswertung der zu integrierenden Funktion f an zufallig

gewahlten Stellen xi fur i=1,...,N im Integrationsbereich, und das arithmetische Mittel

der Funktionswerte f(xi) approximiert den Wert des Integrals. Nach dem Gesetz der

großen Zahlen konvergiert das arithmetische Mittel der Funktionswerte mit zunehmen-

der Anzahl von Auswertungen gegen den wahren Wert des Integrals. Das Monte-Carlo-

Verfahren benutzt bei der Approximation die Pseudo-Zufallsfolgen. Das Wesentliche

ist die Zufalligkeit des Integrationsbereichs. Ein Nachteil ist das langsames Konver-

genzverhalten. Um dieses Verhalten zu verbessern, bietet sich das Quasi-Monte-Carlo-

Verfahren an. Das Quasi-Monte-Carlo-Verfahren verwenden statt Pseudo-Zufallsfolgen

Niederdiskrepanz-Folgen zur Bestimmung der Auswertungsstellen. Die Glieder dieser

Folgen sollen den Integrationsbereich gleichmaßiger ausfullen und dadurch eine schnel-

lere Konvergenz des Verfahrens erzielen.

15

16 KAPITEL 3. QUASI-MONTE-CARLO-VERFAHREN

3.1 Niederdiskrepanz-Folgen

3.1.1 Halton-Folge

Die Halton-Folge ist eine der am leichtesten zu berechnenden Niederdiskrepanz-Folgen.

Zur Berechnung des j-ten Elements des i-ten Glieds x(j)i wird die Zahl i in einer Basis

pj dargestellt mit

i =∞∑k=0

nkpkj .

Die benutzten Basen sind die ersten d Primzahlen (d ist Dimension der xi fur i=1,...,N).

Aus der Entwicklung von i bezuglich pj erhalt man x(j)i durch die radikal-inverse Trans-

formation

x(j)i =

∞∑k=0

nkp−k−1j .

Abbildung 3.1 Halton-2d

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Halton−Folge (N = 1000)

p1 = 2

p2 =

3

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p1 = 2

p2 =

3

In Abbildung 3.1 zeigt die linke Grafik die ersten 1000 und die rechte Grafik die ersten

10000-Punkte aus der 1. und 2. Komponente einer 2-dimensionalen Halton-Folge.

Die Halton-Folge verwendet fur jede Dimension unterschiedliche Basen, da andernfalls

einzelne Komponenten in jedem Glied gleich waren und bestimmte Teile des Integrati-

onsbereichs nicht abgedecken wurden. Die kleinsten moglichen Basen werden benutzt,

um ein sehr uniformes Verhalten auch bei wenigen Folgengliedern zu erreichen, je grosser

die beteiligten Basen sind, desto mehr Folgenglieder werden benotigt, um dem Raum

3.1. NIEDERDISKREPANZ-FOLGEN 17

zu fullen. Benutzt man die Basen p1 < p2 < . . . < pN , dann werden p1 · p2 · · · pN Fol-

genglieder benotigt, um den Raum in den Dimensionen zu fullen, sodass kein großerer

Bereich leer bleibt.

Fur d = 10 werden Beispielsweise (2 · 3 · 5 · 7 · 11 · 13 · 17 · 19 · 23 · 29) = 6469693230

Folgen benotigt. Mit steigender Dimension, wachst die Basen stark auf, so tritt eine

gleichmassige Verteilung erst fur eine sehr grosse Anzahl von Folgen auf. In Abbildung

3.2 wird das Verhalten illustriert.

Abbildung 3.2 Halton-100d

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p1 = 523

p2 =

541

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p1 = 523

p2 =

541


10000-Punkte aus der 99. und 100. Komponente einer 100-dimensionalen Halton-Folge.

Hier erkennt man, dass großen Bereiche leer bleiben. Um dieses Verhalten zu verbessern,

wird die Faure-Folge vorgestellt.

3.1.2 Faure-Folge

Die Faure-Folge beruht wie die Halton-Folge auf der Entwicklung der Folgennummer i

in einer Basis p

i =∞∑k=0

nkpk.

Faure-Folgen verwenden fur alle Komponenten ein einheitliches p, die kleinste Primzahl

mit p ≥ d. Die einzelnen Zahlen nk werden fur jede Dimension j in Ziffern α(j)k wie folgt


geordnet

α(j)k =

[∞∑m≥k

(m

k

)nk(j − 1)m−k

]mod p.

Nach Umformung der Ziffern berechnen sich die Folgekomponenten x(j)i wie bei der

Halton-Folge durch Spiegelung der α(j)k am Dezimalpunkt

x(j)i =

∞∑k=0

α(j)k p−k−1.

Abbildung 3.3 Faure-2d

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Faure−Folge (N = 10000)

p1 = 2

p2 =

2

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p1 = 2

p2 =

2


10000 Punkte der 1. und 2. Komponente einer 2-dimensionalen Faure-Folge.

Im Unterschied zu Halton-Folge fullt die Faure-Folge den d-dimensionalen Einheitswurfel

in einheitlichen Zyklen der Lange p. Bei gleicher Dimension ist p deutlich kleiner als

die großte Basis in der Halton-Folge. Diese wirkt sich auf die Uniformitat in hoheren

Dimension aus. Das Verhalten wird in Abbildung 3.4 illustriert.

3.1. NIEDERDISKREPANZ-FOLGEN 19

Abbildung 3.4 Faure-100d

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p99 = 101

p100

= 1

01

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p99 = 101

p100

= 1

01


10000 Punkte aus der 99. und 100. Komponente einer 100-dimensionalen Faure-Folge.

Hier erkennt man, dass fur wenig Simulationen immer noch großen Bereiche leer bleiben.

Um dieses Verhalten zu verbessern, wird die Sobol-Folge vorgestellt.

3.1.3 Sobol-Folge

Die Sobol-Folge ist eine der am schwierigsten zu generierenden Folgen. Zuerst generiert

man fur jede Komponente j ein primitives Polynom uber F2 der Form

Pj(x) = xdj + α(j)1 xdj−1 + α

(j)2 xdj−2 + . . .+ α

(j)dj−1x+ 1.

Die berechneten α(j)k verwendet man anschliessend zur Generierung von Mengen unge-

rader naturlicher Zahlen

M (j) = m(j)1 , . . . ,m

(j)dlog2Ne.

Hierbei ist 0 < m(j)i < 2i und N bezeichnet die Anzahl der zu generierenden Glieder.

Fur i ≤ dj wahlt man zur Generierung die d kleinsten ungeraden ganzen Zahlen aus.

Fur dj > i erhalt man m(j)i aus den α

(j)k durch die rekursive Formel

m(j)i = 2α

(j)1 m

(j)i−1 ⊕ 22α

(j)2 m

(j)i−2 ⊕ · · · ⊕ 2dj−1α

(j)dj−1m

(j)i−d+1 ⊕ 2djm

(j)i−dj ⊕m

(j)i−dj .


Der Operator ⊕ bezeichnet das bitweise exklusive Oder (XOR).

Beispielsweise ist

1001101⊕ 0101001 = 1100100.

Durch Umwandlung der mji in binare Bruche erhalt man sogenannte Direction-Numbers

vji

v(j)i =

m(j)i

2ifur i = 1, 2, . . . , dlog2Ne.

Zur Generierung der einzelnen Folgenglieder xk geht man fur jede Komponente x(j)k wie

folgt vor

x(j)k−1 =

∞∑i=0

ni2−i−1.

Man bestimmt nun in einer Binardarstellung von k− 1 die Position l der am weitesten

rechts liegenden 0-Ziffer. Anschließend wird die zu l und j zugehorige Direction-Number

binar dargestellt durch

v(j)l =

∞∑k=0

ok2−k−1.

Die Ergebniskomponente x(j)k erhalt man durch XOR aus den Ziffern von x

(j)k−1 und v

(j)l

x(j)k =

∞∑i=0

(ni ⊕ oi)2−i−1.

Abbildung 3.5 Sobol-2d

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1Sobol−Folge (N = 1000)

p1 = 2

p2 =

2

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p1 = 2

p2 =

2

3.2. THEORETISCHE EIGENSCHAFT 21


10000 Punkte aus der 1. und 2. Komponente einer 2-dimensionalen Sobol-Folge.

Die Sobol-Folge verwendet unabhangig von der Dimension die einheitliche Basis 2.

Durch die geringere Lange des Zyklus sollten sich Probleme bei der gleichmaßigen

Ausfullung des Raumes nicht so stark bemerkbar machen wie bei den anderen Folgen.

Abbildung 3.6 Sobol-100d

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p99 = 2

p100

= 2

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9


p99 = 2

p100

= 2


10000 Punkte aus der 99. und 100. Komponente einer 100-dimensionalen Sobol-Folge.

Hier erkennt man, dass die Sobol-Folgen auch bei wenigen Simulationen eine gute

gleichmaßige Auffullung des Raumes aufweißt.

3.2 Theoretische Eigenschaft

Die wesentliche Forderung an die verwendeten Folgen ist die gleichmaßige Verteilung

der Folgenglieder.

Definition 3.1 Sei Xd = (x1, . . . , xN) eine endliche Folge von Punkten, wobei jedes

xi fur i=1,...,N ein Punkt aus [0, 1]d ist. Die Quasi-Monte-Carlo-Approximation mit In-

tegrationsbereich [0, 1]d ist definiert als

∫[0,1]d

f(u) du ≈ 1

N

N∑i=1

f(xi).


Im Idealfall wird die endliche Folge Xd durch eine unendliche Folge Xd∗ = (x1, x2, . . .)

von Punkten aus [0, 1]d ersetzt. Die Mindestanforderung an Xd∗ ist das Konvergenzver-

halten.

Definition 3.2 Sei Xd∗ = (x1, x2, . . .) eine unendliche Folge von Punkten, wobei jedes

xi fur i≥1 ein Punkt aus [0, 1]d ist. Xd∗ heißt gleichverteilt auf [0, 1]d, wenn gilt

limN→∞

1

N

N∑i=1

f(xi) =

∫[0,1]d

f(u) du

fur alle stetigen Funktionen f ∈ [0, 1]d.

Definition 3.3 Sei Xd∗ = (x1, x2, . . .) eine unendliche Folge von Punkte, wobei jedes

xi fur i≥1 ein Punkt aus [0, 1]d ist. Xd∗ heißt gleichverteilt auf [0, 1]d, wenn gilt

limN→∞

1

N

N∑i=1

CM(xi) = λd(M)

fur alle Teilmenge M aus [0, 1]d, wobei CM die charakteristische Funktion von M und

λd das d-dimensionale Lebesgue-Maß ist.

3.2.1 Diskrepanz

Die Diskrepanz ist ein Maß fur die Gleichverteilung.

Definition 3.4 (Zahlfunktion) Sei Xd = (x1, . . . , xN) eine endliche Folge von Punk-

ten, wobei jedes xi fur i=1,...,N ein Punkt aus [0, 1]d ist. Außerdem seiM eine Teilmenge

aus [0, 1]d und CM die charakteristischen Funktion von M.

A(M;Xd) =N∑i=1

CM(xi)

heisst Zahlfunktion, die die Anzahl der xi ∈M fur i=1,...,N zahlt.

Die eindimensionale Diskrepanz

Definition 3.5 (Diskrepanz) Sei X = (x1, . . . , xN) eine endliche Folge von Punkte,

wobei jedes xi fur i=1,...,N ein Punkt aus [0, 1] ist. Außerdem seien a, b zwei Punkte aus


[0, 1] mit 0 ≤ a < b ≤ 1. Die Diskrepanz von der Folge X ist definiert als

DN(X) = sup0≤a<b≤1

∣∣∣∣A((a, b);X)

N− (b− a)

∣∣∣∣ .Definition 3.6 (Sterndiskrepanz) Sei X = (x1, . . . , xN) eine endliche Folge von

Punkte, wobei jedes xi fur i=1,...,N ein Punkt aus [0, 1] ist. Außerdem sei b ein Punkt aus

[0, 1] mit 0 ≤ b ≤ 1. Die Sterndiskrepanz von der Folge X ist definiert als

D∗N(X) = sup0≤b≤1

∣∣∣∣A((0, b);X)

N− b∣∣∣∣ .

Satz 3.1 Die Folge X ist gleichverteilt auf [0, 1] genau dann, wenn

limN→∞

DN(X) = 0.

Beweis: Wahle m ≥ 2. Fur 0 ≤ k ≤ m − 1 sei Ik = ( km, k+1m

). Da X gleichverteilt ist,

existiert ein N0 = N0(m) ∈ N, sodass fur N > N0 und jedes k = 0, 1, . . . ,m− 1 gilt

1

m(1− 1

m) ≤ A(Ik;X)

N≤ 1

m(1 +

1

m). (3.1)

Betrachte nun beliebige halboffene Teilintervalle J = [α, β) aus [0, 1]. Es existieren

Intervalle J1 und J2, die endliche Vereinigungen der Ik sind, mit J1 ⊆ J ⊆ J2, λ(J) −λ(J1) < 2

mund λ(J2)− λ(J) < 2

m. Aus (3.1) folgt fur alle N > N0:

λ(J1)

(1− 1

m

)≤ A(J1;X)

N≤ A(J ;X)

N≤ A(J2;X)

N≤ λ(J2)

(1 +

1

m

)und somit (

λ(J)− 2

m

)(1− 1

m

)<A(J ;X)

N<

(λ(J) +

2

m

)(1 +

1

m

).

Durch λ(J) ≤ 1 ergibt sich

− 3

m− 2

m2<A(J ;X)

N− λ(J) <

3

m+

2

m2∀ N ≥ N0. (3.2)

Da die Schranken in (3.2) unabhangig von J sind, folgt DN(X) ≤ 3m

+ 2m2 ∀ N ≥ N0.

Es kann 3m

+ 2m2 beliebig klein gemacht werden. Damit folgt die Behauptung.


Satz 3.2 Fur jede Folge X = (x1, . . . , xN) aus [0, 1] gilt

1

N≤ DN(X) ≤ 1.

Beweis: Die rechte Ungleichung folgt direkt aus der Definition.

Fur die linke Ungleichung, wahle ε > 0 und sei x ein belibiege Punkt aus [0, 1]. Betrachte

das Intervall J = [x, x+ ε)∩ [0, 1]. Da x ∈ J ist, gilt A(J ;X)N−λ(J) ≥ 1

N−λ(J) ≥ 1

N− ε.

Daraus folgt DN(X) ≥ 1N− ε. Damit ist die gewunschte Ungleichung gezeigt.

Lemma 3.1 Seien X = (x1, . . . , xN) und Y = (y1, . . . , yN) zwei endliche Folgen von

Punkten, wobei jedes xi und yi fur i=1,...,N Punkte aus [0, 1] sind. Wenn |xi− yi| < ε fur

alle i erfullen, dann gelten folgende Abschatzungen

|DN(x1, . . . , xN)−DN(y1, . . . , yN)| ≤ 2ε,

|D∗N(x1, . . . , xN)−D∗N(y1, . . . , yN)| ≤ ε.

Beweis: Betrachte J = [0, b) ⊆ [0, 1). Wenn yi ∈ J ist, dann ist xi ∈ J1 := [0, b + ε) ∩[0, 1].

Daraus folgt

A(J ;Y )

N− λ1(J) ≤ A(J1;X)

N− λ1(J1) + ε ≤ D∗N(X) + ε.

Wenn xi ∈ J2 := [0, b− ε), dann ist yi ∈ J .

Daraus folgt

A(J ;Y )

N− λ1(J) ≥ A(J2;X)

N− λ1(J2)− ε ≥ −D∗N(X)− ε.

Daher gilt D∗N(Y ) ≤ D∗N(X)+ε. Durch Vertauschen von X und Y folgt D∗N ≤ D∗N(Y )+

ε, woraus |D∗N(X)−D∗N(Y )| ≤ ε folgt.

Die erste Ungleichung wird analog gezeigt.


Satz 3.3 Sei X = (x1, . . . , xN) eine endliche Folge von Punkten, wobei jedes xi fur

i=1,...,N ein Punkt aus [0, 1] ist. Außerdem seien die Punkte in aufsteigender Reihenfolge

sortiert, sodass 0 ≤ x1 ≤ x2 ≤ · · · ≤ xN ≤ 1 sind. Dann gilt

DN(X) =1

N+ max

1≤i≤N

(i

N− xi

)− min

1≤i≤N

(i

N− xi

).

Beweis: Setze x0 := 0 und xN+1 := 1. Dann ist

DN(X) = max1≤i≤j≤N

supxi<a≤xi+1xj<b≤xj+1

a<b

∣∣∣∣A([a, b);X)

N− (b− a)

∣∣∣∣= max

1≤i≤j≤Nsup

xi<a≤xi+1xj<b≤xj+1

a<b

∣∣∣∣j − iN− (b− a)

∣∣∣∣= max

1≤i≤j≤N

(∣∣∣∣j − iN− (xj+1 − xi)

∣∣∣∣ , ∣∣∣∣j − iN− (xj − xi+1)

∣∣∣∣) .Setze ri = i

N− xi fur 0≤i≤N+1. Dann ist

DN(X) = max1≤i≤jN

(∣∣∣∣rj+1 − ri −1

N

∣∣∣∣ , ∣∣∣∣rj − ri+1 +1

N

∣∣∣∣)= max

0≤i≤N1≤j≤N+1

∣∣∣∣ 1

N+ ri − rj

∣∣∣∣ .Wenn der letzte Term auf 1 ≤ i, j ≤ N eingeschrankt wird, dann ist DN(X) eindeutig

festgelegt durch

max1≤i≤N

ri ≥ rN ≥ 0, min1≤i≤N

ri ≤ r1 ≤1

N.

Daraus folgt, dass das Maximum entweder fur i = 0 oder j = N + 1 ist und DN(X)

nicht ubertroffen werden kann.

Satz 3.4 Sei X = (x1, . . . , xN) eine endliche Folge von Punkten, wobei jedes xi fur

i=1,...,N ein Punkt aus [0, 1] ist. Außerdem seien die Punkte in aufsteigender Reihenfolge

sortiert, sodass 0 ≤ x1 ≤ x2 ≤ · · · ≤ xN ≤ 1 sind. Dann gilt

D∗N(X) = D∗N(x1, . . . , xN) =1

2N+ max

1≤i≤N

∣∣∣∣xi − 2i− 1

2N

∣∣∣∣ .


Beweis: Setze x0 = 0 und xN+1 = 1. Dann ist

D∗N(X) = max0≤i≤N

supxi<b≤xi+1

∣∣∣∣A([0, b);X)

N− b∣∣∣∣

= max0≤i≤N

supxi<b≤xi+1

∣∣∣∣ iN − b∣∣∣∣

= max0≤i≤N

max

(∣∣∣∣ iN − xi∣∣∣∣ , ∣∣∣∣ iN − xi+1

∣∣∣∣)= max

1≤i≤Nmax

(∣∣∣∣ iN − xi∣∣∣∣ , ∣∣∣∣i− 1

N− xi

∣∣∣∣)=

1

2N+ max

1≤i≤N

∣∣∣∣xi − 2i− 1

2N

∣∣∣∣ .

Mehrdimensionale Diskrepanz

Definition 3.7 (Mehrdimensionale Diskrepanz) Sei Xd = (x1, . . . , xN) eine end-

liche Folge von Punkten, wobei jedes xi fur i=1,...,N ein Punkt aus [0, 1]d ist. Außerdem

sei Q eine d-dimensionale Quader der Form∏d

i=1[ai, bi] mit 0 ≤ ai ≤ bi ≤ 1 fur alle i.

Dann ist die Diskrepanz von der Folge Xd definiert als

DN(Xd) = supQ⊂[0,1]d

∣∣∣∣A(Q;Xd)

N− λd(Q)

∣∣∣∣ .Definition 3.8 (Mehrdimensionale Sterndiskrepanz) Sei Xd = (x1, . . . , xN) ei-

ne endliche Folge von Punkten, wobei jedes xi fur i=1,...,N ein Punkt aus [0, 1]d ist.

Außerdem sei Q∗ eine d-dimensionale Quader der Form∏d

i=1[0, bi] mit 0 ≤ bi ≤ 1 fur

alle i. Dann ist die Sterndiskrepanz von Xd definiert als

DN(Xd) = supQ∗⊂[0,1]d

∣∣∣∣A(Q∗;Xd)

N− λd(Q∗)

∣∣∣∣ .Proposition 3.1 Fur jede endliche Folge Xd, deren Punkte in [0, 1]d liegen, gilt

D∗N(Xd) ≤ DN(Xd) ≤ 2dD∗N(Xd).

Beweis: Die erste Ungleichung folgt direkt aus der Definition.

Fur die zweite Ungleichung, fur d = 1 betrachte A([0, b);Xd) − A([0, a);Xd) mit 0 ≤


a < b ≤ 1 und λ1([a, b)) = λ1([0, b))− λ1([0, a)). Daraus folgt∣∣∣∣A([a, b);Xd)

N− (b− a)

∣∣∣∣ ≤ ∣∣∣∣A([0, b);Xd)

N− b∣∣∣∣+

∣∣∣∣A([0, a);Xd)

N− a∣∣∣∣ .

Bilde das Supremum daruber, folgt die Behauptung.

Fur d = 2 betrachte

Q = (x1, x2 ∈ [0, 1]2 : a1 ≤ x1 < b1 und a2 ≤ x2 < b2

= [a1, b1)× [a2, b2) mit 0 ≤ ai < bi ≤ 1, fur i = 1, 2.

Daraus folgt

Q = ([0, b1)× [0, b2)) \ ([0, a1)× [0, b2)) \ ([0, b1)× [0, a2)) \ ([0, a1)× [0, a2))

= (Q∗1 \ Q∗2) \ (Q∗3 \ Q∗4).

Es gilt

λ(Q) = λ(Q∗1)− λ(Q∗2)− λ(Q∗3) + λ(Q∗4).

Daraus folgt

A(Q;Xd) = A(Q∗1;Xd)− A(Q∗2;Xd)− A(Q∗3;Xd) + A(Q∗4;Xd).

Durch dasselbe Argument wie fur d = 1 folgt die Behauptung.

Fur d ≥ 3 analog.

Fur allgemeine konvexe Integrationsbereiche mussen neue Arten von Diskrepanzen de-

finiert werden.

Definition 3.9 (Isotrope Diskrepanz) Sei Xd = (x1, . . . , xN) eine endliche Folge

von Punkten, wobei jedes xi fur i=1,...,N ein Punkt aus [0, 1]d ist. Außerdem sei C eine

konvexen Teilmengen aus [0, 1]d. Die isotrope Diskrepanz von der Folge Xd ist definiert

als

TN(Xd) = supC⊂[0,1]d

∣∣∣∣A(C;Xd)

N− λd(C)

∣∣∣∣Proposition 3.2 Fur jede endliche Folge Xd, deren Punkte in [0, 1]d liegen, gilt

DN(Xd) ≤ TN(Xd) ≤ 4dDN(Xd)1/d.


Beweis: Siehe H. Niederreiter Diskrepanz und Distanz von Maßen bezuglich konvexer

und Jordanscher Mengen.

Aus Proposition 3.2 folgt, dass eine endliche Folge Xd genau dann auf [0, 1]d gleichver-

teilt ist, wenn limN→∞

TN(Xd) = 0 ist. Wobei TN(Xd) die isotrope Diskrepanz von Xd ist.

Fur Jordan-messbare Teilmengen von [0, 1]d gibt es eine Unterteilung nach der Kom-

plexitat ihrer Rander. Fur B ⊆ [0, 1]d und ε > 0 definiere

Bε = x ∈ [0, 1]d : ‖x− y‖ ≤ ε fur endlich viele y ∈ B,

B−ε = x ∈ [0, 1]d : ‖x− y‖ ≥ ε fur alle y ∈ [0, 1]d \ B.

Sei b eine nicht fallende Funktion fur alle ε > 0 mit limε↓0

b(ε) = 0. Weiter sei Mb die

Familie aller Lebesgue-messbaren Mengen B ⊆ [0, 1]d fur, die folgendes gilt

λd(Bε \ B) ≤ b(ε) und λd(B \ B−ε) ≤ bε fur alle ε > 0.

Jedes B ∈ Mb ist somit Jordan-messbar und umgekehrt gehort jede Jordan-messbare

Teilmenge von [0, 1]d zuMb fur eine passende Funktion b. Wenn die Funktion b fur alle

ε > 0 die Ungleichung b(ε) ≥ ε erfullt, dann kann wie folgt abgeschatzt werden

DN(Mb;Xd) ≤ 4 · b(2

√dDN(Xd)1/d).

In vielen Fallen hat die Funktion b die Form b(ε) = C fur eine Konstante C > 0. Dann

folgt die Abschatzung

DN(Mb;Xd) ≤ (4C

√d+ 2C + 1)DN(Xd)1/d.

Fur jede konvexe Teilmenge von [0, 1]d zu Mb0 mit b0 = 2dε, folgt

TN(Xd) ≤ DN(Mb, Xd).

Wenn die Funktion b die Ungleichung b(ε) ≥ 2dε fur alle ε > 0 erfullt und die Folge

Xd gleichverteilt ist auf [0, 1]d, dann gilt limN→∞

DN(Mb, Xd) = 0. Unter gewissen Vor-

aussetzungen an die Funktion b gilt auch die Umkehrung, dass limN→∞

DN(Mb, Xd) = 0

impliziert, dass Xd gleichverteilt ist auf [0, 1]d.


Bemerkung: Die folgende Aussagen sind aquivalent

• Xd ist gleichverteilt auf [0, 1]d,

• limN→∞

DN(Xd) = 0 und

• limN→∞

D∗N(Xd) = 0.

In diesem Sinn kann die Diskrepanz und die Sterndiskrepanz als Quantifizierung von

gleichverteilten Folgen in [0, 1]d gesehen werden.

3.2.2 Fehlerschranken

Jetzt kann man mit Hilfe der Diskrepanz Fehlerschranke fur die Quasi-Monte-Carlo-

Approximation angeben werden.

Die eindimensionale Fehlerschranken

Satz 3.5 (Koksma-Ungleichung) Besitzt f eine beschrankte Variation V (f) auf [0, 1],

dann gilt fur alle endlichen Folgen X = (x1, . . . , xN), deren Punkte aus [0, 1] sind,∣∣∣∣∣ 1

N

N∑i=1

f(xi)−∫ 1

0

f(u) du

∣∣∣∣∣ ≤ V (f) ·D∗N(X).

Beweis: Sei x1 ≤ x2 ≤ · · · ≤ xN , setze x0 := 0 und xN+1 := 1. Durch die Verwendung

von partieller Summation und partieller Integration folgt

1

N

N∑i=1

f(xi)−∫ 1

0

f(u) du = −N∑i=0

i

N(f(xi+1)− f(xi)) +

∫ 1

0

u df(u)

=N∑i=0

∫ xi+1

xi

(x− i

N

)df(u).

Fur feste i fur i=1,...,N folgt∣∣∣∣x− i

N

∣∣∣∣ ≤ D∗N(X) fur xi ≤ x ≤ xi+1.

Durch Satz 3.4 folgt die Behauptung.


Definition 3.10 Das Stetigkeitsmaß fur eine stetige Funktion f auf [0, 1] ist definiert

als

w(f ; t) = supu,v∈[0,1]|b−a|≤t

|f(b)− f(a)| fur t ≥ 0.

Satz 3.6 Ist f eine stetige Funktion auf [0, 1], dann gilt fur alle endlichen Folgen X =

(x1, . . . , xN), deren Punkte aus [0, 1] sind,∣∣∣∣∣ 1

N

N∑i=1

f(xi)−∫ 1

0

f(u) du

∣∣∣∣∣ ≤ w(f ;D∗N(X)).

Beweis: Sei x1 ≤ x2 ≤ · · · ≤ xN . Aus dem Mittelwertsatz fur Integrale folgt

∫ 1

0

f(u) du =N∑i=1

∫ iN

i−1N

f(u) du =1

N

N∑i=1

f(ti)

mit i−1N< ti <

iN

. Daraus folgt

1

N

N∑i=1

f(xi)−∫ 1

0

f(u) du =1

N

N∑i=1

(f(xi)− f(ti)).

Nun gilt |xi − ti| ≤ D∗N(X) fur alle i, durch Satz 3.4 folgt die Behauptung.

Die mehrdimensionale Fehlerschranken

Die Koksma-Ungleichung lasst sich im mehr-dimensionalen Fall anwenden.

Definition 3.11 Sei f eine Funktion aus [0, 1]d mit d ≥ 2. Unter einer Partition P

von [0, 1]d versteht man eine Menge von d endlichen Folgen (η(j)0 , . . . , η

(j)mj) mit 0 =

η(j)0 ≤ η

(j)1 ≤ · · · ≤ η

(j)mj = 1 fur j=1,...,d. Der ∆j-Operator ist definiert als

∆jf(x(1), . . . , x(j−1), η(j)i , x(j+1), . . . , x(d))

= f(x(1), . . . , x(j−1), η(j)i+1, x

(j+1), . . . , x(d))− f(x(1), . . . , x(j−1), η(j)i , xj+1, . . . , xd).

Definition 3.12 Sei eine Funktion f aus [0, 1]d. Setze

V d(f) := supP

m1−1∑i1=0

· · ·md−1∑id=0

∣∣∣∆1,...,df(η(1)i1, . . . , η

(d)id

)∣∣∣ ,


wobei P alle Partitionen aus [0, 1]d durchlauft. Wenn V d(f) endlich ist, dann sagt man,

f hat eine beschrankte Variation im Sinne von Vitali.

Weitere Formulierungen

V d(f) = supP

∑J∈P

|∆j(f)|,

wobei P alle Partitionen aus [0, 1]d durchlauft. Die Gleichheit

V d(f) =

∫ 1

0

. . .

∫ 1

0

∣∣∣∣ ∂df

∂u1 . . . ∂ud

∣∣∣∣ du1 . . . dud

gilt wenn die partielle Ableitung stetig auf [0, 1]d ist.

Es folgt direkt aus der Definition des ∆j-Operators, dass V d(f) = 0 wenn die Funktion

f aus [0, 1]d von weniger als d Variablen abhangt.

Definition 3.13 Hat f eine Einschrankung auf jeden Bereich F von [0, 1]d mit der

Dimension 1, 2, . . . , d − 1 eine beschrankte Variation im Sinne von Vitali, dann hat f

eine beschrankte Variation auf [0, 1]d im Sinne von Hardy und Krause.

V (f) =d∑

k=1

∑1≤i1≤i2<···<ik≤d

V k(f ; i1, . . . , ik).

Definition 3.14 Sei f eine Funktion aus [0, 1]d mit d ≥ 2. Der ∆∗j -Operator fur j=1,...,d

ist definiert als

∆∗jf(x(1), . . . , x(d))

= f(x(1), . . . , x(j−1), 1, x(j+1), . . . , x(d))− f(x(1), . . . , x(j−1), 0, x(j+1), . . . , x(d))

Sei nun ein Ausdruck F (r, . . . , r+ p− 1; r+ p, . . . , l) gegeben, der von der Partition der

Variation ir, . . . , il in die Teile ir, . . . , ir+p−1 und ir+p, . . . , il abhangt, dann steht∑r,...,s;p

F (r, . . . , r + p− 1; r + p, . . . , l)

fur die Summe aller Ausdrucke, die von F (r, . . . , r + p − 1; r + p, . . . , l) durch das

sukzessive Ersetzen der gegebenen Partitionen ir, . . . , il durch alle anderen Partition

von dieser Menge in eine Menge aus p und eine Menge aus l − r − p + 1 Variablen

abgeleitet wurden. Jede Partition wird genau einmal verwendet. Wenn entweder p = 0


oder p = l− r+ 1 ist, dann interpretiert man die Summe als zu einem Term reduzierte

Summe.

Lemma 3.2 Seien P = (η(j)0 , . . . , η

(j)imj

) und Q = (ξ(j)0 , . . . , ξ

(j)imj

) zwei Partitionen von

[0, 1]d, bestehend jeweils aus d Folgen fur j=1,...,d. Außerdem seien f und g zwei gegebene

Funktionen aus [0, 1]d. Dann gilt

m1−1∑i1=0

· · ·mk−1∑ik=0

f(ξ(1)i1+1, . . . , ξ

(d)id+1)∆1,...,dg(η

(1)i1, . . . , η

(d)id

)

=d∑p=0

(−1)p∑

1,...,d;p

∆p+1,...,d

m1∑i1=0

. . .

mp∑ip=0

g(η1i1, . . . , ηpip , x

p+1, . . . , xd)

∆1,...,pf(ξ1i1, . . . , ξpip , x

p+1, . . . , xd).

Beweis: Siehe L. Kuipers und H. Niederreiter Uniform Distribution of Sequences.

Satz 3.7 (Koksma-Hlawka-Ungleichung) Es besitze die Funktion f die beschrank-

te Variation im Sinne von Hardy und Krause auf [0, 1]d. Außerdem sei Xd = (x1, . . . , xN)

eine endliche Folge von Punkten, sodass deren Punkte aus [0, 1]d sind. Weiter sei

Xp = (xj1 , . . . , xjp) die Projektion der Folge Xd auf den (k − p)-dimensionalen Be-

reich von [0, 1]d. Dann gilt∣∣∣∣∣ 1

N

N∑i=1

f(xi)−∫Idf(u) du

∣∣∣∣∣ ≤d∑p=1

∑1,...,d;p

D∗N(Xdp+1,...,d)V

p(f(. . . , 1, . . . , 1)),

wobei V p(f(. . . , 1, . . . , 1)) die p-dimensionale Variation von f(x(1), . . . , x(p), 1, . . . , 1) auf

[0, 1]p im Sinne von Vitali ist. Wo der Term der Summe p = 1 ist, versteht man

D∗N(Xd)V d(f). Die Diskrepanz D∗N(Xdp+1,...,d) ist berechnet auf dem Bereich von [0, 1]d

wo sie enthalten ist.

Beweis: Fur eine Teilmenge M aus [0, 1]d definiere die Zahlfunktion A(M;Xd) durch

A(M). Definiere eine Funktion g auf [0, 1]d durch

g(x) = g(x(1), . . . , x(d)) =1

NA([0, x(1))× · · · × [0, x(d)))− x(1) . . . x(d).

Setze

D∗N(Xd) = supx∈[0,1]d

|g(x)|


und

D∗N(Xdp+1,...,d) = sup

(x(1),...,x(p))∈[0,1]p

∣∣g(x(1), . . . , x(p), 1, . . . , 1)∣∣ .

Setze xi = (x(1)i , . . . , x

(d)i ) fur i=1,...,N.

Eine zulassige Partitionierung von [0, 1]d durch P und Q wird als Paar (P,Q) aufgefasst.

P besteht aus d Folgen (η(j)0 , . . . , η

(j)imj

) fur j=1,...,d und Q besteht ebenfalls aus d Folgen

(ξ(j)0 , . . . , ξ

(j)imj+1

) fur j=1,...,d, die in folgender Relation zueinander stehen

0 = ξ(j)0 = η

(j)0 ≤ ξ

(j)1 < η

(j)1 ≤ ξ

(j)2 < η

(j)2 ≤ · · · ≤ ξ

(j)imj

< η(j)imj

= ξ(j)imj+1

= 1 fur j=1,...,d.

Außerdem soll fur jedes j fur j=1,...,d die Folge (ξ(j)0 , . . . , ξ

(j)imj+1

) mindestens die (x(j)1 , . . . , x

(j)N )

enthalten. Mit dieser zulassigen 2er-Partition lasst sich das vorherige Lemma mit den

Funktion f und g anwenden und man erhalt

m1−1∑i1=0

· · ·mk−1∑ik=0

f(ξ(1)i1+1, . . . , ξ

(d)id+1)∆1,...,dg(η

(1)i1, . . . , η

(d)id

)

=1

N

m1−1∑i1=0

· · ·mk−1∑ik=0

f(ξ(1)i1+1, . . . , ξ

(d)id+1)∆1,...,dA([0, η

(1)i1

)× · · · × [0, η(d)id

))

−m1−1∑i1=0

· · ·mk−1∑ik=0

f(ξ(1)i1+1, . . . , ξ

(d)id+1)∆1,...,dη

(1)i1. . . η

(d)id.

Nun gilt

∆1,...,dA([0, η(1)i1

)× · · · × [0, η(d)id

))

= ∆1,...,d−1A([0, η(1)i1

)× · · · × [0, η(d)id+1

))− A([0, η(1)i1

)× · · · × [0, η(d)id

))

= ∆1,...,dA([0, η(1)i1

)× · · · × [0, η(d−1)id−1

)× [η(d)id, η

(d)id+1

))

= · · · = A([η(1)i1, η

(1)i1+1)× · · · × [η

(d)id, η

(d)id+1

)).

Somit reduziert der erste Term zu

1

N

m1−1∑i1=0

· · ·mk−1∑ik=0

f(ξ(1)i1+1, . . . , ξ

(d)id+1)A([η

(1)i1, η

(1)i1+1)× · · · × [η

(d)id, η

(d)id+1

)).

Es sind nur diejenigen d-Tupel (i1, . . . , id) von Bedeutung, die ein xi fur 1,...,N im

Quader [η(1)i1, η

(1)i1+1) × · · · × [η

(d)id, η

(d)id+1

) haben. Sollte dieser Fall eintretten und die


Zusatzbedingung an Q, dass xi = ξ(1)i1+1, . . . , ξ

(d)id+1 ist, dann ist die Gleichung gerade

1N

∑Ni=1 f(xi) − f(ξ

(1)i1+1, . . . , ξ

(d)id+1)∆1,...,dg(η

(1)i1, . . . , η

(d)id

). Nun ist g(x) = 0, wenn aber

mindestens eine Koordinate verschwindet und ausserdem ist g(1, . . . , 1) = 0. Wenn

p = 0 ist in diesem Term, dann stimmt er mit dem aus der Gleichung aus Lemma

3.2 uberein. Es ist ∆1,...,dg(x(1), . . . , x(d))f(x(1), . . . , x(d)) daher 0. Außerdem bleiben fur

1 ≤ p ≤ d nur diese Terme uberig, bei denen die Variablen x(p+1), . . . , x(d) durch 1

ersetzt werden. Es folgt aus Lemma 3.2

d∑p=1

(−1)p∑

1,...,d;p

m1∑i1=0

· · ·mp∑ip=0

g(η(1)i1, . . . , η

(p)ip, 1, . . . , 1)∆1,...,pf(ξ

(1)i1+1, . . . , ξ

(p)ip+1, 1, . . . , 1).

Der Term lasst sich nach oben abschatzen mit

p∑p=1

∑1,...,d

m1∑i1=0

· · ·mp∑ip=0

∣∣∣g(η(1)i1, . . . , η

(p)ip, 1, . . . , 1)

∣∣∣ · ∣∣∣∆1,...,pf(ξ(1)i1+1, . . . , ξ

(p)ip+1, 1, . . . , 1)

∣∣∣ .Der Absolutwert bringt mit sich, dass g gleichmaßig beschrankt ist durch D∗N(Xd

p+1,...,d).

Die restliche Summe uber i1, . . . , ip ist dominiert durch V (p)(f(. . . , 1, . . . , 1)).

Daraus folgt∣∣∣∣∣ 1

N

N∑i=1

f(xi)−m1−1∑i1=0

· · ·md−1∑id=0

f(ξ(1)i1+1, . . . , ξ

(d)id+1)∆1,...,dη

(1)i1. . . η

(d)id

∣∣∣∣∣≤

d∑p=1

∑1,...,d;p

D∗N(Xdp+1,...,d)V

(p)(f(. . . , 1, . . . , 1)).

Die ∆1,...,dη(1)i1. . . η

(d)id

= (η(1)i1+1 − η

(1)i1

) . . . (η(1)id+1 − η

(1)id

) und die Summe uber i1, . . . , id

sind eine Riemann-Summe fur∫

[0,1]df(u) du. Die anderen Terme sind unabhangig von

der gewahlten 2er-Partition (P,Q). Der Beweis wird dadurch komplettiert, dass man

(P,Q) durch eine Folge von zulassigen 2er-Partitionen mit

max1≤j≤d

max0≤i<mj

(η(1)i1+1 − η

(1)i1

)→ 0

durchlaufen lasst.

Satz 3.8 Fur alle endlichen Folgen Xd = (x1, . . . , xN) von Punkten, deren Punkte


aus [0, 1]d sind, und jedes beliebige ε > 0 existiert eine Funktion f ∈ C∞([0, 1]d) mit

V (f) = 1, sodas gilt ∣∣∣∣∣ 1

N

N∑i=1

f(xi)−∫

[0,1]df(u) du

∣∣∣∣∣ < D∗N(Xd)− ε.

Beweis: Nach Definition von D∗N(Xd) existiert eine Quader Q∗ =∏d

i=1[0, bi) mit 0 <

bi ≤ 1 fur alle i und es gilt∣∣∣∣A(Q∗;Xd)

N− λd(Q∗)

∣∣∣∣ > D∗N(Xd)− ε

2.

Außerdem existiert eine zweite Quader Q∗2 =∏d

i=1[0, ti] mit 0 ≤ ti ≤ bi und bi− ti < ε2d

fur alle i so, dass Q∗ \ Q∗2 kein xi enthalt. Fur gegebene t und b mit 0 ≤ t < b ≤ 1

sei ft,b ∈ C∞([0, 1]) eine nicht steigende Funktion mit ft,b(a) = 1 fur 0 ≤ a ≤ t und

ft,b(a) = 0 fur t ≤ a ≤ b. Dann ist

f(a) = f(a1, . . . , ad) =d∏i=1

fti,bi(ai)

eine Funktion in C∞([0, 1]d) mit 0 ≤ f(a) ≤ 1 fur alle a ∈ [0, 1]d. Sei f(a) = 1 fur

ein a ∈ Q∗2 und f(a) = 0 fur a /∈ Q∗. Daraus folgt V d(f) = 1, da f(a) = 0 und wenn

eine Koordinate von a den Wert 1 annimmt, folgt V k(f ; i1, . . . , ik) = 0 fur 1 ≤ k < d.

Daraus folgt V (f) = 1. Es ist ersichtlich, dass

1

N

N∑i=1

f(xi) =A(Q∗;Xd)

N

und dass

λd(Q∗2) ≤∫

[0,1]df(u) du ≤ λd(Q∗).

Da λd(Q∗)− λd(Q∗2) ≤∑d

i=1(vi − ti) < ε2

gilt folgt, dass∣∣∣∣∣ 1

N

N∑i=1

f(xi)−∫

[0,1]df(u) du

∣∣∣∣∣ >∣∣∣∣A(Q∗;Xd)

N− λd(Q∗)

∣∣∣∣− ε

2> D∗N(Xd)− ε.


Definition 3.15 Das Stetigkeitsmaß fur eine stetige Funktion f auf [0, 1]d ist definiert

als

w(f ; t) = supa,b∈Id‖b−a‖≤t

|f(b)− f(a)| fur t ≥ 0,

wobei ‖a‖ = max1≤i≤d

|ai| fur a = (a1, . . . , ad) ∈ Rd.

Satz 3.9 Ist f stetig auf Id, dann gilt fur alle Folgen Xd = (x1, . . . , xN) von Punkten

deren Punkte aus [0, 1]d sind∣∣∣∣∣ 1

N

N∑i=1

f(xi)−∫

[0,1]df(u) du

∣∣∣∣∣ ≤ 4 · w(f ;D∗N(Xd)1d ).


Alle diese Fehlerschranken laufen darauf hinaus, dass eine kleine Sterndiskrepanz einen

kleinen Integrationsfehler der Quasi-Monte-Carlo-Verfahren garantiert.

Betrachte nun noch etwas allgemeinere Integrationsbereiche.

Satz 3.10 Es besitze f eine beschrankte Variation V (f) auf [0, 1]d im Sinne von Hardy

und Krause, C sei eine konvexe Menge aus [0, 1]d, dann gilt es fur alle endlichen Folgen

Xd = (x1, . . . , xN) von Punkten, deren Punkte aus [0, 1]d sind∣∣∣∣∣∣ 1

N

∑xi∈Ci

f(xi)−∫Cf(u) du

∣∣∣∣∣∣ ≤ (V (f) + |f(1, . . . , 1)|)TN(Xd).


Satz 3.11 Es besitze f eine beschrankte Variation V (f) auf [0, 1]d im Sinne von Hardy

und Krause, J sei aus einer Familie Jordanm-messbarer Funktionen Mj, dann gilt es

fur alle endlichen Folgen Xd = (x1, . . . , xN) von Punkten, deren Punkte aus [0, 1]d sind,∣∣∣∣∣∣ 1

N

∑xi∈Ji

f(xi)−∫Jf(u) du

∣∣∣∣∣∣ ≤ (V (f) + |f(1, . . . , 1)|)DN(Mj;Xd).

Beweis: Siehe L. Kuipers und H. Niederreiter Uniform Distribution of Sequences. Der

Ausdruck |f(1, . . . , 1)| in den letzten beiden Satzen ist notwendig auch fur konstante

Funktionen f .

3.3. PRAKTISCHE EIGENSCHAFT 37

3.3 Praktische Eigenschaft

Bei experimentellen Untersuchungen weisen dem Quasi-Monte-Carlo-Verfahren zur Be-

wertung pfadunabhangiger Optionen ein deutlich besseres Konvergenzverhalten auf als

Monte-Carlo-Verfahren.

Problematisch ist der Einsatz der Quasi-Monte-Carlo-Verfahren bei der Integration von

hochdimensionalen f , die bei pfadabhangigen Optionen und Optionen auf viele Aktien

auftreten.

In der asymptotischen Betrachtungsweise wachst (logN)d ·N−1 nicht so stark wie N−12 ,

wobei (logN)d · N−1 die Konvergenzordnung von Quasi-Monte-Carlo-Verfahren und

N−12 die Konvergenzordnung von Monte-Carlo-Vefahren darstellt.

3.3.1 Integral-Approximation

Folgendes Integral ist zu berechnen

I(f) =

∫[0,1]10

exp(u1 + . . .+ u10) du1 . . . du10

=

∫ 1

0

exp(u1) du1 · . . . ·∫ 1

0

exp(u10) du10

= (exp(1)− 1)10

= 224.3592464857479

Die Monte-Carlo bzw Quasi-Monte-Carlo-Approximation ist

I(f) =

∫[0,1]10

exp(u1 + . . .+ u10) du1 . . . dud ≈1

N

N∑i=1

exp(x(1)i + . . .+ x

(10)i ).

Das Integral wird mit einer Pseudo-Zufallsfolge und verschiedenen Niederdiskrepanz-

Folgen naherungsweise berechnet

• Pseudo-Zufallsfolge fur Aufwand N · d = 105 · 10

⇒ I(f) ≈ 223.9982885324653

• Halton-Folge fur Aufwand N · d = 105 · 10


⇒ I(f) ≈ 224.0982144317739

• Faure-Folge fur Aufwand N · d = 105 · 10

⇒ I(f) ≈ 224.2106475737543

• Sobol-Folge fur Aufwand N · d = 105 · 10

⇒ I(f) ≈ 224.3359566914221

Die Ergebnisse zeigen, dass bei gleichem Aufwand von 106 die Niederdiskrepanz-Folgen

bessere Approximationen fur das Integral liefern als die Pseudo-Zufallsfolge.

3.3.2 Vergleich Konvergenzverhalten

Um das Konvergenzverhalten des Monte-Carlo-Verfahrens und Quasi-Monte-Carlo-Verfahrens

zu untersuchen, wurden fur N = [103, 104, 105, 106] Simulationen mit den Folgen zu Ap-

proximation Wertes des Integrals durchgefuhrt, dabei ist die Dimension d = 10.

Abbildung 3.7 Integral Konvergenzverhalten

104

105

106

107

10−6

10−5

10−4

10−3

10−2

10−1

Integral Konvergenzverhalten

Aufwand : N x d

Rel

ativ

er F

ehle

r

Monte−CarloHalton−FolgeFaure−FolgeSobol−Folge

In Abbildung 3.7 wurden bei den Berechnung Wertes des Integrals auftretenden relati-

ven Fehler gegen den Aufwand N · d auf einer logarithmischen Skala aufgetragen.

3.3. PRAKTISCHE EIGENSCHAFT 39

Abbildung 3.8 Integral Lineare Regressionsanalyse

104

105

106

107

10−6

10−5

10−4

10−3

10−2

10−1

Integral Lineare Regressionsanalyse

Aufwand : N x d

Rel

ativ

er F

ehle

r

Monte−Carlo : − 0.4988Halton−Folge : − 0.8751Faure−Folge : − 1.0199Sobol−Folge : − 1.0236

In Abbildung 3.8 wurde eine lineare Regressionsanalyse der berechneten Punkte durch-

gefuhrt und die jeweilige Geraden-Steigung (Konvergenz-Geschwindigkeit) abgelesen.

Die Tabelle 3.1 dokumentiert die relatven Fehler und die Konvergenzverhalten, bei der

Berechnung des Wertes des Integrals.

Aufwand Pseudo-Zufallsolge Halton-Folge Faure-Folge Sobol-Folge

103 · 10 0.027437329618 0.050570903864 0.007611077730 0.006184719899

104 · 10 0.008788743788 0.008569960038 0.000662255681 0.000389833861

105 · 10 0.002758624862 0.001163455743 0.000074817839 0.000048145171

106 · 10 0.000877448667 0.000119107581 0.000006271552 0.000004812508

Konvergenz 0.498859452388 0.875111285826 1.019923376267 1.023917233862

Tabelle 3.1: Integral Konvergenzverhalten Tabelle

Die Ergebnisse zeigen, dass bei allen betrachteten Folgen mit zunehmendem Aufwand

die relativen Fehler gegen 0 konvergieren. Dabei weisen die Niederdiskrepanz-Folgen

gegenuber der Pseudo-Zufallsfolge eine deutlich kleinere relative Fehler und schnellere

Konvergenzverhalten auf, eine auch in der Praxis festgestellte Uberlegenheit des Quasi-

Monte-Carlo-Verfahrens zur Bewertung fur niedrig-dimensionale Optionsprobleme.


3.4 Randomisierte Quasi-Monte-Carlo-Verfahren

Niederdiskrepanz-Folgen sind auf die gleichmaßige Ausfullung des Integrationsberei-

ches mit Folgengliedern ausgerichtet. Im Unterschied zur Verwendung von Pseudo-

Zufallsfolgen versucht man durch Niederdiskrepanz-Folgen nicht, Eigenschaften wie die

Unabhangigkeit der einzelnen Glieder zueinander nachzubilden. Daher kann auch kei-

ne Schatzung der Fehlergroße auf Basis einzelner unabhangiger Funktionswerte f(xi)

durchgefuhrt werden.

Randomisierte Quasi-Monte-Carlo-Verfahren kombinieren Niederdiskrepanz-Folgen mit

Pseudo-Zufallsfolgen, um unabhangige Schatzwerte fur den Wert des Integrals zu ge-

nerieren. Mittels dieser Schatzwerte kann dann eine empirische Standardabweichung

berechnet werden.

3.4.1 Randomisierungsansatz von Tuffin

Im Randomisierungsansatz von Tuffin werden die Glieder der gewahlten Niederdiskrepanz-

Folge Xd = (x1, ..., xN) mittels eines in [0, 1]d uniform verteilten Zufallsvektors ξk ver-

schoben. Anstelle der ursprunglichen Niederdiskrepanz-Folge betrachtet man die neue

Folge Xdξk

= (x1 + ξk, . . . , xN + ξk)

Cξk =1

N

N∑i=1

f(xi + ξk)

3.4.2 Erwartungswert und Varianz

Cξk besitzt fur jede in [0, 1]d Riemann-integrierbare Funktion f den Erwartungswert

E[Cξk ] = E

[1

N

N∑i=1

f(xi + ξk)

]= E[f(xi + ξk)] = E[f(ξk)].

Fur die Varianz der Cξk gilt nach dem Satz von Tuffin

V ar[Cξk ] ∈ O(

(logN)2d

N2

),

falls die Stern-Diskrepanz D∗N der verwendeten Folge Xd in O(

(logN)d

N

)liegt.

3.4. RANDOMISIERTE QUASI-MONTE-CARLO-VERFAHREN 41

3.4.3 Fehler des Schatzers

Berechnet man die Zufallsvariablen Cξk fur unabhangige uniform verteilte Vektoren ξk

fur k=1,...,K, kann der Fehler des Schatzers

C =1

K

K∑k=1

Cξk

auf Basis der empirischen Standardabweichung der Cξk

σξk =

√√√√ 1

K − 1

K∑k=1

(Cξk − C)2

und der Normalverteilungsannahme fur C abgeschatzt werden.

Kapitel 4

Numerische Ergebnisse

In diesem Kapitel soll dargestellt werden, auf welche Weise der Quasi-Monte-Carlo-

Verfahren Algorithmus implementiert wird. Dabei werden das Black-Scholes-Modell

betrachtet und fur verschiedene Optionstypen implementiert und zahlreiche Testlaufe

durchgefuhrt, um die Eignung der Verfahren in der Praxis zu verifizieren.

Das Quasi-Mote-Carlo-Verfahren wird in Hinblick auf

• Approximationsfehler des Verfahrens und

• Konvergenzverhalten des Verfahrens

untersucht.

4.1 Implementierung

Das Quasi-Monte-Carlo-Verfahren wird mit MATLAB realisiert, dabei werden die Quasi-

Zufallsfolge Generatoren, Auszahlungsfunktion fur die verschiedene Optionstypen und

Simulationen umfasst. (Quellcode in Anhang)

• Quasi-Zufallsfolgen Generatoren fur

– Halton-Folge

– Faure-Folge

– Sobol-Folge

• Generierung normalverteilter Zahlen

43

44 KAPITEL 4. NUMERISCHE ERGEBNISSE

– Transformation mittels der inversen Normalverteilungsfunktion, Inversions-

methode

• Simulationen von Aktienkurs-Pfadvektoren

– Black-Scholes-Modell mit Ansatz vom Euler-Maruyama-Verfahren

– Black-Scholes-Modell mit Ansatz von Ito-Formel

• Untersuchte pfadunabhangige Optionen

– Europaische Optionen

• Untersuchte pfadabhangige Optionen

– Asiatische Optionen mit diskretem geometrischem Mittel

– Asiatische Optionen mit kontinuierlichem geometrischem Mittel

• Berechnung des exakten Optionspreises

– Black-Scholes-Formel fur Europaische Optionen

– Black-Scholes-Formel fur Asiatische Optionen mit diskretem geometrischem

Mittel

– Black-Scholes-Formel fur Asiatische Optionen mit kontinuierlichem geome-

trischem Mittel

4.2 Bewertung von Optionen

Es wurden pfadunabhangige und pfadabhangige Optionen unter Verwendung verschie-

dener Niederdiskrepanz-Folgen mit verschiedenen Parametern des Black-Scholes-Modells

untersucht und mit den exakten Optionspreisen verglichen. Dabei wurde jede der 243

Kombinationen berucksichtigt, die durch Auswahl von Modellparametern auf den fol-

genden Wertebereichen gebildet werden konnen.

• Aktienkurs S0 ∈ 50, 100, 150

• Ausubungspreis K ∈ 52, 103, 155

• Zinssatz r ∈ 0.02, 0.03, 0.04

4.2. BEWERTUNG VON OPTIONEN 45

• Volatilitat σ ∈ 0.1, 0.25, 0.5

• Verfallszeit T ∈ 0.5, 1, 1.5 Jahre.

Durch die Wahl der Parameter sollen alle Auspragungen von Optionen (Option aus/im

Geld, Langlaufer/Kurzlaufer, niedriges/hohes Risiko), die in einem praktischen Einsatz

des Verfahrens auftreten konnen, berucksichtigt werden.

Jedes der einzelnen numerischen Ergebnisse weist eine Konvergenz des relativen Fehlers

gegen 0 auf, daher werde ich hier nur einige exemplarische Ergebnisse vorfuhren.

4.2.1 Europaische Optionen

Eine Europaische Call-Option mit Aktienkurs S0 = 100, Ausubungspreis K = 103,

Zinssatz r = 0.03, Volatilitat σ = 0.25, Verfallszeit T = 1 Jahr.

Der exakte Optionspreis der mit der Black-Scholes-Formel (Satz 2.3) berechnet wurde,

betragt 9.967521439954112.

Fur die Europaische Call-Option wurden N = [1 · 104, 2 · 104, 4 · 104, 8 · 104, 16 · 104] Si-

mulationen zur Berechnung des Optionspreises durchgefuhrt, dabei wird die Dimension

d = [1, 2, 4, 8, 16] im gleichen Verhaltnis wie N erhoht. Der Grund fur die gleichmaßige

Erhohung ist, den Diskretisierungs-Fehler des Verfahrens zu eliminieren.

Euro-Call Euler-Maruyama-Ansatz

Abbildung 4.1 Euro-Call Euler-Maruyama-Ansatz Konvergenzverhalten

104

105

106

107

10−4

10−3

10−2

10−1

Euro−Call Euler−Maruyama−Ansatz Konvergenzverhalten

Aufwand : N x d

Rel

ativ

er F

ehle

r

Halton−FolgeFaure−FolgeSobol−Folge


In Abbildung 4.1 wurden bei der Berechnung des Wertes des Europaischen Calls mit

Euler-Maruyama-Verfahren auftretenden relativen Fehler gegen den Aufwand N · d auf

einer logarithmischen Skala aufgetragen.

Abbildung 4.2 Euro-Call Euler-Maruyama-Ansatz Lineare Regressionsanalyse

104

105

106

107

10−4

10−3

10−2

10−1

Euro−Call Euler−Maruyama−Ansatz Lineare Regressionsanalyse

Aufwand : N x d

Rel

ativ

er F

ehle

r

Halton−Folge : − 0.5361Faure−Folge : − 0.8180Sobol−Folge : − 0.8017



Die Tabelle 4.1 dokumentiert die relativen Fehler und die Konvergenzverhalten, bei der

Berechnung des Wertes des Europaischen Calls mit dem Euler-Maruyama-Verfahren.

Aufwand Halton-Folge Faure-Folge Sobol-Folge

1 · 104 · 1 0.030424360989676 0.030424360989676 0.030424360989676

2 · 104 · 2 0.009114458292331 0.008550210303867 0.008550210303867

4 · 104 · 4 0.003118451087619 0.002623836864443 0.002005340801068

8 · 104 · 8 0.002472091651737 0.000530888077963 0.000933301778260

16 · 104 · 16 0.001420888689036 0.000420888689036 0.000355451259884

Konvergenz 0.536157385826403 0.818038454094376 0.801720358124229

Tabelle 4.1: Euro-Call Euler-Maruyama-Ansatz relativen Fehler und Konvergenzverhal-

ten Tabelle

Die Ergebnisse zeigen, dass bei allen betrachteten Niederdiskrepanz-Folgen mit zuneh-

mendem Aufwand die relativen Fehler gegen 0 konvergieren.


Euro-Call Ito-Ansatz

Abbildung 4.3 Euro-Call Ito-Ansatz Konvergenzverhalten

104

105

106

107

10−4

10−3

10−2

Euro−Call Ito−Ansatz Konvergenzverhalten

Aufwand : N x d

Rel

ativ

er F

ehle

r


In Abbildung 4.3 wurden bei der Berechnung des Wertes des Europaischen Calls mit Ito-

Formel auftretenden relativen Fehler gegen den Aufwand N ·d auf eine logarithmischen

Skala aufgetragen.

Abbildung 4.4 Euro-Call Ito-Ansatz Lineare Regressionsanalyse

104

105

106

107

10−4

10−3

10−2

Euro−Call Ito−Ansatz Lineare Regressionsanalyse

Aufwand : N x d

Rel

ativ

er F

ehle

r






Berechnung des Wertes des Europaischen Calls mit dem Ito-Formel.


1 · 104 · 1 0.002323057365162 0.002323057365162 0.002323057365162

2 · 104 · 2 0.001699783491951 0.000990012640727 0.000990012640727

4 · 104 · 4 0.001528714599288 0.001088626220606 0.000104452147728

8 · 104 · 8 0.002171297095342 0.000144493660733 0.000507971446410

16 · 104 · 16 0.001432355216490 0.005721186915096 0.000408546514401

Konvergenz 0.052103239010679 0.008793067436457 0.298880230273579

Tabelle 4.2: Euro-Call Ito-Ansatz relativen Fehler und Konvergenzverhalten Tabelle

Die Ergebnisse zeigen, dass das Ito-Verfahren schon bei geringen Aufwand kleine rela-

tiven Fehler aufweist. Bei pfadunabhangige Optionen sollte man daher direkt mit der

Ito-Formel die Option bewerten, dadurch spart man an Aufwand und erhalt auch eine

bessere Genauigkeit bei den Werten.

4.2.2 Asiatische Optionen diskret geometrisch

Eine Asiatische Call-Option mit Aktienkurs S0 = 100, Ausubungspreis K = 103, Zins-

satz r = 0.03, Volatilitat σ = 0.25, Verfallszeit T = 1 Jahr und Zeitpunkt M = 12

(M = d).

Der exakte Optionspreis mit der Black-Scholes-Formel (Satz 2.4) berechnet wurde, be-

tragt 5.181006824978432.

Die Asiatische Call Option mit diskretem geometrischem Mittel wurden N = [5·103, 10·103, 20·103, 40·103, 80·103] Simulationen zur Berechnung des Optionspreis durchgefuhrt,

dabei ist die Dimension d = 12. Hier wird die Dimension nicht erhohrt, da die Zeit-

punkte schon fest sind.


Asia-Call diskret Euler-Maruyama-Ansatz

Abbildung 4.5 Asia-Call diskret Euler-Maruyama-Ansatz Konvergenzverhalten

105

106

107

10−4

10−3

10−2

Asia−Call diskret Euler−Maruyama−Ansatz Konvergenzverhalten

Aufwand : N x d

Rel

ativ

er F

ehle

r


In Abbildung 4.5 wurden bei der Berechnung des Wertes des Asiatischen Calls mit

diskretem geometrischem Mittel mit Euler-Maruyama-Verfahren auftretenden relativen

Fehler gegen den Aufwand N · d auf eine logarithmischen Skala aufgetragen.

Abbildung 4.6 Asia-Call diskret Euler-Maruyama-Ansatz Lineare Regressionsanalyse

105

106

107

10−4

10−3

10−2

Asia−Call diskret Euler−Maruyama−Ansatz Lineare Regressionsanalyse

Aufwand : N x d

Rel

ativ

er F

ehle

r


In Abbildung 4.6 wurde eine Lineare Regressionsanalyse der berechneten Punkte durch-

gefuhrt, und die jeweilige Geraden-Steigung (Konvergenz-Geschwindigkeit) abgelesen.


Die Tabelle 4.3 dokumentiert die relatven Fehler und die Konvergenzverhalten, bei der

Berechnung des Wertes des Asiatischen Calls mit diskretem geometrischem Mittel mit

Euler-Maruyama-Verfahren.


5 · 103 · 12 0.008297610188939 0.004205293852263 0.005442643428294

10 · 103 · 22 0.004197508121534 0.003376623709701 0.005988231933268

20 · 103 · 12 0.004482304605234 0.000437465842291 0.002280001477266

40 · 103 · 12 0.003421711463257 0.000410976385458 0.001503310133530

80 · 103 · 12 0.002707470416369 0.000698139011228 0.001825235745356

Konvergenz 0.352631562479756 0.821969476958533 0.514643684955722

Tabelle 4.3: Asia-Call diskret Euler-Maruyama-Ansatz relativen Fehler und Konver-

genzverhalten Tabelle


mende Aufwand die relativen Fehlern gegen 0 konvergieren, und die Faure-Folge die

kleinsten relativen Fehler und die schnellste Konvergenzverhalten aufweist.

Asia-Call diskret Ito-Ansatz

Abbildung 4.7 Asia-Call diskret Ito-Ansatz Konvergenzverhalten

105

106

107

10−4

10−3

10−2

Asia−Call diskret Ito−Ansatz Konvergenzverhalten

Aufwand : N x d

Rel

ativ

er F

ehle

r


In Abbildung 4.7, wurden bei der Berechnung des Wertes des Asiatischen Calls mit

diskretem geometrischem Mittel mit Ito-Formel auftretenden relativen Fehler gegenuber

den Aufwand N · d auf eine logarithmischen Skala aufgetragen.


Abbildung 4.8 Asia-Call diskret Ito-Ansatz Lineare Regressionsanalyse

105

106

107

10−4

10−3

10−2

Asia−Call diskret Ito−Ansatz Lineare Regressionsanalyse

Aufwand : N x d

Rel

ativ

er F

ehle

r





Berechnung des Wertes des Asiatischen Calls mit diskretem geometrischem Mittel mit

Ito-Formel.


5 · 103 · 12 0.007682205025533 0.003104458614172 0.003824216979071

10 · 103 · 22 0.003283784415296 0.001914685846263 0.004549129018754

20 · 103 · 12 0.003517395409312 0.000924301594618 0.000605733594411

40 · 103 · 12 0.002252138758256 0.000962015619334 0.000099726456442

80 · 103 · 12 0.001406436098152 0.000672484567038 0.000295488317570

Konvergenz 0.544301662368218 0.540651274638860 1.289945301999592

Tabelle 4.4: Asia-Call diskret Ito-Ansatz relativen Fehler und Konvergenzverhalten Ta-

belle


mende Aufwand die relativen Fehler gegen 0 konvergieren, und die Sobol-Folge die

kleinsten relativen Fehler und die schnellste Konvergenzverhalten aufweist.


4.2.3 Asiatische Optionen kontinuierlich geometrisch

Eine Asiatische Call-Option mit Aktienkurs S0 = 100, Ausubungspreis K = 103, Zins-

satz r = 0.03, Volatilitat σ = 0.25, Verfallszeit T = 1 Jahr.

Der exakte Optionspreis mit der Black-Scholes-Formel (Satz 2.5) berechnet wurde, be-

tragt 4.769426376039782.

Fur die Asiatische call-Option mit kontinuierlichem geometrischem Mittel wurden N =

[5 ·103, 10 ·103, 20 ·103, 40 ·103, 80 ·103] Simulationen zur Berechnung des Optionspreises

durchgefuhrt, dabei wird die Dimension d = [5, 10, 20, 40, 80] im gleichen Verhaltnis wie

N erhoht. Der Grund fur die gleichmaßige Erhohung ist, den Diskretisierungs-Fehler

des Verfahrens zu eliminieren.

Asia-Call kontinuierlich Euler-Maruyama-Ansatz

Abbildung 4.9 Asia-Call kontinuierlich Euler-Maruyama-Ansatz Konvergenzverhalten

104

105

106

107

10−3

10−2

10−1

100

Asia−Call kontinuierlich Euler−Maruyama−Ansatz Konvergenzverhalten

Aufwand : N x d

Rel

ativ

er F

ehle

r


In Abbildung 4.9 wurden der bei Berechnung des Wertes des Asiatischen Calls mit

kontinuierlichem geometrischem Mittel mit Euler-Maruyama-Verfahren auftretenden

relativen Fehler gegen den Aufwand N · d auf eine logarithmischen Skala aufgetragen.


Abbildung 4.10 Asia-Call kontinuierlich Euler-Maruyama-Ansatz Lineare Regressi-

onsanalyse

104

105

106

107

10−3

10−2

10−1

100

Asia−Call kontinuierlich Euler−Maruyama−Ansatz Lineare Regressionsanalyse

Aufwand : N x d

Rel

ativ

er F

ehle

r




Die Tabelle 4.5 dokumentiert die relativen Fehler und die Konvergenzverhalten bei

der Berechnung des Wertes des Asiatischen Calls mit kontinuierlichem geometrischem

Mittel mit Euler-Maruyama-Verfahren.

Aufwand : N · d Halton-Folge Faure-Folge Sobol-Folge

5 · 103 · 5 0.192470297333721 0.198345750049843 0.197514959690855

10 · 103 · 10 0.094585807263369 0.097127773339536 0.095893407398357

20 · 103 · 20 0.041485194493920 0.040208601083472 0.046559267605944

40 · 103 · 40 0.007024399741505 0.006562359321326 0.026489625093971

80 · 103 · 80 0.020747447766315 0.005742242255734 0.011767399856740

Konvergenz 0.508921847548554 0.705405817442494 0.499709655042184

Tabelle 4.5: Asia-Call kontinuierlich Euler-Maruyama-Ansatz relativen Fehler und Kon-

vergenzverhalten Tabelle


mende Aufwand die relativen Fehler gegen 0 konvergieren.


Asia-Call kontinuierlich Ito-Ansatz

Abbildung 4.11 Asia-Call kontinuierlich Ito-Ansatz Konvergenzverhalten

104

105

106

107

10−3

10−2

10−1

100

Asia−Call kontinuierlich Ito−Ansatz Konvergenzverhalten

Aufwand : N x d

Rel

ativ

er F

ehle

r


In Abbildung 4.11 wurden der bei Berechnung des Wertes des Asiatischen Calls mit

kontinuierlichem geometrischem Mittel mit Ito-Formel auftretenden relativen Fehler

gegen den Aufwand N · d auf eine logarithmischen Skala aufgetragen.

Abbildung 4.12 Asia-Call kontinuierlich Ito-Ansatz Lineare Regressionsanalyse

104

105

106

107

10−3

10−2

10−1

100

Asia−Call kontinuierlich Ito−Ansatz Lineare Regressionsanalyse

Aufwand : N x d

Rel

ativ

er F

ehle

r





Die Tabelle 4.6 dokumentiert die relativen Fehlern und die Konvergenzverhalten, bei

der Berechnung des Wertes des Asiatischen Calls mit kontinuierlichem geometrischem

Mittel mit Ito-Formel.

Aufwand : N · d Halton-Folge Faure-Folge Sobol-Folge

5 · 103 · 5 0.195869972116154 0.202609956111959 0.201835217204533

10 · 103 · 10 0.095422304411796 0.098040640922101 0.098651870683991

20 · 103 · 20 0.041877486609307 0.041090883211523 0.046832835620440

40 · 103 · 40 0.006697344559224 0.005529652703121 0.026457239078278

80 · 103 · 80 0.021030439279334 0.001997893284779 0.011856437251069

Konvergenz 0.352631562479756 0.821969476958533 0.514643684955722

Tabelle 4.6: Asia-Call kontinuierlich Ito-Ansatz relativen Fehler und Konvergenzverhal-

ten Tabelle


mende Aufwand die relativen Fehler gegen 0 konvergieren.

Kapitel 5

Zusammenfassung und Ausblick

Das Gebiet der Optionsbewertung ist durch die Entwicklungen zu neuen und im-

mer komplexer werdenden Optionstypen und durch Verbesserungen im Bereich der

Aktienkurs-Modelle gepragt. Diese Entwicklung und die gestiegene Leistungsfahigkeit

der Parallelrechner haben das Interesse an den flexiblen Quasi-Monte-Carlo-Verfahren

neu geweckt.

Die experimentellen Untersuchungen bestatigen die Uberlegenheit des Quasi-Monte-

Carlo-Verfahren gegenuber den klassische Monte-Carlo-Verfahren in Bezug auf niedrig-

dimensionale Optionstypen. Dieser Uberlegenheit nimmt aber mit zunehmender Dimen-

sion ab, was eine Nachteil fur das Quasi-Monte-Carlo Verfahren darstellt. Zur Verbesse-

rung des Verfahrens gibt das Dimensions-Reduktions-Prinzip (effective dimension) und

weitere Niederdiskrepanz-Folgen, wie Niederreiter-Folgen, Lattice-Regeln, usw. Weitere

Verbesserungsmoglichkeiten konnten auch durch Wahl von anderen Diskretisierungsver-

fahren mit hoherer starker Ordnung, wie z.B dem Milstein-Verfahren, erreicht werden.

Mit dem Quasi-Monte-Carlo-Verfahren lassen sich auch komplizierte Optionen bewer-

ten, wie z.B.

• Bermuda-Optionen,

• Barrier-Optionen,

• Cap-Optionen,

• Shout-Optionen,

• Lokkback-Optionen,

• Multi-Asset-Optionen,

57

58 KAPITEL 5. ZUSAMMENFASSUNG UND AUSBLICK

• Outperformance-Optionen,

und auch mit weiteren Bewertungs-Modellen kombinieren, wie z.B. dem

• Black-Scholes-Modell mit variabler Verzinsung,

• Black-Scholes-Modell mit zeitabhangiger Volatilitat,

• Heston-Modell fur stochastische Volatilitat,

• Merton-Sprung-Diffusion-Modell und dem

• Libor-Markt Modell fur Zinsderivate,

auf die ich in dieser Bachelorarbeit nicht mehr eingehen werde, mit denen ich mich

jedoch in der Masterarbeit genauer beschafigen werde.

Anhang A

Anhang: Programmcodes

A.1 Niederdiskrepanz-Folgen-Generatoren

Programmcode A.1.1 Halton-Folge

function [qh,Z] = qmch(N,d)

format long;

hf = qrandstream(’halton’,d);

qh = rand(hf,N,d);

Z = icdf(’normal’,qh,0,1);

Programmcode A.1.2 Faure-Folge

function [Z] = qmcf(N,d)

format long;

Z = zeros(N:d);

qs = faure(N,d,13);

Z1 = icdf(’normal’,qs,0,1);

for i = 1 : d

Z(:,i) = Z1(i,1:N)’;

end

59

60 ANHANG A. ANHANG: PROGRAMMCODES

Programmcode A.1.3 Sobol-Folge

function [qs,Z] = qmcs(N,d)

format long;

sf = qrandstream(’sobol’,d);

qs = rand(sf,N,d);

Z = icdf(’normal’,qs,0,1);


function[a] = basexpflip(k,b)

if k

j = fix(log(k)/log(b)) + 1;

a = zeros(1,j);

q = b^(j-1);

for i = 1:j

a(i) = floor(k/q);

k = k - q*a(i);

q = q/b;

end

a = fliplr(a);

else

a = 0;

end

function[c] = comb(n,k)

if n < k

c = 0;

else

c = nchoosek(n,k);

end

function[i] = isint(x)

i = (x == floor(x));

A.1. NIEDERDISKREPANZ-FOLGEN-GENERATOREN 61


function[s] = faure(k,d,b)

if ~(isint(k) && k >= 0)

error(’Input argument "k" must be a non-negative integer’)

end

if ~(isint(d) && d > 0)

error(’Input argument "d" must be a positive integer’)

end

if ~(isint(b) && b > 1)

error(’Input argument "b" must be an integer greater than 1’)

end

s = zeros(d,k+1);

K = k;

D = d;

for k = 0:K

a = basexpflip(k,b);

J = length(a);

L = J - 1;

y = zeros(J,1);

g = b.^(1:J)’;

for d = 1:D

for j = 1:J

S = 0;

for l = 0:L

c = comb(l,j-1);

h = (d-1)^(l-j+1);

if isinf(h)

h = 0;

end

S = S + c*h*a(l+1);

end

y(j) = mod(S,b);

end

s(d,k+1) = sum(y./g);

end

end


A.2 Black-Scholes-Formel

Programmcode A.2.1 Black-Scholes-Formel Euro-Call

function [C0] = bsfec(S0,K,r,sigma,T)

format long;

d1 = (log(S0 / K) + (r + (sigma^2) / 2) * T) / (sigma * T);

d2 = d1 - sigma * sqrt(T);

C0 = S0 * normcdf(d1) - K * exp(-r * T) * normcdf(d2);

Programmcode A.2.2 Black-Scholes-Formel Asia-Call diskret

function [V] = bsacdgm(S0,K,r,sigma,T,M)

format long;

dt = T / M;

T1 = T - (M * (M - 1) * (4 * M + 1) * dt) / (6 * M^2);

T2 = T - ((M - 1) * dt) / 2;

A = exp(-r * (T - T2) - (sigma^2) * (T2 - T1) / 2);

d = (log(S0 / K) + (r - (sigma^2) / 2) * T2) / (sigma * sqrt(T1));

V = S0 * A * normcdf(d + sigma * sqrt(T1)) - K * exp(-r * T) * normcdf(d);

Programmcode A.2.3 Black-Scholes-Formel Asia-Call kontinuierlich

function [V] = bsackgm(S0,K,r,sigma,T)

format long;

d = (log(S0 / K) + 1 / 2 * (r - sigma^2 / 2) * T) / (sigma * sqrt(T / 3));

V = S0 * exp(- 1 / 2 * (r + sigma^2 / 6) * T) * normcdf(d + sigma * sqrt(T / 3))

- K * exp(- r * T) * normcdf(d);

A.3. INTEGRAL APPROXIMATION 63

A.3 Integral Approximation

Programmcode A.3.1 Integral Approximation Pseudo-Zufallsfolgen

function [msm] = kbm

format long;

w = (exp(1) - 1)^10;

N = [1000,10000,100000,1000000];

sm = zeros(10000,length(N));

msm = zeros(1,length(N));

for i = 1 : length(N)

for s = 1 : 1000

mc = rand(N(i),10);

z = zeros(1,N(i));

for j = 1 : N(i)

z(j) = exp(sum(mc(j,:)));

end

ns = mean(z);

sm(s,i)=abs(w-ns)/w;

end

msm(i) = mean(sm(:,i));

end


Programmcode A.3.2 Integral Approximation Halton-Folgen

function [ff] = kbf

format long;

w = (exp(1) - 1)^10;

N = [1000,10000,100000,1000000];

n = length(N);

ff = zeros(1,n);

for i = 1 : n

qf = faure(N(i),10,11);

z = zeros(1,N(i));

for j = 1 : N(i)

z(j) = exp(sum(qf(:,j)));

end

nf = mean(z);

ff(i)=abs(w-nf)/w;

end

Programmcode A.3.3 Integral Approximation Faure-Folgen

function [ff] = kbf

format long;

w = (exp(1) - 1)^10;

N = [1000,10000,100000,1000000];

n = length(N);

ff = zeros(1,n);

for i = 1 : n

qf = faure(N(i),10,11);

z = zeros(1,N(i));

for j = 1 : N(i)

z(j) = exp(sum(qf(:,j)));

end

nf = mean(z);

ff(i)=abs(w-nf)/w;

end

A.3. INTEGRAL APPROXIMATION 65

Programmcode A.3.4 Integral Approximation Sobol-Folgen

function [sf] = kbs

format long;

w = (exp(1) - 1)^10;

N = [1000,10000,100000,1000000];

sf = zeros(1,length(N));


s = qrandstream(’sobol’,10);

qs = rand(s,N(i),10);

z = zeros(1,N(i));

for j = 1 : N(i)

z(j) = exp(sum(qs(j,:)));

end

ns = mean(z);

sf(i) = abs(w - ns) / w;

end

Programmcode A.3.5 Integral Approximation Konvergenzverhalten

N = [1000,10000,100000,1000000];

M = 10;

N2 = N * M;

loglog(N2,kbm,’--ko’,’LineWidth’,2,’MarkerFaceColor’,’k’,’MarkerSize’,8);

hold on

loglog(N2,kbh,’--ro’,’LineWidth’,2,’MarkerEdgeColor’,’r’,’MarkerSize’,8);

hold on

loglog(N2,kbf,’--go’,’LineWidth’,2,’MarkerEdgeColor’,’g’,’MarkerSize’,8);

hold on

loglog(N2,kbs,’--bo’,’LineWidth’,2,’MarkerEdgeColor’,’b’,’MarkerSize’,8)

title(’Konvergenz’)

xlabel(’Aufwand : N x d’)

ylabel(’Relative Fehler’)

legend(’Monte-Carlo’,’Halton-Folge’,’Faure-Folge’,’Sobol-Folge’)

grid


Programmcode A.3.6 Integral Approximation Lineare Regressionsanalyseformat long;

N = [1000,10000,100000,1000000];

N2 = N * 10;

n = log10(N2);

rfmc = log10(kbmc);

pmc = polyfit(n,rfmc,1);

rfh = log10(kbh);

ph = polyfit(n,rfh,1);

rff = log10(kbf);

pf = polyfit(n,rff,1);

rfs = log10(kbs);

ps = polyfit(n,rfs,1);

N3 = [N2(1),N2(length(N))];

regmc = [10^(pmc(1) * log10(N2(1)) + pmc(2)),10^(pmc(1) * log10(N2(length(N))) + pmc(2))];

regh = [10^(ph(1) * log10(N2(1)) + ph(2)),10^(ph(1) * log10(N2(length(N)))

+ ph(2))];

regf = [10^(pf(1) * log10(N2(1)) + pf(2)),10^(pf(1) * log10(N2(length(N)))

+ pf(2))];

regs = [10^(ps(1) * log10(N2(1)) + ps(2)),10^(ps(1) * log10(N2(length(N)))

+ ps(2))];

loglog(N3,regmc,’-k’,’LineWidth’,2);

hold on

loglog(N3,regh,’-r’,’LineWidth’,2);

hold on

loglog(N3,regf,’-g’,’LineWidth’,2);

hold on

loglog(N3,regs,’-b’,’LineWidth’,2)

title(’Lineare Regressionsanalyse Integral’)


ylabel(’Relativer Fehler’)

legend(’Monte-Carlo’,’Halton-Folge’,’Faure-Folge’,’Sobol-Folge’)

grid on

A.4. EUROPAISCHE OPTIONEN 67

A.4 Europaische Optionen

Programmcode A.4.1 Euro-Call Euler-Maruyama-Ansatz Halton-Folgen

function [V] = qmcechalton(S0,K,T,r,Sigma,N,M)

format long;

h = T / M;

DW = qmch(N + 1,M);

S = zeros(N,M + 1);

S(:,1) = S0;

Z = zeros(1,N);

for i = 1 : N

for j = 1 : M

S(i,j + 1) = S(i,j) + r * S(i,j) * h + Sigma * S(i,j) * sqrt(h)

* DW(i + 1,j);

end

end

for i = 1 : N

Z(i) = S(i,M + 1) - K;

if Z(i) < 0

Z(i) = 0;

end

end

V = exp(-r) * mean(Z);


Programmcode A.4.2 Euro-Call Ito-Ansatz Halton-Folgen

function [V] = qmcechaltonito(S0,K,T,r,Sigma,N,M)

format long;

h = T / M;

DW = qmch(N + 1,M);

S = zeros(N,M + 1);

S(:,1) = S0;

Z = zeros(1,N);

for i = 1 : N

for j = 1 : M

S(i,j + 1) = S(i,j) * exp((r - Sigma^2 / 2) * h + Sigma * sqrt(h)

* DW(i + 1,j));

end

end

for i = 1 : N

Z(i) = S(i,M + 1) - K;

if Z(i) < 0

Z(i) = 0;

end

end

Fur Faure-Folge und Sobol-Folge sind die Programm-Codes Analog.

A.4. EUROPAISCHE OPTIONEN 69

Programmcode A.4.3 Euro-Call Euler-Maruyama-Ansatz Konvergenzverhaltenformat long;

N = [10,20,40,80,160];

M = [1,2,4,8,16];

N2 = N.*M;

n = log10(N2);

rfh = zeros(1,length(N));

rff = zeros(1,length(N));

rfs = zeros(1,length(N));


rfh(i) = abs(qmcechalton(100,103,1,0.03,0.25,N(i),M(i))

- 9.967521439954112) / 9.967521439954112;

rff(i) = abs(qmcecfaure(100,103,1,0.03,0.25,N(i),M(i))

- 9.967521439954112) / 9.967521439954112;

rfs(i) = abs(qmcecsobol(100,103,1,0.03,0.25,N(i),M(i))

- 9.967521439954112) / 9.967521439954112;

end

loglog(N2,rfh,’--ro’,’LineWidth’,2,’MarkerFaceColor’,’r’,’MarkerSize’,8);

hold on

loglog(N2,rff,’--go’,’LineWidth’,2,’MarkerFaceColor’,’g’,’MarkerSize’,8);

hold on

loglog(N2,rfs,’--bo’,’LineWidth’,2,’MarkerFaceColor’,’b’,’MarkerSize’,8)


xlabel(’Aufwand : N * d’)

ylabel(’Relative Fehler’)

legend(’Halton-Folge’,’Faure-Folge’,’Sobol-Folge’)

grid on

Fur Ito-Formel ist das Programm-Code Analog.


Programmcode A.4.4 Euro-Call Euler-Maruyama-Ansatz Lineare Regressionsanaly-

seN = [10000,20000,40000,80000,160000];

M = [1,2,4,8,16];

N2 = N.* M;

n = log10(N2);


rfh(i) = abs(qmcechalton(100,103,1,0.03,0.25,N(i),M(i))

- 9.967521439954112) / 9.967521439954112;

rff(i) = abs(qmcecfaure(100,103,1,0.03,0.25,N(i),M(i))

- 9.967521439954112) / 9.967521439954112;

rfs(i) = abs(qmcecsobol(100,103,1,0.03,0.25,N(i),M(i))

- 9.967521439954112) / 9.967521439954112;

end

rrfh = log10(rfh);

ph = polyfit(n,rrfh,1);

rrff = log10(rff);

pf = polyfit(n,rrff,1);

rrfs = log10(rfs);

ps = polyfit(n,rrfs,1);

N3 = [N2(1),N2(length(N))];

regh = [10^(ph(1) * log10(N2(1)) + ph(2)),10^(ph(1) * log10(N2(5)) + ph(2))];

regf = [10^(pf(1) * log10(N2(1)) + pf(2)),10^(pf(1) * log10(N2(5)) + pf(2))];

regs = [10^(ps(1) * log10(N2(1)) + ps(2)),10^(ps(1) * log10(N2(5)) + ps(2))];


hold on


hold on


title(’Euro-Call Lineare Regressionsanalyse Euler-Maruyama’)



legend(’Halton-Folge’, ’Faure-Folge’,’Sobol-Folge’)

A.5. ASIATISCHE OPTIONEN DISKRET GEOMETRISCH 71

Fur Ito-Formel ist das Programm-Code Analog.

A.5 Asiatische Optionen diskret geometrisch

Programmcode A.5.1 Asia-Call diskret Euler-Maruyama-Ansatz Faure-Folgen

function [V] = qmcacfaure(S0,K,r,Sigma,T,N,M)

h = T / M;

DW = qmcf(N + 1,M);

S = zeros(N,M + 1);

S(:,1) = S0;

Z = zeros(1,N);

for i = 1 : N

for j = 1 : M


* DW(i + 1,j);

end

end

ST = zeros(1,N);

for i = 1:N

ST(i) = (prod(S(i,2 : M + 1)))^(1 / M);

end

for i = 1 : N

Z(i) = ST(i) - K;

if Z(i) < 0

Z(i) = 0;

end

end

V = exp(-0.03) * mean(Z);


Programmcode A.5.2 Asia-Call diskret Ito-Ansatz Faure-Folgen

function [V] = qmcacfaure(S0,K,r,Sigma,T,N,M)

h = T / M;

DW = qmcf(N + 1,M);

S = zeros(N,M + 1);

S(:,1) = S0;

Z = zeros(1,N);

for i = 1 : N

for j = 1 : M


* DW(i + 1,j);

end

end

ST = zeros(1,N);

for i = 1:N

ST(i) = (prod(S(i,2 : M + 1)))^(1 / M);

end

for i = 1 : N

Z(i) = ST(i) - K;

if Z(i) < 0

Z(i) = 0;

end

end

V = exp(-0.03) * mean(Z);

Fur Halton-Folge und Sobol-Folge sind die Programm-Codes Analog.

A.5. ASIATISCHE OPTIONEN DISKRET GEOMETRISCH 73

Programmcode A.5.3 Asia-Call diskret Ito-Ansatz Konvergenzverhaltenformat long;

N = [10000,20000,40000,80000,160000];

M = 12;

N2 = N * M;

n = log10(N2);





rfh(i) = abs(qmcachaltonito(100,103,0.03,0.25,1,N(i),M)

- 5.181006824978432) / 5.181006824978432;

rff(i) = abs(qmcacfaureito(100,103,0.03,0.25,1,N(i),M)

- 5.181006824978432) / 5.181006824978432;

rfs(i) = abs(qmcacsobolito(100,103,0.03,0.25,1,N(i),M)

- 5.181006824978432) / 5.181006824978432;

end


hold on

loglog(N2,rff,’--go’,’LineWidth’,2,’MarkerFaceColor’,’g’,’MarkerSize’,8);

hold on

loglog(N2,rfs,’--bo’,’LineWidth’,2,’MarkerFaceColor’,’b’,’MarkerSize’,8);




legend(’Halton-Folge’,’Faure-Folge’,’Sobol-Folge’);

grid on

Fur Euler-Maruyama-Verfahren ist das Programmcode analog.


Programmcode A.5.4 Asia-Call diskret Ito-Ansatz Lineare Regressionsanalyse

N = [10000,20000,40000,80000,160000];

M = 12;

N2 = N * M;

n = log10(N2);


rfh(i) = abs(qmcachaltonito(100,103,0.03,0.25,1,N(i),M)

- 5.181006824978432) / 5.181006824978432;

rff(i) = abs(qmcacfaureito(100,103,0.03,0.25,1,N(i),M)

- 5.181006824978432) / 5.181006824978432;

rfs(i) = abs(qmcacsobolito(100,103,0.03,0.25,1,N(i),M)

- 5.181006824978432) / 5.181006824978432;

end

rrfh = log10(rfh);


rrff = log10(rff);


rrfs = log10(rfs);


N3 = [N2(1),N2(length(N))];

regh = [10^(ph(1) * log10(N2(1)) + ph(2)),10^(ph(1) * log10(N2(length(N))) + ph(2))];

regf = [10^(pf(1) * log10(N2(1)) + pf(2)),10^(pf(1) * log10(N2(length(N))) + pf(2))];

regs = [10^(ps(1) * log10(N2(1)) + ps(2)),10^(ps(1) * log10(N2(length(N))) + ps(2))];


hold on


hold on


title(’Asia-Call diskret Lineare Regressionsanalyse Ito’)



legend(’Halton-Folge : - 0.5443’, ’Faure-Folge : - 0.5406’,’Sobol-Folge : - 1.2899’)


A.6. ASIATISCHE OPTIONEN KONTINUIERLICH GEOMETRISCH 75

A.6 Asiatische Optionen kontinuierlich geometrisch

Programmcode A.6.1 Asia-Call kontinuierlich Euler-Maruyama-Ansatz Sobol-Folgen

function [V] = qmcacsobol(S0,K,r,Sigma,T,N,M)

h = T / M;

DW = qmcs(N + 1,M);

S = zeros(N,M + 1);

S(:,1) = S0;

Z = zeros(1,N);

for i = 1 : N

for j = 1 : M


* DW(i + 1,j);

end

end

ST = zeros(1,N);

for i = 1 : N

ST(i) = (prod(S(i,2 : M + 1)))^(1 / M);

end

for i = 1 : N

Z(i) = ST(i) - K;

if Z(i) < 0

Z(i) = 0;

end

end

V = exp(-0.03) * mean(Z);


Programmcode A.6.2 Asia-Call kontinuierlich Ito-Ansatz Sobol-Folgen

function [V] = qmcacsobolito(S0,K,r,Sigma,T,N,M)

h = T / M;

DW = qmcs(N + 1,M);

S = zeros(N,M + 1);

S(:,1) = S0;

Z = zeros(1,N);

for i = 1 : N

for j = 1 : M

S(i,j + 1) = S(i,j) * exp((r - Sigma^2 / 2) * h + Sigma * sqrt(h)

* DW(i + 1,j));

end

end

ST = zeros(1,N);

for i = 1 : N

ST(i) = (prod(S(i,2 : M + 1)))^(1 / M);

end

for i = 1 : N

Z(i) = ST(i) - K;

if Z(i) < 0

Z(i) = 0;

end

end

V = exp(-0.03) * mean(Z);

A.6. ASIATISCHE OPTIONEN KONTINUIERLICH GEOMETRISCH 77

Programmcode A.6.3 Asia-Call kontinuierlich Ito-Ansatz Konvergenzverhaltenformat long;

N = [5000,10000,20000,40000,80000];

M = [5,10,20,40,80];

N2 = N.*M;

n = log10(N2);





rfh(i) = abs(qmcachaltonito(100,103,0.03,0.25,1,N(i),M(i))

- 4.769426376039782) / 4.769426376039782;

rff(i) = abs(qmcacfaureito(100,103,0.03,0.25,1,N(1),M(1))

- 4.769426376039782) / 4.769426376039782;

rfs(i) = abs(qmcacsobolito(100,103,0.03,0.25,1,N(i),M(i))

- 4.769426376039782) / 4.769426376039782;

end


hold on

loglog(N2,rff,’--bo’,’LineWidth’,2,’MarkerFaceColor’,’b’,’MarkerSize’,8);

hold on

loglog(N2,rfs,’--go’,’LineWidth’,2,’MarkerFaceColor’,’g’,’MarkerSize’,8);



legend(’Halton-Folge’,’Faure-Folge’,’Sobol-Folge’);

grid on



Programmcode A.6.4 Asia-Call kontinuierlich Ito-Ansatz Lineare Regressionsanalyse

N = [5000,10000,20000,40000,80000];

M = [5,10,20,40,80];

N2 = N.* M;

n = log10(N2);


rfh(i) = abs(qmcachaltonito(100,103,0.03,0.25,1,N(i),M(i))

- 4.769426376039782) / 4.769426376039782;

rff(i) = abs(qmcacfaureito(100,103,0.03,0.25,1,N(1),M(1))

- 4.769426376039782) / 4.769426376039782;

rfs(i) = abs(qmcacsobolito(100,103,0.03,0.25,1,N(i),M(i))

- 4.769426376039782) / 4.769426376039782;

end

rrfh = log10(rfh);


rrff = log10(rff);


rrfs = log10(rfs);


N3 = [N2(1),N2(length(N))];

regh = [10^(ph(1) * log10(N2(1)) + ph(2)),10^(ph(1) * log10(N2(length(N))) + ph(2))];

regf = [10^(pf(1) * log10(N2(1)) + pf(2)),10^(pf(1) * log10(N2(length(N))) + pf(2))];

regs = [10^(ps(1) * log10(N2(1)) + ps(2)),10^(ps(1) * log10(N2(length(N))) + ps(2))];


hold on


hold on


title(’Asia-Call kontinuierlich Lineare Regressionsanalyse Ito’)



legend(’Halton-Folge’, ’Faure-Folge’,’Sobol-Folge’)


Literaturverzeichnis

[1] R. Avikainen: Convergence Rates for Approximations of Functionals of SDEs. 2007.

[2] R. Avikainen: On Irregular Functionals of SDEs and the Euler Scheme. Finance

and Stochastics, Vol. 13, No. 3, 2009.

[3] H. Bauer: Wahrscheinlichkeitstheorie. W.de Gruyter, 1991.

[4] V. Bally und D. Talay: The law of the Euler scheme for stochastic differential

equations, I: convergence rate of the distribution function. Probability Theory and

Related Fields, 1995.

[5] C. Bennet und R. Sharpley: Interpolation of Operators. Academic Press, 1988.

[6] K. Bosch: Elementare Einfuhrung in die Wahrscheinlichkeitsrechnung. Vieweg,

1993.

[7] F. Black und M. Scholes: The pricing of option and corporate liabilities. Journal

of Political Economy, 1973.

[8] S. Cambanis und Y. Hu: Exact Convergence Rate of the Euler-Maruyama Scheme,

with Application to Sampling Design. Stochastics: An International Journal of

Probability and Stochastic Processes, Volume 59, Issue 3 and 4, 1996.

[9] L. Clerck: A Proof of Niederreiter’s conjecture concerning error bounds for quasi-

Monte Carlo integration. 1981.

[10] R. Eller und H. Peter: Derivate und Interne Modelle. Schaefer Poeschel, Stuttgart,

1998.

[11] T. Gerstner: Skript zur Vorlesung: Computational Finance. Goethe-Universitat

Frankfurt, 2010.

79

80 LITERATURVERZEICHNIS

[12] T. Gerstner: Skript zur Vorlesung: Monte Carlo. Goethe-Universitat Frankfurt,

2010.

[13] M. Gunther und A. Jungel: Finanzderivate mit MATLAB. VIEWEG, 2010.

[14] M. B. Giles: Multilevel Monte Carlo path simulation. Operations Research, Vol.

56, No. 3, 2008.

[15] M. B. Giles: Improved multilevel Monte Carlo convergence using the Milstein sche-

me, In: Keller, Alexander; Heinrich, Stefan; Niederreiter, Harald: Monte Carlo and

Quasi-Monte Carlo Methods 2006, Springer, Berlin, 2008.

[16] P. Glassermann: Monte Carlo Methods in Financial Engineering. Springer, New

York, 2004.

[17] J. Guyon: Euler Scheme and Tempered Distributions. Stochastic Processes. No. 6,

2006.

[18] J. Hull: Options, Futures and Other Derivatives. Pearson Prentice Hall, 6th edition,

New Jersey, 2005.

[19] C. Joy, P. Boyle und K. S. Tan: Quasi monte carlo methods in numerical finance.

Managment Science, 1996.

[20] I. Karatzas und S. E. Shreve: Brownian Motion and Stochastic Calculus. Springer,

New York, 1988.

[21] A. Klenke: Wahrscheinlichkeitstheorie. Springer, Berlin, 2006.

[22] P. Kloeden: Skript zur Vorlesung: Numerik stochastischer Differenzialgleichungen.

Goethe-Universitat Frankfurt, 2007.

[23] P. Kloeden und E. Platen: Numerical Solution of Stochastic Differential Equations.

Springer, 2nd ed. corr. print. Berlin, 1995.

[24] J. F. Koksma: Een algemeene telling uit de thoridr gelijkmatge verdeeling modlol.

Mathematica B (Zutphen) (1942/1943).

[25] L. Kuipers und H. Niederetr: Uniform Distribution of Sequence. John Wiley and

Sons NY, 1974.

LITERATURVERZEICHNIS 81

[26] W. J. Morokoff und R. E. Caflisch: Quasi-monte carlo simulation of random walks

in finance. In Monte Carlo and Quasi Monte Carlo Methods. University of Salzburg,

Springer, 1996.

[27] H. Niederreiter: Random Number Gereration and Quasi-Monte-Carlo Methods.

Society for Industrial and Applied Mathematics Philadelphia, Pennsyvania, 1992.

[28] H. Niederreiter: Methods for estimationg discrepancy, in Applications of Number

Theory to Numerical Analysis. S.K Zraemba, ed., Acadamic Press, New Yourk,

1972.

[29] H. Niederreiter: Application of diophantine approximations to numerical integra-

tion. In Diophantine Approximation and it’s Applications, C.F. Osgood, ed., Aca-

demic Press, New Yourk, 1973.

[30] H. Niederreiter und J. M. Wills: Diskrepanz und Distanz von Massen bezuglich

konvexer und Jordanscher Mengen. Matz. Z, 1975.

[31] P. D. Proinov: discrepancy and integration of continuous function. J. Approx.

Theory, 1988.

[32] U. Seydel: Tools for Computational Finance. Springer, 2009.

[33] J. Schumacher: Monte-Carlo- und Quasi-Monte-Carlo-Methoden zur Options-

scheinbewertung. Bonn, 2000.

[34] X. Wang und I. H. Sloan: Efficient Weighted Lattice Rules with Applications to

Finance. 2000.

[35] X. Wang und I. H. Sloan: WHY ARE HIGH-DIMENSIONAL FINANCE PRO-

BLEMS OFTEN OF LOW EFFECTIVES DIMENSION. 2000.

[36] G. A. Willard: Calculating prices and sensitivities for pathindependent derivative

securities in multifactor models. Journal of Derivatives, 1997.

[37] M. Zahrnhofer: Einfuhrung in Quasi-Monte Carlo Verfahren. Graz, 2007.

[38] Zaremba: La discrepance isotrope et integration numerique. Ann.Mat.Pura, 1970.

Quasi-Monte-Carlo-Verfahren zur Bewertung von · PDF fileQuasi-Monte-Carlo-Verfahren zur...

Documents

Transcript of Quasi-Monte-Carlo-Verfahren zur Bewertung von · PDF fileQuasi-Monte-Carlo-Verfahren zur...