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);
Y=zeros(1,M);
somcarsig=zeros(1,M);
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');
  X=X.*(1+sqdt*(sig0+Y).*g2);//attention à mettre X à jour avant Y
  Y=exal*Y+betsqal*g1;
  somcarsig=somcarsig+(Y+sig0)^2;
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
sigmoy=sqrt(somcarsig/N);
clear somcarsig;
d1=log(S0/KrT)./sigmoy+sigmoy/2;
BS=S0*cdfnor("PQ",d1,zeros(1,M),ones(1,M))-KrT*cdfnor("PQ",d1-sigmoy, ...
												  zeros(1,M),ones(1,M));
clear sigmoy;
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)