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

clc
clear

% level of expectile
q = 0.05; 
% constraint bounds
b = 0:(0.0013262/99):0.0013262;
% b = 0.00105;
% 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_min_risk, type = minimize',...
        'Objective: objective_expectile_risk',...
        'fnext_nexpectile(headerpoint, valuespoint, matrix_scenarios, q)',...
        ...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_avg_return, lower_bound = ', num2str(b(t))],...
        'avg_g(matrix_scenarios)',...
        ...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('Average return: %f\n',...
        sum(matrix_scenarios * optimal_vector)/size(matrix_scenarios, 1));
    objective_values(t) = tbpsg_objective(solution_str, outargstruc_arr);
    Z(:, t) = optimal_vector;
end

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

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