Replication of Huggett (JME 1997)

aledinola · December 29, 2025, 12:41pm

I wrote a replication of Huggett 1997:

Huggett, Mark, 1997, “The one-sector growth model with idiosyncratic shocks: Steady states and dynamics,” “Journal of Monetary Economics”, Elsevier, volume 39, issue 3, pages 385-403, August.

My replication code is available in this repo.

Problem 1: Bug
Unfortunately, the transition does not run due to a bug in ValueFnIter_InfHorz_TPath_SingleStep_nod_raw. The error message is the following:

Unrecognized function or variable 'n_d'.

Error in ValueFnIter_InfHorz_TPath_SingleStep_nod_raw (line 58)
Policy=UnKronPolicyIndexes_Case1(Policy,n_d,n_a,n_z,vfoptions);
                                        ^^^
Error in ValueFnIter_InfHorz_TPath_SingleStep (line 24)
            [VKron,PolicyKron]=ValueFnIter_InfHorz_TPath_SingleStep_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 144)
        [V, Policy]=ValueFnIter_InfHorz_TPath_SingleStep(Vnext,0,n_a,n_z,[], a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_Case1 (line 349)
            [PricePath,GEcondnPath]=TransitionPath_InfHorz_shooting_nod(PricePathOld, 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 185)
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 bug is quite easy to spot: In line 58 of ValueFnIter_InfHorz_TPath_SingleStep_nod_raw there is the line

%% Policy in transition paths
Policy=UnKronPolicyIndexes_Case1(Policy,n_d,n_a,n_z,vfoptions);

but n_d is not available in the function.

Problem 2: Transition not coded with interpolation?
The above error message happens if vfoptions.gridinterplayer = 0. If I set instead vfoptions.gridinterplayer = 1, there is a toolkit message that says that this case is not implemented yet. Would it be possible to flag this as a priority? I guess that computing a transition in infinite horizon models, with interpolation, is a relatively common case; it could be used for example as a robustness check in Bruggeman 2021. Many thanks!

MichaelTiemann · December 29, 2025, 8:01pm

It looks like UnKronPolicyIndexes_Case1 tolerates passing in 0 (zero) for n_d at line 58. Have you tried that? n_d being zero is effectively the statement that there is no d_grid.

aledinola · December 29, 2025, 8:28pm

Thanks, I will try that (replace n_d with 0 at line 58)

aledinola · December 30, 2025, 12:03pm

I passed n_d=0 in line 58 of ValueFnIter_InfHorz_TPath_SingleStep_nod_raw. Unfortunately there is now a new error:

Unable to perform assignment because the size of the left side is 1-by-1000-by-1 and the size of the right
side is 1-by-1000-by-2.

Error in TransitionPath_InfHorz_shooting_nod (line 149)
            PolicyIndexesPath(:,:,T-ttr)=Policy;
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_Case1 (line 349)
            [PricePath,GEcondnPath]=TransitionPath_InfHorz_shooting_nod(PricePathOld, 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);
                                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Policy is a 3-dim array, so I modified line 149 of TransitionPath_InfHorz_shooting_nod as following:

PolicyIndexesPath(:,:,:,T-ttr)=Policy;

and it is ok, but then I get another error later on, since l_d is not defined.

yechen · December 30, 2025, 1:34pm

I just commented out the relevant lines of code and tested the Aiyagairi1994TransitionPath case—it seems to work.

aledinola · December 30, 2025, 2:06pm

Thanks for checking, @yechen.
This is strange, if I ran the script Aiyagari1994TransitionPath.m, I get this error:

Unrecognized function or variable 'n_d'.

Error in ValueFnIter_InfHorz_TPath_SingleStep_nod_raw (line 58)
Policy=UnKronPolicyIndexes_Case1(Policy,n_d,n_a,n_z,vfoptions);
                                        ^^^
Error in ValueFnIter_InfHorz_TPath_SingleStep (line 24)
            [VKron,PolicyKron]=ValueFnIter_InfHorz_TPath_SingleStep_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 144)
        [V, Policy]=ValueFnIter_InfHorz_TPath_SingleStep(Vnext,0,n_a,n_z,[], a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_Case1 (line 349)
            [PricePath,GEcondnPath]=TransitionPath_InfHorz_shooting_nod(PricePathOld, 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 Aiyagari1994TransitionPath (line 174)
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);
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

I am using the latest version of the toolkit files from github. I suspect that maybe some recent update broke the codes for the transition.

yechen · December 30, 2025, 3:08pm

I copied all of the codes in the ValueFnIter_InfHorz_TPath_SingleStep_nod_raw and commented out the 58th line about the Policy function. Then it can work.

MichaelTiemann · December 30, 2025, 6:13pm

I’m sorry…in which repo is Huggett (JME 1997)? I see Huggett 1996 in the replication repo, but perhaps you are referring to a different repo?

aledinola · December 30, 2025, 9:47pm

Huggett JME 1997 is a replication I am working on but I haven’t updated yet on github; I will soon.

However, the same error happens also if you run the Aiyagari transition example, available here:

github.com/vfitoolkit/VFItoolkit-matlab-examples

HeterogeneousAgentModels/Aiyagari1994TransitionPath.m

master

% Example of computing a general eqm transition path for the model of Aiyagari (1994).
%
% These codes set up and solve the Aiyagari (1994) model for a given
% parametrization. They then show how to solve for the general equilibrium
% transition path in reposonse to a 'surprise' one off change in the
% parameter beta (the time discount parameter).
%
% VFI Toolkit automatically detects hardware (GPU? Number of CPUs?) and
% sets defaults accordingly. It will run without a GPU, but slowly. It is
% indended for use with GPU.

%% Set some basic variables

n_k=2^9;%2^9;
n_z=21; %21;

%Parameters
Params.beta=0.96;
Params.alpha=0.36; % This is actually redundant as declare this below when looking at initial and final eqm
Params.delta=0.08;

This file has been truncated. show original

If you have some time and you have a chance to run the script Aiyagari1994TransitionPath, please let me know how it goes. Thanks!

aledinola · December 30, 2025, 10:06pm

I have uploaded my replication attempt on github: GitHub - aledinola/Huggett_1997

You can run main_huggett_1997.m On my computer, I get an error in line 187.

MichaelTiemann · December 30, 2025, 10:44pm

I’ll take a look after lunchtime (in Wellington).

MichaelTiemann · December 31, 2025, 12:41am

So…there’s some internal VFIToolkit stuff regarding the shapes of PolicyIndexesPath in the no_d / vfoptions.gridinterplayer==0 case which I’m still trying to figure out. But I created a pull request in the VFIToolkit-matlab-examples repo to fix a calling convention error that will get you no matter what. With that done (and the program actually completing) I’m going to see how to make less of a mess dealing with the shapes of PolicyIndexesPath.

MichaelTiemann · December 31, 2025, 12:57am

I cannot report back on Hugget_1997 until I sort the toolkit changes I made for Aiyagari1994TransitionPath. While that seemed to converge, Huggett_1997 almost converged, but it’s been oscillating from quite some time. I’ll report back when I have more info.

yechen · December 31, 2025, 6:16am

@aledinola I first ran Huggett (1997) with a transition horizon of T = 150, but the algorithm hit the iteration limit without converging. When I shortened the horizon to T = 100, the code converged, yet the aggregate capital stock drifted downward in the final decades of the transition, suggesting the economy never fully reaches the new steady state. Moreover, the resulting capital path is the mirror image of Fig. 1 in the paper, which is a wrong pattern.

MichaelTiemann · December 31, 2025, 7:49am

Well, I don’t see an obvious reason why it doesn’t converge. Check out the changes bundled in this commit 294663217e8310464ead222426d410635fd1db38 (in this PR: Resync fixes 2025 alt by MichaelTiemann · Pull Request #29 · vfitoolkit/VFIToolkit-matlab · GitHub ).

aledinola · December 31, 2025, 1:50pm

Thanks for the pull request!

aledinola · December 31, 2025, 2:41pm

Very interesting! I have not managed yet to run the transition using the vfi toolkit, but I have coded myself the Huggett 1997 model. The files are available in my repo here: GitHub - aledinola/Huggett_1997, see folder cpu.

I set T = 100 and the transition converges, getting similar results to what you have shown. I have a number of observations that I summarize below:

I set a high number of grid points for assets, na=2000, to make sure the computation is accurate enough. The model is simple so we can afford many grid points.
I compare the stationary equilibrium in my code and in the toolkit and there are small differences, using exactly the same parameters for tolerance, max number of iterations etc.
Toolkit results:

RESULTS FINAL STEADY STATE, VFI TOOLKIT
Algorithm: pure discretization
No. grid points assets: 2000 
CapitalMarket residual: -0.004739 
Goods market residual:  -0.000195 
Aggregate capital:      4.311640 
Capital-to-labor ratio:  4.311640 
Capital-to-output ratio: 2.547835 
Consumption:             1.261308 
Interest rate:           0.041296 
Wage:                    1.083057 
Run time General Equil: 8.772723

My code results:

FINAL STEADY-STATE, OWN CPU CODE: 
Algorithm: pure discretization
No. grid points assets: 2000 
CapitalMarket residual: 0.001356 
Goods market residual:  0.000052 
Capital stock:           4.311645 
Capital-to-labor ratio:  4.311645 
Capital-to-output ratio: 2.547836 
Consumption:             1.261061 
Interest rate:           0.041296 
Wage:                    1.083057 
Run time General Equil: 6.970329

As I said, I am not able to compute the transition with the toolkit, but I computed it with my code and it converges with tolerance = 10^(-3) and T=100 and dampening parameter equal to 0.9
The transition seems pretty robust but completely different from what Huggett (1997) reports.
As @yechen observed, I also noticed that towards the end, K(t) seems to drift down. Does this mean that the transition is not accurate? I don’t know. I had found a similar problem in this older post.

In the first plot above I show K(t)/L(t) but this is the same as K(t) since average labor is always equal to 1 in this particular model.

Update
Huggett (1997) writes that he used T=1000 while I used T=100. I will recompute the transition with T=1000, but I am sceptical that the results will change substantially. I assume the original code is not available since this is an old paper…

MichaelTiemann · December 31, 2025, 4:47pm

You have a fast CPU! My toolkit results were a fraction faster (8.67s vs 8.77s), but your CPU was almost 2x faster than time (6.97s vs 13.75s). I have replicated all your other numbers.

Happy New Year!

aledinola · December 31, 2025, 5:09pm

On my laptop I have a good cpu but a bad gpu (only 4 gb of ram)

Happy New Year!!

robertdkirkby · December 31, 2025, 10:06pm

Just ran your Huggett (1997) code (toolkit version). It runs and solves just fine. Obviously this is not the same as saying anything about whether or not it is correct [toolkit will be as it is same internal commands that solves the Aiyagari 1994 transition paths].

Those little wiggles at the end are partly because the toolkit never updates period T (the final period values of the PricePath are kept unchanged, precisely to emphasize any jump at the end of the tail). You tend to be able to eliminate them if you want to by just taking the answer, setting a new initial guess based on smoothing out that tail of the answer you got, and then solving again. If you look they are all just differences in K of about 1e-3 or 1e-4, so is just numerical error that can be cleaned up if you want to.

Two minor comments on code. You currently have the stationary eqm command output as [p_eqm_final,~,GeneralEqmCondn_final]=, but I changed it a month or three back so you can now just do the cleaner [p_eqm_final,GeneralEqmCondn_final]=.

I also changed the transition path command so that instead of current output of just PricePath= you can now ask for [PricePath,GeneralEqmCondnPath]=. Which is handy if you are solving lots of transition paths (say in a loop) as you can see later if they converged to a good solution or not.