% This case study investigates two investment strategies for a portfolio of hedge funds.  
% The first strategy rebalances the portfolio by solving an optimization
% problem.
% The second  strategy, called "20-best", every time period selects 20 hedge funds 
% with the highest return (over available data) and give them equal weights. 
% Rebalancing in both strategies is done on monthly basis. 
% Out-of-sample numerical experiments demonstrate  investment strategies in 
% a simulated historical environment.  We consider the portfolio optimization with various
% risk measures: CVaR, CDaR, MAD,... and Variance (code can handle any risk measure 
% included to PSG). Also, we take into account the market-neutrality (i.e., beta-zero) constraint. 
% Performance of  strategies is compared with SP500 benchmark.
% This m-file optimizes problems with Riskprog subroutine. 
% Analogously CS_HedgeFunds_RiskInObjective_psg_solver.m implements 
% the same case study with mpsg_solver optimization subroutine. 

function CS_HedgeFunds_RiskInObjective_riskprog_MAD()

clear;
% set TYPE parameter 
% TYPE='RISK' to optimize portfolio with risk in objective
% TYPE='VARIANCE' to optimize portfolio with variance in objective
TYPE = 'RISK';

%initial settings
InitialPortfolioValue = 10000.0; % initial value of portfolio
BudgetBound = 1.0; % budget bound
K = 0.000001; % K value in risk neutrality constraint
FirstInsampleIndex = 1;
InsampleSize = 12; % in sample set
DefaultLambda = 1.0; % default lambda value
VariablesLowerBound = 0.0; % lower box bounds 
VariablesUpperBound = 0.1; % upper box bounds 
Precision = 4; % default precision
Stages = 30;
BestFundsNumber = 20; % the number of the best returns funds in portfolioEW
SaveSolutionFile = 1; % set this key to write solving statistics to file solution_HedgeFunds.txt

% read funds and benchmark matrices
clear header_main; clear date_vec; clear fundsdata_arr; clear variables_arr;
[header_main, date_vec, fundsdata_arr] = read_mainmatrix('matrix_fundreturns.txt');
clear header_bench; clear datebench_vec; clear benchdata_arr;
NumFunds = length(fundsdata_arr);
[header_bench, datebench_vec, benchdata_arr] = read_mainmatrix('matrix_benchmarkreturns.txt');
if (length(benchdata_arr) ~= length(date_vec))
    Error('Benchmark set is wrong');
end

fundsdata_arr = fundsdata_arr(:, 2:end);
benchdata_arr = benchdata_arr(:, 2:end);

% specify the structure and parameters for equally weighted and benchmark
% portfolios
clear portfolio;
portf_count = 1;
portfolio(portf_count).name = 'portfolioBM'; % benchmark portfolio
portfolio(portf_count).type = 'benchmark';
portfolio(portf_count).lambda = 0; % no risk constraint
portfolio(portf_count).riskfunction = '';  % no risk constraint 
portfolio(portf_count).riskfunction_parameter = [];  % no risk constraint
portfolio(portf_count).modification = @(x) x;  % no risk constraint  
portf_count = portf_count+1;
portfolio(portf_count).name = 'portfolioEW'; % predefined equal weighted portfolio
portfolio(portf_count).type = 'ew';
portfolio(portf_count).lambda = 0; % not risk constrained
portfolio(portf_count).riskfunction = '';  % not risk constrained
portfolio(portf_count).riskfunction_parameter = [];  % not risk constrained
portfolio(portf_count).modification = @(x) x;  % not risk constrained
portf_count = portf_count+1;

% specify the structure and parameters for portfolios with risk measures
if strcmp(TYPE, 'RISK')
    
% Specify parameters for portfolios with some risk measures (which is NOT vriance). 
% You should specify several values of the penalty coefficient in the 
% lambda_list to have several trajectories with these coefficients.
% Currently only three values of the penalty coefficient are specified. 
% According to these values, three optimal portfolios are generated with names "portfolio1", "portfolio2", 
% and "portfolio3".
    lambda_list = [DefaultLambda 1.1*DefaultLambda  1.2*DefaultLambda 1.3*DefaultLambda 1.5*DefaultLambda 1.6*DefaultLambda]; % risk penalty coefficients
    for iportf=1:length(lambda_list)
        portfolio(portf_count).name = ['portfolio' num2str(iportf)]; % portfolio name like portfolio1, portfolio2, ... 
        portfolio(portf_count).type = 'risk';
        portfolio(portf_count).lambda = lambda_list(iportf); % set your value of risk bound
        portfolio(portf_count).riskfunction = 'meanabs_dev';  % set your risk measure  
        %portfolio(portf_count).riskfunction_parameter = 0;  % set parameter of your risk measure  
        portfolio(portf_count).riskfunction_parameter = [];        
        portfolio(portf_count).modification = @(x) x;  % 
        portf_count = portf_count+1;
    end
else

% Specify parameters for portfolios with Variance measure. 
% You should specify several values of the penalty coefficient in the 
% lambda_list to have several trajectories with these coefficients.
% Currently only three values of the penalty coefficient are specified. 
% According to these values, three optimal portfolios are generated with names "portfolio1", "portfolio2", 
% and "portfolio3".
    lambda_list = [0.1*DefaultLambda DefaultLambda 10*DefaultLambda]; % % variance penalty coefficients
    for iportf=1:length(lambda_list)
        portfolio(portf_count).name = ['portfolio' num2str(iportf)]; % portfolio names like portfolio1, portfolio2, ... 
        portfolio(portf_count).type = 'risk';
        portfolio(portf_count).lambda = NumFunds*lambda_list(iportf); % set your value of risk bound
        portfolio(portf_count).riskfunction = 'variance';  % set your risk measure
        portfolio(portf_count).riskfunction_parameter = [];  % there is no parameter in variance 
        portfolio(portf_count).modification = @(x) sqrt(x);  % 
        portf_count = portf_count+1;
    end    
end

% Specify settings for the first pass
pass_count = 1;
outofsampleindex = FirstInsampleIndex+InsampleSize-1;
for iportf=1:length(portfolio)
    portfolio(iportf).returnportfolio(pass_count) = 0.0;
    portfolio(iportf).portfoliovalue(pass_count) = InitialPortfolioValue;
    portfolio(iportf).date(pass_count) = int32(str2num(cell2mat(date_vec{outofsampleindex})));
end
if SaveSolutionFile 
    outfilename = 'solution_HedgeFunds_riskprog.txt';
    fsolid = fopen(outfilename, 'w');    
end

pass_count = pass_count+1;
outofsampleindex = outofsampleindex+1;

mpsg_suppress_message('On'); % PSG messages are swithced off
% Begin generating several passes 
while outofsampleindex <= length(date_vec) % start main in sample - out of sample cycle
    fprintf('\nPass N%8d\n', pass_count);
    
% During each pass the program specifies initial month of the in-sample
% period, FirstInsampleIndex, and its final month, outofsampleindex-1.
% The outofsample index specifies the out of sample month used for testing 
% of the performance of the rebalanced portfolio.
% For in-sample period the program generates the following  data required 
% for solving the optimization problem:    
      
% matrix_fundsreturns containing scenarios of hedge funds returns needed
% for calculating Risk function included in the objective with the penalty
% coefficient    
    clear matrix_fundsreturns;
    matrix_fundsreturns = fundsdata_arr(FirstInsampleIndex:outofsampleindex-1, : );
    
% matrix_expectedreturns needed for calculating expected return of
% portfolio included in the objective      
    clear matrix_expectedreturns;
    matrix_expectedreturns = sum(matrix_fundsreturns)/size(matrix_fundsreturns, 1);

% matrix_budget containing unit coefficients for budget constraint    
    clear matrix_budget;
    matrix_budget = ones(1, size(matrix_fundsreturns, 2));
    
% Calculate matrix_beta containing coefficients for market-neutrality
% constraint    
    clear matrix_beta;
    matrix_funds = fundsdata_arr(FirstInsampleIndex:outofsampleindex-1,:);
    middle_funds = sum(matrix_funds)/size(matrix_funds, 1);
    for i=1:size(matrix_funds, 1)
        matrix_funds(i,:)=matrix_funds(i,:)-middle_funds;
    end
    vector_bench = benchdata_arr(FirstInsampleIndex:outofsampleindex-1);
    vector_bench = vector_bench - sum(vector_bench)/size(vector_bench, 1);
    bench_norm = vector_bench'*vector_bench;
    vector_bench = vector_bench/bench_norm;
    matrix_beta = (matrix_funds'*vector_bench)';
  
    portf_count = 1;
   
 % start portfolio cycle    
    while portf_count <= length(portfolio) 
        current_portfoliovalue = portfolio(portf_count).portfoliovalue(pass_count-1);
        
% generate problem statement for portfolio with risk measure        
        if (strcmp(portfolio(portf_count).type ,'risk')) 
            risk = portfolio(portf_count).riskfunction;
            w = portfolio(portf_count).riskfunction_parameter;
            H = portfolio(portf_count).modification(portfolio(portf_count).lambda)/current_portfoliovalue*matrix_fundsreturns;
            c = []; p = [];
            d = -1.0/current_portfoliovalue*matrix_expectedreturns';
            A = [matrix_beta; -matrix_beta];
            b = [K*current_portfoliovalue; K*current_portfoliovalue];
            Aeq = matrix_budget;
            beq = BudgetBound*current_portfoliovalue;
            lb = current_portfoliovalue*VariablesLowerBound;
            ub = current_portfoliovalue*VariablesUpperBound;
            x0 = [];
            if strcmp(TYPE, 'RISK')
                options.Linearization = 'On';
                options.Solver = 'TANK';
            else
                options.Linearization = 'Off';
                options.Solver = 'VAN';
            end
            options.Precision = Precision; options.Stages = Stages;
% optimize problem  for portfolio with risk measure  
            [optimal_point, fval, solutionstatus, output] = riskprog (risk, w, H, c, p, d, A, b, Aeq, beq, lb, ub, x0, options);
            
% for portfolio with risk measure, determine portfolio return and value corresponding to out of sample month 
            portfolio_return = fundsdata_arr(outofsampleindex, :)*optimal_point;
            portfolio(portf_count).portfoliovalue(pass_count) = current_portfoliovalue+portfolio_return;
            portfolio(portf_count).returnportfolio(pass_count) = portfolio_return;
            portfolio(portf_count).date(pass_count) = int32(str2num(cell2mat(date_vec{outofsampleindex})));
            
% save results for portfolio with risk measure in external file 
            if SaveSolutionFile 
                fprintf(fsolid, 'Portfolio: %s, end date: %s, status: %s\n', portfolio(portf_count).name, ...
                    cell2mat(date_vec{outofsampleindex-1}), solutionstatus);
                fprintf(fsolid, 'Objective: %.6f\n', fval);
                fprintf(fsolid, 'Budjet Constraint: %.6f\nBetaneutrality Constraint: %.6f\n', output.fAeqval, output.fAval(1));
                fprintf(fsolid, '****************************************\n\n');
            end
        else 
            
% Determine benchmark portfolio return and value corresponding to out of sample month            
            if strcmp(portfolio(portf_count).type, 'benchmark') 
                portfolio_return = benchdata_arr(outofsampleindex)*current_portfoliovalue;
                portfolio(portf_count).returnportfolio(pass_count) = portfolio_return;
                portfolio(portf_count).portfoliovalue(pass_count) = current_portfoliovalue+portfolio_return;
                portfolio(portf_count).date(pass_count) = int32(str2num(cell2mat(date_vec{outofsampleindex})));
            else
                
                if strcmp(portfolio(portf_count).type, 'ew')

% Create equally weighted portfolio by including the best portfolio funds 
% having maximum returns for in sample period.  
% The number of the best funds in the portfolio is BestFundsNumber.
                    infundreturns = sum(matrix_fundsreturns);
                    outfundreturns = fundsdata_arr(outofsampleindex, :);
% Determine equally weighted portfolio return and value corresponding to out of sample month                     
                    [locarr, indexarr] = sort(infundreturns, 'descend');
                    portfolio_return = sum(outfundreturns(indexarr(1:BestFundsNumber)))*current_portfoliovalue/BestFundsNumber;
                    portfolio(portf_count).returnportfolio(pass_count) = portfolio_return;
                    portfolio(portf_count).portfoliovalue(pass_count) = current_portfoliovalue+portfolio_return;
                    portfolio(portf_count).date(pass_count) = int32(str2num(cell2mat(date_vec{outofsampleindex})));
                end
            end
        end
        portf_count = portf_count+1;
    end % end portfolio cycle

% Specify settings for the next pass    
    outofsampleindex = outofsampleindex+1;
    pass_count = pass_count+1;
end % end main in sample - out of sample cycle
mpsg_suppress_message('Off'); % PSG messages are swithced on


if SaveSolutionFile 
    fclose(fsolid);
end


% save tables with results of calculations
outfilename = 'CS_HedgeFunds_Output_riskprog.txt';
fid = fopen(outfilename, 'w');
fprintf('Table 1. Portfolios\n%8s', 'Date');
fprintf(fid, 'Table 1. Portfolios\n%8s', 'Date');
for i=1:size(portfolio, 2)
    fprintf('%15s', portfolio(i).name);
    fprintf(fid, '\t%15s', portfolio(i).name);
end
fprintf('\n');
fprintf(fid, '\n');
for i=1:size(portfolio(1).portfoliovalue, 2)
    fprintf('%8d', portfolio(1).date(i));
    fprintf(fid, '%8d', portfolio(1).date(i));
    for ii=1:size(portfolio, 2)
        fprintf('%15.4f', portfolio(ii).portfoliovalue(i)/InitialPortfolioValue);
        fprintf(fid, '\t%15.4f', portfolio(ii).portfoliovalue(i)/InitialPortfolioValue);
    end
    fprintf('\n');
    fprintf(fid, '\n');
end
fprintf(fid, '\nTable 2. Monthly Portfolios Returns\n%8s', 'Date');
for i=1:size(portfolio, 2)
    fprintf(fid, '\t%15s', portfolio(i).name);
end
fprintf(fid, '\n');
for i=1:size(portfolio(1).portfoliovalue, 2)
    fprintf(fid, '%8d', portfolio(1).date(i));
    for ii=1:size(portfolio, 2)
        fprintf(fid, '\t%15.4f', portfolio(ii).returnportfolio(i)/InitialPortfolioValue);
    end
    fprintf(fid, '\n');
end
for i=1:size(portfolio(1).portfoliovalue, 2)
    fprintf(fid, '%8d', portfolio(1).date(i));
    for ii=1:size(portfolio, 2)
        fprintf(fid, '\t%15.4f', portfolio(ii).returnportfolio(i)/InitialPortfolioValue);
    end
    fprintf(fid, '\n');
end
date_vec = (portfolio(1).date)./100;
pos = 1;
year_count = 0;
year_vec = int32([]);
while pos<length(date_vec)
    year = date_vec(pos);
    num_months = sum(date_vec==year);
    if num_months == 12
        year_count = year_count+1;
        year_vec(year_count) = year;
    end
    pos = pos+num_months;
end
fprintf('\nTable 3. Actual Performance of Portfolios\n');
fprintf(fid, '\nTable 3. Actual Performance of Portfolios\n');
fprintf('%15s', 'Portfolio'); fprintf('%9d(%%)', year_vec); 
fprintf('%15s%15s%15s%15s\n', 'Annual(%)', 'Total(%)','Max_DD(%)','Sharpe');
fprintf(fid, '%15s', 'Portfolio'); fprintf(fid, '\t%9d(%%)', year_vec); 
fprintf(fid, '\t%15s\t%15s\t%15s\t%15s\n', 'Annual(%)', 'Total(%)','Max_DD(%)','Sharpe');
for iportf = 1: length(portfolio)
    [year_annual_vec, inception_annual, total, max_dd, sharpe_ratio] = get_portfolio_statistics(year_vec, portfolio(iportf));
    fprintf('%15s', portfolio(iportf).name); fprintf('%12.4f', year_annual_vec*100); 
    fprintf('%15.4f%15.4f%15.4f%15.4f\n', inception_annual*100, total*100, -max_dd*100, sharpe_ratio);
    fprintf(fid, '%15s', portfolio(iportf).name); fprintf(fid, '\t%12.4f', year_annual_vec*100); 
    fprintf(fid, '\t%15.4f\t%15.4f\t%15.4f\t%15.4f\n', inception_annual*100, total*100, -max_dd*100, sharpe_ratio);
end
fclose(fid);

% display Graphs comparing performance of portfolios
YData = zeros(length(portfolio(1).portfoliovalue), length(portfolio));
legent_vec = {};
for iportf = 1: length(portfolio)
    YData(:, iportf) = portfolio(iportf).portfoliovalue./InitialPortfolioValue;
    legent_vec(iportf) = cellstr(portfolio(iportf).name);
end
thisplot = plot(YData);
set(thisplot, 'LineWidth', 2);
set(thisplot, 'MarkerSize', 5);
set(gca,'XTick',1:4:length(portfolio(1).date));
hAxes = gca;
set(hAxes, 'XTickLabel', portfolio(1).date(1:4:length(portfolio(1).date)));
grid on;
xlabel('Month');
ylabel('Portfolio Value');
legend(legent_vec, 'Location', 'NorthWest');
end % end of main function

% utility functions (for this main function only)
function [year_annual_vec, inception_annual, total, max_dd, sharpe] = get_portfolio_statistics(year_vec, portfolio)
    year_annual_vec = [];
    RISKFREE = 0.03;
    for iyear=1:length(year_vec)
        startdate = int32(year_vec(iyear)*100+1);
        enddate = int32(year_vec(iyear)*100+12);
        year_annual_vec(iyear) = 0.0;
        for imonth=1:length(portfolio.date)
            if portfolio.date(imonth) == startdate 
                furstval = portfolio.portfoliovalue(imonth);
            end
            if portfolio.date(imonth) == enddate
                lastval = portfolio.portfoliovalue(imonth);
            end
        end
        year_annual_vec(iyear) = lastval/furstval-1;
    end
    total = portfolio.portfoliovalue(end)/portfolio.portfoliovalue(1)-1.0;
    month_returns_arr = portfolio.returnportfolio./portfolio.portfoliovalue; month_returns_arr = month_returns_arr(2:end);
    month_return = sum(month_returns_arr)/length(month_returns_arr);
    month_deviation = sqrt((month_returns_arr-month_return)*(month_returns_arr-month_return)'/length(month_returns_arr));
    sharpe = sqrt(12)*(month_return-RISKFREE/12)/month_deviation;
    inception_annual = month_return*12;
    dd_arr = zeros(1, length(portfolio.portfoliovalue));
    maxvalue = portfolio.portfoliovalue(1);
    for imonth=2:length(portfolio.portfoliovalue)
        maxvalue = max(maxvalue, portfolio.portfoliovalue(imonth));
        dd_arr(imonth) = (maxvalue-portfolio.portfoliovalue(imonth))/maxvalue;
    end
	max_dd = max(dd_arr);    
end

%=======================================================================|
%American Optimal Decisions, Inc. Copyright                             |
%Copyright ©American Optimal Decisions, Inc. 2007-2015.                 |
%American Optimal Decisions (AOD) retains copyrights to this material.  |
%                                                                       |
%Permission to reproduce this document and to prepare derivative works  |
%from this document for internal use is granted, provided the copyright |
%and “No Warranty”  statements are included with all reproductions       |
%and derivative works.                                                  |
%                                                                       |
%For information regarding external or commercial use of copyrighted    |
%materials owned by AOD, contact AOD at support@aorda.com.              |
%=======================================================================|