function S = carterKohnAR1Restricted(y,yLag,rho,h,P)
%{
Try Carter-Kohn draws for many states with AR1 structures
This version puts in a special restriction about ending up in the last
state

Dick Startz
March 2017
%}
n = length(y);  % number of data points
M = length(rho);  % number of states

pLag = nan(n,M);    %p(S|lagged info)
pCurrent = nan(n,M);       %p(S|current info)
pStar = nan(n,M);   %temp for probability calculation
pSteadyState = steadyStateProbabilities(P);

%forward filter
for t = 1:n
    % Step 1
    if t>1
        for state_t = 1:M
            temp = 0;
            for m=1:M
                temp = temp + P(m,state_t)*pCurrent(t-1,m);
            end
            pLag(t,state_t) = temp;
        end
    else
        for state_t = 1:M
            temp = 0;
            for m=1:M
                temp = temp + P(m,state_t)*pSteadyState(m);
            end
            pLag(t,state_t) = temp;
        end
    end
    
    % Step 2
    for state_t = 1:M
        pStar(t,state_t)...
            = normpdf(y(t),rho(state_t)*yLag(t),1/sqrt(h(state_t)))...
            *pLag(t,state_t);
    end
    % Step 3
    for state_t = 1:M
        pCurrent(t,state_t) = pStar(t,state_t)/sum(pStar(t,:));
    end
end

%backward filter
S = nan(n,1);
S(n) = myRandomState(pCurrent(n,:));
if (y(n) > 0) && S(n) == 3
    S(n) = 1;
end

prob = nan(M,1);

for t = n-1:-1:1
    for state_t = 1:M
        prob(state_t) = ...
            P(state_t,S(t+1))*pCurrent(t,state_t)/pLag(t+1,S(t+1));
    end
    S(t) = myRandomState(prob);
    if (y(t) > 0) && S(t) == 3
        S(t) = 1;
    end
end
end