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(:); 
