clear; close all;

% 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;
mu = mean(M)';
Sigma = cov(M);
N = length(mu);  

% utility elicitation 
alpha = 3;
load(sprintf('EXP1_ALPHA=%d.mat',alpha));
x = utility(2).x;
y = utility(2).y;
port_our = umax(x,y,0,M/4);
u_our = ucal(x,y,M*port_our/4);

% Markowizt portfolios
u_mkw = []; % collecting utility values
d_mkw = []; % collecting distances from our portfolio
r_range = 0:0.01:1;

port_mkw = sdpvar(N,1);
constraints = [port_mkw >= 0, sum(port_mkw) == 1];
ops = sdpsettings('solver','mosek','verbose',0);
    
for i = 1:length(r_range)
    r = r_range(i);
    diagnosis = optimize(constraints, ...
        -r*port_mkw'*mu + (1-r)*sqrt(port_mkw'*Sigma*port_mkw), ops); %#ok<NASGU> 
    u_mkw = [u_mkw, ucal(x,y,M*double(port_mkw)/4)]; %#ok<AGROW>
    d_mkw = [d_mkw, norm(port_our-double(port_mkw))]; %#ok<AGROW> 
end

% Markowitz selection
[~,idx] = min(d_mkw);
r_opt = r_range(idx);
diagnosis = optimize(constraints, ...
    -r_opt*port_mkw'*mu + (1-r_opt)*sqrt(port_mkw'*Sigma*port_mkw), ops);
port_mkw_opt = double(port_mkw);
u_mkw_opt = ucal(x,y,M*port_mkw_opt/4);

fprintf('>> Optimal portfolios [ours, Markowitz] \n');
disp([port_our, port_mkw_opt]);

fprintf('>> Expected utilities [ours, Markowitz] \n');
disp([u_our, u_mkw_opt]);

fprintf('>> Mean-variance objectives [ours, Markowitz]\n');
disp([-r_opt*port_our'*mu + (1-r_opt)*sqrt(port_our'*Sigma*port_our), ...
    -r_opt*port_mkw_opt'*mu + (1-r_opt)*sqrt(port_mkw_opt'*Sigma*port_mkw_opt)]);

fprintf('>> Closest gamma \n');
disp(r_opt);

