%% KRG Cost Index Example.
% The Cost Index determines the Buy-Sell Imbalance for each stock entered
% by the user. The Real Time Cost Index and the Historical Cost Index are
% computed in the same way.
% 
% Real-Time Cost index uses prices from the open through the current time.
% Historical Cost index uses prices over the full day.
% If the user does not enter the Index Prices, e.g., IndexOpen and
% IndexVWAP, then the IndexCost is defined to be 0.
% 
% Historical Costs will need to use the MI Parameters for the region and
% for the stated Date, and this will change from record to record.

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

% Load the example stock data for current analysis
load KRGExampleData TradeDataCurrent
TradeData = TradeDataCurrent;

% To perform historical analysis, use the commands
% load KRGExampleData TradeDataHistorical
% TradeData = TradeDataHistorical;

% Get number of stocks
numStocks = height(TradeData);

% Set optimization options
options = optimoptions('fsolve','MaxIterations',4000);

% Calculate the total cost, index cost, market impact cost, side of micost
% and initial shares values (x0)
totalCost = TradeData.VWAP ./ TradeData.Open - 1;
indexCost = TradeData.Beta .* (TradeData.IndexVWAP ./ TradeData.IndexOpen - 1);
miCost = totalCost - indexCost;
sideIndicator = sign(miCost);
x0 = 0.5 * TradeData.Volume;

% Output results table
costIndexTable = table;
costIndexTable.Symbol = TradeData.Symbol;
costIndexTable.Date = TradeData.Date;
costIndexTable.Side = sideIndicator;
costIndexTable.TotalVolume = TradeData.Volume;
costIndexTable.TotalCost = totalCost;
costIndexTable.IndexCost = indexCost;
    
% Solving for each stock in a loop due to potential for each stock to have
% different market impact code and date 
for i = 1:numStocks
    
  % Set the MiCode and MiDate of the object for each stock
  k.MiCode = TradeData.MICode(i);
  k.MiDate = TradeData.Date(i);
    
  % Solve for Shares for each stock that will result in the target market
  % impact cost.
  
  % In this example, x is the number of shares (imbalance) that will cause
  % the MI impact cost, that is, the number of shares that will result in a market 
  % impact cost of MI. Here we need to use the abs(MI) since market impact 
  % cost is always positive.  If the market impact cost is 0.0050 then 
  % fsolve is trying to find the number of shares x so that the market 
  % impact formula will return 0.0050.
  % Note that fsolve is using the cost in basis points
  x = fsolve(@(x) costIndexExampleEq(x,miCost(i),TradeData(i,:),k), x0(i), options);
  
  % The imbalance must be between 0 and the actual traded volume
  x = max(min(x,TradeData.Volume(i)),0);
  
  % Recalculate the percentage of volume and shares based on x
  TradeData.POV(i) = x/TradeData.Volume(i);
  TradeData.Shares(i) = x;

  % Calculate the new cost as decimal value
  mi = marketImpact(k,TradeData(i,:))/10000;
  
  % imbalance is the share amount denoted as buy or sell by the sideIndicator
  imbalance = sideIndicator(i) * x;
  
  % Calculate the buy and sell volumes
  % Knowing that
  %
  % Volume = buyVolume + sellVolume;
  % Imbalance = buyVolume - sellVolume;
  %
  % Solve for buyVolume and sellVolume
  buyVolume = (TradeData.Volume(i) + imbalance) / 2;
  sellVolume = (TradeData.Volume(i) - imbalance) / 2;
  
  % Fill output table
  costIndexTable.Imbalance(i,1) = imbalance;
  costIndexTable.ImbalancePctADV(i,1) = imbalance/TradeData.ADV(i);
  costIndexTable.ImbalancePctDayVolume(i,1) = imbalance/TradeData.Volume(i); 
  costIndexTable.BuyVolume(i,1) = buyVolume;
  costIndexTable.SellVolume(i,1) = sellVolume;
  costIndexTable.MI(i,1) = mi * sideIndicator(i);
  costIndexTable.ExcessCost(i,1) = totalCost(i) - mi - indexCost(i);
  
end


% This is the equation we are trying to solve
% We want to find the value of x that makes the MI equation equal to
% abs(MICost). This is also the value of x that makes the value of f = 0

function f = costIndexExampleEq(x,MICost,TradeData,k)
% Cost Index example nonlinear equation to solve

  TradeData.POV = x/TradeData.Volume;
  TradeData.Shares = x;
  mi = marketImpact(k,TradeData);
  f = abs(MICost)*10000 - mi;

end