%% KRG Single Stock Optimization Example
% This example calculates optimized single stock strategies by minimizing
% trading costs subject to a user specified risk aversion parameter.
%
% Min ((MI + PA) + Lambda*TR)
%
% given a trading strategy defined by a user input POV, TradeTime or
% TradeSchedule.
%
% To fully use this example, the user will need the Optimization Toolbox.
%
% The user's input data is a table or structure containing Shares, ADV,
% Volatility and Price.  
%
% The strategy is defined in the data as a field or variable named POV,
% TradeTime or TradeSchedule.   When using TradeSchedule, the data will 
% contain a field or variable named VolumeProfile.  The user's data should
% only contain a single strategy.
%
% Lambda is set from 0 to Inf and denotes the level of risk aversion.
%
% References:
% 1) The Science of Algorithmic Trading and Portfolio Management
% 2) Multi-Asset Risk Modeling
% 3) Optimal Trading Strategy
% 4) Algorithmic Trading Strategies

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

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

% Create input data
tradeData.Shares = 100000;
tradeData.ADV = 1000000;
tradeData.Volatility = 0.25;
tradeData.Price = 35;
tradeData.POV = 0.5;
tradeData.Alpha_bp = 50;     

%% Percentage of Volume optimization

% Risk aversion level
Lambda = 1;

% Define the lower and upper bounds of strategy input for optimization
LB = 0;
UB = 1;

% Define the function handle for the objective function
handle = @(pov)krgSingleStockOptimizer(pov,k,tradeData,Lambda);

% Minimize the percentage of volume trade strategy parameter, fminbnd 
% does not require an initial value of pov but calculates it based on the
% lower and upper bound values.
[tradeData.POV,totalcost] = fminbnd(handle,LB,UB);

% Calculate trading costs based on percentage of volume input for display purposes
mi = marketImpact(k,tradeData);
pa = priceAppreciation(k,tradeData);
tr = timingRisk(k,tradeData);
povCosts = [totalcost mi pa tr];

%% Trade Time optimization

% Remove the field POV from the input structure and use the Trade Time code
% path
tradeData = rmfield(tradeData,'POV');
tradeData.TradeTime = .5;

% Redefine the optimization objective function handle because tradeData has
% changed
handle = @(tradetime)krgSingleStockOptimizer(tradetime,k,tradeData,Lambda);

% Minimize the trade time strategy parameter
[tradeData.TradeTime,totalcost] = fminbnd(handle,LB,UB);

% Calculate trading costs based on trade time input for display purposes
mi = marketImpact(k,tradeData);
tr = timingRisk(k,tradeData);
pa = priceAppreciation(k,tradeData);
tradeTimeCosts = [totalcost mi pa tr];

%% Trade Schedule optimization

% Remove the field TradeTime from the input structure
tradeData = rmfield(tradeData,{'TradeTime'});

% Define the number of trades and the volume per trade
numIntervals = 26;
tradeData.VolumeProfile = ones(1,numIntervals) * tradeData.ADV/numIntervals;
tradeData.TradeSchedule = ones(1,numIntervals) .* (tradeData.Shares./numIntervals);

% Define the lower and upper bounds of shares traded per interval
LB = zeros(1,numIntervals);           
UB = ones(1,numIntervals) .* tradeData.Shares;    

% Shares traded over the schedule must match total shares
Aeq = ones(1,numIntervals);
Beq = tradeData.Shares;

% Define the optimization options 
options = optimoptions('fmincon','MaxFunEvals',100000,'MaxIter',100000);

% Redefine the optimization objective function handle
handle = @(tradeschedule)krgSingleStockOptimizer(tradeschedule,k,tradeData,Lambda);

% Optimize the trade schedule strategy parameter
[tradeData.TradeSchedule,totalcost] = fmincon(handle,tradeData.TradeSchedule,[],[],Aeq,Beq,LB,UB,[],options);

% Calculate trading costs based on trade schedule input for display purposes
mi = marketImpact(k,tradeData);
pa = priceAppreciation(k,tradeData);
tr = timingRisk(k,tradeData);
tradeScheduleCosts = [totalcost mi pa tr];