Expected Value Fn

robertdkirkby · September 12, 2026, 10:23pm

Info for those interested in toolkit internals.

Previously, the toolkit always computed the expected next period value fn as EV=V.*pi_z.

I have changed this to use EV=V*pi_z.

The impact is that small models will be 1% slower for the value fn iteration, but large models are substantially faster and much larger models can be run due to avoiding out-of-memory errors.

Currently the change only impacts the markov z in models with standard endogenous states.

You can see a bunch of runtime tests that underlie the decision at:

Most users will notice no appreciable difference, except if you try solving a model with prod(n_z)>1000 in which case it will be much better.

PS. Obviously I am glossing over a bunch of details in how exactly the two are implemented, especially with relation to the way they handle -Inf in the value fn. But these details all just follow from trying to implement the two different .* vs *

robertdkirkby · September 18, 2026, 11:55pm

Another change that I am rolling out throughout the toolkit. It is about ‘post-processing’. Internally, take what user input, reshape everything into form used by toolkit, solve, then reshape everything into form output to user. This is about improving the “then reshape everything into form output to user”. It gives identical output, but just uses a little less memory on the way. [This will be rolling out across the toolkit over the next week or two. Is already in some parts.]

Following is a short explanation from Claude of what is being done.

The grid-interpolation-layer (GI) solvers store the policy as a small stack of rows — a midpoint index on the coarse asset grid, a second-layer index for the fine point between coarse nodes, and a flag recording whether a −Inf return at one of the two bracketing coarse points forced all the weight onto the other. After the age loop, the midpoint/second-layer pair is converted to “lower grid point + offset” and the flag is attached. The first change (Hotspot 1) rewrites the conversion line. It used to be a masked blend, Policy(2,:,:,:)=adjust.Policy(2,:,:,:)+(1-adjust).(Policy(2,:,:,:)-n2short-1);, and is now the algebraically identical Policy(2,:,:,:)=Policy(2,:,:,:)-(n2short+1)*(~adjust); — where adjust is true, nothing changes; where it is false, n2short+1 is subtracted. This is exact on the integer indices (verified over the whole reachable domain), gives a small speedup, and is memory-neutral in practice: MATLAB already fuses the old elementwise expression, so the “extra temporaries” it looked like it created were not the binding allocation.

The second change (Hotspot 2) is the one that actually saves memory. The flag used to live in a separate array, PolicyL2flag=2*ones(1,N_a,N_z,N_j,‘gpuArray’);, written per age as PolicyL2flag(1,:,:,jj)=2+(inLowerStrict&isInfLower)-(inUpperStrict&isInfUpper); and glued on at the end with Policy=[Policy; PolicyL2flag];. Now Policy is preallocated one row taller, Policy=zeros(3,N_a,N_z,N_j,‘gpuArray’); Policy(3,:,:,:)=2;, the per-age write goes straight into that row, Policy(3,:,:,jj)=2+(inLowerStrict&isInfLower)-(inUpperStrict&isInfUpper);, and the concatenation is deleted. The steady-state footprint is unchanged (K+1 rows either way — the flag was always a full row), but the concatenation had to allocate a fresh K+1-row result while its inputs, and the age loop’s large temporaries that are still in scope, were all live; that allocation was the whole solve’s peak. On a full solve at the Guvenen (2007) HIP sizes (N_z=3965, N_j=71, n_a=301–401) removing it lowers the peak GPU memory by about one Policy row (0.55–0.74 GiB) at every size, with output byte-identical to the original in every test bank run so far. The final Policy has the same shape and values as before, so nothing downstream — PolicyInd2Val, simulation, ValueFnFromPolicy — is affected.