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
© Copyright 2024 ExpyDoc