% IndividSolutionQKPREC.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 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
%            "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] = IndividSolutionPREC(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, convergence criterion
                                         % is not satisfied

C_old = zeros(Nk,Nz);                    % Initialize 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; see eq.(69) in JMMT
                                         % (2016)
   for j = 1:Nz  
       aux1 = interp1(kgrid,q(:,j),'linear','pp');
                                         % Use linear interpolation to generate 
                                         % the piecewise polynomial form of 
                                         % the integrand q(:,j).
       q_coef(:,:,j)= aux1.coefs;        % Compute the coefficients of piecewise 
                                         % linear polynomial; (Nk-1)-by-2-by-Nz 
   end
   
   for ii = 1:2                          % There are two coefficients in a 
                                         % local polynomial function
      qp_coef(:,ii,:)= squeeze(q_coef(:,ii,:))*Prob';
                                         % Find the coefficents of precomputed 
                                         % piecewise linear polynomial; see 
                                         % equation (67) in JMMT (2016)
   end
    
   for j = 1:Nz;                         % For today's state j...
       qp_pol = mkpp(kgrid,qp_coef(:,:,j));  
                                         % Make precomputed piecewise 
                                         % polynomial                                 
       qp = ppval(qp_pol,K(:,j));        % Evaluate the expectation function  
       K1(:,j) = (1-delta+R)*kgrid+W*zgrid(j)-(beta*(qp)).^(-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; 