I did write such code but it is not in a good state to be shared right
now.
This is the code I used to for the independent normals:
#include <gsl/gsl_rng.h>
static float recip[] = {1.0, 0.5, 0.333333333333, 0.25, 0.2,
0.166666666666,
0.1428571428, 0.125, 0.111111111, 0.1, 0.09090909,
0.0833333333};
// n<=12 determines precision
static float log(float x, int n){
//if(x==1.0) return 0.0;
float log2 = 0.6931472f;
float inv=1.0f;
int mm;
double y = frexp(x, &mm)-1.0f;
double term = y;
double sum = 0.0;
for(int i=0; i<n; i++){
sum += recip[i]*term;
term *= -y;
}
return inv*(sum+mm*log2);
}
void
my_ran_ugaussian (const gsl_rng * r, double &g1, double &g2) {
double x, y, r2;
do {
x = -1 + 2 * gsl_rng_uniform (r);
y = -1 + 2 * gsl_rng_uniform (r);
r2 = x * x + y * y;
} while (r2 > 1.0 || r2 == 0);
double sq = sqrt (-2.0 * log (r2, 12) / r2);
g1 = x * sq;
g2 = y * sq;
}
void
my_ran_ugaussian (const gsl_rng * r, int n, valarray<double> &rands)
{
int half = n/2;
for(int i=0; i<half; i++){
my_ran_ugaussian (r, rands[i], rands[half+i]);
}
double dummy;
if(2*half != n) my_ran_ugaussian (r, rands[n-1], dummy);
}
where I just used the default uniform random number generator from gsl
(it's very good).
The rejection box-muller is basically a copy of the GSL code (GPL) but
without wasting half of the random gaussians. I also use a faster
logarithm since a lot of compute time seemed to be spent there.
Other things that might be useful:
Antithetics get you about a factor of four in variance reduction.
Using the par CDS levels as control variates also seems like a really good
idea but I havent tried it.
If you have other questions/suggestions please feel free to write.
Regards,
Karsten
On Mon, 19 Jul 2004 00:44:17 +0200, Nguyen Xuan Son
<
[hidden email]> wrote:
>
>
> Hi,
>
> Is there anyone has written the code of a (gaussian) copula? If Yes, can
> he
> share the code to me?
>
> If Not, is there any class which can return me a sequence of "N"
> independent
> Gaussian random variables ( I remember that QuantLib does have this
> class.).
>
> On the numerical side, what is the best random generator for this task
> if N
> is of 10^2 order, and the number of scenario should be 10^4.
>
> My thanks.
>
> Xs
>
Locked
