% IndividSolutionQK.m is a routine for constructing solution to the
% individual problem in the Aiyagari model; for a description of the solution 
% procedure see the article "How to Solve Dynamic Stochastic Models Computing
% Expectations Just Once" by Kenneth L. Judd, Lilia Maliar, Serguei Maliar 
% and Inna Tsener, (2016), Quantitative Economics, forthcoming, (henceforth, 
% JMMT, 2016).
%
% -------------------------------------------------------------------------
% Inputs:    "K"      is an individual decision function for capital for 
%                     each of Nz possible productivity states (initial guess)
%            "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 k
%            "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 © 2016 by Lilia Maliar, Serguei Maliar and Inna Tsener. All 
% rights reserved. The code may be used, modified and redistributed under 
% the terms provided in the file "License_Agreement.txt".
% -------------------------------------------------------------------------

function [K] = IndividSolutionQK(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
q = C.^(-gam)*(1-delta+R);               % Compute the integrand funciton 
                                         % on iteration i
                                         
   for j = 1:Nz;                         % For today's state j...
       q2 = interp1(kgrid,q,K(:,j)','*linear');
                                         % Find the state contingent integrand
                                         % function in t+1
       c2 = (q2./(1-delta+R)).^(-1/gam);
                                         % Find the corresponding
                                         % consumption
       c2 = c2.*(c2>1e-20)+0.0001*(c2<=1e-20);
                                         % Restrict consumption to be 
                                         % nonegative
       q2 = c2.^(-gam).*(1-delta+R);     % Compute the state contingent integrand 
                                         % function                                
                                         
       K1(:,j) = (1-delta+R)*kgrid+W*zgrid(j)-(beta*(q2)*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
                                         % the next iteration using damping

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

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

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