//paramètres
stacksize(40000000);
T=1;//horizon
sig=0.2;//volatilité
r=0.05;//intérêt
S0=100;//valeur initiale du sous-jacent
K=110;// Strike du Put considéré
N=2;//nombre initial de pas de discrétisation
M=2000000;//taille des vecteurs de simulations indépendantes
I=10;//nombre d'itérations
MI=M*I;//nombre de simulations indep

erreul=[];//vecteur des erreurs faibles Euler + Romberg
liceul=[];//largeur des intervalles de confiance 95% Euler
errmil=[];//vecteur des erreurs faibles Milshtein + Romberg
licmil=[];//largeur des intervalles de confiance 95% Milshtein
for j=1:2,//boucle sur le nombre de pas
  
  //variables de stockage
  moyeul=0;
  careul=0;
  moymil=0;
  carmil=0; 
  
  
  
  //paramètres utiles pour la discrétisation avec 2N pas
  sigt=sig*sqrt(T/(2*N));
  rt=r*T/N;
  rt2=rt/2;
  der=rt-sigt^2; 
  
  for i=1:I,//nombre de simulations vectorielles indep
  S=S0*ones(1,M);
  Se=S0*ones(1,M);//schéma d'Euler avec N pas
  Se2=S0*ones(1,M);//schéma d'Euler avec 2N pas
  Sm=S0*ones(1,M);//schéma de Milshtein avec N pas
  Sm2=S0*ones(1,M);//schéma de Milshtein avec 2N pas
	
	for k=1:N,//boucle sur les pas de temps
	  g1=rand(1,M,'g');
	  g2=rand(1,M,'g');
	  g=g1+g2;
	  S=S.*exp(sigt*g+der);
	  Se=Se.*(1+sigt*g+rt);
	  Se2=Se2.*(1+sigt*g1+rt2).*(1+sigt*g2+rt2);
	  Sm=Sm.*(1+sigt*g+rt+sigt^2.*(g^2-2)/2);
	  Sm2=Sm2.*(1+sigt*g1+rt2+sigt^2*(g1^2-1)/2).*(1+sigt*g2+rt2+sigt^2*(g2^2-1)/2);
	end;
	//contribution de la trajectoire de S au Put
	paymc=(S<K).*(K-S);
	//contribution des schémas d'Euler au Put
	payeul=2*(Se2<K).*(K-Se2)-(Se<K).*(K-Se)-paymc;
	//contribution des schémas de Milshtein au Put
	paymil=2*(Sm2<K).*(K-Sm2)-(Sm<K).*(K-Sm)-paymc;

  moyeul=moyeul+sum(payeul);
  careul=careul+sum(payeul.*payeul);
  moymil=moymil+sum(paymil);
  carmil=carmil+sum(paymil.*paymil);
  end;
 moyeul=moyeul/MI;
 careul=careul/MI;
 moymil=moymil/MI;
 carmil=carmil/MI;
 erreul=[erreul,exp(-r*T)*moyeul];
 liceul=[liceul,1.96*exp(-r*T)*sqrt((careul-moyeul^2)/MI)];
 errmil=[errmil,exp(-r*T)*moymil];
 licmil=[licmil,1.96*exp(-r*T)*sqrt((carmil-moymil^2)/MI)];
  
  N=N*2;//multiplication du nombre N de pas par 2
end;

erreul
liceul
errmil
licmil