%slides_3a_210C
%Matlb code for creating figures 14-16
%Uses function PLOT_BANDS.M, CC_BOOT.M
%Economics 210C, Spring 2019, Garey Ramey, UCSD

clear
close all

%Import NIPA, Fernald, NBER data
rawdat = xlsread('data_3a_210C.xlsx');
GDP = log(rawdat(:,1));
TFP = rawdat(:,2);
TFP_util = rawdat(:,3);
Gov = log(rawdat(:,4));
Rec = rawdat(:,5);
T = size(GDP,1);

%Plot annual growth rates
figure(14)
dGDP_A = GDP(5:T)-GDP(1:T-4);
dTFP_A = TFP(5:T)-TFP(1:T-4);
dTFP_util_A = TFP_util(5:T)-TFP_util(1:T-4);
dGov_A = Gov(5:T)-Gov(1:T-4);
plot_bands(Rec(5:T),[-.35,.55],.3,'g')
hold on
h(1) = plot(1:T-4,dGDP_A,'k');
h(2) = plot(1:T-4,dTFP_A,'r');
%h(3) = plot(1:T-4,dTFP_util_A,'r--');
h(4) = plot(1:T-4,dGov_A,'Color',[.9,.5,.1]);
h(5) = plot(1:T-4,zeros(1,T-4),'k');
axis([1,T-4,-.05,.12])
set(gca,'Xtick',13:40:T-4)
set(gca,'Xticklabel',{'1952:1','1962:1','1972:1','1982:1', ...
    '1992:1','2002:1','2012:1'})
title('Annual Growth Rates of GDP, TFP and Gov Spending')
ylabel('Log Difference')
legend(h([1 2 4]),'GDP','TFP','Gov','Location','Northeast')
hold off

%HP filtered cyclical components
GDP_HP = hp(GDP,1600);
TFP_HP = hp(TFP,1600);
TFP_util_HP = hp(TFP_util,1600);
Gov_HP = hp(Gov,1600);

%Band pass filtered cyclical components
K = 12;
GDP_BP = BPF(GDP,6,32,K);
TFP_BP = BPF(TFP,6,32,K);
TFP_util_BP = BPF(TFP_util,6,32,K);
Gov_BP = BPF(Gov,6,32,K);

%Plot cyclical components of quarterly GDP, TFP and Gov for HP filtered 
%   and band pass filtered data
figure(15)
subplot(2,1,1)
plot_bands(Rec(5:T),[-.35,.55],.3,'g')
hold on
h(1) = plot(1:T,GDP_HP,'k');
h(2) = plot(1:T,TFP_HP,'r');
%h(3) = plot(1:T,TFP_util_HP,'r--');
h(4) = plot(1:T,Gov_HP,'Color',[.9,.5,.1]);
h(5) = plot(1:T,zeros(1,T),'k');
axis([1,T,-.055,.075])
set(gca,'Xtick',17:40:T-4)
set(gca,'Xticklabel',{'1952:1','1962:1','1972:1','1982:1', ...
    '1992:1','2002:1','2012:1'})
title('Quarterly Cyclical Components, HP Filtered Log Levels')
legend(h([1 2 4]),'GDP','TFP','Gov','Location','Northeast')
hold off

subplot(2,1,2)
plot_bands(Rec(5:T),[-.35,.55],.3,'g')
hold on
h(1) = plot(1:T,GDP_BP,'k');
h(2) = plot(1:T,TFP_BP,'r');
%h(3) = plot(1:T,TFP_util_BP,'r--');
h(4) = plot(1:T,Gov_BP,'Color',[.9,.5,.1]);
h(5) = plot(1:T,zeros(1,T),'k');
axis([1,T,-.055,.075])
set(gca,'Xtick',17:40:T-4)
set(gca,'Xticklabel',{'1952:1','1962:1','1972:1','1982:1', ...
    '1992:1','2002:1','2012:1'})
title('Quarterly Cyclical Components, Band-Pass Filtered Log Levels')
legend(h([1 2 4]),'GDP','TFP','Gov','Location','Northeast')
hold off

%Estimate shocks to polynomial detrended TFP and Gov
tr = 'lin';     %Choose 'lin', 'quad' or 'cub'
switch tr
    case 'lin'
        X_tr = [ones(T-1,1) (1:T-1)'];
    case 'quad'
        X_tr = [ones(T-1,1) (1:T-1)' ((1:T-1).^2)'];
    case 'cub'
        X_tr = [ones(T-1,1) (1:T-1)' ((1:T-1).^2)' 1e-4*((1:T-1).^3)'];
end
X = [TFP(1:T-1) Gov(1:T-1) X_tr];
Y = [TFP(2:T) Gov(2:T)];
coef = (X'*X)\(X'*Y);
resid = Y-X*coef;
st_dev = std(resid);
para_Z = [coef(1,1) st_dev(1)]
para_G = [coef(2,2) st_dev(2)]
Rsq = ones(2)-(resid'*resid)./((Y-mean(Y))'*(Y-mean(Y)));
Rsq_shocks = diag(Rsq)

%Estimate polynomial detrended GDP
coef_GDP = (X_tr'*X_tr)\(X_tr'*GDP(2:T));
GDP_cyc = GDP(2:T)-X_tr*coef_GDP;

%Calculate cross-correlations and bootstrap error bands
n = 1000;   conf = .95;
cc_eps = cc_boot(resid(:,1),GDP_cyc,4,n,conf);
cc_nu = cc_boot(resid(:,2),GDP_cyc,4,n,conf);

%Plot medians and 95% error bands for cross-correlations of cyclical 
%   component of GDP with estimated shocks to TFP and Gov
figure(16)
subplot(2,1,1)
hold on
h(1) = plot(1:9,cc_eps(:,1),'r:');
h(2) = plot(1:9,cc_eps(:,2),'r');
h(3) = plot(1:9,cc_eps(:,3),'r:');
h(4) = plot(1:9,zeros(1,9),'k');
h(5) = plot(1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k');
axis([1,9,-0.3,1])
title('Corr(TFP shock at t, GDP at t-k), Log-linear Detrended')
%legend(h([2 5]),'Hours','Productivity', ...
%    'Location','Northwest')
set(gca,'Xtick',1:1:9)
set(gca,'Xticklabel',{'-4','-3','-2','-1', ...
    '0','1','2','3','4'})
hold off

subplot(2,1,2)
hold on
h(1) = plot(1:9,cc_nu(:,1),'b:');
h(2) = plot(1:9,cc_nu(:,2),'b');
h(3) = plot(1:9,cc_nu(:,3),'b:');
h(4) = plot(1:9,zeros(1,9),'k');
h(5) = plot(1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k');
axis([1,9,-0.3,1])
title('Corr(Gov shock at t, GDP at t-k), Log-linear Detrended')
%legend(h([2 5]),'Hours','Productivity', ...
%    'Location','Northwest')
set(gca,'Xtick',1:1:9)
set(gca,'Xticklabel',{'-4','-3','-2','-1', ...
    '0','1','2','3','4'})
hold off

%Parallel calculations for utilization-adjusted TFP series
X = [TFP_util(1:T-1) Gov(1:T-1) X_tr];
Y = [TFP_util(2:T) Gov(2:T)];
coef = (X'*X)\(X'*Y);
resid = Y-X*coef;
st_dev = std(resid);
para_Z_u = [coef(1,1) st_dev(1)];
para_G_u = [coef(2,2) st_dev(2)];
Rsq = ones(2)-(resid'*resid)./((Y-mean(Y))'*(Y-mean(Y)));
Rsq_u_shocks = diag(Rsq);

cc_eps_u = cc_boot(resid(:,1),GDP_cyc,4,n,conf);
cc_nu_u = cc_boot(resid(:,2),GDP_cyc,4,n,conf);
