Using the Endogenous Grid Method to solve the Heterogeneous Agent Model.

hetergeneous_agent_model.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Growth Model with Idiosyncratic Shock
% Reference: Huggett, 1997
% Renjie Ge

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% To solve for the optimal decision rule, we use the Endogenous Grid Method.
% This method is legitimate, b/c the policy function is monotone.
% In Step 3 in the main loop, we simulate the economy by generating a million 
% agents and then keep track of these agents' capital holding.
% The mean and variance of earning and consumption are calculated based on 
% the simulation data.

% Note: step 3 here is like the algorithm in HW5_1 


clear all; clc; close all; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Part 1 Initialization %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%---------------------------- Parameters ----------------------------------

alpha=.36; beta=.96; 
sigma=1.5;           % relative risk aversion rate
delta=.1;            % Depreciation rate
N=400;               % No. of grids (first dimension)
S=2;                 % No. of states (second dimension)
gamma=0.1;           % adjustment factor
klow=0;              % lower bound for capital holding

e1=.8; e2=1.2;
e=[e1,e2];                  % endowment distribution
p11=.5; p12=1-p11;
p21=.5; p22=1-p21;
prob=[p11,p12;p21,p22];     % transition matrix

I = ones(N,1)'; J = ones(S,1);
uprime=@(x) x^(-sigma);     % 1st derivative of util func

k0=0; kn=10;  k=linspace(k0,kn,N);
dist=@(x) (x>.2).*10.8104.*(x-0.2)./0.8; % Initial capital distribution (see P395)

%--------------------------- For simulation ------------------------------

M=1000000;                  % No. of agents
agents=linspace(0,1,M);     % Discretized agents
kdist=dist(agents);         % Initial dist of capital
                            % fplot(dist,[0,1]);
K=sum(kdist)/M;             % Initial normalized agg capital stock 

%------------------------- Initialization ---------------------------------

x=zeros(N,1)';
kprime=J*zeros(N,1)';       % Initial guess for policy function. S×N matrix
toler=1e-5;
crit=1;
iter=0;



%%%%%%%%%%%%%%%%%%%%%%%
%%% Part 2 The Loop %%%
%%%%%%%%%%%%%%%%%%%%%%%



tic
while crit>toler
      
      %% Step 1
      % prices solved by zero profit condition
      W=(1-alpha)*K^alpha;
      R=1+alpha*K^(alpha-1)-delta;
      
      % matrix grids for wealth; Dim: S×N
      X=J*k*R+e'*I*W;         
      
      
      %% Step 2  Loop for solving policy functions for state 1 and state 2
      for j=1:200
          
          % b/c the transition rule is the same for both states (p11=p22=0.5),
          % policy function k'(x) (x is wealth) is the same as well, so we only need to solve for one
          % Note: if tran rules are different, we should do another loop using prob(2,s)
          % Note: decision rules for k'(k) here are different though, as we will see later
          
          for i=1:N 
              
              EVprime=0;
              for s=1:S
                  EVprime=prob(1,s)*uprime(R*k(i)+W*e(s)-kprime(s,i))+EVprime; 
              end
              
              x(i)=k(i)+(beta*R*EVprime)^(-1/sigma);  % solve for endogenous grids
              
          end
          
          kprime=interp1(x,k,X,'linear','extrap'); 
          kprime(kprime<klow)=klow;                   % handle the corner solution   
          %plot(k,kprime(1,:),k,kprime(2,:),k,k) 
      end
      policy=interp1(x,k,k,'linear','extrap'); 
      policy(policy<klow)=klow;
      pp=interp1(k,policy,'linear','pp');
      pfun=@(x) ppval(pp,x);        % policy function k'(x)
                                    % fplot(pfun,[0,10])
      
      
      %% Step 3 simulation
      % The idea is that simulate the cross section for a long time until it
      % converges given the policy function from step 2. 
      % Here I just simulated one time because of laziness (wrong)
      
      edist=rand(1,M);
      
      for i=1:M
          if edist(i)<.5
             edist(i)=e(1);
          else
             edist(i)=e(2);
          end
      end
 
      new_kdist=pfun(R*kdist+W*edist);
                                    % plot(agents,new_kdist)
      new_K=sum(new_kdist)/M;
           
      kdist=new_kdist;
      crit=abs(K-new_K);
      K=K+gamma*(new_K-K);
      iter=iter+1;
      
      disp('Iteration')
      disp(iter)
      disp('    Crit      Capital')
      disp([crit,K])
    
end 
toc



%%%%%%%%%%%%%%%%%%%%%%
%%% Part 3 Results %%%
%%%%%%%%%%%%%%%%%%%%%%



disp(['The aggregate capital in stationary equilibrium is: ',num2str(K)])

% decision rules: kt+1 as a function of kt
wealth1=@(x) R*x+W*e1;           % wealth in state 1
pfk1=@(x) pfun(wealth1(x));
wealth2=@(x) R*x+W*e2;           % wealth in state 2
pfk2=@(x) pfun(wealth2(x));
line=@(x) x;                     % 45 degree line

figure (1)
fplot(pfk1,[0,10],'k-.');
hold on;
fplot(pfk2,[0,10],'r--');
fplot(line,[0,10]);
axis([0 10 0 10]);
xlabel('CAPITAL');
ylabel('CAPITAL NEXT PERIOD');
title('OPTIMAL DECISION RULE');
legend('a(k,0.8)', 'a(k,1.2)','45 degree line','location', 'NorthWest');
hold off;

earn=R*kdist+W*edist;
cons=earn-pfun(earn);
disp(['The unconditional mean for earning is: ', num2str(mean(earn))])
disp(['The unconditional variance for earning is: ', num2str(var(earn))])
disp(['The unconditional mean for consumption is: ', num2str(mean(cons))])
disp(['The unconditional variance for consumtpion is: ', num2str(var(cons))])

figure (2)
subplot(1,2,1), hist(earn);
xlabel('EARNING');
ylabel('AGENTS');
subplot(1,2,2), hist(cons);
xlabel('CONSUMPTION');
ylabel('AGENTS');