Awesome state Aiyagari

1

I coded up an example based on an Aiyagari model with an awesome state for the labor efficiency shock. This is in turn inspired by:

Appendix E.3 of Fatih Guvenen, Gueorgui Kambourov, Burhan Kuruscu, Sergio Ocampo, Daphne Chen, Use It or Lose It: Efficiency and Redistributional Effects of Wealth Taxation, The Quarterly Journal of Economics, Volume 138, Issue 2, May 2023, Pages 835–894, Use It or Lose It: Efficiency and Redistributional Effects of Wealth Taxation* | The Quarterly Journal of Economics | Oxford Academic

The purpose of this example was to solve a model with stochastic ageing with the toolkit. The model has a decision variable d, labor supply, and the state variables are:
a: assets, endogenous
e: labor efficiency, exogenous
age={young,retired}, exogenous

The default/standard way of doing this is to set up the exogenous state z as the combo of (e,age) as follows:

n_z    = [n_e,n_age];
z_grid = [e_grid;age_grid];

pi_z = [(1-Params.prob_ret)*pi_e,          Params.prob_ret*eye(n_e);
        Params.prob_death*ones(n_e,1)*G_e', (1-Params.prob_death)*eye(n_e)];

where pi_e is a [n_e*n_e] transition matrix and pi_z is a [2*n_e,2*n_e] transition matrix.
This approach however is “dumb” because it wastes grid points: if age=retired, there is no need to keep track of labor productivity e.

Another approach (thanks for the suggestion @robertdkirkby!) is to use correlated shocks and set up the z variable as:

n_z    = [n_e+1,1];
z_grid = [e_grid,age_grid(1)*ones(n_e,1)
          0,     age_grid(2)];

pi_z = [(1-Params.prob_ret)*pi_e, Params.prob_ret*ones(n_e,1);
        Params.prob_death*G_e',   (1-Params.prob_death)];

Note that I put a 0 for e when age=2 but I could have written any value, since it is not used. Now z_grid is a matrix of size [n_e+1,2] and pi_z is a matrix of size [n_e+1,n_e+1]

I compared these two methods and, as expected, the code with correlated shock is faster. However there are small differences in the results. I tracked down the source of the discrepancy to the distribution (the policy functions are identical)

Compare method 1 to method 2
||Policy_a_m1-Policy_a_m2||
     0

||StationaryDist_m1-StationaryDist_m2||
       0.00307548793025022

My codes are available here: GitHub - aledinola/Awesome-state-Aiyagari

2

Awesome!!!

Together with the Aiyagari with taxes, and Pijona-Mas with taxes it is a nice series of models :smiley:

I think the discrepancy in the distrubution was just a typo in line 99 of main.m, namely where it says StationaryDist(:,1,2) it should instead be sum(StationaryDist(:,:,2),2). Correcting this the discrepancy in the stationary distribution is of the order of 10^(-6). Given this is the max difference between any two points on the whole grid I suspect it is just rounding error (I may be wrong). I tried summing over assets to look at just the distribution over exogenous states,

temp1=sum(StationaryDist_m1,1);
temp2=sum(StationaryDist_m2,1);
max(abs(temp1-temp2))

and this was 10^(-7).

Visually, the cdfs are identical.

1 Like
3

Hi Robert,

Many thanks for checking this. After correcting the typo in main.m, I get the same difference, around 10^(-6), which is also the tolerance for the convergence of the distribution.

I noticed however that this discrepancy generates some differences in a few moments, e.g. the share of wealth owned by the top 1 percent. To make the two methods deliver the same results for the model statistics (“same” here means that the first three decimals are the same, what you would typically show in a paper) you need to set

simoptions.tolerance = 1e-9;

The default value of 1e-6 might be too loose if the focus is on concentration of wealth/income at the top; on the other hand, aggregate variables such as savings or labor are not very sensitive to this.

I would keep 1e-6 as default, and good that there’s an option to change it.

1 Like
4

As a follow up on models with stochastic ageing, I’m trying out a model with the same structure as above, but instead of one labor productivity shock, I have two shocks, called eps and theta (both Markov).

The setup for exogenous z variables is then the following:

  • z = (eps,theta,age)
  • eps, theta: Two markov shocks, I want the full kron product (all combos)
  • age: young or old. If age=young, want all combos of eps,theta. If age=old,
    eps and theta are irrelevant (so can use just one arbitrary value)

In total, the relevant number of points for z should be n_eps*n_theta+1. In order to use the correlated shocks feature of the toolkit, I did the following:

n_z = [n_eps*n_theta+1,1,1];
z_grid = zeros(n_eps*n_theta+1,3);
z_grid(1:n_eps*n_theta,:) = [kron(ones(n_theta,1),eps_grid),kron(theta_grid,ones(n_eps,1)),age_grid(1)*ones(n_eps*n_theta,1)];
z_grid(n_eps*n_theta+1,:) = [0,0,2]; %when old

pi_z = [(1-pr_ret)*kron(pi_theta,pi_eps),pr_ret*ones(n_eps*n_theta,1)
        pr_die*p_star',1-pr_die];

As an illustration, if I have

eps_grid   = [1,2,3]';
theta_grid = [0.5,1.5]';
age_grid   = [1,2]'; %1=young, 2=old

my z_grid looks like this:

1	0.5	1
2	0.5	1
3	0.5	1
1	1.5	1
2	1.5	1
3	1.5 1
0	0	2

The first column of z_grid is the repetition of eps_grid, the second column is theta_grid and the last column is age_grid. The last row of z_grid is when age=old.

Just wanted to double check if this makes sense or I’m off track :slight_smile:
Thanks!

For completeness, I copy and paste here a minimum working example:

clear,clc,close all
format short 

% z = (eps,theta,age)
% eps, theta: Two markov shocks, I want the full kron product (all
% combos)
% age: young or old. If age=young, want all combos of eps,theta. If
% age=old, eps and theta are irrelevant (so can use just one arbitrary
% value)

%% Set up fake data
n_eps      = 3;
n_theta    = 2;
eps_grid   = [1,2,3]';
theta_grid = [0.5,1.5]';
pi_eps     = rand(n_eps,n_eps);
pi_eps     = pi_eps./sum(pi_eps,2);
pi_theta   = rand(n_theta,n_theta);
pi_theta   = pi_theta./sum(pi_theta,2);

% Unconditional distribution of eps and theta
aux = pi_eps^1000;
p_star_eps = aux(1,:)';
aux = pi_theta^1000;
p_star_theta = aux(1,:)';
p_star = kron(p_star_theta,p_star_eps);

age_grid   = [1,2]'; %1=young, 2=old
pr_ret     = 0.2;
pr_die     = 0.05;
pi_age     = [1-pr_ret, pr_ret
              pr_die, 1-pr_die];

%% Move to toolkit notation, correlated shocks

n_z = [n_eps*n_theta+1,1,1];
z_grid = zeros(n_eps*n_theta+1,3);
z_grid(1:n_eps*n_theta,:) = [kron(ones(n_theta,1),eps_grid),kron(theta_grid,ones(n_eps,1)),age_grid(1)*ones(n_eps*n_theta,1)];
z_grid(n_eps*n_theta+1,:) = [0,0,age_grid(2)]; %when old

pi_z = [(1-pr_ret)*kron(pi_theta,pi_eps),pr_ret*ones(n_eps*n_theta,1)
        pr_die*p_star',1-pr_die];

UPDATE: replaced hard-coded 2 with age_grid(2)

1 Like
5

This looks like correctly setup joint-grids.

Only minor thing I see is that you hardcode

z_grid(n_eps*n_theta+1,:) = [0,0,2]; %when old

rather than just

z_grid(n_eps*n_theta+1,:) = [0,0,age_grid(2)]; %when old

(I see in the other part of z_grid you used age_grid(1))

1 Like
6

Thanks! Also for spotting the age_grid(2)