Semi-exogenous shocks with permanent types

1

I might have asked already this question some time ago but I tend to forget :slight_smile:

Suppose I have a semi-exogenous shock that is age and type dependent, so the transition matrix is

\pi(semiz,semiz'| d,j,\theta)

where d is the decision variable affecting the transition probabilities of semiz, j denotes age and \theta is the permanent type.

I find it a bit odd to input the transition probabilities using the function SemiExoStateFn. If semiz and/or \theta have lots of values it is much easier to set up the transition as a matrix (which is also what the code does internally).

So my first question is whether I can do that, i.e. I can set pi_semiz_J as an array and pass it to the toolkit.

Regarding the permanent types, If I have more than 2 or 3 values, I find it odd to have to type pi.e1, pi.e2, pi.e3 etc. What I would like to do is this:

Define pi_semiz_J as a multidimensional array with size [n_semiz,n_semiz,N_d2,N_j,N_ptype]. Can I do this?

Apologies if I have asked this question already, but I haven’t been able to find it here in the discourse. By the way, is there an “official” example in the VFI toolkit that combines these three features:

  1. Semi-exogenous shocks
  2. Transition of semi-exogenous shocks depends also on age
  3. Transition of semi-exogenous shocks depends also on permanent type
2

Looking at the code, there’s support for PTypes in ExoShocks, but not SemiExoShocks. So SubCodes/SemiExoShocks/SemiExogShockSetup_FHorz_PType.m would need to be written.

3

I found an example that has two features (shock that depends on age and fixed type, but not on decision d) in the replication of De Nardi, French and Joines:

The code below shows how to set up a transition matrix for exogenous Markov state h that depends on age and permanent types. The matrix has 3 dimensions: (h,h’,age).

The cumbersome part is that you have to define a different object for each permanent type and there are 10 of them: each type is a separate field of pi_h_J. pi_h_J is a structure, not a numeric array.

I would like to define pi_h_J as a numeric array with 4 dimensions: (h,h’,age,type). Note that type is all combinations of types: so here there are two genders and 5 income groups so in total 10 fixed types.

Question: Can I do this and pass pi_z_J as a 4 dim array to the toolkit?

Question: Can I define pi_semiz_J as a 5 dim array (h,h’,d,age,type) and pass it to the toolkit?

If not, we have to write a lot of boilerplate code. Nowadays is less of an issue given AI tools bit still…

pi_h_J.femaleq1=zeros(n_z(1),n_z(1),N_j);
g_c=1; % female
q_c=1; % 1st income quintile
for hlag_c=1:2
    for jj=1:N_j
        % conditional survival probabilities
        sj=probalive_1yr(hlag_c,jj,g_c,q_c); % current h, age+1, gender, income q
        % health transitions
        pi_h_J.femaleq1(hlag_c,3,jj)=(1-sj); % probability of dying
        for h_c=1:2
            pi_h_J.femaleq1(hlag_c,h_c,jj)=sj*probhealth(hlag_c,h_c,jj,g_c,q_c); % prob of survival*prod of health status
        end
    end
end

The code above shows just one matrix for g=1 and q=1. The user has to copy and paste this block other nine times changing the field names.

4

This is not an authoritative answer, but generally I think that the toolkit expects pi_z_J things to be age-like in the final dimension. When working with PTypes, they are special until they get unpacked, usually into something like:

pi_z_J_this_PType = PTypestructure.(iistr).pi_z_J

where iistr resolves to the name of the PType (usually Names_i{ii}). Thus, from what I’ve seen, PType is not generally a dimension per se, but analogous in a way to the classic split/apply/combine workflow (see Group by: split-apply-combine — pandas 3.0.1 documentation and the following diagram).

image
image1610×857 91 KB

We split the model into PTypes when we construct it, we “apply” the construction of the initial values, policies, and stationary distributions for each PType, and then we “combine” by solving the general equilibria across the various PTypes that all interact with one another.

Of course data science and economics are two different fields, and the tidyverse is not exactly LifeCycle modeling (nor OLG nor Transition Paths). But I think the analogy illustrates how, based on my reading and experience with the toolkit, PTypes are not (yet) dimensional in their nature.

1 Like
5

Thanks for the interesting answer! I agree with you as far as the toolkit outputs are concerned. For example, with permanent types the value function is a structure whose fields are multidimensional arrays:

V.ptype1 is [n_a,n_z,N_j] array
V.ptype2 is [n_a,n_z,N_j] array
Etc.

My point is that if I have a parameter, let’s call it theta, that depends both on age and the permanent type, I would like to input it as a [N_j,n_ptype] matrix instead of

theta.ptype1 vector with N_j elements
theta.ptype2 vector with N_j elements
Etc.

Why? Because the method with field names becomes cumbersome if you have more than 2 or 3 values.

There is actually an example in the intro to lifecycle models that allows the user to input the permanent type as a dimension in an array (see model 24, parameter alpha_i)

6

Very cool! And now that you mention it I “see” the use of a bunch of code within the toolkit that I’ve normally just breezed by.

I suggest that alphadist on line 62 of the Lifecycle example should be alphadist_i.

7

What @MichaelTiemann says above about how things get converted internally (with iistr in a structure, and age as a dimension) is what happens internally to any ‘shape’ of input. This is so all the internal code ‘knows’ what shape the shocks are, but doesn’t tell us about what ‘shapes’ the user is able to input.

I believe the current answer is no. But the answer should be yes. I will see if I can find time next week to make the answer yes (I have a bunch of stuff around calibration and custommodelstats that I have promised a co-author to implement next week, so the ptypes for shocks might get bumped to the following week).

PS. As @aledinola notes. You can already implement parameters that depend on ptype, either by inputting Params.alpha as a 3-by-1 vector (if there are 3 ptypes; or 1-by-3 vector), or by inputting the parameter itself as a structure with fieldnames based on Names_i so something like Params.alpha.ptype1 Params.alpha.ptype2 and Params.alpha.ptype3. This is the intended approach of the toolkit, that you can choose which one suits you as user when setting up the model.

PPS. The “split-apply-combine” is a good description of how most of the ptype code is implemented internally. The main design reason for this choice was simply that it allows most of the “apply” step to just reuse code from the ‘without-ptype’ commands, meaning that implementing ptype commands just requires writing all the “split” and “combine” code.

2 Likes
8

Nice explanation, thanks! Let me know when you introduce this improvement

9

Now works for markov z. You can set up ptype either as fields of a structure, or just as an additional dimension. [I did the main commands, might be I missed a command so if one does error let me know]

Following is code that demos this. Near top is “zbyptype” which you can set =1,2,3 or 4, to see the different setups. Also “zstackcolumn” which you can set to =0 or =1, and which is just different ways to think of z_grid.

% Age
% One endogenous variable: assets
% One stochastic exogenous variables: income shock

J=10; % Ages 1 to 10 inclusive.

% Grid sizes to use
n_d=0;
n_a=501;
n_z=5; % income
N_j=J; % Number of periods in finite horizon

Names_i={'riskyhh','safehh'};

zbyptype=1 % =1 is ptype structure, =2 is ptype as dimension
           % =3 is _J with ptype as structure, and =4 is _J with ptype as dimension
zstackcolumn=0 % =1 z_grid is stacked column, =0 z_grid is joint-grid

%% Declare the model parameters

Params.gamma=2;
% Gamma plays three roles:      
        % gamma is coefficient of relative risk aversion
        % 1/gamma is the intertemporal elasticity of substitution of consumption
        % gamma+1 is the coefficient of relative prudence (Kimball, 1990)

Params.beta=0.96; % rate at which agents discount future
Params.r=0.03; % interest rate on assets
Params.w=1.6;

% Declare the age dependent parameters. This is a simple matter of creating
% the parameter as a row vector of length J (the VFI Toolkit automatically
% detects which parameters depend on age and which do not).
Params.Wj=[1,2,3,5,7,8,8,5,4,4]; % deterministic income depends on age

% Markov z for first ptype:
Params.rho_z1=0.7;
Params.sigma_z1epsilon=0.05;
[z1_grid,pi_z1]=discretizeAR1_FarmerToda(0,Params.rho_z1,Params.sigma_z1epsilon,n_z); % Farmer-Toda discretizes an AR(1) (is better than Tauchen, & better than Rouwenhorst for rho<0.99)
z1_grid=exp(z1_grid);

% Markov z for first ptype:
Params.rho_z2=0.9;
Params.sigma_z2epsilon=0.03;
[z2_grid,pi_z2]=discretizeAR1_FarmerToda(0,Params.rho_z2,Params.sigma_z2epsilon,n_z); % Farmer-Toda discretizes an AR(1) (is better than Tauchen, & better than Rouwenhorst for rho<0.99)
z2_grid=exp(z2_grid);

if zbyptype==1 % =1 is ptype structure
    z_grid.riskyhh=z1_grid;
    pi_z.riskyhh=pi_z1;
    z_grid.safehh=z2_grid;
    pi_z.safehh=pi_z2;
elseif zbyptype==2 % =2 is ptype as dimension
    if zstackcolumn==1
        z_grid(:,1)=z1_grid;
        z_grid(:,2)=z2_grid;
    elseif zstackcolumn==0
        z_grid(:,:,1)=z1_grid;
        z_grid(:,:,2)=z2_grid;
    end
    pi_z(:,:,1)=pi_z1;
    pi_z(:,:,2)=pi_z2;
elseif zbyptype==3 % =3 _J is ptype structure
    if zstackcolumn==1
        z_grid.riskyhh=z1_grid.*ones(1,J);
        z_grid.safehh=z2_grid.*ones(1,J);
    elseif zstackcolumn==0
        z_grid.riskyhh=z1_grid.*ones(1,1,J);
        z_grid.safehh=z2_grid.*ones(1,1,J);
    end
    pi_z.riskyhh=pi_z1.*ones(1,1,J);
    pi_z.safehh=pi_z2.*ones(1,1,J);
elseif zbyptype==4 % =4 _J is ptype as dimension
    if zstackcolumn==1
        z_grid(:,:,1)=z1_grid.*ones(1,J);
        z_grid(:,:,2)=z2_grid.*ones(1,J);
    elseif zstackcolumn==0
        z_grid(:,:,:,1)=z1_grid.*ones(1,1,J);
        z_grid(:,:,:,2)=z2_grid.*ones(1,1,J);
    end
    pi_z(:,:,:,1)=pi_z1.*ones(1,1,J);
    pi_z(:,:,:,2)=pi_z2.*ones(1,1,J);
end

%% Grids
maxa=150;
a_grid=linspace(0,maxa,n_a)'; % Could probably do better by adding more grid points near zero

d_grid=[];

%% Now, create the return function 
DiscountFactorParamNames={'beta'};

ReturnFn=@(aprime,a,z,gamma,r,Wj,w) FiniteHorzStochConsSavings_ReturnFn(aprime,a,z,gamma,r,Wj,w)
% For the return function the first inputs must be (any decision variables), next period endogenous
% state, this period endogenous state (any exogenous shocks). After that come any parameters.

%% Now solve the value function iteration problem

vfoptions.gridinterplayer=1;
vfoptions.ngridinterp=50;
simoptions.gridinterplayer=vfoptions.gridinterplayer;
simoptions.ngridinterp=vfoptions.ngridinterp;

vfoptions.verbose=0;
tic;
[V, Policy]=ValueFnIter_Case1_FHorz_PType(0,n_a,n_z,N_j,Names_i, 0, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, vfoptions);
toc

% max(max(max(max(Policy))))<n_a % Double check that never try to leave top of asset grid.

% V is assets-by-z-by-age

%%
% Just use jequaloneDist that doesn't differ by ptype, you can have it differ by ptype
jequaloneDist=zeros(n_a,n_z,'gpuArray');
jequaloneDist(1,:)=ones(1,n_z)/n_z;

AgeWeightParamNames={'mewj'};
Params.mewj=ones(1,N_j)/N_j;

PTypeDistParamNames={'pmass'};
Params.pmass=[0.6,0.4];

StationaryDist=StationaryDist_Case1_FHorz_PType(jequaloneDist,AgeWeightParamNames,PTypeDistParamNames,Policy,n_d,n_a,n_z,N_j,Names_i,pi_z,Params,simoptions);

%%
FnsToEvaluate.Assets=@(aprime,a,z) a;
FnsToEvaluate.L=@(aprime,a,z) z;

AllStats=EvalFnOnAgentDist_AllStats_FHorz_Case1_PType(StationaryDist,Policy,FnsToEvaluate,Params,n_d,n_a,n_z,N_j,Names_i,d_grid,a_grid,z_grid,simoptions);

AllStats.L.riskyhh.Mean
AllStats.L.riskyhh.StdDeviation
AllStats.L.safehh.Mean
AllStats.L.safehh.StdDeviation


%%
GEPriceParamNames={'r'};
GeneralEqmEqns.capitalmarket=@(r) 0;
heteroagentoptions=struct();

[p_eqm,GECondns]=HeteroAgentStationaryEqm_Case1_FHorz_PType(n_d, n_a, n_z, N_j, Names_i, 0, pi_z, d_grid, a_grid, z_grid,jequaloneDist, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Params, DiscountFactorParamNames, AgeWeightParamNames, PTypeDistParamNames, GEPriceParamNames,heteroagentoptions, simoptions, vfoptions);


%%
CalibParamNames={'w'};
TargetMoments.AllStats.L.Mean=AllStats.L.Mean;
ParametrizeParamsFn=[];
caliboptions=struct();

[CalibParams,calibsummary]=CalibrateLifeCycleModel_PType(CalibParamNames,TargetMoments,n_d,n_a,n_z,N_j,Names_i,d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, jequaloneDist,AgeWeightParamNames, PTypeDistParamNames, ParametrizeParamsFn, FnsToEvaluate, caliboptions, vfoptions,simoptions);

%%
[CalibParams,calibsummary]=CalibrateOLGModel_PType(CalibParamNames,TargetMoments,n_d,n_a,n_z,N_j,Names_i,d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, jequaloneDist, AgeWeightParamNames, GEPriceParamNames, PTypeDistParamNames, ParametrizeParamsFn, FnsToEvaluate, GeneralEqmEqns, heteroagentoptions, caliboptions, vfoptions,simoptions);

The return function file, FiniteHorzStochConsSavings_ReturnFn.m contains


function F=FiniteHorzStochConsSavings_ReturnFn(aprime,a,z,gamma,r,Wj,w)

%jj: age (index 1:J; indicates ages 1 to 10)

W=Wj*z;

c=(1+r)*a+w*W-aprime;

F=-Inf; %(l by aprime)

if c>0
    F=(c.^(1-gamma)-1)/(1-gamma);
end

end
1 Like
10

Next is semiz. But one is done :slight_smile:

1 Like
11

This is a nice addition! I have two comments:

For both comments, I consider a model with the following state variables:

  • a: endogenous state
  • z: exogenous Markov state
  • j: age
  • Permanent types
  1. In the age-dependent case (let j denote age) with ptype as last dimension, we write pi_z as an array with four dimensions (z,z',j,ptype), which is correct. But z_grid should have only three dimensions: (z,j,ptype). Why does your code define z_grid with 4 dimensions (see case zbyptype==4)? If it is a typo, then why doesn’t the code return an error?
  2. Objects such the policy functions or values on grid are always returned with the ptype structure, even if ptype is used as dimensional input. This is not a complaint, just to make sure I understand. Since sometimes I have many types, I would like the policy function for consumption (just as an example) to be defined as pol_c(a,z,j,ptype) instead of pol_c.pt1(a,z,j), pol_c.pt2(a,z,j), …, pol_c.ptN(a,z,j). What I always do in my codes:
ValuesOnGrid=EvalFnOnAgentDist_ValuesOnGrid_FHorz_Case1_PType(Policy,FnsToEvaluate,Params,n_d,n_a,n_z,N_j,Names_i,d_grid,a_grid,z_grid,simoptions);

pol_c = zeros(n_a,n_z,N_j,numel(Names_i));
for ii=1:numel(Names_i)
    name_i = Names_i{ii};
    pol_c(:,:,:,ii) = ValuesOnGrid.C.(name_i);
end

This is really a minor hassle though and changing this convention would break existing code.

12

There is nothing to “1”. I was thinking of z_grid as joint-grid, you can equally do stacked-column grid. I just updated the code in the above post adding in “zstackedcolumn=1” or =0. =0 is the previous functionality where I was thinking of z_grid as a joint-grid, =1 treats z_grid as a stacked column which is how you were thinking of z_grid. Both work equally well.
[looks a little silly as z is one dimesional so the stacked-grid and joint-grid are actually the same thing]

For “2”. I guess this could be done as something that the user can ‘request’ using vfoptions/simoptions. Say, only for the outputs for the ValueFnIter, StationaryDist, and ValuesOnGrid commands. I’ll think about adding this as is probably easy enough.
[Obviously it will just error in more advanced ptype setups where the different agents have different state-spaces, but the user should anyway be aware of that]

1 Like
13

I’m done with the semi-exogenous states for permanent types.

You must set up vfoptions.semiz_grid, and then one of vfoptions.SemiExoShockFn or vfoptions.pi_semiz. You can include age and ptype as dimensions of semiz_grid and pi_semiz, or not as you wish. Dependence on ptype can be done as a (the last) dimension, or as fieldnames.

Here is the code that I was checking all the combos.

% Age
% One endogenous variable: assets
% One stochastic exogenous variables: income shock

J=10; % Ages 1 to 10 inclusive.

% Grid sizes to use
n_d=2;
n_a=501;
n_semiz=2;
n_z=0; % income
N_j=J; % Number of periods in finite horizon

Names_i={'riskyhh','safehh','midhh'};
% Note: second and third types will actually be identical, but I needed to get Names_i to be different size to n_semiz

semizbymatrix=1 % =0 is SemiExogShockFn, =1 is matrix
semizstackcolumn=0 % =1 semiz_grid is stacked column, =0 semiz_grid is joint-grid

% Four approaches to the shape of semi_grid
semizbyptype=1 % =1 is ptype structure, =2 is ptype as dimension
               % =3 is _J with ptype as structure, and =4 is _J with ptype as dimension
% When semizbymatrix=1, these four approaches also apply to the shape of pi_semiz

if semizbymatrix==0
    % Two approaches to SemiExoShockFn
    semizfnbyptype=1 % =0 different parameter values in SemiExoShockFn
                     % =1 SemiExoShockFn is itself a ptype structure
end 

%% Declare the model parameters

Params.gamma=2;
% Gamma plays three roles:      
        % gamma is coefficient of relative risk aversion
        % 1/gamma is the intertemporal elasticity of substitution of consumption
        % gamma+1 is the coefficient of relative prudence (Kimball, 1990)

Params.beta=0.96; % rate at which agents discount future
Params.r=0.03; % interest rate on assets
Params.w=1.6;

Params.agej=1:1:J;

% Declare the age dependent parameters. This is a simple matter of creating
% the parameter as a row vector of length J (the VFI Toolkit automatically
% detects which parameters depend on age and which do not).
Params.Wj=[1,2,3,5,7,8,8,5,4,4]; % deterministic income depends on age


semiz_grid=[0;1];

pi_semiz1A=[0.2,0.8; 0.8,0.2]; % for d=0
pi_semiz2A=[0.8,0.2; 0.2,0.8]; % for d=1
pi_semizA(:,:,1)=pi_semiz1A;
pi_semizA(:,:,2)=pi_semiz2A;

pi_semiz1B=[0.4,0.6; 0.6,0.4]; % for d=0
pi_semiz2B=[0.6,0.4; 0.4,0.6]; % for d=1
pi_semizB(:,:,1)=pi_semiz1B;
pi_semizB(:,:,2)=pi_semiz2B;

Params.ptypeB=[0,1,1]; % so we can input this to SemiExoStateFn, depends on ptype


%% Semi-Exogenous Shock
vfoptions.n_semiz=n_semiz;

if semizbyptype==1 % =1 is ptype structure
    vfoptions.semiz_grid.riskyhh=semiz_grid;
    vfoptions.semiz_grid.safehh=semiz_grid;
    vfoptions.semiz_grid.midhh=semiz_grid;
elseif semizbyptype==2 % =2 is ptype as dimension
    if semizstackcolumn==1
        vfoptions.semiz_grid(:,1)=semiz_grid;
        vfoptions.semiz_grid(:,2)=semiz_grid;
        vfoptions.semiz_grid(:,3)=semiz_grid;
    elseif semizstackcolumn==0
        vfoptions.semiz_grid(:,:,1)=semiz_grid;
        vfoptions.semiz_grid(:,:,2)=semiz_grid;
        vfoptions.semiz_grid(:,:,3)=semiz_grid;
    end
elseif semizbyptype==3 % =3 _J is ptype structure
    if semizstackcolumn==1
        vfoptions.semiz_grid.riskyhh=semiz_grid.*ones(1,J);
        vfoptions.semiz_grid.safehh=semiz_grid.*ones(1,J);
        vfoptions.semiz_grid.midhh=semiz_grid.*ones(1,J);
    elseif semizstackcolumn==0
        vfoptions.semiz_grid.riskyhh=semiz_grid.*ones(1,1,J);
        vfoptions.semiz_grid.safehh=semiz_grid.*ones(1,1,J);
        vfoptions.semiz_grid.midhh=semiz_grid.*ones(1,1,J);
    end
elseif semizbyptype==4 % =4 _J is ptype as dimension
    if semizstackcolumn==1
        vfoptions.semiz_grid(:,:,1)=semiz_grid.*ones(1,J);
        vfoptions.semiz_grid(:,:,2)=semiz_grid.*ones(1,J);
        vfoptions.semiz_grid(:,:,3)=semiz_grid.*ones(1,J);
    elseif semizstackcolumn==0
        vfoptions.semiz_grid(:,:,:,1)=semiz_grid.*ones(1,1,J);
        vfoptions.semiz_grid(:,:,:,2)=semiz_grid.*ones(1,1,J);
        vfoptions.semiz_grid(:,:,:,3)=semiz_grid.*ones(1,1,J);
    end
end

% We can set up pi_semiz or we can set up SemiExoStateFn
if semizbymatrix==0
    if semizfnbyptype==0
        vfoptions.SemiExoStateFn=@(semiz,semizprime,d,ptypeB) Model_semizbyptype_SemiExoShockFn(semiz,semizprime,d,ptypeB);
    elseif semizfnbyptype==1
        vfoptions.SemiExoStateFn.riskyhh=@(semiz,semizprime,d,ptypeB) Model_semizbyptype_SemiExoShockFn(semiz,semizprime,d,ptypeB);
        vfoptions.SemiExoStateFn.safehh=@(semiz,semizprime,d,ptypeB) Model_semizbyptype_SemiExoShockFn(semiz,semizprime,d,ptypeB);
        vfoptions.SemiExoStateFn.midhh=@(semiz,semizprime,d,ptypeB) Model_semizbyptype_SemiExoShockFn(semiz,semizprime,d,ptypeB);
    end
elseif semizbymatrix==1
    if semizbyptype==1 % =1 is ptype structure
        vfoptions.pi_semiz.riskyhh=pi_semizA;
        vfoptions.pi_semiz.safehh=pi_semizB;
        vfoptions.pi_semiz.midhh=pi_semizB;
    elseif semizbyptype==2 % =2 is ptype as dimension
        vfoptions.pi_semiz(:,:,:,1)=pi_semizA;
        vfoptions.pi_semiz(:,:,:,2)=pi_semizB;
        vfoptions.pi_semiz(:,:,:,3)=pi_semizB;
    elseif semizbyptype==3 % =3 _J is ptype structure
        vfoptions.pi_semiz.riskyhh=pi_semizA.*ones(1,1,1,J);
        vfoptions.pi_semiz.safehh=pi_semizB.*ones(1,1,1,J);
        vfoptions.pi_semiz.midhh=pi_semizB.*ones(1,1,1,J);
    elseif semizbyptype==4 % =4 _J is ptype as dimension
        vfoptions.pi_semiz(:,:,:,:,1)=pi_semizB.*ones(1,1,1,J);
        vfoptions.pi_semiz(:,:,:,:,2)=pi_semizB.*ones(1,1,1,J);
        vfoptions.pi_semiz(:,:,:,:,3)=pi_semizB.*ones(1,1,1,J);
    end
    % pi_semiz can be (semiz, semizprime, d) or (semiz, semizprime, d, j) or (semiz, semizprime, d, j, i)
end

% Put them in simoptions as well
simoptions.n_semiz=vfoptions.n_semiz;
simoptions.semiz_grid=vfoptions.semiz_grid;
if semizbymatrix==0
    simoptions.SemiExoStateFn=vfoptions.SemiExoStateFn;
elseif semizbymatrix==1
    simoptions.pi_semiz=vfoptions.pi_semiz;
end
% simoptions.d_grid=d_grid; % Done below after creating d_grid

%% Grids
maxa=150;
a_grid=linspace(0,maxa,n_a)'; % Could probably do better by adding more grid points near zero

d_grid=[0;1];
simoptions.d_grid=d_grid;

z_grid=[];
pi_z=[];

%% Now, create the return function 
DiscountFactorParamNames={'beta'};

ReturnFn=@(d,aprime,a,z,gamma,r,Wj,w,agej) FiniteHorzStochConsSavings2_ReturnFn(d,aprime,a,z,gamma,r,Wj,w,agej)
% For the return function the first inputs must be (any decision variables), next period endogenous
% state, this period endogenous state (any exogenous shocks). After that come any parameters.

%% Now solve the value function iteration problem

vfoptions.gridinterplayer=1;
vfoptions.ngridinterp=50;
simoptions.gridinterplayer=vfoptions.gridinterplayer;
simoptions.ngridinterp=vfoptions.ngridinterp;

vfoptions.verbose=0;
tic;
[V, Policy]=ValueFnIter_Case1_FHorz_PType(n_d,n_a,n_z,N_j,Names_i, d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, vfoptions);
toc

% max(max(max(max(Policy))))<n_a % Double check that never try to leave top of asset grid.

% V is assets-by-z-by-age

%%
% Just use jequaloneDist that doesn't differ by ptype, you can have it differ by ptype
jequaloneDist=zeros(n_a,n_semiz,'gpuArray');
jequaloneDist(1,:)=ones(1,n_semiz)/n_semiz;

AgeWeightParamNames={'mewj'};
Params.mewj=ones(1,N_j)/N_j;

PTypeDistParamNames={'pmass'};
Params.pmass=[0.6,0.2,0.2];

StationaryDist=StationaryDist_Case1_FHorz_PType(jequaloneDist,AgeWeightParamNames,PTypeDistParamNames,Policy,n_d,n_a,n_z,N_j,Names_i,pi_z,Params,simoptions);

%%
FnsToEvaluate.Assets=@(d,aprime,a,semiz) a;
FnsToEvaluate.L=@(d,aprime,a,semiz) semiz;

AllStats=EvalFnOnAgentDist_AllStats_FHorz_Case1_PType(StationaryDist,Policy,FnsToEvaluate,Params,n_d,n_a,n_z,N_j,Names_i,d_grid,a_grid,z_grid,simoptions);

AllStats.L.riskyhh.Mean
AllStats.L.riskyhh.StdDeviation
AllStats.L.safehh.Mean
AllStats.L.safehh.StdDeviation

AllStats.L.Mean
AllStats.L.StdDeviation


%%
GEPriceParamNames={'r'};
GeneralEqmEqns.capitalmarket=@(r) 0;
heteroagentoptions=struct();

[p_eqm,GECondns]=HeteroAgentStationaryEqm_Case1_FHorz_PType(n_d, n_a, n_z, N_j, Names_i, 0, pi_z, d_grid, a_grid, z_grid,jequaloneDist, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Params, DiscountFactorParamNames, AgeWeightParamNames, PTypeDistParamNames, GEPriceParamNames,heteroagentoptions, simoptions, vfoptions);



%%
CalibParamNames={'w'};
TargetMoments.AllStats.L.Mean=AllStats.L.Mean;
ParametrizeParamsFn=[];
caliboptions=struct();

[CalibParams,calibsummary]=CalibrateLifeCycleModel_PType(CalibParamNames,TargetMoments,n_d,n_a,n_z,N_j,Names_i,d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, jequaloneDist,AgeWeightParamNames, PTypeDistParamNames, ParametrizeParamsFn, FnsToEvaluate, caliboptions, vfoptions,simoptions);


%%
[CalibParams,calibsummary]=CalibrateOLGModel_PType(CalibParamNames,TargetMoments,n_d,n_a,n_z,N_j,Names_i,d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, jequaloneDist, AgeWeightParamNames, GEPriceParamNames, PTypeDistParamNames, ParametrizeParamsFn, FnsToEvaluate, GeneralEqmEqns, heteroagentoptions, caliboptions, vfoptions,simoptions);

and here is the content of “Model_semizbyptype_SemiExoShockFn”. I shared the ReturnFn with the above one used for the markov shocks.

function prob=Model_semizbyptype_SemiExoShockFn(semiz,semizprime,d,ptypeB)

prob=-1; % just a placeholder (one that will cause errors if not overwritten)

if ptypeB==0
    if d==0
        if semiz==0
            if semizprime==0
                prob=0.2;
            elseif semizprime==1
                prob=0.8;
            end
        elseif semiz==1
            if semizprime==0
                prob=0.8;
            elseif semizprime==1
                prob=0.2;
            end
        end
    elseif d==1
        if semiz==0
            if semizprime==0
                prob=0.8;
            elseif semizprime==1
                prob=0.2;
            end
        elseif semiz==1
            if semizprime==0
                prob=0.2;
            elseif semizprime==1
                prob=0.8;
            end
        end
    end

elseif ptypeB==1
    if d==0
        if semiz==0
            if semizprime==0
                prob=0.4;
            elseif semizprime==1
                prob=0.6;
            end
        elseif semiz==1
            if semizprime==0
                prob=0.6;
            elseif semizprime==1
                prob=0.4;
            end
        end
    elseif d==1
        if semiz==0
            if semizprime==0
                prob=0.6;
            elseif semizprime==1
                prob=0.4;
            end
        elseif semiz==1
            if semizprime==0
                prob=0.4;
            elseif semizprime==1
                prob=0.6;
            end
        end
    end

end


end
2 Likes
14

I think that covers everything. Tell me if you think I missed something.

1 Like
15

Fantastic! I tested yesterday the code for exogenous Markov and will test today the semiexo Markov.

One quick question: the Model_semizbyptype_SemiExoShockFn is a new thing? To me it seems that this is the approach followed for example in the replication of the Sommers (2016) fertility paper and already well documented in “Introduction to Life Cycle Models”. The new possibility is that now you can input the transition matrix for the semi-exogenous state directly as an array (semiz, semiz’, d) or (semiz, semiz’, j) or (semiz, semiz, j, i). Please correct me if I misstated anything.

16

One small typo in your post:

V is assets-by-Wz-by-age

Replace it with

V is assets-by-z-by-age

1 Like
17

You are right that the ‘vfoptions.SemiExoShockFn’ approach is not new. What is new is the ‘vfoptions.pi_semiz’ approach, and especially the variety of ‘shapes’ it can take. But I wanted to test and make sure that the new was not breaking the old :slight_smile:

1 Like