% IndividSimulation.m is a routine for simulating capital series for a given 
% sequence of productivity states. It is used to simulate 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:    "kgrid"     is a grid for capital stock
%            "kmin"      is a lower bound on capital, coincides with kgrid(1,1)
%            "kmax"      is an upper bound on capital, coincides with kgrid(Nk,1)
%            "K"         is an individual decision function for capital for 
%                        each of Nz possible productivity states
%            "state"     is a vector of simulated series of the states of the Markov chain
%            "ksim_init" is an initial state of the Markov chain
%
% Output:    "ksim"      is a vector of simulated capital series
%            "kmean"     is an average capital stock in the economy
% -------------------------------------------------------------------------
% 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 [kmean,ksim] = IndividSimulation(kgrid,kmin,kmax,K,state,ksim_init)  

T   = size(state,1);     % Simulation length is inferred from the length of  
                         % the given random draw

ksim(1,1) = ksim_init;   % Initial value of capital stock 
for t = 1:T
    ksim(t+1,1) = interp1(kgrid,K(:,state(t,1)),ksim(t,1),'linear');
                         % Find the next-period capital stock by interpolating 
                         % the decision function K
    ksim(t+1,1) = ksim(t+1,1)*(ksim(t+1,1)>kmin).*(ksim(t+1,1)<kmax)+kmin*(ksim(t+1,1)<=kmin)+kmax*(ksim(t+1,1)>=kmax);
                         % Restrict ksim(t+1,1) to be within the admissible
                         % interval
end
                         
kmean = mean(ksim);      % Compute the average capital stock in the economy        