function [Weight,Shares,Value,MI] = KRGPortfolioOptimizationExample(TradeData,CovarianceData,MIOptScenario)
%KRGPortfolioOptimizationExample Optimize portfolio strategies example.
%   [WEIGHT,SHARES,VALUE,MI] = KRGPORTFOLIOOPTIMIZATIONEXAMPLE(T,C,M) runs
%   one of three optimization strategies given the portfolio data, T,
%   covariance data, C, and optimization scenario flag, M.   M is set to 1
%   by default if none is specified.
%
%   For M = 1, the strategy maximizes the return minus lambda multiplied by
%   risk.   
%
%   For M = 2, the strategy minimizes risk subject to net returns greater
%   than or equal to the minimum return.
%
%   For M = 3, the strategy maximizes the net return subject to the risk 
%   less than or equal to the maximum risk.
%
%   For each chosen scenario, the weights, WEIGHT, number of 
%   shares, SHARES, values, VALUE, and market impact costs, MI, are 
%   returned for each member of the portfolio.
%
%   For example, 
%
%   [Weight,Shares,Value,MI] = KRGPortfolioOptimizationExample(TradeDataPortOpt,CovarianceData,1)
%
%   returns
%
%   Weight =
%
%            0.0100
%            0.0100
%            0.0100
%            0.0100
%            0.0100 ...
%
%   Shares =
%
%         24420.02
%         101630.15
%         24995.81
%         5740.20
%         8281.34 ...
%
%   Value =
%
%            1000000
%            1000000
%            1000000
%            1000000
%            1000000 ...
%
%   ans =
%   
%      1.0e-03 *
%   
%       0.1250
%       0.0993
%       0.0758
%       0.1039
%       0.0203 ...
   
%   Copyright 2016-2020 The MathWorks, Inc.

% Load the market impact data provided by the Kissell Research Group
f = ftp('ftp.kissellresearch.com','username','password');
cd(f,'MI_Parameters');
mget(f,'MI_Encrypted_Parameters.csv');

% Create KRG object
miData = readtable('MI_Encrypted_Parameters.csv','Format','auto','delimiter',',','ReadRowNames',false,'ReadVariableNames',true);
k = krg(miData,datetime('today'),1,250);

% Default to first optimization scenario if none specified
if nargin < 3 || isempty(MIOptScenario)
  MIOptScenario = 1;
end

% Specify portfolio optimization parameters 
PortfolioDollars = 100000000;
minReturn = 0.10;
maxRisk = 0.15;
maxTotalMI = 0.0050; 
Lambda = 1;

% Convert covariance data table to array.    
C = table2array(CovarianceData);

% Sum of portfolio weights must sum to 1
numStocks = length(TradeData.Symbol);
Aeq = ones(1,numStocks);
beq = 1;

% Initialize Constraint Matrices
A = [];
b = [];

% Important Note:
%
% 1) The constraints for LB/UB, Min/Max Shares, Size, and Value are all closely related.
% 2) Users may wish to specify constraints in one or more of these values.
% The Market Impact constraints Need to be specified as an extra
% Non-Linear Constraint. 
% Simplify these constraints by converting each to a Dollar Weight
% based on Portfolio Value.
%
% This is as follows:

% LB_Wt
if isempty(TradeData.LB_Wt) || any(isnan(TradeData.LB_Wt))
    TradeData.LB_Wt = -inf(numStocks,1);
end

% UB_Wt
if isempty(TradeData.UB_Wt) || any(isnan(TradeData.UB_Wt))
    TradeData.UB_Wt = inf(numStocks,1);
end

% LB_MinShares
if isempty(TradeData.LB_MinShares) || any(isnan(TradeData.LB_MinShares))
    TradeData.LB_MinShares = -inf(numStocks,1);
else
    TradeData.LB_MinShares = TradeData.LB_MinShares .* TradeData.Price ./ PortfolioDollars;
end

% UB_MaxShares
if isempty(TradeData.UB_MaxShares) || any(isnan(TradeData.UB_MaxShares))
    TradeData.UB_MaxShares = inf(numStocks,1);
else
    TradeData.UB_MaxShares = TradeData.UB_MaxShares .* TradeData.Price ./ PortfolioDollars;
end

% LB_MinPctADV
if isempty(TradeData.LB_MinPctADV) || any(isnan(TradeData.LB_MinPctADV))
    TradeData.LB_MinPctADV = -inf(numStocks,1);
else
    TradeData.LB_MinPctADV = TradeData.LB_MinPctADV .* TradeData.ADV .* TradeData.Price ./ PortfolioDollars;
end

% UB_MaxPctADV
if isempty(TradeData.UB_MaxPctADV) || any(isnan(TradeData.UB_MaxPctADV))
    TradeData.UB_MaxPctADV = inf(numStocks,1);
else
    TradeData.UB_MaxPctADV = TradeData.UB_MaxPctADV .* TradeData.ADV .* TradeData.Price ./ PortfolioDollars;
end

% LB_MinValue
if isempty(TradeData.LB_MinValue) || any(isnan(TradeData.LB_MinValue))
    TradeData.LB_MinValue = -inf(numStocks,1);
else
    TradeData.LB_MinValue = TradeData.LB_MinValue ./ PortfolioDollars;
end

% UB_MaxValue
if isempty(TradeData.UB_MaxValue) || any(isnan(TradeData.UB_MaxValue))
    TradeData.UB_MaxValue = inf(numStocks,1);
else
    TradeData.UB_MaxValue = TradeData.UB_MaxValue ./ PortfolioDollars;
end

% Convert to a Single LB & UB
LB = max([TradeData.LB_Wt, TradeData.LB_MinShares, TradeData.LB_MinPctADV, TradeData.LB_MinValue]')'; %#ok
UB = min([TradeData.UB_Wt, TradeData.UB_MaxShares, TradeData.UB_MaxPctADV, TradeData.UB_MaxValue]')'; %#ok

% specify MI constraints
maxMI = TradeData.UB_MaxMI;

% if there are no constraints specified then set LB to empty matrix LB=[]
if max(LB) == -inf
    LB = [];
end

% if there are no constraints specified then set UB to empty matrix UB=[]
if min(UB) == inf
    UB = [];
end

% Initial Starting Solution
x0 = 1 ./ numStocks(ones(numStocks,1));

% Optimization Options Set
% These can be changed and improved for our objective function
%options = optimoptions('fmincon','Algorithm','interior-point');
options = optimoptions('fmincon','Algorithm','sqp',...
        'TolFun',10E-8,'TolX',10E-16,'TolCon',10E-8,'TolPCG',10E-8,...
        'MaxFunEvals',100000,'MaxIter',100000,'DiffMinChange',10E-8);
      
% MATLAB OPTIMIZATION FUNCTIONS
%
% fmincon function
%
%      fmincon attempts to solve problems of the form:
%      min F(X)  subject to:  A*X  <= B, Aeq*X  = Beq (linear constraints)
%       X                     C(X) <= 0, Ceq(X) = 0   (nonlinear constraints)
%                                LB <= X <= UB        (bounds)

% optimization scenarios
switch MIOptScenario
  
  case 1

    % Scenario#1: Maximize Returns - Lambda*Risk

    % Non-Linear Optimization
    %
    % NetReturns = Returns - MI
    %
    % Max NetReturns-LambdaRisk
    %
    % Max (R'*x - MI'*x) - Lambda*Risk
    %      fmincon attempts to solve problems of the form:
    %      min F(X)  subject to:  A*X  <= B, Aeq*X  = Beq (linear constraints)
    %       X                     C(X) <= 0, Ceq(X) = 0   (nonlinear constraints)
    %                                LB <= X <= UB        (bounds)

    x = fmincon(@(x) krgPortfolioOptimizer(x,Lambda,C,TradeData,PortfolioDollars,k,MIOptScenario),x0,A,b,Aeq,beq,LB,UB,@(x) krgNonLinearConstraint(x,C,TradeData,PortfolioDollars,minReturn,k,MIOptScenario,maxMI,maxTotalMI),options);

  case 2

    % Scenario#2: Minimize Risk subject to: Net Returns >= minReturn

    % Non-Linear Optimization
    %
    % Min:  f=x'*C*x'
    %       s.t. Aeq*x=1
    %            Net Return>= R* (Net Return is NonLinear Function using MI formula)  
    x = fmincon(@(x) krgPortfolioOptimizer(x,[],C,TradeData,PortfolioDollars,k,MIOptScenario),x0,A,b,Aeq,beq,LB,UB,@(x) krgNonLinearConstraint(x,C,TradeData,PortfolioDollars,minReturn,k,MIOptScenario,maxMI,maxTotalMI),options);

  case 3

    % Scenario#3: Max NetReturn s.t. Risk <= maxRisk
     % Non-Linear Optimization
    %
    % Min:  f=x'*C*x'
    %       s.t. Aeq*x=1
    %            Net Return>= R* (Net Return is NonLinear Function using MI formula)  
    x = fmincon(@(x) krgPortfolioOptimizer(x,[],C,TradeData,PortfolioDollars,k,MIOptScenario),x0,A,b,Aeq,beq,LB,UB,@(x) krgNonLinearConstraint(x,C,TradeData,PortfolioDollars,maxRisk,k,MIOptScenario,maxMI,maxTotalMI),options);

end

% Portfolio weight in decimal, shares and dollars
Weight = x;
Shares = x .* PortfolioDollars ./ TradeData.Price;
Value = x .* PortfolioDollars;
    
% Calculate Shares from weights, portfolio value, and stock price
% MI Cost is always positive and based on Positive Shares Quantities
% need to take absolute value
TradeData.Shares = abs(Shares);

% ADV factor in TradeTime is different than tradetime2pov calc.
TradeData.TradeTime = TradeData.TradeTime .* TradeData.ADV;
TradeData.POV = krg.tradetime2pov(TradeData.TradeTime,TradeData.Shares);

% marketImpact calculation, convert from bp to decimal
MI = marketImpact(k,TradeData)/10000;