% Construct Solow residual using BEA and BLS data
%


load ./rawdata/NIPA/NIPAhist/table1_1_5Q.dat
ncon     = table1_1_5Q(7,5:88)'+table1_1_5Q(8,5:88)';
ngdp     = table1_1_5Q(3,5:88)';
load ./rawdata/NIPA/NIPA/table1_1_5Q.dat
ncon     = [ncon;table1_1_5Q(7,1:184)'+table1_1_5Q(8,1:184)'];
ngdp     = [ngdp;table1_1_5Q(3,1:184)'];
load ./rawdata/NIPA/NIPAhist/table1_1_6Q.dat
rgdp     = table1_1_6Q(3,5:88)';
rinv     = table1_1_6Q(9,5:88)'-table1_1_6Q(14,5:88)';
load ./rawdata/NIPA/NIPA/table1_1_6Q.dat
rgdp     = [rgdp;table1_1_6Q(3,1:184)'];
rinv     = [rinv;table1_1_6Q(9,1:184)'-table1_1_6Q(14,1:184)'];
load ./rawdata/NIPA/FA/table1_2A.dat
ka       = [table1_2A(2,24:89)';table1_2A(2,89)^2/table1_2A(2,88)];
kq       = ann2qtr([[1948:2014]',ka],4,0);
kq       = kq(:,2);  

load ./rawdata/Employment/CPU/pop1664.dat
load ./rawdata/Employment/CPU/avghrs.dat
p16q     = pop1664(3:270,2);
hpc      = avghrs(3:270,2);

tfp      = rgdp./(kq.^(.33).*(hpc.*p16q).^(.67));
tfp      = tfp/tfp(1);
lw       = (ncon./ngdp).*(hpc./(1-hpc));
t        = ptime(1948.1,268,4);

%
% Geometrically detrended
%
grate    = 1.0045.^[-267:0]';
gdppc    = rgdp./(grate.*p16q);
invpc    = rinv./(grate.*p16q);
cappc    = kq./(grate.*p16q);
tfppc    = gdppc./(cappc.^(.33).*hpc.^(.67));
zt       = tfppc/tfppc(1);
zt       = tfppc.^(1/.67);
lztlag0  = log(zt(1:267)/mean(zt));
lzt0     = log(zt(2:268)/mean(zt));
beta0    = ols([ones(267,1),lztlag0],lzt0);
eps0     = lzt0-beta0(2)*lztlag0;
lztlag   = log(zt(148:267)/mean(zt(149:268)));
lzt      = log(zt(149:268)/mean(zt(149:268)));
beta     = ols([ones(120,1),lztlag],lzt);
eps      = lzt-beta(2)*lztlag;

figure(1)
plot(t,(rgdp./p16q)/(rgdp(1)/p16q(1)),t,tfp/tfp(1));
title('Real GDP Per Capita and Measured TFP, 1948.1-2014.4')
legend('Real GDP','TFP','Location','Northwest')


%
% HP-filtered
%
hpgdp    = log(rgdp)-hptrend(log(rgdp),1600);
hphpc    = log(hpc)-hptrend(log(hpc),1600);
hphpc    = log(hpc)-hptrend(log(hpc),1600);
hplpr    = log(rgdp./hpc)-hptrend(log(rgdp./hpc),1600);
hptfp    = log(tfp)-hptrend(log(tfp),1600);
hplw     = log(lw)-hptrend(log(lw),1600);

figure(2)
plot(t,hpgdp,t,hptfp,t,hplpr)
title('GDP and Measured Productivity, All logged and detrended, 1948.1-2014.4')
legend('GDP','TFP','GDP per Hour')
figure(3)
plot(t,hphpc,t,hplw)
title('Hours and Measured Labor Wedge, All logged and detrended, 1948.1-2014.4')
legend('Hours','Labor Wedge')

