%% Import data from spreadsheet
% Script for importing data from the following spreadsheet:
%
%    Workbook: /Users/naims/Desktop/Thesis/data/corrdata.xlsx
%    Worksheet: EWJ JAPAN
%
% Auto-generated by MATLAB on 05-Jun-2019 21:20:44

%% Setup the Import Options
opts = spreadsheetImportOptions("NumVariables", 4);

% Specify sheet and range
opts.Sheet = "EWJ JAPAN";
opts.DataRange = "A2:D5790";

% Specify column names and types
opts.VariableNames = ["Date", "LastPrice", "NAVpershare", "sp"];
opts.SelectedVariableNames = ["Date", "LastPrice", "NAVpershare", "sp"];
opts.VariableTypes = ["datetime", "double", "double", "double"];
opts = setvaropts(opts, 1, "InputFormat", "");

% Import the data
corrdata = readtable("/Users/naims/Desktop/Thesis/data/corrdata.xlsx", opts, "UseExcel", false);


%% Clear temporary variables
clear opts
%%
m = table2array(corrdata(:,1));
y = datetime(m,'ConvertFrom','dd-mm-yyyy');
%% Calculating log prices
JAP = table2array(corrdata(:,2:end));
JAP(:,4) = yyyymmdd(y);
JAP(:,5) = log(JAP(:,1)); %etf price
JAP(:,6) = log(JAP(:,2)); %nav price
JAP(:,7) = log(JAP(:,3)); %s&p price
n = length(JAP);


%% Calculating Price deviation between ETF and NAV price
dev = zeros(n,1);
for i = 1:n
    dev(i,1) = abs((JAP(i,1)-JAP(i,2)))/JAP(i,2);
end
% calculating mean and stdev for a monthly horizon
for i = 1:n-21
  md(i,1) = mean(dev(i:i+21,1));
  md(i,2) = std(dev(i:i+21,1));
end
%% log returns of market price, nav and s&p 
dret=zeros(n-1,3); %daily
wret=zeros(n-5,3); %weekly
mret=zeros(n-21,3); %monthly
aret=zeros(n-252,3); %quarterly
%% daily returns
for j = 1:n-1
    dret(j,1)=JAP(j+1,5)-JAP(j,5); %etf return
    dret(j,2)=JAP(j+1,6)-JAP(j,6); %nav return
    dret(j,3)=JAP(j+1,7)-JAP(j,7); %s&p return
    
end
% weekly 
for j = 1:n-5
    wret(j,1)=JAP(j+5,5)-JAP(j,5); %etf return
    wret(j,2)=JAP(j+5,6)-JAP(j,6); %nav return
    wret(j,3)=JAP(j+5,7)-JAP(j,7); %s&p return
end
% monthly
for j = 1:n-22
    mret(j,1)=JAP(j+22,5)-JAP(j,5); %etf return
    mret(j,2)=JAP(j+22,6)-JAP(j,6); %nav return
    mret(j,3)=JAP(j+22,7)-JAP(j,7); %s&p return
end
% Quarterly
for j = 1:n-65
    aret(j,1)=JAP(j+65,5)-JAP(j,5); %etf return
    aret(j,2)=JAP(j+65,6)-JAP(j,6); %nav return
    aret(j,3)=JAP(j+65,7)-JAP(j,7); %s&p return
end


%% descriptive statistics of the monthalized ETF and NAV returns

DS(1,1) = "Japan";
DS(1,2) = "ETF returns";
DS(2,1) = "Mean";
DS(2,2) = mean(mret(:,1));
DS(3,1) = "maxvalue";
DS(3,2) = max(mret(:,1));
DS(4,1) = "minvalue";
DS(4,2) = min(mret(:,1));
DS(5,1) = "standard deviation";
DS(5,2) = std(mret(:,1));
DS(6,1) = "kurtosis";
DS(6,2) = kurtosis(mret(:,1))-3;
DS(7,1) = "skewness";
DS(7,2) = skewness(mret(:,1));
DS(8,1) = "median";
DS(8,2) = median(mret(:,1));


DS(1,1) = "Japan";
DS(1,3) = "NAV returns";
DS(2,3) = "Mean";
DS(2,4) = mean(mret(:,2));
DS(3,3) = "maxvalue";
DS(3,4) = max(mret(:,2));
DS(4,3) = "minvalue";
DS(4,4) = min(mret(:,2));
DS(5,3) = "standard deviation";
DS(5,4) = std(mret(:,2));
DS(6,3) = "kurtosis";
DS(6,4) = kurtosis(mret(:,2));
DS(7,3) = "skewness";
DS(7,4) = skewness(mret(:,2));
DS(8,3) = "median";
DS(8,4) = median(mret(:,2));


% S&P returns
DS(1,5) = "S&P returns";
DS(2,5) = "Mean";
DS(2,6) = mean(mret(:,3));
DS(3,5) = "maxvalue";
DS(3,6) = max(mret(:,3));
DS(4,5) = "minvalue";
DS(4,6) = min(mret(:,3));
DS(5,5) = "standard deviation";
DS(5,6) = std(mret(:,3));
DS(6,5) = "kurtosis";
DS(6,6) = kurtosis(mret(:,3));
DS(7,5) = "skewness";
DS(7,6) = skewness(mret(:,3));
DS(8,5) = "median";
DS(8,6) = median(mret(:,3));

%%
%Correlations at multiple horizons for 252 day rolling period returns

% Daily
n = length(dret);
for i=1:n-252
     corrnav(i,1) = corr(dret(i:i+252,2),dret(i:i+252,3));
     corretf(i,1) = corr(dret(i:i+252,1),dret(i:i+252,3));
     corrnaet(i,1) = corr(dret(i:i+252,1),dret(i:i+252,2));
end

% Weekly
n = length(wret);
for i=1:n-252
    corrnav(i,2) = corr(wret(i:i+252,2),wret(i:i+252,3));
    corretf(i,2) = corr(wret(i:i+252,1),wret(i:i+252,3));
    corrnaet(i,2) = corr(wret(i:i+252,1),wret(i:i+252,2));
end

% Monthly
n = length(mret);
for i=1:n-252
    corrnav(i,3) = corr(mret(i:i+252,2),mret(i:i+252,3));
    corretf(i,3) = corr(mret(i:i+252,1),mret(i:i+252,3));
    corrnaet(i,3) = corr(mret(i:i+252,1),mret(i:i+252,2));
end


%% AR(1) process to fit nav and s&p returns
asynret = dret(1:end,2:3); %asynchronous data as input
ToEstMdl = varm(2,4);
ToEstMdl.Constant = 0;
er = estimate(ToEstMdl,asynret);
summarize(er);
eps = infer(er,asynret);

%%
asynret0 = asynret(1:end-1,:);
asynret1 = asynret(2:end,:);
retdiff = asynret1-asynret0; %error term for the synchronization


%% Calculate synchronous returns from AR(1) parameters
a1 = [-0.086357, 0.4778 ; 0, 0]; %AR(1) parameters adjusted based on tstats and region
w = a1*transpose(retdiff);
synret = asynret1 + transpose(w); % synchronous NAV returns
corr(synret(:,1),synret(:,2))

%% DCC GARCH modelling using kevin sheppards tool box
[PARAMETERS,LL,HT,VCV,SCORES] = dcc(synret,[],1,0,1); %synchronous returns
[PARAMETERS1,LL1,HT1,VCV1,SCORES1] = dcc(asynret,[],1,0,1); %asynchronous returns

%% PLotting and tabulating correlations at multiple horizons
%%
for i = 2:n-1
    wkret(i,1) = ((dret(i,1)+1)^5)-1;
    mnret(i,1) = ((dret(i,1)+1)^22)-1;
    qret(i,1) = ((dret(i,1)+1)^65)-1;
    wkret(i,2) = ((dret(i,2)+1)^5)-1;
    mnret(i,2) = ((dret(i,2)+1)^22)-1;
    qret(i,2) = ((dret(i,2)+1)^65)-1;
    wkret(i,3) = ((dret(i,3)+1)^5)-1;
    mnret(i,3) = ((dret(i,3)+1)^22)-1;
    qret(i,3) = ((dret(i,3)+1)^65)-1;
end
%% correlations
c(1,1) = "nav - sp";
c(2,1) = corr(dret(:,2),dret(:,3));
c(3,1) = corr(wret(:,2),wret(:,3));
c(4,1) = corr(mret(:,2),mret(:,3));
c(5,1) = corr(aret(:,2),aret(:,3));

c(1,3) = "nav - etf";
c(2,3) = corr(dret(:,2),dret(:,1));
c(3,3) = corr(wret(:,2),wret(:,1));
c(4,3) = corr(mret(:,2),mret(:,1));
c(5,3) = corr(aret(:,2),aret(:,1));

c(1,5) = "sp - etf";
c(2,5) = corr(dret(:,1),dret(:,3));
c(3,5) = corr(wret(:,1),wret(:,3));
c(4,5) = corr(mret(:,1),mret(:,3));
c(5,5) = corr(aret(:,1),aret(:,3));

% Asynchronous NAV returns daily
c(1,4) = "snav - etf";
c(2,4) = corr(synret(:,1),dret(2:end,1)); 
c(1,2) = "snav - sp";
c(2,2) = corr(synret(:,1),synret(:,2));
% Lagged S&P returns
c(1,6) = "spt-1 - nav";
c(2,6) = corr(dret(1:end-1,3),dret(2:end,2));



%% regression of S&P returns on NAV returns 
Y = abs(dret(:,1)-dret(:,2));
X = abs(dret(:,3));
mdl = fitlm(X(3:end),Y(3:end)); % NAV at t and S&P at t
mdl1 = fitlm(X(2:end-1),Y(3:end)); % NAV at t and S&P at t-1
mdl2 = fitlm(X(1:end-2),Y(3:end));% NAV at t and S&P at t-1

%% regression with all varibales 

%% Setup the Import Options
opts = spreadsheetImportOptions("NumVariables", 11);

% Specify sheet and range
opts.Sheet = "EWJ";
opts.DataRange = "A2:K5790";

% Specify column names and types
opts.VariableNames = ["Date", "VIX", "FX", "SP500", "Volume", "Volatility", "Bidask", "AUM", "CR", "UScrisis", "TDF"];
opts.SelectedVariableNames = ["Date", "VIX", "FX", "SP500", "Volume", "Volatility", "Bidask", "AUM", "CR", "UScrisis", "TDF"];
opts.VariableTypes = ["datetime", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
opts = setvaropts(opts, 1, "InputFormat", "");

% Import the data
regression = readtable("/Users/nataliiapetromanova/Documents/My Documents/BI/Spring 2019/Master Thesis /regressiondata.xlsx", opts, "UseExcel", false);
%% Clear temporary variables
clear opts
data1 = table2array(regression(:,2:end));
n = length(data1);
%% creating table of variables 

% 1) VIX 
v(:,1) = data1(:,1)/100;

% 2) daily return on FX
FX = log(data1(:,2));
for i = 2:n  
   v(i,2) = abs(FX(i,1)-FX(i-1,1));
end

% 3) log SP500
SP = log(data1(:,3));

for j = 2:n
    v(j,3) = abs(SP(j,1) - SP(j-1,1));
end

% 4) log average daily volume
vol = log (data1(:,4));
for j = 2:n
    v(j,4) = abs(vol(j,1) - vol(j-1,1));
end

% 5) volatility of ETF's daily average trading prices
 v(:,5)= data1(:,5);
 
% 6) bidask spread
 v(:,6) = data1(:,6);

% 7) AUM
size = log(data1(:,7));
for j = 2:n
    v(j,7) = abs(size(j,1) - size(j-1,1));
end

% 8) creation redempion

v(:,8) = data1(:,8);

% 9) US crisis dummy variable

v(:,9) = data1(:,9);

%10) tdf 
 v(:,10) = data1(:,10); 

 %%
tb = array2table(tt,'VariableNames', {'VIX' 'FX' 'SP500' 'VOLUME' 'VOLATILITY' 'Bidask' 'AUM' 'CRP' 'UScrisis' 'TDF'});
%%
fitlm(tb);
disp(ans.Rsquared)
disp(ans.Coefficients)

%%