close all, clear all
% A1 Uni(65,140)
% A2 Norm(120,40)
% A3 Lnorm(4.5,0.9)
% A4 Exp(100)
% A5 Weib(105,2)

% a %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fill the parts XXXXXXX
t = 0:0.5:200;
pd1 = makedist('Uniform','lower',65,'upper',140);
a1 = pdf(pd1,t);
pd2 = makedist('Normal','mu',120,'sigma',40);
a2 = pdf(pd2,t);
pd3 = makedist('Lognormal','mu',4.5,'sigma',0.9);
a3 = pdf(pd3,t);
pd4 = makedist('Exponential','mu',100);
a4 = pdf(pd4,t);
pd5 = makedist('Weibull','a',105,'b',2);
a5 = pdf(pd5,t);

% b %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%create the plot of PDFs
%type 'help pdf' to see what the values a1,...,a5 are and if they can be
%used for plotting
plot(t,a1,t,a2,t,a3,t,a4)
xlabel('monetaty values')
ylabel('probability')
title('PDFs')
legend('A1','A2','A3','A4','location','northwest')

% c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%create samples
samples =zeros(5,5000);
samples(1,:) = random(pd1,1,5000);
samples(2,:) = random(pd2,1,5000);
samples(3,:) = random(pd3,1,5000);
samples(4,:) = random(pd4,1,5000);
samples(5,:) = random(pd5,1,5000);

%vectors to save the results
EV = zeros(5,1);
Var1 = zeros(5,1);
Var5 = zeros(5,1);
Var10 = zeros(5,1);
CVar1 = zeros(5,1);
CVar5 = zeros(5,1);
CVar10 = zeros(5,1);

%calculate the results, some basic idea provided here
%look at slide 10 from lecture 4
%the for loop handles the 5 investment opportunities one at a time
for i = 1:5
    %EV
    EV(i) = mean(samples(i,:));
    %VaRs
    Var1(i) = prctile(samples(i,:),1);
    Var5(i) = prctile(samples(i,:),5);
    Var10(i) = prctile(samples(i,:),10);
    
    %tailindices
    tailIndices1 = find(samples(i,:) <= Var1(i));
    tailIndices5 = find(samples(i,:) <= Var5(i));
    tailIndices10 = find(samples(i,:) <= Var10(i));
    
    %CVaRs
    CVar1(i) = mean(samples(i,tailIndices1));
    CVar5(i) = mean(samples(i,tailIndices5));
    CVar10(i) = mean(samples(i,tailIndices10));
end


% d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%labels for the scatter plots
labels = {'A1','A2','A3','A4','A5'};

%subplot command is useful
% H = subplot(m,n,p), or subplot(mnp), breaks the Figure window
%     into an m-by-n matrix of small axes, selects the p-th axes for
%     the current plot, and returns the axes handle.

%you can uncomment and test the code below
% x = 0:0.01:2*pi;
% subplot(1,3,1)
% plot(x,sin(x))
% subplot(1,3,2)
% plot(x,cos(x))
% subplot(1,3,3)
% plot(x,tan(x))
% ylim([-5,5])
figure
subplot(3,2,1)
plot(EV,Var1,'o')
text(EV,Var1,labels)
xlabel('Expected value')
ylabel('Risk measure')
title('VaR 1%')

subplot(3,2,2)
plot(EV,CVar1,'o')
text(EV,CVar1,labels)
xlabel('Expected value')
ylabel('Risk measure')
title('CVaR 1%')

subplot(3,2,3)
plot(EV,Var5,'o')
text(EV,Var5,labels)
xlabel('Expected value')
ylabel('Risk measure')
title('VaR 5%')

subplot(3,2,4)
plot(EV,CVar5,'o')
text(EV,CVar5,labels)
xlabel('Expected value')
ylabel('Risk measure')
title('CVaR 5%')

subplot(3,2,5)
plot(EV,Var10,'o')
text(EV,Var10,labels)
xlabel('Expected value')
ylabel('Risk measure')
title('VaR 10%')

subplot(3,2,6)
plot(EV,CVar10,'o')
text(EV,CVar10,labels)
xlabel('Expected value')
ylabel('Risk measure')
title('CVaR 10%')

% g %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%first create a matrix that can contain the CVaR results for all the 5
%invesment oppportunities with 100 different alpha values
%then use for loop to calculate the 5*100 cases separately

%plot the CVaRs of each investment opportunity using alpha = 1:1:100
Vars = zeros(5,100);
CVars = zeros(5,100);

for i = 1:5
    for j = 1:100
        Vars(i,j) = prctile(samples(i,:),j);
        tailIndices = find(samples(i,:) <= Vars(i,j));
        CVars(i,j) = mean(samples(i,tailIndices));
    end
end
alpha = 1:1:100;
figure
hold on
for i = 1:5
    plot(alpha,CVars(i,:),'linewidth',1.5) 
end
xlabel('$\alpha$','interpreter','latex')
ylabel('CVaR')
legend('A1','A2','A3','A4','A5','location','northwest')

% h %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% previously we used function pdf(), there is a similar function for
% cumulative distibution functions
t = 0:300;

a1cdf = cdf(pd1,t);
a5cdf = cdf(pd5,t);

figure
plot(t,a1cdf,t,a5cdf)
xlabel('monetary value')
ylabel('probability')
title('CDFs')
legend('A1','A5','location','northwest')

% i %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% analogous to part h
t2 = 0:600;

a2cdf = cdf(pd2,t2);
a3cdf = cdf(pd3,t2);

figure
plot(t2,a2cdf,t2,a3cdf)
xlabel('monetary value')
ylabel('probability')
title('CDFs')
legend('A2','A3','location','northwest')