function S = markov(P,T,S0);
%MARKOV generates a realization of a Markov chain with transition matrix P.
% S = markov(P,T) returns a T x 1 vector S, which is a realization
% of a stationary Markov chain that has transition probabilities
% in P. The elements of S are integers from 1 to n where n is
% the dimension of P.
% S = markov(P,T,S0) uses the initial state S0 for S(1).
% Ellen McGrattan, 9-7-92
% Revised, ERM, 5-23-94
rand('uniform')
n = length(P);
if nargin<3;
S(1) = ceil(rand*n);
else
S(1) = S0;
end;
Q = cumsum(P');
for t = 2:T;
x = find(rand < Q(:,S(t-1)));
S(t) = x(1);
end;
S = S(:);