//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=10000;//nombre de simulations indépendantes

//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 + 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:5,//boucle sur le nombre de pas
   
  
  //variables de stockage
  puteul=0;
  careul=0;
  putmil=0;
  carmil=0;
  
  //paramètres utiles pour la discrétisation avec N et 2N pas
  
  //////////////////////////////////////////////////
  //A compléter avec le calcul des paramètres utiles
  //////////////////////////////////////////////////
  
  for i=1:M,//simulations indépendantes
	S=S0;
	Se=S0;//schéma d'Euler avec N pas
	Se2=S0;//schéma d'Euler avec 2N pas
	Sm=S0;//schéma de Milshtein avec N pas
	Sm2=S0;//schéma de Milshtein avec 2N pas
	
	for k=1:N,//boucle sur les pas de temps
	  g1=rand(1,'g');
	  g2=rand(1,'g');
	  g=g1+g2;
	  
	  ////////////////////////////////////////////////////////////
	  //À 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 des schémas d'Euler au Put
	payeul=2*put(Se2,K)-put(Se,K)-paymc;
	puteul=puteul+payeul;
	careul=careul+payeul^2;
	//contribution des schémas de Milshtein au Put
	paymil=2*put(Sm2,K)-put(Sm,K)-paymc;
	putmil=putmil+paymil;
	carmil=carmil+paymil^2;
  end;
  erreul=[erreul,exp(-r*T)*puteul/M];
  liceul=[liceul,1.96*exp(-r*T)*sqrt((careul/M-(puteul/M)^2)/M)];
  errmil=[errmil,exp(-r*T)*putmil/M];
  licmil=[licmil,1.96*exp(-r*T)*sqrt((carmil/M-(putmil/M)^2)/M)];
  
  N=N*2;//multiplication du nombre N de pas par 2
end;

erreul
liceul
errmil
licmil