Numerische Analyse eines Lucas-Tree

Numerische Analyse eines
Lucas-Tree-Modells mit zwei
heterogenen Investoren
Diplomarbeit
zur Erlangung des akademischen Grades
Diplom-Mathematiker/in
Westfälische Wilhelms-Universität Münster
Fachbereich Mathematik und Informatik
Institut für Numerische und Angewandte Mathematik
Betreuung:
Prof. Dr. Martin Burger
Prof. Dr. Nicole Branger
Eingereicht von:
Jan Philipp Bensch
Münster, Februar 2013
i
Abstract
In this diploma thesis, differential equations for two heterogeneous Epstein-Zin investors are
derived with the help of modern option-price theory. These differential equations describe
the wealth-consumption ratio of the two investors. They are restricted numerically under two
different restrictions. Under the stronger restriction, the equations are simplified to ordinary differential equations. Here, the derivatives of the equations are approximated with the
Upwind-procedure. Then, the ordinary differential equations are solved numerically with the
Newton’s method. Under the weaker restriction, the partial characteristic of the differential
equations is received. In this case, differential equations are dependent on the consumption
share and the long-run risk variable. These partial differential equations are solved numerically with the Forward-Explicit Euler method. The two different solution methods presented
both were tested on their exactness and speed.
ii
Eidesstattliche Erklärung
Hiermit versichere ich, Jan Philipp Bensch, dass ich die vorliegende Arbeit selbstständig
verfasst und keine anderen als die angegebenen Quellen und Hilfsmittel verwendet habe.
Gedanklich, inhaltlich oder wörtlich übernommenes habe ich durch Angabe von Herkunft
und Text oder Anmerkung belegt bzw. kenntlich gemacht. Dies gilt in gleicher Weise für
Bilder, Tabellen, Zeichnungen und Skizzen, die nicht von mir selbst erstellt wurden.
Alle auf der CD beigefügten Programme sind von mir selbst programmiert worden.
Münster, 11. Februar 2013
Vorname Nachname
iii
Danksagung
Ich möchte mich bei einigen Menschen bedanken, die es mir erst ermöglicht haben diese
Diplomarbeit zu erstellen.
• Herrn Professor Dr. Martin Burger, für die gute Begleitung und Beantwortung meiner
Fragen während meiner Arbeit.
• Frau Professorin Dr. Nicole Branger, für die Auswahl und die Einführung in dieses interessante Thema.
• Meiner Freundin Kira Neuschildkamp für die bedingungslose Unterstützung.
• Meiner Schwester Britta Bensch-Horst und ihrem Mann Markus Horst für das Korrekturlesen meiner Arbeit.
• Meinen Eltern Margot und Peter Bensch für die finanzielle und moralische Unterstützung
während des gesamten Studiums.
iv
Inhaltsverzeichnis
1 Einleitung
1
2 Grundlagen
3
2.1
2.2
Betriebswirtschaftliche Grundlagen . . . . . . . . . . . . . . . . . . . . . . . .
3
2.1.1
Optionstheorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.1.2
Marktannahmen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
Mathematische Grundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2.1
Wahrscheinlichkeitstheorie . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2.2
Numerik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
3 Modellherleitung
16
3.1
Heterogene Investoren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
3.2
Konsumanteil und Wohlstandskonsumrate . . . . . . . . . . . . . . . . . . . .
18
3.3
Differentialgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
3.4
Pricing Kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
3.5
Drift und Volatilität . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3.6
Übersicht aller relevanten Gleichungen . . . . . . . . . . . . . . . . . . . . . .
23
4 Wohlstandskonsumquotient als gewöhnliche Differentialgleichung
25
4.1
Gewöhnliche Differentialgleichung . . . . . . . . . . . . . . . . . . . . . . . . .
25
4.2
Randwerte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
4.3
Numerische Berechnung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
5 Wohlstandskonsumquotient als partielle Differentialgleichung
36
5.1
Strukturbetrachtung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
5.2
Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
5.3
Randwerte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
5.4
Ableitungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
5.5
Numerische Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
6 Fazit und Ausblick
57
Abbildungsverzeichnis
59
Tabellenverzeichnis
60
Inhaltsverzeichnis
Literaturverzeichnis
v
61
1
1 Einleitung
In dieser Diplomarbeit wird ein Lucas-Tree-Modell mit zwei heterogenen Investoren numerisch analysiert. Das Lucas-Tree-Modell liefert die Wohlstandskonsumquotientenfunktion als
partielle Differentialgleichung. Diese Gleichung wird untersucht und numerisch berechnet.
Der wichtigste Grundbaustein des Lucas-Tree-Modells ist das 1973 von Black and Scholes
vorgestellte Black-Scholes-Modell. Es gilt als bedeutender Meilenstein in der Finanzwirtschaft und ist die Grundlage der meisten modernen Optionsbewertungen.
Warum ist es sinnvoll heterogene Investoren zu betrachten?
Früher wurde in der Finanzmathematik der Schwerpunkt auf homogene Investoren gelegt.
Modelle spiegeln immer eine Vereinfachung der Wirklichkeit da. Gründe, warum man mit
Vereinfachungen arbeitet, gibt es viele. Unter anderem ist es einfacher (teilweise sogar erst
möglich) das Modell zu berechnen oder es gibt zu vernachlässigende Faktoren. Aber es gibt
auch gute Argumente mit möglichst genauen Modellen zu rechnen - das wichtigste Argument
ist die Realitätsnähe. In der heutigen Zeit rückt die Wichtigkeit der Realität in Modellen immer mehr in den Fokus. So beschäftigen sich immer mehr Modelle mit der Heterogenität der
Investoren. Dadurch muss zwar mehr Arbeit aufgrund der erhöhten Komplexität des Modells
investiert werden, aber dafür liefern diese Modelle genauere Resultate. Heterogene Investoren
unterscheiden sich in ihrer Risikoaversion (risk aversion), in ihrer Zeitpräferenz (time preference) und ihrer Substitutionselastizität (elasticity of intertemporal substituion) voneinander.
Warum werden die Gleichungen numerisch analysiert?
In dieser Arbeit werden mehrere Differentialgleichungen betrachtet. Deren Lösung lassen sich
häufig nur schlecht oder sogar gar nicht analytisch berechnen. Deswegen ist es sinnvoller die
Differentialgleichungen numerisch zu berechen. Dieses anvisierte Ziel (numerisch statt analytisch) ist mit deutlich weniger Rechenaufwand verbunden und kann so einfacher berechnet
werden. Beachtet werden muss aber, dass die numerische Lösung konvergiert und auch stabil
ist.
Diese Diplomarbeit beginnt mit der Einführung der betriebswirtschaftlichen Grundlagen in
Kapitel 2. Dort folgt die Definition des Optionsbegriffs, der die Grundlage für moderne Modelle wie der Epstein-Zin-Nutzenfunktion ist. Weiterhin folgt die Definition für die gegebenen
2
Marktannahmen für diese Diplomarbeit. Im zweiten Teil des Kapitels werden die mathematischen Hilfsmittel der Wahrscheinlichkeitstheorie und Numerik eingeführt. Dabei werden
die Definition eines Wiener Prozesses und das Lemma von Ito wesentliche Bestandteile des
Bereichs Wahrscheinlichkeitstheorie werden. Für die Numerik werden vor allem verschiedene Differenzenquotienten und die Konvergenz von Differentialoperatoren vorgestellt. In
Kapitel 3 folgt die Herleitung der partiellen Differentialgleichung mit Hilfe der Epstein-ZinNutzenfunktion. Diese Gleichungen werden dann unter ein paar Einschränkungen in Kapitel
4 numerisch betrachtet und anschließend für ein konkretes Beispiel berechnet. Unter diesen
Einschränkungen handelt es sich dann um gewöhnliche Differentialgleichungen. Im 5. Kapitel
erfolgt die Betrachtung einer nicht so starken Einschränkung der Differentialgleichungen, so
dass durch die Long-Run-Risk-Variable der partielle Teil der Differentialgleichung erhalten
bleibt. Es folgt eine mathematische Einordnung dieser partiellen Differentialgleichungen und
die Berechnung mit Hilfe des Vorwärts-Euler-Verfahrens. Zum Abschluss der Diplomarbeit
wird ein Fazit über die Differentialgleichungen und ihrer Lösungsverfahren gezogen und Ideen
für mögliche weitere Untersuchungen diskutiert.
3
2 Grundlagen
Im ersten Teil dieses Kapitels erfolgt die Einführung der benötigten grundlegenden Begriffe der Betriebswirtschaftslehre, wobei der Schwerpunkt auf der Optionspreistheorie liegt.
Dabei wird erklärt, was man unter einer Option versteht und welche verschiedenen Optionsarten es gibt. Für tiefere Einblicke in dieses Themengebiet empfiehlt sich Optionen, Futures
und andere Derviate von Hull [Hul09]. Es folgt die Definition der Marktannahmen, die die
Grundlage für das Modell bilden. Im zweiten Teil des Kapitels werden die wichtigen Grundbegriffe und Hilfsmittel der Mathematik eingeführt. Dieses sind Begriffe und Hilfsmittel aus
der Wahrscheinlichkeitstheorie und der Numerik. Der Aufbau ist dabei ähnlich wie in den
Diplomarbeiten [Hä12] und [Dö12].
2.1 Betriebswirtschaftliche Grundlagen
Nachfolgend werden die wichtigsten betriebwirtschaftlichen Grundlagen eingeführt, die für
diese Diplomarbeit benötigt werden. Der Aufbau ist dabei nahe an die Diplomarbeit [Zoc10]
angelehnt. Als erstes wird der Optionsbegriff definiert und der Unterschied zwischen einer
Call- und einer Put-Option erklärt. Anschließend werden noch die Marktannahmen für das
Modell vorgestellt.
2.1.1 Optionstheorie
Der Optionsbegriff
Unter einer Option versteht man das Recht, zu einem späteren Zeitpunkt ein Gut (z.B.: eine
Aktie) zu einem vorher vereinbarten Preis zu kaufen oder zu verkaufen. Dabei unterscheidet
man zwischen verschiedenen Arten von Optionen.
Falls der Optionskäufer das Recht erwirbt zu einem späteren Zeitpunkt eine festgelegte Menge zu kaufen, spricht man von einer Call-Option. Eine Put-Option ist hingegen, das Recht zu
4
einem späteren Zeitpunkt eine festgelegte Menge des Gutes zu verkaufen. Es gibt noch weitere zu unterscheidende Optionen. Zum Beispiel hat der Optionskäufer bei der europäischen
Option nur die Möglichkeit, das Gut am Ende der vorher festgelegten Laufzeit zu kaufen (CallOption) oder zu verkaufen (Put-Option). Bei der amerikanischen Option hingegen kann er
während der gesamten Laufzeit kaufen bzw. verkaufen. Falls es die Möglichkeit gibt, die Option während mehrerer vorher festgelegten Zeiträumen bzw. Zeitpunkten auszuüben, spricht
man von der Bermuda-Option. In dieser Arbeit wird der Schwerpunkt auf die europäische
Call- bzw. Put-Option gelegt. Falls nicht anders erwähnt, wird ab jetzt bei einer Option immer von einer europäischen Option ausgegangen.
Beispiel: Eine europäische Call-Option
Eine Aktie hat zum Zeitpunkt t = 0 den Preis 20 ¤. Ein Spekulant hat die Möglichkeit für
je 0,50 ¤ eine Option auf den Kauf einer Aktie zum Preis von 25 ¤ in einem halben Jahr
(t = 0, 5) zu erwerben.
Da der Optionskäufer nur das Recht, aber nicht die Pflicht erwirbt, die Option auszuüben,
führt er sie nur dann aus, wenn er einen Vorteil durch die Ausübung erlangt. Dafür ist es
entscheidend, wie hoch der Marktwert des Produkts zum Fälligkeitszeitpunkt im Vergleich
zum Kaufpreis ist.
Zurück zu dem Beispiel:
Der Spekulant erwirbt 1000 Optionen zum Preis von 500 ¤
1. Fall: Die Aktie hat zum Fälligkeitszeitpunkt (t = 0, 5) einen Wert von 35 ¤. Der Gewinn bzw. Verlust (B) beim Ausüben der Option berechnet sich wie folgt:
B = 1000 · (35 − 25) − 500 = 9500
Also wird der Spekulant die Option ausüben, um die Aktien am Markt direkt wieder zu verkaufen und macht einen Gewinn von 9500 ¤. Da der Preis der Aktie zum Fälligkeitszeitpunkt
über dem Ausübungspreis (35>25) liegt, spricht man davon, dass ”die Option sich im Geld
(in the money)”befindet.
2. Fall: Die Aktie hat zum Fälligkeitszeitpunkt (t = 0, 5) einen Wert von 25 ¤. Der Gewinn bzw. Verlust (B) beim Ausüben der Option berechnet sich wie folgt
2
Grundlagen
5
Call-Option
Put-Option
St > K
in the money
→ kaufen
out of money
→ nicht verkaufen
St = K
at the money
at the money
St < K
out of money
→ nicht kaufen
in the money
→ verkaufen
Tabelle 2.1: Call-/Put-Option - Vergleich Preis des Gutes und festgelegter Ausübungspreis
B = 1000 · (25 − 25) − 500 = −500
Also ist es für den Spekulanten egal ob er die Option ausübt, da er auch im Falle der
Nichtausübung einen Verlust von 500 ¤ (Kosten für die Option) hat. Da der Ausübungspreis
gleich dem Aktienpreis (25=25) ist, wird hier von ”die Option liegt am Geld (at the money)”gesprochen.
3. Fall: Die Aktie hat zum Fälligkeitszeitpunkt (t = 0, 5) einen Preis von 23 ¤. Falls der
Spekulant die Option ausübt, würde sich der Verlust (B) wie folgt berechnen:
B = 1000 · (23 − 25) − 500 = −2500
In diesem Fall macht es für den Spekulant natürlich keinen Sinn die Option auszuüben, da er
beim Verfallenlassen der Option nur einen Verlust von 500 ¤ verkraften muss. Wenn, wie in
diesem Fall, der Ausübungspreis über dem Aktienpreis zum Fälligkeitszeitpunkt liegt, spricht
man davon, dass ”die Option sich aus dem Geld (out of money)”befindet.
Es muss beachtet werden, dass dieses Beispiel eine Call-Option wiedergibt. Bei einer PutOption gestaltet es sich gegensätzlich wie die Tabelle 2.1 zeigt. Dafür sei St der Preis des
Gutes zum Fälligkeitszeitpunkt und K der laut Option festgelegte Ausübungspreis.
6
Abbildung 2.1: Auszahlungsfunktion einer Call- und einer Put-Option mit Ausübungspreis
K=100
.
Die Auszahlungsfunktionen der beiden Optionsalternativen sind definiert als:
Call : C(S, t) = (St − K)+ := max{0, St − K}
P ut : P (S, t) = (St − K)+ := max{0, K − St }
In der Abbildung 4.1 sieht man die Auszahlungsfunktionen für eine Call- und eine PutOption mit Ausübungspreis K=100. Wie zu sehen ist, sind die Auszahlungsfunktionen stets
nicht negativ. Da allerdings der Optionskäufer den Optionspreis bezahlen muss, kann ein
Optionskäufer trotzdem einen Verlust erwirtschaften (siehe Beispiel).
2
Grundlagen
7
2.1.2 Marktannahmen
Bei diesem Modell wird von einem vollkommenen Markt ausgegangen. Dies bedeutet:
• alle Marktteilnehmer verhalten sich rational und nutzenmaximierend
• keine zeitliche Verzögerung (alle Marktteilnehmer reagieren bei Veränderung der Marktbedingungen sofort)
• keine Steuern, Gebühren oder Transportkosten
• vollkommene Markttransparenz (alle Informationen sind von allen Marktteilnehmern einzusehen)
• einheitlicher Marktzins (Habenzinssatz = Sollzinssatz)
• keine Arbitragemöglichkeiten (es gibt also nicht die Möglichkeit gleiche Güter auf unterschiedlichen Märkten zu unterschiedlichen Preisen zu kaufen oder verkaufen)
Diese Annahmen entsprechen nicht der wirtschaftlichen Realität, bilden aber die Grundlage
vieler Modelle. Häufig können von einem vollkommenen Markt Aussagen auf einen unvollkommenen und damit realistischeren Markt übertragen werden.
8
2.2 Mathematische Grundlagen
Für diese Diplomarbeit werden Grundlagen aus der Wahrscheinlichkeitstheorie und der Numerischen Mathematik benötigt. Die aus betriebswirtschaftlichen Aspekten verwendeten mathematischen Definitionen kann man gut in Financial asset pricing theory von Munk [Mun09]
nachlesen. Tiefere Einblicke in die Wahrscheinlichkeitstheorie bzw. mathematischen Numerik
sind nachzulesen in Wahrscheinlichkeitstheorie von Klenke [Kle08] und in Numerik partieller
Differentialgleichungen von[Bur07]. Zu Beginn werden die wichtigsten Grundlagen zur Wahrscheinlichkeitstheorie dargestellt.
2.2.1 Wahrscheinlichkeitstheorie
Ziel diese Abschnitts wird es sein, einen Wiener Prozess zu definieren und das Lemma von Ito
einzuführen. Für einen Wiener Prozess, der auch Brownsche Bewegung genannt wird, werden
einige grundlegende Begriffe aus der Wahrscheinlichkeitstheorie benötigt. Die Definitionen
und Aussagen von [Mun09] bieten Orientierung bei der Einführung der benötigten Begriffe
der Wahrscheinlichkeitstheorie. Für den Begriff eines Wahrscheinlichkeitsraumes, werden die
Definitionen für einen Ergebnisraum Ω, eine σ-Algebra und eines Wahrscheinlichkeitsmaßes
µ benötigt.
Definition 2.2.1 (Ergebnisraum)
Der Ergebnisraum Ω ist diejenige Menge, die alle möglichen Ergebnisse enthält. Ein Ereignis
ω ist eine Teilmenge von Ω.
Definition 2.2.2 (σ–Algebra)
F ist eine σ - Algebra in Ω, falls F eine Menge von Teilmengen der Grundmenge Ω ist mit
folgenden Eigenschaften:
•Ω∈F
• Für jede Menge F in F gilt, das Kompliment F c ≡ Ω \ F ist wieder in F
S
• Falls F1 , F2 , . . . ∈ F, dann gilt auch für die Vereinigung ∞
n=1 Fn ∈ F
Definition 2.2.3 (Wahrscheinlichkeitsmaß)
P ist ein Wahrscheinlichkeitsmaß, falls P eine Funktion P : F → [0, 1] ∈ (R) mit P(Ω) = 1
S
P∞
ist und für jede Folge (A1 , A2 , . . . ) disjunkter Mengen gilt: P( ∞
n=1 P(An ).
n=1 An ) =
2
Grundlagen
9
Mit diesen drei Definitionen kann ein Wahrscheinlichkeitsraum definiert werden, der die
Grundlage für den Wiener Prozess bietet.
Definition 2.2.4 (Wahrscheinlichkeitsraum)
Falls Ω ein Ergebnisraum, F eine σ-Algebra auf Ω und P ein Wahrscheinlichkeitsmaß ist,
wird das Triple (Ω, F, P) als Wahrscheinlichkeitsraum bezeichnet.
Auf diesem Wahrscheinlichkeitsraum kann eine Zufallsvariable definiert werden. Zufallsvariablen dienen in der Wahrscheinlichkeitstheorie dazu, Ergebnissen aus Zufallsexperimenten
bestimmte Werte zuzuordnen. Formal definiert sind sie al
Definition 2.2.5 (Zufallsvariable)
Sei (Ω0 , F 0 ) ein Messraum und (Ω, F, P) ein Wahrscheinlichkeitsraum. Falls X : Ω → Ω0 eine
F-messbare Funktion ist, so heißt X eine Zufallsvariable.
Definition 2.2.6 (Filtrierung)
Es sei (Ω, F, P) ein Wahrscheinlichkeitsraum und I eine vollständig geordnete Indexmenge.
Sei Ft für t ∈ I eine Folge von Teil-σ-Algebren für die gilt, Ft ⊆ F für alle t ∈ I.
(Ft )t ∈ I heißt Filtration, wenn aus s ≤ t folgt Fs ⊆ Ft .
Definition (Dichtefunktion)
Sei P ein Wahrscheinlichkeitsmaß und X eine reellwertige Zufallsvariable. Sei f : R → [0, ∞)
eine Funktion für die gilt
Z
P(a ≤ X ≤ b) =
b
f (x)dx
a
für alle reelle Zahlen a < b so heißt f Wahrscheinlichkeitsdichte der Zufallsvariable X.
Als bekanntes Beispiel für ein Wahrscheinlichkeitsmaß definieren wir die Normalverteilung.
Definition 2.2.7 (Normalverteilung)
Sei X eine stetige Zufallsvariable mit Wahrscheinlichkeitsdichte f : R → R gegeben durch
1 x−µ 2
1
f (x) = √ e− 2 ( σ )
σ 2π
so heißt X ∼ N (µ, σ 2 )- normalverteilt mit den Parametern µ (Erwartungswert) und σ (Varianz).
10
Abbildung 2.2: Der Pfad eines Wiener Prozesses Zt in der Zeitspanne von t=0 bis t=1.
Definition 2.2.8 (Wiener Prozess/Brownsche Bewegung)
Ein stochastischer Prozess Z auf dem Wahrscheinlichkeitsraum (ω, A, P ) heißt Wiener Prozess
oder auch Brownsche Bewegung, falls er folgende Eigenschaften erfüllt:
(i)Z0 = 0
(ii)für gegebene Zeitpunkte t, t0 ≥ 0 mit t < t0 gilt Zt0 –Zt ∼ N (0, t0 − t)
(iii)für alle gegebenen Zeitpunkte 0 ≤ t0 < t1 , . . . < tn , sind die Zuwächse
Zt1 –Zt0 , . . . , Ztn − Ztn−1 stochastisch unabhängig
(iv)z hat stetige Pfade.
Die Darstellung eines Pfades in Abbildung 2.2 dient als Beispiel für einen Wiener Prozess.
2
Grundlagen
11
Für die Herleitung unserer Differentialgleichungen in Kapitel 3 wird das Lemma von Ito eine
sehr wichtige Rolle spielen. Dafür wird als erstes der stochastische Prozess definiert, da er für
den Ito-Prozess und damit das Lemma von Ito benötigt wird.
Definition 2.2.9 (stochastischer Prozess)
Sei (Ω, F, P) ein Wahrscheinlichkeitsraum. Eine Familie von Zufallsvariablen Xt auf diesem
Wahrscheinlichkeitsraum wird als stochastischer Prozess bezeichnet.
Definition 2.2.10 (Ito-Prozess)
Ein stochastischer Prozess X = (Xt )t≥0 heißt Ito-Prozess, falls seine Veränderung für ein
kleines Zeitintervall [t, t + dt] geschrieben werden kann, als
dXt = µt dt + σt dWt ,
wobei W eine Brownsche Bewegung ist und der Drift µ und die Volatilität σ selbst stochastische Grössen sind. Falls µt und σt als Funktionen von t und Xt geschrieben werden können,
wird von einem Diffusionsprozess gesprochen.
Zum Abschluss dieses Abschnitts führen wir das Lemma von Ito ein. Dieses entspricht der
Kettenregel für stochastische Prozesse.
Satz 2.2.11 (Lemma von Ito) [Bur10]
Sei X = (Xt )t≥0 ein reellwertiger Ito-Prozess mit Dynamiken
dXt = µt dt + σt dWt ,
wobei µ und σ reellwertige Prozesse und W ein Wiener Prozess ist. Sei g(X, t) eine reellwertige
in X zweimal und in t einmal stetig differenzierbare Funktion, dann ist der Prozess y(yt )t≥0
definiert als
yt = g(Xt , t)
ein Ito-Prozess mit Dynamiken
dyt =
∂g
∂g
1 ∂2g
∂g
2
(xt , t) +
(Xt , t)µt +
(Xt , t)σt dt +
(Xt , t)σt dWt .
∂t
∂X
2 ∂X 2
∂X
12
2.2.2 Numerik
Bei den Hilfsmitteln der Numerik dienen uns die Skripte Numerik partieller Differentialgleichungen [Bur07] und Mathematische Modellierung [Bur10] als Orientierung. Für diese Arbeit
benötigen wir verschiedene Formen der Differenzenquotienten, deren Einführung in diesem
Kapitel erfolgten Danach werden wir noch die Begriffe der Konsistenz, Stabilität und Konvergenz betrachten. Zuerst führen wir die Diskretisierung ein.
Diskretisierung
Unter einer Diskretisierung versteht man ein Gitter auf dem man eine Funktionen definieren
möchte. Das Gitter kann je nach Dimension der Funktion definiert werden. Unter einem regulären Gitter versteht man ein Gitter dessen Schrittweiten konstant sind.
Differenzenquotienten
Differenzenquotienten dienen der Approximation von sowohl gewöhnlichen als auch partiellen Ableitungen. Es gibt viele verschiedene Differenzenquotienten mit unterschiedlichen Vorund Nachteilen. Manche haben eine höhere Konsistenzordnung (schneller), manche sind dafür
weniger anfällig für Fehler (stabiler). Wir führen den Vorwärts-, Rückwärts- und zentralen
Differenzenquotienten [Bur07, S.11ff] ein. Dabei sei u(t) die zu betrachtende Funktion. h
beschreibt dabei die Schrittweite für die Differenzenquotienten. uj bezeichnet den Wert der
Funktion u an der Gitterstelle j.
Beim Vorwärtsdifferenzenquotienten gilt
du
1
= − (uj+1 − uj ),
dt
h
(2.1)
beim Rückwärtsdifferenzenquotienten
du
1
= − (uj − uj−1 )
dt
h
(2.2)
und beim zentralen Differenzenquotienten
du
1
= − (uj+1 − uj−1 ).
dt
2h
(2.3)
2
Grundlagen
13
Der zentrale Differenzenquotient hat von diesen dreien mit 2 die höchste Konsistenzordnung.
Je größer die Konsistenzordnung ist, desto größere Schritte können beim Verfahren gewählt
werden. Für die Definition der Konsistenz siehe (2.2.13). Da der zentrale Differenzenquotient
und auch der Vorwärtsdifferenzenquotient aber aufgrund eines wachsenden Fehlers instabil
werden können, wäre der Rückwärtsdifferenzenquotient häufig die bessere Wahl. Da, wie in
Kapitel 4 zu sehen ist, je nach Situation nur uj−1 oder uj+1 zur Verfügung steht, wird eine
Vereinigung der ersten beiden Quotienten eingeführt - das so genannte Upwind-Verfahren.
Dieses Verfahren ist auch im Gegensatz zum zentralen Differenzenquotienten nicht so anfällig
gegen konvektionsdominante Probleme. Das Upwind-Verfahren nutzt, je nach Vorfaktor
von
du
dt ,
den Vorwärts- oder Rückwärtsdifferenzenquotienten. Falls in der betrachteten Glei-
chung a ·
du
dt
steht, gilt für
a·
uj − uj−1
uj+1 − uj
du
= −max(a, 0) ·
− min(a, 0) ·
.
dt
h
h
(2.4)
Falls in der zu untersuchenden Gleichung Ableitungen 2. Ordnung betrachtet werden, wird
die natürliche Approximation [Bur07, S.16] genutzt. Sie approximiert die 2. Ableitung an drei
Gitterpunkten.
uj+1 − 2uj + uj+1
d2 u
=
2
dt
h2
(2.5)
Diese hat als Konsistenzordnung ebenfalls 2.
Alle diese Quotienten lassen sich auch auf partielle Ableitungen übertragen, indem sie jeweils
für die Ableitung nach einer Variablen gelten. Allerdings können bei partiellen Ableitungen
auch Ableitungen in mehrere verschiedene Richtungen vorkommen. Da in dieser Arbeit maximal zwei verschiedene Ableitungsrichtungen vorkommen können, erfolgt die Betrachtung
beschränkt auf diesen Fall. Dafür wird die Fünf-Punkte Stern Approximation [Bur07, S.23]
betrachtet:
4ui,j − ui,j+1 − ui+1,j − ui,j−1 − ui−1,j
d2 u
=
,
dt1 dt2
h1 h2
(2.6)
wobei h1 und h2 die Schrittweiten nach t1 bzw. t2 sind.
Wir werden in Kapitel 5 darauf eingehen, dass am Rand von Definitionsbereichen spezielle
Differenzenquotienten gewählt werden sollten. Es sollte darauf geachtet werden, mit den Differenzenquotienten nur Punkte anzusprechen die innerhalb des Gitters liegen. Zum Abschluss
dieses Kapitels werden wir noch die Begriffe der Konsistenz, Stabilität und der Konvergenz
einführen.
14
Definition 2.2.12 (Diskrete Stabilität)
Sei Lh : Gh → RN die diskrete Approximation eines Differentialoperators. Dann heißt Lh
−1
diskret stabil, wenn L−1
h existiert, für h > 0 hinreichend klein und kLh k gleichmäßig in h
beschränkt ist.
Definition 2.2.13 (Diskrete Konsistenz)
Sei L : C k (Ω) → C 0 (Ω) ein Differentialoperator der Ordnung k und Lh eine diskrete Approximation auf einem Gitter Gh . Die Approximation heisst diskret konsistent, falls
kLh (u|Gh ) − (Lu)|Gh kh → 0
gilt. Die Konsistenzordnung der Approximation ist m, falls kLh (u|Gh ) − (Lu)|Gh kh ≤ Chm
für alle u ∈ C k+m (Ω) gilt.
Aus Konsistenz und Stabilität folgt die Konvergenz, welche im folgenden Satz zum Ausdruck gebracht wird. Für den Beweis und weiteres zum Thema wird verwiesen auf das Skript
Numerik partieller Differentialgleichungen [Bur07].
Satz 2.2.14 (Satz zur Konvergenz) [Bur07, S.18]
e h : C ( Ω) → C 0 (Ω) eine stabile und konsistente Approximation eines DifferentialopeSei L
rators L : C k (Ω) → C 0 (Ω). Sei u die Lösung der Differentialgleichung Lu = f und u
eh die
eh u
Lösung von L
eh = feh , so dass feh → f für h → 0. Dann ist die Approximation konvergent,
d.h. u
eh → u für h → 0.
Zum Abschluss beschäftigen wir uns in diesem Kapitel mit der Betrachtung des NewtonVerfahrens. Es ist ein Beispiel für ein einfaches Dikretisierungsverfahren zur Bestimmung
von Nullstellen einer Funktion [KM74, S.95ff].
Verfahren 2.2.15 (Newton-Verfahren)
Sei x0 der vorgegebene Startwert (Näherung) und n=0,1,2,... die Schrittzahl, dann gilt
xn+1 = xn −
f (xn )
f 0 (xn )
konvergiert für geeigneten Startwert quadratisch gegen eine Nullstelle von f (x).
2
Grundlagen
15
Satz 2.2.16 (Konvergenzsatz für das Newton-Verfahren)
Sei f (x) eine zweimal differenzierbare Funktion und g(x) als g(x) = x −
f (x)
f 0 (x)
definiert. Das
Newton-Verfahren konvergiert in der Umgebung einfacher Nullstellen quadratisch, falls
f (x)f 00 (x) ≤ k < 1.
|g (x)| = [f 0 (x)]2 0
16
3 Modellherleitung
In diesem Kapitel werden wir die zu untersuchenden Differentialgleichungen herleiten. Beim
Herleiten des Modells benutzen wir das Working Paper Preference Heterogeneity and Survival in Long-Run Risk Models von Branger, Dumitrescu, Ivanova und Schlag [BDIS12] und
Explainig Asset Pricing Puzzles Associated with the 1987 Market Crash von Benzone, CollinDufrens und Goldstein [BCDG10]. Die zu bestimmende Funktion ist dabei vi,t (wt , Xt ). Sie
hängt sowohl vom Konsumanteil ωi des Investors i ab, als auch von der Long-Run-RiskVariable X.
Wir nehmen in unserem Modell an, dass zwei Investoren den kompletten Konsum des Marktes
unter sich aufteilen. Deshalb können wir die Gleichung ω1 + ω2 ≡ 1 aufstellen. Der Konsum
des ersten Investors wird als C1 und der des 2. Investors als C2 bezeichnet. Das heißt, wir
können ω1 =
C1
C1 +C2
setzen. Aus Einfachheitsgründen bezeichnen wir ab hier ω ≡ ω1 .
3.1 Heterogene Investoren
Früher wurde bei Investoren immer eine Abhängigkeit der investorenspezifischen Faktoren
Risikoaversion γ und Substitutionsfaktor ψ betrachtet. Es galt die Formel γ =
1
ψ.
Inzwischen
wird immer mehr Wert darauf gelegt, diese beiden Faktoren voneinander zu trennen und sie
als unabhängige Parameter zu betrachten.
Der Risikoaversionsfaktor γ (relative risk aversion) macht deutlich, wieviel Risiko der Investor
eingehen möchte. Je höher er ist, desto höheres Risiko ist der Investor bereit einzugehen. Der
Risikoaversionsfaktor liegt bei unserem Modell im Bereich von [5,10]. Der Substitutionsfaktor
ψ (elasticity of intertemporal substitution EIS) beschreibt, inwieweit der Investor bereit ist,
ein Produkt durch ein anderes zu ersetzen. Je größer dieser Faktor ist, desto eher wäre der
Investor bereit ein Produkt durch ein anderes zu ersetzen. Ein realistischer Wert liegt hierbei
im Bereich [0,2].
Da wir mit ω den Konsumanteil von Investor 1 und mit (1 − ω) den Konsumanteil von
Investor 2 bezeichen können, wird die durchschnittliche Risikoaversion der beiden Investoren
3
Modellherleitung
17
γ als gewichtetes harmonisches Mittel definiert.
γ≡
ω γ11
1
+ (1 − ω) γ12
Der durchschnittliche Substitutionsfaktor ψ ist als arithmetisches Mittel definiert:
ψ ≡ ωψ1 + (1 − ω)ψ2
Neben den beiden Parametern γ und ψ gibt es noch einen dritten Parameter, in dem sich
die beiden Investoren voneinander unterscheiden können, den Zeitpräferenzfaktor β. Er zeigt,
inwiefern der Investor seinen Nutzen lieber früher oder später bekommen möchte. Je größer
der Parameter, desto wichtiger ist es dem Investor den Nutzen möglichst früh zu erhalten.
Ein typischer Wert für β liegt bei ca. 0,02.
Da, wie bereits erwähnt, die Unabhängigkeit der beiden Faktoren γ und ψ benutzt wird
und wir zwei heterogene Investoren betrachten, gilt für die beiden Investoren der rekrusive
Nutzen Ut des Epstein-Zin Typs [BCDG10, S.7]
1
Ut = (1 −
1− 1
e−βdt )Ct ψ
+ e−βdt Et (Ut+dt )1−γ
1− ψ 1−1 1
1−γ
ψ
.
Wie man im Paper Explaining Asset Pricing Puzzles Associated with the 1987 Market Crash
[BCDG10] nachlesen kann, kann durch Umformen und einer Taylorentwicklung der normalisierte Nutzen-Index Jt
Jt =
1 1−γ
U
1−γ t
definiert werden. Außerdem können wir mit Hilfe der normalisierten Sammlerfunktion
Z ∞
Jt = Et
f (Cs , Js )ds
t
den Nutzen-Index für beide Investoren herleiten
f (C, J) =
wobei θi =
1−γi
1− ψ1
i
gilt.
βC
1
1− ψ
1
(1 − ψ1 )[(1 − γ)J] θ −1
− βθJ ,
(3.1)
18
3.2 Konsumanteil und Wohlstandskonsumrate
Wie bereits oben erwähnt, stellt ω den Konsumanteil von Investor 1 dar (1 − ω den Konsumanteil von Investor 2). Die Dynamiken werden durch den stochastischen Prozess
dωt = µω (ωt , Xt )dt + σω (ωt , Xt )0 dWt
beschrieben, wobei als Notation µω ≡ µω (ωt , Xt ) für den Drift und σω ≡ σω (ωt , Xt ) für die
Volatilität vereinbart wird. Sowohl der Drift als auch die Volatilität werden später in den zu
untersuchenden Differentialgleichungen gelöst werden. Die Dynamiken des Konsums Ci sind
definiert durch
dCi,t
0
dWt
= µC1 dt + σC
1
Ci,t
(3.2)
wobei hier nach dem Lemma von Ito [2.2.11] für den Drift
µC1 = µC + Xt +
1
1
µω + σω0 σC
ω
ω
(3.3)
und die Volatilität
σ C1 = σ C +
1
σω
ω
(3.4)
gilt. vi,t ist unsere gesuchte Funktion. Sie stellt die Wohlstandskonsumrate von Investor i zum
Zeitpunkt t da. Sie ist ein stochastischer Prozess mit Dynamiken
0
dvi,t = µvi (ω1,t , Xt )dt + σvi (ωi,t,X
dWt ).
t
(3.5)
Auch hier wird von jetzt an aus Übersichtsgründen für µvi ≡ µvi (ω1,t , Xt ) und σvi ≡ σvi (ωi,t , Xt )
geschrieben. Der Drift und die Volatilität können wie folgt definiert werden, wobei ρωx der
Korrelationsfaktor der beiden Ableitungsrichtungen ist. Für ihn gilt −1 < ρωx < 1.
µvi =
∂vi
∂vi
∂ 2 vi
∂ 2 vi
∂ 2 vi 0
µω −
κXt + 0.5 2 σω0 σω + 0.5 2 σx0 σx + ρωx
σ σx
∂ω
∂x
∂ω
∂x
∂ω∂x ω
(3.6)
∂vi
∂vi
σω +
σx
∂ω
∂x
(3.7)
σvi =
3
Modellherleitung
19
3.3 Differentialgleichung
Die wichtigsten Gleichungen in dieser Diplomarbeit werden die nun hergeleiteten (partiellen)
Differentialgleichungen. Wir werden sie später unter zwei verschiedenen Einschränkungen
analysieren und numerisch berechnen. Wie in Kapitel 2.1.2 schon erwähnt, befinden wir uns
bei unserem Modell in einem vollkommenen Markt. Das heißt unter anderem - dass jeder
Marktteilnehmer - in unserem Fall Investoren, seinen Nutzen bei gegebenem Preis maximiert. Das heißt, wir können annehmen, dass
Et [dJi,t + fi (Ci,t , Ji,t )dt] = 0
(3.8)
gilt (siehe dazu [BCDG10]), wobei Ji,t definiert ist als
Ji,t =
1−γi
Ci,t
1 − γi
βiθi eθi vi,t .
(3.9)
Nun müssen wir noch die erhaltenen Gleichungen zusammenfügen, um die gesuchten Differentialgleichungen zu erhalten. Wir beginnen mit (3.8) und orientieren uns an [Hä12].
Et [dJi,t + fi (Ci,t , Ji,t )dt] = 0
dJi,t fi (Ci,t , Ji,t )dt
⇔ Et
+
=0
J
J
In diese Gleichungen setzen wir (3.1) und (3.9) ein und erhalten nach Vereinfachung
1−γ

Et
Ct i θi θi vi,t
d( 1−γ
βi e
)
i

J

+ e−vi.t θi − βi θi dt = 0.
Künstlich wird eln eingeführt und wir erhalten

Et 
β
θi
i
d( 1−γ
e(1−γi )ln(Ct )+θi vi,t )
i
J

+ e−vi,t θi − βi θi dt = 0,
wobei wir evi,t aus der Gleichung (3.9) erhalten haben. Nun wenden wir auf den ersten Teil
das Lemma von Ito an.
dvi,t 1
(dvi,t )2
dC
1
(dC)2
Et (1 − γi )(
− (1 − γi )γi
) + θi
+ θi (θi − 1) 2
C
2
C
vi,t
2
vi,t
dCt dvi,t
+ (θi )2 (1 − γi )
+ θi e−vi,t dt − βi θi dt = 0
C vi,t
20
dC 2
Jetzt können wir noch die Definitionen für dvi,t (3.5), (dvi,t )2 , dC
C (3.2) und ( C ) einsetzen
und erhalten
1
⇔ Et (1 − γi )(µCi dt + σCi dWt ) − (1 − γi )γi (µCi dt + σCi dWt )2
2
1
+ θi (µvi dt + σvi dWt ) + θi (θi )(µvi dt
2
2
+ σvi dWt ) + (1 − γi )θi (µCi dt + σCi dWt )(µvi dt + σvi dWt ) + θi e
−vi,t
dt − βi θi dt = 0
1
2
⇔ Et (1 − γi )(µCi dt + σCi dWt ) − (1 − γi )γi (µ2Ci (dt)2 + 2µCi σCi dtdWt + σC
(dWt )2 )
i
2
1
+ θi (µvi dt + σvi dWt ) + θi (θi )(µ2vi (dt)2 + 2µvi σvi dtdWt
2
+ σv2i (dWt )2 ) + (1 − γi )θi (µCi µvi (dt)2 + µCi σvi dtdWt + µvi σCi dtdWt
2
−vi,t
+ µvi σvi (dWt ) ) + θi e
dt − βi θi dt = 0
Jetzt nutzen wir noch aus, dass (dWt )2 = dt und wir sowohl (dt)2 als auch dtdWt vernachlässigen können. Da bei einem stochastischen Prozess der Erwartungswert gleich dem
Drift, also in unserem Fall Et [xdt + ydWt ] = x ist, können wir nun den Erwartungswert
ausrechnen. Als letzte Veränderung dividieren wir dann die gesamte Gleichung durch θi . So
erhalten wir unsere gewünschten Differentialgleichungen für i=1,2.
e−vi − βi + (1 −
1
0
0
)(µCi − 0, 5γi σC
σ ) + µvi + 0, 5θi σv0 i σvi + (1 − γi )σC
σ =0
i Ci
i vi
ψi
(3.10)
3.4 Pricing Kernel
Der Pricing Kernel gibt den Wert einer Einheit im jeweiligen volkswirtschaftlichen Zustand
wieder.
ξi,t = βiθi e−βi θi t−(1−θi )
Rt
0
e−vi,s ds (θi −1)vi,t
e
−γi
Ci,t
Wir betrachten in dieser Arbeit nur den Fall γ, ψ 6= 1. Für die anderen Fälle verweisen wir
auf [BCDG10]. Mit Hilfe des Lemmas von Ito folgt:
3
Modellherleitung
dξi,t
ξi,t
21
= −[βi θi + (1 − θi )e−vi + (1 − θi )µvi + γi µCi
0 σ
0 σ ]dt
−0.5(1 − θi )2 σv0 i σvi − 0.5γi (1 + γi )σC
− γi (1 − θi )σC
i Ci
i vi
−[γi σCi + (1 − θi )σvi ]dWt
Aus dieser Gleichung können wir die risikofreie Zinsrate und den Marktpreisrisikosatz erlangen. Der risikofreie Zinssatz ist identisch mit dem negativen Drift und die negative Volatilität
liefert den Marktpreisrisikosatz.
Bevor wir diese beiden Sätze bestimmen, vereinfachen wir die Gleichung. Dafür setzen wir
(3.10), nach e−v aufgelöst, in die Gleichung ein:
dξi,t
1
0
σ )
= − βi θi + (1 − θi ) βi − (1 − )(µCi − 0, 5γi σC
i Ci
ξi,t
ψi
0
− µvi − 0, 5θi σv0 i σvi − (1 − γi )σC
σ + (1 − θi )µvi + γi µCi
i vi
0
0
− 0.5(1 − θi )2 σv0 i σvi − 0.5γi (1 + γi )σC
−
γ
(1
−
θ
)σ
σ
σ
i
i Ci vi dt
i Ci
− [γi σCi + (1 − θi )σvi ]dWt
⇔
dξi,t
1
1 0
0
= − [βi + µCi − 0.5γi (1 + )σC
σ − 0.5(1 − θi )σv0 i σvi − (1 − θi )σC
σ ]dt
i Ci
i vi
ξi,t
ψi
ψi
(3.11)
− [γi σCi + (1 − θi )σvi ]dWt
3.5 Drift und Volatilität
Wie im vorigen Kapitel bereits erwähnt, bildet der negative Drift von (3.11) den risikofreien
Zinssatz für
1 0
1
0
µC1 − 0.5γ1 (1 +
)σ σC − 0.5(1 − θ1 )σv0 1 σv1 − (1 − θ1 )σC
σ
1 v1
ψ1
ψ1 C1 1
1
1 0
0
bzw.r2 = β2 +
µC2 − 0.5γ2 (1 +
)σ σC − 0.5(1 − θ2 )σv0 2 σv2 − (1 − θ2 )σC
σ .
2 v2
ψ2
ψ2 C2 2
r1 = β1 +
22
Da wir in unserem Modell einen vollkommenen Markt betrachten, gilt r1 = r2 . Also folgt
1 0
1
0
σ
µC1 − 0.5γ1 (1 +
)σ σC − 0.5(1 − θ1 )σv0 1 σv2 − (1 − θ1 )σC
1 v1
ψ1
ψ1 C1 1
1
1 0
0
= β2 +
σ .
µC2 − 0.5γ2 (1 +
)σ σC − 0.5(1 − θ2 )σv0 2 σv2 − (1 − θ2 )σC
2 v2
ψ2
ψ 2 C2 2
β1 +
Wir nutzen (3.3) und µC2 = µC + Xt −
1
1−ω µω
−
1
1−ω σω σC
und formen die Gleichung nach
dem gesuchten Drift µω um.
µω = −σω0 σC +
ω(1 − ω)ψ1 ψ2
1
1
−
)µc
{β2 − β1 + (
ψ
ψ
ψ
2
1
1 0
0
− 0, 5γ2 (1 +
σ
)σ σC − 0, 5(1 − θ2 )σv0 2 σv2 − (1 − θ2 )σC
2 v2
ψ2 C2 2
1 0
0
σ }
+ 0, 5γ1 (1 +
)σ σC + 0, 5(1 − θ1 )σv0 1 σv1 + (1 − θ1 )σC
1 v1
ψ1 C1 1
Dieser gilt für ω, also für Investor 1. Für Investor 2 gilt µω2 = −µω .
Bei der Volatilität σω gehen wir recht ähnlich vor. Der gesuchte Marktpreisrisikosatz entspricht der negativen Volatilität von Gleichung (3.11). Da wir einen vollkommenen Markt
haben, gilt hier mp1 = mp2 und damit
γ1 (σC +
σv1 =
dv1
dω σω
+
dv1
dx σx
1
1
σω ) + (1 − θ1 )σv1 = γ2 (σC − σω ) + (1 − θ2 )σv2 .
ω
ω
2
und σv2 = − dv
dω σω +
dv2
dx σx
in die Gleichung eingesetzt und nach σω
aufgelöst, ergibt:
σω =
∂v1
2
(γ2 − γ1 )σC + (1 − θ2 ) ∂v
∂x − (1 − θ1 ) ∂x )σx
Die Gleichung erweitert mit
sion γ =
ω γ1
1
1
+(1−ω) γ1
γ1
ω
+
ω(1−ω)
γ1 ·γ2
γ2
1−ω
∂v2
1
+ [(1 − θ1 ) ∂v
∂ω − (1 − θ2 ) ∂ω ]
.
und mit der Definition der durchschnittlichen Risikoaver-
, ergibt
2
σω =
ω(1−ω)
∂v2
∂v1
1
γ2 )σC + γ1 γ2 [(1 − θ2 ) ∂x − (1 − θ1 ) ∂x ]σx
ω(1−ω)
∂v1
∂v2
1
γ + γ1 γ2 [(1 − θ1 ) ∂ω − (1 − θ2 ) ∂ω ]
ω(1 − ω)( γ11 −
für den 1. Investor. Für den 2. Investor gilt ähnlich wie beim Drift σω2 = −σω .
(3.12)
3
Modellherleitung
23
3.6 Übersicht aller relevanten Gleichungen
Zum Abschluss dieses Kapitels werden wir noch eine kleine Übersicht aller relevanten Gleichungen zusammenstellen. So können wir uns einen besseren Überblick über die Differentialgleichungen und ihre Variablen machen. Die partiellen Differentialgleichungen für i = 1, 2
lauten
0 = e−vi − βi + (1 −
1
0
0
)(µCi − 0, 5γi σC
σ ) + µvi + 0, 5θi σv0 i σvi + (1 − γi )σC
σ
i Ci
i vi
ψi
Drift
σω =
ω(1−ω)
∂v2
∂v1
1
γ2 )σC + γ1 γ2 [(1 − θ2 ) ∂x − (1 − θ1 ) ∂x ]σx
ω(1−ω)
∂v1
∂v2
1
γ + γ1 γ2 [(1 − θ1 ) ∂ω − (1 − θ2 ) ∂ω ]
ω(1 − ω)( γ11 −
σω2 = −σω
Volatilität
µω = −σω0 σC +
1
ω(1 − ω)ψ1 ψ2
1
−
)µc
{β2 − β1 + (
ψ2 ψ1
ψ
1 0
0
− 0, 5γ2 (1 +
)σ σC − 0, 5(1 − θ2 )σv0 2 σv2 − (1 − θ2 )σC
σ
2 v2
ψ2 C2 2
1 0
0
)σ σC + 0, 5(1 − θ1 )σv0 1 σv1 + (1 − θ1 )σC
σ }
+ 0, 5γ1 (1 +
1 v1
ψ1 C1 1
µω2 = −µω
Volatilität des Konsums
µC1 = µC + Xt +
µC2 = µC + Xt −
1
1
µω + σω σC
ω
ω
1
1
µω −
σω σC
1−ω
1−ω
24
Drift des Konsums
σC1 = σC +
σC2 = σC −
1
σω
ω
1
σω
1−ω
Drift der Funktion
σv1 =
∂v1
∂v1
σω +
σx
∂ω
∂x
σv2 = −
∂v2
∂v2
σω +
σx
∂ω
∂x
Volatilität der Funktion
µv1 =
∂v1
∂v1
∂ 2 v1
∂ 2 v1
∂ 2 v1 0
µω −
κXt + 0, 5 2 σω0 σω + 0, 5 2 σx0 σx +
σ σx
∂ω
∂x
∂ω
∂x
∂ω∂x ω
µv2 = −
∂vi
∂ 2 vi 0
∂ 2 v2 0
∂ 2 v2 0
∂vi
µω −
κXt + 0, 5
σ
σ
+
0,
5
σ
σ
−
σ σx
ω
x
∂ω
∂x
∂(ω)2 ω
∂x2 x
∂ω∂x ω
Folgende Hilfsdefinitionen haben wir genutzt:
γ=
θ1 =
1 − γ1
1 − ψ11
θ2 =
1 − γ2
1 − ψ12
ω γ11
1
+ (1 − ω) γ12
ψ = ωψ1 + (1 − ω)ψ2
25
4 Wohlstandskonsumquotient als gewöhnliche
Differentialgleichung
Wir beginnen dieses Kapitel damit beginnen, für unsere Differentialgleichungen ein paar
Einschränkungen zu definieren. Unter diesen Einschränkungen handelt es sich dann um
gewöhnliche Differentialgleichungen. Die Differerentialgleichungen werden wir in 4.1 analysieren. Danach folgt dann die Randwertbetrachtung für ω. Zum Abschluss des Kapitels werden
wir die numerische Berechnung unserer gewöhnlichen Differentialgleichungen mit Hilfe des
Newton-Verfahrens verwirklichen.
4.1 Gewöhnliche Differentialgleichung
Wie im letzten Kapitel hergeleitet, lauten die zu untersuchenden Differentialgleichungen
(3.10)
0 = e−vi (ω) − βi + (1 −
1
0
0
)(µCi − 0, 5γi σC
σ ) + µvi + 0, 5θi σv0 i σvi + (1 − γi )σC
σ
i Ci
i vi
ψi
(4.1)
mit den Definitionen für µCi (3.3), σCi (3.4), µvi (3.6) und σvi (3.7). In diesem Kapitel definieren einige Einschränkungen für unsere Differentialgleichungen. Wir nehmen an, dass sich
sowohl der Zeitpräferenzfaktor der beiden Investoren β ≡ β1 = β2 , als auch ihr Risikoaversionsfaktor γ ≡ γ1 = γ2 , nicht unterscheidet. Außerdem analysieren wir in diesem Kapitel
die Differentialgleichungen unter der Annahme, dass die Long-Run-Risk-Variable X zu vernachlässigen ist. Dadurch vereinfachen sich die Terme einiger Definitionen und wir können
unsere Differentialgleichungen für i = 1, 2 darstellen durch
0 = e−v1 (ω) − β +
ψ1 − 1
ω(1 − ω)(ψ1 − ψ2 ) ∂v1
2
2
· [µC − 0.5γσC
]+
·
[µC − 0.5γσC
]
ωψ1 + (1 − ω)ψ2
ωψ1 + (1 − ω)ψ2
∂ω
(4.2)
0 = e−v2 (ω) − β +
ψ2 − 1
ω(1 − ω)(ψ2 − ψ1 ) ∂v2
2
2
[µC − 0.5γσC
].
· [µC − 0.5γσC
]−
·
ωψ1 + (1 − ω)ψ2
ωψ1 + (1 − ω)ψ2
∂ω
(4.3)
26
Beide Gleichungen hängen damit nur noch von der Funktion vi (w), ihrer Ableitung
∂vi
∂ω
und
den festen Werten β, ψ1 , ψ2 , µC , γ und σC ab.
Zur Analyse schreiben wir beide Differentialgleichungen in der Form:
0 = e−vi + a(ω)i ·
∂vi
+ c(ω)i
∂ω
für i=1,2
(4.4)
mit
ω(1 − ω)(ψ1 − ψ2 )
2
· [µC − 0.5γσC
]
ωψ1 + (1 − ω)ψ2
ψ1 − 1
2
· [µC − 0.5γσC
]
c(ω)1 = −β +
ωψ1 + (1 − ω)ψ2
ψ2 − 1
2
c(ω)2 = −β +
· [µC − 0.5γσC
].
ωψ1 + (1 − ω)ψ2
a(ω) ≡ a(ω)1 = a(ω)2 =
Um zu entscheiden, ob wir zur Approximation der Ableitung
∂vi
∂ω
, den Vorwärts- oder
Rückwärtsdifferenzenquotienten [siehe dazu [Bur07, S.11]] benutzen, untersuchen wir das
Vorzeichen von a(ω). Dazu betrachten wir ein Beispiel.
Beispiel 4.1.1
β = 0, 02
ψ1 = 0, 5
ψ2 = 1, 2
µC = 0, 02
γ=6
σC = 0, 0252
Für dieses Beispiel mit ω = 0, 5 ergibt sich für a(ω) der Wert =-0,0037. Falls wir die beiden
Werte für ψi vertauschen, erhalten wir a(ω) = +0, 0037. Das heißt, es können sowohl positive
als auch negative Werte für a(ω) vorkommen. Deshalb ist es nicht möglich, nur Rückwärtsoder Vorwärtsdifferenzenquotienten zu betrachten. Daher sollten wir zur Approximierung
der vorkommenden Ableitung das Upwind-Verfahren [Bur07, S.13] nutzen. Die Voraussetzung ist jedoch, dass es für feste Paramater für ω von 0 bis 1 keinen Vorzeichenwechsel
von Minus nach Plus von a(ω) gibt. Da allerdings
(ψ1 − ψ2 )[µC −
erwarten.
2]
0.5γσC
ω(1−ω)
ωψ1 +(1−ω)ψ2
> 0 für alle ψ1 , ψ2 6= 0 und
unabhängig von der Variable ω ist, sind keine Schwierigkeiten zu
4
Wohlstandskonsumquotient als gewöhnliche Differentialgleichung
Wir approximieren die Ableitung
∂vi
∂w
27
durch das Upwind-Verfahren, wobei wir jetzt aus
Übersichtsgründen auf den Index i verzichten. Der Index j hängt von der Variablen ω ab, da
wir verschiedene Werte für ω betrachten werden. h ist die Schrittweite für die Ableitungsapproximation. Es gilt
v(ωj ) − v(ωj−1 )
v(ωj+1 ) − v(ωj )
∂v(ωj )
= −max(a(ωj ), 0) ·
− min(a(ωj ), 0) ·
,
∂ω
h
h
(4.5)
wobei wir die beiden Fälle (a(ωj ) positiv oder negativ) trennen können:
∂v(ωj )
v(ωj ) − v(ωj−1 )
=−
∂ω
h
∂v(ωj )
v(ωj+1 ) − v(ωj )
=−
∂ω
h
für a(ωj ) ≥ 0
(4.6)
für a(ωj ) < 0.
(4.7)
Diese Definitionen eingesetzt in (4.4) und etwas vereinfacht liefert uns
a(ωj )
a(ωj )
v(ωj ) +
v(ωj−1 ) + c(ωj )
h
h
a(ωj )
a(ωj )
v(ωj ) −
v(ωj+1 ) + c(ωj )
0 = e−v(ωj ) +
h
h
0 = e−v(ωj ) −
für a(ωj ) ≥ 0
(4.8)
für a(ωj ) < 0,
(4.9)
wobei wir den rechten Teil der Gleichungen als Funktion f + (v(ωj )) bzw. f − (v(ωj )) bezeichnen
können. Das heißt, wir haben jeweils eine Gleichung der Form 0 = f (v(ωj )) zu lösen. Um
diese zu berechnen, nutzen wir das in (2.2.15) vorgestellte Newton-Verfahren.
v(ωn+1 ) = v(ωn ) −
f (v(ωn ))
f 0 (v(ωn ))
(4.10)
Wenn wir nun die beiden Gleichungen f + (v(ωj )) und f − (v(ωj )) in das Newton-Verfahren
einsetzen und weiter zusammenfassen, erhalten wir:
v(ωn+1 ) =
v(ωn )he−v(ωn ) + he−v(ωn ) + a(ωn )v(ωj−1 ) + hc(ωn )
a(ωn ) + he−v(ωn )
für a(ωn ) > 0 ,
(4.11)
v(ωn+1 ) =
v(ωn )he−v(ωn ) + he−v(ωn ) − a(ωn )v(ωj+1 ) + hc(ωn )
−a(ωn ) + he−v(ωn )
für a(ωn ) < 0.
(4.12)
Diese Gleichungen können wir iterativ lösen und erhalten dann v(ωj ). Dafür brauchen wir,
je nach Gleichung, v(ωj−1 ) oder v(ωj+1 ). Da die Randwerte für ω = 0 und ω = 1 berechenbar sind (siehe weiter unten), lassen sich auch mit diesen v(wj ) für alle j berechnen. Dazu
wird je nach Vorzeichen von a(ωj ) der rechte (Rückwärtsdifferenzenquotienten) oder linke
(Vorwärtsdifferenzenquotienten) Randwert der Differentialgleichung verwendet.
28
4.2 Randwerte
In diesem Abschnitt werden wir die Randwerte unserer Differentialgleichungen berechnen.
ω liegt im Bereich von 0 bis 1. Für die Randwerte müssen wir also die Differentialgleichungen
(4.2) und (4.3) für ω = 0 und ω = 1 lösen. Wir betrachten die Ränder für die Differentialgleichung von Investor 1.
ω=0:
ψ1 − 1
2
0 = e−v1 (0) − β +
· [µC − 0.5γσC
]
ψ2
ψ2
⇒ v1 (0) = ln
2]
ψ2 β − (ψ1 − 1)[µC − 0.5γσC
ω=1:
ψ1 − 1
2
· [µC − 0.5γσC
]
0 = e−v1 (1) − β +
ψ1
ψ1
⇒ v1 (1) = ln
2]
ψ1 β − (ψ1 − 1)[µC − 0.5γσC
Die Randwerte für Investor 2 ergeben sich analog, wobei darauf geachtet werden muss, dass
ω2 = 1 − ω gilt.
4
Wohlstandskonsumquotient als gewöhnliche Differentialgleichung
29
4.3 Numerische Berechnung
In diesem Abschnitt werden wir die Realisierung des im letzten Kapitel überlegten Programms
zur Lösung unserer Differentialgleichungen bewirken. Wir nutzen dafür Matlab R2012b.
Wie bereits gezeigt, lässt sich zur numerischen Berechnung der Differentialgleichungen das
Upwind-Verfahren für die Darstellung der Ableitung und das Newton-Verfahren zur Berechnung der Funktionswerte nutzen. In diesem Abschnitt wird gezeigt, wie man anhand der
zugehörigen Matlab-Funktionen ein konkretes Beispiel berechnen kann. Dazu nutzen wir das
Beispiel 4.1.1 mit den Werten β = 0, 02, ψ1 = 0, 5, ψ2 = 1, 2, µC = 0, 02, γ = 6, σC = 0, 0252.
Als erstes übergeben wir dafür an unser Hauptprogramm UpwindNewton die benötigten
konstanten Werte und die Anzahl der Schritte n zur Approximation der Ableitungen. Das
Programm liefert uns dann die Werte für unsere beiden Differentialgleichungen v1 und v2
zurück.
1 [v1,v2]=UpwindNewton(beta,psi 1,psi 2,mu C,gamma,sigma C,n);
Beispiel [v1,v2] = UpwindNewton(0.02,0.5,1.2,0.02,6,0.0252,10000);
Das Hauptprogramm ist wie folgt programmiert, wobei wir uns auf den Fall v1 beschränken,
da v2 analog zu berechnen ist.
1 function [v1,v2]= Name(beta,psi1,psi2,muC,gamma,sigmaC,n,g)
2
3 % Wir erzeugern den Platz fuer unsere Variablen
4 v1=zeros(n+1,1);
6 aw1=zeros(n+1,1);
8 cw1=zeros(n+1,1);
10
11
12 % h ist unsere Schrittweite fuer Omega
13 h=1/n;
14
15 % Omega ist im Bereich von 0 bis 1 mit der Schrittweise h definiert.
16 w=0:h:1;
17
18 % Als naechstes berechnen wir unsere Randwerte fuer v1 (und v2) w=[0,1]
19 % und uebergeben sie an v1 (und v2) an die erste bzw. letzte Stelle
20 v1 0 = log(psi2/(psi2*beta−(psi1−1)*(muC−0.5*gamma*(sigmaC)ˆ2))); %w 1=0
21 v1(1)=v1 0;
22
23 v1 1 = log(psi1/(psi1*beta−(psi1−1)*(muC−0.5*gamma*(sigmaC)ˆ2))); %w 1=1
24 v1(n+1)=v1 1;
30
In unserem Beispiel wären die Randbedingungen für Investor 1: 3,5921 für ω = 0 und 3,2677
für ω = 1. Als nächstes berechnet das Programm die Werte für a(w) und c(w).
37
38 for i=2:n
39 aw1(i)=(w(i)*(1−w(i))*(psi1−psi2))/(w(i)*psi1+(1−w(i))*psi2)*...
(muC−0.5*gamma*sigmaCˆ2);
41
42 cw1(i)=−beta+(psi1−1)/(w(i)*psi1+(1−w(i))*psi2)*(muC−0.5*gamma*sigmaCˆ2);
44 end;
45
46
47
48 for i=2:n
49
% Wir beginnen mit der Berechnung von v1, dafuer machen wir eine
50
% Fallunterscheidung, ob aw1(i) positiv oder negativ ist. Dadurch
51
% entscheiden wir, ob das Programm den Vorwaerts− oder
52
% Rueckwaertsdifferenzenquotienten nutzen soll.
53
if aw1(i) > 0
54
% Nun starten wir das Newton Verfahren mit Startwert 1, wobei wir
55
% damit unsere while−Schleife richtig beginnt, 2 Werte vorgeben muessen.
56
y1(1)=0;
57
y1(2)=1; % Startwert
58
j=2;
59
while abs(y1(j)−y1(j−1)) >0.000001
60
% Das Newtonverfahren wird solange ausgefuehrt, bis der Unterschied
61
% zwischen den beiden letzten Schritten kleiner als 0.000001 ist
62
y1(j+1)= (aw1(i)*v1(i−1)+h*exp(−y1(j))*y1(j)+h*exp(−y1(j))+ ...
h*cw1(i))/(aw1(i)+h*exp(−y1(j)));
63
j=j+1;
64
end;
65
v1(i)=y1(j);
66
% Uebergabe vom errechneten Wert an v1.
67
else
68
% Falls aw1(i) negativ ist, muessen wir "rueckwaerts" rechnen, das heisst
69
% wir berechnen unter Nutzung des rechten Randes den Wert von v1(n)
70
% und dann v1(n−1) usw. Dafuer definieren wir uns als Zaehler k, der
71
% von n bis 2 geht.
72
k=n+2−i;
73
% Newton Verfahren mit Startwert 1
74
y1(1)=0;
75
y1(2)=1;
76
j=2;
77
while abs(y1(j)−y1(j−1)) >0.000001
78
y1(j+1)= (−aw1(k)*v1(k+1)+y1(j)*h*exp(−y1(j))+h*exp(−y1(j))+...
h*cw1(k))/(−aw1(k)+h*exp(−y1(j)));
79
j=j+1;
80
end;
4
Wohlstandskonsumquotient als gewöhnliche Differentialgleichung
31
Abbildung 4.1: Der grüne Graph stellt die Funktion v1 und der rote Graph die Funktion v2
dar.
81
v1(k)=y1(j);
82
end;
83 end
84 % Nun haben wir v1 komplett berechnet und lassen den Funktionsgraphen zeichnen.
85 plot(w,v1,'r')
115 end
Das Programm erfüllt dasselbe nochmal für v2 und liefert die Graphen der beiden Funktionen
v1 und v2 .
Abbildung 4.1 zeigt die für unser Beispiel 4.1.1 ausgerechneten Funktionswert für v1 und v2 .
Was auffällt ist, dass die Randwerte 3,5921 für ω = 0 und 3,2677 für ω = 1 von v1 den
maximale und den minimale Wert der Funktion darstellen.
32
Abbildung 4.2: Die maximale Newton-Schrittanzahl unter den Parametern Genauigkeit und
Schrittanzahl für ω.
Das Newtonverfahren liefert sehr schnell eine gute Annäherung. Es braucht für die Berechnung sehr wenige Schritte, wie man in der Abbildung 4.2 deutlich sehen kann.
Man sieht, dass das Newton-Verfahren selbst bei hoher Genauigkeit in je max. 10 Schritten
zu einer guten Lösung kommt.
4
Wohlstandskonsumquotient als gewöhnliche Differentialgleichung
33
Abbildung 4.3: Der maximale Fehler für v1 roter Graph und für v2 grüner Graph von unserem Beispiel in Abhängigkeit von der Schrittanzahl [2:100] bei einer NewtonGenauigkeit von 10−6 .
Fehler
Wir werden für unser Beispiel noch die Entwicklung des Fehlers betrachten. Dafür berechnen
wir mit Hilfe des zentralen Differenzenquotienten die Ableitungen an den berechneten Stellen
von v1. Wir nehmen an, dass für steigende Schrittanzahl n der maximale Fehler kleiner wird.
Wie man in der Abbildung 4.3 sehen kann, sinkt der Fehler bei einer Schrittanzahl von 2 bis
100 deutlich.
Der maximale Fehler von v1 liegt bei n=2 im Bereich von 1, 0 · 10−4 , bei n=100 im Bereich
von 2, 8 · 10−6 . Wenn wir n noch weiter steigern, fällt der Fehler wie erwartet weiter, bei
n=10000 in einem Bereich von 2, 8 · 10−8 und bei n = 106 liegt er sogar im Bereich von
2, 8 · 10−10 . Der Fehler hängt allerdings nicht nur von der Schrittweite, sondern auch von der
verlangten Genauigkeit des Newtonverfahrens ab. In der Abbildung 4.4 erkennt man deutlich
die Abhängigkeit des maximalen Fehlers von der Schrittweite und der Newtongenauigkeit.
Die gerundeten Werte der Abbildung 4.4 können in der Tabelle 4.1 nachgelesen werden. Wobei
NG die Newtongenauigkeit bezeichnet.
Wie erwartet passiert der kleinste maximale Fehler bei der Schrittanzahl von n=1.000.000
und der Newtongenauigkeit von 10−6 . Allerdings fällt auch auf, dass die Newtongenauigkeit
deutlich größeren Einfluss auf die Fehlerminimierung hat als die Schrittanzahl. Allerdings
kann der Fehler nur bei der Veränderung beider Faktoren wirklich minimiert werden.
34
Abbildung 4.4: Der maximale Fehler für v1 unseres Beispiels in Abhängigkeit von der Newtongenauigkeit und der Schrittanzahl.
101
n=
n = 102
n = 103
n = 104
n = 105
n = 106
N G < 10−1
1, 2 · 10−4
1, 6 · 10−4
1, 4 · 10−4
2, 0 · 10−4
2, 0 · 10−4
2, 0 · 10−4
N G < 10−2
2, 7 · 10−5
2, 8 · 10−6
1, 6 · 10−6
1, 7 · 10−6
1, 9 · 10−6
1, 9 · 10−6
N G < 10−3
2, 7 · 10−5
2, 8 · 10−6
2, 8 · 10−7
2, 8 · 10−8
1, 8 · 10−8
1, 9 · 10−8
N G < 10−4
2, 7 · 10−5
2, 8 · 10−6
2, 8 · 10−7
2, 8 · 10−8
2, 8 · 10−9
2, 8 · 10−10
N G < 10−5
2, 7 · 10−5
2, 8 · 10−6
2, 8 · 10−7
2, 8 · 10−8
2, 8 · 10−9
2, 8 · 10−10
N G < 10−6
2, 7 · 10−5
2, 8 · 10−6
2, 8 · 10−7
2, 8 · 10−8
2, 8 · 10−9
2, 8 · 10−10
Tabelle 4.1: Der maximale Fehler von v1 (von unserem Beispiel) in Abhängigkeit der Newtongenauigkeit und der Schrittanzahl.
4
Wohlstandskonsumquotient als gewöhnliche Differentialgleichung
35
Probleme
Leider kann unsere Vorgehensweise für bestimmte Parameter ein paar Probleme bereiten.
Falls ψ1 und ψ2 zu weit auseinanderliegen (ca. |ψ1 − ψ2 | > 1), kann es zu Problemen bei
der Berechnung der Randwerte kommen. Eine oder sogar beide Randwerte können irregulär
werden, da wir den natürlichen Logarithmus einer negativen Zahlen berechnen wollen. Diese
Irregularität überträgt sich dann auch auf die anderen Werte der Funktion. In Kapitel 6
werden wir nochmals kurz darauf eingehen und eine Möglichkeit zur Verhinderung dieses
Problems vorstellen. Für die meisten und vor allem realistischen Problemstellungen liefert
der Algorithmus allerdings schnell eine gute Lösung für die Differentialgleichungen.
36
5 Wohlstandskonsumquotient als partielle
Differentialgleichung
In Kapitel 4 haben wir die Differentialgleichungen unter verschiedenen Einschränkungen betrachtet und numerisch gelöst. In diesem Kapitel werden wir unsere Differentialgleichungen
nicht so stark einschränken. Wir nehmen zwar weiterhin an, dass der Zeitpräferenzfaktor und
der Risikoaversionsfaktor beider Investoren übereinstimmen (β = β1 = β2 und γ = γ1 = γ2 ),
allerdings werden wir in diesem Kapitel die Long-Run-Risk-Variable X nicht vernachlässigen.
Deswegen wird bei unseren Differentialgleichungen die partielle Eigenschaft weiterhin erhalten bleiben. Das heißt, dass außer den schon in Kapitel 4 auftretenden Ableitungen von v1
nach ω und v2 nach ω, auch die Ableitungen der beiden Funktionen nach x auftreten können.
Außerdem werden auch die Ableitungen 2. Ordnung sowohl nach ω als auch nach x in unseren
Differentialgleichungen vorkommen. Als letztes werden wir auch noch gemischte Ableitungen
nach ω und x betrachten.
Wir beginnen dieses Kapitel mit der allgemeinen Strukturbetrachtung unserer Differentialgleichungen. Als numerische Lösung unserer Differentialgleichungen werden wir den Weg über
ein Anfangswertproblem gehen. Dafür werden wir dann noch die Randwerte beider Variablen
und die Approximationen der Ableitungen betrachten. Wir werden unser Kapitel mit der
Implementierung unseres Lösungsverfahren vervollständigen.
Unsere Differentialgleichungen haben für i = 1, 2 die Form
0 = e−vi (ω) − βi + (1 −
1
0
0
σ ) + µvi + 0, 5θi σv0 i σvi + (1 − γi )σC
)(µCi − 0, 5γi σC
σ
i Ci
i vi
ψi
(5.1)
mit den entsprechenden Definitionen (siehe Kapitel 3.6).
5.1 Strukturbetrachtung
Als erstes möchten wir die allgemeine Struktur unserer Differentialgleichungen betrachten.
Anhand der Zusammenfassung aller relevanten Gleichungen in Kapitel 3.6 ist zu erkennen,
dass es sich bei unseren Differentialgleichungen um Gleichungen 2. Ordnung handelt. Wir unterscheiden 3 verschiedene Grundtypen von Differentialgleichungen: die elliptische, die parabolische und die hyperbolische Differentialgleichung [Bur07, Kapitel 1]. Zur Einteilung unserer
5
Wohlstandskonsumquotient als partielle Differentialgleichung
37
Differentialgleichung von Investor 1, wobei wir den Index weglassen (Investor 2 funktioniert
analog), schreiben wir unsere Gleichung mit Hilfe eines Differentialoperators L um.
Lv = a1 (ω, x)
∂2v
∂2v
∂2v
∂v
∂v
+
a
(ω,
x)
+
a
(ω,
x)
+ b1 (ω, x)
+ b2 (ω, x) +c(ω, x)v
2
3
2
2
∂ω
∂ω∂x
∂x
∂ω
∂x
(5.2)
= f (ω, x)
Wie man leicht sieht, sind die einzigen Ableitungen 2. Ordnung in unserer Differentialgleichung durch µv =
∂v
∂ω µω
−
∂v
∂x κXt
2
2
2
∂ v 0
∂ v 0
∂ v 0
+ 0, 5 ∂ω
2 σω σω + 0, 5 ∂x2 σx σx + ρωx ∂ω∂x σω σx definiert. Wir
können also
a1 (ω, x) = 0, 5σω0 σω
a2 (ω, x) = ρωx σω0 σx
a3 (ω, x) = 0, 5σx0 σx
definieren, wobei σx zwar ein fester Wert ist, aber σω (3.12) in Abhängigkeit von mehreren
Ableitungen nicht linear ist.
Zur Einteilung unserer Differentialgleichung in die 3 Grundtypen betrachten wir die folgende
Definition:
Definition 5.1.1 [Nat01]
Es sei eine Differentialgleichung in der !
Form (5.2) gegeben.
a2 (ω,x)
a1 (ω, x)
2
Falls für A(ω, x) = a2 (ω,x)
a
(ω,
x)
3
2
• A(ω, x) positiv oder negativ definit ist, so nennen wir die Differentialgleichung elliptisch,
• A(ω, x) positiv oder negativ semidefinit, aber nicht definit ist, so bezeichnen wir die Differentialgleichung parabolisch,
• A(ω, x) indefinit (mit genau einem negativen Eigenwert) ist, so nennen wir die Differentialgleichung hyperbolisch.
Zur Untersuchung der Matrix A(ω, x) auf Definitheit betrachten wir die Determinante
D(ω, x) ≡ a1 (ω, x)a3 (ω, x) −
a2 (ω, x) 2
.
2
38
In unserem Fall lautet die Matrix A(ω, x)
A(ω, x) =
0, 5σω0 σω
0 σ
ρωx σω
x
2
0 σ
ρωx σω
x
2
0, 5σx0 σx
!
und es gilt
0, 5σω0 σx 2
2
= 0, 25σω0 σω σx0 σx − 0.25ρωx σω0 σω σx σx
D(ω, x) = 0, 5σω0 σω · 0, 5σx σx −
= (0, 25σω0 σω σx0 σx )(1 − ρωx ),
da für σω0 σω ≥ 0, σx0 σx ≥ 0 und ρωx < 1 (siehe (3.6)) folgt
D(ω, x) > 0
∀(ω, x).
Also handelt es sich bei unserer Gleichung um eine nichtlineare elliptische Differentialgleichung.
Es gibt 2 Arten von Maximumsprinzipien, die wir jetzt näher vorstellen wollen. Dabei orientieren wir uns an [Bur07, S.20ff], da sind auch die Beweise nachzulesen.
Satz 5.1.2 (schwaches Maximumprinzip)
Sei Lv ≤ 0 (Lv ≥ 0) mit L Differentialoperator wie oben. Dann gilt v ≤ 0 (v ≥ 0) oder v
nimmt sein globales Maximum (Minimum) am Rand von Ω an.
Satz 5.1.3 (starkes Maximumprinzip)
Sei Lv < 0 (Lv > 0) mit L wie oben. Dann gilt v ≤ 0 (v ≥ 0) oder u hat kein lokales
Maximum (Minimum) im Inneren von Ω.
Wie wir oben gesehen haben, gilt Lv1 = f (ω, x), wobei f (ω, x) der negative Teil der Differentialgleichung, der unabhängig von der Funktion v1 und ihren Ableitungen, ist. In unserem
Fall gilt also:
1 ) µC + Xt
ψ1
1
(1 − ω)ψ1 ψ2
1
1
1
2
2
2
+
(
−
)µC − 0, 5(1 +
)γσC + 0, 5(1 +
)γσC − 0, 5γσC
ψ2 ψ1
ψ2
ψ1
ψ
f (ω, x) =β − (1 −
5
Wohlstandskonsumquotient als partielle Differentialgleichung
39
Diese Gleichung können wir vereinfachen zu:
f (ω, x) = −β + (1 −
ψ1 − 1
1
2
)Xt +
[µC − 0, 5γσC
]
ψ1
ψ
Wir wissen, dass die Funktion v1 stets nicht negativ ist. Falls für alle ω ∈ [0, 1] und Xt
gilt, dass Lv ≤ 0 folgt daraus, dass unser Differentialoperator das schwache Maximumprinzip
erfüllt. Das heißt v1 nimmst sein globales Maximum am Rand an. Später werden wir für
unser Beispiel sehen, dass v1 genau das erfüllt. Neben dieser wichtigen Eigenschaft können
wir mit Hilfe des Maximumprinzip auch die Eindeutigkeit der Lösung beweisen:
Da für 2 Lösungen v1 , v2 die Differenz v = v1 − v2 die Gleichung Lv = 0 erfüllt, sowie v = 0
am Rand gilt, folgt durch 0 ≤ v ≤ 0 in Ω die Eindeutigkeit der Lösung [Bur07, S.22].
Wie wir gesehen haben, tritt σω in Verbindung mit verschiedenen Termen auf. Die Problematik die sich stellt ist, dass σω sowohl von den Ableitungen des 1. Investors nach x und ω
als auch von den Ableitungen des 2. Investors nach den beiden Argumenten abhängt. Das
heißt, wir können unsere beiden Differentialgleichungen nicht getrennt voneinander lösen man spricht deshalb von gekoppelten partiellen Differentialgleichungen. Deswegen werden wir
den Weg über ein Anfangswertproblem mithilfe des expliziten Eulerverfahrens wählen, wobei
wir erst mit der Diskretisierung beginnen.
5.2 Diskretisierung
In diesem Abschnitt werden wir als erstes das diskrete Gitter, auf dem wir unsere Differentialgleichung lösen möchten einführen. Danach werden wir ein numerisches Verfahren entwickeln,
mit dem wir unsere in ein Anfangswertproblem übertragende Differentialgleichung mithilfe
des Eulerverfahrens lösen können. Dieses Verfahren wird als sogenanntes Vorwärts-EulerVerfahren [Bur07, S.9f] bezeichnet.
Unsere Differentialgleichungen hängen von den Argumenten ω und x ab. Wir werden dafür
ein geeignetes Gitter einführen. Für ω werden wir die natürlichen Grenzen 0 und 1 als Beschränkung des Gitters betrachten. Für x sind zwar die Ränder, wie wir weiter unten sehen
werden, bei −∞ und +∞. Aber da diese Werte nicht realistisch für unser Modell sind, werden
wir unser Gitter auf xmin und xmax , die im Bereich von ±0, 025 liegen, definieren.
G = {(w, x) ∈ U |(w, x) = (wi , xj ), 1 ≤ i ≤ nw, 1 ≤ j ≤ nx},
wobei für die Schrittweite der Diskretisierung nach ω gilt: h1 =
x gilt: h2 =
xmax −xmin
.
nx−1
wmax −wmin
nw−1
und analog für
Wie wir oben bereits erwähnt haben, setzen wir wmin = 0 und
wmax = 1, damit gilt h1 =
1
nw−1 .
Je größer nw bzw. nx (h1 und h2 kleiner) ist, desto genauer
40
ist unser Diskretisierungsgitter. Da wir mit gleichbleibenden h1 und h2 rechnen, handelt es
sich bei unserem Gitter um ein reguläres. Das heißt, wir können unser Gitter als
G = {(w, x) ∈ U |wi = (i − 1) · h1 , xj = xmin + (j − 1) · h2 },
mit 1 ≤ i ≤ nw und 1 ≤ j ≤ nx
definieren. Auf diesem Gitter möchten wir unsere Differentialgleichung lösen. Da wir, wie
oben gesehen, die beiden Differentialgleichungen für Investor 1 und 2 nicht getrennt voneinander betrachten können, möchten wir die Gleichungen mit Hilfe eines Anfangswertproblems
lösen.
Definition 5.2.1 (Anfangswertproblem) [Ohl08, Kapitel 1]
Unter einem Anfangswertproblem verstehen wir ein Gleichungsystem für das gilt
∂y(t)
= f (t, y(t))
∂t
mit
y(t0 ) = y0 ,
wobei die letzte Bedingung unsere Anfangsbedingung darstellt. Das heißt, wir haben zu unserem Startzeitpunkt t0 den Anfangswert y0 .
Wir möchten das ganze mit Hilfe des expliziten Eulerverfahrens lösen. Das explizite Eulerverfahren, oder auch Eulersches Polygonzugverfahren, ist ein einfaches numerisches Verfahren
zur Lösung eines Anfangswertproblems.
Definition 5.2.2 (explizites Eulerverfahren) [Ohl08, Kapitel1]
Wir haben ein Anfangswertproblem gegeben. Nun betrachten wir für die Diskretisierungsschrittweite τ die diskreten Zeitpunkte tk = t0 + i · τ mit i=0,1,2,.... und berechnen die
Werte
yi+1 = yi + τ f (ti , yi )
für i=0,1,2...
Diese Werte stellen Approximationen an die tatsächlichen Werte der exakten Lösung des
Anfangswertproblems da.
5
Wohlstandskonsumquotient als partielle Differentialgleichung
41
Nun möchten wir diese beiden Definitionen verbinden und auf unser Problem anwenden.
Dafür schreiben wir für unsere Differentialgleichungen (5.1) für i=1,2:
1
0
0
σ
σ ) + µv1 + 0, 5θ1 σv0 1 σv1 + (1 − γ1 )σC
)(µC1 − 0, 5γ1 σC
1 v1
1 C1
ψ1
1
0
0
DGLv2 = e−v2 (ω) − β2 + (1 −
σ
)(µC2 − 0, 5γ2 σC
σ ) + µv2 + 0, 5θ2 σv0 2 σv2 + (1 − γ2 )σC
2 v2
2 C2
ψ2
DGLv1 = e−v1 (ω) − β1 + (1 −
Wir beginnen damit unsere DGLv1 bzw. DGLv2 als Anfangswertproblem zu betrachten.
∂v1
= DGLv1
∂t
∂v2
= DGLv2
∂t
Die Ableitung können wir mithilfe des Vorwärtsdifferenzenquotienten
dv1
dt
= − τ1 (v1i+1 − v1i )
approximieren und dann das Ganze nach v1i+1 auflösen
v1i+1 = v1i + τ · DGL(v1)
v2i+1 = v2i + τ · DGL(v2)
Damit haben wir also ein explizites Eulerverfahren. Für genügend große τ strebt v1i bzw.
v2i gegen die Lösung, unabhängig vom Startwert. Dies erfolgt, da DGLv1 und DGLv2 für
genügend große τ gegen 0 gehen und damit die Lösung zur einem stationären Punkt wird.
Durch das explizite Eulerverfahren haben wir, neben der räumlichen Komponente (ω, x),
auch eine zeitliche τ eingeführt. Für diese wählen wir die Diskretisierung τ = 0, ...., T mit
einer regulären Schrittweite.
5.3 Randwerte
Wir untersuchen nun die Randwerte für unsere Variablen ω und x. ω ist im Bereich von 0 und
1 definiert, x hingegen kann Werte im Bereich von −∞ bis +∞ annehmen. Wir betrachten
nur den Rand für Investor 1, da die Randwerte für Investor 2 analog zu berechnen sind.
Randwertbetrachtung für ω
ω=0
Es gilt für ω = 0, dass σω = 0 und µω = 0 ist. Allerdings kommen in unserer Gleichung
sowohl
1
ω
· σω als auch
1
ω
· µω vor. Deswegen müssen wir beim Randwert darauf achten, die
normalisierten Werte zu berechnen. Es gilt
42
1
ω σω
[(1−θ2 )
=
∂v2
∂v
−(1−θ1 ) ∂x1 ]σx
∂x
γ
und
σC1 = σc + ω1 σω
σC2 = σc
∂v1
∂ω σx
∂v2
∂ω σx
σv1 =
σv2 =
Für den normalisierten Wert von µω gilt
1
ω µω
1
ω σω σc
=
1
ψ1 )µC
0 σ
−0, 5γ(1 + ψ12 )σC
2 C2
1
0
+0, 5γ(1 + ψ1 )σC1 σC1
+ ψ1 · ( ψ12 −
0 σ
− 0, 5(1 − θ2 )σv0 2 σv2 − (1 − θ2 )σC
2 v2
0 σ
− 0, 5(1 − θ1 )σv0 1 σv1 − (1 − θ1 )σC
,
v
1
1
außerdem
µC1 = µC + Xt + ω1 µω + ω1 σω σC
2
∂ v1 0
1
µv1 = − ∂v
∂x κXt + 0, 5 ∂x2 σx σx .
Dieses eingesetzt in unsere Differentialgleichung (5.1) und umgeformt nach v1 , liefert unseren
Rand für ω = 0
v1 = ln
−β + (1 −
1
ψ1 )(µC1
−1
− 0.5γσC1 σC1 ) + µv1 + 0.5θ1 σv1 σv1 + (1 − γ)σC1 σv1
!
. (5.3)
ω=1
Der Randwert für ω = 1 ist wesentlich einfacher zu berechnen, da dort sowohl σω = 0, µω = 0
als auch für die normalisierten Werte
1
ω ·σω
=0,
1
ω ·µω
= 0 gelten. An relevanten Gleichungen
bleiben dann
σC1 = σc
∂v1
σx
∂x
= µC + Xt
σv1 =
µC1
µv1 = −κXt
∂v1
∂2v
+ 0, 5 2 σx σx .
∂x
∂x
5
Wohlstandskonsumquotient als partielle Differentialgleichung
43
Unsere nach v1 umgeformte Differentialgleichung lautet:
v1 = ln
−β + (1 −
1
ψ1 )(µC1
−1
− 0, 5γσC1 σC1 ) + µv1 + 0, 5θ1 σv1 σv1 + (1 − γ)σC1 σv1
!
(5.4)
Randwertbetrachtung von x:
x ist im Bereich von [−∞, +∞] definiert. Zwar liegen realistische Werte im Bereich um 0,
aber für die Randbetrachtung betrachten wir sehr kleine und sehr große Werte für x. Für sehr
große oder sehr kleine Werte für x können wir uns auf die Vorfaktoren von x beschränken,
da alle anderen Werte zu vernachlässigen sind.
Also gilt für x → −∞
0=1−
1
∂v1
−
κ,
ψ1
∂x
und analog für sehr große Werte für x (x → +∞) gilt:
0 = −1 +
1
∂v1
+
κ.
ψ1
∂x
Für Investor 2 können wir die Randwerte analog berechnen. Für die numerische Umsetzung
unserer Lösung nehmen wir für x realistischere Werte an, deshalb liefert die Randwertbetrachtung von x für uns keine relevanten Werte.
5.4 Ableitungen
Wir haben in unseren Differentialgleichungen verschiedene Ableitungen, die wir durch passende Differenzenquotienten approximieren müssen. In diesem Kapitel zeigen wir für welche
Ableitung wir welche Approximation benutzen und welche Ableitungen wir für die Ränder
wählen. Wir behandeln nur die Ableitungen von Investor 1, da sie für Investor 2 völlig analog
44
zu berechnen sind.
Als Notation gilt v(i, j) ≡ v1 (ωi , xj ) für i = 1, 2, 3, ..., nw und j = 1, 2, 3, ..., nx. Mit h1
beschreiben wir - wie bei der Diskretisierung gezeigt - die Schrittweite für die Ableitungen
nach ω und entsprechend mit h2 die Schrittweite der Ableitungen nach x.
Wir beginnen mit den inneren Ableitungen, das heißt sowohl i 6= 1 und i 6= nw als
auch j 6= 1 und j 6= nx. Für die einfachen Ableitungen
∂v(i,j)
∂ω
und
∂v(i,j)
∂x
betrachten wir das
in Kapitel 2.3 vorgestellte Upwind-Verfahren. Das heißt, dass wir je nach Vorzeichen des Vorfaktors den Rückwärts- oder den Vorwärtsdifferenzenquotienten nutzen. Also gilt für a· ∂v(i,j)
∂ω
bzw. b ·
∂v(i,j)
∂x
v(i, j) − v(i − 1, j)
∂v(i, j)
= −a ·
∂ω
h1
∂v(i, j)
v(i + 1, j) − v(i, j)
a·
= −a ·
∂ω
h1
a·
für a ≥ 0
(5.5)
für a < 0
(5.6)
für b ≥ 0
(5.7)
für b < 0
(5.8)
v(i, j) − v(i, j − 1)
∂v(i, j)
= −b ·
∂x
h2
∂v(i, j)
v(i, j + 1) − v(i, j)
b·
= −b ·
∂x
h2
b·
Für die reinen Ableitungen der 2. Ordnung nutzen wir den ebenfalls in Kapitel 2.3 vorgestellten Differenzenquotienten, der die Ableitung an drei Gitterpunkten approximiert.
∂ 2 v(i, j)
v(i + 1, j) − 2v(i, j) + v(i − 1, j)
=
∂ω 2
h21
(5.9)
∂ 2 v(i, j)
v(i, j + 1) − 2v(i, j) + v(i, j − 1)
=
∂x2
h22
(5.10)
Als nächstes betrachten wir die gemischten Ableitungen nach ω und x. Dafür nutzen wir die
Fünf-Punkte Stern Approximierung.
4v(i, j) − v(i, j + 1) − v(i + 1, j) − v(i, j − 1) − v(i − 1, j)
∂ 2 v(i, j)
=
∂ω∂x
h1 h2
(5.11)
Damit haben wir alle inneren Ableitungen vorgestellt. Für den Rand müssen wir beachten,
dass wir keine Punkte betrachten die außerhalb unseres Gitters liegen. Allerdings haben wir
auf den Rändern auch deutlich weniger Ableitungen, die wir berechnen müssen. Wir beginnen
mit der Ecke oben links.
5
Wohlstandskonsumquotient als partielle Differentialgleichung
45
i=1,j=1
∂v(1, 1)
v(1, 2) − v(1, 1)
=−
∂x
h2
2
∂ v(1, 1)
v(1, 3) − 2v(1, 2) + v(1, 1)
,
=
∂x2
h22
(5.12)
(5.13)
wobei wir in der ersten Zeile einmal und in der zweiten Zeile zweimal den Vorwärtsdifferenzenquotienten genutzt haben.
i=1,j∈[2,nx-1]
Da am oberen Rand nur Ableitungen nach x vorkommen, können wir sie wie oben mit
dem Upwind-Verfahren (5.7), (5.8) bei Ableitungen 1. Ordnung und mit der 3-PunkteApproximation (5.10) bei Ableitungen 2. Ordnung approximieren.
i=1,j=nx
∂v(1, nx)
v(1, nx) − v(1, nx − 1)
=−
∂x
h2
2
∂ v(1, nx)
v(1, nx) − 2v(1, nx − 1) + v(1, nx − 2)
=
∂x2
h22
(5.14)
(5.15)
Hier haben wir für die erste Ableitung einmal den Rückwärtsdifferenzenquotienten und für
die Ableitung 2. Ordnung zweimal den Rückwärtsdifferenzenquotienten genutzt.
i∈[2,nw-1],j=1
Die Ableitungen 1. Ordnung nach ω approximieren wir mit dem Upwind-Verfahren (5.5),
(5.6), die Ableitung 2. Ordnung nach ω wird wieder mit der 3-Punkte-Approximation (5.9)
dargestellt. Die Ableitung nach x bzw. zweimal nach x approximieren wir wie (5.12) bzw.
(5.13) nach dem Vorwärtsdifferenzenquotienten bzw. zweimal Vorwärtsdifferenzenquotienten.
Für die gemischten Ableitungen können wir nicht die Ableitung wie (5.11) definieren, da
wir sonst einen Punkt außerhalb des Gitters ansprechen. Deshalb nutzen wir für die Ableitung nach x den Vorwärtsdifferenzenquotienten und für die Ableitung nach ω das Upwind-
46
Verfahren.
∂ 2 v(i, 1)
v(i, 2) − v(i, 1) − v(i − 1, 2) + v(i − 1, 1)
=
∂ω∂x
h1 h2
2
∂ v(i, 1)
v(i + 1, 2) − v(i + 1, 1) − v(i, 2) + v(i, 1)
=
∂ω∂x
h1 h2
für positiven Vorfaktor
für negativen Vorfaktor
i∈[2,nw-1],j=nx
Die Ableitungen 1. und 2. Ordnung nach ω können wir wie oben behandeln (UpwindVerfahren (5.5), (5.6), 3-Punkte-Approximation (5.9)). Für die Ableitung 1. Ordnung nach x
bzw. 2. Ordnung nach x nutzen wir den Rückwärtsdifferenzenquotienten einmal bzw. zweimal. Für die gemischte Ableitung nutzen wir das Upwind-Verfahren (Ableitung nach ω) und
den Rückwärtsdifferenzenquotienten (Ableitung nach x).
∂v(i, nx)
v(i, nx) − v(i, nx − 1)
=−
;
∂x
h2
v(i, nx) − 2v(i, nx − 1) + v(i, nx − 2)
∂ 2 v(i, nx)
=
2
∂x
h22
v(i, nx) − v(i, nx − 1) − v(i − 1, nx) + v(i − 1, nx − 1)
∂ 2 v(i, nx)
=
∂ω∂x
h1 h2
(5.16)
(5.17)
für positiven Vorfaktor
(5.18)
∂ 2 v(i, nx)
v(i + 1, nx) − v(i + 1, nx − 1) − v(i, nx) + v(i, nx − 1)
=
∂ω∂x
h1 h2
für negativen Vorfaktor
(5.19)
i=nw,j=1
In dieser Ecke gibt es nur die Ableitungen 1. und 2. Ordnung nach x. Diese approximieren
wir wie (5.12) und (5.13) mit dem Vorwärtsdifferenzenquotienten einmal bzw. zweimal.
i=nw,j∈[2,nx-1]
Die einzigen Ableitungen die an diesem Rand relevant sind, sind wieder die Ableitungen 1.
und 2. Ordnung nach x, diesmal nutzen wir das Upwind-Verfahren (5.7), (5.8) und die 3Punkte-Approximation (5.10)
i=nw,j=nx
Die beiden Ableitungen 1. und 2. Ordnung nach x approximieren wir mit dem Rückwärtsdifferenzenquotienten einmal (5.14) bzw. zweimal (5.15).
Nun haben wir alle Ableitungen approximiert und können jetzt zur Umsetzung unserer Idee
übergehen.
5
Wohlstandskonsumquotient als partielle Differentialgleichung
47
5.5 Numerische Simulation
Die Grundidee unseres Anfangswertproblems haben wir weiter oben bereits vorgestellt. Nun
werden wir mit Hilfe von Matlab R2012b ein Programm entwickeln, mit dem wir konkrete
Beispiele berechnen können. Dafür nutzen wir insgesamt 3 verschiedene Programme. Das
Hauptprogramm EE (explizites Eulerverfahren) liefert die For-Schleife, die uns die benötigte
Anzahl von Schritten liefert. Außerdem aktiviert es die beiden anderen Programme Start und
DGL. Start berechnet unsere DGLv1 und DGLv2 inklusive der Randwerte. Die Randwerte
von ω werden exakt berechnet, die inneren Werte sind erst die erste Näherung. Im Programm
DGL, was in der For-Schleife von EE ausgeführt wird, berechnen wir nur die inneren Werte
von DGLv1 und DGLv2, da die Randwerte ja schon exakt sind.
Wir beginnen damit, an EE die benötigten
Parameter β, ψ1 , ψ2 , µC , γ, σC , σx , κ, ρ, ωmin , ωmax , xmin , xmax , h1 , h2 und τ zu übergeben.
[v1,v2]=EE(beta,psi1,psi2,muC,gamma,sigmaC,sigmax,kappa,rho,
w min,w max,x min,x max,h1,h2,t)
Beispiel:
[v1,v2]=EE(0.02,0.5,1.2,0.02,6,0.0252,0.000085,0.008799,0.5,
0,1,−0.025,0.025,0.1,0.005,0.1)
Unser Programm beginnt mit der Berechnung der maximalen Schrittanzahl.
1 function [v1,v2]=EE(beta,psi1,psi2,muC,gamma,sigmaC,sigmax,kappa,rho,
w min,w max,x min,x max,h1,h2,t)
2
3 nw=(w max−w min)/h1+1
4 nx=(x max−x min)/h2+1
5
6
7 % Startwerte fuer v1 und v2, in unserem Fall beginnen wir mit
8 % v1=1 und v2=1. Die Zahlen koennen wir beliebig waehlen, wir duerfen
9 % aber nicht oberhalb der zu berechnen Werte fuer v1 und v2 liegen.
10 v1=ones(nw,nx);
11 v2=ones(nw,nx);
12
13 % Als erstes berechnen wir mit unseren Startwerten v1 und v2
14 % die DGLv1 und DGLv2 inklusive der Randwerte. Dafuer
15 % uebergeben wir alle benoetigen Parameter an das Programm Start.
48
16 [DGLv1,DGLv2]=Start(v1,v2,beta,psi1,psi2,muC,gamma,sigmaC,sigmax,kappa,rho,
w min,w max,x min,x max,h1,h2);
17
18 % Dieses liefert die erste Annaeherung. Wir setzen die berechneten Werte als
19 % Naeherung fuer v1 und v2 ein. Die Randwerte sind schon korrekt, die anderen
20 % sind noch weit von der richtigen Loesung entfernt.
21 v1=DGLv1;
22 v2=DGLv2;
23
24 % Nun beginnt der wichtigste Teil des Programms. Die For−Schleife laeuft
25 % sehr haeufig durch und loest damit unser Anfangswertproblem. DGLv1 und
26 % DGLv2 streben gegen 0, v1 und v2 streben gegen die korrekte Loesungen
27 % der Differentialgleichungen. Die Anzahl der Schleifendurchlaeufe
28 % darf nicht zu niedrig sein, da wir sonst keine gute Naeherung erreichen.
29 for i=1:1000000
30
[DGLv1,DGLv2]=DGL(v1,v2,beta,psi1,psi2,muC,gamma,sigmaC,sigmax,kappa,rho,
w min,w max,x min,x max,h1,h2);
31
v1=v1+t*DGLv1;
32
v2=v2+t*DGLv2;
33 end;
34
35 end
Damit endet unser Hauptprogramm. Als Ausgabe erhalten wir die Lösung für v1 und v2 . Die
numerisch bestimmte Lösung ist eine gute Annaehrung an die exakte Lösung, falls wir die
For-Schleife oft genug durchlaufen lassen haben - dazu später mehr. Unser Programm Start
wird nur einmal benötigt, das Wichtigste an ihm ist die Berechnung der Randwerte für ω = 0
und ω = 1. Es enthält außerdem die Routine zur Berechnung von DGLv1 und DGLv2 an
allen inneren Punkten. Wir stellen die wichtigsten Punkte des Programms kurz vor.
1 function [DGLv1,DGLv2]= Start(v1,v2,beta,psi1,psi2,muC,gamma,sigmaC,sigmax,
kappa,rho,w min,w max,x min,x max,h1,h2,nw,nx)
2
3 % Wir erzeugen als erstes unser Gitter auf dem wir unsere
4 % Differentialgleichungen loesen moechten.
5 [w,x]=ndgrid(w min:h1:w max,x min:h2:x max);
6
7 % Theta i koennen wir laut Definitionen berechnen.
8 theta1=(1−gamma)/(1−1/psi1);
9 theta2=(1−gamma)/(1−1/psi2);
10
Wir berechnen alle relevanten Ableitungen, genau wie in Kapitel 5.4 erklärt. Aus Übersichtsgründen verzichten wir darauf, diese hier darzustellen.
196 % Die relevanten Raender sind die Raender fuer w=0 und w=1,
5
Wohlstandskonsumquotient als partielle Differentialgleichung
197
198
199
200
201
202
204
205
206
207
208
209
210
211
212
49
% wir beginnen mit w=0 also i=1. Wir berechnen die Randwerte genau
% wie in Kapitel 5.3 dargestellt.
i=1;
for j=1:nx
% D1 ist der normalisierte Wert fuer sigma w also 1/w * sigma w.
D1(i,j) = (((1−theta2)*v2x(i,j)−(1−theta1)*v1x(i,j))*sigmax)/(gamma);
sigmaC1(i,j)=sigmaC+D1(i,j);
sigmaC2(i,j)=sigmaC;
sigmav1(i,j)=v1x(i,j)*sigmax;
sigmav2(i,j)=v2x(i,j)*sigmax;
% Auch bei mu w betrachten wir die normalisierten Werte 1/w* mu w.
V1(i,j) = −D1(i,j)*sigmaC+ psi1*((1/psi2−1/psi1)*muC
−0.5 *gamma *(1+1/psi2)*sigmaC2(i,j)*sigmaC2(i,j)
−0.5*(1−theta2)*sigmav2(i,j)*sigmav2(i,j)
−(1−theta2)*sigmaC2(i,j)*sigmav2(i,j)
+0.5* gamma* (1+1/psi1)*sigmaC1(i,j)*sigmaC1(i,j)
+0.5*(1−theta1)*sigmav1(i,j)*sigmav1(i,j)
+(1−theta1)*sigmaC1(i,j)*sigmav1(i,j));
214
215
muC1(i,j) = muC+x(i,j) +V1(i,j)+D1(i,j)*sigmaC;
216
muC2(i,j) = muC+x(i,j);
217
218
muv1(i,j)=−kappa*x(i,j)*v1x(i,j)+0.5*v1xx(i,j)*sigmax*sigmax;
219
muv2(i,j)=−kappa*x(i,j)*v2x(i,j)+0.5*v2xx(i,j)*sigmax*sigmax;
220
221 DGLv1(i,j)=log(−1/(−beta + (1−1/psi1)*(muC1(i,j)
− 0.5*gamma*sigmaC1(i,j)*sigmaC1(i,j))+muv1(i,j)
+ 0.5* theta1* sigmav1(i,j)*sigmav1(i,j)
+ (1−gamma)* sigmaC1(i,j)*sigmav1(i,j)));
222 DGLv2(i,j)=log(−1/(−beta + (1−1/psi2)*(muC2(i,j)
− 0.5*gamma*sigmaC2(i,j)*sigmaC2(i,j))+muv2(i,j)
+ 0.5* theta2* sigmav2(i,j)*sigmav2(i,j)
+ (1−gamma)* sigmaC2(i,j)*sigmav2(i,j)));
223
224 end;
Anschließend berechnen wir den Rand fuer w=1, also i=nw, recht analog, worauf wir an
dieser Stelle nicht näher eingehen werden.
257
258
259
260
261
262
% Als letztes berechnen wir fuer DGLv1
% und DGLv2 die "inneren" Werte, also fuer i=2,..., nw−1 und
% j=1,...,nx
for i=2:nw−1
for j=1:nx
50
263
264 % Wir starten mit der Berechnung von sigmaw, da wir je nach Vorfaktor
265 % der Ableitungen verschiedene Differenzenquotienten nutzen. Dafuer
266 % berechnen wir zuerst den Nenner von sigmaw:
267
268
Nenner(i,j)= (gamma+w(i,j) * (1−w(i,j))*
(max((1−theta1),0)*v1wp(i,j)
+min((1−theta1),0)*v1wn(i,j)
−max((1−theta2),0)*v2wp(i,j)
−min((1−theta2),0)*v2wn(i,j)));
269
270
sigmaw(i,j) = w(i,j)*(1−w(i,j)) *
(max((1−theta2)/(Nenner(i,j)),0)*v2xp(i,j)
+min((1−theta2)/(Nenner(i,j)),0)*v2xn(i,j)
−max((1−theta1)/(Nenner(i,j)),0)*v1xp(i,j)
−min((1−theta1)/(Nenner(i,j)),0)*v1xn(i,j))*sigmax;
271
sigmaw2(i,j)= −sigmaw(i,j);
272
273
sigmaC1(i,j)=sigmaC+1/(w(i,j))*sigmaw(i,j);
274
sigmaC2(i,j)=sigmaC−1/(1−w(i,j))*sigmaw(i,j);
275
276
sigmav1(i,j)=max(sigmaw(i,j),0)*v1wp(i,j)
+min(sigmaw(i,j),0)*v1wn(i,j)
+max(sigmax,0)*v1xp(i,j)
+min(sigmax,0)*v1xn(i,j);
277
sigmav2(i,j)=−max(sigmaw(i,j),0)*v2wp(i,j)
−min(sigmaw(i,j),0)*v2wn(i,j)
+max(sigmax,0)*v2xp(i,j)
+min(sigmax,0)*v2xn(i,j);
278
279
muw(i,j)= − sigmaw(i,j)*sigmaC +
(w(i,j)*(1−w(i,j))*psi1*psi2)/(w(i,j)*psi1+(1−w(i,j))*psi2)*
((1/psi2 − 1/psi1)*muC
−0.5*gamma*(1+1/psi2)*sigmaC2(i,j)*sigmaC2(i,j)
−0.5*(1−theta2)*sigmav2(i,j)*sigmav2(i,j)
−(1−theta2)*sigmaC2(i,j)*sigmav2(i,j)
+0.5*gamma*(1+1/psi1)*sigmaC1(i,j)*sigmaC1(i,j)
+0.5*(1−theta1)*sigmav1(i,j)*sigmav1(i,j)
+(1−theta1)*sigmaC1(i,j)*sigmav1(i,j));
280
281
muv1(i,j)=max(muw(i,j),0)*v1wp(i,j)+min(muw(i,j),0)*v1wn(i,j)
−max(kappa*x(i,j),0)*v1xp(i,j)
−min(kappa*x(i,j),0)*v1xn(i,j)
+0.5*v1ww(i,j)*sigmaw(i,j)*sigmaw(i,j)
+0.5*v1xx(i,j)*sigmax*sigmax
+rho*v1wx(i,j)*sigmaw(i,j)*sigmax;
282
muv2(i,j)=max(muw(i,j),0)*v2wp(i,j)+min(muw(i,j),0)*v2wn(i,j)
−max(kappa*x(i,j),0)*v2xp(i,j)
−min(kappa*x(i,j),0)*v2xn(i,j)
5
Wohlstandskonsumquotient als partielle Differentialgleichung
51
+0.5*v2ww(i,j)*sigmaw(i,j)*sigmaw(i,j)
+0.5*v2xx(i,j)*sigmax*sigmax
−rho*v2wx(i,j)*sigmaw(i,j)*sigmax;
283
284
285
muC1(i,j) = muC+x(i,j) +1/(w(i,j))*muw(i,j)
+1/(w(i,j))*sigmaw(i,j)*sigmaC;
muC2(i,j) = muC+x(i,j) −1/(1−w(i,j))*muw(i,j)
−1/(1−w(i,j))*sigmaw(i,j)*sigmaC;
286
287
% Nun berechnen wir unsere DGLv1 und DGLv2 fuer die inneren Werte
288
% unseres Gitters.
289 DGLv1(i,j)=exp(−v1(i,j)) −beta + (1−1/psi1)*(muC1(i,j)
− 0.5*gamma*sigmaC1(i,j)*sigmaC1(i,j))+muv1(i,j)
+ 0.5* theta1* sigmav1(i,j)*sigmav1(i,j)
+ (1−gamma)* sigmaC1(i,j)*sigmav1(i,j);
290 DGLv2(i,j)=exp(−v2(i,j)) −beta + (1−1/psi2)*(muC2(i,j)
− 0.5*gamma*sigmaC2(i,j)*sigmaC2(i,j))+muv2(i,j)
+ 0.5* theta2* sigmav2(i,j)*sigmav2(i,j)
+ (1−gamma)* sigmaC2(i,j)*sigmav2(i,j);
298
end;
299 end;
300 end
Unser Programm DGL ist eine verkürzte Version unseres Programms Start. Es vernachlässigt
die Ränder - da diese schon korrekt sind - und berechnet nur die inneren Werte für DGLv1
und DGLv2, genau wie das Programm Start.
52
Abbildung 5.1: Die mit Hilfe unseres Programms berechnete Lösung für v1 auf einem 10x10
Gitter.
Abbildung 5.2: Die mit Hilfe unseres Programms berechnete Lösung für v2 auf einem 10x10
Gitter.
Mit Hilfe des Programms haben wir unser Beispiel berechnet. Dabei haben wir einmal als
Schrittweite für das explizite Eulerverfahren τ = 0, 1 und für das Gitter ein 10x10 Gitter gewählt. Die For-Schleifen Anzahl haben wir auf 10.000 gesetzt. Abbildung 5.1 zeigt den
Graph der Lösung unseres Beispiels von v1 . Abbildung 5.2 zeigt den Graph der Lösung von v2 .
5
Wohlstandskonsumquotient als partielle Differentialgleichung
53
Abbildung 5.3: Die berechnete Funktion v1 bei einer Erhöhung der Diskretisierungschritte
auf ein 100x100 Gitter.
Abbildung 5.4: Die berechnete Funktion v2 bei einer Erhöhung der Diskretisierungschritte
auf ein 100x100 Gitter.
Als zweites haben wir mit dem Programm das selbe Beispiel nochmal berechnet, aber bestimmte Parameter etwas verändert. Wir haben die Schrittweite τ = 1 gesetzt, außerdem
haben wir unsere Funktionen auf einem 100x100 Gitter berechnet. Die For-Schleifen Anzahl
haben wir dagegen auf 2000 gesenkt. Wie wir weiter unten sehen werden, reicht es bei beiden
Berechnungen aus den maximalen Fehler kleiner als 10−6 werden zu lassen. Die Abbildung
5.3 zeigt den Graph der Lösung von v1 auf dem 100x100 Gitter. Die Abbildung 5.4 zeigt den
Graph der Lösung von v2 auf der selbigen Gittergröße.
Wie zu sehen ist, werden die maximalen Werte immer am Rand unseres Gitters angenommen.
Dies entspricht dem Maximumprinzip. Außerdem wird auch der minimale Wert am Rand
angenommen. Mit diesen Gegebenheiten können wir eine kleine Optimierung an unserem
Programm vornehmen: nach dem Berechnen der Randwerte definieren wir alle inneren Werte
der Funktion gleich dem kleineren der beiden Randwerte. Die Optimierung sieht für die erste
Funktion wie folgt aus:
54
29
30
31
32
33
34
35
36
37
38
39
for j=1:nx
if DGLv1(1,j) < DGLv1(nw,j)
for i=2:nw−1
v1(i,j)=DGLv1(1,j);
end;
else
for i=2:nw−1
v1(i,j)=DGLv1(nw,j);
end;
end;
end;
Diese kleine Programmänderung haben wir in dem Programm EEmitbesserenStart verwirklicht. Es spart etwa 10% der benötigten For-Schleifen-Durchläufe um eine spezielle Genauigkeit zu erlangen. Insgesamt ist das Programm recht rechenintensiv, da viele Durchläufe der
For-Schleife nötig sind. Vor allem sollte, wie gleich zu sehen ist, die Schrittweite τ von der
Gitterpunkteanzahl abhängen.
Unser Verfahren ist ein Vorwärts-Euler-Verfahren, das in [Bur07, S.9f] nachzulesen ist. Die
Existenz und Eindeutigkeit unserer Lösung ist aufgrund des Verfahrens schnell gegeben. Allerdings müssen wir zur Stabilität noch folgendes beachten: Der Fehler kann geometrisch
anwachsen. Um dies zu vermeiden, muss unser τ sehr klein im Verhältnis zur Gitterpunkteanzahl gewählt werden. Durch das kleine τ müssen wir dann aber die Anzahl der For-SchleifenDurchläufe deutlich erhöhen. Das heißt, um ein stabiles Verfahren für eine vernünftige Anzahl
an Gitterpunktem zu bekommen, ist die Rechenzeit sehr groß. Diesen Aspekt werden wir nun
genauer untersuchen.
Genauigkeit
Als nächstes beschäftigen wir uns in diesem Abschnitt mit der Genauigkeit, also der Fehlergröße unseres Verfahrens. Dafür haben wir den maximalen Fehler unsers Beispiels für verschiedene Schrittweiten τ und Gitterpunkte berechnet. Wir nehmen an, dass es sich bei unserem
Gitter um ein quadratisches Gitter handelt, das bedeutet nw = nx. In der Tabelle 5.1 sieht
man die Entwicklung des maximalen Fehlers in Abhängigkeit der beiden Größen τ und nw
bei einem For-Schleifen-Durchlauf von 10.000. Hierbei wird folgendes deutlich: je kleiner τ
und je größer nw, desto größer der Fehler. Wobei das kleinere τ größeren Einfluss auf den
maximalen Fehler hat als das größere nw. Das unser Ergebnis keine gute Annäherung liefert,
liegt vor allem daran, dass der For-Schleifen-Durchlauf mit 10.000 nicht groß genug ist.
5
Wohlstandskonsumquotient als partielle Differentialgleichung
τ =1
3x3
< 10−15
6x6
< 10−15
11x11
< 10−15
21x21
< 10−15
τ = 0, 5
< 10−14
< 10−14
< 10−14
< 10−14
τ = 0, 1
< 10−9
< 10−8
< 10−7
< 10−7
τ = 0, 05
< 10−5
< 10−4
< 10−4
< 10−4
τ = 0, 01
< 10−2
< 10−2
< 10−2
< 10−2
55
Tabelle 5.1: Der maximale Fehler in Abhängigkeit der Schrittweite und Gitterdiskretisierung
Um die Stabilität zu garantieren, ist es wichtig τ möglichst klein zu wählen. Im Zusammenhang hiermit machen wir uns jetzt Gedanken darüber welche Schrittanzahl (For-SchleifenDurchläufe) wir benötigen um ein vernünftiges Ergebnis zu erzielen. Wir schreiben unser Programm ein wenig um, indem wir ein Abbruchkriterium einbauen, bei dem wir die gewünschte
Genauigkeit der Lösung verlangen können. Das heißt, unser Programm bricht automatisch
ab, sobald die DGLv1 und DGLv2 kleiner als eine gewünschte Genauigkeit (z.B.: 10−6 ) sind.
Dadurch ist sichergestellt, dass unser Programm auch für kleine τ und große nw noch ein
vernünftiges Ergebnis erzielt. Natürlich erhöht sich dadurch die For-Schleifen-Anzahl deutlich.
Im nachfolgenden Abschnitt haben wir die benötigten For-Schleifen-Durchläufe untersucht.
For-Schleifen-Durchläufe
Mit den berechneten Werten von DGLv1 und DGLv2 in jeder For-Schleife unseres Hauptprogramms EE, können wir unser Programm auch von einem anderen Standpunkt betrachten.
Wir können betrachten, für welche gewünschte Genauigkeit, bei welcher Schrittgröße τ und
einer gewöhnlichen Gittergröße 10x10, die Anzahl der For-Schleifen-Durchläufe wie groß sein
muss. Dafür haben wir unser Programm EEmitbesserenStart um eine Abbruchkomponente
ergänzt.
45
46
47
48
49
50
if max(max((abs(DGLv1))))< 10ˆ(−6);
if max(max(abs(DGLv2))) < 10ˆ(−6);
benoetigteSchritte=i
break
end;
end;
Unser Programm, welches wir nun EEmitabbruch nennen, gibt außer den beiden Lösungen für
v1 und v2 auch noch die benötigte Schrittanzahl aus (strenggenommen müssen wir dazu noch
1 addieren, da wir vor der For-Schleife bereits die erste Näherung mit Hilfe des Programms
56
< 10−3
227
< 10−4
405
< 10−5
562
< 10−6
714
< 10−7
864
< 10−8
1015
τ = 0, 5
453
812
1127
1432
1734
2036
τ = 0, 1
2265
4063
5645
7175
8692
10206
τ = 0, 05
4529
8128
11293
14354
17389
20419
τ = 0, 01
22646
40646
56476
71785
86920*
102060*
τ =1
Tabelle 5.2: Die Schrittanzahl in Abhängigkeit der Schrittweite und Genauigkeit.
Start berechnet haben. Da wir uns aber nicht für die genaue Schrittanzahl, sondern nur für
die Größenordnung interessieren, können wir diese Ungenauigkeit vernachlässigen). In der
Tabelle 5.2 können wir für verschiedene Werte von τ und verschiedene Genauigkeiten die
Anzahl der For-Schleifen-Durchläufe betrachten.
Wie man deutlich sieht, verhält sich die Schrittanzahl bei geringerer Schrittweite τ genau
wie es zu erwarten wäre (halbe Schrittweite = doppelte Schrittanzahl). Auf dieser Tatsache
haben wir die mit * gekennzeichneten Werte in der Tabelle geschätzt. Mit Hilfe unseres
Programms EEmitabbruch können wir auch die benötigten Schritte für z.B.: nw = nx = 100
mit τ = 1 und g < 10−6 berechnen - diese liegen bei 915. Anhand dieser Informationen können
wir abschätzen, in welcher Größenordnung die Anzahl der Schritte bei z.B.: nw = 100 und
τ = 0, 0001 und einer Genauigkeit von < 10−6 liegen müßte: bei ca. 9 · 106 .
Folglich verursacht die Berechnung eines stabilen Verfahrens für ein großes Gitter einen enormen Rechenaufwand.
Probleme
Auch bei diesen Programmen gibt es dieselben Probleme wie in Kapitel 4. Wieder können
die Randwerte von ω, und damit alle Werte von v1 und v2 , irregulär werden. Für die meisten realistischen Werte der Parameter lässt sich dieses Problem vermeiden. Aber neben der
möglichen zu großen Differenz von ψ1 und ψ2 wie in Kapitel 4, kann auch die zusätzliche
Variable x die Irregulärität der Lösung verursachen. Im folgenden Kapitel werden wir kurz
darauf eingehen wie diese Probleme vermieden werden könnten.
57
6 Fazit und Ausblick
In dieser Diplomarbeit wurde ein Lucas-Tree-Modell mit zwei heterogenen Investoren numerisch analysiert. Die Herleitung der Differentialgleichungen erfolgte über die Epstein-ZinNutzenfunktion. Die zwei heterogenen Investoren unterscheiden sich in ihren Substitutionsfaktoren, ihrem Zeitpräferenzfaktor und ihrem Risikoaversionsfaktor.
In Kapitel 4 wurden die beiden letztgenannten Heterogenitätseigenschaften der Investoren
aufgehoben. Außerdem wurden die Differentialgleichungen auf gewöhnliche Differentialgleichungen beschränkt. Diese wurden dann mit Hilfe des Upwind-Verfahrens zur Approximation
von Ableitungen und dem Newton-Verfahren numerisch berechnet. Die numerische Berechnung benutzte dafür die beiden Randwerte der Funktionen. Das numerische Verfahren lieferte
sehr schnell eine gute Lösung für die Differentialgleichungen. Dieses Verfahren wurde in Hinblick auf Geschwindigkeit und maximaler Fehler an einem konkreten Beispiel betrachtet.
Im 5. Kapitel blieb die partielle Eigenschaft der Differentialgleichungen erhalten. Gezeigt
wurde, dass es sich bei den Differentialgleichungen 2. Ordnung um elliptische Differentialgleichungen handelte. Mit Hilfe einer Diskretisierung und verschiedenen Differenzenquotienten
wurde die Berechnung der Differentialgleichungen vorbereitet. Zur numerischen Berechnung
der Differentialgleichungen wurde das Vorwärts-Euler-Verfahren benutzt. Dieses lieferte für
die zweidimensionalen Funktionen eine gute Lösung. Für eine hohe Genauigkeit und ein stabiles Verfahren wurde aber deutlich, dass der Rechenaufwand sehr hoch ist. Es wurden ein
paar kleine Variationen an der numerischen Berechnung vorgestellt.
Beim Vergleich der beiden Verfahren aus Kapitel 4 und Kapitel 5 fällt auf, dass das UpwindNewton-Verfahren deutlich schneller zur gewünschten Lösung führt. Das Vorwärts-EulerVerfahren ist aber aufgrund seiner höheren Variationsmöglichkeiten (Vorgabe des maximalen
Fehlers etc.) und seiner geringeren Einschränkung in den meisten Fällen zu empfehlen.
Sowohl in Kapitel 4, als auch in Kapitel 5 gab es für bestimmte Parameter Probleme bei
der numerischen Berechnung. Diese Probleme enstehen, da ev als Funktion in der Herleitung
festgelegt wurde. Um die Exponentialfunktion zu lösen wurde der natürliche Logarithmus genutzt, welcher bei der Berechnung von negativen Zahlen komplexe Zahlen ergibt. Zur Vermeidung von komplexen Zahlen sollte die partielle Differentialgleichung mit v statt ev hergeleitet
werden. Dadurch könnten diese Probleme vermieden werden.
Für zukünftige Untersuchungen wäre es sinnvoll die Differentialgleichungen nicht mit Hilfe
58
der Exponentialfunktionen herzuleiten. So könnte die Vermeidung der Probleme aus Kapitel
4 und 5 gewährleistet werden. Außerdem könnten in zukünftigen Untersuchungen die Einschränkungen der Gleichheit des Risikoaversionsfaktors und der Zeitpräferenzrate vermieden
werden. Ebenso könnte weiterhin analysiert werden, ob es für die numerische Berechnung aus
Kapitel 5 ein weniger rechenintensives Verfahren gibt.
59
Abbildungsverzeichnis
2.1
Auszahlungsfunktion einer Call- und einer Put-Option mit Ausübungspreis
K=100 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.2
Der Pfad eines Wiener Prozesses Zt in der Zeitspanne von t=0 bis t=1. . . .
10
4.1
Der grüne Graph stellt die Funktion v1 und der rote Graph die Funktion v2 dar. 31
4.2
Die maximale Newton-Schrittanzahl unter den Parametern Genauigkeit und
Schrittanzahl für ω.
4.3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
Der maximale Fehler für v1 roter Graph und für v2 grüner Graph von unserem Beispiel in Abhängigkeit von der Schrittanzahl [2:100] bei einer NewtonGenauigkeit von 10−6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4
Der maximale Fehler für v1 unseres Beispiels in Abhängigkeit von der Newtongenauigkeit und der Schrittanzahl. . . . . . . . . . . . . . . . . . . . . . .
5.1
52
Die berechnete Funktion v1 bei einer Erhöhung der Diskretisierungschritte auf
ein 100x100 Gitter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4
52
Die mit Hilfe unseres Programms berechnete Lösung für v2 auf einem 10x10
Gitter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3
34
Die mit Hilfe unseres Programms berechnete Lösung für v1 auf einem 10x10
Gitter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2
33
53
Die berechnete Funktion v2 bei einer Erhöhung der Diskretisierungschritte auf
ein 100x100 Gitter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
60
Tabellenverzeichnis
2.1
Call-/Put-Option - Vergleich Preis des Gutes und festgelegter Ausübungspreis
4.1
Der maximale Fehler von v1 (von unserem Beispiel) in Abhängigkeit der Newtongenauigkeit und der Schrittanzahl. . . . . . . . . . . . . . . . . . . . . . .
5
34
5.1
Der maximale Fehler in Abhängigkeit der Schrittweite und Gitterdiskretisierung 55
5.2
Die Schrittanzahl in Abhängigkeit der Schrittweite und Genauigkeit. . . . . .
56
61
Literaturverzeichnis
[BCDG10] Luca Benzoni, Pierre Collin-Dufresne, and Robert Goldstein. Explaining Asset
Pricing Puzzles Associated with the 1987 Market Crash. paper, 2010. 16, 17, 19,
20
[BDIS12]
Nicole Branger, Ioana Dumitrescu, Vesela Ivanova, and Christian Schlag. Preference Hetereogeneity and Survival in Long-Run Risk Models paper, 2012. 16
[Bur07]
Martin Burger. Numerik partieller Differentialgleichungen. Skript, 2007. 8, 12,
13, 14, 26, 36, 38, 39, 54
[Bur10]
Martin Burger. Mathematische Modellierung. Skript, 2010. 11, 12
[Dö12]
Michael Dörr. Numerische Analyse von LRR-Modellen mit zwei Lucas-Bäumen
und Heterogenen Investoren. Diplomarbeit, 2012. 3
[Hä12]
Johannes Härtel. Numerische Analyse von Long Run Risk Modellen mit zwei
Bäumen und Sprungrisiko. Diplomarbeit, 2012. 3, 19
[Hul09]
John C. Hull.
Optionen, Futures und andere Derivate.
Pearson Studium,
München, 7th edition, 2009. 3
[Kle08]
Achim Klenke. Wahrscheinlichkeitstheorie. Springer Verlag, 2008. 8
[KM74]
Helmut Kieswetter and Gerhard Maeß. Elementare Methoden der numerischen
Mathematik. 1974. 14
[Mun09]
Claus Munk. Financial asset pricing theory, working paper, 2009. 8
[Nat01]
Frank Natterer. Partielle Differentialgleichungen, Skript, 2001. 37
[Ohl08]
Mario Ohlberger. Höhere Numerische Mathematik. Skript, 2008. 40
[Zoc10]
Jann-Philipp Zocher. Numerische Berechnung von Preis-Konsum Quotienten in
verallgemeinerten Gleichgewichtsmodellen, Diplomarbeit, 2010. 3