function [out] = ms_likfun(M,out_ind)

% This function accepts a parameter vector M and an output indicator where
%   out_ind = 1 returns the value of the (negative) of the log likelihood function
%   out_ind = 2 returns the vector ln(f(x(t);M)), where t=1,...T
%   out_ind = 3 returns the vector of filtered probabilities
% 
% For further details, see Kim and Nelson (1999)
%   - the construction here follows that on pg. 67

global y

[obs col] = size(y);

p = M(1,1);
q = M(1,2);
gam1 = M(1,3);
gam2 = M(1,4);
rho = M(1,5);
sig = M(1,6);

p_00 = exp(p) / (1 + exp(p));  % Markov transition probability
p_11 = (exp(q) / (1+exp(q)));  % Markov transition probability

p_0(1,1) = (1 - p_11) / ((1 - p_00) + (1 - p_11));  % Pr[S(t-1)=0|psi(t-1)]
p_1(1,1) = (1-p_00) / ((1 - p_00) + (1 - p_11));    % Pr[S(t-1)=1|psi(t-1)]

p_00_joint(1,1) = p_00 * p_0(1,1);        %  Pr[S(t)=0,S(t-1)=0|psi(t-1)]
p_10_joint(1,1) = (1 - p_11) * p_1(1,1);  %  Pr[S(t)=0,S(t-1)=1|psi(t-1)]
p_01_joint(1,1) = (1 - p_00) * p_0(1,1);  %  Pr[S(t)=1,S(t-1)=0|psi(t-1)]
p_11_joint(1,1) = p_11 * p_1(1,1);        %  Pr[S(t)=1,S(t-1)=1|psi(t-1)]

f_joint = [p_00_joint(1,1);p_01_joint(1,1);p_10_joint(1,1);p_11_joint(1,1)]; 

lik_vect = zeros(obs,1); 
for k = 2:obs
   ind = 1;
   eta=[(1 / sqrt(2 * pi * sig)) * exp(-(y(k,1) - gam1 - rho * y(k-1,1)) ^ 2 / (2 * sig));
        (1 / sqrt(2 * pi * sig)) * exp(-(y(k,1) - gam2 - rho * y(k-1,1)) ^ 2 / (2 * sig));
        (1 / sqrt(2 * pi * sig)) * exp(-(y(k,1) - gam1 - rho * y(k-1,1)) ^ 2 / (2 * sig));
        (1 / sqrt(2 * pi * sig)) * exp(-(y(k,1) - gam2 - rho * y(k-1,1)) ^ 2 / (2 * sig))];
   f_cond = eta' * f_joint;
   
   lik_vect(k,1) = log(f_cond);
   
   p_i0(1,1) = eta(1,1) * f_joint(1,1) / f_cond;
   p_i0(2,1) = eta(3,1) * f_joint(3,1) / f_cond;
   p_0(k,1) = sum(p_i0);
   
   p_i1(1,1) = eta(2,1) * f_joint(2,1) / f_cond;
   p_i1(2,1) = eta(4,1) * f_joint(4,1) / f_cond;
   p_1(k,1) = sum(p_i1);
   
   p_00_joint(k,1) = (exp(p) / (1 + exp(p))) * p_0(k,1);  
   p_10_joint(k,1) = (1 - (exp(q) / (1 + exp(q)))) * p_1(k,1);  
   p_01_joint(k,1) = (1 - (exp(p) / (1 + exp(p)))) * p_0(k,1);  
   p_11_joint(k,1) = (exp(q) / (1 + exp(q))) * p_1(k,1);   
   
   f_joint = [p_00_joint(k,1);p_01_joint(k,1);p_10_joint(k,1);p_11_joint(k,1)]; 
  
end

if out_ind == 1;
   if sig <= 0 
      out = 10000;
   else
      out = -sum(lik_vect);   % Negative because optimization routine looks for a min
   end
elseif out_ind == 2
   out = lik_vect';
elseif out_ind == 3
   out = [p_0 p_1 p_00_joint p_01_joint p_10_joint p_11_joint];
end





