%% Panel Regression
% TFP growth; Labor hours growth; Top 10% income share; Yield spread

clc
clear all
close all

%% Productivity growth
load("rsatr_Ymean.mat");

US_factor = readtable("data/annual/US_MFPNFBS.csv");
y1 = US_factor{38:74,2};
y1 = diff(log(y1))*100;   % 1986-2021
MFP = readtable("data/annual/MFP.xlsx");
y2 = MFP.UK;
y3 = MFP.GM;
y2 = diff(log(y2))*100;
y3 = diff(log(y3))*100;
Xs = [y1 y2 y3];
% myy = rstars{15:50,2:4};
myy = diff(rstars{14:50,2:4},1);

n = size(myy,1);
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(3), Xs(i,:)'];
end
[b_tfp,~,~,CovB_tfp] = mvregress(myx, myy, 'algorithm','cwls');
t_value = b_tfp(4)/sqrt(CovB_tfp(4,4));    
p_value = 2*(1-tcdf(abs(t_value),n-2))

% tvalue = b_tfp(4)/sqrt(CovB_tfp(4,4)) ;
% tboots_tfp = panel_bootstrap(myy, Xs, 50000, 616);

%% Labor hours growth
us_hours = readtable("data/annual/US_workhours.csv");
uk_hours = readtable("data/annual/UK_workhours.csv");
gm_hours = readtable("data/annual/GM_workhours.csv");

load("rsatr_Qend.mat");
% load("rsatr_Ymean.mat");

y1 = us_hours.Value(90:209);   
y2 = uk_hours.hours(6:125);    
y1 = diff(log(y1))*400;            % 1992Q3-2022Q1
y2 = diff(log(y2))*400;            % 1992Q3-2022Q1
y3 = diff(log(gm_hours.Value));    % 1992-2021
load("rsatr_Ymean.mat");


% myy = rstars{84:end,2:3};
myy = diff(rstars{83:end,2:3});
Xs = [y1 y2];
n = size(myy,1);
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(2), Xs(i,:)'];
end
[b_hour,~,~,CovB_hour] = mvregress(myx, myy, 'algorithm','cwls');
t_value = b_hour(3)/sqrt(CovB_hour(3,3)) ;
p_value = 2*(1-tcdf(abs(t_value),n-2))

% tboots_hour = panel_bootstrap(myy, Xs, 50000, 616);
% findprctile(tvalue,tboots_hour) 

%%  Income inequality, top 10% income share 
load("rsatr_Ymean.mat");
% load("rsatr_Yend.mat");

% Top 10% Income
Income_10 = readtable("data/annual/income_top10.xlsx");
% Top 1% Income
Income_1 = readtable("data/annual/income_top1.xlsx");

% myy = rstars{12:end,2:4};
myy = diff(rstars{12:end,2:4},1);
Xs = diff(Income_10{24:end,3:5}*100);
n = size(myy,1);
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(3), Xs(i,:)'];
end
[b_income,~,~,CovB_income] = mvregress(myx, myy, 'algorithm','cwls');
t_value = b_income(4)/sqrt(CovB_income(4,4)) ;
p_value = 2*(1-tcdf(abs(t_value),n-2))

% tboots_income = panel_bootstrap(myy, Xs, 50000, 616);
% findprctile(tvalue,tboots_income) 

%%  Convenien yield: Yield Spread
% End of quarter r*
load("rsatr_Qend.mat");

% End of quarter yield spreads
Yields_data = readtable("./data/annual/yields.xlsx");
Yields_data.usspread = Yields_data.DBAA - Yields_data.US;
Yields_data.ukspread = Yields_data.DBAA - Yields_data.UK;
Yields_data.gmspread = Yields_data.DBAA - Yields_data.GM;

% y1 = rstars.US(58:202);
% y2 = rstars.UK(58:202);
% y3 = rstars.GM(58:202);
y1 = diff(rstars.US(57:202));
y2 = diff(rstars.UK(57:202));
y3 = diff(rstars.GM(57:202));
myy = [y1 y2 y3];
Xs = [Yields_data.usspread(105:249) Yields_data.ukspread(105:249) Yields_data.gmspread(105:249)];   
n = size(myy,1);
myx = cell(n,1);
for i = 1:n
    myx{i} = [eye(3), Xs(i,:)'];
end
[b_spread,~,~,CovB_spread] = mvregress(myx, myy, 'algorithm','cwls');
t_value = b_spread(4)/sqrt(CovB_spread(4,4)) ;
p_value = 2*(1-tcdf(abs(t_value),n-2))

% tboots_spread = panel_bootstrap(myy, Xs, 50000, 616);
% findprctile(tvalue,tboots_spread)









