%% Program to compute the confidence bands of smoothed estimates in Figures 5 and 6 %% 

clear all

addpath('../result_files')
addpath('../data')

load('ne_3prm_rec_e24_final.mat')
load civpop.mat

global T;global Y;global xi0; global xi1; global P0; global P1; global F;
global d11;global d12;global d13;global d14;global Q;global R;global BD;global RD;
global r11;global r12;global r13;global r14;global ind65;
load indfile
ind65 = ind2(size(U,1)-T+1:size(U,1),1);

sxi_base = sxi;

%program = @ekf_3prm_cfband % for all
program = @ekf_3prm_cfband_tl % for TL

mle = x0;
nn = size(mle,1);
eyemat = eye(nn);

for i=1:nn
mle1(:,i) = mle + (0.1^5)*eyemat(:,i);
end

for i=1:nn
sxi_prime(:,:,i) =feval(program,mle1(:,i));
sxi_diff(:,:,i) = (sxi_prime(:,:,i) -sxi_base);%./(0.1^8);
end

%% Variance of smoothed estimates %%
for i=1:T
    J(:,:,i) = P1(:,:,i)*F'*inv(P0(:,:,i+1));
end

% corresponding mean squared errors %
 SP(:,:,T) = P1(:,:,T);
 SP(:,:,T-1) = P0(:,:,T-1) + J(:,:,T-1)*(P1(:,:,T)-P0(:,:,T))*J(:,:,T-1)';
for i=2:T-1
    SP(:,:,T-i) = P1(:,:,T-i)+J(:,:,T-i)*(SP(:,:,T+1-i)-P0(:,:,T+1-i))*J(:,:,T-i)';
end

%%

for k=1:nn
V(k,k) = white_SE(k,2)^2;
end

for t=1:T
 cfband(:,:,t) = SP(:,:,t) + ...
    reshape(sxi_diff(:,t,:), 192, nn)*V*reshape(sxi_diff(:,t,:), 192, nn)';
% tmp(:,:,t) = reshape(sxi_diff(:,t,:), 192, nn)*V*reshape(sxi_diff(:,t,:), 192, nn)';
end

cfband_wL = sqrt(reshape(cfband(1,1,:),T,1));
cfband_wH = sqrt(reshape(cfband(2,2,:),T,1));
cfband_xL = sqrt(reshape(cfband(3,3,:),T,1));
cfband_xH = sqrt(reshape(cfband(4,4,:),T,1));

for t=1:T
smPL_CI(t,1) = exp(-exp(sxi(3,t)-1.96*cfband_xL(t,1)));
smPL_CI(t,2) = exp(-exp(sxi(3,t)+1.96*cfband_xL(t,1)));
smPH_CI(t,1) = exp(-exp(sxi(4,t)-1.96*cfband_xH(t,1)));
smPH_CI(t,2) = exp(-exp(sxi(4,t)+1.96*cfband_xH(t,1)));

smWL_CI(t,1) = sxi(1,t)-2*cfband_wL(t,1);
smWL_CI(t,2) = sxi(1,t)+2*cfband_wL(t,1);
smWH_CI(t,1) = sxi(2,t)-2*cfband_wH(t,1);
smWH_CI(t,2) = sxi(2,t)+2*cfband_wH(t,1);
end
% dates

for i=1:12
    tt(i,1)=1976+i/13;
end
for i=1:43*12;
    tt(i+12,1)=tt(i,1)+1;
end
% time line 1978.01 - 2013.12
t = tt(48:size(U,1),1); 

%% Figure 5
%{
g1 = figure;
hold on
set(gca,'fontsize',25)
title('New entrants to the labor force')
plot(t,smPU(1:T,1),'-','LineWidth',1,'Color',[1 0 0])
plot(t,smPU(1:T,2),'-','LineWidth',1,'Color',[0 0 1])
plot(t,smPL_CI(1:T,1),'--','LineWidth',1,'Color',[1 0 0])
plot(t,smPL_CI(1:T,2),'--','LineWidth',1,'Color',[1 0 0])
plot(t,smPH_CI(1:T,1),'--','LineWidth',1,'Color',[0 0 1])
plot(t,smPH_CI(1:T,2),'--','LineWidth',1,'Color',[0 0 1])
hleg1 = legend('p_{L,t}(1)','p_{H,t}(1)','Orientation','Horizontal');
%hleg1 = legend('p_{L,t}(1)','p_{H,t}(1)');
set(hleg1,'Location', 'Southeast') 
axis([1979.9,2019.8,0,1])
colorstr=[159 182 205]/256;
shade(t(2,1),t(8,1),colorstr);
shade(t(20,1),t(36,1),colorstr); 
shade(t(128,1),t(136,1),colorstr); 
shade(t(256,1),t(264,1),colorstr); 
shade(t(337,1),t(355,1),colorstr);
hold off
saveas(g1,'g3_ne','epsc');
%}
% Figure 6

g1 = figure;
hold on
set(gca,'fontsize',25)
title('New entrants to the labor force')
plot(t,smW(1:T,1)./pop*100000,'-','LineWidth',1,'Color',[1 0 0])
plot(t,smW(1:T,2)./pop*100000,'-','LineWidth',1,'Color',[0 0 1])
%plot(t,smWL_CI(1:T,1),'--','LineWidth',1,'Color',[1 0 0])
%plot(t,smWL_CI(1:T,2),'--','LineWidth',1,'Color',[1 0 0])
%plot(t,smWH_CI(1:T,1),'--','LineWidth',1,'Color',[0 0 1])
%plot(t,smWH_CI(1:T,2),'--','LineWidth',1,'Color',[0 0 1])
hleg1 = legend('W_L','W_H','Orientation','Horizontal');
set(hleg1,'Location', 'Northeast') 
axis([1979.9,2019.8,0,0.8])
colorstr=[159 182 205]/256;
% shade(t(2,1),t(8,1),colorstr);
% shade(t(20,1),t(36,1),colorstr); 
% shade(t(128,1),t(136,1),colorstr); 
% shade(t(256,1),t(264,1),colorstr); 
% shade(t(337,1),t(355,1),colorstr);
hold off
saveas(g1,'g1_ne_pop','epsc');

%{
g1 = figure;
hold on
set(gca,'fontsize',25)
title('Job Leavers')
plot(t,smW(1:T,1)./sum(smW(1:T,:),2),'-','LineWidth',4,'Color',[1 0 0])
plot(t,smW(1:T,2)./sum(smW(1:T,:),2),'-','LineWidth',4,'Color',[0 0 1])
hleg1 = legend('W_L','W_H','Orientation','Horizontal');
set(hleg1,'Location', 'Northeast') 
axis([1979.9,2019.8,0,1.1])
colorstr=[159 182 205]/256;
shade(t(2,1),t(8,1),colorstr);
shade(t(20,1),t(36,1),colorstr); 
shade(t(128,1),t(136,1),colorstr); 
shade(t(256,1),t(264,1),colorstr); 
shade(t(337,1),t(355,1),colorstr);
hold off
%saveas(g1,'g2_tl','epsc');
%}
%save('tl_3prm_rec_e24_recent.mat')
