View Single Post
Old 16-11-2011, 12:00   #7
mamo139
Senior Member
 
L'Avatar di mamo139
 
Iscritto dal: Sep 2006
Città: Bologna/Milano
Messaggi: 525
Quote:
Originariamente inviato da shinya Guarda i messaggi
Mersenne Twister?
DUE MIGLIORIE:
ieri notte ho separato la generazione dei numeri random immagazzinandoli in un array, quindi ora ci sono due for al posto di uno ma c'è stata una miglioria:

da 0.44 a 0.34

questa mattina ho sostituito il tr1 con:
una versione del Mersenne Twister per generare numeri random e una funzione che pare essere molto ottimizzata per normalizzarli (spero sia fatta abbastanza bene da non crearmi problemi a livello di qualità dei dati )
per esequire 10 simulazioni da 1.000.000 ci mette ora 0.30 al posto di 0.34 del tr1
ma bisogna migliorare ancora

eccovi il codice allo stato attuale:
Codice:
inline double Margrabe_option_MC_price(double r, double t, double s1, double s2, double q1, double q2, 
								double std1, double std2, double corr, long iterations){
	//variables and objects initialization
	register long x,y;
	double price = 0; //the final price will be stored here
	double payoff;
	double *z1; //to store generated random numbers
	double matr21, matr22;
	double sigma1, sigma2;
	double part1, part2, sqrtt;

	//memory allocation
	z1 = new double[iterations*2];

	//seeding the random number generator
#ifdef RANDOM_MODE_TR1
	Myeng eng;
	ndistr stdnorm(0, 1); //creating the standard normal distribution generator
	eng.seed((unsigned long) time(NULL));
#endif
#ifdef RANDOM_MODE_Mersenne_twister	
    init_genrand((int)time(NULL));
#endif
	//calculating sigmas
	matr21 = corr * std2;
	matr22 = sqrt(1-corr*corr) * std2;

	//precalculating fixed part of drive formula
	part1 = exp((r - q1 - 0.5*std1*std1)*t);
	part2 = exp((r - q2 - 0.5*std2*std2)*t);
	sqrtt = sqrt(t);

	for(x=0; x < iterations ;x++){
		y=x+x;

		//generating random numbers
		#ifdef RANDOM_MODE_TR1
		z1[y] = stdnorm(eng);
		z1[y+1] = stdnorm(eng);
		#endif
		#ifdef RANDOM_MODE_Mersenne_twister	
		z1[y] = normsinv(genrand_real3());
		z1[y+1] = normsinv(genrand_real3());
		#endif
	}

	for(x=0; x < iterations ;x++){
		y=x+x;

		//creating sigmas
		sigma1 = z1[y] * std1;
		sigma2 = matr21 * z1[y] + matr22 * z1[y+1];

		//calculating payoff
		payoff = (s1 * part1 * exp(sqrtt * sigma1)) - (s2 * part2 * exp(sqrtt * sigma2));
		if(payoff > 0) price += payoff;
	}
	price /= iterations;
	price *= exp((-r)*t);

	delete [] z1;

	return price;
}
Quote:
Originariamente inviato da starfred Guarda i messaggi
se matlab usa il multithreading sarà molto difficile andar più veloce con un single thread... non credo che il C++ basti, potresti provare soluzioni asm miste c++
si matlab di default usa multithreading in tutte le funzioni matematiche con vettori superiori a qualche migliaio di dati. Ma lo sto facendo andare su un 2 core, quindi possiamo batterlo secondo me

io purtroppo asm non lo conosco, ma se tu hai voglia di farmi vedere come potrei sostituire qualche pezzetto con codice asm apprezzerei molto (e poi così mi metto anche a vedere come funziona l'assembly)
__________________
http://mamo139.altervista.org
mamo139 è offline   Rispondi citando il messaggio o parte di esso