clear; clc;

%addpath c:\dynare\4.4.3\matlab
%addpath c:\dynare\4.6.3\matlab

dynare IRF_as_no_policy_old_calibration.mod
% dynare IRF_as_no_policy_old_calibration_nolog.mod

%var qk Phi_k eps1 eps2 theta  qb  Y  K   N   rk   w   Rb   I   C   BY    Lambda  Rk    spread   omegat   At    L     S   ppit  omgt  haircut theta_original lambda  Phi_k_tilde l_ppi l_omega   V;

%   [ 1    2    3   4     5     6  7  8   9   10  11   12  13   14  15     16    17      18       19     20    21    22   23     24     25         26         27        28       29      30      31 ]

% [K Y l_ppi l_omega]
ghx=oo_.dr.ghx(oo_.dr.inv_order_var,:);

H = ghx([8,7,29,30],:);

save data_H H;

H_bond = ghx(12,:);

save data_H_bond H_bond;



% run data_input_and_save
run data_input_and_save_haoning

% load data_investment_gdp_date
load data_investment_gdp_date_haoning

% investmentestment	gdp	  date	safety_premium	liquidity_premium

investment = investment_gdp_date(:,1);
gdp = investment_gdp_date(:,2);
date = investment_gdp_date(:,3);
safety = investment_gdp_date(:,4);
liquidity = investment_gdp_date(:,5);
bond = investment_gdp_date(:,6);


mean(safety)
std(safety)
mean(liquidity)
std(liquidity)


N = length(investment);
k = zeros(N,1);


delta = 0.025;

k0 = investment(1)/(delta+log(investment(6)/investment(1))/5);


for i = 1:N
    
    if i == 1
        k(i) = investment(1)/(delta+log(gdp(N)/gdp(1))/(N-1));
        
    else
        k(i) = investment(i)+(1-delta)*k(i-1);
    end
    
end
        

g_k = 100*log(k(N)/k(1))/(N-1);

g_y = 100*log(gdp(N)/gdp(1))/(N-1);

g_i = 100*log(investment(N)/investment(1))/(N-1);

growth_rate = [g_k g_y g_i];


figure(1)
plot(date,log(k),'-k',date,log(gdp),'--b','linewidth',2);
xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(N)]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
h = legend('Capital Stock (log)','Real GDP (log)');
set(h,'Interpreter','latex','fontsize',12);

print -depsc2 capital_gdp_1.eps

mean_k_i_grp = [mean(k./gdp) mean(i./gdp)];


% safety = investment_gdp_date(:,4);
% liquidity = investment_gdp_date(:,5);

w = 1600;
dev_log_gdp = log(gdp)-hpfilter(log(gdp),w);
dev_log_k = log(k)-hpfilter(log(k),w);
dev_log_investment = log(investment)-hpfilter(log(investment),w);

%dev_safety = safety-hpfilter(safety,w);

%dev_liquidity = liquidity-hpfilter(liquidity,w);

dev_safety = safety-hpfilter(safety,w);

dev_liquidity = liquidity-hpfilter(liquidity,w);

dev_bond = bond-hpfilter(bond,w);

%plot(date,safety,date,hpfilter(safety,w))

% figure(2)
% plot(lambda,x1,lambda,x2)
% print -depsc2 comparison_1.eps

% figure(2)
% NBERbc(date,[100*dev_log_gdp, 100*dev_log_k, zeros(N,1)], {'-','--','-.'}, [1.5 1.5 1], {'black','red','black'});
% xlabel('Date','Interpreter','latex','fontsize',16);
% xlim([date(1) date(N)]);
% ylim([-6 6]);
% set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
% h = legend('NBER Recession','Real GDP','Capital Stock');
% set(h,'Interpreter','latex','fontsize',12);
% 
% print -depsc2 constructed_capital_hp.eps


max_y = max(max(safety),max(liquidity))*1.08;

% 
% figure(2)
% NBERbc(date,[bond,zeros(length(date),1)], {'-','--'}, [1.5,1], {'blue','black'});
% xlabel('Date','Interpreter','latex','fontsize',16);
% xlim([date(1) date(end)]);
% %ylim([-6 6]);
% set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
% h = legend('NBER Recession','Govt Bond Yield');
% set(h,'Interpreter','latex','fontsize',12);


figure(2)
NBERbc(date,[dev_bond,zeros(length(date),1)], {'-','--'}, [1.5,1], {'blue','black'});
xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(end)]);
%ylim([-6 6]);
%set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
%h = legend('NBER Recession','Govt Bond Yield');
ylabel('Government Bond Yield','Interpreter','latex','fontsize',16);
set(h,'Interpreter','latex','fontsize',12);



figure(3)
NBERbc(date,[safety*4, liquidity*4], {'-','--'}, [1.5 1.5], {'blue','black'});
xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(N)]);
%ylim([0 max_y]);
set(gca,'linewidth',2.5,'fontsize',16,'fontname','helvetica');
h = legend('NBER Recession','Safety Premium','Liquidity Premium');
set(h,'Interpreter','latex','fontsize',16);

print -depsc2 safety_liquidity_premium_cycles.eps



% dynare IRF_as_no_policy_new_2.mod
% 
% %var qk Phi_k eps1 eps2 theta  qb  Y  K   N   rk   w   Rb   I   C   BY    Lambda  Rk    spread   omegat   At    L     S   ppit  omgt  haircut theta_original lambda  Phi_k_tilde l_ppi l_omega   V;
% 
% %   [ 1    2    3   4     5     6  7  8   9   10  11   12  13   14  15     16    17      18       19     20    21    22   23     24     25         26         27        28       29      30      31 ]
% 
% % [K Y l_ppi l_omega]
% ghx=oo_.dr.ghx(oo_.dr.inv_order_var,:);
% H = ghx([8,7,29,30],:);

%H = [0.946739	0.120404	-4.410056	0.005540; 0.180578	1.322872	-19.326121	0.012030; 0.571465	-1.263564	186.029283	0.433555; -0.044135	0.095412  30.264326	-0.797748];


hp_series = [dev_log_k(2:end),dev_log_gdp(1:end-1),dev_safety(1:end-1),dev_liquidity(1:end-1)];
save data_hp_series hp_series;

load data_H;

HH = inv(H);

backout_k_raw = HH(1,1).*hp_series(:,1)+HH(1,2).*hp_series(:,2)+HH(1,3).*hp_series(:,3)+HH(1,4).*hp_series(:,4);
backout_tfp_raw = HH(2,1).*hp_series(:,1)+HH(2,2).*hp_series(:,2)+HH(2,3).*hp_series(:,3)+HH(2,4).*hp_series(:,4);
backout_safety_raw = HH(3,1).*hp_series(:,1)+HH(3,2).*hp_series(:,2)+HH(3,3).*hp_series(:,3)+HH(3,4).*hp_series(:,4);
backout_omega_raw = HH(4,1).*hp_series(:,1)+HH(4,2).*hp_series(:,2)+HH(4,3).*hp_series(:,3)+HH(4,4).*hp_series(:,4);



% load data_pam
% % steady state of omega_bar and ppi
% omega = pam(5);
% ppi = pam(11);


backout_tfp = backout_tfp_raw*100;
backout_safety = backout_safety_raw*100;
backout_omega = backout_omega_raw*100;
% backout_safety = backout_safety_raw*100/ppi;
% backout_omega = backout_omega_raw*100/omega;

%backout_series = [backout_tfp_raw backout_safety_raw/ppi backout_omega_raw/omega];

backout_series = [backout_tfp backout_safety backout_omega];

save data_backout_series backout_series

save 'backout_series.txt' backout_series -ascii 


corr_std_tfp = [corr(backout_tfp(1:end-1),backout_tfp(2:end)) std(backout_tfp)]
corr_std_safety = [corr(backout_safety(1:end-1),backout_safety(2:end)) std(backout_safety)]
corr_std_liquidity =  [corr(backout_omega(1:end-1),backout_omega(2:end)) std(backout_omega)]

corr_std = [corr_std_tfp corr_std_safety corr_std_liquidity];

save data_corr_std corr_std;


figure(4)
NBERbc(date(1:end-1),[backout_tfp, backout_safety, backout_omega], {'-','--','-.'}, [1.5 1.5 1.5], {'blue','red','black'});
xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
h = legend('NBER Recession','A Shock','$\pi$ Shock','$\omega$ Shock');
set(h,'Interpreter','latex','fontsize',12);

print -depsc2 backout_shocks_all_in_one.eps

figure(5)
NBERbc(date(1:end-1),[dev_log_k(1:end-1)*100, backout_k_raw*100], {'-','--','-.'}, [1.5 1.5], {'black','red'});
xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
h = legend('NBER Recession','Data K','Backout K');
set(h,'Interpreter','latex','fontsize',12);

print -depsc2 backout_k_vs_data_k.eps




figure(6)
subplot(3,1,1)
NBERbc(date(1:end-1),[backout_tfp zeros(N-1,1)], {'-','-.'}, [1.5 1], {'blue','black'});
%xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(N-1)]);
ylabel('TFP Shock ($A$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
% h = legend('NBER Recession','TFP Shock (A)');
% set(h,'Interpreter','latex','fontsize',9);
subplot(3,1,2)
NBERbc(date(1:end-1),[backout_safety zeros(N-1,1)], {'-','-.'}, [1.5 1], {'blue','black'});
%xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(N-1)]);
ylabel('Safety Shock ($\pi$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
% h = legend('NBER Recession','Safety Shock ($\pi$)');
% set(h,'Interpreter','latex','fontsize',9);
subplot(3,1,3)
NBERbc(date(1:end-1),[backout_omega zeros(N-1,1)], {'-','-.'}, [1.5 1], {'blue','black'});
xlabel('Date','Interpreter','latex','fontsize',16);
xlim([date(1) date(N-1)]);
ylabel('Liquidity Shock ($\omega$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');
% h = legend('NBER Recession','Liquidity Shock ($\omega$)');
% set(h,'Interpreter','latex','fontsize',9);

print -depsc2 backout_shocks_three_panels.eps



load data_H_bond;

% H_bond = [-2.768995 2.763212 -693.581837 0.567650];

% K(-1)                       0.946876        0.181096        0.564649       -0.044929       -2.768995
% At(-1)                      0.120020        1.321358       -1.247792        0.097068        2.763212
% ppit(-1)                   -4.390167      -19.249617      184.668429       30.090999     -693.581837
% omgt(-1)                    0.005581        0.012206        0.429887
% -0.795983        0.567650




%k_all_0 = backout_k_raw(1);
k_all_0 = 0;
k_all = zeros(N-1,1);
y_all = zeros(N-1,1);
safety_all = zeros(N-1,1);
bond_all = zeros(N-1,1);
liquidity_all = zeros(N-1,1);

for i = 1:(N-1)

    if i ==1
        
        k_all(i) = H(1,1).*k_all_0+H(1,2).*backout_tfp_raw(i)+H(1,3).*backout_safety_raw(i)+H(1,4).*backout_omega_raw(i);
        y_all(i) = H(2,1).*k_all_0+H(2,2).*backout_tfp_raw(i)+H(2,3).*backout_safety_raw(i)+H(2,4).*backout_omega_raw(i);
        safety_all(i) = H(3,1).*k_all_0+H(3,2).*backout_tfp_raw(i)+H(3,3).*backout_safety_raw(i)+H(3,4).*backout_omega_raw(i);
        liquidity_all(i) = H(4,1).*k_all_0+H(4,2).*backout_tfp_raw(i)+H(4,3).*backout_safety_raw(i)+H(4,4).*backout_omega_raw(i);
        bond_all(i) = H_bond(1).*k_all_0+H_bond(2).*backout_tfp_raw(i)+H_bond(3).*backout_safety_raw(i)+H_bond(4).*backout_omega_raw(i);
        
    else
        
        k_all(i) = H(1,1).*k_all(i-1)+H(1,2).*backout_tfp_raw(i)+H(1,3).*backout_safety_raw(i)+H(1,4).*backout_omega_raw(i);
        y_all(i) = H(2,1).*k_all(i-1)+H(2,2).*backout_tfp_raw(i)+H(2,3).*backout_safety_raw(i)++H(2,4).*backout_omega_raw(i);
        safety_all(i) = H(3,1).*k_all(i-1)+H(3,2).*backout_tfp_raw(i)+H(3,3).*backout_safety_raw(i)++H(3,4).*backout_omega_raw(i);
        liquidity_all(i) = H(4,1).*k_all(i-1)+H(4,2).*backout_tfp_raw(i)+H(4,3).*backout_safety_raw(i)++H(4,4).*backout_omega_raw(i);
        bond_all(i) = H_bond(1).*k_all(i-1)+H_bond(2).*backout_tfp_raw(i)+H_bond(3).*backout_safety_raw(i)+H_bond(4).*backout_omega_raw(i);
        
    end
    
end


figure(7)
subplot(3,1,1)
NBERbc(date(1:end-1),[dev_log_gdp(1:end-1),y_all], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Output','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\all$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,2)
NBERbc(date(1:end-1),[dev_safety(1:end-1),safety_all], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Safety Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\all$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,3)
NBERbc(date(1:end-1),[dev_liquidity(1:end-1),liquidity_all], {'-','--'}, [1.5 1.5], {'black','blue'});
xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Liquidity Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\all$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

% 
% NBERbc(date(1:end-1),[dev_bond(1:end-1),bond_all], {'-','--'}, [1.5 1.5], {'black','blue'});
% xlabel('Date','Interpreter','latex','fontsize',12);
% xlim([date(1) date(N-1)]);
% ylim([-6 6]);
% 
% ylabel('Bond Yield','Interpreter','latex','fontsize',16);
% ylim([-6 6]);
% set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

print -depsc2 counterfactual_all_shocks.eps



% 
% figure(7)
% NBERbc(date,[dev_safety, dev_liquidity zeros(N,1)], {'-','--','-.'}, [1.5 1.5 1], {'blue','red','black'});
% xlabel('Date','Interpreter','latex','fontsize',16);
% xlim([date(1) date(N)]);
% %ylim([0 max_y]);
% set(gca,'linewidth',2.5,'fontsize',16,'fontname','helvetica');
% h = legend('NBER Recession','Safety Premium','Liquidity Premium',1);
% set(h,'Interpreter','latex','fontsize',12);
% 
% print -depsc2 safety_liquidity_premium_cycles_hp.eps


% Couterfactual analysis:

% (1) only tfp shock (A)

%y_tfp = H(2,1).*backout_k_raw+H(2,2).*hp_series(:,2)+H(2,3).*hp_series(:,3)+H(2,4).*hp_series(:,4);

k_tfp_0 = 0;
k_tfp = zeros(N-1,1);
y_tfp = zeros(N-1,1);
safety_tfp = zeros(N-1,1);
bond_tfp = zeros(N-1,1);
liquidity_tfp = zeros(N-1,1);


for i = 1:(N-1)

    if i ==1
        
        k_tfp(i) = H(1,1).*k_tfp_0+H(1,2).*backout_tfp_raw(i);
        y_tfp(i) = H(2,1).*k_tfp_0+H(2,2).*backout_tfp_raw(i);
        safety_tfp(i) = H(3,1).*k_tfp_0+H(3,2).*backout_tfp_raw(i);
        liquidity_tfp(i) = H(4,1).*k_tfp_0+H(4,2).*backout_tfp_raw(i);
        bond_tfp(i) = H_bond(1).*k_tfp_0+H_bond(2).*backout_tfp_raw(i);
        
        
    else
        
        k_tfp(i) = H(1,1).*k_tfp(i-1)+H(1,2).*backout_tfp_raw(i);
        y_tfp(i) = H(2,1).*k_tfp(i-1)+H(2,2).*backout_tfp_raw(i);
        safety_tfp(i) = H(3,1).*k_tfp(i-1)+H(3,2).*backout_tfp_raw(i);
        liquidity_tfp(i) = H(4,1).*k_tfp(i-1)+H(4,2).*backout_tfp_raw(i);
        bond_tfp(i) = H_bond(1).*k_tfp(i-1)+H_bond(2).*backout_tfp_raw(i);
        
    end
    
end


figure(8)
subplot(3,1,1)
NBERbc(date(1:end-1),[dev_log_gdp(1:end-1),y_tfp], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
ylabel('Output','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $A$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,2)
NBERbc(date(1:end-1),[dev_safety(1:end-1),safety_tfp], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
ylabel('Safety Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $A$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,3)
NBERbc(date(1:end-1),[dev_liquidity(1:end-1),liquidity_tfp], {'-','--'}, [1.5 1.5], {'black','blue'});
xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
ylabel('Liquidity Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $A$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

print -depsc2 counterfactual_only_tfp.eps

%y_tfp = H(2,1).*backout_k_raw+H(2,2).*hp_series(:,2)+H(2,3).*hp_series(:,3)+H(2,4).*hp_series(:,4);


% (2) only safety shock (pi)

k_ppi_0 = 0;
k_ppi = zeros(N-1,1);
y_ppi = zeros(N-1,1);
bond_ppi = zeros(N-1,1);
safety_ppi = zeros(N-1,1);
liquidity_ppi = zeros(N-1,1);



for i = 1:(N-1)

    if i ==1
        
        k_ppi(i) = H(1,1).*k_ppi_0+H(1,3).*backout_safety_raw(i);
        y_ppi(i) = H(2,1).*k_ppi_0+H(2,3).*backout_safety_raw(i);
        safety_ppi(i) = H(3,1).*k_ppi_0+H(3,3).*backout_safety_raw(i);
        liquidity_ppi(i) = H(4,1).*k_ppi_0+H(4,3).*backout_safety_raw(i);
        bond_ppi(i) = H_bond(1).*k_ppi_0+H_bond(3).*backout_safety_raw(i);
        
        
    else
        
        k_ppi(i) = H(1,1).*k_ppi(i-1)+H(1,3).*backout_safety_raw(i);
        y_ppi(i) = H(2,1).*k_ppi(i-1)+H(2,3).*backout_safety_raw(i);
        safety_ppi(i) = H(3,1).*k_ppi(i-1)+H(3,3).*backout_safety_raw(i);
        liquidity_ppi(i) = H(4,1).*k_ppi(i-1)+H(4,3).*backout_safety_raw(i);
        bond_ppi(i) = H_bond(1).*k_ppi(i-1)+H_bond(3).*backout_safety_raw(i);
        
    end
    
end


figure(9)
subplot(3,1,1)
NBERbc(date(1:end-1),[dev_log_gdp(1:end-1),y_ppi], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Output','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\pi$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,2)
NBERbc(date(1:end-1),[dev_safety(1:end-1),safety_ppi], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
ylabel('Safety Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\pi$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,3)
NBERbc(date(1:end-1),[dev_liquidity(1:end-1),liquidity_ppi], {'-','--'}, [1.5 1.5], {'black','blue'});
xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
ylabel('Liquidity Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\pi$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

print -depsc2 counterfactual_only_safety.eps

% (3) only liquidity shock (omega)

k_omega_0 = 0;
k_omega = zeros(N-1,1);
y_omega = zeros(N-1,1);
bond_omega = zeros(N-1,1);
safety_omega = zeros(N-1,1);
liquidity_omega = zeros(N-1,1);



for i = 1:(N-1)

    if i ==1
        
        k_omega(i) = H(1,1).*k_omega_0+H(1,4).*backout_omega_raw(i);
        y_omega(i) = H(2,1).*k_omega_0+H(2,4).*backout_omega_raw(i);
        safety_omega(i) = H(3,1).*k_omega_0+H(3,4).*backout_omega_raw(i);
        liquidity_omega(i) = H(4,1).*k_omega_0+H(4,4).*backout_omega_raw(i);
        bond_omega(i) = H_bond(1).*k_omega_0+H_bond(4).*backout_omega_raw(i);
        
    else
        
        k_omega(i) = H(1,1).*k_omega(i-1)+H(1,4).*backout_omega_raw(i);
        y_omega(i) = H(2,1).*k_omega(i-1)+H(2,4).*backout_omega_raw(i);
        safety_omega(i) = H(3,1).*k_omega(i-1)+H(3,4).*backout_omega_raw(i);
        liquidity_omega(i) = H(4,1).*k_omega(i-1)+H(4,4).*backout_omega_raw(i);
        bond_omega(i) = H_bond(1).*k_omega(i-1)+H_bond(4).*backout_omega_raw(i);
        
    end
    
end





figure(10)
subplot(3,1,1)
NBERbc(date(1:end-1),[dev_log_gdp(1:end-1),y_omega], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Output','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\omega$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,2)
NBERbc(date(1:end-1),[dev_safety(1:end-1),safety_omega], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Safety Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

subplot(3,1,3)
NBERbc(date(1:end-1),[dev_liquidity(1:end-1),liquidity_omega], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Liquidity Premium','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

print -depsc2 counterfactual_only_liquidity.eps


figure(11)
subplot(3,1,1)
NBERbc(date(1:end-1),[dev_bond(1:end-1),bond_tfp*100], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Bond Yield ($A$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\omega$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,2)
NBERbc(date(1:end-1),[dev_bond(1:end-1),bond_ppi*100], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Bond Yield ($\pi$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\omega$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,3)
NBERbc(date(1:end-1),[dev_bond(1:end-1),bond_omega*100], {'-','--'}, [1.5 1.5], {'black','blue'});
xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

ylabel('Bond Yield ($\omega$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\omega$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

print -depsc2 counterfactual_bond_yield.eps


figure(12)


print -depsc2 counterfactual_bond_yield_all_shocks.eps

%load data_backout_series
%[VAR, VARopt] = VARmodel(backout_series,1)
%vare(backout_series,1)
























figure(13)
title('Output','Interpreter','latex','fontsize',20);
subplot(3,1,1)
NBERbc(date(1:end-1),[dev_log_gdp(1:end-1),y_tfp], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
%ylabel('Output','Interpreter','latex','fontsize',16);
ylabel('TFP Shock ($A$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $A$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,2)
NBERbc(date(1:end-1),[dev_log_gdp(1:end-1),y_ppi], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
%ylabel('Output','Interpreter','latex','fontsize',16);
ylabel('Safety Shock ($\pi$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\pi$ Shock');
% set(h,'Interpreter','latex','fontsize',8);


subplot(3,1,3)
NBERbc(date(1:end-1),[dev_log_gdp(1:end-1),y_omega], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

%ylabel('Output','Interpreter','latex','fontsize',16);
ylabel('Liquidity Shock ($\omega$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

print -depsc2 counterfactual_output.eps





figure(14)

subplot(3,1,1)
NBERbc(date(1:end-1),[dev_safety(1:end-1),safety_tfp], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
%ylabel('Safety Premium','Interpreter','latex','fontsize',16);
ylabel('TFP Shock ($A$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $A$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,2)
NBERbc(date(1:end-1),[dev_safety(1:end-1),safety_ppi], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
%ylabel('Safety Premium','Interpreter','latex','fontsize',16);
ylabel('Safety Shock ($\pi$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\pi$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

% set(gca,'linewidth',1,'fontsize',12,'fontname','helvetica');
% h = legend('NBER Recession','GDP','Simulation by $\omega$ Shock');
% set(h,'Interpreter','latex','fontsize',8);

subplot(3,1,3)
NBERbc(date(1:end-1),[dev_safety(1:end-1),safety_omega], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
%ylabel('Safety Premium','Interpreter','latex','fontsize',16);
ylabel('Liquidity Shock ($\omega$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

print -depsc2 counterfactual_safety_premium.eps



% ylabel('TFP Shock ($A$)','Interpreter','latex','fontsize',16);
% ylabel('Safety Shock ($\pi$)','Interpreter','latex','fontsize',16);
% ylabel('Liquidity Shock ($\omega$)','Interpreter','latex','fontsize',16);



figure(15)

subplot(3,1,1)
NBERbc(date(1:end-1),[dev_liquidity(1:end-1),liquidity_tfp], {'-','--'}, [1.5 1.5], {'black','blue'});
xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
%ylabel('Liquidity Premium','Interpreter','latex','fontsize',16);
ylabel('TFP Shock ($A$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');


subplot(3,1,2)
NBERbc(date(1:end-1),[dev_liquidity(1:end-1),liquidity_ppi], {'-','--'}, [1.5 1.5], {'black','blue'});
xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);
%ylabel('Liquidity Premium','Interpreter','latex','fontsize',16);
ylabel('Safety Shock ($\pi$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');


subplot(3,1,3)
NBERbc(date(1:end-1),[dev_liquidity(1:end-1),liquidity_omega], {'-','--'}, [1.5 1.5], {'black','blue'});
%xlabel('Date','Interpreter','latex','fontsize',12);
xlim([date(1) date(N-1)]);
%ylim([-6 6]);

%ylabel('Liquidity Premium','Interpreter','latex','fontsize',16);
ylabel('Liquidity Shock ($\omega$)','Interpreter','latex','fontsize',16);
%ylim([-6 6]);
set(gca,'linewidth',1,'fontsize',16,'fontname','helvetica');

print -depsc2 counterfactual_liquidity_premium.eps
