% MarkChainSimul.m is a routine for simulating a Markov chain with a given
% transition probability matrix. It is used to construct numerical solution 
% to the Bewley-Aiyagari-Huggett model in the article "The Neoclassical 
% Growth Model with Heterogeneous Quasi-Geometric Consumers", by Lilia Maliar 
% and Serguei Maliar, Journal of Money, Credit, and Banking 38(3), 635-654 
% (henceforth, MM, 2006).  
% -------------------------------------------------------------------------
% Inputs:    "Prob"       is a matrix of transitional probabilities in the 
%                         Markov chain
%            "raux"       is a sequence of uniformly distributed random numbers
%            "state_init" is initial state of the Markov chain
%
% Output:    "state"      is a simulated series of the states of the Markov 
%                         chain
% -------------------------------------------------------------------------
% Copyright © 2002-2016 by Lilia Maliar and Serguei Maliar. All rights 
% reserved. The code may be used, modified and redistributed under the 
% terms provided in the file "License_Agreement.txt".
% -------------------------------------------------------------------------
 
function state = MarkovChainSimul(Prob,raux,state_init)

T   = size(raux,1);     % Simulation length is inferred from the length of  
                        % the given random draw
Nz = size(Prob,2);      % The number of states in the chain is inferred
                        % from transitional probabilities

state(1,1) = state_init;% The initial state is set at state_init

for t = 1:T

    probaux = Prob(state(t,1),:);
                       % Transitional probability from current state to 
                       % all states; 1-by-Nz
    for j = 1:Nz
        if raux(t,1)<= sum(probaux(1:j))        
            state(t+1,1) = j; 
                       % If the random draw raux falls into the interval of 
                       % the probability of the current state to be less
                       % or equal to j, we set the current state to be j
        break
        end        
    end
end