I’m interested in macro-health literature and would like an assessment on whther the VFI-toolkit is able to solve the quantitative models presented in these papers:
Note: by “able to solve” I mean if the toolkit can solve these models as it is now or with some minor modifications that I could program myself. Just basic methods, not grid interpolation.
Hmm,
Kitao (2014) - A life-cycle model of unemployment and disability insurance
Almost, can handle everything but e which is average past earnings, problem is that this is an experienceasset, but that it also depends on another experienceasset g. Toolkit doesn’t yet have a way to handle ‘pensions that are indexed’ in models that also have ‘human capital’, which is what is going on here. Might be possible to implement very easily using a 2-dimensional experienceasset, but I have just never attempted to do this, probably easy enough as an extension of the codes, just needs to be written.
Model has exogenous markov for health, h, and then health expenditures are drawn from a distribution Pi^M(m;j,h) which is i.i.d (condition on j,h) so you can do this with i.i.d. e variable in toolkit, and use ReturnFn (or FnsToEvaluate) to combine the e with j and h to get the medical expense.
Model has one standard endogenous state, assets a. Model is FHorz, and actually has only 4 periods.
There is human capital g, catch is that this depends on employment, and employment is semi-exogenous (there is search effort s which is 0 or 1 and a probability of finding a job). Really this wants an experienceassetsemiz but that doesn’t exist yet (I plan to do experienceassetz at some point, and there are many variations like this). But you can use a trick to get around this, you put in a decision variable d that is 0 or 1, with 1 being ‘employed’ and use the ReturnFn so that household chooses d=1 when employed and d=0 when unemployed (so d just mirrors the semi-exogenous state, by giving -Inf otherwise), this way we can now do g as a standard experienceasset. [Actually, it is a experienceassetu, as it is based on probabilities, which are given in Appendix A, but this doesn’t make much difference.]
Past earnings e, eqn (5) says e'=f(e,gw). Toolkit cannot handle this because it is an experienceasset but it depends on g which is also itself an experienceasset. In principle you could do this as a 2-dimensional experienceasset but I don’t think that works currently in the code (it might do, but I don’t think it does).
d_U is unemployment duration, I expect you can just roll this into the semi-exogenous employment state (e.g., instead of 0 being unemployed, 1 being employed, you could have 0 being first period of unemployment, 1 being second period of unemployment, 2 being employed; have not checked how many periods of unemployment is the max you need to track in the model but guessing not many as there are only four ages in the model anyway and two of them are retirement).
There are also two binaries i_u and i_m that track eligibility for unemployment benefits and medicare benefits, I haven’t looked into exactly how these are determined, but guessing they are easy enough?
I thought this could be easily handled by declaring both human capital g and employment \ell as semi-exogenous. g is restricted to a grid and the transition matrices for g, conditional on employment, are shown in the appendix.
I assumed that the max on accept/reject or stay/quit in the Bellman equations would be the real problem for the toolkit.
Good point. That is a cleaner way to do it.
Also means that just having a experienceassetsemiz would let you do earnings e. Not currently a feature of the toolkit, but is on my vague to-do list somewhere a few bullet points down from the top. (I want experienceassetz for indexed pensions in model the markov, and the experienceassetze and experienceassetsemiz are obvious things you also want for same purpose in models with more shocks)
Oops, I’d totally missed that in my rather quick read, was just looking at the state-space. I suspect that you can still do this, just that you add a decision variable stay/quit and interact it with everything else in the ReturnFn. You also need an i.i.d shock for stay/quit, and if you get the exogenous separation shock then the stay/quit decision variable does nothing. It will be a bit messy in the ReturnFn, but it appears like it should work okay. (Same thing for doing accept/reject, you can reuse the same decision variable since you never have both.) Actually, since these accept/reject and stay/quit influence the employment status the binary decision variable for them has to be part of influencing the semi-exogenous state, but that should be okay (the i.i.d. separation shock should then be rolled into the semi-exo state too).
Also, do you need the accept/reject decision, aren’t all jobs identical so why are you exerting search effort if you would reject the job? Or did I miss something that makes jobs differ in the model?
One quick question: Can you input into the toolkit the transition matrix for semi-exogenous shocks as a 4-dimensional (semiz,semiz',d,j) array? You said recently that you can’t do it for permanent type but without permanent types? In the Kitao (2014) model that would be nice because g has many points and inputting the transition probabilities for g, \ \ell in a function with if conditions is impractical.
The answer is: yes you can input directly the array.
I found the answer to my own question by reading this code:
As far as I understand, the user can directly input
vfoptions.pi_semiz_J
which is quite handy. But then, it shouldn’t be a problem to input
vfoptions.pi_semiz_J.type1
And
vfoptions.pi_semiz_J.type2
Correct?
The only missing feature that would help me is the possibility to input directly the type as 5-th dimension in pi_semiz_J.
Regarding the accept/reject and stay/quit decisions:
I think it depends on the model’s assumptions.
There are several models in the macro literature that employ this feature: for example, a series of papers by Krusell, Mukoyama, Rogerson, Sahin and coauthors:
The one linked above is cast in infinite-horizon but there are also life-cycle versions with a similar structure (see for example Gross Worker Flows over the Life Cycle). It would be nice to have an example that explains how to model this feature with the toolkit.
De Nardi, Pashchenko & Porapakkarm (2025) - The Lifetime Costs of Bad Health
Short answer: They have two sub-periods, in the first sub-period you buy health insurance and choose labor supply, in the second you do everything else. Toolkit doesn’t really have ability to handle subperiods (conceptually would be easy to code, but has not been done). You could however lump everything together in one period, and just make is so that instead of choosing insurance purchase this period, you switch it to an endogenous state you chose last period, and you would get much the same model. You then get something toolkit can solve easy enough. [Not really clear in model why you would want labor supply chosen without knowing medical expense shocks, seems more plausible that people choose their current work amount based on their current health? Obvs this is an empirical question.]
Section 4.6 on page 211: The state space for the first sub-period is:
(a,h,tau,upsilon,g;j,p1,p2,p3)
where a is assets a standard endogenous state, h is health status a 3-valued markov, tau is current-health duration a 3-valued markov, upsilon is markov (AR(1) on labor productivity), g is binary valued for ‘insurance offer’ and is a markov (explained below).
j is age, p1,p2,p3 are three dimensions of permanent type (explained below, there are 3x2x3=18 ptypes)
While state-space for the second sub-period also includes:
(i_H, l, x^H)
where i_H is four-valued (essentially a standard endogenous state), l is binary labor supply decision, x^H is an medical expense shock that depends on age and health, so essentially it is an extra dimension to the markov on (h,tau).
Retirees have a simplified version of this.
Health-status, h, is a (permanent-type dependent) markov process. There are three health levels (poor, bad, fair) but you also have to allow the transition probabilities to depend on ‘length of time in current health state’ and they allow up to three periods for this [this is my reading of eqns (2), (4) and (7) together with first paragraph in section 2.4.2, where they say “In our life-cycle model, we use our estimated health process with T = 3, which means 4 years of lagged dependence.”]. They implement this as markov h with 3 states (poor,fair,good health) and a markov tau with 3 states (1,2,3 years in current health). This part is easy enough (although the exact equations to create them are quite complex, they are easy enough from a compute perspective).
g is a markov for ‘can buy insurance’. It depends on age (j), labor producitivity (upsilon) and health (h). This is trivial as just means treating all the markovs (g,upsilon,h,tau) as one big age-dependent markov, but will take a lot care to set it up correctly.
Medical expense shocks x^H are a markov that depends on age, and health. So yet another dimension on our big markov which is now (g,upsilon,h,tau,x^h).
Model period is two-years.
Permanent-type has three dimensions, health type, discount factor and labour productivity, and these are correlated (this is easy for toolkit). They have 3x2x3=18 permanent types total.
Labor supply is 0 or 1.
Earnings has deterministic age component, and AR(1) component, and a permanent type. All easy.
There is standard consumption-savings decision.
Survival probabilities depend on age and health, this won’t be trivial but is do-able (add a state to joint-markov on (h,tau) that is ‘dead’, then use simoptions.conditionalrestrictions so it is omitted from all the model stats).
Summing up: This is a very complicated version of a very simple model. You have one standard endogenous state (two if you count i_H, but the second is just 4-values), and then a big age-dependent markov (g,upsilon,h,tau,x^H), and one decision variable (l), together with permanent types (p1,p2,p3). If it wasn’t for the two-subperiods the toolkit would solve this easily enough, as it is you could ‘merge’ the subperiods and get a very similar model the toolkit would solve just fine (by making i_H something chosen previous period so it becomes a second endogenous state, and allowing choice of l to depend on medical expense shock so it becomes a decision variable). That said, the complexity of setting up the markov correctly should not be underestimated, in principle it is trivial (you are just creating a grid and transition matrix for a discrete first-order markov process), but in practice it involves keeping track of a large number of equations and writing everything out very carefully.
PS. I guess the issue with merging sub-periods and making i_H an endogenous state is that it then depends on previous period upsilon and h rather than this period. But seems unlikely to be an important difference. Not a priori obvious that sub-periods is better/worse than merging the periods. Empirical fit is only thing that would distinguish the two and hard to know whether that would be better or worse.
Tha thanks for the very detailed explanation!
I agree that if it wasn’t for the two sub-periods structure it would be a relatively simple model: one endogenous state with an age-dependent Markov process and permanent types. The Markov states are not semi-exogenous but purely exogenous. The Markov states are not independent so you need care in setting up the transition matrix. But compitationally is trivial right?
Say you have two exogenous Markov z1 and z2 and pi1(z1,z1’) but pi2(z2,z2’|z1), so that the transitions of z2 depend on current value of z1. In the toolkit you set
n_z = [n_z1,n_z2];
z_grid = [z1;z2];
pi_z4 = zeros(n_z1,n_z1,n_z2,n_z2p);
for z1=1:n_z1
for z2=1:n_z2
for z1p=1:n_z1
for z2p=1:n_z2
pi_z(z1,z1p,z2,z2p) = pi1(z1,z1p)*pi2(z2,z2p|z1);
end
end
end
end
pi_z = reshape(pi_z4,prod(n_z),prod(n_z));
I was thinking about
max_{x,y} f(x,y)
versus
max_{x} g(x); where g(x)\equiv max_y f(x,y)
for CPU vs GPU. So merged period vs two sub-periods.
By doing two subperiods instead of one, you eliminate some of the ability to interact the x and y (not in the above simple equations as long as all the functions are bijections, but yes in model because there is a shock that is only known in between the two sub periods, versus known prior to the merged period). So this has a computational gain. On the CPU this is the end of the story. But on the GPU, we went from parallel to serial which comes with a computational cost. So I kind of suspect that the reason their paper does subperiods is just to cut computation, whereas on the GPU this is instead computationally costly because you loose parallelization (as long as you don’t run out of memory solving the merged period problem).
PS. Not really ‘suspect’, I am highly confident even without asking them that they used subperiods to cut computational runtimes. They are a competent and experienced bunch of authors and I would be very surprised if this wasn’t deliberate.
Interesting observation!
I have another partially unrelated question: when I input age-dependent parameters, does the toolkit check the dimension internally? For example, if I set params.ej, does the toolkit check that ej has N_j elements? I am not sure if it should, just asking.
It could be that ej is relevant only during working age, so it should have only Jr-1 elements. I always set the remaining elements equal to 0.
Regarding pi_z_J: does the toolkit check that all elements are positive and that each row sums to 1 by age?
Not only does the toolkit check the dimensions of parameters and infer that anything by-N_j or N_j-by is age-dependent, but in the case of Transition Paths, it pays attention to by-N_T and N_T-by. Which means you want your ages and transition periods to be distinct in length, lest it transpose or fail to transpose your parameters when seeing square inputs. On the flipside, it’s not the name of the parameter that makes the toolkit assume it’s age-like (other than N_j, which is defined by the toolkit and which you use as you need it). It knows nothing about ej other than what you tell it.
You can set an unused age to NaN, just to ensure that it’s really not used.
I haven’t looked into the pi_z_J question, but there are a number of places where summing to one is checked (age weights, PType masses).
As Michael says toolkit essentially checks if any dimension of the parameter is N_j (or N_i or T when using ptypes and TransPath, respectively). It then parses all the values of that parameter appropriately to extract the relevant j (or i or t) and passes that (now scalar) value to any functions.
As long as you don’t tell the ReturnFn to use ej for the retirement ages, while the toolkit will internally parse them and hand them to ReturnFn, and they then just get ignored by the code in your ReturnFn. So, e.g., if you ReturnFn contains
if agej<Jr
c=(1+r)*a+w*kappa_j*l*z-aprime;
else
c=(1+r)*a + pension;
end
Then while the toolkit is parsing kappa_j for all the different j=1,…,N_j, the ReturnFn is just ignoring all the values for j>=Jr. Putting zeros or NaN so that it would be very obvious to you as user if you accidently write code that uses them is good practice.
For pi_z (and pi_z_J, pi_e, pi_e_J):
There is a check in ValueFnIter for InfHorz that the elements are all positive and that rows sum to one (to be precise, it checks rows summing to within 1e-13 of one, otherwise numerical rounding errors caused it problems).
There is not currently a check in ValueFnIter for FHorz. But there should be, I will try add this same check sometime today.
[The check is done only in ValueFnIter on the grounds you probably always run this as the first thing you do after creating the pi_z. And so the other commands just assume that if it worked for ValueFnIter then it presumably still works for all the others.]
PS. Michael rightly points out if you ever write a model where N_j, N_i or T are the same lengths, you would be relying on complete fluke for the toolkit to do the correct thing. DO NOT do this. [I figure it is crazy low odds that anyone ever does, but I should probably just add error messages for this.]
Now has tests for pi_z and pi_e (elements are positive, things sum to one when they should).
There are already tests for this with semi-exogenous shocks. I changed them from error to warning in a subsequent push [maybe someone wants to set pi_z that doesn’t sum to one for some reason, not obvious why but I can see it being some kind of thing people wanted to play with for some weird behavioural/preference thingy].