% This example maximizes average returns with a negative expectile
% constraint using data from "Case Study: Basic CVaR Optimization Problem,
% Beyong Black-Litterman."

clc
clear

% level of expectile
q = 0.05;
% constraint bound
b = 0.025:((0.039-0.025)/99):0.039;
% solver = VAN, TANK, CAR, BULDOZER, VANGRB, CARGRB, HELI
solver = 'VANGRB';
% Load data
load('.\expectiles_problem_data.mat')

D = size(matrix_scenarios, 2);
headerpoint = ['x1' sprintf('\tx%d', 2:D)];
valuespoint = ones(1, D);
matrix_linear = ones(1, D);

objective_values = zeros(1, length(b));
Z = zeros(D, 100);
for t = 1:length(b)
problem_statement = sprintf('%s\n',...
    ...define problem name and type
    'Problem: expectiles_test_max_return, type = maximize',...
    'Objective: objective_avg_return',...
    'avg_g(matrix_scenarios)',...
    ...define linear constraint
    'Constraint: constraint_linear, lower_bound = 1, upper_bound = 1',...
    ...define linear function that implements linear constraint
    'linear_balance(matrix_linear)',...
    ...define external constraint
    ['Constraint: constraint_expectile_risk, upper_bound = ', num2str(b(t))],...
    'fnext_nexpectile(headerpoint, valuespoint, matrix_scenarios, q)',...
    ...define bounds on variables
    'Box_of_Variables: lowerbounds = 0, upperbounds = 1',...
    ...define solver type, precision
    ['Solver: ' solver ', precision = 9']);

 iargstruc_arr(1) = matrix_pack('matrix_linear', matrix_linear,headerpoint);
 iargstruc_arr(2) = matrix_pack('matrix_scenarios', matrix_scenarios,headerpoint);
 
 [solution_str, outargstruc_arr] = mpsg_solver(problem_statement, iargstruc_arr);
 
fprintf('mpsg_solver function using:\nSolution status: %s\n%s\nOptimal point:\n',...
    get_solution(solution_str, outargstruc_arr, 'solution_status'), ...
    get_solution(solution_str, outargstruc_arr, 'objective'));
optimal_point = get_optimalpoint(solution_str, iargstruc_arr, outargstruc_arr);
optimal_vector = zeros(D, 1);
for i=1:D
    fprintf('%s\t=\t%f\n', optimal_point{i, 1}, optimal_point{i, 2});
    optimal_vector(i) = optimal_point{i, 2};
end
fprintf('Sum of x: %f\n', sum(optimal_vector));
fprintf('Negative expectile: %f\n',...
    fnext_nexpectile(headerpoint, optimal_vector, matrix_scenarios, q));
    objective_values(t) = tbpsg_objective(solution_str, outargstruc_arr);
    Z(:, t) = optimal_vector;
end

figure
plot(b, objective_values)
title('Efficient Frontier')
xlabel('Upper Bound on Negative Expectile Risk')
ylabel('Average Returns')

figure
[X, Y] = meshgrid(b, 1:4);
mesh(X,Y,Z,'FaceAlpha',0.8);
title('Optimal Portfolio as Function of Negative Expectile Constraint')
xlabel('Upper Bound on Negative Expectile Risk')
ylabel('Instruments')
zlabel('Percentage Invested')