Yes, I noticed the paper reports the periods for the transition are 1000. I did test the code at T = 150, 200, 500, and 1000 yesterday. Convergence became markedly slow for T > 150, and none of the runs reached the default tolerance before hitting the 1,000-iteration limit.
1 Like
Just to comment that convergence of the transition path is calculated with the distance being measured as
% See how far apart the price paths are
PricePathDist=max(abs(reshape(PricePathNew(1:T-1,:)-PricePathOld(1:T-1,:),[numel(PricePathOld(1:T-1,:)),1])));
that is, distance is in the L_{\infty} norm. This PricePathDist is what gets compared to transpathoptions.tolerance (default=10^(-5))to determine when convergence is reached.
That the PricePath does not converge before hitting the transpathoptions.maxiter (default=1000) iteration limit before reaching convergence can also be related to how the damping parameters are set in the update rules (obviously more damping means more iterations will be required for any given convergence tolerance).
Run-time for doing a single iteration should be essentially linear in T (but changing T potentially means you need more/less iterations).
2 Likes
I checked just now and I get this error:
Error using gpuArray/reshape
Sparse gpuArray matrices are not supported for this function.
Error in StationaryDist_InfHorz_TPath_SingleStep (line 9)
AgentDist=reshape(Gammatranspose*AgentDist,[N_a,N_z]); %No point checking distance every single iteration. Do 100, then check.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_InfHorz_shooting_nod (line 217)
AgentDistnext=StationaryDist_InfHorz_TPath_SingleStep(AgentDist,PolicyaprimezPath(:,tt),II1,IIones,N_a,N_z,pi_z_sparse);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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 238)
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);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Is it related to using R2024b where reshape of sparse gpu arrays does not work? I thought that the toolkit would not use this feature. If it does, it is going a problem for me because the HPC server at my university has updated Matlab only to R2024b (after my insistence) but said they would not update to R2025b for a while, due to some technical difficulties I haven’t completely understood (apparently Matlab switched from Java to something else and this is creating problems for the server guys) 
I realised while working on a different model that the runtime gain from moving Tan improvement (agent distribution iteration) to the sparse GPU matrices was minor (decent speedup, >2x even approaching 5 or 10x sometimes, but agent dist iteration is often only taking 1/100th of what the value fn iteration takes, so it is a good speed improvement on something that is already almost none of the runtime). But it also creates a new memory bottleneck, as the gpu just has way less memory than cpu. As a result I have decided not to move any of the Tan iterations to the GPU, I will leave them all on CPU, as the runtime gain is minor and the bottleneck is real.
Clearly though this one (InfHorz TransPath) is currently on GPU. I will move it back to CPU. Hopefully I can find time early next week to do this.
[The model in question was an OLG with two endogenous states, one standard and one experienceasset for human capital. We run it on a server with big grids and with Tan improvement on the GPU we get out of memory errors there even for grids that can be solved for the value fn.]
In an ideal world I would move the Tan improvement for one endogenous state models to GPU, and keep it on CPU for anything with two or more endo states; you are unlikely to out-of-memory in one endo state models. But as it is currently done these both use the same command, so easier just to do everything on CPU.
Just pushed. I think that the agent dist for InfHorz TPath is now all doing Tan improvement on CPU. @aledinola I don’t really have an easy way to test it though, so let me know if it doesn’t run for you.
1 Like
Hi Rob, the error is still there (see update below):
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
Error using gpuArray/reshape
Sparse gpuArray matrices are not supported for this function.
Error in AgentDist_InfHorz_TPath_SingleStep (line 12)
AgentDist=reshape(AgentDist*pi_z_sparse,[N_a*N_z,1]);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_InfHorz_shooting_nod (line 211)
AgentDistnext=AgentDist_InfHorz_TPath_SingleStep(AgentDist,PolicyaprimezPath(:,tt),II1,IIones,N_a,N_z,pi_z_sparse);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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 238)
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);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Update
The same error happens even if I try to do the transition with gridinterplayer=1 and divide and conquer =1. But now notice that the error happens in a different function: AgentDist_InfHorz_TPath_SingleStep_nProbs
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
Error using gpuArray/reshape
Sparse gpuArray matrices are not supported for this function.
Error in AgentDist_InfHorz_TPath_SingleStep_nProbs (line 13)
AgentDist=reshape(AgentDist*pi_z_sparse,[N_a*N_z,1]);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in TransitionPath_InfHorz_shooting_nod (line 213)
AgentDistnext=AgentDist_InfHorz_TPath_SingleStep_nProbs(AgentDist,PolicyaprimezPath(:,:,tt),II2,PolicyProbsPath(:,:,tt),N_a,N_z,pi_z_sparse);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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 239)
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);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
I’ll have another go on Monday and use whos to check that things are on CPU when they should be.
1 Like
@aledinola, hopefully fixed now (agent dist on CPU). Let me know if it still errors.
1 Like
Ran the transition successfully. You can remove or comment out line 211 in TransitionPath_InfHorz_shooting_nod:
whos AgentDist PolicyaprimezPath II1 IIones pi_z_sparse
1 Like
I’m trying to understand why the transition dynamics display a reversal of the pattern described in the paper, effectively contradicting the author’s conclusions. Initially, 20 % of households are credit-constrained and hold zero capital, exactly as the author assumes. Along the transition to the stationary distribution, however, the share of borrowing-constrained households falls. Consequently, should aggregate capital first rise and then decline toward its steady-state level? Any thoughts on this?
1 Like
My thought is that it will depend on how you set the other 80%. If you set to other 80% to have same assets as before, then initial capital is lower once you put the 20% with zero assets. But if you set the other 80% to have enough assets that total assets is unchanged (so scale up the 80% to compensate for the 20% with zero) you are going to get something rather different.
Does Huggett (1997) talk about what he does with the other 80%? Because this will tell you if aggregate capital has to rise over the path (which is what the graph Alessandro put in an earlier post shows).
2 Likes
Huggett (1997): The initial distribution puts 20% of the agents exactly at zero asset holdings and equal numbers of agents at all capital levels between 0 and 10.8104. Thus, the aggregate capital stock is 4.3242 which is the steadystate capital stock that I calculate for this economy.
I checked my replication and I set the initial distribution in this way:
a_cut_init = find(a_grid>10.81, 1 );
pdf_assets = zeros(n_a,1);
pdf_assets(1) = 0.2;
pdf_assets(2:a_cut_init) = 0.8/numel(2:a_cut_init);
StationaryDist_init = zeros(n_a,n_z);
StationaryDist_init(:,1) = 0.5*pdf_assets;
StationaryDist_init(:,2) = 0.5*pdf_assets;
StationaryDist_init = StationaryDist_init/sum(StationaryDist_init,"all");
Total assets using this initial distribution are
sum(a_grid.*StationaryDist_init,"all") = 2.17
which are lower than the value in the final equilibrium (this confirms @robertdkirkby intuition) In my code I set up the initial distribution exactly as Huggett describes, but I don’t know how he defines the grid for assets (is it equally spaced? What is the upper bound? If it is not equally spaced, is it double exp?). I think he chose the cutoff 10.8104 so that the initial equilibrium replicates the aggregate assets in the stationary equilibrium. I do replicate aggregate assets in the stationary equilibrium: my values is 4.3116, his value is 4.3242.
@yechen If you are interested in this paper, you can recompute the transition using 21.5 as cutoff. With this value, I get
sum(a_grid.*StationaryDist_init,“all”)
Aggregate capital in initial distribution:
ans =
4.31078954174934
1 Like
To be more precise, here is the updated code to calculate the correct cutoff:
%% Define initial distribution via cutoff on asset grid
% Huggett-style initialization:
% - 20% of agents are exactly at zero assets (a = 0)
% - Remaining 80% are spread uniformly over assets up to a cutoff
%
% We choose the cutoff by bisection so that aggregate capital in the
% initial distribution is close to the stationary value p_eqm_final.K.
c_low = a_grid(1);
c_high = a_grid(end);
maxiter = 50;
tol = 1e-6; % tolerance on aggregate capital
for it = 1:maxiter
% Midpoint cutoff
c_mid = 0.5 * (c_low + c_high);
% First grid point strictly above the cutoff
a_cut_init = find(a_grid > c_mid, 1);
% Cross-sectional distribution over assets
pdf_assets = zeros(n_a,1);
pdf_assets(1) = 0.20; % 20% at zero assets
if ~isempty(a_cut_init)
% Spread remaining 80% uniformly across 2:a_cut_init
pdf_assets(2:a_cut_init) = 0.80 / numel(2:a_cut_init);
end
% Joint distribution over (a,z)
StationaryDist_init = zeros(n_a,n_z);
StationaryDist_init(:,1) = 0.5 * pdf_assets;
StationaryDist_init(:,2) = 0.5 * pdf_assets;
% Normalize (should already sum to 1 but safe)
StationaryDist_init = StationaryDist_init ./ sum(StationaryDist_init,"all");
% Aggregate capital implied by current cutoff
K_mid = sum(a_grid .* StationaryDist_init,"all");
fprintf('Iter %2d: cutoff c = %.6f, K = %.6f (target %.6f)\n', ...
it, c_mid, K_mid, p_eqm_final.K);
% Check stopping rule
if abs(K_mid - p_eqm_final.K) < tol
fprintf('Tolerance reached. Stopping bisection.\n');
break
end
% Bisection update
if K_mid > p_eqm_final.K
% Too much capital -> reduce cutoff
c_high = c_mid;
else
% Too little capital -> increase cutoff
c_low = c_mid;
end
end %end iter bisection
Sadly, the transition does not converge 
@robertdkirkby: What is the default tolerance for transition paths? Can I change it by setting
transpathoptions.tolerance = 1e-4;
?
Also, I don’t understand the output on the screen
Number of iterations on transition path: 12
Current distance between old and new price path (in L-Infinity norm): 0.561249
Current distance to convergence: 5612.49 (convergence when reaches 1)
What does Current distance to convergence mean? I think the algorithm stops when the distance between two successive iterations falls below the tolerance, no?
Note that 0.2*0+0.8*(10.8104/2)=4.3242. So if you put evenly spaced grid points from 0 to 10.8104 and evenly weight a mass of 0.8 agents onto them, then you will get the same as Huggett. [Use an odd number of points, including the end points, so that there is a point exactly in the middle, otherwise the aggregate will just be approximate.] It is a safe bet this is what Huggett did given how cleanly the numbers play out.
Currently you use:
a_grid = a_min + (a_max - a_min) * (linspace(0, 1, n_a)'.^3);
with a_min=0 so as long as a_max>10.8104 you can just switch this to
a_grid=[linspace(0,10.8104,floor(n_a/2)),linspace(10.8104,a_max),n_a-floor(n_a/2)+1]';
a_grid=unique(a_grid); % eliminate the double value at 10.8104
and then it will be easy enough to set up the agent dist for this exercise. [now n_a/2 should be odd, so n_a should be even; otherwise you will get close to but not exact 4.3242]
1 Like
Thanks @robertdkirkby, this is a great suggestion. There is a small typo in the construction of a_grid:
a_grid = [linspace(0,10.8104,floor(n_a/2)),linspace(10.8104,a_max,n_a-floor(n_a/2)+1)]';
Even with this change, though, the transition does not converge: keeps oscillating.
Default for InfHorz transition paths is
transpathoptions.tolerance=1e-5;
Number of iterations on transition path: 12
Current distance between old and new price path (in L-Infinity norm): 0.561249
Current distance to convergence: 5612.49 (convergence when reaches 1)
The second line of this is the distance, and it will converge when this reaches simoptions.tolerance.
The third line is the ratio of the distance to the tolerance, and will converge when this reaches 1. I just put this there as it is much easier to read and see how close, or not, you are to convergence (for one, you don’t need to remember what the tolerance is set at). But as you highlight it is not well named. I will think a bit and give it a new name, along the lines of “Current ratio of distance to convergence tolerance: XX (convergence when reaches 1)”
1 Like
In my experience the most common reason for oscillating close to convergence is that you have simoptions.tolerance to tight for the number of grid points. So either just loosen tolerance or increase grid points. [Obviously this is not the only reason you can see oscillations, but it is a common one.]
1 Like
Not sure… I have
n_a = 2000; % number of asset grid points
n_z = 2; % number of idiosyncratic productivity states
T = 500;
transpathoptions.GEnewprice = 3;
transpathoptions.tolerance = 1e-4;
% Specify how to update GE prices from GE conditions (for price scheme 3).
% A row has the form: {GEcondnName, priceName, add, factor}
% Here: CapitalMarket > 0 means K is too big, so we want to reduce K.
transpathoptions.GEnewprice3.howtoupdate = …
{‘CapitalMarket’, ‘K’, 0, 0.1};
1 Like
@aledinola Do you really need T=500? I ran your code and it was quite slow, then I realized that it is slow since there are so many time periods in the transition.