%% Figures 5/6/7/8 %%

clear all
load result_baseline_3prm

nt=48;PS = smPU; P = ones(T,2)-PS;W = smW;

tau = zeros(nt,nt);
 for i=1:nt
 for j=2:nt
     k=i-j+1;
     if k>0
     tau(i,j-1)=k;
     else 
     tau(i,j-1)=0;
     end
 end
 end
 
for t=1:T
    
 for i=1:nt
 W1(i,t) = sxi(1+4*(i-1),t);
 W2(i,t) = sxi(2+4*(i-1),t);
 X1(i,t) = sxi(3+4*(i-1),t);
 X2(i,t) = sxi(4+4*(i-1),t);
 end

 for i=2:nt
 for j=1:nt            
     if tau(i,j)>0
 MW1(i,j,t) = exp(-exp(X1(j,t)+D(tau(i,j),1)));
 MW2(i,j,t) = exp(-exp(X2(j,t)+D(tau(i,j),1)));
     else
 MW1(i,j,t) = 1;
 MW2(i,j,t) = 1;
 end
 end
 end
 
for j=1:nt
 D_MW1(j,j,t)= prod(MW1(j,:,t),2);
 D_MW2(j,j,t)= prod(MW2(j,:,t),2);
end
 
 D_MW1(1,1,t)= 1;
 D_MW2(1,1,t)= 1;

 UL(1,t)=W1(1,t)*D_MW1(1,1,t);
 UH(1,t)=W2(1,t)*D_MW2(1,1,t);
 TU(1,t)=UL(1,t)+UH(1,t);

for k=2:nt
 UL(k,t)=W1(k,t)*D_MW1(k,k,t);
 UH(k,t)=W2(k,t)*D_MW2(k,k,t);
 TU(k,t)=UL(k,t)+UH(k,t);
end
end

TUH = sum(UH,1)';TUL = sum(UL,1)';
clearvars TS;
TS(:,1) = TUH'./(TUH'+TUL');
TS(:,2) = TUL'./(TUH'+TUL');

%% confidence intevals calculation %%

sxi_base = sxi;

program = @ekf_chv_3prm_cfband

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

for i=1:nn
mle1(:,i) = mle + (0.1^8)*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

 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

% Ansley and Kohn's Variance

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)';
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));

% confidence intervals for smoothed inflows and probabilities

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

%% confidence intervals for GDD

se11 = table4_se(5,1);
se12 = table4_se(6,1);
se13 = table4_se(7,1);

CI_D = zeros(nt,1); % matrix of duration dependence parameter

CI_D(1,1)=0;

CI_D(2,1)=se11*(1/48) + se12*[2*(1/48)^2-1] + se13*[4*(1/48)^3-3*(1/48)];
CI_D(3,1)=se11*(2/48) + se12*[2*(2/48)^2-1] + se13*[4*(2/48)^3-3*(2/48)];
CI_D(4,1)=se11*(3/48) + se12*[2*(3/48)^2-1] + se13*[4*(3/48)^3-3*(3/48)];
CI_D(5,1)=se11*(4/48) + se12*[2*(4/48)^2-1] + se13*[4*(4/48)^3-3*(4/48)];
CI_D(6,1)=se11*(5/48) + se12*[2*(5/48)^2-1] + se13*[4*(5/48)^3-3*(5/48)];
CI_D(7,1)=se11*(6/48) + se12*[2*(6/48)^2-1] + se13*[4*(6/48)^3-3*(6/48)];
CI_D(8,1)=se11*(7/48) + se12*[2*(7/48)^2-1] + se13*[4*(7/48)^3-3*(7/48)];
CI_D(9,1)=se11*(8/48) + se12*[2*(8/48)^2-1] + se13*[4*(8/48)^3-3*(8/48)];
CI_D(10,1)=se11*(9/48) + se12*[2*(9/48)^2-1] + se13*[4*(9/48)^3-3*(9/48)];
CI_D(11,1)=se11*(10/48) + se12*[2*(10/48)^2-1] + se13*[4*(10/48)^3-3*(10/48)];

for i=12:nt
CI_D(i,1)=se11*((i-2)/48) + se12*[2*((i-2)/48)^2-1] + se13*[4*((i-2)/48)^3-3*((i-2)/48)];
end %duration 13+ m

% dates

for i=1:12
    tt(i,1)=1976+i/13;
end
for i=1:37*16;
    tt(i+12,1)=tt(i,1)+1;
end
% time line 1978.01 - 2013.12
t = tt(25:T,1); 

%% Figure 5
g1 = figure;
hold on
plot(t,smPU(25:T,1),'-','LineWidth',2,'Color',[1 0 0])
plot(t,smPU(25:T,2),'--','LineWidth',2,'Color',[0 0 1])
plot(t,smPL_CI(25:T,1),':','LineWidth',2,'Color',[1 0 0])
plot(t,smPL_CI(25:T,2),':','LineWidth',2,'Color',[1 0 0])
plot(t,smPH_CI(25:T,1),':','LineWidth',2,'Color',[0 0 1])
plot(t,smPH_CI(25:T,2),':','LineWidth',2,'Color',[0 0 1])
hleg1 = legend('p_{L,t}(1)','p_{H,t}(1)');
set(hleg1,'Location', 'Southwest') 
axis([1978.1,2017.5,0.2,1])
colorstr=[159 182 205]/256;
shade(t(25,1),t(31,1),colorstr);
shade(t(43,1),t(59,1),colorstr); 
shade(t(151,1),t(159,1),colorstr); 
shade(t(279,1),t(287,1),colorstr); 
shade(t(360,1),t(378,1),colorstr);
hold off
saveas(g1,'Figure5','epsc');

% Figure 6
g2 = figure;
hold on
plot(t,smW(25:T,1),'-','LineWidth',2,'Color',[1 0 0])
plot(t,smW(25:T,2),'--','LineWidth',2,'Color',[0 0 1])
plot(t,smWL_CI(25:T,1),':','LineWidth',2,'Color',[1 0 0])
plot(t,smWL_CI(25:T,2),':','LineWidth',2,'Color',[1 0 0])
plot(t,smWH_CI(25:T,1),':','LineWidth',2,'Color',[0 0 1])
plot(t,smWH_CI(25:T,2),':','LineWidth',2,'Color',[0 0 1])
hleg1 = legend('W_L','W_H');
set(hleg1,'Position', [0.79,0.82,0.05,0.05]) 
axis([1978.1,2017.5,0.2,3.5])
colorstr=[159 182 205]/256;
shade(t(25,1),t(31,1),colorstr);
shade(t(43,1),t(59,1),colorstr); 
shade(t(151,1),t(159,1),colorstr); 
shade(t(279,1),t(287,1),colorstr); 
shade(t(360,1),t(378,1),colorstr);
hold off
saveas(g2,'Figure6','epsc');

% Figure 7
g3 = figure;
hold on
plot(t,TS(25:T,1),'-','LineWidth',2,'Color',[1 0 0])
plot(t,TS(25:T,2),':','LineWidth',2,'Color',[0 0 1])
hleg1 = legend('Share_H','Share_L');
set(hleg1,'Location','Southwest') 
axis([1978.1,2017.5,0,1])
colorstr=[159 182 205]/256;
shade(t(25,1),t(31,1),colorstr);
shade(t(43,1),t(59,1),colorstr); 
shade(t(151,1),t(159,1),colorstr); 
shade(t(279,1),t(287,1),colorstr); 
shade(t(360,1),t(378,1),colorstr);
hold off
saveas(g3,'Figure7','epsc');

% Figure 8: graphs for duration dependence parameters

dd =  figure;
hold on
plot(D(1:24,1),'-','LineWidth',2,'Color',[0 0 0]) % black
plot([D(1:24,1)+2*CI_D(1:24,:)],':','LineWidth',2,'Color',[0 0 0]) % black
plot([D(1:24,1)-2*CI_D(1:24,:)],':','LineWidth',2,'Color',[0 0 0]) % black
hold off
saveas(dd,'Figure8_a','epsc');

for i=1:24
DH(i,1) = exp(-exp(D(i,1)+mean(sxi(4,:),2)));
CI_DH(i,1) = exp(-exp(D(i,1)+2*CI_D(i,:)+mean(sxi(4,:),2)));
CI_DH(i,2) = exp(-exp(D(i,1)-2*CI_D(i,:)+mean(sxi(4,:),2)));
DL(i,1) = exp(-exp(D(i,1)+mean(sxi(3,:),2)));
CI_DL(i,1) = exp(-exp(D(i,1)+2*CI_D(i,:)+mean(sxi(3,:),2)));
CI_DL(i,2) = exp(-exp(D(i,1)-2*CI_D(i,:)+mean(sxi(3,:),2)));
end

dd_all =  figure;
hold on
plot(DL(1:24,:),'-','LineWidth',2,'Color',[1 0 0]) % red
plot(DH(1:24,:),'--','LineWidth',2,'Color',[0 0 1]) % red
plot(CI_DL(1:24,1),':','LineWidth',2,'Color',[1 0 0]) % red
plot(CI_DL(1:24,2),':','LineWidth',2,'Color',[1 0 0]) % red
plot(CI_DH(1:24,1),':','LineWidth',2,'Color',[0 0 1]) % red
plot(CI_DH(1:24,2),':','LineWidth',2,'Color',[0 0 1]) % red
axis([1,24,0.1,1])
hleg1 = legend('type L','type H');
set(hleg1,'Location','SouthWest') 
hold off
saveas(dd_all,'Figure8_b','epsc');

%% Figure 12 %%
% time line 1978.01 - 2013.12
load udata_rs;

t_rs = tt(12:T,1); 

g41 = figure;
hold on
plot(t,smW(25:T,1),'-','LineWidth',2,'Color',[1 0 0]) % reason duration data avaialble from 1976.06
plot(t,u_rs(25:T,1)/1000,':','LineWidth',2,'Color',[0 0 1])
ylabel('(Millions)');
hleg1 = legend('W_L','Permanent Separations');
set(hleg1,'Position', [0.32,0.80,0.1,0.1]) 
axis([1977.385,2017.5,0,2])
colorstr=[159 182 205]/256;
shade(t(25,1),t(31,1),colorstr);
shade(t(43,1),t(59,1),colorstr); 
shade(t(151,1),t(159,1),colorstr); 
shade(t(279,1),t(287,1),colorstr); 
shade(t(360,1),t(378,1),colorstr);
hold off
saveas(g41,'Figure14_c','epsc');

g42 = figure;
hold on
plot(t,TUL(25:T,1),'-','LineWidth',2,'Color',[1 0 0])
plot(t,u_rs(25:T,2)/1000,':','LineWidth',2,'Color',[0 0 1])
ylabel('(Millions)');
hleg1 = legend('Type L','Permanent Separations');
set(hleg1,'Position', [0.32,0.80,0.1,0.1]) 
axis([1977.385,2017.5,0,14])
colorstr=[159 182 205]/256;
shade(t(25,1),t(31,1),colorstr);
shade(t(43,1),t(59,1),colorstr); 
shade(t(151,1),t(159,1),colorstr); 
shade(t(279,1),t(287,1),colorstr); 
shade(t(360,1),t(378,1),colorstr);
hold off
saveas(g42,'Figure14_d','epsc');