%slides_2_210C_v2
%Matlb code for creating figures 11-13
%Economics 210C, Spring 2019, Garey Ramey, UCSD

clear all
close all

%Import NIPA, NBER, CPS data
rawdat = xlsread('data_2_210C.xlsx');
dGDP_A = log(rawdat(5:284,1))-log(rawdat(1:280,1));
dHours_A = log(rawdat(5:284,2))-log(rawdat(1:280,2));
Rec = rawdat(1:284,3);
dProd_A = dGDP_A-dHours_A;

%Plot GDP, Hours, Productivity
figure(11)
plot_bands(Rec(5:280),[-.06,.14],.3,'g')
hold on
h(1) = plot(1:280,dGDP_A,'k');
h(2) = plot(1:280,dHours_A,'b');
h(3) = plot(1:280,dProd_A,'r');
h(4) = plot(1:280,zeros(1,280),'k');
axis([1,284,-.06,.14])
set(gca,'Xtick',13:40:280)
set(gca,'Xticklabel',{'1952:1','1962:1','1972:1','1982:1', ...
    '1992:1','2002:1','2012:1'})
title('Annual Growth Rates of GDP, Hours and Productivity')
ylabel('Log Difference')
legend(h([1 2 3]),'GDP','Hours','Productivity', ...
    'Location','Northeast')
hold off

%Quarterly growth rates
dGDP = log(rawdat(2:284,1))-log(rawdat(1:283,1));
dHours = log(rawdat(2:284,2))-log(rawdat(1:283,2));
dProd = dGDP-dHours;

%HP filtered cyclical components
GDP_HP = hp(log(rawdat(:,1))',1600)';
Hours_HP = hp(log(rawdat(:,2))',1600)';
Prod_HP = hp(log(rawdat(:,1))'-log(rawdat(:,2))',1600)';

%Band pass filtered cyclical components
K = 12;
GDP_BP = BPF(log(rawdat(:,1)),6,32,K);
Hours_BP = BPF(log(rawdat(:,2)),6,32,K);
Prod_BP = BPF(log(rawdat(:,1))-log(rawdat(:,2)),6,32,K);

%Plot cyclical components of quarterly GDP, Hours, Productivity for
%   HP filtered and band pass filtered data
figure(12)
subplot(2,1,1)
plot_bands(Rec(1:284),[-.06,.09],.3,'g')
hold on
h(1) = plot(1:284,GDP_HP,'k');
h(2) = plot(1:284,Hours_HP,'b');
h(3) = plot(1:284,Prod_HP,'r');
h(4) = plot(1:284,zeros(1,284),'k');
axis([1,284,-.06,.09])
set(gca,'Xtick',17:40:284)
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 3]),'GDP','Hours','Productivity', ...
    'Location','Northeast')
hold off

subplot(2,1,2)
plot_bands(Rec(1:284),[-.06,.09],.3,'g')
hold on
h(1) = plot(1:284,GDP_BP,'k');
h(2) = plot(1:284,Hours_BP,'b');
h(3) = plot(1:284,Prod_BP,'r');
h(4) = plot(1:284,zeros(1,284),'k');
axis([1,284,-.06,.09])
set(gca,'Xtick',17:40:284)
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 3]),'GDP','Hours','Productivity', ...
    'Location','Northeast')
hold off

%Calculate cross-correlations and bootstrap error bands
n = 1000;   conf = .95;
cc_Hours_HP = cc_boot(Hours_HP,GDP_HP,4,n,conf);
cc_Prod_HP = cc_boot(Prod_HP,GDP_HP,4,n,conf);
cc_Hours_BP = cc_boot(Hours_BP,GDP_BP,4,n,conf);
cc_Prod_BP = cc_boot(Prod_BP,GDP_BP,4,n,conf);

%Plot medians and 95% error bands for cross-correlations of GDP with
%   Hours and Prod for HP filtered and band pass filtered data
figure(13)
subplot(2,1,1)
hold on
h(1) = plot(1:9,cc_Hours_HP(:,1),'b:');
h(2) = plot(1:9,cc_Hours_HP(:,2),'b');
h(3) = plot(1:9,cc_Hours_HP(:,3),'b:');
h(4) = plot(1:9,cc_Prod_HP(:,1),'r:');
h(5) = plot(1:9,cc_Prod_HP(:,2),'r');
h(6) = plot(1:9,cc_Prod_HP(:,3),'r:');
h(7) = plot(1:9,zeros(1,9),'k');
h(8) = plot(1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k');
axis([1,9,-0.3,1.2])
title('Corr(Variable at t, GDP at t-k), HP Filtered')
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_Hours_BP(:,1),'b:');
h(2) = plot(1:9,cc_Hours_BP(:,2),'b');
h(3) = plot(1:9,cc_Hours_BP(:,3),'b:');
h(4) = plot(1:9,cc_Prod_BP(:,1),'r:');
h(5) = plot(1:9,cc_Prod_BP(:,2),'r');
h(6) = plot(1:9,cc_Prod_BP(:,3),'r:');
h(7) = plot(1:9,zeros(1,9),'k');
h(8) = plot(1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k');
axis([1,9,-0.3,1.2])
title('Corr(Variable at t, GDP at t-k), HP Filtered')
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

%Standard deviations of cyclical components
SD_Diff = std([dGDP dHours dProd]);
SD_HP = std([GDP_HP Hours_HP Prod_HP]);
SD_BP = std([GDP_BP Hours_BP Prod_BP]);
SDs = [SD_Diff; SD_HP; SD_BP]

%Correlations of cyclical components of Hours and Prod with GDP
corr_Diff = [corr(dGDP,dHours) corr(dGDP,dProd)];
corr_HP = [corr(GDP_HP,Hours_HP) corr(GDP_HP,Prod_HP)];
corr_BP = [corr(GDP_BP,Hours_BP) corr(GDP_BP,Prod_BP)];
corrs = [corr_Diff(1) corr_Diff(2); corr_HP(1) corr_HP(2); ...
    corr_BP(1) corr_BP(2)]

%Betas of cyclical components of Hours and Prod with respect to GDP
cov_Diff = [cov(dGDP,dHours) cov(dGDP,dProd)];
cov_HP = [cov(GDP_HP,Hours_HP) cov(GDP_HP,Prod_HP)];
cov_BP = [cov(GDP_BP,Hours_BP) cov(GDP_BP,Prod_BP)];
beta_Diff = [cov_Diff(1,2) cov_Diff(1,4)]./var(dGDP);
beta_HP = [cov_HP(1,2) cov_HP(1,4)]./var(GDP_HP);
beta_BP = [cov_BP(1,2) cov_BP(1,4)]./var(GDP_BP);
betas = [beta_Diff; beta_HP; beta_BP]
