function [hist_Ulevel] = histdcmp(hh,sxi,F,D)

global T;global F;
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);
 TT = T-nt+1;
 nn = 48;
 Xi=dssxi;
 
 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

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

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

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 

 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

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

end