## Likelihood Function ##########################################

scr_fcn <- function(prmtr){
  
  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]
  
  Fmat = matrix(0,6,6) # Transition matrix
  Fmat[1,] = c(0, 0, 0, 0, 0, 0)
  Fmat[2,] = c(1, 0, 0, 0, 0, 0)
  Fmat[3,] = c(0, 0, 0, 0, 0, 0)
  Fmat[4,] = c(0, 0, 0, 1, 0, 0)
  Fmat[5,] = c(0, 0, 0, 0, 1, 0)
  Fmat[6,] = c(0, 0, 0, 0, 0, 1)
  
  Fstar = Fmat[1:3,1:3] # Transition matrix of I(0) part
  
  # Measurement equation
  Hy = matrix(c(1, theta, 0, 1, 0, 0),1,6)
  Hc = matrix(c(0, 0, 1, gamma_ctau, gamma_ceps, 1),1,6)
  H = rbind(Hy, Hc)
  
  Q = matrix(0,6,6) # Covariance matrix
  Q[1,] = c(sig_eps[1]^2, 0, 0, 0, sig_eps[1]^2, 0)
  Q[2,] = c(0, 0, 0, 0, 0, 0)
  Q[3,] = c(0, 0, sig_v[1]^2, 0, 0, 0)
  Q[4,] = c(0, 0, 0, sig_eta[1]^2, 0, 0)
  Q[5,] = c(sig_eps[1]^2, 0, 0, 0, sig_eps[1]^2, 0)
  Q[6,] = c(0, 0, 0, 0, 0, sig_u^2)
  
  Qstar = Q[1:3,1:3] # Covariance matrix of I(0) part
  
  Ay = 0
  Ac = 0
  A = matrix(c(Ay, Ac),2,1)
  
  # Vectorization of covariance matrix for stationary components
  vecQstar = matrix(Qstar, ncol = 1)
  vecP_LL = solve(diag(nrow(Qstar)^2) - Fstar%x%Fstar)%*%vecQstar 
  # Covariance matrix
  P_LL = matrix(vecP_LL, nrow(Qstar), nrow(Qstar))
  P_LL = cbind(P_LL, matrix(0,nrow(Qstar),nrow(Qstar)))
  P_LL = rbind(P_LL, c(matrix(0,1,nrow(Qstar)), prior, 0, 0),
                c(0, 0, 0, 0, prior, 0),
                c(0, 0, 0, 0, 0, prior))
  
  val <- vector()
  
  i_iter = 1
  while(i_iter <= N){ # Loop over individuals
  	
    yin = y[1 + (i_iter - 1)*Ts, 1]

    
    for(z in 1:Ts){
      if(yin == 999){
        yin = y[1 + z + (i_iter - 1)*Ts, 1]
      }
    } 
    
    cin = y[1 + (i_iter - 1)*Ts, 2]
         
    for(z in 1:Ts){
      if(cin == 999){
        cin = y[1+ z + (i_iter - 1)*Ts, 2]
      }
    } 
     
      
  beta_LL = matrix(c(0, 0, 0, yin, 0, cin - gamma_ctau*yin),6,1)
  P_LL = matrix(vecP_LL, nrow(Qstar), nrow(Qstar))
  P_LL = cbind(P_LL, matrix(0,nrow(Qstar),nrow(Qstar)))
  P_LL = rbind(P_LL, c(matrix(0,1,nrow(Qstar)), prior, 0, 0),
                c(0, 0, 0, 0, prior, 0),
                c(0, 0, 0, 0, 0, prior))
  mu = matrix(c(0, 0, 0, 0, 0, 0),6,1)
  
    
    eta_iter = 1
    eps_iter = 1
    v_iter = 1
    
    j_iter = 1
    while(j_iter <= Ts){ # Loop over time periods
      Q[1,] = c(sig_eps[eps_iter]^2, 0, 0, 0, sig_eps[eps_iter]^2, 0)
      Q[2,] = c(0, 0, 0, 0, 0, 0)
      Q[3,] = c(0, 0, sig_v[v_iter]^2, 0, 0, 0)
      Q[4,] = c(0, 0, 0, sig_eta[eta_iter]^2, 0, 0)
      Q[5,] = c(sig_eps[eps_iter]^2, 0, 0, 0, sig_eps[eps_iter]^2, 0)
      Q[6,] = c(0, 0, 0, 0, 0, sig_u^2)
      
      beta_TL = mu + Fmat%*%beta_LL
      P_TL = Fmat%*%P_LL%*%t(Fmat) + Q
      
      if(y[j_iter + (i_iter - 1)*Ts, 1] == 999 && y[j_iter + (i_iter - 1)*Ts, 2] == 999){
        beta_TT = beta_TL
        P_TT = P_TL
        
      } else if(y[j_iter + (i_iter - 1)*Ts, 1] == 999){
        vt = y[j_iter + (i_iter - 1)*Ts, 2] - Hc%*%beta_TL - Ac # Prediction error
        ft = Hc%*%P_TL%*%t(Hc) # Variance of forecast error
        
        invft = solve(ft)
        
        beta_TT = beta_TL + P_TL%*%t(Hc)%*%invft%*%vt
        P_TT = P_TL - P_TL%*%t(Hc)%*%invft%*%Hc%*%P_TL
        
        if (j_iter >= start){
        val <- rbind(val, -0.5*log(((2*pi)^1)*det(ft)) - 0.5*t(vt)%*%solve(ft)%*%vt)}
        
        
      } else if(y[j_iter + (i_iter - 1)*Ts, 2] == 999){
        vt = y[j_iter + (i_iter - 1)*Ts, 1] - Hy%*%beta_TL - Ay # Prediction error
        ft = Hy%*%P_TL%*%t(Hy) # Variance of forecast error 
        
        invft = solve(ft)
        
        beta_TT = beta_TL + P_TL%*%t(Hy)%*%invft%*%vt
        P_TT = P_TL - P_TL%*%t(Hy)%*%invft%*%Hy%*%P_TL
        
        if (j_iter >= start){
        val <- rbind(val, -0.5*log(((2*pi)^1)*det(ft)) - 0.5*t(vt)%*%solve(ft)%*%vt)}
        
      } else if(y[j_iter + (i_iter - 1)*Ts, 1] != 999 && y[j_iter + (i_iter - 1)*Ts, 2] != 999){
        vt = y[j_iter + (i_iter - 1)*Ts, 1:2] - H%*%beta_TL - A # Prediction error
        ft = H%*%P_TL%*%t(H) # Variance of forecast error 
        
        invft = solve(ft)
        
        beta_TT = beta_TL + P_TL%*%t(H)%*%invft%*%vt
        P_TT = P_TL - P_TL%*%t(H)%*%invft%*%H%*%P_TL
        
        if (j_iter >= start){
        val <- rbind(val, -0.5*log(((2*pi)^2)*det(ft)) - 0.5*t(vt)%*%solve(ft)%*%vt)}
      }
      
      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
        
      } else if((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      
        
      } else if(j_iter == 10 || j_iter == 11){
        eps_iter = eps_iter + 1
        eta_iter = eta_iter + 1
      }
    } # End of loop over time periods
    
    i_iter = i_iter + 1
  } # End of loop over individuals
   
  return(val) # Negative log likelihood value
  
} # End of likelihood function