
clear all
pwd

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

load tl_3prm_rec_e24_final;

U1 = U(:,1);
TU = sum(U,2);
T = size(TU,1);
for t=2:T
p(t,1) = (TU(t,1)-U1(t,1))/TU(t-1,1);
end

DU1 = U1;
T1 = TU;
p1 = p;

keep DU1 T1 p1

load os_3prm_rec_e24_final;

U1 = U(:,1);
TU = sum(U,2);
T = size(TU,1);
for t=2:T
p(t,1) = (TU(t,1)-U1(t,1))/TU(t-1,1);
end

DU2 = U1;
T2 = TU;
p2 = p;

keep DU1 T1 p1 ...
     DU2 T2 p2
 
load qt_3prm_rec_e24_final;

U1 = U(:,1);
TU = sum(U,2);
T = size(TU,1);
for t=2:T
p(t,1) = (TU(t,1)-U1(t,1))/TU(t-1,1);
end

DU3 = U1;
T3 = TU;
p3 = p;

keep DU1 T1 p1 ...
     DU2 T2 p2 ...
     DU3 T3 p3 

 
load re_3prm_rec_e24_final;

U1 = U(:,1);
TU = sum(U,2);
T = size(TU,1);
for t=2:T
p(t,1) = (TU(t,1)-U1(t,1))/TU(t-1,1);
end

DU4 = U1;
T4 = TU;
p4 = p;

keep DU1 T1 p1 ...
     DU2 T2 p2 ...
     DU3 T3 p3 ... 
     DU4 T4 p4 
 
load ne_3prm_rec_e24_final;

U1 = U(:,1);
TU = sum(U,2);
T = size(TU,1);
for t=2:T
p(t,1) = (TU(t,1)-U1(t,1))/TU(t-1,1);
end

DU5 = U1;
T5 = TU;
p5 = p;

keep DU1 T1 p1 ...
     DU2 T2 p2 ...
     DU3 T3 p3 ... 
     DU4 T4 p4 ...
     DU5 T5 p5 

%% ongoing distribution calculation

w(:,1) = DU1;
w(:,2) = DU2;
w(:,3) = DU3;
w(:,4) = DU4;
w(:,5) = DU5;

p(:,1) = p1;
p(:,2) = p2;
p(:,3) = p3;
p(:,4) = p4;
p(:,5) = p5;

T = size(p,1);

for n=1:5
for j=1:48
for i=j:T
U(i,j,n) = w(i-j+1,n)*prod(p(i-j+2:i,n),1);
end
end
end

U(:,49,:) = sum(U(:,1:48,:),2);
TU = sum(U,3);
% distribution

for nn=1:5
for jj=1:48
for t=1:T
dis(t,jj,nn) = U(t,jj,nn)/U(t,49,nn);
end
end
end


for nn=1:5
for jj=1:48
for t=1:T
Tdis(t,jj,nn) = U(t,jj,nn)/TU(t,49);
end
end
end

% moment calculation

for n=1:5
for t=1:T
for jj=1:48
D_dis(t,jj,n) = jj*dis(t,jj,n);
D2_dis(t,jj,n) = (jj^2)*dis(t,jj,n);
end
end
end

MD_dis = sum(D_dis,2);
MD2_dis = sum(D2_dis,2);

for n=1:5
for t=1:T
for jj=1:48
TD_dis(t,jj,n) = jj*Tdis(t,jj,n);
TD2_dis(t,jj,n) = (jj^2)*Tdis(t,jj,n);
end
end
end

TMD_dis = sum(sum(TD_dis,3),2);
TMD2_dis = sum(sum(TD2_dis,3),2);

for n=1:5
for t=1:T
mean_ig(t,n) = MD_dis(t,1,n);
var_ig(t,n) = MD2_dis(t,1,n)- MD_dis(t,1,n)^2;
end
end

for t=1:T
mean_group(t,1) = TMD_dis(t,1);
var_group(t,1) = TMD2_dis(t,1)- TMD_dis(t,1)^2;
end

mean_group = mean_group(1:T,1);
var_group = var_group(1:T,1);

keep mean_group var_group

%% MEMO: section 1 of A1 is calculated from this %
%clear all;
%load data_agg
% 1948:m01 - 2016:m02 => p and w

[nun] = xlsread('udur_fred');
data = nun(:,2)/(52/12); % mean duration of unemployment
wdata = nun(:,2);

w = nun(:,3);
U = nun(:,4);

T =size(U,1);
for t=1:T-1
p(t+1,1) = (U(t+1,1)-w(t+1,1))./U(t,1);
end

U(1,1) = w(1,1);

for j=1:48
for i=j:T
U(i,j) = w(i-j+1,1)*prod(p(i-j+2:i,1),1);
end
end
U(:,49) = sum(U(:,1:48),2);


for j=1:48
for i=j:T
U(i,j) = w(i-j+1,1)*prod(p(i-j+2:i,1),1);
end
end

U(:,49,:) = sum(U(:,1:48,:),2);

% distribution

for jj=1:48
for t=j:T
dis(t,jj,1) = U(t,jj)/U(t,49);
end
end

% moment calculation

for t=1:T
for jj=1:48
D_dis(t,jj) = jj*dis(t,jj);
D2_dis(t,jj) = (jj^2)*dis(t,jj);
end
end

MD_dis = sum(D_dis,2);
MD2_dis = sum(D2_dis,2);

for t=1:T
mean(t,1) = MD_dis(t,1);
var(t,1) = sqrt(MD2_dis(t,1)- MD_dis(t,1)^2);
end

mean = mean(1:T,:);
var = var(1:T,:);

%% data: distribution of unemployment duration

[tmp0] = xlsread('udur_week_dist');

dis_dur = tmp0(1:end,1:end);
tt = size(dis_dur,1);
J = size(dis_dur,2);

for t=1:tt
    for j=1:J
tmp(t,j) = dis_dur(t,j)*(j-wdata(528-tt+1,1))^2;
    end    
end

data_stdev = sqrt(sum(tmp,2))/(52/12);

[tmp1] = xlsread('udur_stdev_1978_2013');
tmp1 = tmp1(24:end,2)

data_stdev = [tmp1;data_stdev];

%% Figure 2
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

beg = 48;
t = tt(beg:T-1,1); 

g1 = figure;
hold on
set(gca,'fontsize',15)
plot(t,mean(beg:T-1,1),'-','LineWidth',4,'Color',[1 0 0])
plot(t,mean_group(beg:T-1,1),'--','LineWidth',4,'Color',[1 0 1])
plot(t,data(beg:T-1,1),'-','LineWidth',2,'Color',[0 0 1])
hleg1 = legend('No heterogeneity','Observed heterogeneity','Data');
set(hleg1,'Location', 'Northwest') 
colorstr=[159 182 205]/256;
axis([1979.9,2019.8,0,12])
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,'/figures/mean_agg','epsc');

g2 = figure;
hold on
set(gca,'fontsize',15)
plot(t,var(beg:T-1,1),'-','LineWidth',4,'Color',[1 0 0])
plot(t,sqrt(var_group(beg:T-1,1)),'--','LineWidth',4,'Color',[1 0 1])
plot(t,data_stdev(1:end-1,1),'-','LineWidth',2,'Color',[0 0 1])
hleg1 = legend('No heterogeneity','Observed heterogeneity','Data');
set(hleg1,'Location', 'Northwest') 
axis([1979.9,2019.8,0,12])
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(g2,'/figures/std_agg','epsc');