%% KRG Portfolio Liquidity Optimization Example.
% The user runs an optimization strategy to calculate the portfolio value,
% target number of shares held, shares to be traded and market impact cost
% of trades.

%   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);

% Load the example data
load KRGExampleData.mat TradeDataPortOpt CovarianceData
n = 10;
TradeDataPortOpt = TradeDataPortOpt(1:n,:);
CovarianceData = CovarianceData(1:n,1:n);

% Value to withdraw from the portfolio
PortLiquidationValue = 100000000; 

% Portfolio risk needs to remain within 90% and 110% of its current value
PortRiskBounds = [0.9 1.10]; 

% Maximum market impact cost for all trades
maxTotalMI = 0.0050;              

% Number of stocks 
numPortStocks = length(TradeDataPortOpt.Symbol);

% Convert Covariance Table to Covariance Matrix for Stock Sample Set J
C = table2array(CovarianceData);

% Individual Market Impact Cost Constraint
maxMI = TradeDataPortOpt.UB_MaxMI;

% Portfolio values
PortfolioValue = sum(TradeDataPortOpt.Value);
absPortValue = abs(TradeDataPortOpt.Value);
PortfolioAbsValue = sum(absPortValue);
PortfolioTargetValue = PortfolioValue - PortLiquidationValue;

% Current portfolio weight
w = sign(TradeDataPortOpt.Shares) .* absPortValue / PortfolioAbsValue;

% Sum of weights = 1
Aeq = ones(1,numPortStocks);
beq = 1;

A = [];
b = [];

% Bounded constraints
LB = TradeDataPortOpt.LB_Wt;
UB = TradeDataPortOpt.UB_Wt;

% User enters constraint on Shares, Value, PctADV (Optional)
% If any of these constraints are selected, the constraints can be
% simplified into one major constraint based on SharesToTrade
lbShares = max([TradeDataPortOpt.LB_MinShares,TradeDataPortOpt.LB_MinValue./TradeDataPortOpt.Price,TradeDataPortOpt.LB_MinPctADV.*TradeDataPortOpt.ADV],[],2);
ubShares = min([TradeDataPortOpt.UB_MaxShares,TradeDataPortOpt.UB_MaxValue./TradeDataPortOpt.Price,TradeDataPortOpt.UB_MaxPctADV.*TradeDataPortOpt.ADV],[],2);

% Initial starting solution
x0 = TradeDataPortOpt.Value ./ sum(TradeDataPortOpt.Value);
x = x0;

% Optimization options
options = optimoptions('fmincon','Algorithm','sqp',...
        'TolFun',10E-8,'TolX',10E-16,'TolCon',10E-8,'TolPCG',10E-8,...
        'MaxFunEvals',50000,'MaxIter',50000,'DiffMinChange',10E-8);

% MATLAB nonlinear optimization setup
% 
% 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)


% Liquidation Non-Linear Optimization Problem
x = fmincon(@(x) krgLiquidityFunction(x,TradeDataPortOpt,PortfolioTargetValue,k),x0,A,b,Aeq,beq,LB,UB,...
    @(x) krgLiquidityConstraint(x,w,C,TradeDataPortOpt,PortfolioTargetValue,PortRiskBounds,lbShares,ubShares,maxMI,maxTotalMI,k),...
    options);

% Optimized values

% Portfolio value 
x1 = x .* PortfolioTargetValue / PortfolioValue;

% Number of shares held
TargetShares = x * PortfolioTargetValue ./ TradeDataPortOpt.Price;

% Number of shares to be traded to meet Number of shares held value
SharesToTrade = TradeDataPortOpt.Shares - TargetShares;

% Portfolio value denoted in currency
TargetValue = x * PortfolioTargetValue;
    
% Calculate shares from weights, portfolio value, and stock price.
% Market Impact cost is always positive and based on positive shares 
% quantities
TradeDataPortOpt.Shares = abs(SharesToTrade);

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

% marketImpact calculation, returned as decimal values
MI = marketImpact(k,TradeDataPortOpt)/10000;
