Seminar: Grundlagen der Simulation und Statistik
von dynamischen Systemen
SS 2012
Maximum-Likelihood-Schätzung für das
Black-Scholes-Merton-Modell und das
Cox-Ingersoll-Ross-Modell
Thema 10
Philipp Probst
19. Juni 2012
Dozentin:
Prof. Dr. Christine Müller
Inhaltsverzeichnis
1 Einleitung 3
2 Statistische Modelle und Methoden 3
2.1 Die stochastische Differentialgleichung . . . . . . . . . . . . . . . . . . . . 4
2.2 Das Black-Scholes-Merton-Modell . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 Das Cox-Ingersoll-Ross-Modell . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Die Maximum-Likelihood-Schätzung . . . . . . . . . . . . . . . . . . . . . 5
2.5 Die ML-Schätzung bei stochastischen DGLs . . . . . . . . . . . . . . . . . 6
3 Simulationen 7
3.1 Black-Scholes-Merton-Modell . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2 Cox-Ingersoll-Ross-Modell . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4 Zusammenfassung 13
A Anhang 15
Literatur 25
1 Einleitung 3
1 Einleitung
Stochastische Integral- und Differentialrechnung spielen in den Wirtschaftswissenschaften
eine immer größer werdende Rolle. Sie kann zum einen zur Modellierung von Finanz-
märkten eingesetzt werden, zum anderen baut die moderne Zeitreihenökonometrie darauf
auf. Das von Black und Scholes aufgestellte Modell zur Lösung einer stochastischen Dif-
ferentialgleichung brachte den beiden sogar den Nobelpreis und fand Anwendung auf
Aktienmärkten. Es sollte dabei nicht außer acht gelassen werden, dass Finanzmärkte sehr
schwierig zu modellieren sind und bestimmte Annahmen wie die, dass die Fluktuationen
normalverteilt sind, nicht mit der Realität übereinstimmen. Solche Fehlannahmen haben
in der Vergangenheit zu großen Verlusten im Derivatenhandel geführt (Ball, 2004, S.251).
Lösungen von anderen stochastischen Differentialgleichungen haben aber nicht nur An-
wendungen in der Wirtschaft, das Cox-Ingersoll-Ross-Modell wird beispielsweise auch zur
Modellierung von Bevölkerungswachstum oder Stickstoffoxidemissionen verwendet.
In der folgenden Arbeit werden in Kapitel 2 zunächst das Black-Scholes-Modell, das Cox-
Ingersoll-Ross-Modell sowie die Maximum-Likelihood-Schätzung vorgestellt. Daraufhin
werden in Kapitel 3 jeweils konkrete Möglichkeiten zur Maximum-Likelihood-Schätzung
der Parameter der beiden Modelle besprochen. Bei den darauf folgenden Simulationen zu
den beiden Modellen soll vor allem die Genauigkeit der Schätzungen bei unterschiedlicher
Anzahl an Realisierungen oder unterschiedlichen Zeitabständen der aufeinander folgen-
den Realisierungen ins Auge gefasst werden. Außerdem soll die asymptotische Normal-
verteilungsannahme der Maximum-Likelihood-Schätzung überprüft werden. In Kapitel 4
werden die wichtigsten Ergebnisse der Simulationen zusammengefasst.
2 Statistische Modelle und Methoden
In diesem Kapitel werden die statistischen Modelle und Methoden vorgestellt, die bei den
darauf folgenden Simulationen benutzt werden.
2 Statistische Modelle und Methoden 4
2.1 Die stochastische Differentialgleichung
Die Differentialgleichung dxt
dt= f(t, x) kann als Integralgleichung xt = xa +
∫ ta f(s, x(s))ds
geschrieben werden (Kuo, 2006, S.185). Wird die Funktion durch eine brownsche Bewe-
gung Bt beeinflusst, so ergibt sich die Differentialgleichung dXt
dt= f(t, X) + σ(t, X)B′
t.
Diese kann umgeschrieben werden zu dXt = f(t, X)dt + σ(t, X)dBt, die als stochastische
Integralgleichung Xt = Xa +∫ t
a f(s, Xs)ds +∫ t
a σ(s, Xs)dB(s), a ≤ t ≤ b interpretiert
werden kann. Lösungen Xt dieser Gleichungen werden Diffusionsprozesse genannt.
In den beiden folgenden Kapiteln werden zwei wichtige stochastische Differentialgleichun-
gen vorgestellt, die sich aus der allgemeinen Form
dXt = f(Xt, θ)dt + σ(Xt, θ)dB(t) (1)
ableiten lassen. θ ∈ Θ ∈ Rp ist dabei ein unbekannter p-dimensionaler Vektor, der die
Funktionen f und σ näher beschreibt und f : R × Θ → R und σ : R × Θ → (0, ∞) seien
bekannt und so gewählt, dass (1) eine Lösung besitzt.
2.2 Das Black-Scholes-Merton-Modell
Das Black-Scholes-Merton-Modell wird je nach Kontext auch Black-Scholes-Modell oder
“Modell der brownschen geometrischen Bewegung” genannt. Es löst die spezielle stochas-
tische Differentialgleichung
dXt = θ1Xtdt + θ2XtdBt, X0 = x0, θ1 ∈ R, θ2 ∈ R+. (2)
(vgl. Iacus, 2008, S.117f). Die zugehörige bedingte Dichtefunktion pθ(t, .|x) zum Zeitpunkt
t bei gegebenem Wert x ist logarithmisch normalverteilt mit Erwartungswert m(t, x) =
xeθ1t und Varianz v(t, x) = x2e2θ1t(eθ22t − 1). Daraus folgt die bedingte Dichte
pθ(t, y|x) =1
θ2y√
2πtexp
{
−(log y − (log x + (θ1 − 12θ2
2)t))2
2θ22t
}
. (3)
2 Statistische Modelle und Methoden 5
2.3 Das Cox-Ingersoll-Ross-Modell
Das Cox-Ingersoll-Ross-Modell löst die spezielle stochastische Differentialgleichung
dXt = (θ1 − θ2Xt)dt + θ3
√
XtdBt, X0 = x0 > 0, (4)
(vgl. Iacus, 2008, S.119), wobei θ1, θ2, θ3 ∈ R+ sind. Im Falle von 2θ1 > θ23 ist der Prozess
strikt positiv, ansonsten nur nichtnegativ. Die bedingte Dichte Yt|Y0 der Zufallsvariablen
Yt = 2cXt mit c = 2θ2/(θ23(1−e−θ2t)) verhält sich wie eine nichtzentrale χ2-Verteilung mit
4θ1/θ23 Freiheitsgraden und Nichtzentralitätsparameter y0e
−θ2t (Iacus, 2008, S.47). Daraus
lässt sich leicht die bedingte Dichte von Xt|X0 = x0 berechnen. Das Ergebnis ist
pθ(t, y|x) = ce−u−v(
v
u
)q/2
Iq(2√
uv), x, y ∈ R+ (5)
mit den Parametern c = 2θ2
θ23(1−e−θ2t) , q = 2θ1
θ23
− 1, u = cxe−θ2t, v = cy.
Iq(.) ist dabei die modifizierte Bessel-Funktion ersten Grades q, die sich als
Iq(x) =∞
∑
k=0
(x
2)2k+q 1
k!Γ(k + q + 1), x ∈ R (6)
aufschreiben lässt. Γ(.) bezeichnet die Gammafunktion Γ(z) =∫ ∞
0 xz−1e−xdx, z ∈ R.
2.4 Die Maximum-Likelihood-Schätzung
Seien X1, ..., Xn unabhängig identisch verteilte Zufallsvariablen mit gemeinsamer Dichte
f(x1, ..., xn|θ) = f(x1|θ) · · · f(xn|θ). θ ∈ Rk sei dabei ein Parametervektor, der eine be-
stimmte X1, ..., Xn zugrunde liegende Verteilung charakterisiert. Bei einer Realisierung
von Werten x1, ..., xn und unbekanntem Parametervektor θ kann die Dichte als eine Funk-
tion in Abhängigkeit von θ aufgefasst werden. Diese Funktion L(θ) = f(x1, ..., xn|θ) wird
als Likelihoodfunktion bezeichnet (Fahrmeir et al., 1999, S.372f). Sinnvoll ist es dann
als Schätzung für θ denjenigen Parameter θ̂ zu wählen, der die Likelihoodfunktion maxi-
miert, also L(θ̂) = maxθ f(x1, ..., xn|θ). Letztendlich liefert uns diese Funktion also für die
2 Statistische Modelle und Methoden 6
Zufallsvariablen X1, ..., Xn denjenigen Parameter θ̂, der am plausibelsten erscheint. Die
Schätzfunktion θ̂(x1, ..., xn), die uns in Abhängigkeit der Realisierungen eine Schätzung
für θ liefert, nennt man Maximum-Likelihood-Schätzer. Das Maximum dieser Funktion
kann, falls die Funktion differenzierbar ist, durch partielles Ableiten nach θ1, ..., θn und
Nullsetzen der Ableitungen berechnet werden. In vielen Fällen ist es sinnvoll die Like-
lihoodfunktion zunächst zu logarithmieren und diese dann zu maximieren, da die Terme
oft einfacher abzuleiten sind und das Maximum einer Funktion wegen der strikten Mo-
notonie des Logarithmus auch das Maximum seiner Logarithmierung ist. Diese Funktion
l(θ) = ln L(θ) =n∑
i=1ln f(xi|θ) wird Log-Likelihood-Funktion genannt.
2.5 Die ML-Schätzung bei stochastischen DGLs
Seien x1, ..., xn zeitdiskrete Beobachtungen eines stochastischen Prozesses, der über ei-
ne stochastische Differentialgleichung wie in (1) definiert ist. Sei außerdem für jede Be-
obachtung X1, ..., Xn die auf die vorhergehende Beobachtung bedingte Dichtefunktion
pθ(t, xi, |xi−1) mit unbekanntem θ bekannt. Dann kann der Parametervektor θ der sto-
chastischen Differentialgleichung mittels einer Maximum-Likelihood-Schätzung geschätzt
werden. Aufgrund der Markov-Eigenschaft von X hat nur die Beobachtung xi−1 einen
Einfluss auf xi und die Likelihoodfunktion ist dann gegeben durch
L(θ) =n
∏
i=1
pθ(t, xi, |x1, ..., xi−1)pθ(x0) =n
∏
i=1
pθ(t, xi, |xi−1)pθ(x0). (7)
Mittels Logarithmierung lässt sich daraus leicht die Log-Likelihoodfunktion l(θ) = ln(L(θ))
bilden. Problematisch ist hierbei der Wert pθ(x0), da die Dichte zu diesem Zeitpunkt oft
nicht gegeben ist. Für eine große Anzahl n an Beobachtungen wird der Einfluss von pθ(x0)
so gering, dass er vernachlässigt und auf 1 gesetzt werden kann. Um das Maximum zu
erhalten können wie bei der normalen Maximum-Likelihood-Schätzung die partiellen Ab-
leitungen nach θ1, ..., θn gleich Null gesetzt werden. Unter bestimmten Bedingungen gilt,
dass die Maximum-Likelihood-Schätzung konsistent und asymptotisch normalverteilt ist
(Iacus, 2008, S.112). Statt die Bedingungen zu überprüfen kann man dies jedoch auch
mittels Empirie feststellen. Bei bestimmten bedingten Dichten ist es nur schwer möglich
3 Simulationen 7
mittels analytischer Methoden ein Maximum zu finden. In diesen Fällen können numeri-
sche Verfahren wie das L−BFGS−B angewendet werden, das bereits im vorhergehenden
Vortrag beschrieben wurde.
3 Simulationen
In den folgenden Simulationen soll von Interesse sein, die Parameter der bedingten Dich-
ten der stochastischen Differentialgleichungen bei gegebenen Realisierungen möglichst gut
zu schätzen. Die Realisierungen eines bestimmten stochastischen Prozesse können mit der
Funktion sde.sim in R (R Development Core Team 2012) aus dem Paket sde (Iacus,
2009) simuliert werden. Dazu müssen der Drift, der Diffusionskoeffizient und ein Simulati-
onsschema oder die bedingte Dichtefunktion einer bestimmten stochastischen Differential-
gleichung angegeben werden. Außerdem sollten noch die Anzahl an Realisierungen n und
die Zeitabstände t zwischen zwei aufeinander folgenden Realisierungen angegeben wer-
den. Anhand der simulierten Daten soll daraufhin eine Maximum-Likelihood-Schätzung
(ML-Schätzung) durchgeführt werden, um die Parameter der stochastischen Differenti-
algleichungen zu schätzen. Als letztes soll evaluiert werden, wie gut die durchgeführten
Schätzungen sind. Unter bestimmten Bedingungen ist die Konsistenz und die asymptoti-
sche Normalverteilung der ML-Schätzung gegeben. Ob die Schätzungen dies erfüllen, soll
hier empirisch überprüft werden. Außerdem soll das Verhalten der Schätzer bei Variierung
der Zeitabstände t untersucht werden.
3.1 Black-Scholes-Merton-Modell
Zunächst wird das Black-Scholes-Merton-Modell betrachtet. Als erstes werden mit der
Funktion sde.sim jeweils n1 = 100, n2 = 1000 und n3 = 5000 Realisierungen einer geo-
metrischen Brownschen Bewegung mit Parametern θ1 = 0.5, θ2 = 0.2 und Zeitabständen
∆t = 0.01 erzeugt, der Anfangswert ist x0 = 1. Es werden dann zwei Funktionen zur
Berechnung des ML-Schätzers programmiert (vgl. Anhang). Es wird dazu die bedingte
Dichte dieses Modells benutzt, sie ist lognormalverteilt mit Dichte wie in (3) auf Seite 4.
3 Simulationen 8
Jede Beobachtung wird dabei auf die vorhergehende Beobachtung bedingt. Allerdings ist
das Maximum der Likelihoodfunktion von lognormalverteilten Zufallsvariablen analytisch
nicht einfach berechenbar, weswegen die numerische Methode “L-BFGS-B” verwendet
wird. Diese ist in der Funktion mle enthalten, die in R ohne Paketinstallation aufrufbar
ist. Da diese standardmäßig das Minimum einer Funktion sucht, wird das Vorzeichen der
Likelihood-Funktion umgedreht. Der Parametervektor θ̂ = (θ̂1, θ̂2) der möglichen Lösun-
gen wird auf den Raum (0.01, 0.01) ≤ θ̂ ≤ (1, 1) beschränkt. Außerdem wird mit Hilfe der
Hesse-Matrix die Standardabweichung sowie das 95%-Konfidenzintervall berechnet. Die
jeweils erhaltenen Parameterschätzer sind in Tabelle 1 zu sehen. Für eine größer werdende
Anzahl an Realisierungen liegen die Schätzungen immer näher am wahren θ = (0.5, 0.2).
Die Konfidenzintervalle überdecken bei θ1 immer den wahren Parameter, bei θ2 wird er
erst bei der letzten Schätzung mit n3 = 5000 Beobachtungen überdeckt.
Tabelle 1: ML-Schätzer, Standardabweichung und 95% Konfidenzintervalle für θ
Anzahl Realisierungen θ̂1 θ̂2
n1 = 100 ML-Schätzung 0.7311 0.2340Standardabweichung 0.2340 0.016595%-Konfidenzintervall [0.2690, 1.1953] [0.2050, 0.2706]
n2 = 1000 ML-Schätzung 0.4493 0.2097Standardabweichung 0.06632 0.004795%-Konfidenzintervall [0.3193, 0.5795] [0.2008, 0.2192]
n3 = 5000 ML-Schätzung 0.4704 0.2025Standardabweichung 0.0286 0.002095%-Konfidenzintervall [0.4143, 0.5266] [0.1986, 0.2066]
Um genauere Aussagen treffen zu können, werden für n = 100, 200, 300, ..., 5000 Beobach-
tungen einer brownschen Bewegung mit gleichen Einstellungen wie bei obigen Schätzun-
gen jeweils 1000 Prozesse simuliert. Für diese Prozesse wird jeweils eine ML-Schätzung
durchgeführt. Von diesen ML-Schätzung wird das arithmetische Mittel sowie das 0.025
bzw. 0.975-Quantil gebildet. Das Ergebnis für θ̂1 ist in Abbildung 1 zu sehen, für θ̂2 ist
eine analoge Grafik im Anhang zu sehen (Abbildung 3). Es ist erkennbar, dass bei θ1 und
θ2 die Mittelwerte der ML-Schätzungen immer nahe an den wahren θ1 und θ2 liegen und
die Abstände der 0.025 bzw. 0.975-Quantile immer enger werden. Es kann deshalb davon
3 Simulationen 9
ausgegangen werden, dass die Schätzungen im Mittel immer richtig sind und dass bessere
Schätzungen resultieren, wenn die Anzahl der Beobachtungen erhöht wird, da sich dann
der Abstand der 0.025-Quantilen zu den 0.975-Quantilen und damit auch die Varianz
verkleinert.
0 1000 2000 3000 4000 5000
0.2
0.4
0.6
0.8
Anzahl Realisierungen
Dur
chsc
hnitt
übe
r 10
00 W
iede
rhol
unge
n
ML−Schätzung0.025 und 0.975 Quantilewahrer Wert
Abbildung 1: Black-Scholes-Modell: Durchschnitt und 0.025 bzw. 0.975-Quantile der ML-Schätzungen von θ1 bei wachsender Anzahl an Realisierungen
Von Interesse ist außerdem, ob die asymptotische Normalverteilungsannahme richtig ist.
Dafür werden für eine große Anzahl an Beobachtungen n (n = 5000) 1000 Wiederholungen
einer ML-Schätzung für θ1 und θ2 durchgeführt, die restlichen Einstellungen bleiben gleich.
Diese Schätzungen werden dann mittels QQ-Plots, die in Abbildung 2 zu sehen sind, auf
Normalverteilung überprüft. Es ist hier keine markante Abweichung von einer Geraden
erkennbar, weshalb von asymptotischer Normalverteilung ausgegangen werden kann.
Darüber hinaus kann betrachtet werden, wie sich die ML-Schätzung verhält, wenn der
Zeitabstand t zwischen den Beobachtungen verändert wird. Für ∆t ≥ 0.3 ist die nume-
rische Berechnung des ML-Schätzers nicht mehr möglich, da die geometrische Bewegung
am Ende Werte erzeugt, die für R unendlich sind. Um das allgemeine Verhalten des ML-
Schätzers bei verschiedenen Zeitabständen ∆t zu betrachten, wird ∆t < 0.3 variiert, die
3 Simulationen 10
−3 −2 −1 0 1 2 3
0.40
0.45
0.50
0.55
Theoretische Quantile
Em
pris
che
Qua
ntile
−3 −2 −1 0 1 2 3
0.19
40.
198
0.20
2Theoretische Quantile
Em
pris
che
Qua
ntile
Abbildung 2: QQ-Plots der 1000 ML-Schätzungen für geometrische brownsche Bewegun-gen mit n = 5000 Realisierungen und θ1 (links) und θ2 (rechts)
Grenzen der möglichen Lösungen werden neu gesetzt ((−6, 0.01) ≤ θ ≤ (1, 1)). Es werden
für ∆t = e−11.4, e−11.2, e−11, ..., e−1.6 wieder jeweils 1000 Wiederholungen einer brownschen
Bewegung mit n = 1000 Beobachtungen durchgeführt und das arithmetische Mittel sowie
die 0.025 und 0.975-Quantile der ML-Schätzungen berechnet. Das Ergebnis ist in Abbil-
dung 4 im Anhang zu sehen. Die x-Achse ist logarithmiert, um eine bessere Übersicht
über die Werte zu bekommen. Für θ1 ist die Schätzung im Mittel immer richtig, mit grö-
ßer werdenden Zeitabständen ∆t wird die Schätzung immer besser, da der Abstand der
Quantile, die den wahren Wert 0.05 einschließen immer kleiner wird. Für θ2 dagegen sind
die Schätzungen unabhängig von ∆t überall gleich gut. Dies liegt vermutlich daran, dass
θ1 der Parameter für die deterministische Steigung ist. Diese wird bei größeren Zeitabstän-
den offensichtlicher, da die Schätzung nicht mehr so stark durch den zufälligen Verlauf der
Brownschen Bewegung beeinflusst wird. Der Parameter der brownschen Bewegung θ2 ist
im Gegensatz dazu unabhängig von den Zeitabständen gut berechenbar, da die Brownsche
Bewegung auf beliebig kleinen Räumen die gleiche Struktur aufweist.
3 Simulationen 11
3.2 Cox-Ingersoll-Ross-Modell
Als nächstes soll eine Schätzung der Parameter des Cox-Ingersoll-Ross-Modells (CIR-
Modell) durchgeführt werden. Dazu wird auch hier die bedingte Dichte verwendet um
eine Maximum Likelihood-Schätzung für θ durchzuführen. Problematisch ist hierbei die
in der bedingten Dichte enthaltene Bessel-Funktion Iq. Die in R implementierte Funkti-
on besselI(x,q) führt zu Overflows bei großem x und q. Overflow, auch arithmetischer
Überlauf genannt, bedeutet, dass das Ergebnis einer Rechnung für den gültigen Zahlen-
bereich (in R Absolutwerte von 2 · 10−308 bis 2 · 10308) zu groß ist, um noch richtig inter-
pretiert werden zu können. Ähnliche Probleme ergeben sich, wenn die χ2 -Verteilung zur
Likelihood-Schätzung benutzt wird, auch wenn hier seltener Overflows entstehen. Um bes-
sere Ergebnisse zu erzielen, wird die exponentiell skalierte Bessel-Funktion e−2√
uvIq(2√
uv)
verwendet, statt die normale in besselI(x,q) implementierte Funktion zu verwenden.
Diese kann mittels asymptotischer Expansion berechnet werden (siehe expBes im An-
hang). Die Log-Likelihood des CIR Modells ergibt sich zu
li(θ) = log c − (u + v) +q
2log
(
v
u
)
+ log Iq(2√
uv). (8)
Um die exponentiell skalierte Bessel-Funktion anzuwenden wird der Rechte Term von
log Iq(2√
uv) = log Iq(2√
uv) − 2√
uv + 2√
uv = log(e−2√
uvIq(2√
uv)) + 2√
uv (9)
benutzt, der sich aus dem letzten Term der Loglikelihood ergibt.
Trotz dieser Modifikation ist bei allen Alternativen die Berechnung des ML-Schätzers
aufgrund numerischer Overflows oftmals nicht möglich. Im weiteren Verlauf wird deshalb
sowohl die Variante mit der exponentiell skalierten Bessel-Funktion (CIR.lik) und die
Variante mit der χ2-Verteilung (CIR.lik2) benutzt (vgl. Anhang). Generell kommt es
wohl bei kleinen Beobachtungsanzahlen n und bei zu großen oder zu kleinen Zeitabständen
∆t vermehrt zu Fehlermeldungen. Zunächst werden 1000 CIR-Prozesse mit N = 1000
Realisierungen, Parametern θ1 = 0.2, θ2 = 0.06 und θ3 = 0.15 sowie Zeitabständen
∆t = 0.01 simuliert. Der Parametervektor θ̂ = (θ̂1, θ̂2, θ̂3) der möglichen Lösungen wird
3 Simulationen 12
auf den Raum (0.01, 0.01, 0.01) ≤ θ̂ ≤ (1, 1, 1) beschränkt. Bei CIR.lik konnte genau wie
bei CIR.lik2 132-mal keine Schätzung durchgeführt werden, da ein Fehler auftrat. Leider
traten die Fehler jeweils bei den selben CIR-Prozessen auf. Bei anderen Einstellungen
konnte es aber durchaus sein, dass eine Funktion ein Ergebnis lieferte und die andere nicht,
weswegen dennoch beide Möglichkeiten im Auge behalten werden sollten. Ein Beispiel
für eine dieser ML-Schätzungen mit den berechneten Werten, Standardabweichungen und
Konfidenzintervallen ist in Tabelle 2 zu finden. Beide Schätzungen sind hier in etwa gleich,
der Parameter θ1 wird stark überschätzt.
Tabelle 2: Maximum Likelihood-Schätzer für den Parameter θ
Funktion θ̂1 θ̂2 θ̂3
(CIR.lik) ML-Schätzung 0.2938 0.0991 0.1525Standardabweichung 0.0699 0.0801 0.003495%-Konfidenzintervall [NA,0.4311] [-0.0579, 0.2563] [0.1461, 0.1595]
(CIR.lik2) ML-Schätzung 0.2938 0.0991 0.1525Standardabweichung 0.0699 0.0801 0.003495%-Konfidenzintervall [NA,0.4311] [-0.0579, 0.2563] [0.1461, 0.1595]
Ein weiterer interessanter Aspekt ist die Dauer der Berechnungen. Einen Schätzung mit-
tels der χ2-Verteilung zu berechnen benötigt mehr als zehnmal so viel Zeit wie die Funk-
tion, welche die exponentiell skalierte Bessel-Funktion verwendet.
Eine saubere Analyse der programmierten ML-Schätzungen ist wegen der häufig auftre-
tenden Fehlermeldungen schwerlich möglich, da die Fehlermeldungen wohl bei einer spe-
ziellen Klasse von Daten auftreten. Dies kann anhand Abbildung 5 im Anhang nach-
vollzogen werden. Hier wurden mit der Funktion CIR.lik ML-Schätzungen für n =
100, 200, 300, ..., 5000 jeweils 1000 mal durchgeführt, die Parameter wurden wie oben ge-
wählt. Der Parametervektor wurde auf den Raum (−2, 0.001, 0.001) ≤ θ̂ ≤ (1, 5, 1) ein-
geschränkt. Die arithmetischen Mittel und die 0.975-Quantile über die ML-Schätzungen
von θ1 und θ2 bei kleinen n, die keine Fehlermeldung aufweisen (ca. 85%), liegen weit über
dem wahren Wert. Dies könnte daran liegen, dass die Berechnung vor allem für Daten,
die eine kleine ML-Schätzung ergeben würden, scheitert. θ3 scheint dagegen für alle n
erwartungstreu geschätzt werden zu können. Für größer werdende n werden die Schätzun-
4 Zusammenfassung 13
gen besser für alle θ, dies ist an der Konvergenz der Mittelwerte gegen die wahren Werte
sowie den kleineren Abständen zwischen den Quantilen zu sehen. Auch die Anzahl an
Fehlermeldungen reduziert sich, so tritt beispielsweise für n = 5000 keine Fehlermeldung
mehr auf.
Um wieder die asymptotische Normalverteilungsannahme zu überprüfen werden für die
Ergebnisse bei n = 5000 für alle Parameter QQ-Plots erstellt, die im Anhang in den
Abbildungen 6 und 7 zu finden sind. Da auch keine Fehlermeldungen auftreten, kön-
nen Verzerrungen aus diesem Grund ausgeschlossen werden. Bei θ1 und θ2 sind leichte
Krümmungen zu erkennen. Einzig die ML-Schätzung für θ3 scheint die asymptotische
Normalverteilungsannahme zu erfüllen.
Auf eine grafische Analyse des Verhaltens bei unterschiedlichen Zeitabständen ∆t wird,
wegen der häufig auftretenden Fehlermeldungen bei zu großen und zu kleinen ∆t, verzich-
tet. Diese machen eine Überprüfung der Erwartungstreue und Varianz der ML-Schätzer
bei unterschiedlichen ∆t unmöglich.
4 Zusammenfassung
In dieser Ausarbeitung geht es darum, Methoden zur Schätzung von Parametern des
Black-Scholes-Modells und des Cox-Ingersoll-Ross-Modells vorzustellen und zu bewerten.
Beide Modelle bieten Lösungen zu speziellen stochastischen Differentialgleichungen. Die
Dichten zu einer Realisierung, bedingt auf die vorhergehende Realisierung, bei bekann-
ten Zeitabständen sind in beiden Fällen bekannt, weswegen die Parameter jeweils mit der
Maximum-Likelihood-Methode geschätzt werden können. Da eine analytische Berechnung
des Maximums wegen der Unkenntnis der partiellen Ableitungen nicht möglich ist, wer-
den numerische Methoden verwendet. Dazu muss der Lösungsraum des Parametervektors
vorab beschränkt werden.
Beim Black-Scholes-Modell werden mit dieser Methode für eine ausreichend große Anzahl
an Beobachtungen recht gute Ergebnisse erzielt, die Schätzung für beide Parameter ist
durchgehend erwartungstreu. Zusätzlich wird die ML-Schätzung bei wachsender Anzahl
4 Zusammenfassung 14
an Beobachtungen immer genauer. Für eine große Anzahl an Beobachtungen kann auch
festgestellt werden, dass die ML-Schätzer normalverteilt um den wahren Wert θ streu-
en, dass also eine asymptotische Normalverteilung vorliegt. Bei konstant gehaltener An-
zahl an Beobachtungen und größer werdenden Zeitabständen wird die Schätzung für den
Steigungsparameter θ1 immer besser. Die Schätzung für den Parameter der brownschen
Bewegung θ2 bleibt dagegen konstant gut.
Die ML-Schätzung beim Cox-Ingersoll-Ross-Modell ist problematischer. Bei Benutzung
der normalen bedingten Dichte dieses Modells ergeben sich große Probleme bei der ML-
Schätzung, da es zu numerischen Overflows kommt. Deshalb werden modifizierte Versio-
nen dieser bedingten Dichte verwendet, die jedoch auch nicht immer einwandfrei funktio-
nieren und oft zu Fehler führen. Um Fehler zu vermeiden sollte die Beobachtungsanzahl
n ausreichend groß und der Zeitabstand in etwa um ∆t=0.01 gewählt werden. Die Schät-
zungen für die deterministischen Parameter θ1 und θ2 sind nicht erwartungstreu, aber
wenigstens konsistent. Auch die asymptotische Normalverteilungsannahme der Schätzer
ist bei den beiden Parametern nicht sicher gewährleistet. Die Schätzung für θ3 ist dagegen
erwartungstreu und asymptotisch normalverteilt. Auf eine Analyse bei unterschiedlichen
Zeitabständen ∆t wird verzichtet, da eine zu große Variation von ∆t nach oben oder nach
unten zu einer großen Anzahl an Fehlermeldungen führt, welche die Analyse verzerrt.
Weitergehende Untersuchungen könnten beim Cox-Ingersoll-Ross-Modell versuchen die
Gefahr von numerischen Overflows bei der Maximum-Likelihood-Schätzung einzudämmen.
Dazu existieren bereits einige Überlegungen, beispielsweise könnte die bedingte Dichte
mittels einer Poisson mixing-Gamma Charakterisierung approximiert werden (vgl. Iacus,
2008, S.121).
A Anhang 15
A Anhang
A.1 Zusätzliche Grafiken und Tabellen
0 1000 2000 3000 4000 5000
0.17
0.19
0.21
0.23
Anzahl Realisierungen
Dur
chsc
hnitt
übe
r 10
00 W
iede
rhol
unge
n
ML−Schätzung0.025 und 0.975 Quantilewahrer Wert
Abbildung 3: Black-Scholes-Modell:Durchschnitt und 0.025 bzw. 0.975-Quantile der ML-Schätzungen von θ2 bei wachsender Anzahl an Realisierungen
A Anhang 16
−10 −8 −6 −4 −2
−2
02
4
logarithmierte Zeitabstände
Dur
chsc
hnitt
übe
r 10
00 W
iede
rhol
unge
n
ML−Schätzung0.025 und 0.975 Quantilwahrer Wert
−10 −8 −6 −4 −2
0.19
00.
195
0.20
00.
205
0.21
0
logarithmierte Zeitabstände
Dur
chsc
hnitt
übe
r 10
00 W
iede
rhol
unge
n
ML−Schätzung0.025 und 0.975 Quantilwahrer Wert
Abbildung 4: Black-Scholes-Modell: Durchschnitt und 0.025 bzw. 0.975-Quantile der ML-Schätzungen von θ1 (oben) und θ2 (unten) bei wachsenden Zeitabständen
A Anhang 17
0 1000 2000 3000 4000 5000
0.2
0.3
0.4
0.5
Anzahl Realisierungen
Dur
chsc
hnitt
übe
r 10
00 W
iede
rhol
unge
n
ML−Schätzung0.025 und 0.975 Quantilwahrer Wert
0 1000 2000 3000 4000 5000
0.0
0.2
0.4
0.6
0.8
1.0
Anzahl Realisierungen
Dur
chsc
hnitt
übe
r 10
00 W
iede
rhol
unge
n
ML−Schätzung0.025 und 0.975 Quantilwahrer Wert
0 1000 2000 3000 4000 5000
0.13
50.
145
0.15
50.
165
Anzahl Realisierungen
Dur
chsc
hnitt
übe
r 10
00 W
iede
rhol
unge
n
ML−Schätzung0.025 und 0.975 Quantilwahrer Wert
Abbildung 5: Cox-Ingersoll-Ross-Modell: Durchschnitt und 0.025 bzw. 0.975-Quantile derML-Schätzungen von θ1 (oben), θ2 (mittig) und θ3 (unten) bei wachsendenZeitabständen
A Anhang 18
−3 −2 −1 0 1 2 3
0.15
0.25
0.35
Theoretische Quantile
Em
pris
che
Qua
ntile
−3 −2 −1 0 1 2 3
0.05
0.10
0.15
Theoretische Quantile
Em
pris
che
Qua
ntile
Abbildung 6: QQ-Plots der 1000 ML-Schätzungen des Cox-Ingersoll-Ross-Modells mitn = 5000 Realisierungen für θ1 (links) und θ2 (rechts)
−3 −2 −1 0 1 2 3
0.14
60.
150
0.15
4
Theoretische Quantile
Em
pris
che
Qua
ntile
Abbildung 7: QQ-Plot der 1000 ML-Schätzungen des Cox-Ingersoll-Ross-Modells mit n =5000 Realisierungen für θ3
A Anhang 19
A.2 R-Code
# Black-Scholes Model
dcBS <- function(x,t,x0,theta,log=TRUE){
m1 <- log(x0) + (theta[1]-theta[2]^2/2)*t
s1 <- sqrt(t)*theta[2]
lik <- dlnorm(x,meanlog=m1,sdlog=s1,log=TRUE)
if(!log)
lik <-exp(lik)
lik
}
BS.lik <- function(theta1,theta2){
n <- length(X)
dt <- deltat(X)
-sum(dcBS(x=X[2:n], t=dt, x0=X[1:(n-1)], theta=c(theta1,theta2)))
}
#install.packages("sde")
library(sde)
# ex3.03.R
set.seed(216)
X <- sde.sim(model="BS", theta=c(0.5,0.2), delta=0.01)
fit <- mle(BS.lik, start=list(theta1=1, theta2=1), method="L-BFGS-B", lower=c(0.01,0.01))
coef(fit)
summary(fit)
confint(fit)
length(X)*deltat(X)
set.seed(216)
X <- sde.sim(model="BS", theta=c(0.5,0.2),N=1000, delta=0.01)
fit <- mle(BS.lik,start=list(theta1=1, theta2=1), method="L-BFGS-B", lower=c(0.01,0.01))
coef(fit)
summary(fit)
confint(fit)
length(X)*deltat(X)
set.seed(216)
X <- sde.sim(model="BS", theta=c(0.5,0.2), N=5000, delta=0.01)
fit <- mle(BS.lik,start=list(theta1=1, theta2=1), method="L-BFGS-B", lower=c(0.01,0.01))
coef(fit)
summary(fit)
confint(fit)
length(X)*deltat(X)
#extra (KI’s, Normalverteilung und Erwartungstreue)
# Funktion
set.seed(216)
k <- seq(100,5000,100)
e1 <- e2 <- q11 <- q12 <- q21 <- q22 <- numeric(length(k))
for (j in 1:length(k)){
f1 <- f2 <- numeric(1000)
for (i in 1:1000){
X <- sde.sim(model="BS", theta=c(0.5,0.2), N=k[j], delta=0.01)
fit <- mle(BS.lik,start=list(theta1=1, theta2=1), method="L-BFGS-B", lower=c(0.01,0.01))
f1[i] <- coef(fit)[[1]]
f2[i] <- coef(fit)[[2]]}
# Erwartungswert
e1[j] <- mean(f1)
e2[j] <- mean(f2)
# Quantile
q11[j] <- quantile(f1,0.025)
q12[j] <- quantile(f1,0.975)
q21[j] <- quantile(f2,0.025)
q22[j] <- quantile(f2,0.975)
}
setwd("/media/DATA/Studium/6. Semester/Nichtlineare dynamische Systeme")
# save(list(e1,e2,q11,q12,q21,q22,f1,f2),file="Realisationsanzahl.RData")
load("Realisationsanzahl.RData")
A Anhang 20
# Erwartungswerte
# theta1
par(mfrow=c(1,1))
plot(k,a[[1]],type="l",ylim=range(c(a[[1]],a[[3]],a[[4]])),main="",xlab="Anzahl Realisierungen",
ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(k,a[[3]],col="blue")
lines(k,a[[4]],col="blue")
lines(k, rep(0.5,50),lty="dashed")
legend(1000,0.32,c("ML-Schätzung","0.025 und 0.975 Quantile","wahrer Wert"),col=c("red","blue",
"black"),lty=c(1,1,2))
#theta2
plot(k,a[[2]],type="l",ylim=range(c(a[[2]],a[[5]],a[[6]])),main="",xlab="Anzahl Realisierungen",
ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(k,a[[5]],col="blue")
lines(k,a[[6]],col="blue")
lines(k, rep(0.2,50),lty="dashed")
legend(1000,0.189,c("ML-Schätzung","0.025 und 0.975 Quantile","wahrer Wert"),col=c("red","blue",
"black"),lty=c(1,1,2))
# asymptotisch normalverteilte Daten? (laut Buch)
par(mfrow=c(1,2))
qqnorm(f1,main="",ylab="Emprische Quantile",xlab="Theoretische Quantile")
qqline(f1)
qqnorm(f2,main="",ylab="Emprische Quantile",xlab="Theoretische Quantile")
qqline(f2)
# Bandbreitenveränderung
# l muss kleiner als 0.3 sein
set.seed(216)
l <- exp(seq(-11.4,-1.6,0.2))
log(l)
e1 <- e2 <- q11 <- q12 <- q21 <- q22 <- numeric(length(l))
for (j in 1:length(l)){
f1 <- f2 <- numeric(1000)
for (i in 1:1000){
X <- sde.sim(model="BS", theta=c(0.5,0.2), N=1000, delta=l[j])
fit <- mle(BS.lik,start=list(theta1=1, theta2=1), method="L-BFGS-B", lower=c(-6,0.01))
f1[i] <- coef(fit)[[1]]
f2[i] <- coef(fit)[[2]]}
# Erwartungswert
e1[j] <- mean(f1)
e2[j] <- mean(f2)
# Quantile
q11[j] <- quantile(f1,0.025)
q12[j] <- quantile(f1,0.975)
q21[j] <- quantile(f2,0.025)
q22[j] <- quantile(f2,0.975)
}
setwd("/media/DATA/Studium/6. Semester/Nichtlineare dynamische Systeme")
# save(list(e1,e2,q11,q12,q21,q22,f1,f2),file="Bandbreiten.RData")
load("Bandbreiten.RData")
# Erwartungswerte
# theta1
par(mfrow=c(2,1))
plot(log(l),b[[1]],type="l",ylim=range(c(b[[1]],b[[3]],b[[4]])),main="",xlab="logarithmierte
Zeitabstände",ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(log(l),b[[3]],col="blue")
lines(log(l),b[[4]],col="blue")
lines(log(l), rep(0.5,50),lty="dashed")
legend(log(l[20]),4,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),
col=c("red","blue","black"),lty=c(1,1,2))
#theta2
plot(log(l),b[[2]],type="l",ylim=range(c(b[[2]],b[[5]],b[[6]])),main="",
xlab="logarithmierte Zeitabstände",ylab="Durchschnitt über 1000 Wiederholungen",col="red")
A Anhang 21
lines(log(l),b[[5]],col="blue")
lines(log(l),b[[6]],col="blue")
lines(log(l), rep(0.2,50),lty="dashed")
legend(log(l[20]),0.207,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),
col=c("red","blue","black"),lty=c(1,1,2))
############################ Cox-Ingersoll-Ross model ######################################
expBes <- function(x,nu){
mu <- 4*nu^2
A1 <- 1
A2 <- A1*(mu-1)/(1*(8*x))
A3 <- A2*(mu-9)/(2*(8*x))
A4 <- A3*(mu-25)/(3*(8*x))
A5 <- A4*(mu-49)/(4*(8*x))
A6 <- A5*(mu-81)/(5*(8*x))
A7 <- A6*(mu-121)/(6*(8*x))
1/sqrt(2*pi*x)*(A1-A2+A3-A4+A5-A6+A7)
}
dcCIR <- function(x, t, x0, theta, log=FALSE){
c <- 2*theta[2]/((1-exp(-theta[2]*t))*theta[3]^2)
ncp <- 2*c*x0*exp(-theta[2]*t)
df <- 4*theta[1]/theta[3]^2
u <- c*x0*exp(-theta[2]*t)
v <- c*x
q <- 2*theta[1]/theta[3]^2-1
lik <- (log(c)- (u+v)+q/2*log(v/u)+log(expBes(2*sqrt(u*v),q))+2*sqrt(u*v))
if(!log)
lik <- exp(lik)
lik
}
CIR.lik <- function(theta1, theta2, theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcCIR(x=X[2:n],t=dt,x0=X[1:(n-1)],theta=c(theta1,theta2,theta3),log=TRUE))
}
# inefficient version based on noncentral chi^2 density
dcCIR2 <- function(x,t,x0,theta,log=FALSE){
c <- 2*theta[2]/((1-exp(-theta[2]*t))*theta[3]^2)
ncp <- 2*c*x0*exp(-theta[2]*t)
df <- 4*theta[1]/theta[3]^2
lik <- (dchisq(2*x*c,df=df,ncp=ncp,log=TRUE)+log(2*c))
if(!log)
lik <- exp(lik)
lik
}
CIR.lik2 <- function(theta1,theta2,theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcCIR2(x=X[2:n], t=dt, x0=X[1:(n-1)], theta=c(theta1,theta2,theta3),log=TRUE))
}
# Simulationen CIR
# Fehlersuche
set.seed(216)
f1 <- f2 <- f3 <- f21 <- f22 <- f23 <- numeric(1000)
for (i in 1:1000){
X <- sde.sim(X0=0.1, model="CIR", theta=c(0.2,0.06,0.15),N=1000,delta=0.01)
try(fit <- mle(CIR.lik, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(0.001,0.001,0.001),upper=c(1,1,1)))
try(fit2 <- mle(CIR.lik2, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(0.001,0.001,0.001),upper=c(1,1,1)))
f1[i] <- coef(fit)[[1]]
f2[i] <- coef(fit)[[2]]
A Anhang 22
f3[i] <- coef(fit)[[3]]
f21[i] <- coef(fit2)[[1]]
f22[i] <- coef(fit2)[[2]]
f23[i] <- coef(fit2)[[3]]}
for (i in 2:1000){
if((is.na(f1[i])|is.na(f1[i-1]))==F){
if(f1[i]==f1[i-1]){
f1[i] <- NA
f2[i] <- NA
f3[i] <- NA
}}}
for (i in 2:1000){
if((is.na(f21[i])|is.na(f21[i-1]))==F){
if(f21[i]==f21[i-1]){
f21[i] <- NA
f22[i] <- NA
f23[i] <- NA
}}}
setwd("/media/DATA/Studium/6. Semester/Nichtlineare dynamische Systeme")
# c <- list(f1,f21)
# save(c,file="FehlerCIR.RData")
load("FehlerCIR.RData")
length(which(is.na(c[[1]])))
length(which(is.na(c[[2]])))
which(is.na(c[[1]]))==which(is.na(c[[2]]))
set.seed(216)
X <- sde.sim(X0=0.1, model="CIR", theta=c(0.2,0.06,0.15),N=1000,delta=0.01)
fit <- mle(CIR.lik, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(0.001,0.001,0.001),upper=c(1,1,1))
coef(fit)
summary(fit)
confint(fit)
length(X)*deltat(X)
set.seed(216)
fit2 <- mle(CIR.lik2, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(0.001,0.001,0.001),upper=c(1,1,1))
coef(fit2)
summary(fit2)
confint(fit2)
length(X)*deltat(X)
system.time(fit <- mle(CIR.lik, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(0.001,0.001,0.001),upper=c(1,1,1)))
system.time(fit <- mle(CIR.lik2, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(0.001,0.001,0.001),upper=c(1,1,1)))
#extra (KI’s, Normalverteilung und Erwartungstreue)
# Funktion
set.seed(216)
k <- seq(100,5000,100)
e1 <- e2 <- e3 <- q11 <- q12 <- q21 <- q22 <- q31 <- q32<- numeric(length(k))
for (j in 1:length(k)){
f1 <- f2 <- f3 <- numeric(1000)
for (i in 1:1000){
X <- sde.sim(X0=0.1, model="CIR", theta=c(0.2,0.06,0.15),N=k[j],delta=0.01)
try(fit <- mle(CIR.lik, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(-2,0.001,0.001),upper=c(1,5,1)))
f1[i] <- coef(fit)[[1]]
f2[i] <- coef(fit)[[2]]
f3[i] <- coef(fit)[[3]]}
for (i in 2:1000){
if((is.na(f1[i])|is.na(f1[i-1]))==F){
if(f1[i]==f1[i-1]){
f1[i] <- NA
f2[i] <- NA
f3[i] <- NA
A Anhang 23
}}}
# Erwartungswert
e1[j] <- mean(f1,na.rm=T)
e2[j] <- mean(f2,na.rm=T)
e3[j] <- mean(f3,na.rm=T)
# Quantile
q11[j] <- quantile(f1,0.025,na.rm=T)
q12[j] <- quantile(f1,0.975,na.rm=T)
q21[j] <- quantile(f2,0.025,na.rm=T)
q22[j] <- quantile(f2,0.975,na.rm=T)
q31[j] <- quantile(f3,0.025,na.rm=T)
q32[j] <- quantile(f3,0.975,na.rm=T)
}
setwd("/media/DATA/Studium/6. Semester/Nichtlineare dynamische Systeme")
d <- list(e1,e2,e3,q11,q12,q21,q22,q31,q32,f1,f2,f3)
# save(d,file="RealisationsanzahlCIR.RData")
load("RealisationsanzahlCIR.RData")
# Erwartungswerte
par(mfrow=c(3,1))
# theta1
plot(k,d[[1]],type="l",ylim=range(c(d[[1]],d[[4]],d[[5]])),main="",xlab="Anzahl Realisierungen",
ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(k,d[[4]],col="blue")
lines(k,d[[5]],col="blue")
lines(k, rep(0.2,50),lty="dashed")
legend(3000,0.53,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),
col=c("red","blue","black"),lty=c(1,1,2))
#theta2
plot(k,d[[2]],type="l",ylim=range(c(d[[2]],d[[6]],d[[7]])),main="",xlab="Anzahl Realisierungen",
ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(k,d[[6]],col="blue")
lines(k,d[[7]],col="blue")
lines(k, rep(0.06,50),lty="dashed")
legend(3000,1,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),col=c("red","blue","black"),
lty=c(1,1,2))
#theta3
plot(k,d[[3]],type="l",ylim=range(c(d[[3]],d[[8]],d[[9]])),main="",xlab="Anzahl Realisierungen",
ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(k,d[[8]],col="blue")
lines(k,d[[9]],col="blue")
lines(k, rep(0.15,50),lty="dashed")
legend(3000,0.17,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),
col=c("red","blue","black"),lty=c(1,1,2))
# asymptotisch normalverteilte Daten? (laut Buch)
par(mfrow=c(1,2))
qqnorm(d[[10]],main="",ylab="Emprische Quantile",xlab="Theoretische Quantile")
qqline(d[[10]])
qqnorm(d[[11]],main="",ylab="Emprische Quantile",xlab="Theoretische Quantile")
qqline(d[[11]])
par(mfrow=c(1,1))
qqnorm(d[[12]],main="",ylab="Emprische Quantile",xlab="Theoretische Quantile")
qqline(d[[12]])
# Bandbreitenveränderung
set.seed(216)
l <- exp(seq(-12,-2.2,0.2))
#l <- exp(seq(-10,-0.2,0.2))
e1 <- e2 <- e3 <- q11 <- q12 <- q21 <- q22 <- q31 <- q32<- numeric(length(l))
for (j in 1:length(l)){
f1 <- f2 <-f3 <- numeric(1000)
for (i in 1:1000){
X <- sde.sim(X0=0.1, model="CIR", theta=c(0.2,0.06,0.15),N=1000,delta=l[j])
try(fit <- mle(CIR.lik, start=list(theta1=0.1,theta2=0.1,theta3=0.3),method="L-BFGS-B",
lower=c(-2,0.001,0.001),upper=c(1,5,1)))
f1[i] <- coef(fit)[[1]]
A Anhang 24
f2[i] <- coef(fit)[[2]]
f3[i] <- coef(fit)[[3]]}
for (i in 2:1000){
if((is.na(f1[i])|is.na(f1[i-1]))==F){
if(f1[i]==f1[i-1]){
f1[i] <- NA
f2[i] <- NA
f3[i] <- NA
}}}
# Erwartungswert
e1[j] <- mean(f1,na.rm=T)
e2[j] <- mean(f2,na.rm=T)
e3[j] <- mean(f3,na.rm=T)
# Quantile
q11[j] <- quantile(f1,0.025,na.rm=T)
q12[j] <- quantile(f1,0.975,na.rm=T)
q21[j] <- quantile(f2,0.025,na.rm=T)
q22[j] <- quantile(f2,0.975,na.rm=T)
q31[j] <- quantile(f3,0.025,na.rm=T)
q32[j] <- quantile(f3,0.975,na.rm=T)
}
setwd("/media/DATA/Studium/6. Semester/Nichtlineare dynamische Systeme")
e <- list(e1,e2,e3,q11,q12,q21,q22,q31,q32,f1,f2,f3)
# save(e,file="BandbreitenCIR.RData")
load("BandbreitenCIR.RData")
# Erwartungswerte
# theta1
plot(log(l),e[[1]],type="l",ylim=range(c(e[[1]],e[[4]],e[[5]])),main="theta_1",xlab="Bandbreite",
ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(log(l),e[[4]],col="blue")
lines(log(l),e[[5]],col="blue")
legend(100,0.19,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),
col=c("red","blue","black"),lty=c(1,1,2))
#theta2
plot(log(l),e[[2]],type="l",ylim=range(c(e[[2]],e[[6]],e[[7]])),
main="theta_2",xlab="Bandbreite",ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(log(l),e[[6]],col="blue")
lines(log(l),e[[7]],col="blue")
legend(100,0.19,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),col=c("red","blue","black"),
lty=c(1,1,2))
#theta3
plot(log(l),e[[3]],type="l",ylim=range(c(e[[3]],e[[8]],e[[9]])),main="theta_3",xlab="Bandbreite",
ylab="Durchschnitt über 1000 Wiederholungen",col="red")
lines(log(l),e[[8]],col="blue")
lines(log(l),e[[9]],col="blue")
legend(100,0.19,c("ML-Schätzung","0.025 und 0.975 Quantil","wahrer Wert"),
col=c("red","blue","black"),lty=c(1,1,2))
A Anhang 25
Literatur
Ball, P. (2005): Critical Mass: How one thing leads to another, Arrow Books, London.
Fahrmeir, L., Künstler, R., Pigeot, I. und Tutz, G. (1999): Statistik - Der Weg zur Daten-
analyse, 2. Auflage, Springer, Berlin.
Iacus, S. M. (2008): Simulation and Inference for Stochastic Differential Equations, Sprin-
ger, New York.
Iacus, S. M. (2009): sde: Simulation and Inference for Stochastic Differential Equations,
R package version 2.0.10. URL: http://CRAN.R-project.org/package=sde
Kuo, H.-H. (2006): Introduction to Stochastic Integration, Springer, New York.
R Development Core Team (2012): A Language and Environment for Statistical Compu-
ting, R Foundation for Statistical Computing, Wien, URL: http://www.R-project.org/.
Top Related