
% 1980 recession: hh=1
% 1981 recession: hh=16
% 1990 recession: hh=126
% 2001 recessoin: hh=253
% 2007 recession: hh=336 (2007m11)

T=498;
hh=336;
s=48; b=hh+s; e=b+s; 

nt=48; qnt=48*4; 

for nn=1:nt
    tmp = 4*(nn-1);
    ssxi(tmp+1,nn:nn+T-1) = sxi(1,:);
    ssxi(tmp+2,nn:nn+T-1) = sxi(2,:);
    ssxi(tmp+3,nn:nn+T-1) = sxi(3,:);
    ssxi(tmp+4,nn:nn+T-1) = sxi(4,:);
end

 dssxi = ssxi(:,nt:T);
 dsxi = sxi(:,nt:T);
 dxi = xi0(:,nt:T);
 TT = T-nt+1;
 nn = 48;
 Xi=dssxi;
 
F = zeros(qnt,qnt);
F(1:4,1:4)=eye(4);
for q=1:nt-1
F(1+4*q:4+4*q,1+4*(q-1):4+4*(q-1))=eye(4);
end
 
%% baseline structure specification %%

W1 = zeros(nt,TT);
W2 = zeros(nt,TT);
X1 = zeros(nt,TT);
X2 = zeros(nt,TT);
MW1 = ones(nt,nt,TT);
MW2 = ones(nt,nt,TT);
D_MW1 = zeros(nt,nt,TT);
D_MW2 = zeros(nt,nt,TT);
MP1 = zeros(nt,nt,TT);
MP2 = zeros(nt,nt,TT);
D_MP1 = zeros(nt,nt,TT);
D_MP2 = zeros(nt,nt,TT);
MC1 = zeros(nt,nt,TT);
MC2 = zeros(nt,nt,TT);
M_JH1 = zeros(nt,qnt,TT);
UL=zeros(nt,TT);
UH=zeros(nt,TT);
TU=zeros(nt,TT);
yhat = zeros(5,TT);
fJH1 = zeros(5,qnt,TT);
fJH = zeros(qnt,5,TT);
CFJH=zeros(4,nt);
FJH = zeros(4,nt,TT);
FJH1=zeros(TT,qnt);
fxi0 = zeros(qnt,TT);

for b=hh:hh
    e=b+s;
for kk = 1:s+1
fxi0(:,b+kk-1) = (F^(kk-1))*Xi(:,b);
end

for t=b:e
      
 for i=1:nt
 W1(i,t) = fxi0(1+4*(i-1),t);
 W2(i,t) = fxi0(2+4*(i-1),t);
 X1(i,t) = fxi0(3+4*(i-1),t);
 X2(i,t) = fxi0(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
 
 D_MW1(1,1,t)= 1;
 D_MW2(1,1,t)= 1;
 
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
  
 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

b_fcst(b,:) = sum(TU(:,b:e),1); % baseline U forecast
b_fcst1(b,:) = sum(TU(1,b:e),1); % baseline U1 forecast
b_fcst23(b,:) = sum(TU(2:3,b:e),1); % baseline U23 forecast
b_fcst46(b,:) = sum(TU(4:6,b:e),1); % baseline U46 forecast
b_fcst712(b,:) = sum(TU(7:12,b:e),1); % baseline U712 forecast
b_fcst13(b,:) = sum(TU(13:48,b:e),1); % baseline U13+ forecast

%% nonlinear forecasting %%
% baseline structure specification %%

W1 = zeros(nt,TT);
W2 = zeros(nt,TT);
X1 = zeros(nt,TT);
X2 = zeros(nt,TT);
MW1 = ones(nt,nt,TT);
MW2 = ones(nt,nt,TT);
D_MW1 = zeros(nt,nt,TT);
D_MW2 = zeros(nt,nt,TT);
MP1 = zeros(nt,nt,TT);
MP2 = zeros(nt,nt,TT);
D_MP1 = zeros(nt,nt,TT);
D_MP2 = zeros(nt,nt,TT);
MC1 = zeros(nt,nt,TT);
MC2 = zeros(nt,nt,TT);
M_JH1 = zeros(nt,qnt,TT);
UL=zeros(nt,TT);
UH=zeros(nt,TT);
TU=zeros(nt,TT);
yhat = zeros(5,TT);
fJH1 = zeros(5,qnt,TT);
fJH = zeros(qnt,5,TT);
CFJH=zeros(4,nt);
FJH = zeros(4,nt,TT);
FJH1=zeros(TT,qnt);

fxi = fxi0;
fxi_lx = zeros(qnt,TT);
fxi_lx1 = zeros(qnt,TT);
fxi_lx2 = zeros(qnt,TT);
fxi_lx3 = zeros(qnt,TT);
fxi_lx4 = zeros(qnt,TT);

ups_lx = zeros(qnt,TT);
ups_lx1 = zeros(qnt,TT);
ups_lx2 = zeros(qnt,TT);
ups_lx3 = zeros(qnt,TT);
ups_lx4 = zeros(qnt,TT);

%% smoothed estimates from b to e %%

for t=b+1:e
    ss = t-b+1;
ups(1:4,ss) = Xi(1:4,t)-Xi(1:4,b);
end

ups_lx = zeros(qnt,s+1);
ups_lx1 = zeros(qnt,s+1);
ups_lx2 = zeros(qnt,s+1);
ups_lx3 = zeros(qnt,s+1);
ups_lx4 = zeros(qnt,s+1);

for nn=1:nt
    tmp = 4*(nn-1);
    ups_lx(tmp+1,nn:nn+s) = ups(1,:);
    ups_lx(tmp+2,nn:nn+s) = ups(2,:);
    ups_lx(tmp+3,nn:nn+s) = ups(3,:);
    ups_lx(tmp+4,nn:nn+s) = ups(4,:);
    ups_lx1(tmp+1,nn:nn+s) = ups(1,:);
    ups_lx2(tmp+2,nn:nn+s) = ups(2,:);
    ups_lx3(tmp+3,nn:nn+s) = ups(3,:);
    ups_lx4(tmp+4,nn:nn+s) = ups(4,:);
end

ups_lx = ups_lx(1:qnt,1:s+1);
ups_lx1 = ups_lx1(1:qnt,1:s+1);
ups_lx2 = ups_lx2(1:qnt,1:s+1);
ups_lx3 = ups_lx3(1:qnt,1:s+1);
ups_lx4 = ups_lx4(1:qnt,1:s+1);

for t=b+1:e %TT

 for i=1:nt
 W1(i,t) = fxi0(1+4*(i-1),t);
 W2(i,t) = fxi0(2+4*(i-1),t);
 X1(i,t) = fxi0(3+4*(i-1),t);
 X2(i,t) = fxi0(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
 
 D_MW1(1,1,t)= 1;
 D_MW2(1,1,t)= 1;

 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
 
 MC1(1,:,t)=zeros(1,nt);
 MC2(1,:,t)=zeros(1,nt);

 for i=2:nt
 for j=1:nt 
     if tau(i,j)>0
 MC1(i,j,t) = -exp(X1(j,t)+D(tau(i,j),1));
 MC2(i,j,t) = -exp(X2(j,t)+D(tau(i,j),1));
     else
 MC1(i,j,t) = 0;
 MC2(i,j,t) = 0;
 end
 end
 end
 
 MP1(1,:,t)=zeros(1,nt);
 MP2(1,:,t)=zeros(1,nt);

 for i=2:nt
 MP1(i,1:i-1,t)=W1(i,t)*D_MW1(i,i,t);
 MP2(i,1:i-1,t)=W2(i,t)*D_MW2(i,i,t);
 end

 for i=1:nt
 for j=1:nt 
 D_MP1(i,j,t)= MP1(i,j,t).*MC1(i,j,t);
 D_MP2(i,j,t)= MP2(i,j,t).*MC2(i,j,t);
 end
 end

for i=1:nt
M_JH1(:,1+4*(i-1),t)=D_MW1(:,i,t);
M_JH1(:,2+4*(i-1),t)=D_MW2(:,i,t);
M_JH1(:,3+4*(i-1),t)=D_MP1(:,i,t);
M_JH1(:,4+4*(i-1),t)=D_MP2(:,i,t);
end

fJH1(1,:,t)= M_JH1(1,:,t);
fJH1(2,:,t)= M_JH1(2,:,t)+M_JH1(3,:,t);
fJH1(3,:,t)= M_JH1(4,:,t)+M_JH1(5,:,t)+M_JH1(6,:,t);
fJH1(4,:,t)= sum(M_JH1(7:12,:,t),1);
fJH1(5,:,t)= sum(M_JH1(13:48,:,t),1);

FJH1(t,:) = ones(1,5)*fJH1(:,:,t);
FJH(:,:,t) = reshape(FJH1(t,:)',4,nt);

n=t-b;
if n<nt
FJH(:,n+1:nt,t) = zeros(4,nt-n);
end
MFJH(:,n+1,b)=reshape(FJH(:,:,t),192,1);
end
FJH = FJH1';

ups_m(:,:,b)= ups_lx.*MFJH(:,:,b);
ups_m_lx1(:,:,b)= ups_lx1.*MFJH(:,:,b);
ups_m_lx2(:,:,b)= ups_lx2.*MFJH(:,:,b);
ups_m_lx3(:,:,b)= ups_lx3.*MFJH(:,:,b);
ups_m_lx4(:,:,b)= ups_lx4.*MFJH(:,:,b);


seps_m(b,:) = sum(ups_m(:,:,b),1);
seps_m_lx1(b,:) = sum(ups_m_lx1(:,:,b),1);
seps_m_lx2(b,:) = sum(ups_m_lx2(:,:,b),1);
seps_m_lx3(b,:) = sum(ups_m_lx3(:,:,b),1);
seps_m_lx4(b,:) = sum(ups_m_lx4(:,:,b),1);

lu_fcst_all(b,:) = b_fcst(b,:)+seps_m(b,:); %green
lu_fcst_lx1(b,:) = b_fcst(b,:)+seps_m_lx1(b,:); %red
lu_fcst_lx2(b,:) = b_fcst(b,:)+seps_m_lx2(b,:); %cyan
lu_fcst_lx3(b,:) = b_fcst(b,:)+seps_m_lx3(b,:); %purple
lu_fcst_lx4(b,:) = b_fcst(b,:)+seps_m_lx4(b,:); %mustard

% U1 %
seps_m1(b,:) = sum(ups_m(1:2,:,b),1);
seps_m1_lx1(b,:) = sum(ups_m_lx1(1:2,:,b),1);
seps_m1_lx2(b,:) = sum(ups_m_lx2(1:2,:,b),1);
seps_m1_lx3(b,:) = sum(ups_m_lx3(1:2,:,b),1);
seps_m1_lx4(b,:) = sum(ups_m_lx4(1:2,:,b),1);

lu1_fcst_all(b,:) = b_fcst1(b,:)+seps_m1(b,:); %green
lu1_fcst_lx1(b,:) = b_fcst1(b,:)+seps_m1_lx1(b,:); %red
lu1_fcst_lx2(b,:) = b_fcst1(b,:)+seps_m1_lx2(b,:); %cyan
lu1_fcst_lx3(b,:) = b_fcst1(b,:)+seps_m1_lx3(b,:); %purple
lu1_fcst_lx4(b,:) = b_fcst1(b,:)+seps_m1_lx4(b,:); %mustard

lu1_fcst_ex1(b,:) = b_fcst1(b,:)+seps_m1_lx2(b,:)+seps_m1_lx3(b,:)+seps_m1_lx4(b,:); %red
lu1_fcst_ex2(b,:) = b_fcst1(b,:)+seps_m1_lx1(b,:)+seps_m1_lx3(b,:)+seps_m1_lx4(b,:); %cyan
lu1_fcst_ex3(b,:) = b_fcst1(b,:)+seps_m1_lx1(b,:)+seps_m1_lx2(b,:)+seps_m1_lx4(b,:); %purple
lu1_fcst_ex4(b,:) = b_fcst1(b,:)+seps_m1_lx1(b,:)+seps_m1_lx2(b,:)+seps_m1_lx3(b,:); %mustard

% U23 %
seps_m23(b,:) = sum(ups_m(3:10,:,b),1);
seps_m23_lx1(b,:) = sum(ups_m_lx1(3:10,:,b),1);
seps_m23_lx2(b,:) = sum(ups_m_lx2(3:10,:,b),1);
seps_m23_lx3(b,:) = sum(ups_m_lx3(3:10,:,b),1);
seps_m23_lx4(b,:) = sum(ups_m_lx4(3:10,:,b),1);

lu23_fcst_all(b,:) = b_fcst23(b,:)+seps_m23(b,:); %green
lu23_fcst_lx1(b,:) = b_fcst23(b,:)+seps_m23_lx1(b,:); %red
lu23_fcst_lx2(b,:) = b_fcst23(b,:)+seps_m23_lx2(b,:); %cyan
lu23_fcst_lx3(b,:) = b_fcst23(b,:)+seps_m23_lx3(b,:); %purple
lu23_fcst_lx4(b,:) = b_fcst23(b,:)+seps_m23_lx4(b,:); %mustard

lu23_fcst_ex1(b,:) = b_fcst23(b,:)+seps_m23_lx2(b,:)+seps_m23_lx3(b,:)+seps_m23_lx4(b,:); %red
lu23_fcst_ex2(b,:) = b_fcst23(b,:)+seps_m23_lx1(b,:)+seps_m23_lx3(b,:)+seps_m23_lx4(b,:); %cyan
lu23_fcst_ex3(b,:) = b_fcst23(b,:)+seps_m23_lx1(b,:)+seps_m23_lx2(b,:)+seps_m23_lx4(b,:); %purple
lu23_fcst_ex4(b,:) = b_fcst23(b,:)+seps_m23_lx1(b,:)+seps_m23_lx2(b,:)+seps_m23_lx3(b,:); %mustard

% U46 %
seps_m46(b,:) = sum(ups_m(11:22,:,b),1);
seps_m46_lx1(b,:) = sum(ups_m_lx1(11:22,:,b),1);
seps_m46_lx2(b,:) = sum(ups_m_lx2(11:22,:,b),1);
seps_m46_lx3(b,:) = sum(ups_m_lx3(11:22,:,b),1);
seps_m46_lx4(b,:) = sum(ups_m_lx4(11:22,:,b),1);

lu46_fcst_all(b,:) = b_fcst46(b,:)+seps_m46(b,:); %green
lu46_fcst_lx1(b,:) = b_fcst46(b,:)+seps_m46_lx1(b,:); %red
lu46_fcst_lx2(b,:) = b_fcst46(b,:)+seps_m46_lx2(b,:); %cyan
lu46_fcst_lx3(b,:) = b_fcst46(b,:)+seps_m46_lx3(b,:); %purple
lu46_fcst_lx4(b,:) = b_fcst46(b,:)+seps_m46_lx4(b,:); %mustard

lu46_fcst_ex1(b,:) = b_fcst46(b,:)+seps_m46_lx2(b,:)+seps_m46_lx3(b,:)+seps_m46_lx4(b,:); %red
lu46_fcst_ex2(b,:) = b_fcst46(b,:)+seps_m46_lx1(b,:)+seps_m46_lx3(b,:)+seps_m46_lx4(b,:); %cyan
lu46_fcst_ex3(b,:) = b_fcst46(b,:)+seps_m46_lx1(b,:)+seps_m46_lx2(b,:)+seps_m46_lx4(b,:); %purple
lu46_fcst_ex4(b,:) = b_fcst46(b,:)+seps_m46_lx1(b,:)+seps_m46_lx2(b,:)+seps_m46_lx3(b,:); %mustard

% U712 %
seps_m712(b,:) = sum(ups_m(23:46,:,b),1);
seps_m712_lx1(b,:) = sum(ups_m_lx1(23:46,:,b),1);
seps_m712_lx2(b,:) = sum(ups_m_lx2(23:46,:,b),1);
seps_m712_lx3(b,:) = sum(ups_m_lx3(23:46,:,b),1);
seps_m712_lx4(b,:) = sum(ups_m_lx4(23:46,:,b),1);

lu712_fcst_all(b,:) = b_fcst712(b,:)+seps_m712(b,:); %green
lu712_fcst_lx1(b,:) = b_fcst712(b,:)+seps_m712_lx1(b,:); %red
lu712_fcst_lx2(b,:) = b_fcst712(b,:)+seps_m712_lx2(b,:); %cyan
lu712_fcst_lx3(b,:) = b_fcst712(b,:)+seps_m712_lx3(b,:); %purple
lu712_fcst_lx4(b,:) = b_fcst712(b,:)+seps_m712_lx4(b,:); %mustard

lu712_fcst_ex1(b,:) = b_fcst712(b,:)+seps_m712_lx2(b,:)+seps_m712_lx3(b,:)+seps_m712_lx4(b,:); %red
lu712_fcst_ex2(b,:) = b_fcst712(b,:)+seps_m712_lx1(b,:)+seps_m712_lx3(b,:)+seps_m712_lx4(b,:); %cyan
lu712_fcst_ex3(b,:) = b_fcst712(b,:)+seps_m712_lx1(b,:)+seps_m712_lx2(b,:)+seps_m712_lx4(b,:); %purple
lu712_fcst_ex4(b,:) = b_fcst712(b,:)+seps_m712_lx1(b,:)+seps_m712_lx2(b,:)+seps_m712_lx3(b,:); %mustard

% U13 %
seps_m13(b,:) = sum(ups_m(47:192,:,b),1);
seps_m13_lx1(b,:) = sum(ups_m_lx1(47:192,:,b),1);
seps_m13_lx2(b,:) = sum(ups_m_lx2(47:192,:,b),1);
seps_m13_lx3(b,:) = sum(ups_m_lx3(47:192,:,b),1);
seps_m13_lx4(b,:) = sum(ups_m_lx4(47:192,:,b),1);

lu13_fcst_all(b,:) = b_fcst13(b,:)+seps_m13(b,:); %green
lu13_fcst_lx1(b,:) = b_fcst13(b,:)+seps_m13_lx1(b,:); %red
lu13_fcst_lx2(b,:) = b_fcst13(b,:)+seps_m13_lx2(b,:); %cyan
lu13_fcst_lx3(b,:) = b_fcst13(b,:)+seps_m13_lx3(b,:); %purple
lu13_fcst_lx4(b,:) = b_fcst13(b,:)+seps_m13_lx4(b,:); %mustard

lu13_fcst_ex1(b,:) = b_fcst13(b,:)+seps_m13_lx2(b,:)+seps_m13_lx3(b,:)+seps_m13_lx4(b,:); %red
lu13_fcst_ex2(b,:) = b_fcst13(b,:)+seps_m13_lx1(b,:)+seps_m13_lx3(b,:)+seps_m13_lx4(b,:); %cyan
lu13_fcst_ex3(b,:) = b_fcst13(b,:)+seps_m13_lx1(b,:)+seps_m13_lx2(b,:)+seps_m13_lx4(b,:); %purple
lu13_fcst_ex4(b,:) = b_fcst13(b,:)+seps_m13_lx1(b,:)+seps_m13_lx2(b,:)+seps_m13_lx3(b,:); %mustard

PLOT_L = [b_fcst(b,:)',lu_fcst_all(b,:)',...
    lu_fcst_lx1(b,:)',lu_fcst_lx2(b,:)',lu_fcst_lx3(b,:)',lu_fcst_lx4(b,:)'];

PLOT_L1 = [b_fcst1(b,:)',lu1_fcst_all(b,:)',...
    lu1_fcst_lx1(b,:)',lu1_fcst_lx2(b,:)',lu1_fcst_lx3(b,:)',lu1_fcst_lx4(b,:)'];

PLOT_X1 = [b_fcst1(b,:)',lu1_fcst_all(b,:)',...
    lu1_fcst_ex1(b,:)',lu1_fcst_ex2(b,:)',lu1_fcst_ex3(b,:)',lu1_fcst_ex4(b,:)'];

PLOT_L23 = [b_fcst23(b,:)',lu23_fcst_all(b,:)',...
    lu23_fcst_lx1(b,:)',lu23_fcst_lx2(b,:)',lu23_fcst_lx3(b,:)',lu23_fcst_lx4(b,:)'];

PLOT_X23 = [b_fcst23(b,:)',lu23_fcst_all(b,:)',...
    lu23_fcst_ex1(b,:)',lu23_fcst_ex2(b,:)',lu23_fcst_ex3(b,:)',lu23_fcst_ex4(b,:)'];

PLOT_L46 = [b_fcst46(b,:)',lu46_fcst_all(b,:)',...
    lu46_fcst_lx1(b,:)',lu46_fcst_lx2(b,:)',lu46_fcst_lx3(b,:)',lu46_fcst_lx4(b,:)'];

PLOT_X46 = [b_fcst46(b,:)',lu46_fcst_all(b,:)',...
    lu46_fcst_ex1(b,:)',lu46_fcst_ex2(b,:)',lu46_fcst_ex3(b,:)',lu46_fcst_ex4(b,:)'];

PLOT_L712 = [b_fcst712(b,:)',lu712_fcst_all(b,:)',...
    lu712_fcst_lx1(b,:)',lu712_fcst_lx2(b,:)',lu712_fcst_lx3(b,:)',lu712_fcst_lx4(b,:)'];

PLOT_X712 = [b_fcst712(b,:)',lu712_fcst_all(b,:)',...
    lu712_fcst_ex1(b,:)',lu712_fcst_ex2(b,:)',lu712_fcst_ex3(b,:)',lu712_fcst_ex4(b,:)'];

PLOT_L13 = [b_fcst13(b,:)',lu13_fcst_all(b,:)',...
    lu13_fcst_lx1(b,:)',lu13_fcst_lx2(b,:)',lu13_fcst_lx3(b,:)',lu13_fcst_lx4(b,:)'];

PLOT_X13 = [b_fcst13(b,:)',lu13_fcst_all(b,:)',...
    lu13_fcst_ex1(b,:)',lu13_fcst_ex2(b,:)',lu13_fcst_ex3(b,:)',lu13_fcst_ex4(b,:)'];

end

APLOT = [PLOT_L,PLOT_L1,PLOT_L23,PLOT_L46,PLOT_L712,PLOT_L13,PLOT_X1,PLOT_X23,PLOT_X46,PLOT_X712,PLOT_X13];

%%
U1980r = PLOT_L; %h=1 4param
U1981r = PLOT_L; %h=1 4param
U1990r = PLOT_L; %h=1 4param
U2001r = PLOT_L; %h=1 4param
U2007r = PLOT_L; %h=1 4param

%% JDH exercise %

J1 = [Y(hh+47:hh+47+23,1),PLOT_L1(1:24,1:6)];
J23 = [Y(hh+47:hh+47+23,2),PLOT_L23(1:24,1:6)];
J46 = [Y(hh+47:hh+47+23,3),PLOT_L46(1:24,1:6)];
J712 = [Y(hh+47:hh+47+23,4),PLOT_L712(1:24,1:6)];
J13 = [Y(hh+47:hh+47+23,5),PLOT_L13(1:24,1:6)];

JJ = [J1,J23,J46,J712,J13]; 

for i=1:24
JD(i,:)=JJ(i,:)-JJ(1,:)
end
