Krusell and Smith (1998): Heterogeneous Agent Models with Aggregate Uncertainty

The original Krusell and Smith (1998) algorithm can be implemented with the toolbox, by transforming the optimization problem of the households to a system of first order conditions and complementarity-slackness conditions.

Though the toolbox is not designed for solving the equilibrium of this type of model directly, since the decision problem is characterized by an equation system (the Euler equation

\[u'(c_t) = \beta \mathbb{E}_t[(1-\delta+r_{t+1})u'(c_{t+1})] + \lambda_t,\]

where \(\lambda_t\) is the Lagrange multiplier on the borrowing constraint, and the complementary-slackness condition, \(\lambda_t k_{t+1}=0\)) with state transition functions, it readily fits in the toolbox’s general framework. One just needs an extra fixed-point loop to update the equilibrium object, which can be coded in MATLAB. In Krusell and Smith (1998), the equilibrium objects are the aggregate state transition coefficients.

The model can be solved within 100 lines of GDSGE code plus 100 lines of MATLAB code. We present the heterogeneous-beta version of the model below.

The gmod File and MATLAB File

ks1998.gmod (you also need the file for grid points from Krusell and Smith (1998) points.kgd)

% Toolbox options
INTERP_ORDER = 4; EXTRAP_ORDER = 4;
SIMU_RESOLVE = 0; SIMU_INTERP = 1;
SaveFreq = inf; PrintFreq = 100; TolEq = 1e-5;
% Parameters
parameters alpha delta phi kMin;
alpha = 0.36;   % capital share in production
delta = 0.025;  % depreciation rate
% phi_k_1 phi_k_0 phi_l_1 phi_l_0; for z in {bad, good}
phi = [
      0.96053     0.095014            0      -1.2228
      0.96072     0.099212            0      -1.1583
    ]';
% States, use the ones as in Krusell and Smith (1998)
var_state kbar k;
kbar = [11.1 11.52 11.94 12.36];
k = [0,1];      % A place holder, will be overwritten
kMin = k(1);    % A place holder, will be overwritten
% Shock process in Krusell and Smith (1998)
var_shock beta e unemp z z_idx;
shock_num = 12;
ez_trans = [
    0.5250 0.3500 0.0312 0.0938
    0.0389 0.8361 0.0021 0.1229
    0.0938 0.0312 0.2917 0.5833
    0.0091 0.1159 0.0243 0.8507
    ];
beta_trans = [
    0.9950 0.0050 0.0000
    0.000625 0.99875 0.000625
    0.0000 0.0050 0.9950
    ];
shock_trans = kron(ez_trans,beta_trans);
z_idx_grid = [0, 1];
z_grid = [0.99, 1.01];
e_grid = [0.00, 0.3271];
unemp_grid = [0.07, 0.00];  % unemployement transfer
beta_grid = [0.9858 0.9894 0.9930];
[beta,e,z] = ndgrid(beta_grid,e_grid,z_grid);
[~,unemp,z_idx] = ndgrid(beta_grid,unemp_grid,z_idx_grid);
beta=beta(:); e = e(:); unemp = unemp(:);
z = z(:); z_idx = z_idx(:);

% State transition functions
var_interp Evp_interp Ev_interp;
initial Evp_interp 1./(k+e+unemp-kMin);
initial Ev_interp log(k+e+unemp-kMin);
% Update
Evp_interp = shock_trans * vp;
Ev_interp = shock_trans * v;

% Tensor
var_tensor kbar_next lbar Y r w budget;
% Forecasting rule
kbar_next = exp(log(kbar).*phi(z_idx*4+1) + phi(z_idx*4+2));
lbar = exp(log(kbar).*phi(z_idx*4+3) + phi(z_idx*4+4));
% Interst rate and wage
Y = z .* kbar.^alpha .* lbar.^(1-alpha)
r = alpha*Y./kbar - delta;
w = (1-alpha)*Y./lbar;
budget = k.*(1+r) + e.*w + unemp;

% Endogenous variables
var_policy k_next lambda;
inbound k_next kMin budget;
inbound lambda 0 1.0;
% Other variables out of the block
var_aux c vp v l;
% Used in simulation
var_output c k_next l;

model;
  c = budget - k_next;
  u_prime = 1/c;
  [Evp_future,Ev_future] = GDSGE_INTERP_VEC(shock,kbar_next,k_next);
  Evp_future = Evp_future;
  euler_residual = -1 + beta*Evp_future/u_prime + lambda;
  vp = (1+r)*u_prime;
  v = log(c) + beta*Ev_future;
  l = e;
  equations;
    euler_residual;
    lambda*(k_next-kMin);
  end;
end;

simulate;
  num_periods = 11000;
  num_samples = 10000;
  initial k 1;
  initial kbar 1;
  initial shock 1;
  var_simu c l;
  k' = k_next;
  kbar' = k_next;    % A place holder
end;
post_simulate_iter;
  SimuRslt.kbar(:,GDSGE_t+1) = mean(SimuRslt.k(:,GDSGE_t+1));
  SimuRslt.lbar(:,GDSGE_t) = mean(SimuRslt.l(:,GDSGE_t));
end;

The MATLAB file that calls the toolbox codes and manually update equilibrium objects main.m

%% Solve the warm-up problems
options.k = csvread('points.kgd'); options.kMin = options.k(1);
IterRslt = iter_ks1998(options);

%% Generate random numbers
% Parameters
num_periods = 11000;
num_samples = 10000; % Need around 8GB to run the simulation

% Transition matrix for aggregate shocks
shock_trans = IterRslt.shock_trans;
% Reduction along the z dimension (beta,e,z) -> (beta',e',z')
shock_trans_zp = sum(reshape(shock_trans,[6,2,6,2]),3);
z_trans = reshape(shock_trans_zp(1,:,:),[2,2]);

% Aggregate shock
rng(0729);
z_idx = gen_discrete_markov_rn(z_trans,1,num_periods,1);
% Idiosyncratic shock conditional on aggregate shock
% Permute to (beta,e) -> (beta',e') | (z,z')
shock_trans_conditional_on_z = permute(reshape(shock_trans,[6,2,6,2]),[1,3,2,4]) ./ reshape(z_trans,[1,1,2,2]);
ind_shock_idx = gen_discrete_markov_rn_with_agg(shock_trans_conditional_on_z, num_samples, num_periods, 6*ones(num_samples,1), z_idx);

simuOptions.num_samples = num_samples;
simuOptions.GEN_SHOCK_START_PERIOD = num_periods;
simuOptions.init.shock = ind_shock_idx + (z_idx-1)*6;
simuOptions.init.k = 11*ones(num_samples,1);
simuOptions.init.kbar = 11*ones(num_samples,1);
simuOptions.EnforceSimuStateInbound = 0;

%% Iterate transition coefficients phi
phi = [
      0.96053     0.095014            0      -1.2228
      0.96072     0.099212            0      -1.1583
    ]';
metric_phi = inf;
tol_phi = 1e-4;
update_speed = 0.5;
iter = 0;
while metric_phi>tol_phi
    options.phi = phi;
    options.WarmUp = IterRslt;
    IterRslt = iter_ks1998(options);
    SimuRslt = simulate_ks1998(IterRslt,simuOptions);
    
    % Collect samples and run regression
    burnPeriods = 1000;
    periodIndex = 1:num_periods;
    r2_z = zeros(1,2);
    rmse_z = zeros(1,2);
    phi_new = phi;
    for i_z=1:2
        periodSet = periodIndex(z_idx==i_z & periodIndex>burnPeriods);
        sample_lbar = SimuRslt.lbar(1,periodSet);
        sample_kbar = SimuRslt.kbar(1,periodSet);
        sample_kbar_next = SimuRslt.kbar(1,periodSet+1);
        mdl = fitlm(log(sample_kbar),log(sample_kbar_next));
        
        phi_new(1,i_z) = mdl.Coefficients.Estimate(2);
        phi_new(2,i_z) = mdl.Coefficients.Estimate(1);
        phi_new(4,i_z) = log(mean(sample_lbar));
        r2_z(i_z) = mdl.Rsquared.Ordinary;
        rmse_z(i_z) = mdl.RMSE;
    end
    
    iter = iter+1;
    metric_phi = max(abs(phi(:) - phi_new(:)));
    fprintf('iter, metric_phi: %d, %g\n',iter, metric_phi);
    % Update
    phi = phi_new*update_speed + phi*(1-update_speed);
end

%% Post analysis
disp('Converged transition coefficients (phi_k_1 phi_k_0 phi_l_1 phi_l_0):');
disp(phi')
% fraction negative wealth
sum(SimuRslt.k(:,num_periods)<0) / numel(SimuRslt.k(:,num_periods))
% Gini coefficients
samplePeriods = num_periods;
k = SimuRslt.k(:,samplePeriods);
% ginicoeff(k(k>=0)) % Uncomment this to calculate gini if you have ginicoeff routine

%%%%%%%%%%%%%%% MAIN ENDS HERE %%%%%%%%%%%%%%%%%%
function shock = gen_discrete_markov_rn_with_agg(trans,numPaths,lenPath,initShock,aggIdx)
% Generate idiosyncratic shocks conditional on aggregate shocks
numStates = size(trans,1);
shock = zeros(numPaths,lenPath);
shock(:,1) = initShock;
cumTrans = cumsum(trans, 2);
for t=1:lenPath-1
    for j=1:numStates
        % Find samples that has shock j
        idxOfShockJ = find(shock(:,t)==j);
        % Fill in numbers based on distribution weights
        un = rand(length(idxOfShockJ),1);
        % Look up trans_j
        [~,shock(idxOfShockJ,t+1)] = histc(un, [0 cumTrans(j,:,aggIdx(t),aggIdx(t+1))]);
    end
end
end