function f = ekf_chv_wk(x)

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

 f1=1;f2=1;f3=1;f4=1;f5=1;
 nr=5; nt=48;qnt=nt*4;n=5;

% duration dependence matrix : D%
% D starts from t=0 to t=47 : 48 months to consider %
% change this for linear spline (1)%

d11= -1+2*exp(x(5))/(1+exp(x(5)));
d12= -1+2*exp(x(6))/(1+exp(x(6)));
d13= -1+2*exp(x(7))/(1+exp(x(7)));

D = zeros(nt+1,1); % matrix of duration dependence parameter

%{
D(1,1)=0;%0
D(2,1)=0;%1
D(3,1)=d11;%2
D(4,1)=d11+d12;%3
D(5,1)=d11+2*d12;%4
D(6,1)=d11+3*d12;%5

D(7,1)=d11+3*d12+d13;%6
D(8,1)=d11+3*d12+2*d13;%7
D(9,1)=d11+3*d12+3*d13;%8
D(10,1)=d11+3*d12+4*d13;%9
D(11,1)=d11+3*d12+5*d13;%10
D(12,1)=d11+3*d12+6*d13;%11

for i=13:nt+1
D(i,1)=d11+3*d12+6*d13+(i-12)*d14;
end
%}

D = zeros(nt+1,1); % matrix of duration dependence parameter

D(1,1)=0;
D(2,1)=0;

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

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

%(3) variance for v: 6 numbers for q
 
  q1=exp(10*x(1))/(1+exp(10*x(1)))+0.000001;
  q2=exp(10*x(2))/(1+exp(10*x(2)))+0.000001;
  q3=exp(10*x(3))/(1+exp(10*x(3)))+0.000001; 
  q4=exp(10*x(4))/(1+exp(10*x(4)))+0.000001;

 R(1,1)=0.0000000001+0.1*exp(10*x(8))/(1+exp(10*x(8)));
 R(2,2)=0.0000000001+0.1*exp(10*x(9))/(1+exp(10*x(9)));
 R(3,3)=0.0000000001+0.1*exp(10*x(10))/(1+exp(10*x(10)));
 R(4,4)=0.0000000001+0.1*exp(10*x(11))/(1+exp(10*x(11)));
 R(5,5)=0.0000000001+0.1*exp(10*x(12))/(1+exp(10*x(12)));
  
ff=eye(4);
F=zeros(4*nt,4*nt);

F(1:4,1:4)=ff;

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;

%P0(:,:,1) = eye(4*nt,4*nt)*10000;
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);

for i=1:nt
    xx(1,4*(i-1)+1) = i-1;
    xx(1,4*(i-1)+2) = i-1;
    xx(1,4*(i-1)+3) = i-1;
    xx(1,4*(i-1)+4) = i-1;
end

for i=1:qnt-1
tt(i+1,i:i+4*(nt)-1)=xx;
end

tau = tt(1:qnt,1:qnt);

for t=1:T

for i=1:nt
WL(i,t) = xi0(1+4*(i-1),t);
WH(i,t) = xi0(2+4*(i-1),t);
XL(i,t) = xi0(3+4*(i-1),t);
XH(i,t) = xi0(4+4*(i-1),t);
end

for i=1:nt
VXL(4*(i-1)+1,t)=XL(i,t);
VXL(4*(i-1)+2,t)=XL(i,t);
VXL(4*(i-1)+3,t)=XL(i,t);
VXL(4*(i-1)+4,t)=XL(i,t);
VXH(4*(i-1)+1,t)=XH(i,t);
VXH(4*(i-1)+2,t)=XH(i,t);
VXH(4*(i-1)+3,t)=XH(i,t);
VXH(4*(i-1)+4,t)=XH(i,t);
VWL(4*(i-1)+1,t)=WL(i,t);
VWL(4*(i-1)+2,t)=WL(i,t);
VWL(4*(i-1)+3,t)=WL(i,t);
VWL(4*(i-1)+4,t)=WL(i,t);
VWH(4*(i-1)+1,t)=WH(i,t);
VWH(4*(i-1)+2,t)=WH(i,t);
VWH(4*(i-1)+3,t)=WH(i,t);
VWH(4*(i-1)+4,t)=WH(i,t);
end

for i=1:qnt
    for j=1:qnt
        if i>j
      PML(i,j,t) = 1;
      PMH(i,j,t) = 1;
        else
        PML(i,j,t) = exp(-exp(VXL(i,t)+D(tau(i,j)+1,1)));
        PMH(i,j,t) = exp(-exp(VXH(i,t)+D(tau(i,j)+1,1)));
        end
    end
end

PML(1,:,t)=ones(1,qnt);
PMH(1,:,t)=ones(1,qnt);

% (4) computing weekly number unemployed: WU1, WU2 : in row "t", we have the weekly number unemployed at t%
% WU1, WU2: T by qnt matrix%

WUL(t,:) = VWL(:,t)'.*prod(PML(:,:,t),1); % weekly number unemployed of group 1
WUH(t,:) = VWH(:,t)'.*prod(PMH(:,:,t),1); % weekly number unemployed of group 2

for i=1:nt
MUL(t,i) = WUL(t,4*(i-1)+1)+WUL(t,4*(i-1)+2)+WUL(t,4*(i-1)+3)+WUL(t,4*(i-1)+4);% monthly number unemployed of group 1
MUH(t,i) = WUH(t,4*(i-1)+1)+WUH(t,4*(i-1)+2)+WUH(t,4*(i-1)+3)+WUH(t,4*(i-1)+4);% monthly number unemployed of group 2
end

TU(t,:) = MUL(t,:)+MUH(t,:);

%%%%%%%%%%%%%%%%%%%%%%%%
% Measurement equation %
%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%
% jacobian %
%%%%%%%%%%%%

% jacobian wrt w1,w2 %

PPML(t,:) = prod(PML(:,:,t),1); % product of MPP1 matrix
PPMH(t,:) = prod(PMH(:,:,t),1); % product of MPP2 matrix

for i=1:nt
JPML(t,i) = PPML(t,4*(i-1)+1)+PPML(t,4*(i-1)+2)+PPML(t,4*(i-1)+3)+PPML(t,4*(i-1)+4);
JPMH(t,i) = PPMH(t,4*(i-1)+1)+PPMH(t,4*(i-1)+2)+PPMH(t,4*(i-1)+3)+PPMH(t,4*(i-1)+4);
end

% JH_W1/W2 matrix: rows indicate the index of measurement eqn, column indicate
% the variables %

JH_WL = zeros(nt,nt,T);
JH_WH = zeros(nt,nt,T);

for i=1:nt
JH_WL(i,i,t)=JPML(t,i);
JH_WH(i,i,t)=JPMH(t,i);
end

% JH_x1/x2 matrix %
% baseline matrix : nt by dnt matrix (nt is monthly inflows, dnt is weekly number unemployed) % 

% matrix of derivatives (FML, FMH) : f(x,tau) = 1-p(x,tau)
% FML/FMH  is qnt(# of x's) by qnt(# of weekly meas eqn) matrices.

for i=1:qnt
    for j=1:qnt
        if i>j
      FML(i,j,t) = 0;
      FMH(i,j,t) = 0;
        else
        FML(i,j,t) = -exp(VXL(i,t)+D(tau(i,j)+1,1));
        FMH(i,j,t) = -exp(VXH(i,t)+D(tau(i,j)+1,1));
        end
    end
end

FML(1,:,t)=zeros(1,qnt);
FMH(1,:,t)=zeros(1,qnt);

% convert FML and FMH into nt by qnt matrices (sum 1~4 rows, 5~8 rows, ...
% etc) => MFML, MFMH

for i=1:nt
    MFML(i,:,t)=FML(4*(i-1)+1,:,t)+FML(4*(i-1)+2,:,t)+FML(4*(i-1)+3,:,t)+FML(4*(i-1)+4,:,t);
    MFMH(i,:,t)=FMH(4*(i-1)+1,:,t)+FMH(4*(i-1)+2,:,t)+FMH(4*(i-1)+3,:,t)+FMH(4*(i-1)+4,:,t);
end

% nt by qnt matrix
for i=1:nt
JW_XL(i,:,t) = MFML(i,:,t).*WUL(t,:);
JW_XH(i,:,t) = MFMH(i,:,t).*WUH(t,:);
end

% nt(variable) by nt(meas, number unemployed) matrix 

for j=1:nt
JM_XL(:,j,t) = JW_XL(:,4*(j-1)+1,t)+JW_XL(:,4*(j-1)+2,t)+JW_XL(:,4*(j-1)+3,t)+JW_XL(:,4*(j-1)+4,t); 
JM_XH(:,j,t) = JW_XH(:,4*(j-1)+1,t)+JW_XH(:,4*(j-1)+2,t)+JW_XH(:,4*(j-1)+3,t)+JW_XH(:,4*(j-1)+4,t); 
end

% now the row is monthly number unemployed, while column is variables %
JH_XL(:,:,t) = JM_XL(:,:,t)'; 
JH_XH(:,:,t) = JM_XH(:,:,t)';


for i=1:nt
    for j=1:nt
        M_JH1(i,4*(j-1)+1,t) = JH_WL(i,j,t);
        M_JH1(i,4*(j-1)+2,t) = JH_WH(i,j,t);
        M_JH1(i,4*(j-1)+3,t) = JH_XL(i,j,t);
        M_JH1(i,4*(j-1)+4,t) = JH_XH(i,j,t);
    end
end

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

%(3)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

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

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

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

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

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

W(t,1)=xi0(1,t);
W(t,2)=xi0(2,t);
PP(t,1) = exp(-exp(xi0(3,t)));
PP(t,2) = exp(-exp(xi0(4,t)));
end
f =-sum(ll(1:T,1),1); % finally, we want to find the initial values maximizing this
end

