% compute long-run correlation relying on the method in Muller and Watson (2018)
% TFP growth; Labor hours growth; Top 10% income share; Yield spread
% rstars from HLW downloaded from Holston's website

clc
clear all
close all

HLWdata = readtable("./data/rstar_world.xlsx");

% annualized rstars 1972 - 2019
% use euro rstar as a proxy of GM rstar
tempus = HLWdata.US(HLWdata.Date>datetime(1971,12,1) & HLWdata.Date<datetime(2020,1,1),:);
tempuk = HLWdata.UK(HLWdata.Date>datetime(1971,12,1) & HLWdata.Date<datetime(2020,1,1),:);
tempeuro = HLWdata.Euro(HLWdata.Date>datetime(1971,12,1) & HLWdata.Date<datetime(2020,1,1),:);

n = size(tempus,1)/4;
usr = nan(n,1);  
ukr = nan(n,1);  
euror = nan(n,1); 

for i = 1:length(tempus)
    if mod(i,4)==0
        usr(i/4) = mean(tempus(i-3:i));
        ukr(i/4) = mean(tempuk(i-3:i));
        euror(i/4) = mean(tempeuro(i-3:i));
    end
end

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

T1 = 62-14;  % 1972-2019
T2 = 62-22;  % 1980-2019
T3 = 62-22;  % 1980-2019

x1 = usr;
x2 = ukr(9:end);
x3 = euror(9: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
q1 = floor(2*T1/lowest_period_years);
q2 = floor(2*T2/lowest_period_years);
q3 = floor(2*T3/lowest_period_years);

% 10-percent inequality
y1 = Income_10.US(13:end-2);
y2 = Income_10.UK(21:end-2);
y3 = Income_10.GM(21:end-2);
% ADF_Individual(diff(y1));  % Income_inequality is I(1)
% ADF_Individual(diff(y2));
% ADF_Individual(diff(y3));
[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);

% 1-percent inequality
y1 = Income_1.US(13:end-2);
y2 = Income_1.UK(21:end-2);
y3 = Income_1.GM(21:end-2);
[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);




path(p);  % Restore original path

