%slides_1b_210C
%Matlb code for creating figures 9 and 10
%Economics 210C, Spring 2019, Garey Ramey, UCSD

clear all
close all

%Import NIPA, NBER, CPS data
rawdat = xlsread('data_1_210C.xlsx');
GDP = rawdat(1:288,1);
Cons = rawdat(1:288,2);
Inv = rawdat(1:288,3);
Cons_Def = rawdat(1:288,10);
EmpPop = rawdat(1:288,16);  %First four quarters are missing

%Import Shiller data
rawdat2 = xlsread('Shiller_ie_data.xls');
m_price = rawdat2(913:1776,2);
N = size(m_price,1)/3;
q_price = zeros(N,1);
for i = 1:N
    q_price(i) = m_price(3*i);
end

%Compute quarterly growth rates
dGDP = log(GDP(2:288))-log(GDP(1:287));
dCons = log(Cons(2:288))-log(Cons(1:287));
dInv = log(Inv(2:288))-log(Inv(1:287));
dEmpPop = log(EmpPop(2:288))-log(EmpPop(1:287));
dSandP = log(q_price(2:288))-log(q_price(1:287));

%Initialize bootstrap simulation
n_sim = 1000;
t = size(dGDP,1);
conf = .95;

%Simulation for Cons
C_sim = zeros(9,n_sim);
for i = 1:n_sim
    for k = 0:4
        draw_ld = randi(t-k,t-k,1);
        draw_lg = randi(t-k,t-k,1);
        C_sim(5-k,i) = corr(dGDP(k+draw_ld),dCons(draw_ld));
        C_sim(5+k,i) = corr(dGDP(draw_lg),dCons(draw_lg+k));
    end
end
C_sim = sort(C_sim,2);
low = round(n_sim*(1-conf));   high = round(n_sim*conf);
plot_C_sim = zeros(9,3);
plot_C_sim(:,1) = C_sim(:,low);
plot_C_sim(:,2) = median(C_sim,2);
plot_C_sim(:,3) = C_sim(:,high);

%Simulation for Inv
I_sim = zeros(9,n_sim);
for i = 1:n_sim
    for k = 0:4
        draw_ld = randi(t-k,t-k,1);
        draw_lg = randi(t-k,t-k,1);
        I_sim(5-k,i) = corr(dGDP(k+draw_ld),dInv(draw_ld));
        I_sim(5+k,i) = corr(dGDP(draw_lg),dInv(draw_lg+k));
    end
end
I_sim = sort(I_sim,2);
low = round(n_sim*(1-conf));   high = round(n_sim*conf);
plot_I_sim = zeros(9,3);
plot_I_sim(:,1) = I_sim(:,low);
plot_I_sim(:,2) = median(I_sim,2);
plot_I_sim(:,3) = I_sim(:,high);

%Simulation for EmpPop
E_sim = zeros(9,n_sim);
for i = 1:n_sim
    for k = 0:4
        draw_ld = randi(t-k-4,t-k-4,1);
        draw_lg = randi(t-k-4,t-k-4,1);
        E_sim(5-k,i) = corr(dGDP(k+4+draw_ld),dEmpPop(4+draw_ld));
        E_sim(5+k,i) = corr(dGDP(4+draw_lg),dEmpPop(4+draw_lg+k));
    end
end
E_sim = sort(E_sim,2);
low = round(n_sim*(1-conf));   high = round(n_sim*conf);
plot_E_sim = zeros(9,3);
plot_E_sim(:,1) = E_sim(:,low);
plot_E_sim(:,2) = median(E_sim,2);
plot_E_sim(:,3) = E_sim(:,high);

%Simulation for SandP
S_sim = zeros(9,n_sim);
for i = 1:n_sim
    for k = 0:4
        draw_ld = randi(t-k,t-k,1);
        draw_lg = randi(t-k,t-k,1);
        S_sim(5-k,i) = corr(dGDP(k+draw_ld),dSandP(draw_ld));
        S_sim(5+k,i) = corr(dGDP(draw_lg),dSandP(draw_lg+k));
    end
end
S_sim = sort(S_sim,2);
low = round(n_sim*(1-conf));   high = round(n_sim*conf);
plot_S_sim = zeros(9,3);
plot_S_sim(:,1) = S_sim(:,low);
plot_S_sim(:,2) = median(S_sim,2);
plot_S_sim(:,3) = S_sim(:,high);

%Plot medians and 95% error bands for cross-correlations of dGDP with
%   dCons, dInv, dEmpPop, dSandP
figure(9)
subplot(2,2,1)
plot(1:9,plot_C_sim(:,1),'b:',1:9,plot_C_sim(:,2),'b', ...
    1:9,plot_C_sim(:,3),'b:',1:9,zeros(1,9),'k', ...
    1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k')
axis([1,9,-0.3,1])
title('Corr(dCons at t, dGDP at t-k)')
set(gca,'Xtick',1:1:9)
set(gca,'Xticklabel',{'-4','-3','-2','-1', ...
    '0','1','2','3','4'})

subplot(2,2,2)
plot(1:9,plot_I_sim(:,1),'r:',1:9,plot_I_sim(:,2),'r', ...
    1:9,plot_I_sim(:,3),'r:',1:9,zeros(1,9),'k',...
    1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k')
axis([1,9,-0.3,1])
title('Corr(dInv at t, dGDP at t-k)')
set(gca,'Xtick',1:1:9)
set(gca,'Xticklabel',{'-4','-3','-2','-1', ...
    '0','1','2','3','4'})

subplot(2,2,3)
plot(1:9,plot_E_sim(:,1),'g:',1:9,plot_E_sim(:,2),'g', ...
    1:9,plot_E_sim(:,3),'g:',1:9,zeros(1,9),'k', ...
    1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k')
axis([1,9,-0.3,1])
title('Corr(dEmpPop at t, dGDP at t-k)')
set(gca,'Xtick',1:1:9)
set(gca,'Xticklabel',{'-4','-3','-2','-1', ...
    '0','1','2','3','4'})

subplot(2,2,4)
hold on
plot(1:9,plot_S_sim(:,1),'r:',1:9,plot_S_sim(:,2),'r', ...
    1:9,plot_S_sim(:,3),'r:','Color',[.9,.5,.1])
plot(1:9,zeros(1,9),'k',...
    1:9,[-1000*ones(1,4),0,1e20*ones(1,4)],'k')
axis([1,9,-0.3,1])
title('Corr(dS&P at t, dGDP at t-k)')
set(gca,'Xtick',1:1:9)
set(gca,'Xticklabel',{'-4','-3','-2','-1', ...
    '0','1','2','3','4'})
hold off

%Simulation for GDP autocorrelations
G_sim = zeros(8,n_sim);
for i = 1:n_sim
    for k = 1:8
        draw = randi(t-k,t-k,1);
        G_sim(k,i) = corr(dGDP(k+draw),dGDP(draw));
    end
end
G_sim = sort(G_sim,2);
low = round(n_sim*(1-conf));   high = round(n_sim*conf);
plot_G_sim = zeros(8,3);
plot_G_sim(:,1) = G_sim(:,low);
plot_G_sim(:,2) = median(G_sim,2);
plot_G_sim(:,3) = G_sim(:,high);

%Plot median autocorrelations and 95% error bands of dGDP
figure(10)
plot(1:8,plot_G_sim(:,1),'k:',1:8,plot_G_sim(:,2),'k', ...
    1:8,plot_G_sim(:,3),'k:',1:8,zeros(1,8),'k')
%axis([1,9,-0.3,1])
title('Autocorrelations of GDP Growth Rates')
set(gca,'Xtick',1:1:8)
set(gca,'Xticklabel',{'1','2','3','4', ...
    '5','6','7','8'})

