____________________________________________________________________
Inhalt....: Empirische Untersuchungen zur Waldschen Identität
Kategorie.: Arbeitsblatt
Mathematik: Stochastik, Statistik
MuPAD.....: 3.1.1
Datum.....: 2005-07-04
Autoren...: Artur Foitzik <arturfoitzik@web.de>
Autoren...: Thorsten Fraas <thorstenfraas@web.de>
Autoren...: Gudrun Freitag <freitag@math.uni-goettingen.de>
Autoren...: Zeljko Golic <zeljkogo@aol.com>
Autoren...: Susanne Koch <skoch@math.uni-goettingen.de>
Autoren...: Katrin Kreh <katrin.kreh@stud.uni-goettingen.de>
Autoren...: Thorsten Mehlich <mehlicht@web.de>
Funktionen: stats::binomialRandom, stats::exponentialRandom, stats::normalPDF
Funktionen: stats::gammaRandom, stats::poissonRandom, stats::normalCDF
Funktionen: stats::mean, stats::variance
Funktionen: case, op, numeric::solve, stringlib::remove
Funktionen: plot::Histogram2d, plot::Point2d, plot::Function2d
____________________________________________________________________
Empirische Untersuchungen zur Waldschen
Identität
Die Waldsche Identität taucht bei der Behandlung bedingter Verteilungen und
bedingter Wahrscheinlichkeiten auf: Hier wird eine Summe S(N) von N
Zufallsvariablen, die identisch verteilt sind wie eine Zufallsgröße X mit Erwartungswert
E(X), betrachtet. Dabei ist N eine von den Summanden unabhängige Zufallsvariable mit
Werten in den natürlichen Zahlen. Die Waldsche Identität sagt aus, dass der
Erwartungswert von S(N) gerade das Produkt aus den Erwartungswerten E(X) und E(N)
ist, d.h. E(S(N))=E(N)*E(X). Hier ist die Unabhängigkeit von N und X entscheidend.
Sind auch die Summanden unabhängig voneinander, so lässt sich auch die Varianz von
S(N) als Funktion der ersten beiden Momente von X und N einfach darstellen (s.u.).
Im vorliegenden Notebook werden diese Identitäten für eine Beispielsituation aus der
Versicherungsmathematik empirisch überprüft. Hierbei werden unterschiedliche
Verteilungen zum einen für die Summanden und zum anderen für die Anzahl N
angenommen. Abschließend wird eine praktische Anwendung demonstriert.
Voraussetzungen:
(1) Zufallsvariablen und ihre Verteilung, insbesondere die Binomial-, die Poisson-, die
negative Binomial-, die Exponential- , die Normalverteilung und die Log-Normalverteilung,
(2) Momente (theoretisch und empirisch),
(3) Zentraler Grenzwertsatz,
(4) Histogramm.
Zunächst wollen wir die folgende Situation simulieren:
Bei einer Versicherung gehen jährlich N Schadensmeldungen ein; dabei ist N
eine N-wertige Zufallsvariable. Jeder der Schäden habe eine zufällige
Schadenshöhe, deren Betrag unabhängig von N ist.
Für die Simulation schreiben wir eine Prozedur Gesamtschaden. Hierbei nehmen
wir an, dass N eine binomialverteilte Zufallsvariable (mit Parametern n und p) ist
und die Schadenshöhen exponentialverteilte Zufallsvariablen (mit Parametern 0
und lambda) sind. Als Argumente übergeben wir der Prozedur Gesamtschaden
die Verteilungsparameter. Der Rückgabewert der Prozedur ist ein möglicher
Gesamtschaden eines Jahres. Zur Generierung der einzelnen Zufallsvariablen
werden die entsprechenden Prozeduren der stats-Bibliothek verwendet.
Zunächst schreiben wir aber eine Prozedur Euro, mittels der im Folgenden
alle Schadenshöhen auf zwei Nachkommastellen gerundet werden. Auch mit
den Schadenshöhen in Zusammenhang stehende Mittelwerte werden im
gesamten Notebook immer auf zwei Nachkommastellen gerundet.
Euro:=proc(x)
begin
return(float(round(x*100)/100))
end_proc:
Für eine einheitliche Ausgabe von Euro-Beträgen am Bildschirm (ggf. mit
'aufgefüllten Nullen') verwenden wir die Prozedur Eupr.
Eupr:=proc(x)
begin
if round(100*(float(trunc(10*Euro(x)))-10*Euro(x)))=0 then
return(expr2text(Euro(x))."0")
else return(expr2text(Euro(x)))
end_if
end_proc:
Gesamtschaden:=proc(n,p,lambda)
local i;
begin
N:=stats::binomialRandom(n,p)(); /*Schadensanzahl*/
X:=[Euro(stats::exponentialRandom(0,lambda)()) $ i=1..N];
/*Liste der einzelnen Schadenshöhen*/
return(_plus(X[i] $ i=1..N)) /*Gesamtschaden*/
end_proc:
Exemplarisch simulieren wir nun die Gesamtschadenshöhen für m=10
aufeinander folgende Jahre unter Annahme der Parameter n=1000, p=1/100,
lambda=1/50. Dazu berechnen wir den Mittelwert (im Folgenden als "empirischer
Erwartungswert" bezeichnet) und die empirische Varianz.
(Den Parameter m muss man nicht als Zeitdauer in Jahren betrachten; als
alternative Interpretationen bieten sich Zeitdauer in Monaten, Zeitdauer in
Wochen oder Anzahl von Versicherungsgesellschaften an.)
m:=10:
l:=[Gesamtschaden(1000,1/100,1/50) $ i=1..m]:
/*Liste von m Gesamtschäden*/
mitt:=stats::mean(l): var:=stats::variance(l):
print(Unquoted, "Gesamtschadenswerte: ". expr2text(op(l))):
print(Unquoted, "empirischer Erwartungswert: ". Eupr(mitt)):
print(Unquoted, "empirische Varianz: ". Eupr(var)):
Gesamtschadenswerte: 825.78, 67.92, 607.62, 349.42, 778.13, 691.79, 281.36, 614\
.25, 410.72, 251.53
empirischer Erwartungswert: 487.85
empirische Varianz: 63367.87
In diesem Beispiel erhielten wir Werte, die stark um ihren Mittelwert streuen,
was sich in einer sehr großen empirischen Varianz niederschlägt. Wie kommt
das zustande?
Zunächst betrachten wir die ersten beiden Momente der zugrunde liegenden
Verteilungen. Mit den oben angegebenen Parametern (n=1000, p=1/100,
lambda=1/50) erhalten wir für die theoretischen Erwartungswerte und die
theoretischen Varianzen:
E(N)=n*p=10, Var(N)=n*p*(1-p)=9.9,
E(X)=1/lambda=50, Var(x)=(1/lambda)^2=2500.
Für den theoretischen Erwartungswert von S(N) erhalten wir mit der
Waldschen Identität:
E(S(N))=E(N)*E(X)
=10*50=500.
Die Formel für die theoretische Varianz von S(N) findet sich z.B. in
Schmidt (2002, Satz 7.1.4) und in Klugman et al. (2004, Kapitel 6):
Var(S(N))=E(N)*Var(X)+Var(N)*(E(X))^2.
In unserem Beispiel erhalten wir damit
Var(S(N))=10*50^2+9.9*50^2=49750.
Wir sehen, dass diese theoretischen Werte ungefähr den oben erhaltenen
empirischen Werten entsprechen.
Im Folgenden wollen wir verschiedene Konstellationen von Verteilungen
von N und X untersuchen. Dazu schreiben wir eine Prozedur SN, die als
Argumente die Anzahl m der betrachteten Jahre (bzw. Monate / Wochen oder
Anzahl der Versicherungsgesellschaften) und für N und X jeweils eine Liste
(vertN bzw. vertX) mit einer Abkürzung (Zeichenkette) für den Verteilungstyp
und mit den zugehörigen Parameterwerten erhält.
Für N stehen die Binomialverteilung (Abkürzung "bi", Parameter n und p),
die Poissonverteilung (Abkürzung "po", Parameter lambda) und die negative
Binomialverteilung (Abkürzung "nb", Parameter r und beta; vgl. z.B. Gentle,
1998, Abschnitt 3.1.10) zur Auswahl.
Für X steht neben der Binomial- und der Poissonverteilung auch die
Exponentialverteilung zur Verfügung, mit Abkürzung "ex" und Parametern
theta (das ist der Verschiebungsparameter) und lambda (das ist der
Skalierungsparameter; es gilt lambda = 1/E(X) - theta).
Die Prozedur SN simuliert wie die obige Prozedur Gesamtschaden für m
Jahre die jeweilige Schadensanzahl N und die zugehörigen N Schadenshöhen
mit Verteilung von X. Dazu werden die einzelnen Listenelemente aus vertN
bzw. vertX verwendet. Die case-Anweisung wird zur Unterscheidung der
einzelnen Konstellationen herangezogen. Der Aufruf op(vertN,1) liefert das
erste Element der Liste vertN, in diesem Fall also den Verteilungstyp der
Schadensanzahl N (d.h. "bi", "po" oder "nb"). Analog spricht man mit op(vertN,2)
bzw. op(vertN,3) die zugehörigen Parameterwerte an.
Die Ausgabe der Prozedur SN umfasst Informationen zur zugrunde liegenden
Situation, sowie den ermittelten empirischen Erwartungswert (Mittelwert) und
die empirischeVarianz (über die m Jahre) der simulierten Werte von S(N) und
die entsprechenden theoretischen Momente aus den Waldschen Identitäten.
Der Rückgabewert der Prozedur ist eine Folge mit fünf Gliedern, bestehend aus
der Liste der m Gesamtschäden und den ersten beiden empirischen und
theoretischen Momenten.
SN:=proc(m, vertN , vertX)
local N, X, sn, k, i, EN, varN, EX, varX, mitt, var, wald1, wald2;
begin
sn:=[0 $ m]; print();
print(NoNL, "Anzahl m: " .expr2text(m). ", ");
for k from 1 to m do
case op(vertN,1)
of "bi" do
N:=stats::binomialRandom(op(vertN,2),op(vertN,3))();
EN:=op(vertN,2)*op(vertN,3);
varN:=op(vertN,2)*op(vertN,3)*(1-op(vertN,3));
if k=1 then
print(NoNL,
"Verteilung v. N: Bi(".expr2text(op(vertN,2)).",".expr2text(op(vertN,3))."), ");
end_if;
break;
of "po" do
N:=stats::poissonRandom(op(vertN,2))();
EN:=op(vertN,2);
varN:=EN;
if k=1 then
print(NoNL, "Verteilung v. N: Po(".expr2text(op(vertN,2))."), ");
end_if;
break;
of "nb" do
L:=stats::gammaRandom(op(vertN,2),op(vertN,3))();
N:=stats::poissonRandom(L)();
EN:=op(vertN,2)*op(vertN,3);
varN:=EN*(1+op(vertN,3));
if k=1 then
print(NoNL, "Verteilung v. N: Nb(".expr2text(op(vertN,2)).",".expr2text(op(vertN,3))."), ");
end_if;
break;
otherwise print(Unquoted, "Falsche Verteilungseingabe für N. ");return();
end_case;
case op(vertX,1)
of "bi" do
X:=[stats::binomialRandom(op(vertX,2),op(vertX,3))() $ i=1..N];
EX:=op(vertX,2)*op(vertX,3);
varX:=op(vertX,2)*op(vertX,3)*(1-op(vertX,3));
if k=1 then
print(NoNL, "Verteilung v. X: Bi(".expr2text(op(vertX,2)).",".expr2text(op(vertX,3)).")");
end_if;
break;
of "po" do
X:=[stats::poissonRandom(op(vertX,2))() $ i=1..N];
EX:=op(vertX,2); varX:=EX;
if k=1 then
print(NoNL, "Verteilung v. X: Po(".expr2text(op(vertX,2)).")");
end_if;
break;
of "ex" do
X:=[Euro(stats::exponentialRandom(op(vertX,2),(op(vertX,3)))()) $ i=1..N];
EX:=op(vertX,2)+1/op(vertX,3);
varX:=op(vertX,3)^(-2);
if k=1 then
print(NoNL, "Verteilung v. X: Ex(".expr2text(op(vertX,2)).",".expr2text(op(vertX,3)).")");
end_if;
break;
otherwise print(Unquoted, "Falsche Verteilungseingabe für X!");return();
end_case;
sn[k]:=_plus(X[i] $ i=1..N); /*Das ist der Gesamtschaden im k-ten Jahr*/
end_for;
mitt:=Euro(stats::mean(sn)); var:=Euro(stats::variance(sn)); /*empirische Momente von S(N)*/
wald1:=Euro(EX*EN); wald2:=Euro(EN*varX+varN*EX^2); /* theoretische Momente von S(N)*/
print();
print(NoNL, "theor. Erwartungswert v. S(N): ". Eupr(wald1). " theor. Varianz von S(N): ". Eupr(wald2));
print();
print(NoNL, "empir. Erwartungswert v. S(N): ". Eupr(mitt). " empir. Varianz von S(N): ". Eupr(var));
print();
return(sn,wald1,wald2,mitt,var):
end_proc:
Exemplarisch verwenden wir die Prozedur SN für die Situation von m=100
Jahren, binomialverteilten Schadensanzahlen mit wachsendem n und festem
p=0.05, sowie poissonverteilten Schadenshöhen mit erwarteter Schadenshöhe
lambda=100.
for n from 100 to 500 step 100 do
SN(100,["bi",n,0.05],["po",100]):
end_for:
Anzahl m: 100, Verteilung v. N: Bi(100,0.05), Verteilung v. X: Po(100)theor\
. Erwartungswert v. S(N): 500.00 theor. Varianz von S(N): 48000.00empir. Erwa\
rtungswert v. S(N): 441.06 empir. Varianz von S(N): 40788.38Anzahl m: 100, \
Verteilung v. N: Bi(200,0.05), Verteilung v. X: Po(100)theor. Erwartungswert \
v. S(N): 1000.00 theor. Varianz von S(N): 96000.00empir. Erwartungswert v. S(\
N): 949.41 empir. Varianz von S(N): 90725.15Anzahl m: 100, Verteilung v. N:\
Bi(300,0.05), Verteilung v. X: Po(100)theor. Erwartungswert v. S(N): 1500.00\
theor. Varianz von S(N): 144000.00empir. Erwartungswert v. S(N): 1506.41 e\
mpir. Varianz von S(N): 126003.90Anzahl m: 100, Verteilung v. N: Bi(400,0.05)\
, Verteilung v. X: Po(100)theor. Erwartungswert v. S(N): 2000.00 theor. Var\
ianz von S(N): 192000.00empir. Erwartungswert v. S(N): 1954.81 empir. Varianz\
von S(N): 194908.62Anzahl m: 100, Verteilung v. N: Bi(500,0.05), Verteilun\
g v. X: Po(100)theor. Erwartungswert v. S(N): 2500.00 theor. Varianz von S(N)\
: 240000.00empir. Erwartungswert v. S(N): 2468.54 empir. Varianz von S(N): 24\
7257.73
Wieder entsprechen die empirischen Werte der ersten beiden Momente
ungefähr den theoretischen Werten. Wir veranschaulichen nun die empirische
Verteilung der Gesamtschäden von m Jahren mit Hilfe eines Histogramms.
Dazu schreiben wir eine Prozedur Hist, welche als Argumente genau die
gleichen wie die Prozedur SN übergeben bekommt. Die Ausgabe der
Prozedur Hist besteht zum einen aus den Informationen, die die Prozedur SN
zu den entsprechenden Argumenten ausgibt und zum anderen aus einem
Histogramm der Gesamtschäden, wobei der theoretische Erwartungswert
E(S(N)) des Gesamtschadens durch einen grünen Punkt markiert wird.
Hist:=proc(m, vertN , vertX)
local sn,h,p,simul,wald1,wald2,mitt,var;
begin
simul:=SN(m,vertN,vertX):
sn:=op(simul,1);
wald1:=op(simul,2);
wald2:=op(simul,3);
mitt:=op(simul,4);
var:=op(simul,5);
h:=plot::Histogram2d(sn, Area=1, Cells=round(sqrt(m)),
XTicksLabelStyle=Vertical,XAxisTitle="",YAxisTitle="");
p:=plot::Point2d([wald1,0],PointSize=3, Color=RGB::SeaGreen);
plot(h,p,Header="Relative Häufigkeit der Gesamtschäden");
end_proc:
Im Folgenden lassen wir uns das Histogramm für einen Zeitraum von m=1000
Jahren in der Situation poissonverteilter Schadensanzahl (hier ist der Parameter
zunächst 400) und exponentialverteilter Schadenshöhe (mit Parametern 0 und
1/5) zeichnen.
Hist(1000,["po",400],["ex",0,1/5])
Anzahl m: 1000, Verteilung v. N: Po(400), Verteilung v. X: Ex(0,1/5)\
theor. Erwartungswert v. S(N): 2000.00 theor. Varianz von S(N): 20000.00empir\
. Erwartungswert v. S(N): 2000.54 empir. Varianz von S(N): 21888.63
Das Histogramm suggeriert, dass S(N) annähernd normalverteilt ist. Dies ist
nicht sofort einzusehen, da S(N) eine Summe von zufällig vielen, unabhängigen
und identisch verteilten Zufallsvariablen ist. Der klassische zentrale Grenzwertsatz
(vgl. Dehling & Haupt, 2003, Kapitel 10) sagt nur aus, dass die Verteilung einer
Summe von n (fest!) unabhängigen, identisch verteilten Zufallsvariablen für große
n gegen eine Normalverteilung strebt.
Dennoch ist unsere Beobachtung nicht völlig überraschend: Denn im gewählten
Fall der Poissonverteilung mit Parameter 400 für N ist die durchschnittliche Anzahl
von Schadensfällen pro Jahr, also E(N), groß: In diesem Fall gilt nämlich gerade
E(N)=400. Und für große Werte von E(N) ist auch S(N) näherungsweise
normalverteilt (hierzu sei auf Klugman et al., 2004, Kapitel 6, verwiesen).
Wir erweitern nun die Prozedur Hist um den Plot der Normalverteilungsdichte mit
Erwartungswert E(S(N)) und Varianz Var(S(N)). Dies geschieht in der Prozedur
Hist2:
Hist2:=proc(m, vertN , vertX)
local sn,h,p,simul,wald1,wald2,mitt,var,norm;
begin
simul:=SN(m,vertN,vertX);
sn:=op(simul,1); wald1:=op(simul,2); wald2:=op(simul,3); mitt:=op(simul,4); var:=op(simul,5);
h:=plot::Histogram2d(sn, Area=1, Cells=round(sqrt(m)),
XTicksLabelStyle=Vertical,XAxisTitle="",YAxisTitle="");
p:=plot::Point2d([wald1,0],PointSize=3, Color=RGB::SeaGreen);
norm:=plot::Function2d(stats::normalPDF(wald1,wald2)(x),
x=mitt-3.5*sqrt(var)..mitt+3.5*sqrt(var), Color=RGB::Blue,LineWidth=0.5);
plot(h,p,norm,Header="Relative Häufigkeit der Gesamtschäden");
end_proc:
In der gleichen Situation wie eben ergibt sich dann folgendes Bild:
Hist2(1000,["po",400],["ex",0,1/5])
Anzahl m: \
1000, Verteilung v. N: Po(400), Verteilung v. X: Ex(0,1/5)theor. Erwartungs\
wert v. S(N): 2000.00 theor. Varianz von S(N): 20000.00empir. Erwartungswert \
v. S(N): 1994.78 empir. Varianz von S(N): 19852.79
Da die eben gewählte Exponentialverteilung für die Schadenshöhen rechtsschief
(das heißt auch linkssteil) ist, funktioniert der zentrale Grenzwertsatz für kleinere
E(N) noch nicht so gut (die Verteilung von S(N) ist dann ebenfalls rechtsschief):
Dies können wir gut erkennen, wenn wir die mittlere Schadensanzahl durch einen
kleineren Parameter bei der Poissonverteilung (10 statt 400) reduzieren:
Hist2(1000,["po",10],["ex",0,1/5])
Anzahl m: 1000, Verteilun\
g v. N: Po(10), Verteilung v. X: Ex(0,1/5)theor. Erwartungswert v. S(N): 50.0\
0 theor. Varianz von S(N): 500.00empir. Erwartungswert v. S(N): 50.08 empir\
. Varianz von S(N): 484.48
Für ungefähr symmetrisch verteilte Schadenshöhen wird der zentrale
Grenzwertsatz auch schon für kleinere N 'greifen', z.B. in der folgenden
Situation:
Hist2(1000,["po",10],["bi",20,0.45])
Anzahl m: 1000, Verteilung v. N: Po(10), Verteilu\
ng v. X: Bi(20,0.45)theor. Erwartungswert v. S(N): 90.00 theor. Varianz von S\
(N): 859.50empir. Erwartungswert v. S(N): 88.38 empir. Varianz von S(N): 853.\
75
Wir wollen uns nun einer praktischen Anwendung widmen. In Klugman et
al. (2004, Examples 6.4, 6.5), wird folgende Situation vorgestellt:
Empirische Mittelwerte und Standardabweichungen der Schadensanzahlen
bzw. Schadenshöhen aus vergangenen 10 Monaten sind gegeben.
Zunächst (vgl. Example 6.4) sollen daraus mit Hilfe der Waldschen Formeln
die Werte E(S(N)) und Var(S(N)) geschätzt werden (hierbei ist S(N) nun der
monatliche Gesamtschaden).
Die folgende Prozedur EMPIR1 dient dazu, diese Berechnungen für
simulierte Daten (wie oben in der Prozedur SN) für eine Anzahl m von
Monaten durchzuführen. Als Argumente enthält die Prozedur EMPIR1 also
eine Anzahl m von Monaten, den Verteilungstyp der Schadensanzahl (mit
Parametern) und den Verteilungstyp der Schadenshöhen (mit Parametern).
Ausgegeben werden dann die beiden Verteilungen, die theoretischen Werte
der ersten beiden Momente von S(N) und die geschätzten Werte dieser beiden
Größen. Letztere gehen aus den simulierten Daten durch Anwendung der
Waldschen Identitäten hervor. Der Rückgabewert der Prozedur ist die Folge
der geschätzten Werte für die ersten beiden Momente.
EMPIR1:=proc(m, vertN , vertX)
local N, X, k, i, l, Nlist, sumN, sumN2, mittN, varN, empvarN,
sumX, sumX2, mittX, varX, empvarX, empwald1, empwald2;
begin
Nlist:=[0 $ m];
/*Initialisierung der Liste der m Schadensanzahlen*/
sumX:=0;
sumX2:=0;
print();
print(NoNL, "Anzahl m: " .expr2text(m). ", ");
for k from 1 to m do
case op(vertN,1)
of "bi" do
N:=stats::binomialRandom(op(vertN,2),op(vertN,3))();
EN:=op(vertN,2)*op(vertN,3);
varN:=op(vertN,2)*op(vertN,3)*(1-op(vertN,3));
if k=1 then
print();
print(NoNL, "Verteilung von N: Bi(".expr2text(op(vertN,2)).",".expr2text(op(vertN,3))."), ");
end_if;
break;
of "po" do N:=stats::poissonRandom(op(vertN,2))();
EN:=op(vertN,2);
varN:=EN;
if k=1 then
print(NoNL, "Verteilung von N: Po(".expr2text(op(vertN,2))."), ");
end_if;
break;
of "nb" do
L:=stats::gammaRandom(op(vertN,2),op(vertN,3))();
N:=stats::poissonRandom(L)();
EN:=op(vertN,2)*op(vertN,3);
varN:=EN*(1+op(vertN,3));
if k=1 then
print(NoNL, "Verteilung von N: Nb(".expr2text(op(vertN,2)).",".expr2text(op(vertN,3))."), ");
end_if;
break;
otherwise print(Unquoted, "Falsche Verteilungseingabe für N.");return();
end_case;
case op(vertX,1)
of "bi" do
X:=[stats::binomialRandom(op(vertX,2),op(vertX,3))() $ i=1..N];
EX:=op(vertX,2)*op(vertX,3);
varX:=op(vertX,2)*op(vertX,3)*(1-op(vertX,3));
if k=1 then
print(NoNL, "Verteilung von X: Bi(".expr2text(op(vertX,2)).",".expr2text(op(vertX,3)).")");
end_if;
break;
of "po" do
X:=[stats::poissonRandom(op(vertX,2))() $ i=1..N];
EX:=op(vertX,2);
varX:=EX;
if k=1 then
print(NoNL, "Verteilung von X: Po(".expr2text(op(vertX,2)).")");
end_if;
break;
of "ex" do
X:=[Euro(stats::exponentialRandom(op(vertX,2),op(vertX,3))()) $ i=1..N];
EX:=op(vertX,2)+1/op(vertX,3);
varX:=op(vertX,3)^(-2);
if k=1 then
print(NoNL, "Verteilung von X: Ex(".expr2text(op(vertX,2)).",".expr2text(op(vertX,3)).")");
end_if;
break;
otherwise print(Unquoted, "Falsche Verteilungseingabe für X!");return();
end_case;
Nlist[k]:=N;
sumX:=sumX+_plus(X[i] $ i=1..N);
sumX2:=sumX2+_plus((X[i])^2 $ i=1..N);
end_for:
wald1:=Euro(EX*EN); wald2:=Euro(EN*varX+varN*EX^2); /*theoretische Momente von S(N)*/
print();
print(NoNL, " E(S(N)) theoretisch: ". Eupr(wald1). " Var(S(N)) theoretisch: ". Eupr(wald2));
sumN:=_plus(Nlist[i] $ i=1..m); sumN2:=_plus((Nlist[i])^2 $ i=1..m);
mittN:=sumN/m; empvarN:=(sumN2-m*(mittN)^2)/(m-1); /*empirische Momente von N*/
mittX:=sumX/sumN; empvarX:=(sumX2-sumN*(mittX)^2)/(sumN-1);/*empirische Momente von X*/
empwald1:=Euro(mittX*mittN);
empwald2:=Euro(mittN*empvarX+empvarN*mittX^2);/*geschätzte und gerundete Momente von S(N)*/
print();
print(NoNL, " E(S(N)) geschätzt: ". Eupr(empwald1). " Var(S(N)) geschätzt: ". Eupr(empwald2));
print();
return(empwald1, empwald2):
end_proc:
EMPIR1(10, ["po", 100], ["ex",0,1.5]):
Anzahl m: 10, Verteilung von N: Po(100), Verteilung von X: Ex(0,1.5) E(S(\
N)) theoretisch: 66.67 Var(S(N)) theoretisch: 88.89 E(S(N)) geschätzt: 6\
9.97 Var(S(N)) geschätzt: 98.81
Nun soll, wie in Example 6.5 in Klugman et al. (2004), die Wahrscheinlichkeit
dafür berechnet werden, dass S(N) mehr als das 1.4-fache des geschätzten
Wertes E(S(N)) beträgt.
Dazu wird zum einen die Normalverteilungsapproximation für S(N) verwendet.
Zum anderen wird in der angegebenen Literatur (ohne theoretische Grundlage)
vorgeschlagen, für kleine Werte von E(N) die Log-Normalverteilung zur
Approximation zu nutzen, da die Verteilung von S(N) dann üblicherweise schief
ist.
Die folgende Prozedur EMPIR2 ermöglicht diese Berechnungen für simulierte
Daten aus m Monaten. Der Parameter f gibt an, um welchen Faktor der
Erwartungswert überschritten wird (im obigen Beispiel ist f=1.4). Dieser wird,
zusammen mit der Anzahl m von Monaten und den Verteilungen für N und X,
der Prozedur als Argument übergeben. Ausgegeben werden dann zum einen
die bereits in der Prozedur EMPIR1 ausgegebenen Informationen und zum
anderen die Wahrscheinlichkeiten, dass S(N) den f-fachen Wert des geschätzten
Wertes von E(S(N)) überschreitet; einmal unter der Annahme, dass S(N)
normalverteilt ist (die Parameter werden dabei gerade durch die geschätzten
Werte der ersten beiden Momente gegeben) und einmal unter der Annahmne,
dass S(N) log-normalverteilt ist. Der Rückgabewert der Prozedur ist die Folge der
Überschreitungswahrscheinlichkeiten (einmal unter Normalverteilungs- und einmal
unter Log-Normalverteilungsannahme).
Es sei angemerkt, dass eine Zufallsvariable Y log-normalverteilt zu den
Parametern mu und sigma heißt, wenn log(U) normalverteilt zu den Parametern
mu und sigma ist. Dann gilt
E(Y)=exp(mu+0.5*sigma^2), E(Y^2)=exp(2*mu+2*sigma^2) . (*)
EMPIR2:=proc(m, vertN , vertX, f)
local estim, erwS, varS, fh1, fh2, h1, h2, h3, mu, sig2,
p1, pn1, pn2, p2, pln1, pln2;
begin
estim:=EMPIR1(m, vertN , vertX);
erwS:=op(estim,1); varS:=op(estim,2); /* geschätzte, gerundete Momente von S(N)*/
fh1:=floor(f); fh2:=stringlib::remove(expr2text(frac(f)),"0.");
print(NoNL, "Wahrscheinlichkeit f. Überschreitung d. ".expr2text(fh1).
".".fh2."-fachen Erwartungswertes v. S(N):");
print();
h1:=(erwS*f-erwS)/sqrt(varS); /*f-facher Erwartungswert, standardisiert*/
p1:=1-stats::normalCDF(0,1)(h1); /*Überschreitungswahrscheinlichkeit, mit NV-Annahme*/
pn1:=floor(p1); pn2:=stringlib::remove(expr2text(frac(p1)),"0.");
DIGITS:=10;
print(NoNL, " - m. Normalverteilungs-Approximation: ".expr2text(pn1).".".pn2); print();
DIGITS:=100;
h2:=[x+0.5*y=log(E,erwS),2*x+2*y=log(E,varS+erwS^2)]; /* Ausnutzen von (*) */
h3:=numeric::solve(h2,[x,y=0.0001..infinity]);
if h3={} then
print(NoNL, "- mit Log-Normalverteil.-Approximation: keine Lösung");
DIGITS:=10; break;
end_if;
mu:=text2expr(stringlib::remove(expr2text(op(op(h3,1),1)),"x = "));
sigma:=sqrt(text2expr(stringlib::remove(expr2text(op(op(h3,1),2)),"y = ")));
h4:=(log(E,erwS*f)-mu)/sigma; /*f-facher Erwartungswert, standardisiert*/
DIGITS:=10;
p2:=1-stats::normalCDF(0,1)(h4); /*Überschreitungswahrscheinlichkeit, mit Log-NV*/
pln1:=floor(p2); pln2:=stringlib::remove(expr2text(frac(p2)),"0.");
print(NoNL, " - mit Log-Normalvert.-Approximation: ".expr2text(pln1).".".pln2);
print();
return(p1,p2):
end_proc:
EMPIR2(12, ["nb", 4000, 0.9], ["bi",40,0.5], 1.1):
Anzahl m: 12, Verteilung von N: Nb(4000,\
0.9), Verteilung von X: Bi(40,0.5) E(S(N)) theoretisch: 72000.00 Var(S(N)\
) theoretisch: 2772000.00 E(S(N)) geschätzt: 72341.83 Var(S(N)) geschätz\
t: 2389358.31Wahrscheinlichkeit f. Überschreitung d. 1.1-fachen Erwartungswer\
tes v. S(N): - m. Normalverteilungs-Approximation: 0.0000014341 - mit Log-No\
rmalvert.-Approximation: 0.0000038793
______________________________________________________________________
Anmerkungen:
1. Literatur:
Dehling, H. & Haupt, B.: "Einführung in die Wahrscheinlichkeitstheorie und Statistik", Springer-Verlag
Berlin Heidelberg 2003
Gentle, J. E.: "Random Number Generation and Monte Carlo Methods", Springer-Verlag New York
1998.
Klugman, S. A., Panjer, H. H. & Willmot, G. E.: "Loss models. From data to decisions", 2. Auflage,
John Wiley & Sons Hoboken New Jersey 2004
Schmidt, K. D.: "Versicherungsmathematik", Springer-Verlag Berlin Heidelberg 2002
2. Aufgaben:
Variieren Sie die Verteilungstypen und Verteilungsparameter in den Prozeduren SN,
Hist, Hist2, EMPIR1 und EMPIR2 und
(a) vergleichen Sie theoretische und empirische Mittelwerte;
(b) studieren Sie die Gültigkeit des zentralen Grenzwertsatzes.
_____________________________________________________________________