//paramètres
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=1;//nombre initial de pas de discrétisation
M=100000;//nombre de simulations indépendantes

//formule Black-Scholes pour le prix du Put
d1=(log(S0/K)+r*T)/(sig*sqrt(T))+sig*sqrt(T)/2;
d2=d1-sig*sqrt(T);
BS=K*exp(-r*T)*cdfnor("PQ",-d2,0,1)-S0*cdfnor("PQ",-d1,0,1);




erreul=[];//vecteur des erreurs faibles Euler
liceul=[];//demi-largeur des intervalles de confiance 95% Euler
conterreul=[];//vecteur des erreurs faibles Euler avec variable de contrôle
contliceul=[];//demi-largeur des intervalles de confiance 95%
Npas=[];// vecteur des nombres de pas



for j=1:6,//boucle sur le nombre de pas
  
  
  //paramètres utiles pour la discrétisation avec N pas
  
  //////////////////////////////////////////////////
  //A compléter avec le calcul des paramètres utiles
  //////////////////////////////////////////////////

  
	S=S0*ones(1:M);
	Se=S0*ones(1:M);
	for k=1:N,//boucle sur les pas de temps
	  g=rand(1:M,'g');//génération d'un vecteur de gaussiennes centrées réduites
	  
	  ////////////////////////////////////////////////////////////
	  //À compléter avec l'evolution du sous-jacent et du schéma
	  ////////////////////////////////////////////////////////////
	  
	end;
	//contribution des trajectoires de S au Put
	paymc=(K-S).*(K-S);
	//contribution des trajectoires du schéma d'Euler au Put
	payeul=(K-Se).*(K-Se);
	puteul=sum(payeul);
	careul=sum(payeul^2);
	
	////////////////////////////////////////////////////////////
	//À compléter en stockant dans contputeul (resp. contcareul)
	//la somme (resp. somme des carrés) des différences entre payoff 
	//avec Se et payoff avec S
	////////////////////////////////////////////////////////////
	
	

  erreul=[erreul,exp(-r*T)*puteul/M-BS];
  liceul=[liceul,1.96*exp(-r*T)*sqrt((careul/M-(puteul/M)^2)/M)];
  conterreul=[conterreul,exp(-r*T)*contputeul/M];
  contliceul=[contliceul,1.96*exp(-r*T)*sqrt((contcareul/M-(contputeul/M)^2)/M)];
  Npas=[Npas,N];
  N=N*2;//multiplication du nombre N de pas par 2
end;

erreul
liceul
conterreul
contliceul

//representation graphique de (h=T/N,conterreul)
clf;
plot2d(T./Npas',[conterreul;conterreul-contliceul;conterreul+contliceul]');