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 log likelihood of the sample 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,2);
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(d==1);
if isnull(j1)
   lhout=-1e80;sp=[];
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,:)=(H*pst).*(1../sig').*exp(-.5*(y(t)-x(t,:)*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

function smsp=mcsmooth(H,sp)
% provides for each t p(S_t|data_T)
smsp=zeros(size(sp));
smsp(T,:)=sp(T,:);
for it=T-1:-1:1
   X=H.*(repmat(sp(it,:)',1,ns));
   smsp(it,:)=smsp(it+1,:)*(X./repmat(sum(X')',1,ns);
end

function spdraw=mcdraw(H,sp)
% provides a draw from the pdf p({S_t}|data_T)
T=size(y,1);
spdraw=zeros(T,1);
u=rand(1,1);
p=spdraw(T,:);
for it=T:-1:1
   cdf=cumsum(p);
   if u==0
      spdraw(it)=1;
   else
      spdraw(it)=find(u<=cdf & u> [0 cdf(1:end-1)]);
   end
   p= H(spdraw(it))*sp(it,:)';
   p=p/sum(p);
end

% Code above all works for any parameterization of b, sig, and H.  
% Next function isspecialized to models in which b is not dependent 
% on S (all its columns the same) and sig(S) is constrained only to be monotone
% increasing in S.

function sig=mcsdraw(y,x,b,sdraw)
% draw sig from posterior (assuming flat prior on the log sig's, except for the 
% constraint that sig is monotone increasing in S) conditional 
% on b and sdraw.  This code does not
% assume that b is necessarily the same for all S.
ns=size(b,2);
T=size(y,1);
sdum=find(repmat(sdraw,1,ns)==repmat([1:ns],T,1));
u=repmat(y,1,ns)-x*b;
rsss=sum(u.^2.*sdum);
sigdraw=1:ns;
while any(sig(1:ns-1)>sigdraw(2:ns))
   gv=repmat(randn(T,1).^2,1,ns);
   gv=sum(gv.*sdum); % This is probably not the most efficient way to generate these 
                     % chi-squared random draws.
   sig=sqrt(rsss./gv);
end  % This way of enforcing monotonicity might be very inefficient if the sig values
     % are close together or ill-determined and ns is large.
     
function b = mcbdraw(y,x,sig,sdraw)
% produces a draw from the normal posterior of b conditional on everything else.
ns=size(sig,2);
T=size(y,1);
nx=size(x,2);
ywt=y./sig(sdraw);
xwt=x./repmat(sig(sdraw),nx);
bhat=xwt.\ywt;
[u,d,v]=svd(xwt);
b=bhat+(v/d)*randn(nx,1);

function H=mchdraw(sdraw,oldH)
% produces a Metropolis-Hastings draw from the Dirichlet, modified 
% by the initial steady-state density, that is the posterior on the elements of H.
% The draw is made from the Dirichlet itself, and then a retain/repeat decision is
% made using the initial steady-state distribution as weighting function, so that the
% the target distribution becomes the actual posterior.
ns=size(sig,2);
T=size(y,1);
H=zeros(ns,ns);
for it=2:T
   H(sdraw(it,it-1))=H(sdraw(it,it-1))+1;
end
for is=1:m
   for js=1:m
      H(is,js)=sum(randn(2*H(is,js),1).^2)/2;
   end
end
H=H./repmat(sum(H),ns,1);      
[v,d]=eig(H');
j1=find(d==1);
pbar=v(:,j1)./sum(v(:,j1);
if size(pbar,2)>1
   pbar = mean(pbar,2);% in case of non-ergodic H matrix.
end
[v,d]=eig(oldH');
j1=find(d==1);
oldpbar=v(:,j1)./sum(v(:,j1);
if size(oldpbar,2)>1
   oldpbar = mean(oldpbar,2);% in case of non-ergodic H matrix.
end
r=pbar(sdraw(1))/oldpbar(sdraw(1));
u=rand;
if u>r
   H=oldH;
end
