clear;
stacksize(50000000);
//paramètres
T=2;//horizon
r=0.05;//intérêt
S0=100;//valeur initiale du sous-jacent
sig=0.2;//volatilité
K=90;//strike
N=1;//nombre initial de pas de discrétisation
M=2000000;//nombre de simulations indépendantes


trapprix=[];//vecteur des prix méthode des trapèzes
traplic=[];//largeur des intervalles de confiance 95%
rectprix=[];//vecteur des prix méthode des rectangles
rectlic=[];//largeur des intervalles de confiance 95%
corprix=[];//vecteur des prix avec correction
corlic=[];//largeur des intervalles de confiance 95% Pont
Npas=[];// vecteur des nombres de pas

for j=1:5,//boucle sur les nombres de pas
  
 //paramètres utiles pour la discrétisation avec N pas
  sigt=sig*sqrt(T/N);
  sigtt=sigt/(sqrt(12));
  rt=r*T/N;
  der=rt-sigt^2/2; 
  S=S0*ones(1,M);
  trapsomS=0.5*S0*ones(1,M);
  rectsomS=zeros(1,M);
  somW=zeros(1,M);
  varKV=0;
  corsomS=zeros(1,M);
  for k=0:N-1, //boucle sur les pas de temps
	varKV=varKV+(N+0.5-k)^2;
	g1=rand(1,M,'g');
	somW=somW+(N+0.5-k)*g1;
	g2=sigtt*rand(1,M,'g');
	corsomS=corsomS+S.*(1+rt/2+sigt*g1/2+g2);
	S=S.*exp(sigt*g1+der);
	rectsomS=S+rectsomS;
  end;
  clear g1;
trapmoyS=(rectsomS+0.5*(S0-S))/N;
rectmoyS=rectsomS/N;
cormoyS=corsomS/N;
KV=S0*exp(sigt*somW/N+der*N/2);
sigKV=sig*sqrt(T*varKV/N^3);
clear S;
clear somS;
clear rectsomS;
clear corsomS;

//espérance de la variable de contrôle Kemma et Vorst
d2=(log(S0/K)+der*N/2)/sigKV;
d1=d2+sigKV;
espKV=exp(-r*T)*(S0*exp(der*N/2+sigKV^2/2)*cdfnor("PQ",d1,0,1)-K*cdfnor("PQ",d2,0,1));

//contribution de la moyenne trapèzes au prix du Call
trappay=exp(-r*T)*((trapmoyS>K).*(trapmoyS-K)-(KV>K).*(KV-K));
trappr=sum(trappay)/M;
trapprix=[trapprix,trappr+espKV];
traplic=[traplic,1.96*sqrt((sum(trappay^2)/M-trappr^2)/M)];
  
//contribution de la moyenne rectangles au prix du Call
rectpay=exp(-r*T)*((rectmoyS>K).*(rectmoyS-K)-(KV>K).*(KV-K));
rectpr=sum(rectpay)/M;
rectprix=[rectprix,rectpr+espKV];
rectlic=[rectlic,1.96*sqrt((sum(rectpay^2)/M-rectpr^2)/M)];

//contribution de la moyenne corrigée au prix du Call
corpay=exp(-r*T)*((cormoyS>K).*(cormoyS-K)-(KV>K).*(KV-K));
corpr=sum(corpay)/M;
corprix=[corprix,corpr+espKV];
corlic=[corlic,1.96*sqrt((sum(corpay^2)/M-corpr^2)/M)];
  
  Npas=[Npas,N];
  N=N*2;//multiplication du nombre de pas par 2
end;


trapprix
traplic
rectprix
rectlic
corprix
corlic

//representation graphique de (N,prixeul) et (N,prixmil) et de prixBS
clf();
subplot(2,2,1);
plot2d(Npas',[rectprix;rectprix-rectlic;rectprix+rectlic]');
xtitle(['rectangles'],'N','prix')
subplot(2,2,2);
plot2d(Npas',[trapprix;trapprix-traplic;trapprix+traplic]');
xtitle(['trapèzes'],'N','prix')
subplot(2,2,3);
plot2d(Npas',[corprix;corprix-corlic;corprix+corlic]');
xtitle(['correction'],'N','prix')
subplot(2,2,4);
plot2d(Npas',[trapprix;trapprix-traplic;trapprix+traplic;corprix;corprix-corlic;corprix+corlic]');
xtitle(['comparaison';'trapèzes/correction'],'N','prix')