//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=2;//nombre initial de pas de discrétisation
M=10000;//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);



//fonction de payoff du put
function y=put(S,K)
  y=0;
  if (S<K)
	y=K-S;
  end;
endfunction;




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



for j=1:5,//boucle sur le nombre de pas
  
  //variables de stockage
  puteul=0;
  careul=0;
  putmil=0;
  carmil=0;
  puteul=0;
  contputeul=0;
  contcareul=0;
  contputmil=0;
  contcarmil=0;
  
  //paramètres utiles pour la discrétisation avec N pas
  
  //////////////////////////////////////////////////
  //A compléter avec le calcul des paramètres utiles
  //////////////////////////////////////////////////

  for i=1:M,//simulations indépendantes
	S=S0;
	Se=S0;
	Sm=S0;
	for k=1:N,//boucle sur les pas de temps
	  g=rand(1,'g');//génération d'une gaussienne centrée réduite
	  
	  ////////////////////////////////////////////////////////////
	  //À compléter avec l'evolution du sous-jacent et des schémas
	  ////////////////////////////////////////////////////////////
	  
	end;
	//contribution de la trajectoire de S au Put
	paymc=put(S,K);
	//contribution de la trajectoire schéma d'Euler au Put
	payeul=put(Se,K);
	puteul=puteul+payeul;
	careul=careul+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
	////////////////////////////////////////////////////////////
	
	//contribution de la trajectoire schéma de Milshtein au Put
	
	////////////////////////////////////////////////////////////
	//Reprendre ce qui précède en remplaçant Euler par Milshtein
	////////////////////////////////////////////////////////////
	
  end;
  erreul=[erreul,exp(-r*T)*puteul/M-BS];
  liceul=[liceul,1.96*exp(-r*T)*sqrt((careul/M-(puteul/M)^2)/M)];
  errmil=[errmil,exp(-r*T)*putmil/M-BS];
  licmil=[licmil,1.96*exp(-r*T)*sqrt((carmil/M-(putmil/M)^2)/M)];
  conterreul=[conterreul,exp(-r*T)*contputeul/M];
  contliceul=[contliceul,1.96*exp(-r*T)*sqrt((contcareul/M-(contputeul/ ...
												  M)^2)/M)];
  conterrmil=[conterrmil,exp(-r*T)*contputmil/M];
  contlicmil=[contlicmil,1.96*exp(-r*T)*sqrt((contcarmil/M-(contputmil/ ...
												  M)^2)/M)];
  Npas=[Npas,N];
  N=N*2;//multiplication du nombre N de pas par 2
end;

erreul
liceul
errmil
licmil
conterreul
contliceul
conterrmil
contlicmil

//representation graphique de (1/N,conterreul) et (1/N,conterrmil)
xbasc();
subplot(2,1,1);
plot2d(1../Npas',[conterreul;conterreul-contliceul;conterreul+contliceul]');
subplot(2,1,2);
plot2d(1../Npas',[conterrmil;conterrmil-contlicmil;conterrmil+contlicmil]');