%% Likelihood Function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = scr_fcn_FD(prmtr,N,T,y)

    prmtr = trans(prmtr);

    theta = prmtr(1);         % MA parameter
    gamma_ctau = prmtr(2);    % Consumption response to permanent income
    gamma_ceps = prmtr(3);    % Consumption response to transitory income
    sig_u = prmtr(4);         % Standard deviation of permanent consumption shock
    sig_eta = prmtr(5:14);    % Standard deviation of permanent income shock
    sig_eps = prmtr(15:26);   % Standard deviation of transitory income shoc
    sig_v = prmtr(27:36);     % Standard deviation of transitory consumption shock (1989 and 1988 missing consumption)

    % Transition matrix
    
      F = [0, 0, 0, 0, 0, 0, 0;
         1, 0, 0, 0, 0, 0, 0;
         0, 1, 0, 0, 0, 0, 0;
         0, 0, 0, 0, 0, 0, 0;
         0, 0, 0, 1, 0, 0, 0;
         0, 0, 0, 0, 0, 0, 0;
         0, 0, 0, 0, 0, 0, 0]; 
   

    % Transition matrix of I(0) part
    Fstar = F;

    % Measurement equation
    Hy = [1, (theta-1), -theta, 0, 0, 1, 0];
    Hc = [gamma_ceps,0,0,1,-1,gamma_ctau,1];
    H = [Hy; Hc];

    % Covariance matrix
    Q = [sig_eps(1)^2, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, sig_v(1)^2, 0, 0, 0;
        0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, sig_eta(1)^2, 0;
        0, 0, 0, 0, 0, 0, sig_u^2];

    % Covariance matrix of I(0) part
    Qstar = Q;

    Ay = 0;
    Ac = 0;
    A = [Ay; Ac];
    
    rowsQstar = size(Qstar, 1);
    % Vectorization of covariance matrix for stationary components
    vecQstar = reshape(Qstar, [numel(Qstar),1]);
    vecP_LL = inv(eye(rowsQstar^2) - kron(Fstar, Fstar))*vecQstar; 
    % Covariance matrix
    P_LL = reshape(vecP_LL, rowsQstar, rowsQstar);
    
    
    val = []; 

    i_iter = 1;
    while i_iter <= N % Loop over individuals
      
        
        beta_LL = [0, 0, 0, 0, 0, 0, 0]';
        P_LL = reshape(vecP_LL, rowsQstar, rowsQstar);
        mu = [0, 0, 0, 0, 0, 0, 0]';

        eta_iter = 1;
        eps_iter = 1;
        v_iter = 1;

        j_iter = 2;
        while j_iter <= T % Loop over time periods
            
             Q = [sig_eps(eps_iter)^2, 0, 0, 0, 0, 0, 0;
                0, 0, 0, 0, 0, 0, 0;
                0, 0, 0, 0, 0, 0, 0;
                0, 0, 0, sig_v(v_iter)^2, 0, 0, 0;
                0, 0, 0, 0, 0, 0, 0;
                0, 0, 0, 0, 0, sig_eta(eta_iter)^2, 0;
                0, 0, 0, 0, 0, 0, sig_u^2];

            beta_TL = mu + F*beta_LL;
            P_TL = F*P_LL*F' + Q;

           if ((y(j_iter + (i_iter - 1)*T, 1) == 999) || (y((j_iter-1) + (i_iter - 1)*T, 1) == 999))  && ((y(j_iter + (i_iter - 1)*T, 2) ~= 999) && (y((j_iter-1) + (i_iter - 1)*T, 2) ~= 999))
       
                vt = (y(j_iter + (i_iter - 1)*T, 2)- y((j_iter-1) + (i_iter - 1)*T, 2))' - Hc*beta_TL - Ac; % Prediction error
                ft = Hc*P_TL*Hc'; % Variance of forecast error 

                invft=inv(ft); 

                beta_TT = beta_TL + P_TL*Hc'*invft*vt;
                P_TT = P_TL - P_TL*Hc'*invft*Hc*P_TL;
                
                val = [val;
                          -0.5*log(((2*pi)^1)*det(ft)) - 0.5*vt'*invft*vt];
                
            elseif ((y(j_iter + (i_iter - 1)*T, 2) == 999) || (y((j_iter-1) + (i_iter - 1)*T, 2) == 999))  && ((y(j_iter + (i_iter - 1)*T, 1) ~= 999) && (y((j_iter-1) + (i_iter - 1)*T, 1) ~= 999))
                vt = (y(j_iter + (i_iter - 1)*T, 1)- y((j_iter-1) + (i_iter - 1)*T, 1))' - Hy*beta_TL - Ay; % Prediction error
                ft = Hy*P_TL*Hy'; % Variance of forecast error 

                invft=inv(ft); 

                beta_TT = beta_TL + P_TL*Hy'*invft*vt;
                P_TT = P_TL - P_TL*Hy'*invft*Hy*P_TL;

                val = [val;
                          -0.5*log(((2*pi)^1)*det(ft)) - 0.5*vt'*invft*vt];
                

            elseif (y(j_iter + (i_iter - 1)*T, 1) ~= 999) && (y((j_iter-1) + (i_iter - 1)*T, 1) ~= 999)  && (y(j_iter + (i_iter - 1)*T, 2) ~= 999) && (y((j_iter-1) + (i_iter - 1)*T, 2) ~= 999)
                vt = (y(j_iter + (i_iter - 1)*T, 1:2)- y((j_iter-1) + (i_iter - 1)*T, 1:2))' - H *beta_TL - A; % Prediction error % Prediction error
                ft = H*P_TL*H'; % Variance of forecast error 
                
                invft=inv(ft); 

                beta_TT = beta_TL + P_TL*H'*invft*vt;
                P_TT = P_TL - P_TL*H'*invft*H*P_TL;

                val = [val;
                          -0.5*log(((2*pi)^2)*det(ft)) - 0.5*vt'*invft*vt];    
                
            else 
                
            beta_TT = beta_TL;
            P_TT = P_TL;
                
            end
            
            
            
 

            beta_LL = beta_TT;
            P_LL = P_TT;

            j_iter = j_iter + 1;

            if j_iter == 3 || j_iter == 4
                eps_iter = eps_iter + 1;
                v_iter = v_iter + 1;
                
            elseif (j_iter > 4 && j_iter < 10) || j_iter == 12 || j_iter == 13
                eps_iter = eps_iter + 1;
                v_iter = v_iter + 1;
                eta_iter = eta_iter + 1;
 
            elseif  j_iter == 10 || j_iter == 11
                eps_iter = eps_iter + 1;
                eta_iter = eta_iter + 1;
            end
        end % End of loop over time periods
               
    i_iter = i_iter + 1;
    end % End of loop over individuals    
    
end % End of likelihood function