% compute long-run correlation relying on the method in Muller and Watson (2018)
% TFP growth; Labor hours growth; Top 10% income share; Yield spread

clc
clear all
close all

%% Total factor productivity

load("rsatr_Ymean_LW_rr2.mat");

US_factor = readtable("data/annual/US_MFPNFBS.csv");
MFP = readtable("data/annual/MFP.xlsx");
y1 = diff(log(US_factor{24:74,2}));  % 1972-2021 TFP growth rate
y2 = diff(log(MFP.UK));  % 1986-2021
y3 = diff(log(MFP.GM));  % 1986-2021

T1 = length(y1);  % 1972-2021
T = length(y2);   % 1985-2021

parmdir = 'lrcov_parms_dir/';        % Location of various pre-computed parameters and covariance matrices (these are supplied with the replication files)
progdir = 'matlab_lrcov_programs/';
p = path; % Change path so that Matlab accesses utility programs
path(progdir,p);

lowest_period_years = 10;  % This is low-period cut-off .. here in years
q1 = floor(2*T1/lowest_period_years);
q = floor(2*T/lowest_period_years);

x1 = rstars.US;   % 1972-2021
x2 = rstars.UK(15:end);  % 1986-2021
x3 = rstars.GM(15:end);  % 1986-2021
 
[corr1,bet,sig,mis] = lrcov(x1,y1,q1,1,0,0,parmdir);
[corr2,bet,sig,mis] = lrcov(x2,y2,q,1,0,0,parmdir);
[corr3,bet,sig,mis] = lrcov(x3,y3,q,1,0,0,parmdir);

%% Labor hours

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


y1 = diff(log(us_hours.Value(7:209))); % 1971Q4-2022Q1, labor hour growth
y2 = diff(log(uk_hours.hours(6:125))); % 1992Q3-2022Q1

T1 = length(y1);
T2 = length(y2);
load("rsatr_Qend_LW.mat");
x1 = rstars.US;
x2 = rstars.UK(84:end);


y3 = diff(log(gm_hours.Value/52));      % 1992-2021,labor hour growth
T3 = length(y3);   
load("rsatr_Ymean_LW.mat");
x3 = rstars.GM(21:end);

parmdir = 'lrcov_parms_dir/';        % Location of various pre-computed parameters and covariance matrices (these are supplied with the replication files)
progdir = 'matlab_lrcov_programs/';
p = path; % Change path so that Matlab accesses utility programs
path(progdir,p);

lowest_period_years = 10;  % This is low-period cut-off .. here in years
nper = 4; 
lowest_period = nper*lowest_period_years;
q1 = floor(2*T1/lowest_period);
q2 = floor(2*T2/lowest_period);
q3 = floor(2*T3/lowest_period_years);

[corr1,bet,sig,mis] = lrcov(x1,y1,q1,1,0,0,parmdir);
[corr2,bet,sig,mis] = lrcov(x2,y2,q2,1,0,0,parmdir);
[corr3,bet,sig,mis] = lrcov(x3,y3,q3,1,0,0,parmdir);


%% Income inequality 
 
load("rsatr_Ymean_LW.mat");

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

T1 = 62-12;  % 1972-2021
T2 = 62-23;  % 1983-2021
T3 = 62-20;  % 1980-2021

parmdir = 'lrcov_parms_dir/';        % Location of various pre-computed parameters and covariance matrices (these are supplied with the replication files)
progdir = 'matlab_lrcov_programs/';

p = path; % Change path so that Matlab accesses utility programs
path(progdir,p);

lowest_period_years = 10;  % This is low-period cut-off .. here in years
q1 = floor(2*T1/lowest_period_years);
q2 = floor(2*T2/lowest_period_years);
q3 = floor(2*T3/lowest_period_years);

y1 = Income_10.US(13:end);
y2 = Income_10.UK(24:end);
y3 = Income_10.GM(21:end);
ADF_Individual(diff(y1));  % Income_inequality is I(1)
ADF_Individual(diff(y2));
ADF_Individual(diff(y3));

x1 = rstars.US;
x2 = rstars.UK(12:end);
x3 = rstars.GM(9:end);

[corr1,bet,sig,mis] = lrcov(x1,y1,q1,1,0,0,parmdir);
[corr2,bet,sig,mis] = lrcov(x2,y2,q2,1,0,0,parmdir);
[corr3,bet,sig,mis] = lrcov(x3,y3,q3,1,0,0,parmdir);


y1 = Income_1.US(13:end);
y2 = Income_1.UK(24:end);
y3 = Income_1.GM(21:end);

[corr1,bet,sig,mis] = lrcov(x1,y1,q1,1,0,0,parmdir);
[corr2,bet,sig,mis] = lrcov(x2,y2,q2,1,0,0,parmdir);
[corr3,bet,sig,mis] = lrcov(x3,y3,q3,1,0,0,parmdir);

%% Convenience Yield:   r* and yield spread
% End of quarter r*
load("rsatr_Qend_LW_RR2.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;

% ADF test using the code from Stock and Watson (2012) Disentangling the
% channel of the 2007-2009 recession
ADF_Individual(diff(rstars.US));
ADF_Individual(diff(rstars.UK(46:end)));


parmdir = 'lrcov_parms_dir/';        % Location of various pre-computed parameters and covariance matrices (these are supplied with the replication files)
progdir = 'matlab_lrcov_programs/';

p = path; % Change path so that Matlab accesses utility programs
path(progdir,p);

lowest_period_years = 10;  % This is low-period cut-off .. here in years
nper = 4; 

T_obs = 145; % 1986Q1 - 2022Q1
lowest_period = nper*lowest_period_years;
q = floor(2*T_obs/lowest_period);
lowest_period_years_q = 2*T_obs/(nper*q);
highest_period_years_q = 2*T_obs/(nper);

x1 = rstars.US(58:202);
x2 = rstars.UK(58:202);
x3 = rstars.GM(58:202);

y1 = Yields_data.usspread(105:249);
y2 = Yields_data.ukspread(105:249);
y3 = Yields_data.gmspread(105:249);


% call the lrcov.m function to compute long-run correlation
[corr1,bet,sig,mis] = lrcov(x1,y1,q,1,0,0,parmdir);
[corr2,bet,sig,mis] = lrcov(x2,y2,q,1,0,0,parmdir);
[corr3,bet,sig,mis] = lrcov(x3,y3,q,1,0,0,parmdir);

path(p);  % Restore original path






