function total_Y = gdd_counterfactual(D,sxi)
% notation : handbook of econometrics "state space model" by Hamilton.
global T;global Y;global P0; global P1; global F;
global Q;global R;

 nn=4; nr=5; nt=48;qnt=192;
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 t=1:T
 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

 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

 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

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

yhat(1,t)=TU(1,t);
yhat(2,t)=TU(2,t)+TU(3,t);
yhat(3,t)=TU(4,t)+TU(5,t)+TU(6,t);
yhat(4,t)=sum(TU(7:12,t),1);
yhat(5,t)=sum(TU(13:nt,t),1);

total_Y = sum(yhat,1);
end
