%% This code calculates the contribution to the distribution of unemployment duration %%
% Note: The related part is page 21, and footnotes 17 and 18. 
% This code considers various cases of 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

%% Simulation
% Option/baseline set-up

% A. Base set-up
% 1.RD=BD; 2. BD=RD; 3. nochange 
base = 0;

% B. Option set-up
% 0. no GDD; 1. RD=BD; 2. original GDD; 3. net-GDD(LTU); 
% 4. -GDD(LTU); 5. +GDD(LTU)

%option = 1
option = 5;

% contribution of negative GDD in LTU (duration in option 4-option 3): option = 3;
% contribution of positive GDD in LTU (duration in option 2-option 4): option = 4;
% contribution of GDD in LTU on net (duration in option 2 - option 3)
% contribution of short-term positive GDD (duration in option 5 - option 2): option = 5;

% 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 base == 1
    RD = BD;
elseif base == 2
    BD = RD;
end

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

% inflows by reason for unemployment

w_o = reshape(sum(est_TU,1),T,5); 
w_uh = reshape(sum(est_UH,1),T,5); 
w_ul = reshape(sum(est_UL,1),T,5); 

in = sum(w_o,2);

for n=1:5
reason_share(:,n) = w_o(:,n)./sum(w_o,2);
end

tl_u = [w_uh(:,1), w_ul(:,1)];
os_u = [w_uh(:,2), w_ul(:,2)];
qt_u = [w_uh(:,3), w_ul(:,3)];
re_u = [w_uh(:,4), w_ul(:,4)];
ne_u = [w_uh(:,5), w_ul(:,5)];

for n=1:2
tl_share(:,n) = tl_u(:,n)./sum(tl_u,2);
os_share(:,n) = os_u(:,n)./sum(os_u,2);
qt_share(:,n) = qt_u(:,n)./sum(qt_u,2);
re_share(:,n) = re_u(:,n)./sum(re_u,2);
ne_share(:,n) = ne_u(:,n)./sum(ne_u,2);
end

%% count of unemployed individuals by reason for unemployment

% load temp_all
% load temp_dur

% 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
%beg = 1
%beg = find(tt == (2000 + 1/13)); % Gr%eat Recession

% (1) counterfactual duration with the inflows fixed at December 2007
creason_share = reason_share;

ctl_share = tl_share;
cos_share = os_share;
cqt_share = qt_share;
cre_share = re_share;
cne_share = ne_share;

% date December 2007 : line number beg 
% cin(beg:end,1) = in(beg,1)*ones(409-beg+1,size(cin,2));

creason_share(beg:end,:) = ones(T-beg+1,1)*reason_share(beg,:);
ctl_share(beg:end,:) = ones(T-beg+1,1)*tl_share(beg,:);
cos_share(beg:end,:) = ones(T-beg+1,1)*os_share(beg,:);
cqt_share(beg:end,:) = ones(T-beg+1,1)*qt_share(beg,:);
cre_share(beg:end,:) = ones(T-beg+1,1)*re_share(beg,:);
cne_share(beg:end,:) = ones(T-beg+1,1)*ne_share(beg,:);
cdur_PH = dur_PH;
cdur_PL = dur_PL;

for n=1:5
cdur_PH(beg:end,n) = dur_PH(beg,n)*ones(1,T-beg+1);
cdur_PL(beg:end,n) = dur_PL(beg,n)*ones(1,T-beg+1);
end

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

% data, no changes, observed share, unobserved share, type H, and type L
hist_dur = [D(beg:end,1), c0_D(beg:end,1), co_D(beg:end,1), cu_D(beg:end,1), ...
     cl_D(beg:end,1), ch_D(beg:end,1)];
    
plot(hist_dur)

hist_LR1 = [D, c0_D, co_D, cu_D, cl_D, ch_D]
hist_LR2 = [D-c0_D, co_D-c0_D, cu_D-co_D, cl_D-cu_D, ch_D-cl_D]

%% plot hist_LR2

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_LR2(plot_beg:plot_end,:);

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') 
%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,'/figures/hist_dur','epsc');

%clearvars -except hist_dur
%open hist_dur

% compute counter-factural mean duration
A1_14_counterfactual
