% based on bootstrapping
clear; close all;

% handle for computing a discrete distribution's entropy
entropy = @(prob) -sum(prob.*log10(max(prob,1e-8)));

% compute stats of yearly total returns based on 1927-2023
M = csvread('10_Industry_Portfolios.csv', 0,1); %#ok<CSVRD> 
M = 1 + M(end-19:end,1:10)/100;

% utility elicitation 
alpha = 5;
load(sprintf('EXP1_ALPHA=%d.mat',alpha));
x = utility(2).x;
y = utility(2).y;
u = @(x) mean((1-exp(-alpha*x))/(1-exp(-alpha)));

% simulate returns (R) of length N
SIM = 1e2;  % # simulations

radius_all = 0:5e-4:5e-3;
entropy_all = zeros(1,length(radius_all));
expu_all = zeros(1,length(radius_all));
entropy_all_all = zeros(SIM,length(radius_all));
expu_all_all = zeros(SIM,length(radius_all));

for sim = 1:SIM
    rng(sim);
    I = randperm(size(M,1));
    I = I(1:size(M,1)/2);
    R_is = M(I,:); % in-sample returns
    R_os = M(setdiff(1:size(M,1),I),:); % out-of-sample returns
    for i = 1:length(radius_all)
        port = umax(x,y,radius_all(i),R_is/4);
        entropy_all(i) = entropy(port);
        u_approx = ucal(x,y,R_os*port/4);
        u_exact = u(R_os*port/4);
        expu_all(i) = u_exact;
    end
    entropy_all_all(sim,:) = entropy_all;
    expu_all_all(sim,:) = expu_all;
end

f1 = figure; hold on;
multiplot_combined(radius_all, ...
    (expu_all_all./repmat(expu_all_all(:,1),[1 length(radius_all)])-1)'*100, ...
    [0.3010 0.7450 0.9330]/2, f1);
plot(radius_all, zeros(size(radius_all)), '-.', 'color', [0 0 0]/2, 'linewidth', 2);
grid on;
box on;
xlabel('$$\delta$$', 'Interpreter','latex');
ylabel('Gain in Expected Utility (%)');
axis([-inf +inf -0.8 2.0]);
title("$\alpha = $ " + alpha, 'Interpreter', 'Latex');

f2 = figure; hold on;
multiplot_combined(radius_all, ...
    entropy_all_all', ...
    [0.4660 0.6740 0.1880]/2, f2);
grid on;
box on;
axis([-inf +inf 0 1]);
xlabel('$$\delta$$', 'Interpreter','latex');
ylabel('Portfolio''s Entropy');
title("$\alpha = $ " + alpha, 'Interpreter', 'Latex');


function [y_output] = multiplot_combined(x_data, y_data, color, f)

% select or create a figure to draw ---------------------------------------
if nargin <= 3
    figure;
    hold on;
else
    figure(f);
    hold on;
end
pastel = min(1.6*color, 1);
% -------------------------------------------------------------------------
% remove NaN columns
y_data = y_data(:,all(~isnan(y_data)));
% calculate means and percentiles
y_m = mean(y_data, 2);
% y_m = mean(y_data, 50, 2);
y_u = prctile(y_data, 90, 2);
y_l = prctile(y_data, 10, 2);
% -------------------------------------------------------------------------
fill([x_data(1:1:end), x_data(end:-1:1)], [y_u; y_l(end:-1:1)], pastel, 'FaceAlpha', 1/2);
plot(x_data, y_u, 'Color', pastel, 'Linewidth', 2);
plot(x_data, y_l, 'Color', pastel,  'Linewidth', 2);
plot(x_data, y_m, '-o', 'Color', color, 'Linewidth', 2);
y_output = [y_l y_m y_u];
% -------------------------------------------------------------------------
grid on;
box on;
axis([min(x_data), max(x_data), -inf, + inf]);

end

