%RBC_evaluate
%Matlb code for creating figures 14-18
%Calculates empirical moments and compares them to simulated moments
%Uses functions HP.M, SD_BOOT.M, CC_BOOT.M, AC_BOOT.M and C_COR.M, and 
%   script RBC_SOLVE_SIM.M
%Economics 210C, Spring 2019, Garey Ramey, UCSD

clear
close all

%IMPORT DATA, CALCULATE EMPIRICAL MOMENTS
%Bootstrapping parameters for data
n = 1000;   conf = .95;

%Import NIPA, NBER, BLS data, log and HP filter it
rawdat = xlsread('data_3_210C.xlsx');
sample = 1:size(rawdat,1);
%Cooley-Prescott sample 1954:I-1991_IV
%sample = 25:174;

GDP = log(rawdat(sample,1));
Cons = log(rawdat(sample,2));
Inv = log(rawdat(sample,3));
Gov = log(rawdat(sample,4));
Hours = log(rawdat(sample,5));
Prod = GDP-Hours;

GDP_HP = hp(GDP,1600);
Cons_HP = hp(Cons,1600);
Inv_HP = hp(Inv,1600);
Gov_HP = hp(Gov,1600);
Hours_HP = hp(Hours,1600);
Prod_HP = hp(Prod,1600);

%Calculate moments
sd_GDP = sd_boot(GDP_HP,n,conf);
sd_Cons = sd_boot(Cons_HP,n,conf);
sd_Inv = sd_boot(Inv_HP,n,conf);
sd_Gov = sd_boot(Gov_HP,n,conf);
sd_Hours = sd_boot(Hours_HP,n,conf);
sd_Prod = sd_boot(Prod_HP,n,conf);
sd_all = [sd_GDP; sd_Cons; sd_Inv; sd_Hours; sd_Prod; sd_Gov]';

cc_Cons = cc_boot(Cons_HP,GDP_HP,5,n,conf);
cc_Inv = cc_boot(Inv_HP,GDP_HP,5,n,conf);
cc_Gov = cc_boot(Gov_HP,GDP_HP,5,n,conf);
cc_Hours = cc_boot(Hours_HP,GDP_HP,5,n,conf);
cc_Prod = cc_boot(Prod_HP,GDP_HP,5,n,conf);

ac_dGDP = ac_boot(GDP(2:end)-GDP(1:end-1),5,n,conf);

sd_all_sim = zeros(6,6);
cc_c_sim_all = zeros(6,11);
cc_i_sim_all = zeros(6,11);
cc_l_sim_all = zeros(6,11);
cc_p_sim_all = zeros(6,11);
ac_dy_sim_all = zeros(6,5);

eta_vec = [.25 1 10];
para = 'nogov';
for i = 1:3
    eta = eta_vec(i);
    RBC_solve_sim
    y_sim_HP = hp(y_sim,1600);
    c_sim_HP = hp(c_sim,1600);
    i_sim_HP = hp(i_sim,1600);
    g_sim_HP = hp(g_sim,1600);
    l_sim_HP = hp(l_sim,1600);
    p_sim_HP = hp(p_sim,1600);
    sd_y_sim = std(y_sim_HP);
    sd_c_sim = std(c_sim_HP);
    sd_i_sim = std(i_sim_HP);
    sd_g_sim = std(g_sim_HP);
    sd_l_sim = std(l_sim_HP);
    sd_p_sim = std(p_sim_HP);
    if i == 1                   %Parameterization used for histograms
        sd_y_sim_0 = sd_y_sim;
        sd_c_sim_0 = sd_c_sim;
        sd_i_sim_0 = sd_i_sim;
        sd_l_sim_0 = sd_l_sim;
        sd_p_sim_0 = sd_p_sim;
    end
    Xsd = mean([sd_y_sim; sd_c_sim; sd_i_sim; sd_l_sim; sd_p_sim; ...
        sd_g_sim],2)';
    sd_all_sim(i,:) = Xsd; 
    cc_c_sim_all(i,:) = mean(c_cor(c_sim_HP,y_sim_HP,5),2)';
    cc_i_sim_all(i,:) = mean(c_cor(i_sim_HP,y_sim_HP,5),2)';
    cc_l_sim_all(i,:) = mean(c_cor(l_sim_HP,y_sim_HP,5),2)';
    cc_p_sim_all(i,:) = mean(c_cor(p_sim_HP,y_sim_HP,5),2)';
    dy_sim = y_sim(2:end,:)-y_sim(1:end-1,:);
    Xac = mean(c_cor(dy_sim,dy_sim,5),2)';
    ac_dy_sim_all(i,:) = Xac(1,7:11);
end
para = 'gov';
for i = 4:6
    eta = eta_vec(i-3);
    RBC_solve_sim
    y_sim_HP = hp(y_sim,1600);
    c_sim_HP = hp(c_sim,1600);
    i_sim_HP = hp(i_sim,1600);
    g_sim_HP = hp(g_sim,1600);
    l_sim_HP = hp(l_sim,1600);
    p_sim_HP = hp(p_sim,1600);
    sd_y_sim = std(y_sim_HP);
    sd_c_sim = std(c_sim_HP);
    sd_i_sim = std(i_sim_HP);
    sd_g_sim = std(g_sim_HP);
    sd_l_sim = std(l_sim_HP);
    sd_p_sim = std(p_sim_HP);
    Xsd = mean([sd_y_sim; sd_c_sim; sd_i_sim; sd_l_sim; sd_p_sim; ...
        sd_g_sim],2)';
    sd_all_sim(i,:) = Xsd; 
    cc_c_sim_all(i,:) = mean(c_cor(c_sim_HP,y_sim_HP,5),2)';
    cc_i_sim_all(i,:) = mean(c_cor(i_sim_HP,y_sim_HP,5),2)';
    cc_l_sim_all(i,:) = mean(c_cor(l_sim_HP,y_sim_HP,5),2)';
    cc_p_sim_all(i,:) = mean(c_cor(p_sim_HP,y_sim_HP,5),2)';
    dy_sim = y_sim(2:end,:)-y_sim(1:end-1,:);
    Xac = mean(c_cor(dy_sim,dy_sim,5),2)';
    ac_dy_sim_all(i,:) = Xac(1,7:11);
end

sd_table = [sd_all(2,:)' sd_all_sim']

%PLOT EMPIRICAL AND SIMULATED MOMENTS
figure(14)
subplot(2,2,1)
hold on
h(1) = plot(1:11,cc_c_sim_all(1,:),'b');
h(2) = plot(1:11,cc_c_sim_all(2,:),'b--');
h(3) = plot(1:11,cc_c_sim_all(3,:),'b:');
h(4) = plot(1:11,cc_Cons(:,2)','k');
h(5) = plot(1:11,zeros(1,11),'k');
h(6) = plot(1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k');
axis([1,11,-0.3,1])
title({'Corr(Cons at t, GDP at t-k)';'No Government Sector'})
legend(h([1 2 3]),{'eta = .25','eta = 1','eta = 10'},'Location','Southeast')
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})
hold off

subplot(2,2,2)
plot(1:11,cc_Cons(:,2)','k',1:11,cc_c_sim_all(4,:),'b', ...
    1:11,cc_c_sim_all(5,:),'b--',1:11,cc_c_sim_all(6,:),'b:', ...
    1:11,zeros(1,11),'k',1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k')
axis([1,11,-0.3,1])
title({'Corr(Cons at t, GDP at t-k)';'With Government Sector'})
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})

subplot(2,2,3)
plot(1:11,cc_Inv(:,2)','k',1:11,cc_i_sim_all(1,:),'r', ...
    1:11,cc_i_sim_all(2,:),'r--',1:11,cc_i_sim_all(3,:),'r:', ...
    1:11,zeros(1,11),'k',1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k')
axis([1,11,-0.3,1])
title({'Corr(Inv at t, GDP at t-k)';'No Government Sector'})
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})

subplot(2,2,4)
plot(1:11,cc_Inv(:,2)','k',1:11,cc_i_sim_all(4,:),'r', ...
    1:11,cc_i_sim_all(5,:),'r--',1:11,cc_i_sim_all(6,:),'r:', ...
    1:11,zeros(1,11),'k',1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k')
axis([1,11,-0.3,1])
title({'Corr(Inv at t, GDP at t-k)';'With Government Sector'})
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})


figure(15)
subplot(2,2,1)
plot(1:11,cc_Hours(:,2)','k',1:11,cc_l_sim_all(1,:),'b', ...
    1:11,cc_l_sim_all(2,:),'b--',1:11,cc_l_sim_all(3,:),'b:', ...
    1:11,zeros(1,11),'k',1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k')
axis([1,11,-0.3,1])
title({'Corr(Hours at t, GDP at t-k)';'No Government Sector'})
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})

subplot(2,2,2)
plot(1:11,cc_Hours(:,2)','k',1:11,cc_l_sim_all(4,:),'b', ...
    1:11,cc_l_sim_all(5,:),'b--',1:11,cc_l_sim_all(6,:),'b:', ...
    1:11,zeros(1,11),'k',1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k')
axis([1,11,-0.3,1])
title({'Corr(Hours at t, GDP at t-k)';'With Government Sector'})
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})

subplot(2,2,3)
plot(1:11,cc_Prod(:,2)','k',1:11,cc_p_sim_all(1,:),'r', ...
    1:11,cc_p_sim_all(2,:),'r--',1:11,cc_p_sim_all(3,:),'r:', ...
    1:11,zeros(1,11),'k',1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k')
axis([1,11,-0.3,1])
title({'Corr(Prod at t, GDP at t-k)';'No Government Sector'})
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})

subplot(2,2,4)
plot(1:11,cc_Prod(:,2)','k',1:11,cc_p_sim_all(4,:),'r', ...
    1:11,cc_p_sim_all(5,:),'r--',1:11,cc_p_sim_all(6,:),'r:', ...
    1:11,zeros(1,11),'k',1:11,[-1000*ones(1,5),0,1e20*ones(1,5)],'k')
axis([1,11,-0.3,1])
title({'Corr(Prod at t, GDP at t-k)';'With Government Sector'})
set(gca,'Xtick',1:1:11)
set(gca,'Xticklabel',{'-5','-4','-3','-2','-1', ...
    '0','1','2','3','4','5'})


figure(16)
subplot(2,1,1)
plot(1:5,ac_dGDP(:,2)','k',1:5,ac_dy_sim_all(1,:),'b', ...
    1:5,ac_dy_sim_all(2,:),'b--',1:5,ac_dy_sim_all(3,:),'b:', ...
    1:5,zeros(1,5),'k')
title({'Autocorrelations of GDP Growth Rates';'No Government Sector'})
set(gca,'Xtick',1:1:5)
set(gca,'Xticklabel',{'1','2','3','4','5'})

subplot(2,1,2)
plot(1:5,ac_dGDP(:,2)','k',1:5,ac_dy_sim_all(4,:),'b', ...
    1:5,ac_dy_sim_all(5,:),'b--',1:5,ac_dy_sim_all(6,:),'b:', ...
    1:5,zeros(1,5),'k')
title({'Autocorrelations of GDP Growth Rates';'With Government Sector'})
set(gca,'Xtick',1:1:5)
set(gca,'Xticklabel',{'1','2','3','4','5'})

%PLOT HISTOGRAMS ACROSS SIMULATION RELICATIONS
%Uses parameterization eta = .25, para = 'nogov'
figure(17)
histogram(sd_y_sim_0,100)
title({'Histogram of St Dev of GDP Across 1000 Simulations'; ...
    '\eta  = .25, No Government Sector'})


figure(18)
subplot(2,2,1)
histogram(sd_c_sim_0./sd_y_sim_0,100)
title('SD of Cons/SD of GDP')

subplot(2,2,2)
histogram(sd_i_sim_0./sd_y_sim_0,100)
title('SD of Inv/SD of GDP')

subplot(2,2,3)
histogram(sd_l_sim_0./sd_y_sim_0,100)
title('SD of Hours/SD of GDP')

subplot(2,2,4)
histogram(sd_p_sim_0./sd_y_sim_0,100)
title('SD of Prod/SD of GDP')
