% PORTFOLIO_FUNCTION.m 
%
% Solves a portfolio problem with CVaR objective using 
%
% 1. Sampling Method
% 2. Moments Method
% 3. Segregated Moments Method
%

function [x_sol, objval] = portfolio_function(cur_data, beta)

% extract moments and parameters from historical data
mean_ret  = mean(cur_data)';
covar_ret = cov(cur_data);
tau   = mean(mean_ret);
N = length(mean_ret);
T = size(cur_data, 1);

% partitioned statisitcs
part_data = [pos(cur_data), neg(cur_data)];
part_mean = mean(part_data)';
part_covar = cov(part_data);

% output arrays
x_sol = zeros(N, 3);
objval = zeros(1, 3);

%%%%%%%%%%%%%%%%%%%%%%
% 1. Sampling Method %
%%%%%%%%%%%%%%%%%%%%%%
h = rome_begin('Portfolio Example (Sampling)'); % make model

% Decisions
newvar v;           % cvar aux variable
newvar x(N) nonneg; % asset weights
newvar y(T) nonneg; % aux variable

% Objective
rome_minimize(v + (1./((1-beta)*T)) * sum(y));

% Constraints
rome_constraint(mean_ret' * x >= tau); % Target Mean Return
rome_constraint(sum(x) == 1);          % Sum of Weights
rome_constraint(y >= -cur_data*x - v); % Auxilliary constraint

% Completing the Model
h.solve;                     % solve
x_sol(:, 1) = h.eval(x); % get portfolio weights
objval(1)   = h.objective;   % get objective value
rome_end;                    % clear memory

%%%%%%%%%%%%%%%%%%%%%
% 2. Moments Method %
%%%%%%%%%%%%%%%%%%%%%
h = rome_begin('Portfolio Example (Moments)');  % make model

% Uncertainties
newvar r(N) uncertain;  % represents uncertain returns
r.set_mean(mean_ret);   % specify mean
r.Covar = covar_ret;    % specify covariance

% Decisions
newvar v;               % cvar aux variable
newvar x(N) nonneg;     % asset weights
newvar y(1, r) linearrule nonneg;  % auxilliary variable

% Objective
rome_minimize(v + (1./(1-beta)) * mean(y));

% Constraints
rome_constraint(mean_ret' * x >= tau); % Target Mean Return
rome_constraint(sum(x) == 1);          % Sum of Weights
rome_constraint(y >= -r'*x - v); % Auxilliary constraint

% Complete 
h.solve_deflected;           % solve
x_sol(:, 2) = h.eval(x); % get portfolio weights
objval(2) = h.objective;     % get objective value
rome_end;                    % clear memory

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Partitioned Moments Method %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h = rome_begin('Portfolio Example (Partitioned)'); % make model

% Uncertainties
newvar s(2*N) uncertain nonneg;  % Uncertainties: partitioned returns
s.set_mean(part_mean);           % Specify partitioned mean
s.Covar = part_covar;            % Specify partitioned covariance

% Actual returns: function of uncertainties
r = [eye(N), -eye(N)] * s;

% Decisions
newvar v;               % cvar aux variable
newvar x(N) nonneg;     % asset weights
newvar y(1, s) linearrule nonneg;  % auxilliary variable

% Objective
rome_minimize(v + (1./(1-beta)) * mean(y));

% Constraints
rome_constraint(mean_ret' * x >= tau); % Target Mean Return
rome_constraint(sum(x) == 1);          % Sum of Weights
rome_constraint(y >= -r'*x - v); % Auxilliary constraint

% Complete 
h.solve_deflected;           % solve
x_sol(:, 3) = h.eval(x); % get portfolio weights
objval(3) = h.objective;     % get objective value
rome_end;                    % clear memory