clear;
stacksize(50000000);
//paramètres
T=1;//horizon
r=0.05;//intérêt
S0=100;//valeur initiale du sous-jacent
sig0=0.2;//valeur moyenne de la volatilité
rho=-0.5;//corrélation des browniens
alpha=2;//vitesse de retour à la moyenne de la volatilité
beta=0.1;//coefficient multiplicatif du bruit de la vol
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-sig0^2/2)*T)/(sig0*sqrt(T))
  d2m=(log(S0/b)+(r-sig0^2/2)*T)/(sig0*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 sigma0 (cas beta=0) par la méthode des images 
prixBS=pr(S0,K,b,r,sig0,T)-(b/S0)^(2*r/sig0^2-1)*pr(b^2/S0,K,b,r,sig0,T)


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

for j=1:4,//boucle sur les nombres de pas
  S=S0*ones(1,M);
  Y=zeros(1,M);
  indb=(S0<b)*ones(1,M);//indicatrice non franchissement
  probab=(S0<b)*ones(1,M);//proba condit non franchissement
  dt=T/N;
  sqdt=sqrt(dt);
  exal=exp(-alpha*dt);
  betsqal=beta*sqrt((1-exal^2)/(2*alpha));
  c1=rho*sqrt(2*(1-exal)^2/(alpha*dt*(1-exal^2)));
  c2=sqrt(1-c1^2);
  KrT=K*exp(-r*T);
	
	
  for k=0:N-1, //boucle sur les pas de temps
	g1=rand(1,M,'g');
	g=c1*g1+c2*rand(1,M,'g');
	Sold=S;
	S=Sold.*(1+sqdt*(sig0+Y).*g+r*dt);//attention à mettre S et probab à jour avant Y
	indpas=(S<b);
	indb=indb.*indpas;
	probab=probab.*indpas.*(1-exp(-2*(b-Sold).*(b-S)./(dt*Sold.*Sold.*(sig0+Y)^2)));
	Y=exal*Y+betsqal*g1;
  end;
	
	
  clear g1;
  clear g2;
  clear Sold;
  clear Y;
  //contribution des trajectoires au prix du Call
  pay=exp(-r*T)*(S>K).*(S-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)];
  Npas=[Npas,N];
  N=N*5;//multiplication du nombre de pas par 5
end;


prixeul
liceul
prixpont
licpont

//representation graphique de (N,prixeul) et (N,prixmil) et de prixBS
xbasc();
subplot(2,1,1);
plot2d(Npas',[prixeul;prixeul-liceul;prixeul+liceul]');
if (beta==0) 
  plot2d(Npas',prixBS*ones(Npas)');
end;
subplot(2,1,2);
plot2d(Npas',[prixpont;prixpont-licpont;prixpont+licpont]');
if (beta==0) 
  plot2d(Npas',prixBS*ones(Npas)');
end;