Box-Müller-Methode in C++

Hallo,
habe mir eine Monte-Carlo-Simulation in C++ programmiert zwecks Bewertung von Optionen. Anstatt der Polar-Methode habe ich mal die (instabilere) Box-Müller-Methode verwendet. Und nach nur 3119 Simulationen der pseudo-random-numbers wird als Preis bereits 1.#INF angegeben. (es ist immer die Anzahl 3119 weil ich die Zahlen bewusst nicht an die uhr koppele damit ich den Fehler finden kann.)
Hat jemand von euch eine Idee woran der Fehler liegen könnte? Wenn ich richtig programmiert habe, wäre die letzte random number = -1,09295.

Danke im voraus, für eventuelle Hilfe!!!

Die Programmierung sieht wie folgt aus:

// MONTE CARLO SIMULATION MIT KONSTANTER VOLATILITÄT

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <ctime>
#define PI 3.141592653589793 // PI muss erst definiert werden

using namespace std;

// Eingabe der Parameter für die Optionsberechnung:
double T=1; // Laufzeit in Jahren
double K=1000; // Strike
double S0=1000; // Startkurswert des Underlying
double Vola=0.40; // Konstante Volatilität
double r=0.05; // Drift bzw. risikolose Zins
unsigned long Pfade=3119; // Fehler bei 3119
double KnotenProJahr=12; // Kurs-Betrachtungszeiträume pro Jahr
// Ende Parametereingabe

double KnotenGesamt=KnotenProJahr*T; // Anzahl der gesamten Betrachtungszeitpunkten
double t=(1/KnotenProJahr); // Dauer zwischen den Betrachtungszeitpunkten (in Jahren)

// 1.1 Normalverteilte Zufallsvariable nach der Polar-Methode
double pNVZV()
{
double z;
double x;
double y;
double sizeSquared;
do
{
x = 2.0000rand()/static_cast<double>(RAND_MAX)-1; // erstellt gleichverteilte ZVen von -1 bis +1
y = 2.0000rand()/static_cast<double>(RAND_MAX)-1; // erstellt gleichverteilte ZVen von -1 bis +1
sizeSquared = x*x + y*y;
}
while
( sizeSquared >= 1.0000);
z = x*sqrt(-2*log(sizeSquared)/sizeSquared);
return z;
}
// Ende 1.1 Normalverteilte Zufallsvariable nach der Polar-Methode

// 1.2 Normalverteilte Zufallsvariable nach der Box-Müller-Methode
double NVZV()
{
double result;
double x = rand()/static_cast<double>(RAND_MAX); // erstellt gleichverteilte ZVen von 0 bis 1
double y = rand()/static_cast<double>(RAND_MAX); // erstellt gleichverteilte ZVen von 0 bis 1
result = cos(2 * PI * x)*sqrt(-2*log(y)); // erstellt normalverteilte ZVen von 0 bis 1
return result;
}
// Ende 1.2 Normalverteilte Zufallsvariable nach der Box-Müller-Methode

double z;
// 2. Monte-Carlo-Simulation
double MonteCarlo
(double T,
double K,
double S0, // Achtung: 0=NULL! Nicht der Buchstabe "o". Startaktienkurs
double Vola,
double r,
double KnotenGesamt,
unsigned long Pfade)
{
double Var= pow(Vola,2);
double Summe=0;

double ST; // Aktienendkurs
// Beginn äußere Schleife
for(unsigned long i=0; i <Pfade; i++) // Anzahl der Simulationen
{
double PerformanceKnoten;
double PerformanceGesamt=1;
// Beginn innere Schleife
for(unsigned long i=0; i <KnotenGesamt; i++) // Anzahl der Knoten
{
z=NVZV();
PerformanceKnoten=exp(((r- 0.5*Var)*t)+(sqrt(t)*z*Vola));
PerformanceGesamt = PerformanceKnoten;
ST=S0PerformanceGesamt;
}
// Ende innere Schleife

double Auszahlung = ST - K;
Auszahlung = Auszahlung >0 ? Auszahlung : 0;
Summe += Auszahlung;
}
// Ende äußere Schleife
double average = Summe / Pfade;
average = exp(-rT);
return average;
}
// Ende 2. Monte-Carlo-Simulation

// 3. Das Programm
int main()
{

// Ausführung der MCS Simulation und Ausgabe als Ergebnis
double Ergebnis = MonteCarlo(T, K, S0, Vola, r, KnotenGesamt, Pfade);
cout << "Der Preis der Call-Option beträgt: " <<Ergebnis << " EUR" << endl << endl;
// Ausführung der MCS Simulation und Ausgabe als Ergebnis

cout << z;

system("PAUSE");
}
// Emde 3. Das Programm
// ENDE PROGRAMM MONTE CARLO SIMULATION MIT KONSTANTER VOLATILITÄT

SeppJ

Deine Zufallszahlen, so wie du sie ziehst, können auch 0 sein. Bei Box-Müller ist aber Voraussetzung, dass dies nicht passiert.

Über den Codestil sage ich mal lieber nix.

zum inhalt:
also die WK mit der ZV=0 kenne ich. aber statistisch ist dieser fehler/risiko eher zu vernachlässigen, da die WK bei einer zahl mit 15komma-stellen (double-datentyp) bei 10^(-15) liegt. also das wird nicht der grund sein. den wie bereits erwähnt taucht der fehler ca. alle 3000-4000 ZV auf. also DEUTLICH überrepräsentativ. und wenn ichs korrekt programiert haben sollte, dann ist die letzte ZV ja was mit -1,xxx

zum code-stil:
ich bin bwler und kann KEINE EINZIGE programmiersprache und C++ hab ich mir in den letzten wochen selber beigebracht. von daher ist der stil für mich erstmal nebensächlich solange es läuft.^^

SeppJ

Was heißt hier statistisch? Du ziehst doch so wie das sehe immer die gleichen Zahlen. Also wird das immer die gleiche 0 sein und deswegen tritt der Fehler auch immer bei der gleichen Pfadtiefe auf. Guck doch einfach mal, welche Zahl da gezogen wird!

Also ich denke, daß hier "double" an seine Grenzen kommt. Es liegt nur indirekt an den Zufallszahen. Der Darstellungsbereich von "double" ist zwar hoch 2 hoch 64 Bits. Es wäre ein Versuch wert "double" durch "long double" zu ersetzen. Tritt der Fehler immer noch bei der selben Pfadtiefe auf?

SeppJ

Exhumed schrieb:

Also ich denke, daß hier "double" an seine Grenzen kommt. Es liegt nur indirekt an den Zufallszahen. Der Darstellungsbereich von "double" ist zwar hoch 2 hoch 64 Bits. Es wäre ein Versuch wert "double" durch "long double" zu ersetzen. Tritt der Fehler immer noch bei der selben Pfadtiefe auf?

Nein.

@SeppJ:
das ziehen der immer gleichen ZV hab ich doch jetzt nur gemacht UM den fehler zu finden. (hatte ich doch glaube ich auch zuvor geschrieben.)
bei ca. 5000 simulierungen (mit NICHT immer den gleichen ZV!) wird das ergebnis jedes mal "zerschossen". bei 3000 nur ca. jedes zweite mal. also wenn ich "statistisch" schreibe, dann bezieht sich das natürlich NICHT auf immer die gleichen ZV-reihen!
und wie auch bereits im ausgangs-thread geschrieben:
die letzte ZV (bei "festen" ZV an 3119. stelle) ist NICHT null sondern -1,09295!!!

SeppJ schrieb:

Exhumed schrieb:

Also ich denke, daß hier "double" an seine Grenzen kommt. Es liegt nur indirekt an den Zufallszahen. Der Darstellungsbereich von "double" ist zwar hoch 2 hoch 64 Bits. Es wäre ein Versuch wert "double" durch "long double" zu ersetzen. Tritt der Fehler immer noch bei der selben Pfadtiefe auf?

Nein.

=> Doch.
also mit nur 3000 simulationen (mit den nicht immer identischen ZV) ist die Warscheinlichkeit, dass das "Ergebnis" 1.#INF kommt bei gefühlten. 50%. und das trotz des "long double" datentypes. es hat also wie ich mir schon dachte nichts mit dem datentyp zu tun. weil da hätte double eigentlich ausreichen müssen. denn wie auch schon erwähnt kann es nicht an der null liegen!
die wahrscheinlichkeit dafür wäre bei nur 3000 simulationen bei über eins zu einer milliarde. (wenn es sich wirklich um eine zahl mit 15 kommastellen handelt wie es beim double-typ eigentlich sein sollte)
und dann bei mehreren Simulationen, wäre die WK das jedes mal eine 0,000000000000000 auftaucht weit höher als 1:Trillionen....

SeppJ

Tobias K. schrieb:

=> Doch.

Ich sage, es liegt nicht am Datentyp, dann sagst du doch und dann erklärst du, warum es nicht am Datentyp liegt?

Außerdem hat die WK dass eine 0 kommt, nichts mit deinem Datentyp zu tun, sondern mit deinem Zufallsgenerator. Guck dir mal dessen RAND_MAX an.

Und irgendwie weiß ich nicht, wie man dir helfen soll. Du argumentierst bei einem deterministischen Programm mit irgendwelchen Wahrscheinlichkeiten und so. Wenn dein Programm einen Fehler hat, dann guck das Programm an, wie ich dir gesagt habe. Guck dir an, welche Zufallszahl gezogen wird, wenn NAN kommt. Reparier den genannten Fehler in der Box-Müller Methode. Wechsel den Seed. Wechsel den Zufallsgenerator. Wenn du einfach alles ablehnst und (zu unrecht) meinst dass es nicht sein kann, dann ist dir nicht zu helfen.

@SeppJ:
sry, da habe ich dich dann falsch verstanden.
Denn der letzte Satz war
"Tritt der Fehler immer noch bei der selben Pfadtiefe auf?"
und deine antwort war "NEIN". also nahm ich an, dein "nein" bezog sich auf den letzten satz und nicht auf seinen vorschlag.

SeppJ

tobias kempkens schrieb:

@SeppJ:
sry, da habe ich dich dann falsch verstanden.
Denn der letzte Satz war
"Tritt der Fehler immer noch bei der selben Pfadtiefe auf?"
und deine antwort war "NEIN". also nahm ich an, dein "nein" bezog sich auf den letzten satz und nicht auf seinen vorschlag.

Ahh, ok. Tschuldige dass das missverständlich war. Ich wollte Exhumed darauf hinweisen, dass dies ganz offensichtlich kein Genauigkeitsproblem ist, aber auch nicht seitenlang über Numerik schreiben um ihm zu erklären, warum.

edit: Falls möglich, wäre auch die rand-Implementierung deiner Standardbibliothek hilfreich. Bei meiner (libstdc++ auf 64-Bit Debian) tritt der Fehler nämlich nicht auf, selbst bei vielen Millionen Pfaden. Da spielt dann eventuell die Wahrscheinlichkeit der 0 eine Rolle, weil ich vielleicht Glück hatte und einfach keine bekommen habe.

Schau dir mal den Wert von RAND_MAX an. Dieser sagt dir wie viele verschiedene Zahlen zwischen 0 und 1 du höchstens haben kannst. Wahrscheinlich gibt es da nur 2^16 Zahlen. Was die Genauigkeit von double bei weitem unterschreitet.

also das könnte natürlich der grund sein. ich bin natürlich absoluter informatik-laie und dachte, dass mit dem typen double dann auch eben wirklich 15 kommastellen zur verfügung ständen und alle werte im intervall [0;1] entsprechend gleichverteilt wären.

da ich leider nicht weiss wo ich sowas für den bereich RAND_MAX nachgucken kann oder evtl ändern hab ich einfach mal folgendes gemacht:

die gleichverteilte ZV (welche später logarithmiert wird) so bestimmt, dass wenn sie =0 ist, stattdessen den wert 0,000001 annimmt.
und ihr hattet recht. selbst mit 10 Millionen Simulation bleibt das Ergebnis stabil bzw. errechenbar. (so wie es eigentlich seien sollte.)

also wäre das problem an der stelle gelöst. (warum es ZU instabil war).

noch eine lästige frage:
könntet ihr jdm der (quasi) NULL PLAN von der ganzen materie hat (ausser dem zusammenbasteln) erklären wie ich den befehl
rand()/static_cast<double>(RAND_MAX)
dann entsprechend umwandeln könnte, damit er in etwa ZV im double format zieht?
gibts einen anderen befehl in der standard-bib? kann man den befehl iwie erweitern? muss ich mir erst ne andere bibliothek-datei besorgen?

der GRUND warum ich eigentlich die suboptimal Box-Müller-Methode genommen habe anstatt die bessere Polar-Methode ist, dass ich gehört habe, dass BM-Methode für C++ instabil wäre bei bereits zahlen NAHE null!!! und ich gerne rausfinden möchte inwiefern.

auf auf jeden fall schonmal DANKE für eure mühe!

JaykopX

Um Zufallszahlen zu bekommen die größer sind als RAND_MAX hab ich mal ein Beispiel gebastelt für mein RAND_MAX = 32767 = 2^15 -1:

#include <iostream>

using namespace std;

#define MAX_UINT64 18446744073709551615 //2^64-1
double RandomNumber() {
	unsigned long long retVal=0;//Rückgabewert
	retVal = rand();			//Ersten 15 bits füllen, da RAND_MAX bei mir 2^15
	for(int i=0; i<3; ++i) {		//3*15 bits nachfüllen
		retVal <<=15;			//15 Bits verschieben
		retVal |= rand();		//nächsten 15 Bits füllen
	}
	return (double)retVal/(double)(MAX_UINT64>>4);	//Wir nutzen nur 4*15=60 bits, also 4 bits "abziehen"
}

int main() {
	cout<<"Wert von RAND_MAX:"<<RAND_MAX<<'\n';
	for(int i=0; i<100; ++i) {
		cout<<RandomNumber()<<" ## ";
	}
}

Das ganze ohne Gewähr! Da könnte noch was faul dran sein, hab es nicht ausgiebig getestet. Das ganze sollte nun eine 60 bit Zufallszahl liefern, ich habe aber nicht nachgeprüft ob da durch Rundungsfehler was verloren geht. Die Idee dahinter ist einfach das man eine Zufallszahl aus 4 Zufallszahlen zusammenbastelt.

Michael E.

Das verkürzt die Periode auf ein Viertel und hat wahrscheinlich bescheidene Eigenschaften bzgl. Zufälligkeit. Nimm lieber nen PRNG, der für größere Zahlen ausgelegt ist.

SeppJ

Also für Monte Carlo ist rand() sowieso nicht gut, weil man nie weiß was dahinter steckt. Normalerweise ist das ein ganz einfacher LCG und der ist für ernsthaftes MC absolut ungeeignet. Such dir am besten einen guten, schnellen Zufallsgenerator (multiply and carry, TAUS, ...) und implementier den selbst oder nimm eine fertige Implementierung, z.B. aus der GSL. Da ist der Zahlenraum dann auch mindestens 32 Bit, du findest bestimmt auch welche mit 64-Bit.

o.k., danke für den tipp / die tipps.
(wobei rand() für diese monte carlo simulation noch auszureichen schien. denn bei 50Mio. Simulationen war das numerische Ergebnis identisch mit dem analytischen aus dem Black-Scholes-Modell zur Optionspreisbewertung.)

kurze frage nochmal des informatik-laien. hab in der library hinsichtlich des befehls RAND_MAX folgendes gefunden:

#define RAND_MAX 0x7FFF

was genau sagt jetzt dieses 0x7FFF aus?

SeppJ

Tobias K. schrieb:

was genau sagt jetzt dieses 0x7FFF aus?

32767

Und das erklärt auch, warum bei dir die 0 so wahrscheinlich ist und warum ich auf dem Computer an dem ich getestet habe (RAND_MAX = 2147483647) viele Millionen Schritte ohne 0 machen konnte.

(Und es sagt, dass ich mit meiner Antwort ganz zu Anfang zu 100% Recht hatte )

Das ist auf jeden Fall schonmal interessant zu wissen.

was genau müsste ich denn da statt 0x7FFF eingeben um auch 2147483647 mögliche zahlen zu erhalten? oder kann ich da nicht einfach so rumfummeln?

SeppJ

Tobias K. schrieb:

oder kann ich da nicht einfach so rumfummeln?

Dies. Der Wert ist dazu da, dich darüber zu informieren, wie rand() implementiert ist. Der Wert wird von der rand()-Implementierung selbst aber nicht benutzt. D.h. du kannst gerne was ändern, aber das Verhalten von rand() ändert sich nicht. Aber die Programme die du compilierst würden dann falsche Annahmen über RAND_MAX machen.