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é
rho=-0.8;//corrélation des browniens
alpha=2;//vitesse de retour à la moyenne de la volatilité
beta=0.1;//coefficient multiplicatif du bruit de la vol
K=100;//strike
KrT=K*exp(-r*T);//strike actualisé
N=10;//nombre de pas de discrétisation
M=100000;//nombre de simulations indépendantes

//prix du Call Black-Scholes vol sig0
sigT=sig0*sqrt(T);
d1=log(S0/KrT)/sigT+sigT/2;
BS=S0*cdfnor("PQ",d1,0,1)-KrT*cdfnor("PQ",d1-sigT,0,1);


X=S0*ones(1,M);//sous-jacent actualisé
Y=zeros(1,M);//processus dirigeant la volatilité
Xbs=S0*ones(1,M);//modèle de Black-Scholes actualisé
dt=T/N;
sqdt=sqrt(dt);
exal=exp(-alpha*dt);
betsqal=beta*sqrt((1-exal^2)/(2*alpha));
c1=rho*sqrt(2*(1-exal)^2/(alpha*dt*(1-exal^2)));
c2=sqrt(1-c1^2);
	
	
for k=0:N-1, //boucle sur les pas de temps
	g1=rand(1,M,'g');
	g=c1*g1+c2*rand(1,M,'g');
	
	//////////////////////////////////
    //mettre à jour X et Xbs puis Y //
    //////////////////////////////////
end;
clear Y;
clear g1;
clear g2;

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


//contribution des trajectoires au prix du Call avec variable de controle

//////////////////////////////////////////////////////////////////////
//stockez dans paycont le vecteur pay moins les payoff Black-Scholes//
//////////////////////////////////////////////////////////////////////


clear Xbs;

//calcul des prix et des largeurs des intervalles de confiance
//sans variable de controle
prix=sum(pay)/M
licprix=1.96*sqrt((sum(pay^2)/M-(prix)^2)/M)
//avec variable de controle
contprix=BS+sum(paycont)/M
liccontprix=1.96*sqrt((sum(paycont^2)/M-(sum(paycont)/M)^2)/M)