function [chain,state] = simulate_markov(x,P,pi0,T);

%  x     = the quantity corresponding to each state, typical element x(i)
%  P     = Markov transition matrix, typical element p(i,j) i,j=1,...n
%  pi0   = probability distribution over initial state
%  T     = number of periods to simulate
%  chain = sequence of realizations from the simulation

n = length(x); %% what is the size of the state vector?
E = rand(T,1); %% T-vector of draws from independent uniform [0,1]  

cumsumP = P*triu(ones(size(P)));
               %% creates a matrix whose rows are the cumulative sums of
               %% the rows of P               
               
%%%%% SET INITIAL STATE USING pi0               

E0   = rand(1,1);
ppi0 = [0,cumsum(pi0)];
s0   = ((E0<=ppi0(2:n+1)).*(E0>ppi0(1:n)))';
s    = s0; 

%%%%% ITERATE ON THE CHAIN

for t=1:T,
    state(:,t) = s;
    ppi        = [0,s'*cumsumP];
    s          = ((E(t)<=ppi(2:n+1)).*(E(t)>ppi(1:n)))';
end

chain = x'*state;