clc; clear; close all; format bank
%%
tic
%% Set default figure properties
set(0,'DefaultFigureWindowStyle','default') % Set to 'default' to scale figures correctly, else 'docked'
set(0,'DefaultLineLineWidth',1.8)
set(0,'DefaultAxesColorOrder',[0 0 0],...
    'DefaultAxesLineStyleOrder','-|--|:|-.')
set(0,'defaultAxesFontSize',10)
set(0,'defaultfigurecolor',[1 1 1])
%% Load data
addpath '/Users/MikkelFoss/Library/Mobile Documents/com~apple~CloudDocs/Documents/Skole og utdanning/Master/4th Semester/Master Thesis/Data/Analysis/Preliminary_analysis'
addpath '/Users/MikkelFoss/Library/Mobile Documents/com~apple~CloudDocs/Documents/MATLAB'
load('clean_data.mat');
%% Convert data to quarterly frequency 
    freq = 12; % For annualization % DO NOT CHANGE 
flag = 0; % Set 1 for quarterly, 0 for monthly

if flag == 1   
for j = 1:1
    data.HML_devil = [];
    freq = 4; % For annualization 
data2 = data; data2.Date = []; data2 = ret2tick(data2); data2(1,:) = [];
data2.Date = data.Date; data = data2;
data.Date = datestr(data.Date,'dd.mm.yyyy'); data.Date = datenum(data.Date,'dd.mm.yyyy');
startdate = data.Date(1,1); enddate = data.Date(end,1);

busday_quart = busdays(startdate,enddate,'quarterly'); busday_quart = double(ismember(data.Date,busday_quart));
date_dummy = busday_quart.*data.Date; data.Date = date_dummy;
idx = data.Date == 0; data(idx,:) = [];

one = array2table((ones(1,23)));
data = movevars(data,24,'Before',1);

dd = datestr(datenum('31.01.1996','dd.mm.yyyy'),'dd.mm.yyyy'); dates = data.Date; data = data(:,2:end);
one.Properties.VariableNames = data.Properties.VariableNames; data = [one;data];



ret = tick2ret(data); dates = table(dates); data = [dates,ret];

data.dates = datestr(data.dates,'dd.mm.yyyy'); data.Properties.VariableNames(1) = {'Date'};

clear busday_quart data2 date_dummy dates dd enddate idx j one ret startdate
end
else 
end

data.Date = datetime(data.Date,'InputFormat','dd.MM.yyyy');
clear flag 
%% Limit data to Jurek&Stafford sample size (January 1996 - June 2012)
flag = 1; % 1 = full sample, 0 = Jurek and Stafford sample 

if flag == 0
last_busdate = datestr(lbusdate(2012,6)); last_date = find(data.Date == last_busdate);
data = data(1:last_date,:);
clear last_date last_busdate
else 
end
clear flag
%% Descriptive statistics
annual_ret = mean(data{:,2:end},1).*freq; annual_ret(1,18) = mean(data.VIX);   % Create annualized data
annual_vol = std(data{:,2:end}).*sqrt(freq); annual_vol(1,18) = std(data.VIX);
annual_kurt = kurtosis(data{:,2:end},1);
annual_skew = skewness(data{:,2:end},1);

%% Jarque-Bera test for normality (H0 is of normality)
warning('off'); % Hide warning: P is less than the smallest tabulated value, returning 0.001.  
jbstat = zeros(1,21);
p_jb = zeros(1,21);
n = width(data)-1;
x = data{:,2:end};

for i = 1:n
[~,p_jb(1,i),jbstat(1,i)] = jbtest(x(:,i));
end

warning('on');
clear i n
%% Sharpe-ratio
SR = (annual_ret)./annual_vol;

%% Drawdown 
data.Datedum = datenum(data.Date);

% 1998
idx = data.Datedum <= datenum('01.01.1998','dd.mm.yyyy')|...
    data.Datedum >= datenum('01.01.1999','dd.mm.yyyy');
data1998 = data; data1998(idx,:) = []; 

% 2008
idx = data.Datedum <= datenum('01.01.2008','dd.mm.yyyy')|...
    data.Datedum >= datenum('01.01.2009','dd.mm.yyyy');
data2008 = data; data2008(idx,:) = [];

data1998.Datedum = []; data1998.Date = []; data1998 = table2array(data1998);
data2008.Datedum = []; data2008.Date = []; data2008 = table2array(data2008);
data.Datedum = [];

x = data; x.Date = []; x = table2array(x);
[data1998, ~] = ret2tick(data1998);
[data2008, ~] = ret2tick(data2008);
[x, ~] = ret2tick(x);

drawdown = maxdrawdown(x);
drawdown1998 = maxdrawdown(data1998);
drawdown2008 = maxdrawdown(data2008);

clear idx x
%% CAPM for descriptive statistics
X = table2array(data(:,2:end));
[~,n] = size(X);
alpha = zeros(n,1); beta = zeros(n,1);
warning('off', 'stats:pvaluedw:ExactUnavailable');

for i = 1:n
   capm = regstats(X(:,i),data.SP500,'linear');
   alpha(i,1) = capm.beta(1,1);
   beta(i,1) = capm.beta(2,1);
   
end

alpha = alpha*100*12;
clear i n X capm

%% AR(1) model
X = table2array(data(:,2:end));
[~,n] = size(X);
mdl = arima(1,0,0);

for i = 1:n
    AR1(i,:) = estimate(mdl,X(:,i),'Display','off');
end 

for i = 1:n
AR1_Coeff = summarize(AR1(i,:));
AR1_coeff(i,:) = AR1_Coeff.Table(2,1);
AR1_Tstat(i,:) = AR1_Coeff.Table(2,3);
end

AR1_coeff = table2array(AR1_coeff); AR1_Tstat = table2array(AR1_Tstat);
clear X n mdl AR1_Coeff AR1 i 

%% Create table 
Table_descr_stats = table(annual_ret'*100,annual_vol'*100,annual_skew',...
    annual_kurt',jbstat',p_jb',SR',alpha,beta,...
    AR1_coeff,AR1_Tstat,drawdown1998'*100,drawdown2008'*100,drawdown'*100);
Table_descr_stats.Properties.VariableNames = {'Mean','Vol','Skew','Kurt.','JB',...
    'P_JB','SR','Alpha','Beta','AR1_Coeff','AR1_tstat','Drawdown_1998','Drawdown_2008','Drawdown'}; 
Table_descr_stats.Properties.RowNames = data.Properties.VariableNames(2:end);
fprintf('\n ----------------------------------------------------\n')
fprintf('\n ---------- Summary Statistics ---------- \n')
disp(Table_descr_stats)

%% Export table
filename = 'table2_sum_stats.xlsx';
writetable(Table_descr_stats,filename,'Sheet',1,'Range','A1','WriteRowNames',true)
clear filename

%% Regressions for models (HFRI)
% Carhart 4-factor 
x = [data.MKTRF, data.SMB, data.HML, data.MOM];
carhart4f_reg_HFRI = regstats(data.HFRI_pre_fee,x,'linear');
carhart4f_infeasible_ret_HFRI = (carhart4f_reg_HFRI.yhat);
carhart4f_feasible_ret_HFRI = (carhart4f_reg_HFRI.yhat-carhart4f_reg_HFRI.beta(1,1));

% CAPM
capm_reg_HFRI = regstats(data.HFRI_pre_fee,data.MKTRF,'linear');
capm_feasible_ret_HFRI = (capm_reg_HFRI.yhat-capm_reg_HFRI.beta(1,1));
capm_infeasible_ret_HFRI = capm_reg_HFRI.yhat;

% Fung-Hsieh
x = [data.SP500,data.Size_SP500, data.Treasury, data.Credit,...
        data.PTFSBD, data.PTFSFX, data.PTFSCOM, data.PTFSIR, data.PTFSSTK];
fung_hsieh_reg_HFRI = regstats(data.HFRI_pre_fee,x,'linear');
fung_hsieh_feasible_ret_HFRI = (fung_hsieh_reg_HFRI.yhat-fung_hsieh_reg_HFRI.beta(1,1));
fung_hsieh_infeasible_ret_HFRI = (fung_hsieh_reg_HFRI.yhat);

% [Z = -1, L = 2]
S1_reg_HFRI = regstats(data.HFRI_pre_fee, data.S1,'linear');
S1_feasible_ret_HFRI = (S1_reg_HFRI.yhat-S1_reg_HFRI.beta(1,1));
S1_infeasible_ret_HFRI = S1_reg_HFRI.yhat;

% [Z = -2, L = 4]
S2_reg_HFRI = regstats(data.HFRI_pre_fee, data.S2,'linear');
S2_feasible_ret_HFRI = S2_reg_HFRI.yhat-S2_reg_HFRI.beta(1,1);
S2_infeasible_ret_HFRI = S2_reg_HFRI.yhat;

%% Summary tables for regressions
for i = 1:1
% HFRI
fprintf('----------\n')
fprintf('The following tables report coefficents from monthly excess return regressions under \n several risk models for the HFRI \n')
fprintf('\n')

% CAPM
table_capm_reg_HFRI = table([capm_reg_HFRI.beta(1,1)*100,capm_reg_HFRI.beta(2:end,1),capm_reg_HFRI.adjrsquare]',[capm_reg_HFRI.tstat.t',0]'); 
table_capm_reg_HFRI.Properties.RowNames = {'Alpha','RMRF','Adj. R2'}; 
table_capm_reg_HFRI.Properties.VariableNames = {'CAPM','T-stat'};
disp(table_capm_reg_HFRI)
fprintf('\n')

% Carhart 4 Factor model
table_carhart_reg_HFRI = table([carhart4f_reg_HFRI.beta(1,1)*100;carhart4f_reg_HFRI.beta(2:end,1);carhart4f_reg_HFRI.adjrsquare],[carhart4f_reg_HFRI.tstat.t',0]'); 
table_carhart_reg_HFRI.Properties.RowNames = {'Alpha','RMRF','SMB','HML','MOM','Adj. R2'}; 
table_carhart_reg_HFRI.Properties.VariableNames = {'Carhart 4 Factor','T-stat'};
disp(table_carhart_reg_HFRI)
fprintf('\n')

% Fung-Hsieh
table_FH_reg_HFRI = table([fung_hsieh_reg_HFRI.beta(1,1)*100;fung_hsieh_reg_HFRI.beta(2:end,1);fung_hsieh_reg_HFRI.adjrsquare],[fung_hsieh_reg_HFRI.tstat.t',0]'); 
table_FH_reg_HFRI.Properties.RowNames = {'Alpha','RMRF','SIZE','TSY','CREDIT','TF-BD','TF-FX','TF-COM','TF-IR','TF-STK','Adj. R2'}; 
table_FH_reg_HFRI.Properties.VariableNames = {'Fung-Hsieh','T-stat'};
disp(table_FH_reg_HFRI)
fprintf('\n')

% S1
table_s1_reg_HFRI = table([S1_reg_HFRI.beta(1,1)*100,S1_reg_HFRI.beta(2:end,1),S1_reg_HFRI.adjrsquare]',[S1_reg_HFRI.tstat.t',0]');
table_s1_reg_HFRI.Properties.RowNames = {'Alpha','[Z = -1, L = 2]','Adj. R2'};
table_s1_reg_HFRI.Properties.VariableNames = {'S1','T-stat'};
disp(table_s1_reg_HFRI)
fprintf('\n')

% S2
table_s2_reg_HFRI = table([S2_reg_HFRI.beta(1,1)*100,S2_reg_HFRI.beta(2:end,1),S2_reg_HFRI.adjrsquare]',[S2_reg_HFRI.tstat.t',0]');
table_s2_reg_HFRI.Properties.RowNames = {'Alpha','[Z = -2, L = 4]','Adj. R2'};
table_s2_reg_HFRI.Properties.VariableNames = {'S2','T-stat'};
disp(table_s2_reg_HFRI)
fprintf('\n')
%% Export tables to excel
filename = 'table_reg_HFRI.xlsx';
writetable(table_capm_reg_HFRI,filename,'Sheet',1,'Range','A1','WriteRowNames',true)
writetable(table_carhart_reg_HFRI,filename,'Sheet',1,'Range','D1','WriteRowNames',true)
writetable(table_FH_reg_HFRI,filename,'Sheet',1,'Range','A6','WriteRowNames',true)
writetable(table_s1_reg_HFRI,filename,'Sheet',1,'Range','J1','WriteRowNames',true)
writetable(table_s2_reg_HFRI,filename,'Sheet',1,'Range','G1','WriteRowNames',true)
clear filename

%% Regressions for models (DJCS)
% Carhart 4-factor 
x = [data.MKTRF, data.SMB, data.HML, data.MOM];
carhart4f_reg_DJCS = regstats(data.DJCS_pre_fee,x,'linear');
carhart4f_feasible_ret_DJCS = (carhart4f_reg_DJCS.yhat-carhart4f_reg_DJCS.beta(1,1));
carhart4f_infeasible_ret_DJCS = carhart4f_reg_DJCS.yhat;

% CAPM
capm_reg_DJCS = regstats(data.DJCS_pre_fee,data.MKTRF,'linear');
capm_feasible_ret_DJCS = capm_reg_DJCS.yhat-capm_reg_DJCS.beta(1,1);
capm_infeasible_ret_DJCS = capm_reg_DJCS.yhat;

% Fung-Hsieh
x = [data.SP500,data.Size_SP500, data.Treasury, data.Credit,...
        data.PTFSBD, data.PTFSFX, data.PTFSCOM, data.PTFSIR, data.PTFSSTK];
fung_hsieh_reg_DJCS = regstats(data.DJCS_pre_fee,x,'linear');
fung_hsieh_feasible_ret_DJCS = (fung_hsieh_reg_DJCS.yhat-fung_hsieh_reg_DJCS.beta(1,1));
fung_hsieh_infeasible_ret_DJCS = fung_hsieh_reg_DJCS.yhat;

% [Z = -1, L = 2]
S1_reg_DJCS = regstats(data.DJCS_pre_fee, data.S1,'linear');
S1_feasible_ret_DJCS = S1_reg_DJCS.yhat-S1_reg_DJCS.beta(1,1);
S1_infeasible_ret_DJCS = S1_reg_DJCS.yhat;

% [Z = -2, L = 4]
S2_reg_DJCS = regstats(data.DJCS_pre_fee, data.S2,'linear');
S2_feasible_ret_DJCS = S2_reg_DJCS.yhat-S2_reg_DJCS.beta(1,1);
S2_infeasible_ret_DJCS = S2_reg_DJCS.yhat;

%% Summary tables for regressions on DJCS

fprintf('---------- \n')
fprintf('The following tables report coefficents from monthly excess return regressions under \n several risk models for the DJCS \n')
fprintf('\n')

% CAPM
table_capm_reg_DJCS = table([capm_reg_DJCS.beta(1,1)*100,capm_reg_DJCS.beta(2:end,1),capm_reg_DJCS.adjrsquare]',[capm_reg_DJCS.tstat.t',0]'); 
table_capm_reg_DJCS.Properties.RowNames = {'Alpha','RMRF','Adj. R2'}; 
table_capm_reg_DJCS.Properties.VariableNames = {'CAPM','T-stat'};
disp(table_capm_reg_DJCS)
fprintf('\n')

% Carhart 4 Factor model
table_carhart_reg_DJCS = table([carhart4f_reg_DJCS.beta(1,1)*100;carhart4f_reg_DJCS.beta(2:end,1);carhart4f_reg_DJCS.adjrsquare],[carhart4f_reg_DJCS.tstat.t',0]'); 
table_carhart_reg_DJCS.Properties.RowNames = {'Alpha','RMRF','SMB','HML','MOM','Adj. R2'}; 
table_carhart_reg_DJCS.Properties.VariableNames = {'Carhart 4 Factor','T-stat'};
disp(table_carhart_reg_DJCS)
fprintf('\n')

% Fung-Hsieh
table_FH_reg_DJCS = table([fung_hsieh_reg_DJCS.beta(1,1)*100;fung_hsieh_reg_DJCS.beta(2:end,1);fung_hsieh_reg_DJCS.adjrsquare],[fung_hsieh_reg_DJCS.tstat.t',0]'); 
table_FH_reg_DJCS.Properties.RowNames = {'Alpha','RMRF','SIZE','TSY','CREDIT','TF-BD','TF-FX','TF-COM','TF-IR','TF-STK','Adj. R2'}; 
table_FH_reg_DJCS.Properties.VariableNames = {'Fung-Hsieh','T-stat'};
disp(table_FH_reg_DJCS)
fprintf('\n')

% S1
table_s1_reg_DJCS = table([S1_reg_DJCS.beta(1,1)*100,S1_reg_DJCS.beta(2:end,1),S1_reg_DJCS.adjrsquare]',[S1_reg_DJCS.tstat.t',0]');
table_s1_reg_DJCS.Properties.RowNames = {'Alpha','[Z = -1, L = 2]','Adj. R2'};
table_s1_reg_DJCS.Properties.VariableNames = {'S1','T-stat'};
disp(table_s1_reg_DJCS)
fprintf('\n')

% S2
table_s2_reg_DJCS = table([S2_reg_DJCS.beta(1,1)*100,S2_reg_DJCS.beta(2:end,1),S2_reg_DJCS.adjrsquare]',[S2_reg_DJCS.tstat.t',0]');
table_s2_reg_DJCS.Properties.RowNames = {'Alpha','[Z = -2, L = 4]','Adj. R2'};
table_s2_reg_DJCS.Properties.VariableNames = {'S2','T-stat'};
disp(table_s2_reg_DJCS)
fprintf('\n')

%% Export tables to excel
filename = 'table_reg_DJCS.xlsx';
writetable(table_capm_reg_DJCS,filename,'Sheet',1,'Range','A1','WriteRowNames',true)
writetable(table_carhart_reg_DJCS,filename,'Sheet',1,'Range','D1','WriteRowNames',true)
writetable(table_FH_reg_DJCS,filename,'Sheet',1,'Range','A6','WriteRowNames',true)
writetable(table_s1_reg_DJCS,filename,'Sheet',1,'Range','J1','WriteRowNames',true)
writetable(table_s2_reg_DJCS,filename,'Sheet',1,'Range','G1','WriteRowNames',true)

clear filename
%% Create cumulative returns using HFRI  %% Returns are NOT excess %%

carhart4f_feasible_cumret = ret2tick(carhart4f_feasible_ret_HFRI+data.rf);
carhart4f_infeasible_cumret = ret2tick(carhart4f_infeasible_ret_HFRI+data.rf);
capm_feasible_cumret = ret2tick(capm_feasible_ret_HFRI+data.rf);
capm_infeasible_cumret = ret2tick(capm_infeasible_ret_HFRI+data.rf);
fung_hsieh_feasible_cumret = ret2tick(fung_hsieh_feasible_ret_HFRI+data.rf);
fung_hsieh_infeasible_cumret = ret2tick(fung_hsieh_infeasible_ret_HFRI+data.rf);

S1_feasible_cumret = ret2tick(data.S1+data.rf); 
S2_feasible_cumret = ret2tick(data.S2+data.rf);
HFRI_cumret = ret2tick(data.HFRI_pre_fee+data.rf);

end
%% Create plots
for i = 1:1
dates = [datetime('30.01.1996','InputFormat','dd.MM.yyyy');data.Date];
width_fig = 4;
height_fig = 4;
titlesize = 13;

f = figure('Units','inches',...
        'Position',[2 5 width_fig height_fig],...
        'PaperPositionMode','auto',...
        'Name','Feasible Linear Replicating Portfolio');

plot_feasible_linear = plot(dates,[HFRI_cumret,capm_feasible_cumret,...
    carhart4f_feasible_cumret,fung_hsieh_feasible_cumret]);
hold on
set(gca,...
'Units','normalized',... 
'YTick',1:1:20,...
'Position',[.15 .2 .75 .7],...
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

xlim([dates(1) dates(end)])
ylim([1 16])
legend({'HFRI','CAPM','Carhart','Fung-Hsieh'},... 
    'interpreter','latex',... 
    'FontSize',11,... 
    'FontName','Times',... 
    'Location','NorthWest');

title('Feasible Linear Replicating Portfolio',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times');

print -deps fig_feasible_linear_rp.eps
%%
figure('Units','inches',...
        'Position',[6 5 width_fig height_fig],...
        'PaperPositionMode','auto',...
        'Name','Infeasible Linear Replicating Portfolio');

plot_infeasible_linear = plot(dates,[HFRI_cumret,capm_infeasible_cumret,...
    carhart4f_infeasible_cumret,fung_hsieh_infeasible_cumret]);
hold on
set(gca,...
'Units','normalized',... 
'YTick',1:1:20,...
'Position',[.15 .2 .75 .7],... 
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

xlim([dates(1) dates(end)])
ylim([1 16])

legend ({'HFRI','CAPM','Carhart','Fung Hsieh'},... 
    'interpreter','latex',... 
    'FontSize',11,... 
    'FontName','Times',... 
    'Location','NorthWest');

title('Infeasible Linear Replicating Portfolio',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times');

print -deps fig_infeasible_linear_rp.eps
%% Plot 'Feasible Nonolinear Replicating Portfolio'
figure('Units','inches',...
        'Position',[10 5 width_fig height_fig],...
        'PaperPositionMode','auto',...
        'Name','Feasible Nonolinear Replicating Portfolio');

plot_feasible_nonlinear = plot(dates,[HFRI_cumret, S1_feasible_cumret, S2_feasible_cumret]);
hold on 
set(gca,...
'Units','normalized',... 
'YTick',1:1:20,...
'Position',[.15 .2 .75 .7],... 
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

ylim([1 16])
legend ({'HFRI','[Z = -1, L = 2]','[Z = -2, L = 4]'},...
     'interpreter','latex',... 
    'FontSize',11,... 
    'FontName','Times',... 
    'Location','NorthWest');
title('Feasible Nonolinear Replicating Portfolio',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times');

print -deps fig_feasible_nonlinear_rp.eps
%% Drawdown plot 
    % Created using HFRI

% Set up variables 
n = height(data)+1;
HFRI_drawdown_plot = zeros(n,1);
fh_drawdown_plot_feasible = zeros(n,1);
capm_drawdown_plot_feasible = zeros(n,1);
carhart_drawdown_plot_feasible = zeros(n,1);
s1_drawdown_plot_feasible = zeros(n,1);
s2_drawdown_plot_feasible = zeros(n,1);

fh_drawdown_plot_infeasible = zeros(n,1);
capm_drawdown_plot_infeasible = zeros(n,1);
carhart_drawdown_plot_infeasible = zeros(n,1);

% Run loop to generate drawdown path
for i = 2:n
   HFRI_drawdown_plot(i,1) = HFRI_cumret(i,1)/max(HFRI_cumret(1:i,1));
   fh_drawdown_plot_feasible(i,1) = fung_hsieh_feasible_cumret(i,1)/max(fung_hsieh_feasible_cumret(1:i,1));
   capm_drawdown_plot_feasible(i,1) = capm_feasible_cumret(i,1)/max(capm_feasible_cumret(1:i,1));
   carhart_drawdown_plot_feasible(i,1) = carhart4f_feasible_cumret(i,1)/max(carhart4f_feasible_cumret(1:i,1));
   s1_drawdown_plot_feasible(i,1) = S1_feasible_cumret(i,1)/max(S1_feasible_cumret(1:i,1));
   s2_drawdown_plot_feasible(i,1) = S2_feasible_cumret(i,1)/max(S2_feasible_cumret(1:i,1));
   
   fh_drawdown_plot_infeasible(i,1) = fung_hsieh_infeasible_cumret(i,1)/max(fung_hsieh_infeasible_cumret(1:i,1));
   capm_drawdown_plot_infeasible(i,1) = capm_infeasible_cumret(i,1)/max(capm_infeasible_cumret(1:i,1));
   carhart_drawdown_plot_infeasible(i,1) = carhart4f_infeasible_cumret(i,1)/max(carhart4f_infeasible_cumret(1:i,1));
   
end

%% Plot Drawdown - Linear portfolios - Feasible
figure('Units','inches',...
        'Position',[2 0 width_fig height_fig],...
        'PaperPositionMode','auto',...
        'Name','Drawdown feasible linear replicating portfolios');

x = [HFRI_drawdown_plot,fh_drawdown_plot_feasible,capm_drawdown_plot_feasible, carhart_drawdown_plot_feasible]; x(1,:) = [];
plot(data.Date, x-1);
hold on
set(gca,...
'Units','normalized',... 
'YTick',-0.3:0.05:0,...
'Position',[.15 .2 .75 .7],... 
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

ylim([-0.3 0])
legend ({'HFRI','Fung-Hsieh','CAPM','Carhart'},...
    'interpreter','latex',... 
    'FontSize',11,... 
    'FontName','Times',... 
    'Location','SouthWest');
title('Drawdown',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times');
clear x

print -deps fig_drawdown_feasible_linear_rp.eps
%% Plot Drawdown - Linear portfolios - Infeasible
figure('Units','inches',...
        'Position',[6 0 width_fig height_fig],...
        'PaperPositionMode','auto',...
        'Name','Drawdown infeasible linear replicating portfolios');
    
x = [HFRI_drawdown_plot,fh_drawdown_plot_infeasible,capm_drawdown_plot_infeasible, carhart_drawdown_plot_infeasible]; x(1,:) = [];
plot(data.Date, x-1);
hold on
set(gca,...
'Units','normalized',... 
'YTick',-0.3:0.05:0,...
'Position',[.15 .2 .75 .7],... 
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

ylim([-0.3 0])
legend ({'HFRI','Fung-Hsieh','CAPM','Carhart 4F'},...
    'interpreter','latex',... 
    'FontSize',11,... 
    'FontName','Times',... 
    'Location','SouthWest');

title('Drawdown',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times');
clear x

print -deps fig_drawdown_infeasible_linear_rp.eps
%% Plot Drawdown - Nonlinear portfolios
figure('Units','inches',...
        'Position',[10 0 width_fig height_fig],...
        'PaperPositionMode','auto',...
        'Name','Drawdown feasible nonlinear replicating portfolios');
    
x = [HFRI_drawdown_plot,s1_drawdown_plot_feasible,s2_drawdown_plot_feasible]; x(1,:) = [];
plot(data.Date, x-1);
hold on
set(gca,...
'Units','normalized',... 
'YTick',-0.3:0.05:0,...
'Position',[.15 .2 .75 .7],... 
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

ylim([-0.3 0])
legend ({'HFRI','[Z = -1, L = 2]','[Z = -2, L = 4]'},...
     'interpreter','latex',... 
    'FontSize',11,... 
    'FontName','Times',... 
    'Location','SouthWest');
title('Drawdown',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times');
clear x

print -deps fig_drawdown_feasible_nonlinear_rp.eps
end 


%% Joint hypothesis test (Alpha = 0, Beta = 1)
for i = 1:1
ss_tot_HFRI = sum((data.HFRI_pre_fee-(mean(data.HFRI_pre_fee))).^2);
ss_tot_DJCS = sum((data.DJCS_pre_fee-(mean(data.DJCS_pre_fee))).^2);
    
    % Run regression of feasible returns on HFRI
jh_capm_HFRI = regstats(data.HFRI_pre_fee,capm_feasible_ret_HFRI,'linear');
jh_carhart4f_HFRI = regstats(data.HFRI_pre_fee,carhart4f_feasible_ret_HFRI,'linear');
jh_fh_HFRI = regstats(data.HFRI_pre_fee,fung_hsieh_feasible_ret_HFRI,'linear');
jh_s1_HFRI = regstats(data.HFRI_pre_fee,S1_feasible_ret_HFRI,'linear');
jh_s2_HFRI = regstats(data.HFRI_pre_fee,S2_feasible_ret_HFRI,'linear');

% Compute residuals exclusive of intercept 
ss_res_capm_HFRI = sum((data.HFRI_pre_fee-(jh_capm_HFRI.yhat-jh_capm_HFRI.beta(1,1))).^2);
ss_res_carhart4f_HFRI = sum((data.HFRI_pre_fee-(jh_carhart4f_HFRI.yhat-jh_carhart4f_HFRI.beta(1,1))).^2);
ss_res_fh_HFRI = sum((data.HFRI_pre_fee-(jh_fh_HFRI.yhat-jh_fh_HFRI.beta(1,1))).^2);
ss_res_s1_HFRI = sum((data.HFRI_pre_fee-(jh_s1_HFRI.yhat-jh_s1_HFRI.beta(1,1))).^2);
ss_res_s2_HFRI = sum((data.HFRI_pre_fee-(jh_s2_HFRI.yhat-jh_s2_HFRI.beta(1,1))).^2);

% Compute adj.R2
adjr2_feasible_capm_HFRI = (1-(1-((1-(ss_res_capm_HFRI/ss_tot_HFRI))))*((n-1)/(n-2)));
adjr2_feasible_carhart4f_HFRI = (1-(1-((1-(ss_res_carhart4f_HFRI/ss_tot_HFRI))))*((n-1)/(n-2)));
adjr2_feasible_fh_HFRI = (1-(1-((1-(ss_res_fh_HFRI/ss_tot_HFRI))))*((n-1)/(n-2)));
adjr2_feasible_s1_HFRI = (1-(1-((1-(ss_res_s1_HFRI/ss_tot_HFRI))))*((n-1)/(n-2)));
adjr2_feasible_s2_HFRI = (1-(1-((1-(ss_res_s2_HFRI/ss_tot_HFRI))))*((n-1)/(n-2)));

% Run test
h = [1 0;0 1]; c = [0;1];
[p_capm_HFRI,~] = linhyptest(jh_capm_HFRI.beta, jh_capm_HFRI.covb, c, h, jh_capm_HFRI.fstat.dfe);
[p_carhart4f_HFRI,~] = linhyptest(jh_carhart4f_HFRI.beta, jh_carhart4f_HFRI.covb, c, h, jh_carhart4f_HFRI.fstat.dfe);
[p_fh_HFRI,~] = linhyptest(jh_fh_HFRI.beta, jh_fh_HFRI.covb, c, h, jh_fh_HFRI.fstat.dfe);
[p_s1_HFRI,~] = linhyptest(jh_s1_HFRI.beta, jh_s1_HFRI.covb, c, h, jh_s1_HFRI.fstat.dfe);
[p_s2_HFRI,~] = linhyptest(jh_s2_HFRI.beta, jh_s2_HFRI.covb, c, h, jh_s2_HFRI.fstat.dfe);

% Table
jh_HFRI = table([p_capm_HFRI,p_carhart4f_HFRI,p_fh_HFRI,p_s1_HFRI,p_s2_HFRI]',100.*[adjr2_feasible_capm_HFRI,adjr2_feasible_carhart4f_HFRI,adjr2_feasible_fh_HFRI,adjr2_feasible_s1_HFRI,adjr2_feasible_s2_HFRI]');
jh_HFRI.Properties.RowNames = {'CAPM','Carhart','FH','S1','S2'};
jh_HFRI.Properties.VariableNames ={'P Value','Adj.R2 [feasible]'};
fprintf('\n ----------------------------------------------------\n')
fprintf('P-value for joint hypothesis test and feasible adj.R2 \n - HFRI \n')
disp(jh_HFRI)

    % Run regression of feasible returns on DJCS
jh_capm_DJCS = regstats(data.DJCS_pre_fee,capm_feasible_ret_DJCS,'linear');
jh_carhart4f_DJCS = regstats(data.DJCS_pre_fee,carhart4f_feasible_ret_DJCS,'linear');
jh_fh_DJCS = regstats(data.DJCS_pre_fee,fung_hsieh_feasible_ret_DJCS,'linear');
jh_s1_DJCS = regstats(data.DJCS_pre_fee,S1_feasible_ret_DJCS,'linear');
jh_s2_DJCS = regstats(data.DJCS_pre_fee,S2_feasible_ret_DJCS,'linear');

% Compute sum squared residuals exclusive of intercept 
ss_res_capm_DJCS = sum((data.DJCS_pre_fee-(jh_capm_DJCS.yhat-jh_capm_DJCS.beta(1,1))).^2);
ss_res_carhart4f_DJCS = sum((data.DJCS_pre_fee-(jh_carhart4f_DJCS.yhat-jh_carhart4f_DJCS.beta(1,1))).^2);
ss_res_fh_DJCS = sum((data.DJCS_pre_fee-(jh_fh_DJCS.yhat-jh_fh_DJCS.beta(1,1))).^2);
ss_res_s1_DJCS = sum((data.DJCS_pre_fee-(jh_s1_DJCS.yhat-jh_s1_DJCS.beta(1,1))).^2);
ss_res_s2_DJCS = sum((data.DJCS_pre_fee-(jh_s2_DJCS.yhat-jh_s2_DJCS.beta(1,1))).^2);

% Compute adj.R2
adjr2_feasible_capm_DJCS = (1-(1-((1-(ss_res_capm_DJCS/ss_tot_DJCS))))*((n-1)/(n-2)));
adjr2_feasible_carhart4f_DJCS = (1-(1-((1-(ss_res_carhart4f_DJCS/ss_tot_DJCS))))*((n-1)/(n-2)));
adjr2_feasible_fh_DJCS = (1-(1-((1-(ss_res_fh_DJCS/ss_tot_DJCS))))*((n-1)/(n-2)));
adjr2_feasible_s1_DJCS = (1-(1-((1-(ss_res_s1_DJCS/ss_tot_DJCS))))*((n-1)/(n-2)));
adjr2_feasible_s2_DJCS = (1-(1-((1-(ss_res_s2_DJCS/ss_tot_DJCS))))*((n-1)/(n-2)));

% Run test 
h = [1 0;0 1]; c = [0;1];
[p_capm_DJCS,~] = linhyptest(jh_capm_DJCS.beta, jh_capm_DJCS.covb, c, h, jh_capm_DJCS.fstat.dfe);
[p_carhart4f_DJCS,~] = linhyptest(jh_carhart4f_DJCS.beta, jh_carhart4f_DJCS.covb, c, h, jh_carhart4f_DJCS.fstat.dfe);
[p_fh_DJCS,~] = linhyptest(jh_fh_DJCS.beta, jh_fh_DJCS.covb, c, h, jh_fh_DJCS.fstat.dfe);
[p_s1_DJCS,~] = linhyptest(jh_s1_DJCS.beta, jh_s1_DJCS.covb, c, h, jh_s1_DJCS.fstat.dfe);
[p_s2_DJCS,~] = linhyptest(jh_s2_DJCS.beta, jh_s2_DJCS.covb, c, h, jh_s2_DJCS.fstat.dfe);

% Table
jh_DJCS = table([p_capm_DJCS,p_carhart4f_DJCS,p_fh_DJCS,p_s1_DJCS,p_s2_DJCS]',100.*[adjr2_feasible_capm_DJCS,adjr2_feasible_carhart4f_DJCS,adjr2_feasible_fh_DJCS,adjr2_feasible_s1_DJCS,adjr2_feasible_s2_DJCS]');
jh_DJCS.Properties.RowNames = {'CAPM','Carhart','FH','S1','S2'};
jh_DJCS.Properties.VariableNames ={'P Value','Adj.R2 [feasible]'};
fprintf('\n ----------------------------------------------------\n')
fprintf('P-value for joint hypothesis test and feasible adj.R2 \n - DJCS \n')
disp(jh_DJCS)
end

%% Correlation between factors and Put-writing 
x = [data.HML, data.SMB, data.MOM, data.MKTRF, data.QMJ,...
    data.BAB, data.PTFSBD, data.PTFSFX, data.PTFSCOM,...
    data.PTFSIR, data.PTFSSTK, data.Treasury, data.Credit, data.VIX];

% Correlation matrix 
[corr_mat,corr_mat_p] = corrcoef([data.S1,data.S2,x]);
corr_mat = array2table(corr_mat);
corr_mat.Properties.RowNames = {'S1','S2','HML','SMB','MOM','MKTRF','QMJ','BAB','PTFSBD','PTFSFX','PTFSCOM','PTFSIR','PTFSSTK','Treasury','Credit','VIX'};
corr_mat.Properties.VariableNames = {'S1','S2','HML','SMB','MOM','MKTRF','QMJ','BAB','PTFSBD','PTFSFX','PTFSCOM','PTFSIR','PTFSSTK','Treasury','Credit','VIX'};

corr_mat_p = array2table(corr_mat_p);
corr_mat_p.Properties.RowNames = {'S1','S2','HML','SMB','MOM','MKTRF','QMJ','BAB','PTFSBD','PTFSFX','PTFSCOM','PTFSIR','PTFSSTK','Treasury','Credit','VIX'};
corr_mat_p.Properties.VariableNames = {'S1','S2','HML','SMB','MOM','MKTRF','QMJ','BAB','PTFSBD','PTFSFX','PTFSCOM','PTFSIR','PTFSSTK','Treasury','Credit','VIX'};

fprintf('\n ----------------------------------------------------\n')
fprintf('Correlation Matrix \n')
fprintf('\n')
disp(corr_mat);
fprintf('\n')
fprintf('\n ----------------------------------------------------\n')
fprintf('P-Values from Correlation Matrix \n')
fprintf('\n')
disp(corr_mat_p);

clear x
%% Rolling window Correlation
window = 30;
x = [data.HFRI_pre_fee, data.DJCS_pre_fee, data.S1, data.S2,...
    data.HML, data.SMB, data.MOM, data.MKTRF, data.QMJ,...
    data.BAB];

r_s1 = zeros(size(x));
p_s1 = zeros(size(x));

r_s2 = zeros(size(x));
p_s2 = zeros(size(x));

r_credit = zeros(size(x));
p_credit = zeros(size(x));

for i = 1:size(x,2)
[r_s1(:,i), p_s1(:,i)] = movcorr(data.S1, x(:,i),window);
[r_s2(:,i), p_s2(:,i)] = movcorr(data.S2, x(:,i),window);
[r_credit(:,i), p_credit(:,i)] = movcorr(data.Credit, x(:,i),window);
end

table_rolling_corr_s1 = array2table(r_s1);
table_rolling_corr_s2 = array2table(r_s2);
table_rolling_corr_credit = array2table(r_credit);

table_rolling_corr_s1.Properties.VariableNames = {'HFRI','DJCS','S1','S2','HML','SMB','MOM','MKTRF','QMJ','BAB'};
table_rolling_corr_s2.Properties.VariableNames = {'HFRI','DJCS','S1','S2','HML','SMB','MOM','MKTRF','QMJ','BAB'};
table_rolling_corr_credit.Properties.VariableNames = {'HFRI','DJCS','S1','S2','HML','SMB','MOM','MKTRF','QMJ','BAB'};
%% Plot Rolling correlations
figure('Name','Correlation between Variables')
hold on
tiledlayout(3,1,'Padding','compact');
nexttile 
plot(data.Date,[table_rolling_corr_credit.S1]);
title('Correlation between [Z=-1, L=2] and Credit')
yline(0,'r','Linewidth',2);
hold on 
ylim([-1 1])
legend({'[Z=-1,L=2]'})

nexttile 
plot(data.Date,[table_rolling_corr_credit.QMJ,table_rolling_corr_credit.BAB]);
yline(0,'r','Linewidth',2);
title('Correlation between Credit and BAB/QMJ factors ')
hold on 
ylim([-1 1])
legend({'QMJ','BAB'})

nexttile 
plot(data.Date,[table_rolling_corr_s1.HFRI,table_rolling_corr_s1.DJCS]);
yline(0,'r','Linewidth',2);
title('Correlation between [Z=-1, L=2] and DJCS and HFRI')
hold on 
ylim([-1 1])
legend({'HFRI','DJCS'})

%% Regress factors on S1/S2
y = [data.HML, data.SMB, data.MOM, data.MKTRF, data.QMJ, data.BAB];

for i = 1:size(y,2)
S1_reg_factors(i) = regstats(y(:,i),data.S1,'linear');
S2_reg_factors(i) = regstats(y(:,i),data.S2,'linear');

tstat_s1(:,i) = S1_reg_factors(1,i).tstat.t;
tstat_s2(:,i) = S2_reg_factors(1,i).tstat.t;
end

%% Create table of regressions 
A1 = [S1_reg_factors.beta];
B1 = A1(2,:);
A1 = A1(1,:).*100;

table_S1_reg_factors = table([A1;B1]',tstat_s1',100*[S1_reg_factors.adjrsquare]');
table_S1_reg_factors.Properties.RowNames = {'HML','SMB','MOM','MKTRF','QMJ','BAB'};
table_S1_reg_factors.Properties.VariableNames = {'Alpha/beta','T-stats','Adj.R2'};
fprintf('\n ----------------------------------------------------\n')
fprintf('Regressions of Factors on S1 (Factor = beta*Put-writing) \n')
disp(table_S1_reg_factors)

A = [S2_reg_factors.beta];
B = A(2,:);
A = A(1,:).*100;

table_S2_reg_factors = table([A;B]',tstat_s2',100*[S2_reg_factors.adjrsquare]');
table_S2_reg_factors.Properties.RowNames = {'HML','SMB','MOM','MKTRF','QMJ','BAB'};
table_S2_reg_factors.Properties.VariableNames = {'Alpha/beta','T-stats','Adj.R2'};
fprintf('\n ----------------------------------------------------\n')
fprintf('Regressions of Factors on S2 (Factor = beta*Put-writing) \n')
disp(table_S2_reg_factors)

%% BAB/QMJ = FF + Put-writing 
% S1
x = [data.MKTRF,data.SMB,data.HML,data.MOM,data.S1];
bab_ff_reg_s1 = regstats(data.BAB,x,'linear');
qmj_ff_reg_s1 = regstats(data.QMJ,x,'linear');

fprintf('\n ----------------------------------------------------\n')

fprintf('Regressing BAB on FF factors and S1 \n')
tab_bab_ff_reg_s1 = table([100*bab_ff_reg_s1.beta(1,1);bab_ff_reg_s1.beta(2:end,:);100*bab_ff_reg_s1.adjrsquare],[bab_ff_reg_s1.tstat.t;0]);
tab_bab_ff_reg_s1.Properties.RowNames = {'Intercept x100','MKTRF','SMB','HML','MOM','S1','Adj.R2'};
tab_bab_ff_reg_s1.Properties.VariableNames = {'Coeff','T-stat'};
disp(tab_bab_ff_reg_s1)

% QMJ
fprintf('Regressing QMJ on FF factors and S1 \n')
tab_qmj_ff_reg_s1 = table([100*qmj_ff_reg_s1.beta(1,1);qmj_ff_reg_s1.beta(2:end,:);100*qmj_ff_reg_s1.adjrsquare],[qmj_ff_reg_s1.tstat.t;0]);
tab_qmj_ff_reg_s1.Properties.RowNames = {'Intercept x100','MKTRF','SMB','HML','MOM','S1','Adj.R2'};
tab_qmj_ff_reg_s1.Properties.VariableNames = {'Coeff','T-stat'};
disp(tab_qmj_ff_reg_s1)

fprintf('\n ----------------------------------------------------\n')

% S2
x = [data.MKTRF,data.SMB,data.HML,data.MOM,data.S2];
bab_ff_reg_s2 = regstats(data.BAB,x,'linear');
qmj_ff_reg_s2 = regstats(data.QMJ,x,'linear');

fprintf('Regressing BAB on FF factors and S2 \n')
tab_bab_ff_reg_s2 = table([100*bab_ff_reg_s2.beta(1,1);bab_ff_reg_s2.beta(2:end,:);100*bab_ff_reg_s2.adjrsquare],[bab_ff_reg_s2.tstat.t;0]);
tab_bab_ff_reg_s2.Properties.RowNames = {'Intercept x100','MKTRF','SMB','HML','MOM','S2','Adj.R2'};
tab_bab_ff_reg_s2.Properties.VariableNames = {'Coeff','T-stat'};
disp(tab_bab_ff_reg_s2)

% QMJ
fprintf('Regressing QMJ on FF factors and S2 \n')
tab_qmj_ff_reg_s2 = table([100*qmj_ff_reg_s2.beta(1,1);qmj_ff_reg_s2.beta(2:end,:);100*qmj_ff_reg_s2.adjrsquare],[qmj_ff_reg_s2.tstat.t;0]);
tab_qmj_ff_reg_s2.Properties.RowNames = {'Intercept x100','MKTRF','SMB','HML','MOM','S2','Adj.R2'};
tab_qmj_ff_reg_s2.Properties.VariableNames = {'Coeff','T-stat'};
disp(tab_qmj_ff_reg_s2)

%% BAB/QMJ = FF+Carhart
x = [data.MKTRF,data.SMB,data.HML,data.MOM];
bab_ff_reg = regstats(data.BAB,x,'linear');
qmj_ff_reg = regstats(data.QMJ,x,'linear');

fprintf('\n ----------------------------------------------------\n')

% BAB
fprintf('Regressing BAB on FF factors \n')
tab_bab_ff_reg = table([100*bab_ff_reg.beta(1,1);bab_ff_reg.beta(2:end,:);100*bab_ff_reg.adjrsquare],[bab_ff_reg.tstat.t;0]);
tab_bab_ff_reg.Properties.RowNames = {'Intercept x100','MKTRF','SMB','HML','MOM','Adj.R2'};
tab_bab_ff_reg.Properties.VariableNames = {'Coeff','T-stat'};
disp(tab_bab_ff_reg)

% QMJ
fprintf('Regressing QMJ on FF factors\n')
tab_qmj_ff_reg = table([100*qmj_ff_reg.beta(1,1);qmj_ff_reg.beta(2:end,:);100*qmj_ff_reg.adjrsquare],[qmj_ff_reg.tstat.t;0]);
tab_qmj_ff_reg.Properties.RowNames = {'Intercept x100','MKTRF','SMB','HML','MOM','Adj.R2'};
tab_qmj_ff_reg.Properties.VariableNames = {'Coeff','T-stat'};
disp(tab_qmj_ff_reg)

%% Non-normality plots
width_fig = 4;
height_fig = 4;
titlesize = 13;

%% S1 QQ-Plot
figure('Units','inches',...
    'Position',[2 5 width_fig height_fig],...
    'PaperPositionMode','auto',...
    'Name','QQ-plot HFRI');

a = qqplot(data.HFRI_pre_fee);
set(a, 'Color', 'k','MarkerEdgeColor', 'k')
hold on
set(gca,...
'Units','normalized',... 
'Position',[.10 .10 .80 .80],...
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

ylim([-0.17 0.07])
title('Quartile-Quartile plot',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times',...
    'interpreter','latex');

xlabel('')

print -deps fig_qq_HFRI_pre_fee.eps
matlab2tikz('fig_qq_HFRI_pre_fee.tex','width','\figW','height','\figH','showInfo', false);
%% S2 - QQ-Plot
figure('Units','inches',...
    'Position',[2 5 width_fig height_fig],...
    'PaperPositionMode','auto',...
    'Name','QQ-plot [Z=-2,L=4]');

a = qqplot(data.S2);
set(a, 'Color', 'k','MarkerEdgeColor', 'k')
hold on
set(gca,...
'Units','normalized',... 
'Position',[.10 .10 .80 .80],...
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')

ylim([-0.17 0.07])
title('Quartile-Quartile plot',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times',...
    'interpreter','latex');

xlabel('')

print -deps fig_qq_s2.eps
matlab2tikz('fig_qq_s2.tex','width','\figW','height','\figH','showInfo', false);
%% HFRI
figure('Units','inches',...
    'Position',[2 5 width_fig height_fig],...
    'PaperPositionMode','auto',...
    'Name','Histogram HFRI');
a = histfit(data.HFRI_pre_fee);
set(a(1),'facecolor',([0.5 0.5 0.5]),'edgecolor','k'); 
set(a(2),'color','k')
hold on
set(gca,...
'Units','normalized',... 
'Position',[.15 .15 .80 .75],...
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')
ylim([0 160])
xlim([-0.20 0.10])

title('Distribution of returns',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times',...
    'interpreter','latex');

print -deps fig_histfit_HFRI_pre_fee.eps
matlab2tikz('fig_histfit_HFRI_pre_fee.tex','width','\figW','height','\figH','showInfo', false);
%% S2
figure('Units','inches',...
    'Position',[2 5 width_fig height_fig],...
    'PaperPositionMode','auto',...
    'Name','Histogram [Z=-2,L=4]');
a = histfit(data.S2);
set(a(1),'facecolor',([0.5 0.5 0.5]),'edgecolor','k'); 
set(a(2),'color','k')
hold on
set(gca,...
'Units','normalized',... 
'Position',[.15 .15 .80 .75],...
'FontUnits','points',... 
'FontWeight','normal',... 
'FontSize',11,... 
'FontName','Times')
ylim([0 160])
xlim([-0.20 0.10])


title('Distribution of returns',...
    'FontUnits','points',... 
    'FontWeight','normal',... 
    'FontSize',titlesize,... 
    'FontName','Times',...
    'interpreter','latex');

print -deps fig_histfit_s2.eps
matlab2tikz('fig_histfit_s2.tex','width','\figW','height','\figH','showInfo', false);

%%
toc