% This MATLAB software solves a version of the Bewley-Aiyagari-Huggett model 
% with uninsurable idiosyncratic risk and occasionally binding borrowing 
% constraints by using an Euler equation method that parameterizes capital 
% function jointly with Q function, as described in 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).
%
% This software is based on that of Lilia  Maliar and Serguei Maliar for 
% solving a version of the Bewley-Aiyagari-Huggett model with uninsurable 
% idiosyncratic risk and occasionally binding borrowing constraints by using 
% an Euler equation method, as described 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). 
%
% First version: May 15, 2015
% This version:  November 15, 2016
% ------------------------------------------------------------------------
% This software uses the following files: 
% ------------------------------------------------------------------------
% 1. "Main_AiyagariQK.m"    solves the model and evaluates the accuracy by 
%                           computing the residuals in the Euler equations
% 2. "IndividSolutionQK.m"  computes a solution to the individual problem
%                           satisfying the Euler equation for the given 
%                           aggregate capital (and interest rate and wage) 
% 3. "IndividSimulation.m"  simulates a time-series solution of the individual 
%                           problem for a given sequence of shocks
% 4. "MarkovProcessDiscr.m" discretizes a Markov process with continuous 
%                           states into a Markov chain with discrete 
%                           states using a method of Tauchen (1986), 
%                           Economics Letters 20/2, 177-181. 
% 5. "MarkovChainSimul.m"   simulates a Markov chain for a given sequence
%                           of shocks 
% 6. "epsi2.mat"            contains two fixed sequences of uniformly 
%                           distributed random numbers of lenghts 10,000 
%                           and 10,200 for computing a solution and for 
%                           evaluating accuracy, respectively
% -------------------------------------------------------------------------
% 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".
% -------------------------------------------------------------------------

clc;
clear all;
tic;

% 1. Model parameters
% -------------------
alpha   = 0.36;      % Capital share in production 
beta    = 0.96;      % Discount factor
delta   = 0.08;      % Depreciation rate
gam     = 3;         % Risk aversion parameter

% 2. Discretizing a first-order Markov process
% --------------------------------------------
Nz    =  7;         % Number of states in the Markov chain
rho     = .6;        % Persistence of the log of the productivity level
                     % of the Markov process
sigma   = .2*(1-rho^2)^0.5; 
                     % Standard deviation of shocks to the log of the 
                     % productivity level   
Nsigma  =  5;        % A constant that determines the size of the interval
                     % [-sigma*Nsigma, sigma*Nsigma] on which shocks can 
                     % take values                                                     
[zgrid, Prob,PrUnc] = MarkovProcessDiscr(rho,sigma,Nsigma,Nz);
                     % Construct a grid for productivity levels, 
                     % transition probability and stationary 
                     % distribution of productivity levels

% 3. Simulating Markov chain
% --------------------------
T       = 10000;     % The length of simulation in the solution procedure
%epsi    = rand(T,1);% A draw of uniformly distributed random numbers of 
                     % length T for constructing a solution
load epsi2;          % Load previously stored sequences of uniformly 
                     % distributed random numbers "epsi" and "epsi_test"
                     % for constructing solution and evaluating the
                     % residuals, respectively
state_init  = 4;     % Initial state for simulation is 4                     
state   = MarkovChainSimul(Prob,epsi,state_init);
                     % Produce a simulated series of productivity states 
                     % for constructung solution for the given draw "epsi"
                                          
% 4. Algorithm's parameters
% -------------------------
Nk      = 1000;      % Number of grid points in assets
criter  = 1e-10;     % Convergence criterion
kdamp   = 0.5;       % Damping parameter for a (fixed-point) iteration on 
                     % the capital policy function; a default value "1" means 
                     % no damping and a smaller value in the interval (0,1)
                     % means that capital policy function is updated only 
                     % partially from one iteration to another to enhance 
                     % convergence 

% 5. Grid for capital stock
% -------------------------
kmin    = 0;         % Lower bound on capital is zero (no borrowing)    
kmax    = delta^(1/(alpha-1));
                     % Upper bound on capital is maximum sustainable 
                     % capital 
kgrid(1:900,1) = linspace(kmin,kmin+3,900)';      
                     % Dense grid for assets near the kink [kmin,kmin+3]
kgrid(900:1000,1) = linspace(kmin+3,kmax,101)';      
                     % Less dense grid for assets away from the kink
                     % [kmin+3,kmax]

                     
% 6. Initial guess and initial condition for capital
% --------------------------------------------------
K = kgrid*ones(1,Nz);  % Initial guess about the decision function for capital;
                       % Nk-by-Nz 
kss = ((1/beta-1+delta)/alpha)^(1/(alpha-1));
                       % Steady state level of capital
ksim_init = kss;       % Initial individual capital for simulation is steady 
                       % state      
R2 = 1/beta-1+delta;   % Upper bound on the interest rate is constructed from
                       % the steady state of the deterministic economy
R1 = delta;            % Lower bound on the interest rate is depreciation

dif = inf;             % Initially, convergence criterion is not satisfied
                                         
i = 1;                 % The current iteration is 1                                         

% 7. Computing a solution to the model
% ------------------------------------
while dif > 1e-10;	   % To compute the equilibrium interest rate, we use a 
                       % bisection method as in Aiyagari (1994)       
                    
    R = (R1+R2)/2;     % Set the interest rate to the average of the upper and 
                       % lower bounds
                    
    fprintf(1,'Iteration = %i:, Lower bound, R1 = %.4f, Current interest rate, R = %.4f, Upper bound, R2 = %.4f\n\n',i,R1*100,R*100,R2*100);
                       % Display the results on the current iteration                
                    
    [K] = IndividSolutionQK(K,R,Prob,criter,kdamp,alpha,Nk,Nz,kgrid,zgrid,delta,gam,beta,kmin,kmax);
                       % Compute a solution to the individual problem
    [kmean] = IndividSimulation(kgrid,kmin,kmax,K,state,ksim_init); 
                       % Simulate the model to find aggregate capital 
    if kmean>=(R/alpha)^(1/(alpha-1));
        R2 = R;	       % If aggregate capital is too high, reduce the 
                       % upper bound on the interest rate        
    else
        R1 = R;        % Otherwise, increase the lower bound on the 
                       % interest rate
    end;
    
   dif = abs(1-R2/R1); % Compute the difference between upper and lower bounds
                       % on the interest rate
    
    i = i+1;           % Go to the next iteration 
    
end;

CPU = toc;             % Store running time


% 8. Simulating a Markov chain to produce a set of points for evaluating the
% residuals in the Euler equation
% -------------------------------------------------------------------------
T_test = 10200;     % The length of simulation for accuracy evaluation
% epsi_test = rand(T_test,1); 
                    % A draw of uniformly distributed random numbers of 
                    % length T_test for constructing solution
state_init  = 4;    % Initial state for simulation is 4                     
state_test   = MarkovChainSimul(Prob,epsi_test,state_init);
                    % Produce a simulated series of productivity states 
                    % for evaluating the residuals for the draw "epsi_test"
ksim_init = kss;    % Initial individual capital for simulation is steady 
                    % state    
                    
[kmean_test,ksim_test] = IndividSimulation(kgrid,kmin,kmax,K,state_test,ksim_init); 
                    % Simulate the model to produce capital series and
                    % aggregate capital 
                      
W = (1-alpha)*(R/alpha)^(alpha/(alpha-1)); 
                    % Compute wage corresponding to the given interest rate
CSS = kgrid*ones(1,Nz)*(1-delta+R)+ones(Nk,1)*zgrid'*W-K;
                    % Compute the consumption decision function; Nk-by-Nz
QSS = CSS.^(-gam).*(1-delta+R);  
                    % Compute the decision function for the integrand; 
                    % Nk-by-Nz
                    
for j = 1:Nz  
    aux1 = interp1(kgrid,QSS(:,j),'linear','pp');
                    % Use linear interpolation to generate the piecewise 
                    % polynomial form of the integrand QSS.
    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 t = 1:T_test; 
    c1 = ksim_test(t,1)*(1-delta+R)+W*zgrid(state_test(t,1),1)-ksim_test(t+1,1);
                    % Compute consumption today
    qp_pol = mkpp(kgrid,qp_coef(:,:,state_test(t,1)));  
                    % Make precomputed piecewise polynomial                                 
    qp = ppval(qp_pol,ksim_test(t+1,1)); 
                    % Evaluate the expectation function of the Euler
                    % equation 
    residuals(t,1) = ((beta*qp)/c1^(-gam)-1).*(ksim_test(t+1,1)>kmin);
                    % Construct the unit-free residuals in the Euler equation 
                    % for those periods in which the borrowing limit is 
                    % not binding
end

discard = 200;      % The number of observations to discard 

residuals = residuals(1+discard:T_test,:);                           
                    % We discard the first 200 observations to eliminate
                    % the effect of the initial conditions
                    
residuals =  residuals(find(residuals~=0));   
                    % Drop the observations in which the residual is zero
                    % (i.e., the borrowing constraint is binding)

Mean_Residuals   =  mean(log10(abs(residuals)));
                    % The average residuals on the simulated path
Max_Residuals    =  max(log10(abs(residuals)));                    
                    % The maximum residuals on the simulated path

fprintf(1,'ACCURACY EVALUATION AND RUNNING TIME:\n\n');
fprintf(1,'Running time = %.2f, Mean residuals = %.2f, Max residuals = %.2f\n\n',CPU,Mean_Residuals,Max_Residuals);
                    % Display the running time and residuals


% 9. Display the output on the graphs 
% ------------------------------------
figure(1)

subplot(2,2,1)
mesh(kgrid,zgrid,K');
xlabel('k')
ylabel('z')
zlabel('K(k,z)')
title('Asset function')

subplot(2,2,2)
mesh(kgrid,zgrid,CSS');
xlabel('k')
ylabel('z')
zlabel('C(k,z)')
title('Consumption function') 
