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

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

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);


out = zeros(22,6);


tstart = 3;
tend = 127;


clear y1 n control1 x

'Latest RGDP on initial RGDP, RGDI 1978-present'
y1 = RGDP_latest(tstart:tend,:);
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' 0 0 rsq];
out(2,:) = [se_newey' 0 0 0];

'Latest RGDP on initial RGDP(t,t-1), RGDI(t,t-1) 1978-present'
y1 = RGDP_latest(tstart:tend,:);
n = size(y1,1);
control1 = ones(size(y1));
x = [RGDP_adv(tstart:tend,:) RGDP_adv(tstart-1:tend-1,:) RGDI_adv(tstart:tend,:) 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];


'Latest RGDP on initial RGDP(t-1), RGDI(t-1) 1978-present'
y1 = RGDP_latest(tstart:tend,:);
n = size(y1,1);
control1 = ones(size(y1));
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(5,:) = [beta' 0 0 rsq];
out(6,:) = [se_newey' 0 0 0];



'Initial RGDP on initial RGDP(t-1), RGDI(t-1) 1978-present'
y1 = RGDP_fin(tstart:tend,:);
n = size(y1,1);
control1 = ones(size(y1));
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(17,:) = [beta' 0 0 rsq];
out(18,:) = [se_newey' 0 0 0];

'Initial RGDP on initial RGDP(t-2), RGDI(t-2) 1978-present'
y1 = RGDP_fin(tstart:tend,:);
n = size(y1,1);
control1 = ones(size(y1));
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(23,:) = [beta' 0 0 rsq];
out(24,:) = [se_newey' 0 0 0];

out_row1 = ['$\left(GDP(E)_t/GDP(E)_{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,6),'%3.2f') dslash 'latest' 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 bf num2str(out(3,1),'%3.2f') anp bf num2str(out(3,2),'%3.2f') anp anp bf num2str(out(3,3),'%3.2f') anp bf num2str(out(3,4),'%3.2f') anp anp bf num2str(out(3,5),'%3.2f') anp num2str(out(3,6),'%3.2f') dslash anp leftp num2str(out(4,1),'%3.2f') rightp anp leftp num2str(out(4,2),'%3.2f') rightp anp anp leftp num2str(out(4,3),'%3.2f') rightp anp leftp num2str(out(4,4),'%3.2f') rightp anp anp leftp num2str(out(4,5),'%3.2f') rightp anp dslash]
out_row3 = [ anp anp bf num2str(out(5,1),'%3.2f') anp anp anp bf num2str(out(5,2),'%3.2f') anp anp bf num2str(out(5,3),'%3.2f') anp num2str(out(5,6),'%3.2f') dslash anp anp leftp num2str(out(6,1),'%3.2f') rightp anp anp anp leftp num2str(out(6,2),'%3.2f') rightp anp anp leftp num2str(out(6,3),'%3.2f') rightp anp dslash]
out_row4 = ['$\left(GDP(E)_t/GDP(E)_{t-1}\right)^4,$'  anp anp bf num2str(out(17,1),'%3.2f') anp anp anp bf num2str(out(17,2),'%3.2f') anp anp bf num2str(out(17,3),'%3.2f') anp num2str(out(17,6),'%3.2f') dslash 'initial (3rd)'  anp anp leftp num2str(out(18,1),'%3.2f') rightp anp anp anp leftp num2str(out(18,2),'%3.2f') rightp anp anp leftp num2str(out(18,3),'%3.2f') rightp anp dslash]
out_row5 = [ anp anp anp bf num2str(out(23,1),'%3.2f') anp anp anp bf num2str(out(23,2),'%3.2f') anp bf num2str(out(23,3),'%3.2f') anp num2str(out(23,6),'%3.2f') dslash anp anp anp leftp num2str(out(24,1),'%3.2f') rightp anp anp anp leftp num2str(out(24,2),'%3.2f') rightp anp leftp num2str(out(24,3),'%3.2f') rightp anp dslash]
out_table = [out_row1 out_row2 out_row3 hline out_row4 out_row5 hline]
dlmwrite('table_GDP_forecasts_2adv.txt',out_table,'');
