%% KRG Trade Schedule Optimization.
% To Compute Market Impact Cost and Optimal Trading Strategies:
% Determine the optimal order slicing strategy for the basket or program.
% Optimize trading schedules based on risk aversion parameters, and trading needs.
% Perform real-time re-optimization based on changing conditions.
% Uses Trademetrics Optimization process
% Based on techniques presented in:
%   The Science of Algorithmic Trading and Portfolio Management, Robert Kissell 2013
%   Optimal Trading Strategies, Robert Kissell & Mort Glantz, 2003
%   Multi-Period Trade Schedule Optimization, Roberto Malamut & Robert Kissell, 2002

%   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');
close(f)

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

% Import data
load KRGExampleData TradeDataTradeOpt CovarianceTradeOpt

% Convert table to array
CovarianceTradeOpt = table2array(CovarianceTradeOpt);

% Set the example optimization type
% For trading schedule optimization, optType = 1
% For POV optimization, optType = 2
optType = 1;

% Use total trading time of 1 day with 13 trading periods 
totalDays = 1;
periodsPerDay = 13;

% Set risk aversion level
Lambda = 0.5;

% Set minimum and maximum percentage of volume
minPOV = 0.00;
maxPOV = 0.60;

% total number of trading periods
totalNumberPeriods = totalDays * periodsPerDay;

% Portfolio Value
PortfolioValue = TradeDataTradeOpt.Price'*TradeDataTradeOpt.Shares;

% Number of stocks
numberStocks = height(TradeDataTradeOpt);

% Covariance matrix is annualized covariance matrix in decimals.
% Convert to ($/Shares)^2 units for the trade period, this matrix is for a
% two-sided portfolio, buys and sells or long and short.
CC = diag(TradeDataTradeOpt.Price) * CovarianceTradeOpt * diag(TradeDataTradeOpt.Price);        

% Scale to one trading period
CC = CC / periodsPerDay / k.TradeDaysInYear;       

% Convert to buy-only problem (e.g., one-sided problem)
CC = TradeDataTradeOpt.SideIndicator * TradeDataTradeOpt.SideIndicator' .* CC;      


% Convert Alpha_bp from basis points per day to cents/share per period
TradeDataTradeOpt.Alpha_bp = TradeDataTradeOpt.Alpha_bp / 1000 .* TradeDataTradeOpt.Price / totalNumberPeriods;

% Run trade schedule or POV optimization depending on optType flag

% Set the initial trade schedule or POV values
theta0 = rand(numberStocks,1);

% Optimization settings for fmincon
optionsold = optimset;
options = optimset(optionsold,'LargeScale','on','GradObj','off','DerivativeCheck',...
    'off','FinDiffType','central','FinDiffRelStep',1E-12,...
    'TolFun',10E-5,'TolX',10E-12,'TolCon',10E-12,'TolPCG',10E-12,...
    'MaxFunEvals',20000,'MaxIter',20000,'DiffMinChange',10E-04);

% Calculate optimal trading schedule or percentage of volume
if optType == 1
  % Set lower and upper bounds for percentage of shares in the trade
  % schedule
  LB = zeros(numberStocks,1);
  UB = 100 * ones(numberStocks,1);
  [theta,fval,exitflag,output] = fmincon(@optimizeTradingSchedule,theta0,[],[],[],[],LB,UB,[],options,totalNumberPeriods,numberStocks,periodsPerDay,TradeDataTradeOpt,CC,Lambda,k); 
elseif optType == 2
  % Set the upper and lower bounds of the POV values
  UB_POV = maxPOV.*ones(numberStocks,1);
  LB_POV = max(minPOV,TradeDataTradeOpt.Shares./(TradeDataTradeOpt.Shares+TradeDataTradeOpt.ADV.*totalDays));
  [POV,fval,exitflag,output] = fmincon(@optimizePOV,theta0,[],[],[],[],LB_POV,UB_POV,[],options,totalNumberPeriods,TradeDataTradeOpt,CC,Lambda,k); 
end

% Calculation Transaction Costs and Schedule

if optType == 1
    
  % Trade Schedule optimization
  numPeriods = 1:totalNumberPeriods;
  K_Matrix = repmat(numPeriods,numberStocks,1);
  Theta_Matrix = repmat(theta,1,totalNumberPeriods);
  Volume_Matrix = repmat(TradeDataTradeOpt.ADV/periodsPerDay,1,totalNumberPeriods);
  TradeDataTradeOpt.VolumeProfile = Volume_Matrix;
  Shares_Matrix = repmat(TradeDataTradeOpt.Shares,1,totalNumberPeriods);

  % X = Shares to trade in period i
  Xpct = (exp(-K_Matrix .* Theta_Matrix) .* (exp(Theta_Matrix)-1)) ./ (1 - exp(-totalNumberPeriods * Theta_Matrix));
  X = repmat(TradeDataTradeOpt.Shares,1,totalNumberPeriods) .* Xpct;
  TradeDataTradeOpt.TradeSchedule = X;
  
  % R = Residual Shares at beginning of period i
  Rpct = (exp(-(K_Matrix-1).*Theta_Matrix) - exp(-totalNumberPeriods.*Theta_Matrix)) ./ (1-exp(-totalNumberPeriods.*Theta_Matrix));
  R = repmat(TradeDataTradeOpt.Shares,1,totalNumberPeriods) .* Rpct;

  % Price Appreciations in Dollars
  PA = sum(R,2) .* TradeDataTradeOpt.Alpha_bp;
  
  % Market Impact in Dollars
  MI = marketImpact(k,TradeDataTradeOpt) .* TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Price ./10000;
  
  % Timing Risk in Dollars
  TR = sqrt(sum(R.^2,2) .* diag(CC));
  TR_bp = TR ./ (TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Price) * 10000;
    
  % Avg POV Rate
  kTR = ((TR_bp/10000*1./TradeDataTradeOpt.Volatility).^2).*(k.TradeDaysInYear*3./(TradeDataTradeOpt.Shares./TradeDataTradeOpt.ADV));
  POV = 1./(1+kTR);
  POV = max(POV,TradeDataTradeOpt.Shares./(TradeDataTradeOpt.Shares+totalDays.*TradeDataTradeOpt.ADV));
    
  % TradeTime
  TradeDataTradeOpt.TradeTime = TradeDataTradeOpt.Shares./TradeDataTradeOpt.ADV.*(1-POV)./POV;
    
elseif optType == 2
  
  % Percentage of volume optimization
  numPeriods = 1:totalNumberPeriods;
  k1 = 1:totalNumberPeriods+1;
  k2 = 0:totalNumberPeriods-1;
  TradeDataTradeOpt.TradeTime = (TradeDataTradeOpt.Shares ./ TradeDataTradeOpt.ADV) .* (1-POV) ./ POV;
  T = TradeDataTradeOpt.TradeTime * totalNumberPeriods;
  phi = 2.50 * T.^(-1.25) + 0.01;

  % Trade Schedule
  R = repmat(TradeDataTradeOpt.Shares,1,totalNumberPeriods+1).*exp(-phi*(k1-1));
  X = R(:,1:end-1)-R(:,2:end);  
  
  % Price Appreciations in Dollars
  PA = 1/2 * TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Alpha_bp .* TradeDataTradeOpt.TradeTime;

  % Market Impact in Dollars
  MI = marketImpact(k,TradeDataTradeOpt) .* TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Price ./ 10000;
  
  % Timing Risk in Dollars
  TR = TradeDataTradeOpt.Volatility .* sqrt(1/k.TradeDaysInYear*1/3*(TradeDataTradeOpt.Shares./TradeDataTradeOpt.ADV) .* (1-POV)./POV) .* TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Price;
  
end    

TotMI = sum(MI) / (TradeDataTradeOpt.Shares' * TradeDataTradeOpt.Price) .* 10000;       % bp
TotPA = sum(PA) / (TradeDataTradeOpt.Shares' * TradeDataTradeOpt.Price) .* 10000;   % bp
TotTR = sqrt(trace(R'*CC*R)) ./ (TradeDataTradeOpt.Shares' * TradeDataTradeOpt.Price) * 10000;  

function TotCost = optimizeTradingSchedule(theta,NumberPeriods,NumberStocks,periodsPerDay,TradeDataTradeOpt,CC,Lambda,kObj)
% TotCost = optimizeTradingSchedule(theta,NumberPeriods,NumberStocks,periodsPerDay,TradeDataTradeOpt,CC,Lambda,kObj)
%   TradeMetrics Trade Schedule
%   - Exponential Optimization Function - following:
%   - "The Science of Algorithmic Trading and Portfolio Management," Robert Kissell 2013
%   - "Multi-Asset Risk Modeling," Mort Glantz and Robert Kissell, 2013
%   - "Algorithmic Trading Strategies," Robert Kissell, 2006
%   - "Optimal Trading Strategies," Robert Kissell & Morton Glantz, 2003
%   - "Multi-Period Trade Schedule Optimization," Roberto Malamut & Robert Kissell, 2002

% Notations:
% Theta = Optimization Parameter (1 parameter per stock)
% POV = Optimized POV Rate (1 per stock)
% X = Shares to trade in period i
% R = Residual Shares at beginning of period i


% Matrices Used to Help Speed up the Calculations
kPeriods = 1:NumberPeriods;
K_Matrix = repmat(kPeriods,NumberStocks,1);
Theta_Matrix = repmat(theta,1,NumberPeriods);
TradeDataTradeOpt.VolumeProfile = repmat(TradeDataTradeOpt.ADV/periodsPerDay,1,NumberPeriods);

% Trade Schedule Matrices
% Shares
% X = Shares to Trade in each trade period
Xpct = (exp(-K_Matrix.*Theta_Matrix).*(exp(Theta_Matrix)-1)) ./ (1-exp(-NumberPeriods*Theta_Matrix));
X = repmat(TradeDataTradeOpt.Shares,1,NumberPeriods) .* Xpct;
TradeDataTradeOpt.TradeSchedule = X;

% Residuals
% R = Residual Shares (unexecuted Shares at beginning of each trade period)
Rpct = (exp(-(K_Matrix-1).*Theta_Matrix) - exp(-NumberPeriods.*Theta_Matrix)) ./ (1-exp(-NumberPeriods.*Theta_Matrix));
R = repmat(TradeDataTradeOpt.Shares,1,NumberPeriods) .* Rpct;

% Compute Cost Components

% Price Appreciation
pa = sum(R,2) .* TradeDataTradeOpt.Alpha_bp;
PA = sum(pa);

% Market Impact in dollars
mi= marketImpact(kObj,TradeDataTradeOpt) .* TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Price ./ 10000;
MI = sum(mi);

% Timing Risk
% Compute TR if Lambda (Risk Aversion) is not zero to save on calculations
% The Trace Function in MatLab is much quicker than using a "loop"
TR = 0;
if Lambda ~= 0
    TR = sqrt(trace(R'*CC*R));
end

% Objective Function Calculations
TotCost = PA + MI + Lambda * TR;

end


function TotCost = optimizePOV(POV,NumberPeriods,TradeDataTradeOpt,CC,Lambda,kobj)
% TotCost = optimizePOV(POV,NumberPeriods,TradeDataTradeOpt,CC,Lambda,kobj)
% Derive Trade Schedule using POV parameter
%    T = t*periodsPerDay
%    phi = 2.5*T.^(1.25)+0.01
%
%    phi is set so that the residual schedule at time t expressed in trading intervals j is 10%
%    t is the completion time expressed in days
%    T is the completion time expressed in trading intervals
%    N is the total number of trading intervals

% Matrices Used to Help Speed up the Calculations
TradeDataTradeOpt.POV = POV;
k1 = 1:NumberPeriods + 1;
TradeDataTradeOpt.TradeTime = (TradeDataTradeOpt.Shares ./ TradeDataTradeOpt.ADV) .* (1-TradeDataTradeOpt.POV) ./ TradeDataTradeOpt.POV;
T = TradeDataTradeOpt.TradeTime * NumberPeriods;
phi = 2.50 * T.^(-1.25) + 0.01;

% Residual Trade Schedule
% R = Residual Shares (unexecuted Shares at beginning of each trade period)
R = repmat(TradeDataTradeOpt.Shares,1,NumberPeriods+1) .* exp(-phi*(k1-1));

% Compute Cost Components
% Price Appreciation
pa = 1/2 * TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Alpha_bp .* TradeDataTradeOpt.TradeTime;
PA = sum(pa);

% Market Impact

% Calculate market impact in dollars
mi = marketImpact(kobj,TradeDataTradeOpt) .* TradeDataTradeOpt.Shares .* TradeDataTradeOpt.Price ./ 10000;
MI = sum(mi);

% Timing Risk
% Compute TR if Lambda (Risk Aversion) is not zero to save on calculations
% The Trace Function in MatLab is much quicker than using a "loop"
TR = 0;
if Lambda ~= 0
    TR = sqrt(trace(R'*CC*R));
end

% Objective Function Calculations
TotCost = PA + MI + Lambda * TR;

end