%The case study demonstrates optimization setting 
%for a portfolio replication problem with the
%replication error measured by Mean Absolute Penalty.
%Underperformance of the portfolio compared 
%to S&P100 index is measures by CVaR.
%Distribution of residuals is shaped with a CVaR constraint
%(several constraints can be specified, if of interest).

clc
clear;

%Define point_prices values for stocks 1- 30
point_prices = [19.78129959 42.5625 73.625 64 12.6875 47.0625...
    38.85430145 34.8125 68.125 73 50.1875 20.5625...
    8.3125 38 26.1406002 59.1875 32.625 17.6875 8.875...
    69.875 29.90629959 17.3125 19.5 42.21879959 11.1875...
    14.7031002 19.15629959 19.92000008 26.25 31.8125]; 

%Specify a set of values of parameter (upper bound on risk), for which runs should be conducted:
Parameters = [0.0002 0.0005 0.001 0.003 0.005 0.01 0.02];

%Load the problem:
load('problem_CS_Portfolio_Replication_CVaR_0p001_data.mat');

warning('off', 'all');
%Set solver options:
options.Solver = 'VAN';
options.Precision  = 3; 
 
%Set graph options:
stroptions.PlotGraph = 'On';
stroptions.PlotType = 1;

%Set the Scale parameter
stroptions.ScaleObjective = -10000;
stroptions.ScaleParam = 10000;

%Solve the optimization problem for specified  values of the parameter
[Objectives, Parameters, Points, GraphHandle]=riskconstrparam (risk1, risk2, [], H1, c1, p1, w2, H2, c2, p2,...
[], r, [], [], Aeq, beq, lb, [], 'r', [],[],[],...
[], [], Parameters, [], stroptions);

%Sort components in lexicigraphical order
component_name_arr_sorted = sort (component_name_arr);

%Sort Points and Matrices in accordance with components order
for j=1:size(component_name_arr) 
    for i=1:size(component_name_arr) 
        k = strcmp(component_name_arr_sorted(i), component_name_arr(j));
        if k == true
                Points_sorted(:,i)= Points (:,j);
                H2_sorted (:,i) = H2 (:,j);
                H1_sorted (:,i) = H1 (:,j);
        end
    end
end
    
% Display results as table
fprintf('\n%30s','Bound on the underperformance');
for i=1:1:length(Parameters)
    fprintf('%12.6f', Parameters(i));
end

fprintf('\n%30s','Replication error');
for i=1:1:length(Parameters)
    fprintf('%12.4f', Objectives(i));
end

fprintf('\n%30s','Optimal value of components');
for j=1:size(component_name_arr)
    fprintf('\n%30s', component_name_arr_sorted {j});
   for i=1:1:length(Parameters)
    fprintf('%12f', Points_sorted(i,j));
   end
end 

disp('  ')

%--------------------------------------------------------------------------
% Create plot for meanabs_pen_Portfolio_Replication vs. cvar_risk_Portfolio_Replication

clear x;
clear y;
for i=1:7
  x(i) = functionvalue('cvar_risk', w2, H2, c2, p2, Points(i,:)')*10000;
  y(i) = functionvalue('meanabs_pen', [], H1, c1, p1, Points (i,:)')*10000;
end 
  figure
plot(...
  x,y,...
  'Color',[1 0 0],...
  'LineStyle','none',...
  'Marker','square',...
  'MarkerFaceColor',[1 0 0]);

xlabel('cvar\_risk\_Portfolio\_Replication');
ylabel('meanabs\_pen\_Portfolio\_Replication');
title('meanabs\_pen\_Portfolio\_Replication vs. cvar\_risk\_Portfolio\_Replication');
grid on;

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Create plot for relative contribution to CVaR

% Define increment of the var_risk_Portfolio_Replication for point
% point_problem_CS_Portfolio_Replication_CVaR_0p001

for i = 1:30
    x(i) = i;
end

[fval] = functionvalue('cvar_risk', w2, H2_sorted, c2, p2, Points_sorted(3,:)');
[Increments] = functionincrement('cvar_risk', w2, H2_sorted, c2, p2, Points_sorted (3,:)');
y = Increments/fval;

    figure
    plot(...
  x,y,...
  'Color',[1 0 0],...
  'LineStyle','none',...
  'Marker','square',...
  'MarkerFaceColor',[1 0 0]);
grid on;
set(gca,'XTick',x);
hAxes = gca;
set (hAxes, 'YLim',([-155 35]))
xlabel('# of component');
ylabel('relative contribution to CVaR');
title('relative contribution to CVaR');
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Create plot for relative contribution to meanabs_pen

% Define increment of the meanabs_pen_Portfolio_Replication for point
% point_problem_CS_Portfolio_Replication_CVaR_0p001

for i = 1:30
    x(i) = i;
  end

[fval] = functionvalue('meanabs_pen', [], H1_sorted, c1, p1, Points_sorted(3,:)');
[Increments] = functionincrement('meanabs_pen', [], H1_sorted, c1, p1, Points_sorted (3,:)');
y = Increments/fval;

    figure
    plot(...
  x,y,...
  'Color',[1 0 0],...
  'LineStyle','none',...
  'Marker','square',...
  'MarkerFaceColor',[1 0 0]);
grid on;
set(gca,'XTick',x);

hAxes = gca;
set (hAxes, 'YLim',([-8 2]))
xlabel('# of component');
ylabel('relative contribution to meanabs');
title('relative contribution to meanabs\_pen');
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
% Create plot for meanabs_pen_marginal vs. exposure

% Define increment of the function meanabs_pen_marginal for point
% point_problem_CS_Portfolio_Replication_CVaR_0p001
format short e;
[Increments] = functionincrement('meanabs_pen', [], H1_sorted, c1, p1, Points_sorted (3,:)');
    clear x;
    clear y;
    j=1;
    for i=1:30
       if abs(Points_sorted (3,i)) > 1.0E-6
           x(j) = Points_sorted(3,i)*point_prices(i);
           y(j) = 10000*Increments(i)/x(j);
           j = j+1;
       end
    end
  figure
  plot(...
  x,y,...
  'Color',[1 0 0],...
  'LineStyle','none',...
  'Marker','square',...
  'MarkerFaceColor',[1 0 0]);
   grid on;
   xlabel('exposure');
   ylabel('marginal');
   title('meanabs\_pen\_marginal vs. exposure');
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Create plot for CVaR_marginal vs. exposure

% Define increment of the function CVaR_pen_marginal for point
% point_problem_CS_Portfolio_Replication_CVaR_0p001
format short e;
[Increments] = functionincrement('cvar_risk', w2, H2_sorted, c2, p2, Points_sorted (3,:)');
    clear x;
    clear y;
    j=1;
    for i=1:30
       if abs(Points_sorted (3,i)) > 1.0E-6
           x(j) = Points_sorted(3,i)*point_prices(i);
           y(j) = 10000*Increments(i)/x(j);
           j = j+1;
       end
    end
  figure
  plot(...
  x,y,...
  'Color',[1 0 0],...
  'LineStyle','none',...
  'Marker','square',...
  'MarkerFaceColor',[1 0 0]);
   grid on;
   xlabel('exposure');
   ylabel('marginal');
   title('CVaR\_marginal vs. exposure');
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Create histogram for Probability Density Function of Losses

loss = 10000*(c1 - H1_sorted*Points_sorted(3,:)');
figure
[bins,xout] = hist(loss,40);
bins = bins/sum(bins);
bar(xout,bins,'r');
grid on;
xlabel('Losses (USD)');
ylabel ('Probability');
tit =sprintf('%s%','Probability Density Function of Losses: matrix\_portfolio\_replication\_scenarios');
title(tit);

%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Create pie diagram for exposure of the point
% point_problem_CS_Portfolio_Replication_CVaR_0p001
% The graph shows n_plot maximal values of the exposures. 
% Exposures with less values are combined into category "Others".
n_plot = 10;
clear u;
u = sort(abs(Points_sorted(3,:).*point_prices), 'descend');
i = 0;
other = 0;
for j=1:1:size(Points_sorted,2)
    if abs(Points_sorted(3,j)*point_prices(j)) >= u(n_plot)
         i = i + 1;
         point_variables(i) = component_name_arr_sorted(j);
         Point1(i) = Points_sorted(3,j)*point_prices(j);
    else
        other = other + Points_sorted(3,j)*point_prices(j);
    end
end

Point1(n_plot+1) = other;
point_variables(n_plot+1)={'other'};

figure;
pie(Point1);
legend(point_variables, 'Location','EastOutside');
title(sprintf('%0.f%s',n_plot,' Largest Exposures of  the Point point\_problem\_CS\_Portfolio\_Replication\_CVaR\_0p001'));

%=======================================================================|
%American Optimal Decisions, Inc. Copyright                             |
%Copyright ©American Optimal Decisions, Inc. 2007-2014.                 |
%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.              |
%=======================================================================|