function f = ekf_chv_rec_jac(x)
% notation : handbook of econometrics "state space model" by Hamilton.
global T;global Y;global xi0; global xi1; global P0; global P1; global F;
global d11;global d12;global d13;global Q;global R;
global r11;global r12;global r13;

 nn=4; nr=5; nt=48;qnt=192;
load indfile;
IND=ind2;
% parameters for the dynamics of latent variables
% variance for epsilon: innovation of the latent variables

  q1=x(1)^2;
  q2=x(2)^2;
  q3=x(3)^2;
  q4=x(4)^2;
  
F=zeros(4*nt,4*nt);
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

qq = [q1,0,0,0;...
      0,q2,0,0;...
      0,0,q3,0;...
      0,0,0,q4];
  
Q=zeros(4*nt,4*nt);
Q(1:4,1:4)=qq;

d11= x(5); %duration 2-3 m
d12= x(6); %duration 4-6 m
d13= x(7); %duration 7-12 m

r11= x(13); %duration 2-3 m
r12= x(14); %duration 4-6 m
r13= x(15); %duration 7-12 m

% variance of the measurement errors

 R(1,1)=x(8)^2;
 R(2,2)=x(9)^2;
 R(3,3)=x(10)^2;
 R(4,4)=x(11)^2;
 R(5,5)=x(12)^2;

% size specification: Matrix for filter 
c=10;
sigma2 = 0.1;
I4 = c*eye(4,4);
S = zeros(192,192,48);
S(1:192,1:192,1)=kron(ones(48,48),I4);
for j=2:48
    S(4*(j-1)+1:192,4*(j-1)+1:192,j)=kron(ones(48-(j-1),48-(j-1)),(sigma2)*eye(4,4));
end

P0(:,:,1)=sum(S,3);

ll=ones(T,1);
PP=ones(T,2);
W=ones(T,2);
JH1 = zeros(nr,4*nt,T);
JH = zeros(4*nt,nr,T);
VY = zeros(nr,nr,T);
yhat = zeros(nr,T);
xi1 = zeros(4*nt,T);
P1 = zeros(4*nt,4*nt,T);

% Matrix of GDD by duration (column)
BD = zeros(nt,1);

BD(1,1)=0;
BD(2,1)=d11*(1/48) + d12*[2*(1/48)^2-1] + d13*[4*(1/48)^3-3*(1/48)];
BD(3,1)=d11*(2/48) + d12*[2*(2/48)^2-1] + d13*[4*(2/48)^3-3*(2/48)];
BD(4,1)=d11*(3/48) + d12*[2*(3/48)^2-1] + d13*[4*(3/48)^3-3*(3/48)];
BD(5,1)=d11*(4/48) + d12*[2*(4/48)^2-1] + d13*[4*(4/48)^3-3*(4/48)];
BD(6,1)=d11*(5/48) + d12*[2*(5/48)^2-1] + d13*[4*(5/48)^3-3*(5/48)];
BD(7,1)=d11*(6/48) + d12*[2*(6/48)^2-1] + d13*[4*(6/48)^3-3*(6/48)];
BD(8,1)=d11*(7/48) + d12*[2*(7/48)^2-1] + d13*[4*(7/48)^3-3*(7/48)];
BD(9,1)=d11*(8/48) + d12*[2*(8/48)^2-1] + d13*[4*(8/48)^3-3*(8/48)];
BD(10,1)=d11*(9/48) + d12*[2*(9/48)^2-1] + d13*[4*(9/48)^3-3*(9/48)];
BD(11,1)=d11*(10/48) + d12*[2*(10/48)^2-1] + d13*[4*(10/48)^3-3*(10/48)];

for i=12:nt
BD(i,1)=d11*((i-1)/48) + d12*[2*((i-1)/48)^2-1] + d13*[4*((i-1)/48)^3-3*((i-1)/48)];
end %duration 13+ m

RD(1,1)=0;
RD(2,1)=r11*(1/48) + r12*[2*(1/48)^2-1] + r13*[4*(1/48)^3-3*(1/48)];
RD(3,1)=r11*(2/48) + r12*[2*(2/48)^2-1] + r13*[4*(2/48)^3-3*(2/48)];
RD(4,1)=r11*(3/48) + r12*[2*(3/48)^2-1] + r13*[4*(3/48)^3-3*(3/48)];
RD(5,1)=r11*(4/48) + r12*[2*(4/48)^2-1] + r13*[4*(4/48)^3-3*(4/48)];
RD(6,1)=r11*(5/48) + r12*[2*(5/48)^2-1] + r13*[4*(5/48)^3-3*(5/48)];
RD(7,1)=r11*(6/48) + r12*[2*(6/48)^2-1] + r13*[4*(6/48)^3-3*(6/48)];
RD(8,1)=r11*(7/48) + r12*[2*(7/48)^2-1] + r13*[4*(7/48)^3-3*(7/48)];
RD(9,1)=r11*(8/48) + r12*[2*(8/48)^2-1] + r13*[4*(8/48)^3-3*(8/48)];
RD(10,1)=r11*(9/48) + r12*[2*(9/48)^2-1] + r13*[4*(9/48)^3-3*(9/48)];
RD(11,1)=r11*(10/48) + r12*[2*(10/48)^2-1] + r13*[4*(10/48)^3-3*(10/48)];

for i=12:nt
RD(i,1)=r11*((i-1)/48) + r12*[2*((i-1)/48)^2-1] + r13*[4*((i-1)/48)^3-3*((i-1)/48)];
end %duration 13+ m

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

% size specification: Derivatives for EKF
 MW1(:,:,:) = ones(nt,nt,T);
 MW2(:,:,:) = ones(nt,nt,T);
 MC1(:,:,:) = zeros(nt,nt,T);
 MC2(:,:,:) = zeros(nt,nt,T);
 MP1(:,:,:) = zeros(nt,nt,T);
 MP2(:,:,:) = zeros(nt,nt,T);
 D_MW1(:,:,:)=zeros(nt,nt,T);
 D_MW2(:,:,:)=zeros(nt,nt,T);
 D_MP1(:,:,:)=zeros(nt,nt,T);
 D_MP2(:,:,:)=zeros(nt,nt,T);
 M_JH1 = zeros(nt,qnt,T);

 % Filtering %
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) = xi0(1+4*(i-1),t);
 W2(i,t) = xi0(2+4*(i-1),t);
 X1(i,t) = xi0(3+4*(i-1),t);
 X2(i,t) = xi0(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

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

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

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

JH(:,:,t)=JH1(:,:,t)';

% Variance of Y: E[(Y-EY)(Y-EY)'] 

VY(:,:,t)=(JH(:,:,t)')*P0(:,:,t)*(JH(:,:,t))+R; %+ 0.1*eye(5,5); % Variance of Y

% xi(t|t)
xi1(:,t)=xi0(:,t)+P0(:,:,t)*JH(:,:,t)/(VY(:,:,t))*(Y(t,:)'-yhat(:,t));
 
% P(t|t)

P1(:,:,t) = P0(:,:,t)-P0(:,:,t)*JH(:,:,t)/(VY(:,:,t))*JH(:,:,t)'*P0(:,:,t);

% Forecast of xi
xi0(:,t+1) = F*xi1(:,t);

% P0(t+1|t)
P0(:,:,t+1) =  F*P1(:,:,t)*F' + Q;

f(t,1) = -(nr/2)*log(2*pi)-(1/2)*log(det(VY(:,:,t)))...
    -(1/2)*((Y(t,:)'-yhat(:,t))'/(VY(:,:,t))*(Y(t,:)'-yhat(:,t)));

end
end
