//- Profil pertes cdo square
clear 
xdel(winsid());
stacksize(100000000)
//---------------------------------------------------------------
//- Variation de la moyenne, de la variance et du quantil selon
//- differents tranches et probabilites de defaut 
//----------------------------------------------------------------

p=[0.005 0.01 0.03 0.05 0.1 0.2 0.3 0.5];
tranches=[0,0.03,0.03,0.06,0.06,0.09,0.15,1];

//- Quantile
alpha = 0.01;

//- Nb d'emprunteur dans un panier
N_total = 500; 

//- Nombre de CDOs consideres dans un CDO squared
M= 100;
//- Nb d'iterations Monte Carlo
Nb_mc = 10000;

a_sq = 0.03;
b_sq = 0.06;
for j=1:length(tranches)/2
  a=tranches(2*(j-1)+1);
  b=tranches(2*(j-1)+2);
  for i=1:length(p)
    N = grand(M,Nb_mc,'bin',N_total,p(i));
    loss=1/(N_total*(b-a))*min(max(N-a*N_total,0),(b-a)*N_total);
    loss=sum(loss,'r')/M; //- Perte normalisée entre 0 et 1
    loss=1/(b_sq-a_sq)*min(max(loss-a_sq,0),b_sq-a_sq);
    m_loss(i,j)=mean(loss);
    v_loss(i,j)=stdev(loss);
    [xk,k]=gsort(-loss,'g','i');
    VaR(i,j)=-xk(floor(alpha*Nb_mc));
    CVaR(i,j)=-sum(xk(1:floor(alpha*Nb_mc)))/floor(alpha*Nb_mc);
  end
end

xset("window",1)
xtitle("Moyenne des pertes","Probabilite de defaut","Pertes");
plot2d(p',m_loss(:,1),rect=[0,0,p(length(p)),1.1],style=2)
plot2d(p',m_loss(:,2),style=3)
plot2d(p',m_loss(:,3),style=4)
plot2d(p',m_loss(:,4),style=5)
legends(['0%-3%','3%-6%','6%-9%','15%-100%'],[2,3,4,5],'ur')
xset("window",2)
xtitle("Ecart type des pertes","Probabilite de defaut","Ecart type");
plot2d(p',v_loss(:,1),style=2)
plot2d(p',v_loss(:,2),style=3)
plot2d(p',v_loss(:,3),style=4)
plot2d(p',v_loss(:,4),style=5)
legends(['0%-3%','3%-6%','6%-9%','15%-100%'],[2,3,4,5],'ur')
xset("window",3)
xtitle("VaR","Probabilite de defaut","Quantile");
plot2d(p',VaR(:,1),rect=[0,0,p(length(p)),1.1],style=2)
plot2d(p',VaR(:,2),style=3)
plot2d(p',VaR(:,3),style=4)
plot2d(p',VaR(:,4),style=5)
legends(['0%-3%','3%-6%','6%-9%','15%-100%'],[2,3,4,5],'ur')
xset("window",4)
xtitle("CVaR","Probabilite de defaut","CVaR");
plot2d(p',CVaR(:,1),rect=[0,0,p(length(p)),1.1],style=2)
plot2d(p',CVaR(:,2),style=3)
plot2d(p',CVaR(:,3),style=4)
plot2d(p',CVaR(:,4),style=5)
legends(['0%-3%','3%-6%','6%-9%','15%-100%'],[2,3,4,5],'ur')