close all

fonts1=10;
fonts2=5;

ynewsshock=irf_fevd_out_2.eps_structural(:,3);
rvolshock=irf_fevd_out_4.eps_structural(:,1);
bcnewsshock=irf_fevd_out_5.eps_structural(:,1);

shocks_mat=[ynewsshock,...
    rvolshock,...
    bcnewsshock];

shock_name_vec={'Yield News Shock'
    'Realized Volatility Shock'
    'BC News Shock'};

shocks_ma_rho90=[];
for i=1:size(shocks_mat,2)
    shocki=shocks_mat(:,i);
    shocks_ma_rho90=[shocks_ma_rho90,weighted_ma(shocki,.9)];    
end

%-- plot shocks and exponentially weighted moving averages
r=0.9;
y=data_raw(:,ismember(name_vec,{'INDPRO'}));
y=100*(log(y)-log(lag(y,11))); % 12-month growth in IP
dates=datenum([calds_m,ones(size(shocks_ma_rho90,1),1)]);


% *********************************************************************
% Figure 5 -- Time series of estimated yield news shock
% *********************************************************************
namevectmp={'Yield News Shock'
    'Realized Volatility Shock'
    'BC News Shock'};

shocknameveccolor={'Yield News'
    'Realized Volatility'
    'BC News'};
NM=length(namevectmp);
indshocks=nan(1,NM);
for i=1:NM
    indshocks(i)=find(ismember(shock_name_vec,namevectmp{i}));
end

ynewsshock2=packr(ynewsshock);
n_sims=500; 
st=std(ynewsshock2);
ksk=(prctile(ynewsshock2,90)+prctile(ynewsshock2,10)-2*prctile(ynewsshock2,50))./(prctile(ynewsshock2,90)-prctile(ynewsshock2,10));
sk=skewness(ynewsshock2);
kur=kurtosis(ynewsshock2);

n_obs=length(ynewsshock2);
bs = datasample(ynewsshock2,n_obs*n_sims);
bs2 = reshape(bs,n_obs,n_sims);

se_st=std(std(bs2));
ksk_bs = (prctile(bs2,90)+prctile(bs2,10)-2*prctile(bs2,50))./(prctile(bs2,90)-prctile(bs2,10));
se_ksk = std(ksk_bs);
se_sk = std(skewness(bs2));
se_kur = std(kurtosis(bs2));

fig5_1=figure;
bar(dates,ynewsshock,'FaceColor',[1, 0, 0],'EdgeColor',[1, 0, 0],'LineWidth',.5);
ax = gca;
ax.XTick = dates(1:48:end);
datetick('x','yyyy','keeplimits')
xlim([dates(1),dates(end)])
recessionplot;
title('Yield News Shock','fontsize',fonts1)
text(725180,-3,sprintf('Skewness, Kelley Skewness = %.2f, %.2f',sk,ksk),'fontsize',fonts1)
text(733850,-3.5,sprintf('(%.2f) (%.2f)',se_sk,se_ksk),'fontsize',fonts1)
text(725180,-4.5,sprintf('Kurtosis                                  = %.2f',kur),'fontsize',fonts1)
text(733850,-5,sprintf('(%.2f)',se_kur),'fontsize',fonts1)


corrmat=corr(packr([shocks_ma_rho90(:,indshocks),y]));

fig5_2=figure;
plot(dates,shocks_ma_rho90(:,1),'-','LineWidth',1,'Color', [1, 0, 0]);
ax = gca;
ax.XTick = dates(1:48:end);
datetick('x','yyyy','keeplimits')
xlim([dates(1),dates(end)])
recessionplot;
title(['Exponentially Weighted Moving Average of Yield News Shock, \rho = ',sprintf('%.2f',r)],'fontsize',fonts1)
text(728686,-7,sprintf('Correlation with IP = %.2f',corrmat(end,1)),'fontsize',fonts1)


% *********************************************************************
% Figure 12 -- Exponentially weighted moving average of estimated 
% realized volatility and business cycle news shocks
% *********************************************************************
fig12_1=figure;
hold on
p_ynews=plot(dates,shocks_ma_rho90(:,1),'-','LineWidth',.1,'Color', [1, 0, 0]);
p_rvol=plot(dates,shocks_ma_rho90(:,2),'-','LineWidth',1,'Color', [0, 0.4470, 0.7410]);
hold off
ax = gca;
ax.XTick = dates(1:48:end);
datetick('x','yyyy','keeplimits')
xlim([dates(1),dates(end)])
recessionplot;
title(['Exponentially Weighted Moving Average of Realized Volatility Shock, \rho = ',sprintf('%.2f',r)],'fontsize',fonts1)
text(728686,-7,sprintf('Correlation with IP = %.2f',corrmat(end,2)),'fontsize',fonts1)
legend([p_rvol,p_ynews],{'Realized Volatility','Yield News'},...
            'location','northwest','fontsize',fonts2)


fig12_2=figure;
hold on
p_ynews=plot(dates,shocks_ma_rho90(:,1),'-','LineWidth',.1,'Color', [1, 0, 0]);
p_bcnews=plot(dates,shocks_ma_rho90(:,3),'-','LineWidth',1,'Color', [0.3922, 0.8314, 0.0745]);
hold off
ax = gca;
ax.XTick = dates(1:48:end);
datetick('x','yyyy','keeplimits')
xlim([dates(1),dates(end)])
recessionplot;
title(['Exponentially Weighted Moving Average of Business Cycle News Shock, \rho = ',sprintf('%.2f',r)],'fontsize',fonts1)
text(728686,-7,sprintf('Correlation with IP = %.2f',corrmat(end,3)),'fontsize',fonts1)
legend([p_bcnews,p_ynews],{'Business Cycle News','Yield News'},...
            'location','northwest','fontsize',fonts2)


         
%-------- save figures
saveas(fig5_1,sprintf('results/Fig5_yield_news_shock.eps'),'epsc')
saveas(fig5_1,sprintf('results/Fig5_yield_news_shock.pdf'),'pdf')
saveas(fig5_2,sprintf('results/Fig5_yield_news_shock_filtered.eps'),'epsc')
saveas(fig5_2,sprintf('results/Fig5_yield_news_shock_filtered.pdf'),'pdf')
saveas(fig12_1,sprintf('results/Fig12_rvol_shock_filtered.eps'),'epsc')
saveas(fig12_1,sprintf('results/Fig12_rvol_shock_filtered.pdf'),'pdf')
saveas(fig12_2,sprintf('results/Fig12_bc_news_shock_filtered.eps'),'epsc')
saveas(fig12_2,sprintf('results/Fig12_bc_news_shock_filtered.pdf'),'pdf')