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=1 % =1 semiz_grid is stacked column, =0 semiz_grid is joint-grid
% Four approaches to the shape of semi_grid
semizbyptype=2 % =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.05; % 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;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_semiz1A=[0,1; 0,1]; % for d=0
pi_semiz2A=[1,0; 1,0]; % 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_semiz1B=[1,0; 1,0]; % for d=0
% pi_semiz2B=[0,1; 0,1]; % 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-Wz-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);
temp1=squeeze(sum(StationaryDist.riskyhh,1)) % riskyhh semi-exo shocks have very specific distribution following the ReturnFn forcing d=1 every even period, and the pi_semiz1A and pi_semiz1B [when using matrix setup for pi_semiz]
%%
FnsToEvaluate.Assets=@(d,aprime,a,semiz) a;
FnsToEvaluate.L=@(d,aprime,a,semiz) semiz;
FnsToEvaluate.decision=@(d,aprime,a,semiz) d;
simoptions.conditionalrestrictions.positivea=@(d,aprime,a,semiz) (a>0);
simoptions.conditionalrestrictions.alla=@(d,aprime,a,semiz) (a>=0);
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
%%
AllStats.Assets.Mean
%%
AgeConditionalStats=LifeCycleProfiles_FHorz_Case1_PType(StationaryDist,Policy,FnsToEvaluate,Params,n_d,n_a,n_z,N_j,Names_i,d_grid,a_grid,z_grid,simoptions);
AgeConditionalStats.decision.Mean % 0 at odd ages, 1 at even ages, following the ReturnFn
AgeConditionalStats.L.riskyhh.Mean
AgeConditionalStats.L.safehh.Mean
AgeConditionalStats.L.midhh.Mean
AgeConditionalStats.L.Mean
%%
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);
%%
SimPanelValues=SimPanelValues_FHorz_Case1_PType(jequaloneDist,PTypeDistParamNames,Policy,FnsToEvaluate,Params,n_d,n_a,n_z,N_j,Names_i,d_grid,a_grid,z_grid,pi_z, simoptions);
mean(SimPanelValues.decision,2) % same as AgeConditionalStats.decision.Mean
mean(SimPanelValues.L,2) % same as AgeConditionalStats.L.Mean
%%
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”.
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
and here is the content of “FiniteHorzStochConsSavings2_ReturnFn”, which has unusual content to force the choice of d each period so can easily see if semi-exo shock seems to be working.
function F=FiniteHorzStochConsSavings2_ReturnFn(d,aprime,a,z,gamma,r,Wj,w,agej)
F=-Inf;
%jj: age (index 1:J; indicates ages 1 to 10)
W=Wj*z;
c=(1+r)*a+w*W-aprime;
if c>0
F=(c.^(1-gamma)-1)/(1-gamma);
end
if rem(agej,2)==0 % when agej is even
if d==0
F=-Inf; % choose d=1 when agej is even
end
elseif rem(agej,2)==1 % when agej is odd
if d==1
F=-Inf; % choose d=0 when agej is odd
end
end
end