% MarkProcessDiscr.m is a routine for discretizing a first-order Markov 
% process with continuous states into a Markov chain with discrete states 
% using the method of Tauchen (1986), Economics Letters 20/2, 177-181. It 
% is used to construct a 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:    "Nz"      is a number of states in the Markov chain
%            "Nsigma"  is a constant that determines the size of the interval
%                      on which shocks can take values
%            "rho"     is a persistence of the log of the productivity level 
%                      in the Markov process
%            "sigma"   is a standard deviation of shocks to the log of the 
%                      productivity level    
%
% Output:    "zgrid"   is a grid for productivity level in the Markov
%                      chain
%            "Prob"    is a matrix of transitional probabilities in the 
%                      Markov chain
%            "PrUnc"   is an unconditional probability of the productivity 
%                      levels in 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 [zgrid, Prob, PrUnc] = MarkovProcessDiscr(rho,sigma,Nsigma,Nz)

zgridlog = linspace(-sigma*Nsigma,sigma*Nsigma,Nz)';
          % Construct a uniformly spaced grid of Nz productivity levels
          % between -sigma*Nsigma and sigma*Nsigma
for i = 1:Nz-1
   zbound(i) = (zgridlog(i)+zgridlog(i+1))/2;
end;      % Define intervals for computing probabilities      

for i = 1:Nz
   Prob(i,1) = normcdf(zbound(1)-rho*zgridlog(i),0,sigma);
          % Probability that a Normally distributed variable is smaller than
          % zbound(1)
   Prob(i,Nz) = 1-normcdf(zbound(Nz-1)-rho*zgridlog(i),0,sigma); 
          % Probability that a Normally distributed variable is larger than
          % zbound(Nz-1)
   for j = 2:Nz-1   
      Prob(i,j) = normcdf(zbound(j)-rho*zgridlog(i),0,sigma)-normcdf(zbound(j-1)-rho*zgridlog(i),0,sigma);
          % Probability that a Normally distributed variable belongs to the 
          % interval zbound(j) and zbound(j-1)
   end
end

% Finding an unconditional probability distribution of the productivity
% levels in the Markov chain
% ---------------------------------------------------------------------
% Consider an example with 3 states. The matrix of transitional probabilities 
% in the Markov chain is:
%         [p11 p12 p13]
% Prob =  [p21 p22 p23]
%         [p31 p32 p33]
% where pij is the probability to be in a state i at t and to move to state
% j at t+1
%
% We must have:
% [p11 p21 p31]   [pie1] = [pie1]            [pie1] = [pie1]
% [p12 p22 p32] x [pie2] = [pie2] = Prob' x  [pie2] = [pie2]
% [p13 p23 p33]   [pie3] = [pie3]            [pie3] = [pie3]
% where [pie1]
%       [pie2]  
%       [pie3]  is an unconditional probability distribution of the 
%               productivity levels in the Markov chain, PrUnc
%
% Since the sum of probabilities is 1, we have: 
% p11*pie1 + p21*pie2 + p31*(1-pie1-pie2) = pie1
% p12*pie1 + p22*pie2 + p32*(1-pie1-pie2) = pie2
% pie3 = 1-pie1-pie2;
%
% Combine the terms to obtain:
% [(p11-p13-1)*pie1 + (p12-p13)*pie2 + p13] = 0
% [(p21-p23)*pie1 + (p22-p23-1)*pie2 + p23] = 0
%
% In a matrix form, we have:
% [(p11-p31-1) (p21-p31)]   `[pie1] = [-p31]
%                          x
% [(p12-p32) (p22-p32-1)]    [pie2] = [-p32] 
%  
% The first matrix is denoted by M:
% [(p11-p31-1) (p21-p31)]
% [(p12-p32)   (p22-p32-1)]
%
ProbT = Prob'; % Transpose the matrix of transitional probabilities in the 
               % Markov chain
M = ProbT(1:Nz-1,1:Nz-1) - ProbT(1:Nz-1,Nz)*ones(1,Nz-1) - eye(Nz-1,Nz-1);
               % Compute the matrix M, defined above
PrUnc(1:Nz-1,:) = -inv(M)*ProbT(1:Nz-1,Nz);
               % Find the first Nz-1 elements of the stationary probability
               % distribution
PrUnc(Nz,:) = 1 - sum(PrUnc(1:Nz-1,:));
               % Find the last element of the unconditional probability
               % distribution of the productivity levels in the Markov chain

% Constructing the grid of the productivity levels
%-------------------------------------------------
zgrid = exp(zgridlog)/(PrUnc'*exp(zgridlog));
               % The mean of productivity levels on the grid is normalized 
               % to one