Issues with Infinite Horizon Transition

[aledinola]

Dear Rob,

Playing around with toolkit with inifinite horizon transitions, I have found that these two toolkit examples do not work (i.e. they return some error messages):

1. Aiyagari (1994) transition: VFItoolkit-matlab-examples/HeterogeneousAgentModels/Aiyagari1994TransitionPath.m at master · vfitoolkit/VFItoolkit-matlab-examples · GitHub

2. Guerrieri and Lorenzoni (2017) transition: VFItoolkit-matlab-examples/HeterogeneousAgentModels/GuerrieriLorenzoni2017_Example.m at master · vfitoolkit/VFItoolkit-matlab-examples · GitHub

3. Additionally, I wrote a replication of Huggett (JME, 1997) using the toolkit (and also my own CPU code). The toolkit code also returns an error. Huggett_1997/vfi_toolkit at main · aledinola/Huggett_1997 · GitHub

It would be great if you could check these three codes whenever you have time, and verify that the transition runs both with pure discretization and with interpolation (I guess grid interp layer plus divide and conquer is supposed to work). Then I’ll move on to the models with entrepreneurs

@MichaelTiemann has kindly made a pull request on Github to fix some of the problems with the transition, please see here: Pull requests · vfitoolkit/VFIToolkit-matlab · GitHub

[MichaelTiemann]

I am very interested in modeling entrepreneurs vs. rent-seeking alternatives. I started down this rabbit hole trying to learn about brain drain vs. brain gain. While there were plenty of papers to read on the subject, the one that I found most interesting was this one by Murphy et al: The Allocation of Talent: Implications for Growth | NBER. On the lifecycle front, Estimation of a Life-Cycle Model with Human Capital, Labor Supply and Retirement | NBER is the next paper I’m reading. I’m excited to have tools that can help me explore these questions!

[aledinola]

If you are interested in models with entrepreneurs, you might find this replication of Buera and Shin (2013) interesting:

GitHub - aledinola/BueraShin2013_v6

[robertdkirkby]

I think that

Aiyagari transition path: Does not run

Fixed (and uploaded a new cleaner version of Aiyagari1994TransitionPath file, including correcting the comments). I was doing a few changes to infinite horizon transition paths in December, substantial improvements. But all my examples use grid interpolation layer and divide and conquer so I didn’t realize I had some typos when those are not used. Cleaned them up now

likely fixes all three of these. I have not yet tried to run the Huggett (1997), nor done a complete run of the GL2017 (I am going to clean it while I’m there), but I hopefully can do these later today.

[aledinola]

Great, thanks!
About GL2017 example: I think it makes sense to update this example to allow for grid interpolation.

(1) Guerrieri and Lorenzoni use EGM which delivers very accurate results. Grid interpolation would increase the accuracy of the replication relative to the original paper. With pure discretization the policy function for consunption has some wiggles that should not be there.

(2) Good opportunity to test grid interpolation algorithm with transition in the context of a paper that replicates well.

I tried to add myself

vfoptions.gridinterplayer=1;
Pass to simoptions

to the GL2017 example: the value function and the distribution commands work as intended but other commands do not. Probably you have to include simoptions as extra argument when you call some functions, otherwise the code does not know that you are using grid interpolation?

[MichaelTiemann]

See Use vfoptions,simoptions,etc properly in calls by MichaelTiemann · Pull Request #4 · vfitoolkit/VFItoolkit-matlab-examples · GitHub

The caller fixups I made likely indicate similar cleanups may be needed in other examples.

[robertdkirkby]

Merged the pull request. Thanks Michael!

Quick comment: simoptions and vfoptions were traditionally optional inputs that you could simply omit. Because of how the toolkit has changed over the past few years, they have become something you will almost always want to use. Many of these older examples were written back before divide-and-conquer and grid interpolation layer existed (the later has only existed since mid-2025, the former since late 2024). Hence many older examples simply omit these options from the inputs to functions, but would benefit from being updated to use them.

[aledinola]

Fantastic!
In the GL2017 example, this line

StationaryDist_initial=StationaryDist_Case1(Policy_initial,n_d,n_a,n_z,pi_z)

does not show simoptions being passed to the stationary distribution command. Is that ok?

[robertdkirkby]

No. It should get passed.

[aledinola]

I will run tomorrow the examples and report back

Update

I was faster Surprisingly, the line StationaryDist_initial=StationaryDist_Case1(Policy_initial,n_d,n_a,n_z,pi_z); runs without any error! But how does the toolkit know that we are using grid interpolation?

[robertdkirkby]

It doesn’t know, it just ignores the interpolation part and treats the lower grid index as if it was just a standard grid index. Really I need to add a check to make it throw an error when this happens (you can tell because Policy is the wrong shape, but this is not tested for). Because it is the ‘lower grid index’ it will anyway be close to correct, but it is not good that it does so.

[robertdkirkby]

Just pushed an update of GL2017 example. Is now good as far as line 400. I want to overhaul the SimPanelValues for InfHorz TPath so will fix the second half later.

[robertdkirkby]

Alright, done. Just uploaded GL2017, runs right through nice and easy now. First transition takes a fair while, but that is just because the update step size is so small. Grid interpolation layer and divide-and-conquer means that it is now easy to run on my desktop.

[aledinola]

Hi Robert,
Could you please check also the case for infinite horizon with transition, no d code. Running my Huggett 1997 code with R2025b I do not get the reshape gpu error anymore, as expected, but I get a new error:

transpathoptions = 

  struct with fields:

                 verbose: 1
            weightscheme: 0
               tolerance: 1.0000e-03
              GEnewprice: 3
             GEnewprice3: [1×1 struct]
    updateaccuracycutoff: 1.0000e-09
                parallel: 2
           oldpathweight: 0
                  Ttheta: 1
                 maxiter: 1000
          graphpricepath: 0
        graphaggvarspath: 0
           graphGEcondns: 0
      historyofpricepath: 0
               stockvars: 0
          weightsforpath: [100×1 double]
          tanimprovement: 1
     useintermediateEqns: 0

Unrecognized function or variable 'n_d'.

Error in ValueFnIter_InfHorz_TPath_SingleStep_DC1_nod_raw (line 65)
Policy=UnKronPolicyIndexes_Case1(Policy,n_d,n_a,n_z,vfoptions);
                                        ^^^
Error in ValueFnIter_InfHorz_TPath_SingleStep (line 32)
                [VKron,PolicyKron]=ValueFnIter_InfHorz_TPath_SingleStep_DC1_nod_raw(VKron,n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_InfHorz_shooting_nod (line 138)
        [V, Policy]=ValueFnIter_InfHorz_TPath_SingleStep(V,0,n_a,n_z,[], a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_Case1 (line 360)
            [PricePath,GEcondnPath]=TransitionPath_InfHorz_shooting_nod(PricePath0, PricePathNames, PricePathSizeVec, ParamPath, ParamPathNames, ParamPathSizeVec, T, V_final, AgentDist_initial, n_a, n_z, pi_z, a_grid,z_gridvals, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, ReturnFnParamNames, GEeqnNames, vfoptions, simoptions,transpathoptions);
                                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in main_huggett_1997 (line 264)
PricePath=TransitionPath_Case1(PricePath0, ParamPath, T, V_final, StationaryDist_init, n_d, n_a, n_z, pi_z, d_grid,a_grid,z_grid, ReturnFn, FnsToEvaluate, TransPathGeneralEqmEqns, Params, DiscountFactorParamNames, transpathoptions,vfoptions, simoptions);
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The problem happens here:

TransitionPaths/InfHorz/subcodes/ValueFnSingleStep/DivideAndConquer/ValueFnIter_InfHorz_TPath_SingleStep_DC1_nod_raw.m

If you look at this function, there is no n_d variable. I think we already found the same error with @MichaelTiemann. The solution is simply to set n_d=0 in this function, as Michael had pointed out.

[aledinola]

Bug fixes for Infinite Horizon Transition No d Variable

1. Could you please update line 129 of Aiyagari1994TransitionPath.m from

StationaryDist_final=StationaryDist_Case1(Policy_final,n_d,n_a,n_z,pi_z);

to

StationaryDist_final=StationaryDist_Case1(Policy_final,n_d,n_a,n_z,pi_z,simoptions);

Note: Adding simoptions as last argument should be done not only here but for all instances of StationaryDist_Case1.

2. There is a mistake here plot(0:1:T, [p_eqm_init.r;PricePath.r]) since now PricePath.r is a row vector, replace ; with , (we are concatenating horizontally and not vertically).

3. Please add the following lines at the beginning of the script:

% The user can experiment with gridinterplayer=0 (pure discretization) or gridinterplayer=1 (linear interpolation b/w grid points). If gridinterplayer=1, then you must set vfoptions.divideandconquer=1 (required for transition).
vfoptions.gridinterplayer  = 1;
vfoptions.ngridinterp      = 20;
vfoptions.divideandconquer = 1;
simoptions.gridinterplayer = vfoptions.gridinterplayer;
simoptions.ngridinterp     = vfoptions.ngridinterp;

Item 3 is a suggestion, not a bug fix, obviously. However, I find it useful to test transition codes **with** and **without** grid interp layer to find out bugs.

P.S. I should learn at some point how to send pull requests :slight_smile: Whenever I try, I mess up my repo :slight_smile:

[MichaelTiemann]

Kia ora Alessandro,

I have felt the pain of doing pull requests wrong…and it is a pain. The easiest way I have found is this:

1. Start with a clean branch of upstream/main (or whatever the upstream branch is named). The upstream might be vfitoolkit or it might be your own repository, depending on the scope of changes you want to make. If you make changes that require multiple repositories to be updated, say so in your descriptions (later).
2. Create a new branch/branches for the fix/feature you want to implement it. It helps to give branches a descriptive name. When creating branches that need to sync across repos, a name that makes it clear to each repo owner what’s going on can be helpful. A branch named “minor bugfix” doesn’t give enough details to do cross-repo coordination. (It also doesn’t really support multiple people working on the source code at the same time, because it doesn’t convey any causal information.)
3. Your branches are private until you publish them. When published, github will make it easy for you to create a pull request. But if you haven’t made any changes, there’s nothing to pull (yet).
4. In your new branch, make and test the changes you want to make. If you need to test a change that depends on other repos being changed, all repos need to be at their self-consistent branch level. It is important to note that if one of the repos you depend upon is being modified by others, it can require you to do an update (either fast-forward or merge). This is the beauty and challenge of git and collaboration. When people are making small changes elsewhere in the code–no problem. When 50 files are being renamed and moved around in an upstream repo, it can require considerable care and patience. Best to know ahead of time about such changes, and avoid starting a fresh branch until things are settled. But as you gain experience (and as the team works out its own patterns of teamwork), git can make handling even complex development relatively painless.
5. When your changes are behaving as expected, commit your changes to your branch. If you haven’t yet published your branch, publish it and your commits will become visible to others.
6. When your branch is published, you can create a pull request (either using GitHub Desktop or visiting the repo on GitHub, where there should be a mention up top of new branches and their readiness for creating a pull request). The creation of the pull request will cause github to do some quick checks as to whether there are merge conflicts or not.
7. If there are merge conflicts, you must resolve them. This is the pain of collaborative development. If not, then you can happily wait for the maintainer (which may be you in your repo) to accept the pull request. You might not want to merge pull your own changes until other repos you depend upon do their merges. Version numbers can help you communicate what levels of which libraries are tested with your software.
8. When a pull request has been merged, you can switch to the upstream repo (which has just been updated), delete your feature branch, and be ready to start some new development, back at step 1.

Like anything else, it takes practice. And like anything, failure is the best teacher. Doing one single feature at a time is easiest. The more complicated your feature work (especially stacking different features on different branches all at the same time) can become a nightmare. Wishing you good luck,

Michael

[aledinola]

GL2017

Hi Rob, I found an issue with the policy function for consumption plotted in Figure 1. It has a strange wiggle and, more importantly, it is negative for low values of wealth.

I show in the table below a few values of the policy function of consumption, size [n_a,n_z] = [401,13]

You can see that consumption<0 occurs for \theta_2 and \theta_3.

The same problem arises with pure discretization. Below I show the policy function for consumption with pure discretization.

I’m really puzzled by this: negative consumption cannot happen in equilibrium, unless the budget constraint is empty, but in that case the problem is not specified correctly. Also, notice that in the states where optimal consumption is negative, optimal labor supply is zero, which is also strange: for very low wealth, agents would optimally want to work a lot in these models.

[javier_fernan]

Hi @aledinola, your finding is a bit concerning. I don’t understand how this can happen: if consumption is negative, F is set to -inf, reading the function GuerrieriLorenzoni2017_ReturnFn. There is another -inf later on, maybe this is the problem?

@robertdkirkby Can you please shed some light on this?

[robertdkirkby]

Is all just places with V=-Inf and zero agent mass.

Here is a version of the graph.

First panel is much the same thing as in Alessandro’s post.
Second panel is same, but now only showing parts of the state space with V is finite (so get rid of the V=-Inf)
Third panel is same, but now only showing parts of the state space where the mass of agents in the stationary distribution is positive (so not zero).

As can be seen, all the ‘oddities’ are just from the parts with V=-Inf and mass=0.

Following is the code generating the figure, put it straight after solving the initial general eqm.

FnsToEvaluate.Consumption = @(d, aprime,a,z,r, v, B, Bprime) GuerrieriLorenzoni2017_ConsumptionFn(d, aprime, a, z,r, v, B, Bprime); % Consumption
ValuesOnGrid_initial=EvalFnOnAgentDist_ValuesOnGrid_Case1(Policy_initial,FnsToEvaluate,Params_initial,[],n_d,n_a,n_z,d_grid,a_grid,z_grid,simoptions);
figure(1)
subplot(3,1,1); plot(a_grid,ValuesOnGrid_initial.Consumption(:,1),a_grid,ValuesOnGrid_initial.Consumption(:,6),a_grid,ValuesOnGrid_initial.Consumption(:,13))
legend('z1','z6','z13')
title('Consumption')
Vfinite=isfinite(V_initial);
subplot(3,1,2); plot(a_grid,Vfinite(:,1).*ValuesOnGrid_initial.Consumption(:,1),a_grid,Vfinite(:,6).*ValuesOnGrid_initial.Consumption(:,6),a_grid,Vfinite(:,13).*ValuesOnGrid_initial.Consumption(:,13))
legend('z1','z6','z13')
title('Consumption where isfinite(V)')
positivemass=(StationaryDist_initial>0);
subplot(3,1,3); plot(a_grid,positivemass(:,1).*ValuesOnGrid_initial.Consumption(:,1),a_grid,positivemass(:,6).*ValuesOnGrid_initial.Consumption(:,6),a_grid,positivemass(:,13).*ValuesOnGrid_initial.Consumption(:,13))
legend('z1','z6','z13')
title('Consumption where mass>0')

PS. There was also one typo, namely I used the Params from final eqm with the Policy from initial eqm to create ValuesOnGrid that was used draw the graph. Have updated GL2017 example to fix this.