# 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)

  1% Toolbox options
2INTERP_ORDER = 4; EXTRAP_ORDER = 4;
3SIMU_RESOLVE = 0; SIMU_INTERP = 1;
4SaveFreq = inf; PrintFreq = 100; TolEq = 1e-5;
5% Parameters
6parameters alpha delta phi kMin;
7alpha = 0.36;   % capital share in production
8delta = 0.025;  % depreciation rate
9% phi_k_1 phi_k_0 phi_l_1 phi_l_0; for z in {bad, good}
10phi = [
11      0.96053     0.095014            0      -1.2228
12      0.96072     0.099212            0      -1.1583
13    ]';
14% States, use the ones as in Krusell and Smith (1998)
15var_state kbar k;
16kbar = [11.1 11.52 11.94 12.36];
17k = [0,1];      % A place holder, will be overwritten
18kMin = k(1);    % A place holder, will be overwritten
19% Shock process in Krusell and Smith (1998)
20var_shock beta e unemp z z_idx;
21shock_num = 12;
22ez_trans = [
23    0.5250 0.3500 0.0312 0.0938
24    0.0389 0.8361 0.0021 0.1229
25    0.0938 0.0312 0.2917 0.5833
26    0.0091 0.1159 0.0243 0.8507
27    ];
28beta_trans = [
29    0.9950 0.0050 0.0000
30    0.000625 0.99875 0.000625
31    0.0000 0.0050 0.9950
32    ];
33shock_trans = kron(ez_trans,beta_trans);
34z_idx_grid = [0, 1];
35z_grid = [0.99, 1.01];
36e_grid = [0.00, 0.3271];
37unemp_grid = [0.07, 0.00];  % unemployement transfer
38beta_grid = [0.9858 0.9894 0.9930];
39[beta,e,z] = ndgrid(beta_grid,e_grid,z_grid);
40[~,unemp,z_idx] = ndgrid(beta_grid,unemp_grid,z_idx_grid);
41beta=beta(:); e = e(:); unemp = unemp(:);
42z = z(:); z_idx = z_idx(:);
43
44% State transition functions
45var_interp Evp_interp Ev_interp;
46initial Evp_interp 1./(k+e+unemp-kMin);
47initial Ev_interp log(k+e+unemp-kMin);
48% Update
49Evp_interp = shock_trans * vp;
50Ev_interp = shock_trans * v;
51
52% Tensor
53var_tensor kbar_next lbar Y r w budget;
54% Forecasting rule
55kbar_next = exp(log(kbar).*phi(z_idx*4+1) + phi(z_idx*4+2));
56lbar = exp(log(kbar).*phi(z_idx*4+3) + phi(z_idx*4+4));
57% Interst rate and wage
58Y = z .* kbar.^alpha .* lbar.^(1-alpha)
59r = alpha*Y./kbar - delta;
60w = (1-alpha)*Y./lbar;
61budget = k.*(1+r) + e.*w + unemp;
62
63% Endogenous variables
64var_policy k_next lambda;
65inbound k_next kMin budget;
66inbound lambda 0 1.0;
67% Other variables out of the block
68var_aux c vp v l;
69% Used in simulation
70var_output c k_next l;
71
72model;
73  c = budget - k_next;
74  u_prime = 1/c;
75  [Evp_future,Ev_future] = GDSGE_INTERP_VEC(shock,kbar_next,k_next);
76  Evp_future = Evp_future;
77  euler_residual = -1 + beta*Evp_future/u_prime + lambda;
78  vp = (1+r)*u_prime;
79  v = log(c) + beta*Ev_future;
80  l = e;
81  equations;
82    euler_residual;
83    lambda*(k_next-kMin);
84  end;
85end;
86
87simulate;
88  num_periods = 11000;
89  num_samples = 10000;
90  initial k 1;
91  initial kbar 1;
92  initial shock 1;
93  var_simu c l;
94  k' = k_next;
95  kbar' = k_next;    % A place holder
96end;
97post_simulate_iter;
98  SimuRslt.kbar(:,GDSGE_t+1) = mean(SimuRslt.k(:,GDSGE_t+1));
99  SimuRslt.lbar(:,GDSGE_t) = mean(SimuRslt.l(:,GDSGE_t));
100end;


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

  1%% Solve the warm-up problems
2options.k = csvread('points.kgd'); options.kMin = options.k(1);
3IterRslt = iter_ks1998(options);
4
5%% Generate random numbers
6% Parameters
7num_periods = 11000;
8num_samples = 10000; % Need around 8GB to run the simulation
9
10% Transition matrix for aggregate shocks
11shock_trans = IterRslt.shock_trans;
12% Reduction along the z dimension (beta,e,z) -> (beta',e',z')
13shock_trans_zp = sum(reshape(shock_trans,[6,2,6,2]),3);
14z_trans = reshape(shock_trans_zp(1,:,:),[2,2]);
15
16% Aggregate shock
17rng(0729);
18z_idx = gen_discrete_markov_rn(z_trans,1,num_periods,1);
19% Idiosyncratic shock conditional on aggregate shock
20% Permute to (beta,e) -> (beta',e') | (z,z')
21shock_trans_conditional_on_z = permute(reshape(shock_trans,[6,2,6,2]),[1,3,2,4]) ./ reshape(z_trans,[1,1,2,2]);
22ind_shock_idx = gen_discrete_markov_rn_with_agg(shock_trans_conditional_on_z, num_samples, num_periods, 6*ones(num_samples,1), z_idx);
23
24simuOptions.num_samples = num_samples;
25simuOptions.GEN_SHOCK_START_PERIOD = num_periods;
26simuOptions.init.shock = ind_shock_idx + (z_idx-1)*6;
27simuOptions.init.k = 11*ones(num_samples,1);
28simuOptions.init.kbar = 11*ones(num_samples,1);
29simuOptions.EnforceSimuStateInbound = 0;
30
31%% Iterate transition coefficients phi
32phi = [
33      0.96053     0.095014            0      -1.2228
34      0.96072     0.099212            0      -1.1583
35    ]';
36metric_phi = inf;
37tol_phi = 1e-4;
38update_speed = 0.5;
39iter = 0;
40while metric_phi>tol_phi
41    options.phi = phi;
42    options.WarmUp = IterRslt;
43    IterRslt = iter_ks1998(options);
44    SimuRslt = simulate_ks1998(IterRslt,simuOptions);
45
46    % Collect samples and run regression
47    burnPeriods = 1000;
48    periodIndex = 1:num_periods;
49    r2_z = zeros(1,2);
50    rmse_z = zeros(1,2);
51    phi_new = phi;
52    for i_z=1:2
53        periodSet = periodIndex(z_idx==i_z & periodIndex>burnPeriods);
54        sample_lbar = SimuRslt.lbar(1,periodSet);
55        sample_kbar = SimuRslt.kbar(1,periodSet);
56        sample_kbar_next = SimuRslt.kbar(1,periodSet+1);
57        mdl = fitlm(log(sample_kbar),log(sample_kbar_next));
58
59        phi_new(1,i_z) = mdl.Coefficients.Estimate(2);
60        phi_new(2,i_z) = mdl.Coefficients.Estimate(1);
61        phi_new(4,i_z) = log(mean(sample_lbar));
62        r2_z(i_z) = mdl.Rsquared.Ordinary;
63        rmse_z(i_z) = mdl.RMSE;
64    end
65
66    iter = iter+1;
67    metric_phi = max(abs(phi(:) - phi_new(:)));
68    fprintf('iter, metric_phi: %d, %g\n',iter, metric_phi);
69    % Update
70    phi = phi_new*update_speed + phi*(1-update_speed);
71end
72
73%% Post analysis
74disp('Converged transition coefficients (phi_k_1 phi_k_0 phi_l_1 phi_l_0):');
75disp(phi')
76% fraction negative wealth
77sum(SimuRslt.k(:,num_periods)<0) / numel(SimuRslt.k(:,num_periods))
78% Gini coefficients
79samplePeriods = num_periods;
80k = SimuRslt.k(:,samplePeriods);
81% ginicoeff(k(k>=0)) % Uncomment this to calculate gini if you have ginicoeff routine
82
83%%%%%%%%%%%%%%% MAIN ENDS HERE %%%%%%%%%%%%%%%%%%
84function shock = gen_discrete_markov_rn_with_agg(trans,numPaths,lenPath,initShock,aggIdx)
85% Generate idiosyncratic shocks conditional on aggregate shocks
86numStates = size(trans,1);
87shock = zeros(numPaths,lenPath);
88shock(:,1) = initShock;
89cumTrans = cumsum(trans, 2);
90for t=1:lenPath-1
91    for j=1:numStates
92        % Find samples that has shock j
93        idxOfShockJ = find(shock(:,t)==j);
94        % Fill in numbers based on distribution weights
95        un = rand(length(idxOfShockJ),1);
96        % Look up trans_j
97        [~,shock(idxOfShockJ,t+1)] = histc(un, [0 cumTrans(j,:,aggIdx(t),aggIdx(t+1))]);
98    end
99end
100end