#include <stdlib.h>
#include <time.h>

const int M=365; // nombre de pas de temps
const double two_pi =6.28318530718;
double gauss(){
double x,y;
x= (1.+rand())/(1.+RAND_MAX);
y= (1.+rand())/(1.+RAND_MAX);
return sqrt( -2 *log(x) )*cos(two_pi*y);
}
double EDOstoch(double S0, double dt, double sdt, double rdt){
double S= S0;
int i; for(i=1;i<M;i++)
S= S*(1.+gauss()*sdt+rdt);
return S;
}

int main(argc, argv) int argc; char *argv[];
{
const double K=110, S0=100, T=1., r=0.03, sigma=0.2;
const int k, kmax=20000; // nb de realisations
double dt=T/M, sdt=sigma*sqrt(dt), rdt=r*dt, P0=0, t0=clock();
srand(time(NULL));
for(k=0; k<kmax;k++){
double Sa= EDOstoch(S0, dt, sdt, rdt);
if(K>Sa) P0 += K-Sa;
}
t0=(clock()-t0)/CLOCKS_PER_SEC;
printf("P_0 = %f CPUtime=%f \n",P0*exp(-r*T)/kmax, t0);
return 0;
}