%% This code calculates the contribution to the distribution of unemployment duration %%
% [Final] The contribution when GDD = E during the Great Recession
% separately reports the contribution of changes in GDD

clear all

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

load tl_3prm_rec_e24_final
sxi_tl = sxi;
BD_tl = BD;
RD_tl = RD;
clearvars -except sxi_tl BD_tl RD_tl

load os_3prm_rec_e24_final
sxi_os = sxi;
BD_os = BD;
RD_os = RD;
clearvars -except sxi_tl BD_tl RD_tl ...
                  sxi_os BD_os RD_os

load qt_3prm_rec_e24_final
sxi_qt = sxi;
BD_qt = BD;
RD_qt = RD;
clearvars -except sxi_tl BD_tl RD_tl ...
                  sxi_os BD_os RD_os ...
                  sxi_qt BD_qt RD_qt

load re_3prm_rec_e24_final
sxi_re = sxi;
BD_re = BD;
RD_re = RD;
clearvars -except sxi_tl BD_tl RD_tl ...
                  sxi_os BD_os RD_os ...
                  sxi_qt BD_qt RD_qt ...
                  sxi_re BD_re RD_re

load ne_3prm_rec_e24_final
sxi_ne = sxi;
BD_ne = BD;
RD_ne = RD;
clearvars -except sxi_tl BD_tl RD_tl ...
                  sxi_os BD_os RD_os ...
                  sxi_qt BD_qt RD_qt ...
                  sxi_re BD_re RD_re ...
                  sxi_ne BD_ne RD_ne U

est_sxi(:,:,1) = sxi_tl;
est_sxi(:,:,2) = sxi_os;
est_sxi(:,:,3) = sxi_qt;
est_sxi(:,:,4) = sxi_re;
est_sxi(:,:,5) = sxi_ne;

est_BD(:,:,1) = BD_tl;
est_BD(:,:,2) = BD_os;
est_BD(:,:,3) = BD_qt;
est_BD(:,:,4) = BD_re;
est_BD(:,:,5) = BD_ne;

est_RD(:,:,1) = RD_tl;
est_RD(:,:,2) = RD_os;
est_RD(:,:,3) = RD_qt;
est_RD(:,:,4) = RD_re;
est_RD(:,:,5) = RD_ne;

T = size(sxi_re,2)
load indfile
IND = ind2(size(U,1)-T+1:size(U,1),1);
load data_dur

%% 

% option set-up
% 0. no GDD; 1. same GDD (BD=RD); 2. original GDD

option = 1

% calendar set-up : 

for i=1:12
    tt(i,1)=1980+i/13;
end

for i=1:39*12
    tt(i+12,1)=tt(i,1)+1;
end

nt=48; W1 = []; W2 = []; X1 =[]; X2 = [];

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 nn = 1:5
 
sxi = est_sxi(:,:,nn);
BD = est_BD(:,:,nn);
RD = est_RD(:,:,nn);

if option == 0 
    BD = zeros(48,1);
    RD = zeros(48,1);
elseif option == 1
    RD = BD;
elseif option == 3
    BD(6:end) = ones(48-6+1,1)*BD(5,1);
    RD(6:end) = ones(48-6+1,1)*RD(5,1);
elseif option == 4
    BD(11:end) = ones(48-11+1,1)*BD(11,1);
    RD(11:end) = ones(48-11+1,1)*RD(11,1);
elseif option == 5
   BD(1:5) = ones(5,1)*BD(5,1);
   RD(1:5) = ones(5,1)*RD(5,1);   
end

for t=1:T
    
if IND(t,1)==0;
 D=BD;
else IND(t,1)==1;
 D=RD;
end

 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

for kk = 1:nt
est_PL(kk,t,nn) = kk*UL(kk,t);
est_PH(kk,t,nn) = kk*UH(kk,t);
end

end

est_TU(:,:,nn) = TU;
est_UL(:,:,nn) = UL;
est_UH(:,:,nn) = UH;

dur_PL(1,:,nn) = sum(est_PL(:,:,nn),1)./sum(est_UL(:,:,nn),1);
dur_PH(1,:,nn) = sum(est_PH(:,:,nn),1)./sum(est_UH(:,:,nn),1);

 end
 
  dur_PL = reshape(dur_PL, T, 5);
  dur_PH = reshape(dur_PH, T, 5);
 
  clearvars -except for dur_PL dur_PH tt T
  
%% count of unemployed individuals by reason for unemployment

% share and duration for the baseline case (option =2)
load baseline

% total inflows * share of obs. * share pf unobs * cumulative probability

for t=1:T
for n=1:5    
    
D_tl_H(t,1) = reason_share(t,1) * tl_share(t,1) * dur_PH(t,1);
D_tl_L(t,1) = reason_share(t,1) * tl_share(t,2) * dur_PL(t,1);

D_os_H(t,1) = reason_share(t,2) * os_share(t,1) * dur_PH(t,2);
D_os_L(t,1) = reason_share(t,2) * os_share(t,2) * dur_PL(t,2);

D_qt_H(t,1) = reason_share(t,3) * qt_share(t,1) * dur_PH(t,3);
D_qt_L(t,1) = reason_share(t,3) * qt_share(t,2) * dur_PL(t,3);

D_re_H(t,1) = reason_share(t,4) * re_share(t,1) * dur_PH(t,4);
D_re_L(t,1) = reason_share(t,4) * re_share(t,2) * dur_PL(t,4);

D_ne_H(t,1) = reason_share(t,5) * ne_share(t,1) * dur_PH(t,5);
D_ne_L(t,1) = reason_share(t,5) * ne_share(t,2) * dur_PL(t,5);

end
end

 D = (D_tl_H + D_tl_L ...
    + D_os_H + D_os_L ...
    + D_qt_H + D_qt_L ...
    + D_re_H + D_re_L ...
    + D_ne_H + D_ne_L);

 %% historical decomposition 

% calendar set-up 
 
beg = find(tt == (2007 + 12/13)); % Gr%eat Recession

% (1) counterfactual duration with the inflows fixed at December 2007
creason_share = zeros(T,5);
creason_share(beg:end,:) = ones(T-beg+1,1)*reason_share(beg,:);

% no changes

for t=1:T
    
c0_D_tl_H(t,1) = creason_share(t,1) * ctl_share(t,1) * cdur_PH(t,1);
c0_D_tl_L(t,1) = creason_share(t,1) * ctl_share(t,2) * cdur_PL(t,1);

c0_D_os_H(t,1) = creason_share(t,2) * cos_share(t,1) * cdur_PH(t,2);
c0_D_os_L(t,1) = creason_share(t,2) * cos_share(t,2) * cdur_PL(t,2);

c0_D_qt_H(t,1) = creason_share(t,3) * cqt_share(t,1) * cdur_PH(t,3);
c0_D_qt_L(t,1) = creason_share(t,3) * cqt_share(t,2) * cdur_PL(t,3);

c0_D_re_H(t,1) = creason_share(t,4) * cre_share(t,1) * cdur_PH(t,4);
c0_D_re_L(t,1) = creason_share(t,4) * cre_share(t,2) * cdur_PL(t,4);

c0_D_ne_H(t,1) = creason_share(t,5) * cne_share(t,1) * cdur_PH(t,5);
c0_D_ne_L(t,1) = creason_share(t,5) * cne_share(t,2) * cdur_PL(t,5);

end

c0_D = (c0_D_tl_H + c0_D_tl_L + ...
        c0_D_os_H + c0_D_os_L + ...
        c0_D_qt_H + c0_D_qt_L + ...
        c0_D_re_H + c0_D_re_L + ...
        c0_D_ne_H + c0_D_ne_L);

% observed groups

for t=1:T
    
co_D_tl_H(t,1) = reason_share(t,1) * ctl_share(t,1) * cdur_PH(t,1);
co_D_tl_L(t,1) = reason_share(t,1) * ctl_share(t,2) * cdur_PL(t,1);

co_D_os_H(t,1) = reason_share(t,2) * cos_share(t,1) * cdur_PH(t,2);
co_D_os_L(t,1) = reason_share(t,2) * cos_share(t,2) * cdur_PL(t,2);

co_D_qt_H(t,1) = reason_share(t,3) * cqt_share(t,1) * cdur_PH(t,3);
co_D_qt_L(t,1) = reason_share(t,3) * cqt_share(t,2) * cdur_PL(t,3);

co_D_re_H(t,1) = reason_share(t,4) * cre_share(t,1) * cdur_PH(t,4);
co_D_re_L(t,1) = reason_share(t,4) * cre_share(t,2) * cdur_PL(t,4);

co_D_ne_H(t,1) = reason_share(t,5) * cne_share(t,1) * cdur_PH(t,5);
co_D_ne_L(t,1) = reason_share(t,5) * cne_share(t,2) * cdur_PL(t,5);

end

co_D = (co_D_tl_H + co_D_tl_L + ...
        co_D_os_H + co_D_os_L + ...
        co_D_qt_H + co_D_qt_L + ...
        co_D_re_H + co_D_re_L + ...
        co_D_ne_H + co_D_ne_L);
      
% unobserved types

for t=1:T
    
cu_D_tl_H(t,1) = reason_share(t,1) * tl_share(t,1) * cdur_PH(t,1);
cu_D_tl_L(t,1) = reason_share(t,1) * tl_share(t,2) * cdur_PL(t,1);

cu_D_os_H(t,1) = reason_share(t,2) * os_share(t,1) * cdur_PH(t,2);
cu_D_os_L(t,1) = reason_share(t,2) * os_share(t,2) * cdur_PL(t,2);

cu_D_qt_H(t,1) = reason_share(t,3) * qt_share(t,1) * cdur_PH(t,3);
cu_D_qt_L(t,1) = reason_share(t,3) * qt_share(t,2) * cdur_PL(t,3);

cu_D_re_H(t,1) = reason_share(t,4) * re_share(t,1) * cdur_PH(t,4);
cu_D_re_L(t,1) = reason_share(t,4) * re_share(t,2) * cdur_PL(t,4);

cu_D_ne_H(t,1) = reason_share(t,5) * ne_share(t,1) * cdur_PH(t,5);
cu_D_ne_L(t,1) = reason_share(t,5) * ne_share(t,2) * cdur_PL(t,5);

end

cu_D = (cu_D_tl_H + cu_D_tl_L + ...
        cu_D_os_H + cu_D_os_L + ...
        cu_D_qt_H + cu_D_qt_L + ...
        cu_D_re_H + cu_D_re_L + ...
        cu_D_ne_H + cu_D_ne_L);
    
% type L duration

for t=1:T
    
cl_D_tl_H(t,1) = reason_share(t,1) * tl_share(t,1) * cdur_PH(t,1);
cl_D_tl_L(t,1) = reason_share(t,1) * tl_share(t,2) * dur_PL(t,1);

cl_D_os_H(t,1) = reason_share(t,2) * os_share(t,1) * cdur_PH(t,2);
cl_D_os_L(t,1) = reason_share(t,2) * os_share(t,2) * dur_PL(t,2);

cl_D_qt_H(t,1) = reason_share(t,3) * qt_share(t,1) * cdur_PH(t,3);
cl_D_qt_L(t,1) = reason_share(t,3) * qt_share(t,2) * dur_PL(t,3);

cl_D_re_H(t,1) = reason_share(t,4) * re_share(t,1) * cdur_PH(t,4);
cl_D_re_L(t,1) = reason_share(t,4) * re_share(t,2) * dur_PL(t,4);

cl_D_ne_H(t,1) = reason_share(t,5) * ne_share(t,1) * cdur_PH(t,5);
cl_D_ne_L(t,1) = reason_share(t,5) * ne_share(t,2) * dur_PL(t,5);

end

cl_D = (cl_D_tl_H + cl_D_tl_L + ...
        cl_D_os_H + cl_D_os_L + ...
        cl_D_qt_H + cl_D_qt_L + ...
        cl_D_re_H + cl_D_re_L + ...
        cl_D_ne_H + cl_D_ne_L);

    % type H duration
for t=1:T
    
ch_D_tl_H(t,1) = reason_share(t,1) * tl_share(t,1) * dur_PH(t,1);
ch_D_tl_L(t,1) = reason_share(t,1) * tl_share(t,2) * dur_PL(t,1);

ch_D_os_H(t,1) = reason_share(t,2) * os_share(t,1) * dur_PH(t,2);
ch_D_os_L(t,1) = reason_share(t,2) * os_share(t,2) * dur_PL(t,2);

ch_D_qt_H(t,1) = reason_share(t,3) * qt_share(t,1) * dur_PH(t,3);
ch_D_qt_L(t,1) = reason_share(t,3) * qt_share(t,2) * dur_PL(t,3);

ch_D_re_H(t,1) = reason_share(t,4) * re_share(t,1) * dur_PH(t,4);
ch_D_re_L(t,1) = reason_share(t,4) * re_share(t,2) * dur_PL(t,4);

ch_D_ne_H(t,1) = reason_share(t,5) * ne_share(t,1) * dur_PH(t,5);
ch_D_ne_L(t,1) = reason_share(t,5) * ne_share(t,2) * dur_PL(t,5);

end

ch_D = (ch_D_tl_H + ch_D_tl_L + ...
        ch_D_os_H + ch_D_os_L + ...
        ch_D_qt_H + ch_D_qt_L + ...
        ch_D_re_H + ch_D_re_L + ...
        ch_D_ne_H + ch_D_ne_L); 
    
hist_dur = [D(beg:end,1), co_D(beg:end,1), cu_D(beg:end,1) - co_D(beg:end,1), ...
     cl_D(beg:end,1) - cu_D(beg:end,1), ch_D(beg:end,1)-cl_D(beg:end,1)];

% this step is needed because the starting points are different due to GDD
hist_dur2 = hist_dur - ones(size(hist_dur,1),1)*hist_dur(1,:); 

%plot(hist_dur2)

%% plot : duration decomposition

nm = 48;
plot_beg = beg ;
plot_end = plot_beg + nm ;

t2007r=tt(plot_beg:plot_end,1); % timeline 2007M11-2011M11 (Recession: 2007M12-2009M06)
hist_dur = hist_dur2(1:nm+1,:);

load data_duration

b = t2007r(1,1);
e = t2007r(end,1);

%{
g2007r = figure;
hold on
set(gca,'fontsize',15)
title('Great Recession and its recovery (Dec.2007-Dec.2011)')
plot(t2007r,(duration(plot_beg:plot_end,1)-duration(plot_beg,1))/4.33,'--','LineWidth',1,'Color',[0 0 0]) % black
plot(t2007r,hist_dur(:,1),'-','LineWidth',3,'Color',[0 0 0]) % black
plot(t2007r,hist_dur(:,2),'-','LineWidth',1,'Color',[0 0 0]) % black
plot(t2007r,hist_dur(:,3),'-','LineWidth',3,'Color',[1 0 0]) % red
plot(t2007r,hist_dur(:,4),'--','LineWidth',2,'Color',[0 0 1]) % magenta
plot(t2007r,hist_dur(:,5),'-x','LineWidth',2,'Color',[1 0 1]) % blue 
hleg1 = legend('Data','D_{t}^{full}','D_t^{o}', ...
               'D_{t}^{u}','D_{Lt}','D_{Ht}');
set(hleg1,'Location', 'NorthWest') 
colorstr=[159 182 190]/205;
hold off
%}

G = sum(hist_dur_full(:,4:5),2) - sum(hist_dur(:,4:5),2);

hist_dur = hist_dur;

g2007r = figure;
hold on
set(gca,'fontsize',15)
title('Great Recession and its recovery (Dec.2007-Dec.2011)')
plot(t2007r,hist_dur_full(:,1),'-','LineWidth',3,'Color',[0 0 0]) % black
plot(t2007r,hist_dur(:,2),'-','LineWidth',1,'Color',[0 0 0]) % black
plot(t2007r,hist_dur(:,3),'-','LineWidth',3,'Color',[1 0 0]) % red
plot(t2007r,hist_dur(:,4),'--','LineWidth',2,'Color',[0 0 1]) % magenta
plot(t2007r,hist_dur(:,5),'-x','LineWidth',2,'Color',[1 0 1]) % blue 
plot(t2007r,G(:,1),'--','LineWidth',2,'Color',[0 0 0]) % blue 
hleg1 = legend('D_{t}^{full}','D_t^{o}', ...
               'D_{t}^{u}','D_{Lt}^{E}','D_{Ht}^{E}','G_t');
set(hleg1,'Location', 'NorthWest') 
axis([b,e,-0.5,4.5])
%axis([b,e,3.5,9.0])
set(gca,'XTick',[2008.1, 2009.1, 2010.1, 2011.1])
colorstr=[159 182 190]/205;
shade(t2007r(1,1),t2007r(19,1),colorstr); % Dec 2007 - June 2009
hold off
saveas(g2007r,'/cmc/ahn-projects/A1_final/figures/hist_dur','epsc');

