clear;
stacksize(50000000);
//paramètres
T=1;//horizon
r=0.05;//intérêt
S0=100;//valeur initiale du sous-jacent   
sig=0.2;//valeur moyenne de la volatilité
K=100;//strike
b=130;//barrière up and out : l'option est désactivée 
      //lorsque le sous-jacent dépasse b
N=1;//nombre initial de pas de discrétisation
M=100000;//nombre de simulations indépendantes


//expression explicite de E^*(1_{S_T<b}(S_T-K)^+) dans le modèle de Black-Scholes
function pr=pr(S0,K,b,r,sig0,T)
  d2p=(log(S0/K)+(r-sig^2/2)*T)/(sig*sqrt(T))
  d2m=(log(S0/b)+(r-sig^2/2)*T)/(sig*sqrt(T))
  d1p=d2p+sig0*sqrt(T)
  d1m=d2m+sig0*sqrt(T)
  pr=S0*(cdfnor("PQ",d1p,0,1)-cdfnor("PQ",d1m,0,1))-K*exp(-r*T)*(cdfnor("PQ",d2p,0,1)-cdfnor("PQ",d2m,0,1))
endfunction;

//prix Black-Scholes vol sigma par la méthode des images 
prixBS=pr(S0,K,b,r,sig,T)-(b/S0)^(2*r/sig^2-1)*pr(b^2/S0,K,b,r,sig,T)


prixeul=[];//vecteur des prix Euler
liceul=[];//demi-largeur des intervalles de confiance 95% Euler
prixpont=[];//vecteur des prix avec Pont
licpont=[];//demi-largeur des intervalles de confiance 95% Pont
prixdec=[];//vecteur des prix avec decalage
licdec=[];//demi-largeur des intervalles de confiance 95% decalage
Npas=[];// vecteur des nombres de pas

for j=1:4,//boucle sur les nombres de pas
  Se=S0*ones(1,M);//schéma d'Euler
  indb=(S0<b);//indicatrice non franchissement barriere
  probab=(S0<b)*ones(1,M);//proba condit non franchissement barrière
  sigt=sig*sqrt(T/N);
  indbdec=(S0<b-c0*sigt*S0);//indicatrice non franchissement barriere decalee
  rt=r*T/N;
  der=rt-sigt^2/2;
  
	
	
  for k=0:N-1, //boucle sur les pas de temps
	g=rand(1,M,'g');
	Seold=Se;
	///////////////////
	//mettre à jour Se//
	///////////////////
	indb=indb.*(Se<b);
	////////////////////////
	//mettre à jour probab//
	////////////////////////
	////////////////////////
	//mettre à jour indbdec//
	////////////////////////
  end;
	
	
  clear g;
  clear Sold;
   //contribution des trajectoires au prix du Call
  pay=exp(-r*T)*(Se>K).*(Se-K);
  payeul=pay.*indb;
  prix=sum(payeul)/M;
  prixeul=[prixeul,prix];
  liceul=[liceul,1.96*sqrt((sum(payeul^2)/M-prix^2)/M)];
  
//contribution des trajectoires au prix du Call avec technique de pont
  paypont=pay.*probab;
  prix=sum(paypont)/M;
  prixpont=[prixpont,prix];
  licpont=[licpont,1.96*sqrt((sum(paypont^2)/M-prix^2)/M)];

  //contribution des trajectoires au prix du Call avec technique de decalage
  paydec=pay.*indbdec;
  prix=sum(paydec)/M;
  prixdec=[prixdec,prix];
  licdec=[licdec,1.96*sqrt((sum(paydec^2)/M-prix^2)/M)];

Npas=[Npas,N];
  N=N*5;//multiplication du nombre de pas par 5
end;


prixeul
liceul
prixpont
licpont
prixdec
licdec
//representation graphique de (N,prixeul), (N,prixpont), (N,prixdec) et de prixBS
clf;
subplot(3,1,1);
plot2d(Npas',[prixeul;prixeul-liceul;prixeul+liceul]');
plot2d(Npas,prixBS*ones(Npas),style=5);
subplot(3,1,2);
plot2d(Npas',[prixpont;prixpont-licpont;prixpont+licpont]');
plot2d(Npas,prixBS*ones(Npas),style=5);
subplot(3,1,3);
plot2d(Npas',[prixdec;prixdec-licdec;prixdec+licdec]');
plot2d(Npas,prixBS*ones(Npas),style=5);