# 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 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 % 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

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 %% 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