clear all
global y t;
OPTIONS = optimset('MaxFunEvals',600,'Maxiter',200,'TolFun',1e-12);

anp = ' & ';
dslash = ' \\ ';
hline = ' \hline ';
bf = ' \bf ';
leftp = '(';
rightp = ')';

% 1978 to 2005 data
dataRGDP = dlmread('vintage_growth_currqtrly.csv',',');
dataRGDPa = dlmread('vintage_growth_advcurrqtrly.csv',',');
dataRGDI = dlmread('vintage_realGDI_growth_currqtrly.csv',',');
dataRGDIa = dlmread('vintage_realGDI_predgrowth_latest.csv',',');


RGDP_adv = dataRGDPa(:,5);
RGDI_adv = dataRGDIa(:,8);
RGDP_fin = dataRGDP(:,7);
RGDI_fin = dataRGDI(:,5);
RGDP_latest = dataRGDP(:,8);
RGDI_latest = dataRGDI(:,6);

tstart = 3;
tend = 127;

dataUR = dlmread('UR.csv',',');
t_u = size(dataUR,1);

UR_diff = 4*(dataUR(53:t_u,2)-dataUR(52:t_u-1,2));

clear y1 n control1 x
'UR_diff on initial RGDP, RGDI 1978-present, intercept'
y1 = UR_diff(tstart+68:tend+68,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(1,:) = [beta' rsq];
out(2,:) = [se_newey' 0];

' UR_diff on initial RGDP(t-1), RGDI(t-1) 1978-present'
x = [RGDP_adv(tstart-1:tend-1,:) RGDI_adv(tstart-1:tend-1,:) control1];
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(3,:) = [beta' rsq];
out(4,:) = [se_newey' 0];


' UR_diff on initial RGDP(t-2), RGDI(t-2) 1978-present'
x = [RGDP_adv(tstart-2:tend-2,:) RGDI_adv(tstart-2:tend-2,:) control1];
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(5,:) = [beta' rsq];
out(6,:) = [se_newey' 0];



dataLFCE = dlmread('lfce_fixbreaks.csv',',');
t_u = size(dataLFCE,1)-1;

lfce_diff = (((dataLFCE(53:t_u,2)./dataLFCE(52:t_u-1,2)).^4)-1)*100;

clear y1 n control1 x
' HH hours growth on initial RGDP, RGDI 1978-present, intercept'
y1 = lfce_diff(tstart+68:tend+68,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(7,:) = [beta' rsq];
out(8,:) = [se_newey' 0];

' HH hours growth on initial RGDP(t-1), RGDI(t-1) 1978-present'
x = [RGDP_adv(tstart-1:tend-1,:) RGDI_adv(tstart-1:tend-1,:) control1];
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(9,:) = [beta' rsq];
out(10,:) = [se_newey' 0];


' HH hours growth on initial RGDP(t-2), RGDI(t-2) 1978-present'
x = [RGDP_adv(tstart-2:tend-2,:) RGDI_adv(tstart-2:tend-2,:) control1];
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(11,:) = [beta' rsq];
out(12,:) = [se_newey' 0];




dataISM = dlmread('NBER_ISM.csv',',');
t_u = size(dataISM,1)-1;

ISM = dataISM(53:t_u,3);
[dataISM(53:t_u,1)/1000 ISM]

clear y1 n control1 x
' ISM on initial RGDP, RGDI 1978-present, intercept'
y1 = ISM(tstart+68:tend+68,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(13,:) = [beta' rsq];
out(14,:) = [se_newey' 0];

' ISM on initial RGDP(t-1), RGDI(t-1) 1978-present'
x = [RGDP_adv(tstart-1:tend-1,:) RGDI_adv(tstart-1:tend-1,:) control1];
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(15,:) = [beta' rsq];
out(16,:) = [se_newey' 0];

' ISM on initial RGDP(t-2), RGDI(t-2) 1978-present'
x = [RGDP_adv(tstart-2:tend-2,:) RGDI_adv(tstart-2:tend-2,:) control1];
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(17,:) = [beta' rsq];
out(18,:) = [se_newey' 0];



data_assets = dlmread('stock_price_changes_qtrly_52q3.csv',',');
t_assets = size(data_assets,1);
%sp500 = data_assets(53:t_assets,2);
sp500_diff1 = data_assets(53:t_assets,3);
sp500_diff3 = (data_assets(53:t_assets,3)+data_assets(52:t_assets-1,3))./2;
sp500_diff3 = (data_assets(53:t_assets,3)+data_assets(52:t_assets-1,3)+data_assets(51:t_assets-2,3))./3;
sp500_diff4 = (data_assets(53:t_assets,3)+data_assets(52:t_assets-1,3)+data_assets(51:t_assets-2,3)+data_assets(50:t_assets-3,3))./4;
sp500_diff5 = (data_assets(53:t_assets,3)+data_assets(52:t_assets-1,3)+data_assets(51:t_assets-2,3)+data_assets(50:t_assets-3,3)+data_assets(49:t_assets-4,3))./5;
sp500_diff6 = (data_assets(53:t_assets,3)+data_assets(52:t_assets-1,3)+data_assets(51:t_assets-2,3)+data_assets(50:t_assets-3,3)+data_assets(49:t_assets-4,3)+data_assets(48:t_assets-5,3))./6;
sp500_diff7 = (data_assets(53:t_assets,3)+data_assets(52:t_assets-1,3)+data_assets(51:t_assets-2,3)+data_assets(50:t_assets-3,3)+data_assets(49:t_assets-4,3)+data_assets(48:t_assets-5,3)+data_assets(47:t_assets-6,3))./7;
term_2yr = data_assets(53:t_assets,5) - data_assets(53:t_assets,4);
term_2yr_4 = data_assets(49:t_assets-4,5) - data_assets(49:t_assets-4,4);
term_2yr_5 = data_assets(48:t_assets-5,5) - data_assets(48:t_assets-5,4);
term_2yr_6 = data_assets(47:t_assets-6,5) - data_assets(47:t_assets-6,4);
term_2yr_7 = data_assets(46:t_assets-7,5) - data_assets(46:t_assets-7,4);
term_2yr_8 = data_assets(45:t_assets-8,5) - data_assets(45:t_assets-8,4);
term_2yr_9 = data_assets(44:t_assets-9,5) - data_assets(44:t_assets-9,4);
t_x = size(sp500_diff1,1);
[data_assets(53:t_assets,1)/1000 sp500_diff1 sp500_diff4 sp500_diff7 term_2yr term_2yr_4 term_2yr_8]

'SP500 on initial RGDP, RGDI 1978-present, intercept'
y1 = sp500_diff4(tstart+68:tend+68,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(19,:) = [beta' rsq];
out(20,:) = [se_newey' 0];



'Term Spread on initial RGDP, RGDI 1978-present, intercept'
y1 = term_2yr_8(tstart+68:tend+68,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(21,:) = [beta' rsq];
out(22,:) = [se_newey' 0];


data_assets = dlmread('stock_price_changes_qtrly_52q3.csv',',');
t_assets = size(data_assets,1);
spf_currqtr = data_assets(103:t_assets,6);
spf_1qtr_ahead = data_assets(102:t_assets-1,7);
spf_2qtr_ahead = data_assets(101:t_assets-2,8);
spf_3qtr_ahead = data_assets(100:t_assets-3,9);
[data_assets(103:t_assets,1)/1000 spf_currqtr spf_1qtr_ahead spf_2qtr_ahead spf_3qtr_ahead]

'SPF, current qtr on initial RGDP, RGDI 1978-present, intercept'
y1 = spf_currqtr(tstart+18:tend+18,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(23,:) = [beta' rsq];
out(24,:) = [se_newey' 0];

'SPF, 1 qtr ahead on initial RGDP, RGDI 1978-present, intercept'
y1 = spf_1qtr_ahead(tstart+18:tend+18,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(25,:) = [beta' rsq];
out(26,:) = [se_newey' 0];

'SPF, 2 qtr ahead on initial RGDP, RGDI 1978-present, intercept'
y1 = spf_2qtr_ahead(tstart+18:tend+18,:)
n = size(y1,1)
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDI_adv(tstart:tend,:) control1]
kv = size(x,2);
beta = inv(x'*x)*x'*y1;
pred = x*beta; 
e = y1 - x*beta; 
ss = y1'*(eye(n) - control1*inv(control1'*control1)*control1')*y1;
rsq = 1 - ((e'*e/(n-kv))/(ss/(n-1)))
rmse = (e'*e/(n-kv))^(0.5);

% Computing Newey-West Standard Errors for beta
 % Calculating main diagonal of var-cov. matrix, correcting for cross correlation
   P0 = x'*diag(diag(e*e'))*x;
 % Calculating off-diagonals of var-cov. matrix
   h_1 = diag(diag(e*e',1)',1);
   h_1_ = diag(diag(e*e',-1)',-1);
   P1a = x'*h_1*x;
   P1b = x'*h_1_*x;

   h_2 = diag(diag(e*e',2)',2);
   h_2_ = diag(diag(e*e',-2)',-2);
   P2a = x'*h_2*x;
   P2b = x'*h_2_*x;

   h_3 = diag(diag(e*e',3)',3);
   h_3_ = diag(diag(e*e',-3)',-3);
   P3a = x'*h_3*x;
   P3b = x'*h_3_*x;

   h_4 = diag(diag(e*e',4)',4);
   h_4_ = diag(diag(e*e',-4)',-4);
   P4a = x'*h_4*x;
   P4b = x'*h_4_*x;

   h_5 = diag(diag(e*e',5)',5);
   h_5_ = diag(diag(e*e',-5)',-5);
   P5a = x'*h_5*x;
   P5b = x'*h_5_*x;

   h_6 = diag(diag(e*e',6)',6);
   h_6_ = diag(diag(e*e',-6)',-6);
   P6a = x'*h_6*x;
   P6b = x'*h_6_*x;

   h_7 = diag(diag(e*e',7)',7);
   h_7_ = diag(diag(e*e',-7)',-7);
   P7a = x'*h_7*x;
   P7b = x'*h_7_*x;
    
   h_8 = diag(diag(e*e',8)',8);
   h_8_ = diag(diag(e*e',-8)',-8);
   P8a = x'*h_8*x;
   P8b = x'*h_8_*x;

omega = P0 + ((8/9)*P1a) + ((8/9)*P1b) + ((7/9)*P2a) + ((7/9)*P2b) + ((6/9)*P3a) + ((6/9)*P3b) + ((5/9)*P4a) + ((5/9)*P4b) + ((4/9)*P5a) + ((4/9)*P5b) + ((3/9)*P6a) + ((3/9)*P6b) + ((2/9)*P7a) + ((2/9)*P7b) + ((1/9)*P8a) + ((1/9)*P8b);
asyvarc = inv(x'*x)*omega*inv(x'*x);
se_newey = diag(asyvarc).^(1/2);
asyvar = inv(x'*x)*P0*inv(x'*x);
se_white = diag(asyvar).^(1/2);
se_ols = (diag(sum(e'*e)/(n-kv)*inv(x'*x))).^(1/2);
[beta'; se_ols'; se_white'; se_newey']
plot(y1,x(:,1),'x')

out(27,:) = [beta' rsq];
out(28,:) = [se_newey' 0];


out_row1 = ['$\left(UR_t - UR_{t-1}\right)*4$' anp bf num2str(out(1,1),'%3.2f') anp anp anp bf num2str(out(1,2),'%3.2f') anp anp anp bf num2str(out(1,3),'%3.2f') anp num2str(out(1,4),'%3.2f') dslash anp leftp num2str(out(2,1),'%3.2f') rightp anp anp anp leftp num2str(out(2,2),'%3.2f') rightp anp anp anp leftp num2str(out(2,3),'%3.2f') rightp anp dslash]
out_row2 = [ anp anp bf num2str(out(3,1),'%3.2f') anp anp anp bf num2str(out(3,2),'%3.2f') anp anp bf num2str(out(3,3),'%3.2f') anp num2str(out(3,4),'%3.2f') dslash anp anp leftp num2str(out(4,1),'%3.2f') rightp anp anp anp leftp num2str(out(4,2),'%3.2f') rightp anp anp leftp num2str(out(4,3),'%3.2f') rightp anp dslash]
out_row3 = [ anp anp anp bf num2str(out(5,1),'%3.2f') anp anp anp bf num2str(out(5,2),'%3.2f') anp bf num2str(out(5,3),'%3.2f') anp num2str(out(5,4),'%3.2f') dslash anp anp anp leftp num2str(out(6,1),'%3.2f') rightp anp anp anp leftp num2str(out(6,2),'%3.2f') rightp anp leftp num2str(out(6,3),'%3.2f') rightp anp dslash]
out_row4 = ['$\left(E_t^{household}/E_{t-1}^{household}\right)^4$' anp bf num2str(out(7,1),'%3.2f') anp anp anp bf num2str(out(7,2),'%3.2f') anp anp anp bf num2str(out(7,3),'%3.2f') anp num2str(out(7,4),'%3.2f') dslash anp leftp num2str(out(8,1),'%3.2f') rightp anp anp anp leftp num2str(out(8,2),'%3.2f') rightp anp anp anp leftp num2str(out(8,3),'%3.2f') rightp anp dslash]
out_row5 = [ anp anp bf num2str(out(9,1),'%3.2f') anp anp anp bf num2str(out(9,2),'%3.2f') anp anp bf num2str(out(9,3),'%3.2f') anp num2str(out(9,4),'%3.2f') dslash anp anp leftp num2str(out(10,1),'%3.2f') rightp anp anp anp leftp num2str(out(10,2),'%3.2f') rightp anp anp leftp num2str(out(10,3),'%3.2f') rightp anp dslash]
out_row6 = [ anp anp anp bf num2str(out(11,1),'%3.2f') anp anp anp bf num2str(out(11,2),'%3.2f') anp bf num2str(out(11,3),'%3.2f') anp num2str(out(11,4),'%3.2f') dslash anp anp anp leftp num2str(out(12,1),'%3.2f') rightp anp anp anp leftp num2str(out(12,2),'%3.2f') rightp anp leftp num2str(out(12,3),'%3.2f') rightp anp dslash]
out_row7 = ['$ISM_t^{manuf.}$' anp bf num2str(out(13,1),'%3.2f') anp anp anp bf num2str(out(13,2),'%3.2f') anp anp anp bf num2str(out(13,3),'%3.2f') anp num2str(out(13,4),'%3.2f') dslash anp leftp num2str(out(14,1),'%3.2f') rightp anp anp anp leftp num2str(out(14,2),'%3.2f') rightp anp anp anp leftp num2str(out(14,3),'%3.2f') rightp anp dslash]
out_row8 = [ anp anp bf num2str(out(15,1),'%3.2f') anp anp anp bf num2str(out(15,2),'%3.2f') anp anp bf num2str(out(15,3),'%3.2f') anp num2str(out(15,4),'%3.2f') dslash anp anp leftp num2str(out(16,1),'%3.2f') rightp anp anp anp leftp num2str(out(16,2),'%3.2f') rightp anp anp leftp num2str(out(16,3),'%3.2f') rightp anp dslash]
out_row9 = [ anp anp anp bf num2str(out(17,1),'%3.2f') anp anp anp bf num2str(out(17,2),'%3.2f') anp bf num2str(out(17,3),'%3.2f') anp num2str(out(17,4),'%3.2f') dslash anp anp anp leftp num2str(out(18,1),'%3.2f') rightp anp anp anp leftp num2str(out(18,2),'%3.2f') rightp anp leftp num2str(out(18,3),'%3.2f') rightp anp dslash]
out_row10 = ['$log\left(SP500_t - SP500_{t-4}\right)/4$' anp bf num2str(out(19,1),'%3.2f') anp anp anp bf num2str(out(19,2),'%3.2f') anp anp anp bf num2str(out(19,3),'%3.2f') anp num2str(out(19,4),'%3.2f') dslash anp leftp num2str(out(20,1),'%3.2f') rightp anp anp anp leftp num2str(out(20,2),'%3.2f') rightp anp anp anp leftp num2str(out(20,3),'%3.2f') rightp anp dslash]
out_row11 = ['$r_{t-8}^{Treas(10yr)} - r_{t-8}^{Treas(2yr)}$' anp bf num2str(out(21,1),'%3.2f') anp anp anp bf num2str(out(21,2),'%3.2f') anp anp anp bf num2str(out(21,3),'%3.2f') anp num2str(out(21,4),'%3.2f') dslash  anp leftp num2str(out(22,1),'%3.2f') rightp anp anp anp leftp num2str(out(22,2),'%3.2f') rightp anp anp anp leftp num2str(out(22,3),'%3.2f') rightp anp dslash]
out_row12 = ['$\widehat{\Delta GDP(E)_{t,t}}$, SPF forecast' anp bf num2str(out(23,1),'%3.2f') anp anp anp bf num2str(out(23,2),'%3.2f') anp anp anp bf num2str(out(23,3),'%3.2f') anp num2str(out(23,4),'%3.2f') dslash 'current quarter' anp leftp num2str(out(24,1),'%3.2f') rightp anp anp anp leftp num2str(out(24,2),'%3.2f') rightp anp anp anp leftp num2str(out(24,3),'%3.2f') rightp anp dslash]
out_row13 = ['$\widehat{\Delta GDP(E)_{t,t-1}}$, SPF forecast' anp bf num2str(out(25,1),'%3.2f') anp anp anp bf num2str(out(25,2),'%3.2f') anp anp anp bf num2str(out(25,3),'%3.2f') anp num2str(out(25,4),'%3.2f') dslash '1 quarter ahead' anp leftp num2str(out(26,1),'%3.2f') rightp anp anp anp leftp num2str(out(26,2),'%3.2f') rightp anp anp anp leftp num2str(out(26,3),'%3.2f') rightp anp dslash]
out_row14 = ['$\widehat{\Delta GDP(E)_{t,t-2}}$, SPF forecast' anp bf num2str(out(27,1),'%3.2f') anp anp anp bf num2str(out(27,2),'%3.2f') anp anp anp bf num2str(out(27,3),'%3.2f') anp num2str(out(27,4),'%3.2f') dslash '2 quarters ahead' anp leftp num2str(out(28,1),'%3.2f') rightp anp anp anp leftp num2str(out(28,2),'%3.2f') rightp anp anp anp leftp num2str(out(28,3),'%3.2f') rightp anp dslash]
out_table = [out_row1 out_row2 out_row3 hline out_row4 out_row5 out_row6 hline out_row7 out_row8 out_row9 hline out_row10 hline out_row11 hline out_row12 out_row13 out_row14]
dlmwrite('table_aux_forecasts_adv.txt',out_table,'');
