Optionsbewertung - Heston Modell - Zhus Transformed Volatility Scheme

Hallo,

ich weiss nicht, ob mir hier jdm weiterhelfen kann, da es schon ein wenig speziell ist. aber ich habe folgendes problem:

Ich habe in C++ (programmierung siehe unten) das Heston Modell zwecks Bewertung einer einfachen Call-Option programmiert mit einer Monte Carlo Simulation. Diskretisiert habe ich es mit "the transformed volatility scheme" von Zhu mit der moment-matching methode.
Nun mein Problem: Ich habe die gleichen Parameter gewählt wie Zhu in seinem Paper. aber während bei ihm preise von rund 47 rauskommen, erhalte ich immer preise um die 56. ich find aber den fehler in meiner programmierung nicht und wäre für eure hilfe sehr dankbar!!!!

Um Kritik an der (schlechten) Generierung von Zufallszahlen vorweg zu nehmen:
für eine einfache Monte Carlo Simulation mit konstanter Volatilität waren diese zumindest ausreichend. (Heisst: identische Preise wie über den analytischen Weg des Black&Scholes Modelles) Und vor allem so große Preisunterschiede dürften darin nicht begründet sein.

zum mir: ich bin bwl'er und habe daher für gewöhnlich nichts mit programmiersprachen am hut. daher bitte ich um nachsicht für meinen vllt. schlechten programmier-stil. mir gings bisher nur darum, dass es funktioniert.^^ (was es bei einfacheren Modellen ja auch tat.)
ich hab mir C++ nur selber in ein paar wochen beibringen müssen (ohne irgend eine andere sprache zu können).

VIELEN DANK im voraus!!!
wäre klasse, falls jdm den Fehler finden könnte!

grüße
tobias

----------------------------------------------------------------------------
die c++ programmierung:

// HESTON MODELL - Transformed Volatility Scheme - Moment Matching Method

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <ctime>

using namespace std;

// Eingabe der Parameter für die Optionsberechnung:
double T=6; // Laufzeit der Option
double K=70; // Strike
double S0=100; // Startkurs Basiswert
double Vola=0.15; // Start-Vola (sigma)
double r=0.04; // risikoloser Zins (c.c.)
double Speed=2; // Stärke des Mean Reverse (kappa)
double meanVar=0.16; // durchscnittliche langfristige Varianz (theta)
double VolVol=0.3; // Die Volatilität der Varianz des Underlying (xi)
unsigned long Pfade=20000; // Anzahl Akitenpfad-Simulationen
double KnotenProJahr=32; // Anzahl der unterjährigen Kursbeobachtungen (time steps)
double p=(-0.5); // Korrelation der ZVen zwischen Basiswert und Volatilität (rho)
// Ende Parametereingabe

double Var; // Varianz des Basiswertes
double KnotenGesamt=KnotenProJahr*T; // Anzahl der gesamten Kursbeobachtungen über gesamte Laufzeit
double t=(1/KnotenProJahr); // Dauer zwischen den Knoten (in Jahren)

// 1. Generierung normalverteilter Zufallszahlen
double NVZV()
{
double result;
double x;
double y;
double sizeSquared;
do
{
x = 2.0000rand()/static_cast<double>(RAND_MAX)-1; // erstellt GVZV von -1 bis +1
y = 2.0000rand()/static_cast<double>(RAND_MAX)-1;
sizeSquared = x*x + y*y;
}
while
( sizeSquared >= 1.0000);
while
( sizeSquared <= 0.0000);
result = x*sqrt(-2*log(sizeSquared)/sizeSquared);
return result;
}
// Ende 1.

// 2. Monte-Carlo-Simulation
double MonteCarlo
(double T,
double K,
double S0,
double Vola,
double r,
double KnotenGesamt,
unsigned long Pfade)
{

double Summe=0;

double ST; // Aktienendkurs
//xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
for(unsigned long i=0; i <Pfade; i++) // Anzahl der Simulationen
{
double PerformanceKnoten;
double PerformanceGesamt=1;

//xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
for(unsigned long i=0; i <KnotenGesamt; i++) // Anzahl der Knoten
{
Var=VolaVola;
double Z=NVZV();
double q=1-(pp);
double W= sqrt(q)*NVZV() + p*Z;

PerformanceKnoten=exp(((r- 0.5*Var)*t)+(sqrt(t)*Vola*Z));
PerformanceGesamt *= PerformanceKnoten;

double alpha = meanVar+(Var-meanVar)*exp(-Speed*t)-(VolVol*VolVol/(4*Speed)(1-exp(-Speedt)));
alpha = alpha >0 ? alpha : 0;
double beta = sqrt(alpha);
double mmm= (beta-Vola*exp(-0.5*Speed*t))/(1-exp(-0.5*Speed*t)); // mmm = moment matched mean level

Vola += (0.5*Speed*(mmm - Vola)*t) + (0.5*VolVol*sqrt(t)*W);

}
//xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
ST=S0*PerformanceGesamt;
double Auszahlung = ST - K;
Auszahlung = Auszahlung >0 ? Auszahlung : 0;
Summe += Auszahlung;

}
//xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
double average = Summe / Pfade;
average = exp(-rT);
return average;
}
// Ende 2.

// Hauptprogramm:
int main()
{
// START UHR
$ clock_t $ start, end;
double dif;
start = clock();

// Start: für zeitabhängig generierte Zufallsvariablen
long n = time(0);
srand(n);
// Ende: für zeitabhängig generierte Zufallsvariablen

double Ergebnis = MonteCarlo(T, K, S0, Vola, r, KnotenGesamt, Pfade);
cout << "Der Preis der Call-Option beträgt: " <<Ergebnis << " EUR" << endl << endl;

// Ende Uhr
end = clock();
dif = double (end - start) / CLOCKS_PER_SEC;
cout << "Benoetigte Zeit: " << dif << "s" << endl;
cout <<endl<<endl<<endl<<endl<<endl;

system("PAUSE");

Hallo, mir scheint beim Durchlesen, daß deine Zufallsroutine schon falsch ist.

Die Bedingung sollte eher 'do {..} while(sizeSquared <= 0 || sizeSquared >= 1);' oder so heißen.

Aber nach der Änderung bekomme ich auch ~ 56EUR raus. Ich hab mir dann mal http://www.l-p-a.de/_files/cms/files/HestonSimulation.pdf reingezogen (das ist die Veröffentlichung, auf die Du dich beziehst, oder?).

Da wird der Zufallsprozess ganz anders mit (V_0, Theta, sigma, rho) parametriert und wenn ich die einsetze, kommt auch einigermaßen das Erwartete raus.