%% 2.3, 2.4
clc, clf, clear all, close all;

mean = [-3.0588 -2.5371 -2.3961]; %Means
std = [0.7202 0.9899 1.5508]; %Standard deviations
alpha = [0.95 0.99 0.999];
m = 1000; %Number of obligors
l = 0.6; %Default loss

ind_def_prob = [];      %individual default probability
pair_def_corr = [];     %pairwise default correlation
rho = [];               %mutual dependence of obligors

for (i = 1:3)
    mix_prob = @(x)  ((1 ./ (1 + exp(-(mean(i) + std(i).*x)))) .* (1./sqrt(2*pi)).*exp((-x.^2)/2)); %Mixing probability p(Z) * F(Z)
    ind_def_prob(i) = integral(mix_prob, -Inf, Inf); %E[p(Z)]
end

for (i = 1:3)
    mix_prob = @(x)  ((1 ./ (1 + exp(-(mean(i) + std(i).*x)))).^2) .* (1./sqrt(2*pi)).*exp((-x.^2)/2); %Mixing probability p(Z) * F(Z)
    EpZ2 = integral(mix_prob, -Inf, Inf); %E[p(Z)^2]
    pair_def_corr(i) = (EpZ2 - (ind_def_prob(i)^2)) / (ind_def_prob(i)*(1-ind_def_prob(i)));
end

for (i = 1:3)
    fun = pair_def_corr(i) * ind_def_prob(i) * (1 - ind_def_prob(i)) + ind_def_prob(i)^2;
    EpZ2 = @(x, Z) (((normcdf((norminv(ind_def_prob(i)) - Z .* sqrt(x)) ./ sqrt(1 - x))).^2) .* (normcdf((1 ./ sqrt(x)) .* (sqrt(1 - x) .* norminv(Z) - norminv(ind_def_prob(i))))) - fun);
    g = @(x) integral(@(Z) EpZ2(x,Z), 0, 1);
    rho(i) = lsqnonlin(g, 0.21, 0, 1);
end

VaR_95 = [];
ES_95 = [];

for (i = 1:length(mean))
    VaR_95(i) = l * m * normcdf(sqrt(rho(i))*norminv(alpha(1))+norminv(ind_def_prob(i))/(sqrt(1-rho(i))));
    ES_95(i) = (1/(1-alpha(1)))*integral(@(x) l * m * normcdf(sqrt(rho(i))*norminv(x)+norminv(ind_def_prob(i))/(sqrt(1-rho(i)))), alpha(1), 1);
end

VaR_99 = [];
ES_99 = [];

for (i = 1:length(mean))
    VaR_99(i) = l * m * normcdf(sqrt(rho(i))*norminv(alpha(2))+norminv(ind_def_prob(i))/(sqrt(1-rho(i))));
    ES_99(i) = (1/(1-alpha(1)))*integral(@(x) l * m * normcdf(sqrt(rho(i))*norminv(x)+norminv(ind_def_prob(i))/(sqrt(1-rho(i)))), alpha(2), 1);
end

VaR_999 = [];
ES_999 = [];

for (i = 1:length(mean))
    VaR_999(i) = l * m * normcdf(sqrt(rho(i))*norminv(alpha(3))+norminv(ind_def_prob(i))/(sqrt(1-rho(i))));
    ES_999(i) = (1/(1-alpha(1)))*integral(@(x) l * m * normcdf(sqrt(rho(i))*norminv(x)+norminv(ind_def_prob(i))/(sqrt(1-rho(i)))), alpha(3), 1);
end