% IndividSolution.m is a routine for simulating capital series for a given 
% interest rate and a sequence of productivity states. 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:    "K"      is an initial guess on the individual decision function 
%                     for capital for each of Nz possible productivity states 
%            "R"      is an interest rate
%            "Prob"   is a matrix of transitional probabilities in the 
%                     Markov chain
%            "criter" is a convergence criterion
%            "kdamp"  is an updating parameter for decision function
%            "alpha"  is a capital share in Cobb-Douglas production function
%            "Nk"     is a number of grid points in capital (assets)
%            "Nz"     is a number of productivity states
%            "kgrid"  is a grid for capital stock
%            "zgrid"  is a set of possible productivity states
%            "delta"  is a depreciation rate
%            "gam"    is a coefficient of risk aversion
%            "beta"   is a discount factor
%            "kmin"   is a lower bound on capital, coincides with kgrid(1,1)
%            "kmax"   is an upper bound on capital, coincides with kgrid(Nk,1)
%
% Output:    "K"      is an individual decision function for capital for 
%                     each of Nz possible productivity states
% -------------------------------------------------------------------------
% 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 [K] = IndividSolution(K,R,Prob,criter,kdamp,alpha,Nk,Nz,kgrid,zgrid,delta,gam,beta,kmin,kmax)

W = (1-alpha)*(R/alpha)^(alpha/(alpha-1)); 
                                         % Compute wage corresponding to 
                                         % the given interest rate
                                         
dif = inf;                               % Initially, the convergence criterion
                                         % is not satisfied

C_old = zeros(Nk,Nz);                    % Initialize the consumption function
                                         % which is used to check the 
                                         % convergence

i = 1;                                   % The current iteration is 1                                         
                                         
while dif>criter;
   
C = kgrid*ones(1,Nz)*(1-delta+R)+ones(Nk,1)*zgrid'*W-K;
                                         % Compute the consumption function 
                                         % on iteration i
      
   for j = 1:Nz;                         % For today's state j...
       
       k2 = interp1(kgrid,K,K(:,j)','linear');
                                         % Find the state contingent capital
                                         % choice in period t+1
       c2 = K(:,j)*ones(1,Nz)*(1-delta+R)+W*ones(Nk,1)*zgrid'-k2;
                                         % Find the corresponding
                                         % consumption
       c2 = c2.*(c2>1e-20)+0.0001*(c2<=1e-20);
                                         % Restrict consumption to be 
                                         % nonegative
       K1(:,j) = (1-delta+R)*kgrid+W*zgrid(j)-((beta*c2.^(-gam).*(1-delta+R))*Prob(j,:)').^(-1/gam);
                                         % Compute a new capital decision
                                         % function from the Euler 
                                         % equation
   end;    
   
K1 = K1.*(K1>kmin).*(K1<kmax)+kmin*(K1<=kmin)+kmax*(K1>=kmax);
                                         % Restrict the decision to be
                                         % within the given range
   
K = K1*kdamp+K*(1-kdamp);                % Compute the capital function for
                                         % next iteration using damping

dif = max(max(abs(1-C./C_old)));         % Compute the maximum difference
                                         % between the consumption functions
                                         % on current and previous
                                         % iterations

C_old = C;                               % Set the old consumption 
                                         % function to the current one

i = i+1;                                 % Go to the next iteration 
                                   
end; 