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
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é
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);
KrT=K*exp(-r*T);
	
	
for k=0:N-1, //boucle sur les pas de temps
	g1=rand(1,M,'g');
	g=c1*g1+c2*rand(1,M,'g');
	X=X.*(1+sqdt*(sig0+Y).*g);//attention à mettre X à jour avant Y
	Y=exal*Y+betsqal*g1;
end;
clear Y;
clear g1;
clear g2;

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

//contribution des trajectoires au prix du Call avec variable de controle
//////////////////////////////////////////////
//calculez dans coef un estimateur de gamma*//
//////////////////////////////////////////////
coef=

paycont=pay-coef*(X-S0);


//calcul des prix et des largeurs des intervalles de confiance
//sans variable de controle
prix
varpay=(sum(pay^2)/M-(prix)^2);
licprix=1.96*sqrt(varpay/M)
//avec variable de controle
contprix=sum(paycont)/M
varcont=varpay-coef^2*(sum(X^2)/M-S0^2);
liccontprix=1.96*sqrt(varcont/M)