clear;
stacksize(50000000);
//paramètres
T=1;//horizon
r=0.05;//intérêt
S0=100;//valeur initiale du sous-jacent
sig0=0.2;//valeur moyenne de la volatilité
alpha=2;//vitesse de retour à la moyenne de la volatilité
beta=0.1;//coefficient multiplicatif du bruit de la vol
K=100;//strike
N=10;//nombre de pas de discrétisation
M=100000;//nombre de simulations indépendantes

X=S0*ones(1,M);//sous-jacent actualisé
Y=zeros(1,M);//processus dirigeant la volatilité
somcarsig=zeros(1,M);//somme des carrés des volatilité
dt=T/N;
sqdt=sqrt(dt);
exal=exp(-alpha*dt);
betsqal=beta*sqrt((1-exal^2)/(2*alpha));
KrT=K*exp(-r*T);
	
for k=0:N-1, //boucle sur les pas de temps
  g1=rand(1,M,'g');
  g2=rand(1,M,'g');
  
  /////////////////////////////////////////
  //mettre à jour X puis Y et somcarsig//
  /////////////////////////////////////////
  
end;
clear Y;
clear g2;

//contribution des trajectoires au prix du Call
pay=(X>KrT).*(X-KrT);
clear X;


//contribution des trajectoires à l'espérance conditionnée
sigT=sqrt(somcarsig/N);
clear somcarsig;
d1=log(S0/KrT)./sigT+sigT/2;
BS=S0*cdfnor("PQ",d1,zeros(1,M),ones(1,M))-KrT*cdfnor("PQ",d1-sigT, ...
												  zeros(1,M),ones(1,M));
clear sigT;
clear d1;

//calcul des prix et des largeurs des intervalles de confiance
//sans conditionnement
prix=sum(pay)/M
licprix=1.96*sqrt((sum(pay^2)/M-(prix)^2)/M)
clear pay;
//avec conditionnement
condprix=sum(BS)/M
liccondprix=1.96*sqrt((sum(BS^2)/M-(condprix)^2)/M)