function [lh,sp]=mclh(y,x,sig,b,H)
% function [lh,sp]=mclh(y,x,sig,b,H)
% The model is 
%              y(t)=x(t)b(S(t))+eps(t)sig(S(t)), t=1,...,T
%              P(S(t)=i|S(t-1)=j)= h(i,j), t=1,...,T
%              eps(t) i.i.d. N(0,1), independent of current and past x, past y
%              [P(S(0)=i)]=pbar, the steady-state distribution implied by H.
%              If x contains lagged y's, those occurring in x(1) are conditioned on, i.e.
%              treated as non-random.
% This function evaluates the likelihood of the sample (actually, a vector whose product is
% the likelhihood) and, for each date t in 
% the sample, p(S_t|data_t) as a function of S_t.  
% The b argument is a kxm matrix, where k is the number of x's and m is the number of states.
% The sig argument is 1xm.
ns=size(H,1);
T=size(y,1);
m=size(b,1);
if size(b,2)~=ns, error('b/H dimension mismatch'),end
if size(sig,2)~=ns, error('sig/H dimension mismatch'),end
if size(x,1)~=T, error('x/y length mismatch'),end
[v,d]=eig(H');
j1=find(abs(1-diag(d))<eps*100);
if isempty(j1) | any(sig<eps*100)
   lh=-1e80*ones(T,1);sp=.5*ones(T,2);
else
   pbar=v(:,j1)/sum(v(:,j1));
   if size(pbar,2)>1
   pbar = mean(pbar,2);% in case of non-ergodic H matrix (more than one unit root).
   end
   pst=pbar';
   lh=zeros(T,1); % time series of likelihood increments, like standardized residuals
   sp= zeros(T,ns); % Filtered state probabilities.
   sr2pfac=1/sqrt(2*pi);
   for it=1:T
      sp(it,:)=(pst*H').*(1../sig).*exp(-.5*(y(it)-x(it,:)*b).^2../sig.^2)*sr2pfac;
      lh(it)=sum(sp(it,:)); % i.e., lh is a weighted sum of the Gaussian likelihoods for
                            % for the separate states.
      sp(it,:)=sp(it,:)/sum(sp(it,:));
      pst=sp(it,:);
   end
end
