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

function [shocks_mat_v,shocks_mat_u,shocks_mat_eps,shocks_mat_eta] = shocks(prmtr,N,T,y,prior,smoother)  %need to update to produce four outputs

    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, 0;
        1, 0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 1, 0, 0, 0, 0;
        0, 0, 0, 0, 1, 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, 0];

    % Transition matrix of I(0) part
    Fstar = [0, 0, 0, 0, 0;
            1, 0, 0, 0, 0;
            0, 0, 0, 0, 0;
            0, 0, 0, 0, 0;
            0, 0, 0, 0, 0];

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

    % Covariance matrix
    Q = [sig_eps(1)^2, 0, 0, 0, sig_eps(1)^2, 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, sig_eta(1)^2, 0, 0, sig_eta(1)^2, 0;
        sig_eps(1)^2, 0, 0, 0, sig_eps(1)^2, 0, 0, 0;
        0, 0, 0, 0, 0, sig_u^2, 0, sig_u^2;
        0, 0, 0, sig_eta(1)^2, 0, 0, sig_eta(1)^2, 0;
        0, 0, 0, 0, 0, sig_u^2, 0, sig_u^2];

    % Covariance matrix of I(0) part
    Qstar = [sig_eps(1)^2, 0, 0, 0, 0;
            0, 0, 0, 0, 0;
            0, 0, sig_v(1)^2, 0, 0;
            0, 0, 0, sig_eta(1)^2, 0;            
            0, 0, 0, 0, sig_u^2];

    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);
    P_LL = [P_LL(1:3,1:3), zeros(3, 5);
            zeros(1,3), prior, zeros(1,4);
            zeros(1,4), prior, zeros(1,3);
            zeros(1,5), prior, zeros(1,2);
            zeros(1,6), sig_eta(1)^2, 0;
            zeros(1,7), sig_u^2];
        
shocks_mat_v=999*ones(T,N);
shocks_mat_u=999*ones(T,N);
shocks_mat_eta=999*ones(T,N);
shocks_mat_eps=999*ones(T,N);

    i_iter = 1;
    while i_iter <= N % Loop over individuals

        yin = y(1 + (i_iter - 1)*T, 1);
        cin = y(1 + (i_iter - 1)*T, 2);
        
        for z = 1:T 
            if yin == 999
                yin = y(1 + z + (i_iter - 1)*T, 1);
            end
        end
        if yin == 999
            yin = 0;
        end    

        for z = 1:T 
            if cin == 999
                cin = y(1+ z + (i_iter - 1)*T, 2);
            end
        end
        if cin == 999
            cin = 0;
        end    
        
        beta_LL = [0, 0, 0, yin, 0, cin-gamma_ctau*yin, 0, 0]';
            P_LL = reshape(vecP_LL, rowsQstar, rowsQstar);
            P_LL = [P_LL(1:3,1:3), zeros(3, 5);
            zeros(1,3), prior, zeros(1,4);
            zeros(1,4), prior, zeros(1,3);
            zeros(1,5), prior, zeros(1,2);
            zeros(1,6), sig_eta(1)^2, 0;
            zeros(1,7), sig_u^2];
        mu = [0, 0, 0, 0, 0, 0, 0, 0]';

        eta_iter = 1;
        eps_iter = 1;
        v_iter = 1;
        
   BTT_MAT=zeros(T,8);
   PTL_MAT=zeros(T,8^2);
   PTT_MAT=zeros(T,8^2);

        j_iter = 1;
        while j_iter <= T % Loop over time periods
        Q = [sig_eps(eps_iter)^2, 0, 0, 0, sig_eps(eps_iter)^2, 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, sig_eta(eta_iter)^2, 0, 0, sig_eta(eta_iter)^2, 0;
            sig_eps(eps_iter)^2, 0, 0, 0, sig_eps(eps_iter)^2, 0, 0, 0;
            0, 0, 0, 0, 0, sig_u^2, 0, sig_u^2;
            0, 0, 0, sig_eta(eta_iter)^2, 0, 0, sig_eta(eta_iter)^2, 0;
            0, 0, 0, 0, 0, sig_u^2, 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 + (i_iter - 1)*T, 2) == 999)
                beta_TT = beta_TL;
                P_TT = P_TL;

            elseif y(j_iter + (i_iter - 1)*T, 1) == 999  
                vt = y(j_iter + (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;

    shocks_mat_v(j_iter,i_iter) = beta_TT(3)/sig_v(v_iter);
    shocks_mat_u(j_iter,i_iter) = beta_TT(8)/sig_u;

            elseif y(j_iter + (i_iter - 1)*T, 2) == 999
                vt = y(j_iter + (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;

    shocks_mat_eps(j_iter,i_iter) = beta_TT(1)/sig_eps(eps_iter);
    shocks_mat_eta(j_iter,i_iter) = beta_TT(7)/sig_eta(eta_iter);	


            elseif y(j_iter + (i_iter - 1)*T, 1) ~= 999 && y(j_iter + (i_iter - 1)*T, 2) ~= 999
                vt = y(j_iter + (i_iter - 1)*T, 1:2)' - H*beta_TL - A; % 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;

    shocks_mat_v(j_iter,i_iter) = beta_TT(3)/sig_v(v_iter);
    shocks_mat_u(j_iter,i_iter) = beta_TT(8)/sig_u;                
    shocks_mat_eps(j_iter,i_iter) = beta_TT(1)/sig_eps(eps_iter);
    shocks_mat_eta(j_iter,i_iter) = beta_TT(7)/sig_eta(eta_iter);
    
            end

        vecP_TL = reshape(P_TL, [numel(P_TL),1]);
        vecP_TT = reshape(P_TT, [numel(P_TT),1]);            
            
    BTT_MAT(j_iter,1:8)=beta_TT';
    PTL_MAT(j_iter,1:8^2)=vecP_TL';
    PTT_MAT(j_iter,1:8^2)=vecP_TT';            

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

    BETA_PALL=beta_TT;
    P_PALL=P_TT;

j_iter = T-1;
while 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    
                  
    BETA_TALL = BTT_MAT(j_iter,1:8)' + reshape(PTT_MAT(j_iter,1:8^2)',8,8)*F'*inv(reshape(PTL_MAT(j_iter+1,1:8^2)',8,8))*(BETA_PALL - F*BTT_MAT(j_iter,1:8)');
    P_TALL = reshape(PTT_MAT(j_iter,1:8^2)',8,8) + reshape(PTT_MAT(j_iter,1:8^2)',8,8)*F'*inv(reshape(PTL_MAT(j_iter+1,1:8^2)',8,8))*(P_PALL - reshape(PTL_MAT(j_iter+1,1:8^2)',8,8))*inv(reshape(PTL_MAT(j_iter+1,1:8^2)',8,8))'*F*reshape(PTT_MAT(j_iter,1:8^2)',8,8)';

if (smoother==1)
    if (shocks_mat_v(j_iter,i_iter) ~= 999)     
    shocks_mat_v(j_iter,i_iter) = BETA_TALL(3)/sig_v(v_iter); 
    end
    if (shocks_mat_u(j_iter,i_iter) ~= 999)     
    shocks_mat_u(j_iter,i_iter) = BETA_TALL(8)/sig_u; 
    end
    if (shocks_mat_eps(j_iter,i_iter) ~= 999)   
    shocks_mat_eps(j_iter,i_iter) = BETA_TALL(1)/sig_eps(eps_iter); 
    end
    if (shocks_mat_eta(j_iter,i_iter) ~= 999)   
    shocks_mat_eta(j_iter,i_iter) = BETA_TALL(7)/sig_eta(eta_iter); 
    end
end

    BETA_PALL=BETA_TALL;   P_PALL=P_TALL;

j_iter = j_iter - 1;
end        
        
        i_iter = i_iter + 1;
    end % End of loop over individuals
  
j_iter = 1;

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

while j_iter <= T 
    
if (sum(shocks_mat_v(j_iter,1:N)')==999*N) 
    shocks_mat_v(j_iter,1:N)=999*ones(N,1)';   
else   
    v_sub=selif(shocks_mat_v(j_iter,1:N)',shocks_mat_v(j_iter,1:N)' ~= 999);
    if (size(v_sub,1)==N) 
        miss_obs=[]; 
    else
        miss_obs=999*ones(N-size(v_sub,1),1);
    end
    shocks_mat_v(j_iter,1:N)=[sig_v(v_iter)*(v_sub-mean(v_sub))./std(v_sub);miss_obs]';
end

if (sum(shocks_mat_u(j_iter,1:N)')==999*N) 
    shocks_mat_u(j_iter,1:N)=999*ones(N,1)';   
else   
    u_sub=selif(shocks_mat_u(j_iter,1:N)',shocks_mat_u(j_iter,1:N)' ~= 999);
    if (size(u_sub,1)==N) 
        miss_obs=[]; 
    else
        miss_obs=999*ones(N-size(u_sub,1),1);
    end
    shocks_mat_u(j_iter,1:N)=[sig_u*(u_sub-mean(u_sub))./std(u_sub);miss_obs]';
end
   
if (sum(shocks_mat_eta(j_iter,1:N)')==999*N) 
    shocks_mat_eta(j_iter,1:N)=999*ones(N,1)';   
else  
    eta_sub=selif(shocks_mat_eta(j_iter,1:N)',shocks_mat_eta(j_iter,1:N)' ~= 999);
    if (size(eta_sub,1)==N) 
        miss_obs=[]; 
    else
        miss_obs=999*ones(N-size(eta_sub,1),1); 
    end
    shocks_mat_eta(j_iter,1:N)=[sig_eta(eta_iter)*(eta_sub-mean(eta_sub))./std(eta_sub);miss_obs]';
end

if (sum(shocks_mat_eps(j_iter,1:N)')==999*N) 
    shocks_mat_eps(j_iter,1:N)=999*ones(N,1)';   
else   
    eps_sub=selif(shocks_mat_eps(j_iter,1:N)',shocks_mat_eps(j_iter,1:N)' ~= 999);
    if (size(eps_sub,1)==N) 
        miss_obs=[]; 
    else
        miss_obs=999*ones(N-size(eps_sub,1),1); 
    end
    shocks_mat_eps(j_iter,1:N)=[sig_eps(eps_iter)*(eps_sub-mean(eps_sub))./std(eps_sub);miss_obs]';
end
 
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 % End of likelihood function