% Reproducing example in section 3.2.2, Hickman e Wollman (2002) 

clear
nsim = 20000;
PD = [.0001 .0002 .0006 .0018 .0106 .0494 .1914];
s = [.4 3.6 21.1 71.9 128.9 197.6];            % average exposure
rho = .5;
LGD = .3;
e = [.3 66.6 681.5 461 385 200;
1.8 236.3 1067.9 619.5 264.9 611.6;
237.7 1493.4 2027.3 68 243.9 207.8;
364.6 2709.9 2664.3 713 239.8 0;
232.2 1717.8 2389.5 480.7 141.9 163.8
61.5 563.5 1760.9 637 142.5 200;
64.5 386.5 309.1 182.6 0 0];
L_y = zeros(length(PD),length(s),nsim);
y = randn(nsim,1);
for i = 1:length(PD)
    for j = 1:length(s)
        n = round(e(i,j)/s(j));            % number of loans
        p_y = normcdf((norminv(PD(i)) - rho*y) ./ ((1-rho^2)^.5));      % Vasicek formula
        L_y(i,j,:) = binornd(n,p_y,nsim,1)*s(j)*LGD;
    end
end

% computing subportfolio loss distribution

L_tot = zeros(nsim,1);
for i = 1:length(PD)
    for j = 1:length(s)
        L_tot = L_tot + squeeze(L_y(i,j,:));
    end
end
quantili = quantile(L_tot,[.5 .75 .9 .95 .99 .999 .9999])/sum(sum(e));
quantili
mean(L_tot)
skewness(L_tot)
kurtosis(L_tot)