%slides_3bc_210C_simpl
%Simplified Matlab code for creating RBC tables and figures
%Uses script RBC_SOLVE_SIM_SIMPL.M and functions RE_LIN_SOLVE.M, HP.M, 
%   and C_COR.M
%Economics 210C, Spring 2019, Garey Ramey, UCSD

clear
close all

RBC_solve_sim_simpl

%CALCULATE MOMENTS OF SIMULATED DATA
%HP filter simulated variables, after discarding initial observations
y_sim_HP = hp(y_sim(T_init+1:end),1600);
c_sim_HP = hp(c_sim(T_init+1:end),1600);
i_sim_HP = hp(i_sim(T_init+1:end),1600);
g_sim_HP = hp(g_sim(T_init+1:end),1600);
l_sim_HP = hp(l_sim(T_init+1:end),1600);
p_sim_HP = hp(p_sim(T_init+1:end),1600);

%Calculate standard deviations
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);
sd_vec = [sd_y_sim; sd_c_sim; sd_i_sim; sd_l_sim; sd_p_sim; sd_g_sim]

%Calculate hours/productivity statistics
hrs_prd = [mean(diag(corr(p_sim_HP,l_sim_HP)))
        mean(std(l_sim_HP)./std(p_sim_HP))]
    
%Calculate cross-correlations
cc_c_sim = mean(c_cor(c_sim_HP,y_sim_HP,5),2)';
cc_i_sim = mean(c_cor(i_sim_HP,y_sim_HP,5),2)';
cc_l_sim = mean(c_cor(l_sim_HP,y_sim_HP,5),2)';
cc_p_sim = mean(c_cor(p_sim_HP,y_sim_HP,5),2)';
cc_vec = [cc_c_sim; cc_i_sim; cc_l_sim; cc_p_sim]

%Calculate autocorrelations of output growth rates
dy_sim = y_sim(2:end,:)-y_sim(1:end-1,:);
ac = mean(c_cor(dy_sim,dy_sim,5),2)'

%GENERATE AND PLOT IMPULSE RESPONSES
%Initialize IRFs for z and g shocks
T_ir = 26;
IR_S_z = zeros(T_ir,3);
IR_U_z = zeros(T_ir,1);
IR_v_z = zeros(T_ir,2);
IR_S_g = zeros(T_ir,3);
IR_U_g = zeros(T_ir,1);
IR_v_g = zeros(T_ir,2);

%IRF for z shock - shock occurs at t = 2
z_shock = .01;
IR_S_z(2,2) = z_shock;
IR_U_z(2,:) = PHI_U*IR_S_z(2,:)';
IR_v_z(2,:) = PHI_v*[IR_U_z(2,:); IR_S_z(2,:)']; 
for t = 3:T_ir
    IR_S_z(t,:) = PHI_S*IR_S_z(t-1,:)';
    IR_U_z(t,:) = PHI_U*IR_S_z(t,:)';
    IR_v_z(t,:) = PHI_v*[IR_U_z(t,:); IR_S_z(t,:)'];
end
IR_y_z = IR_S_z(:,2)+alpha*IR_S_z(:,1)+(1-alpha)*IR_v_z(:,2);
IR_i_z = (Y*IR_y_z-C*IR_v_z(:,1)-G*IR_S_z(:,3))/(delta*K);
IR_p_z = IR_y_z-IR_v_z(:,2);

%IRF for g shock - shock occurs at t = 2
g_shock = .01;
IR_S_g(2,3) = g_shock; 
IR_U_g(2,:) = PHI_U*IR_S_g(2,:)';
IR_v_g(2,:) = PHI_v*[IR_U_g(2,:); IR_S_g(2,:)']; 
for t = 3:T_ir
    IR_S_g(t,:) = PHI_S*IR_S_g(t-1,:)';
    IR_U_g(t,:) = PHI_U*IR_S_g(t,:)';
    IR_v_g(t,:) = PHI_v*[IR_U_g(t,:); IR_S_g(t,:)']; 
end
IR_y_g = IR_S_g(:,2)+alpha*IR_S_g(:,1)+(1-alpha)*IR_v_g(:,2);
IR_i_g = (Y*IR_y_g-C*IR_v_g(:,2)-G*IR_S_g(:,3))/(delta*K);
IR_p_g = IR_y_g-IR_v_g(:,2);

%Plot impulse-responses for TFP shock
figure(1)
subplot(2,1,1)
hold on
plot(1:T_ir,IR_y_z,'k',1:T_ir,IR_v_z(:,1),'b',1:T_ir,IR_i_z,'r')
plot(1:T_ir,zeros(T_ir,1),'k')
title('Response to One Percent TFP Shock, \eta = 10')
legend('GDP','Cons','Inv')
axis([1 T_ir -.015 .10])
hold off

subplot(2,1,2)
hold on
plot(1:T_ir,IR_y_z,'k',1:T_ir,IR_v_z(:,2),'b',1:T_ir,IR_p_z,'r')
plot(1:T_ir,zeros(T_ir,1),'k')
title('Response to One Percent TFP Shock, \eta = 10')
legend('GDP','Hours','Prod','Location','Southwest')
axis([1 T_ir -.015 .02])
hold off

%Plot impulse-responses for government spending shock
figure(2)
subplot(2,1,1)
hold on
plot(1:T_ir,IR_y_g,'k',1:T_ir,IR_v_g(:,1),'b',1:T_ir,IR_i_g,'r')
plot(1:T_ir,zeros(T_ir,1),'k')
title('Response to One Percent Gov Spending Shock, \eta = 10')
legend('GDP','Cons','Inv')
axis([1 T_ir -.015 .10])
hold off

subplot(2,1,2)
hold on
plot(1:T_ir,IR_y_g,'k',1:T_ir,IR_v_g(:,2),'b',1:T_ir,IR_p_g,'r')
plot(1:T_ir,zeros(T_ir,1),'k')
title('Response to One Percent Gov Spending Shock, \eta = 10')
legend('GDP','Hours','Prod','Location','Southwest')
axis([1 T_ir -.015 .02])
hold off
