//- Profil pertes cdo square
clear 
xdel(winsid());
stacksize(100000000)
//- La probabilite de defaut depend du taux d interet r
p = 0.03;
//- Tranches du CDO de base
a = 0.03;
b = 0.09;
//- Nombre total d'emprunteur dans un panier
N_total = 500;

//- Tranches du CDO squared
a_sq = 0.03;
b_sq = 0.09;
//- Nombre de CDOs consideres dans un CDO squared
M = 100;
//- Nombre de simulations de Monte Carlo
Nb_mc = 10000;
//- Quantile
alpha=0.01;

//-----------------------------------
//- Incertitude sur p
//-----------------------------------


sg=0.1; // sg entre 0 et 1
p_values=[ p*(1-sg), p*(1+sg)];
p_proba=[0.5 0.5];
p_cumul=cumsum([0,p_proba]);
U=ones(M,1)*rand(1,Nb_mc);
N=0;
for i=1:length(p_proba)
  N=N+((U<p_cumul(i+1))&(U>p_cumul(i))).*grand(M,Nb_mc,'bin',N_total, ...
                                             p_values(i));
end
loss=1/(N_total*(b-a))*min(max(N-a*N_total,0),(b-a)*N_total);
loss=sum(loss,'r')/M; //- Perte normailsée entre 0 et 1
loss=1/(b_sq-a_sq)*min(max(loss-a_sq,0),b_sq-a_sq);

//---------------------------
//- Recherche d un quantile
[xk,k]=gsort(-loss,'g','i');
VaR=-xk(floor(alpha*Nb_mc));
CVaR=-sum(xk(1:floor(alpha*Nb_mc)))/floor(alpha*Nb_mc);
disp(VaR,'La VaR est')
disp(CVaR,'La CVaR est')
//---------------------------
histplot(20,loss)