The Sequential Endogenous Grid Method

1Introduction¶

Solving heterogeneous-agent models with multiple simultaneous decisions has traditionally required either joint optimization, which is slow, or restrictive assumptions about preferences that limit economic relevance. These heterogeneous-agent models, in the tradition of Huggett (1993), Aiyagari (1994), and Krusell & Smith (1998), are the workhorses of macroeconomic analysis Krueger et al., 2016, and increasingly they ask the household to make several decisions at once, such as labor supply, portfolio choice, or illiquid saving, the setting in which joint optimization grows costly. The Endogenous Grid Method (EGM) transformed computation of consumption-saving problems by inverting first-order conditions, but extending that insight to multiple continuous choices has proved difficult. The leading extensions either constrain the problem class or keep a numerical solve inside the iteration Iskhakov, 2015Druedahl & Jørgensen, 2017Druedahl, 2021.

A household choosing consumption and labor supply makes both decisions with the same information. Yet the solver need not treat them jointly. I introduce Sequential EGM (EGM${}^n$, whose superscript records the chain of stagewise inversions), which orders the problem into stages and applies an EGM inversion wherever the stage permits one. Stages without an inversion fall back to lower-dimensional optimization.

Each stage of an EGM${}^n$ decomposition applies EGM when its first-order condition and its recovery map, the composite of the inverted condition with the budget identities, satisfy the inversion and optimality conditions developed in the paper. The recovery map carries each post-decision gridpoint to its endogenous pre-decision state, as $𝔪(a)=a+𝔠(a)$ does in one-dimensional EGM, where $a$ is end-of-period saving, $𝔠$ the consumption the inverted Euler equation assigns to it, and $𝔪$ the market resources from which that choice is optimal. Throughout, choices are continuous and satisfy standard constraint qualifications. EGM${}^n$ requires stronger stage structure than nested approaches that embed optimization within the EGM iteration, but it avoids optimization at qualifying stages and parallelizes across grid cells because the closed-form inversions have no cross-cell dependencies.

On the two-asset model of Auclert et al. (2021), with both solvers in my own matched implementation, EGM${}^n$ is faster than the Nested Endogenous Grid Method (NEGM) at every grid and thread count tested, by roughly an order of magnitude in serial and by about 3.7 to 7.9 times even at sixteen matched threads. Threads narrow the gap because the nested search has per-node work to spread across them, while the closed-form inversion has almost none. For the curvilinear grids the inversions produce, I develop ENGINE (ENdogenous Grid INterpolation and Extrapolation), with which the full solve runs roughly three to four and a half times faster than with the curvilinear interpolation previously proposed by White (2015) for such grids, at matched accuracy on the two-asset application.

As a direct product of the backward step, the decomposition also returns intermediate-stage marginal values, quantities single-block solvers never construct. In the two-asset model, their ratio minus one is the adjustment-cost wedge that prices converting illiquid wealth into liquid.

The Endogenous Grid Method was built to avoid the cost of searching over choices: for a three-choice problem with $G$ grid points per dimension, that search costs $O(G^3)$ evaluations at every state, a price paid across the whole state space.^[1] Carroll (2006) introduces it for the consumption-saving problem: starting from a grid of post-decision states, EGM inverts the first-order condition to recover the optimal policy and the pre-decision state that implies it, so a numerical search becomes a function evaluation. Barillas & Fernández-Villaverde (2007) extends it to a second control, the labor-leisure choice, in a neoclassical growth model.^[2]

In the decade after EGM’s introduction, its scope broadened along two fronts. On the structural side, Fella (2014) relaxes continuity requirements for non-convex problems and Hintermaier & Koeniger (2010) handles occasionally binding constraints among endogenous variables; Druedahl & Jørgensen (2017) combines both advances into G2EGM, which admits multiple controls and states. On the discrete-choice side, Iskhakov et al. (2017) incorporates jumps between alternatives through extreme value errors. Their method, however, remains restricted to a single control and state variable.^[3] Separately, Clausen & Strub (2020) develops Reverse Calculus for nested and non-convex optimization; G2EGM Druedahl & Jørgensen, 2017 builds its first-order-condition arguments on an earlier working-paper version of these differentiability results.

EGM${}^n$ builds most directly on the multivariate extensions. White (2015) formalizes the conditions under which EGM applies to multidimensional problems and builds interpolation methods for the curvilinear grids that result; Ludwig & Schön (2018) proposes Delaunay triangulation for the same grids, though triangulation costs can offset the time saved by avoiding grid search. Iskhakov (2015) characterizes a triangular structure in the system of first-order conditions that permits sequential solution by substitution, with no root-finding.^[4] Triangular EGM asks for triangularity of the entire system; EGM${}^n$ asks only for sequential invertibility, stage by stage, which admits the mixed structures described above. Druedahl (2021) develops NEGM, which pairs an upper envelope with a one-dimensional numerical solve nested in the outer EGM iteration; nesting cuts complexity relative to joint optimization but keeps the embedded solve, which EGM${}^n$ removes.^[5]

Section 2 develops the decomposition on a two-choice consumption-labor problem, the smallest example that shows the chaining. Section 3 takes it to the canonical two-asset liquid-illiquid household block of Auclert et al. (2021), the household side of a typical HANK model, in which two asset states generate two-dimensional curvilinear grids and the benchmarks above are measured. ENGINE exists to interpolate those grids, and Section 4 develops and tests it. Section 5 gives practical guidelines for splitting a problem into stages and two sufficient structures, separable utility and invertible transition marginals, under which a stage inverts, and Section 6 marks where the method stops applying.

2The Sequential Endogenous Grid Method¶

A household that simultaneously chooses consumption and labor supply presents the joint solver with a two-dimensional search at every state space point. But when decisions are separable in stages, though not necessarily globally, sequential decomposition lets EGM inversion solve each stage efficiently.

2.1Problem Setup¶

I demonstrate EGM${}^n$ on a consumption-labor model in which the household chooses leisure and consumption under a risky wage. It is the smallest example that shows the chaining: two simultaneous decisions, a sequential ordering in which each stage sheds a state variable, and an EGM inversion at both stages, only one of which is the consumption margin EGM was built for. The household maximizes the present discounted lifetime utility of consumption and leisure over a finite horizon of $T+1$ periods, where $t=0$ denotes the beginning of the planning horizon (the first period of economic life) and $t=T$ the terminal period:

where ${\mathrm{V}}_0$ is the lifetime value function at the start of period zero, $a_0$ is the stock of liquid assets the household starts with, ${\theta }_0$ is the transitory wage shock, and ${𝔼}_0$ denotes the expectation conditional on information available at the start of the planning horizon.

Period utility is additively separable in consumption and leisure:

where $\rho >0$ is the coefficient of relative risk aversion, $\zeta >0$ determines the curvature of leisure utility (and thus the Frisch elasticity of labor supply) and $\nu >0$ scales leisure utility. I chose additive separability because it gives the leisure margin a marginal utility of its own, and therefore something to invert. The household solves

where $\beta$ is the discount factor and $𝖱$ is the risk-free gross return. The terminal condition is ${\mathrm{v}}_{T+1}\equiv 0$, so the household in the final period $T$ consumes all remaining resources; leisure then solves the static first-order condition (or its corners). Throughout, I assume standard constraint qualifications hold, so that interior solutions satisfy first-order conditions.^[6] The intermediate quantities the constraints define are market resources $m_t$ (incoming assets with their return, plus labor income) and end-of-period liquid assets $a_{t+1}$.

Although the two decisions are economically simultaneous, the dependence structure permits sequential solution. The labor-leisure choice determines market resources, and given those resources the consumption-saving choice determines liquid assets.

The sequential decomposition starts the period with the labor-leisure problem. To distinguish value functions at different stages, I introduce stage superscripts: the original problem has ${\mathrm{v}}_t\equiv v_t^0$, the value at the first decision stage (labor-leisure), and $v_t^i$ is the value after the decisions at stages $0,1,\dots ,i-1$. At this stage, the household observes the liquid assets $a_t$ it carries into the period and the wage offer ${\theta }_t$, choosing leisure to maximize the sum of current leisure utility and the continuation value from market resources $m_t$:

Once market resources are realized, the pure consumption-saving problem determines how to allocate $m_t$ between current consumption and liquid assets. The state space shrinks to a single dimension because the wage offer no longer matters:

No choice remains after consumption, so the period closes with a post-decision value: the discounted expectation of next period’s problem, taken over the wage draw the household does not yet know,

The sequential formulation follows the nested approaches of Clausen & Strub (2020) and Druedahl (2021) but composes EGM inversions without embedding optimization. Each choice is self-contained in a subproblem, with the structure chosen to minimize state variables at each stage. The sequential formulation is mathematically equivalent to the original joint problem by Bellman’s principle of optimality Bellman, 1957: the expectation operator appears only in the final subproblem, so no uncertainty resolves between stages and every decision uses the same information set. Were shocks to resolve between stages, households would gain information unavailable in the original formulation, changing the problem’s solution.

Limiting cases confirm that the decomposition nests simpler models. When $\nu \to 0$, leisure utility vanishes and the labor-leisure stage degenerates: the household works full time regardless of wages, recovering the standard consumption-saving EGM of Carroll (2006). When $\beta \to 0$, the fully impatient household consumes all resources immediately, recovering the static problem.

2.2Sequential Solution¶

EGM${}^n$ inverts both first-order conditions of the period in closed form, without grid search or root-finding, applying the standard EGM inversion to the consumption Euler equation and an analogous one to the leisure condition.

The consumption-saving stage is the standard EGM of Carroll (2006). I take it as given and pause only to fix notation: bracketed variables ($[\mathrm{a}]$, $[\mathrm{m}]$) denote exogenous grids, gothic letters ($𝔠$, $𝔪$) denote endogenous quantities built by inverting a first-order condition, and primes are ordinary derivatives. Inverting the Euler equation $u^{\prime }(c_t)={(v_t^2)}^{\prime }(a_{t+1})$ for consumption, which is possible because the marginal-utility map $u^{\prime }$ is strictly decreasing and hence one-to-one, gives consumption on the exogenous post-decision grid, and the budget identity returns the endogenous grid of market resources:

The pair $({𝔪}_t,{𝔠}_t)$ supports an interpolating consumption function for market-resources values off the grid. Through the envelope condition ${(v_t^1)}^{\prime }(m_t)=u^{\prime }(c_t)$, the same pair gives the labor-leisure stage the marginal value it needs without additional computation.

The leisure choice has a first-order condition of its own. Substituting $n_t=1-z_t$ and $m_t=𝖱a_t+{\theta }_t{n}_t$ puts market resources inside $v_t^1$, so the whole trade-off between leisure and labor income runs through that one function:

This equates the marginal utility of leisure to the opportunity cost of forgone labor income ${\theta }_t$, scaled by the marginal utility of consumption ${(v_t^1)}^{\prime }(m_t)=u^{\prime }(c_t)$.

The marginal utility of leisure inverts in closed form,^[7] so on an exogenous grid of market resources $[\mathrm{m}]$ and wage shocks $[\theta ]$ leisure follows directly:

Nothing in that step is special to consumption. The inverse marginal utility of leisure ${({\mathrm{h}}^{\prime })}^{-1}$ carries an exogenous grid of market resources and wages to an endogenous grid of incoming assets, exactly as ${(u^{\prime })}^{-1}$ carries post-decision assets to market resources one stage later. EGM${}^n$ chains the two.

Households with high market resources $m_t$ and low wage offers ${\theta }_t$, however, may find the unconstrained optimum violates the feasibility constraint $z_t\in [\mathrm{0,1}]$. When this occurs, projecting the solution onto the constraint boundary defines the constrained optimal function ${\widehat{𝔷}}_t([\mathrm{m}],[\theta ])$ as

The projection keeps leisure within the time endowment; in practice only the upper bound binds.^[8] Interpolation handles the resulting kink by sorting the endogenous grid and fitting piecewise segments on either side of it.

Labor follows as $1-{\widehat{𝔷}}_t(m_t,{\theta }_t)$, and for each ${\theta }_t$ and $m_t$ on the exogenous grid the endogenous grid of incoming assets is ${𝔞}_t(m_t,{\theta }_t)=(m_t-{\theta }_t(1-{\widehat{𝔷}}_t(m_t,{\theta }_t)\left)\right)/𝖱$.

By the envelope theorem, the marginal value the previous period needs is $\partial v_t^0/\partial a=𝖱{(v_t^1)}^{\prime }(m_t)$, which completes the backward step. At interior solutions this marginal value equals $𝖱{\mathrm{h}}^{\prime }(z_t)/{\theta }_t$. At a corner, however, only the first equality survives, because a Kuhn-Tucker inequality replaces the first-order condition.

The resulting endogenous grid for the labor-leisure problem is curvilinear rather than rectilinear.^[9]

Figure 1 collects the full solution; no numerical optimization produced any of the policies it shows.^[10] In the first period, the marginal propensity to consume falls from about 0.48 at the borrowing constraint to about 0.024 near $m_t=4.4$, then rises slowly to about 0.06 at higher wealth. The slow rise plausibly mirrors the decline in the household’s expected future labor, which reaches zero near $m_t=18$, shortly before the propensity levels off near $m_t=20$. Below $m_t=18$, saving most of an extra unit lets the household work less later in life. The labor-leisure inversion already produces a curvilinear grid, the simplest case of the meshes Section 4 interpolates. The two-asset problem of Section 3, where two state variables persist across stages, provides the more demanding test.

3A Two-Asset Liquid-Illiquid Household¶

The consumption-labor problem in Section 2 carries at most one post-decision state variable per stage. The canonical stress test for a sequential EGM step is the two-asset household block at the center of the heterogeneous-agent New Keynesian (HANK) literature Kaplan et al., 2018Auclert et al., 2021: a household holding liquid wealth alongside an illiquid account, moving funds between them at a convex adjustment cost. Two continuous state variables persist across the EGM stages, so both coordinates of the endogenous grid warp. The same structure, a second persistent stock alongside liquid wealth, appears in durable goods, human capital, and health investment problems; Online Appendix A develops the health investment model of Ludwig & Schön (2018) and White (2015) as a second worked example. I solve the exact household block of Auclert et al. (2021) at its published calibration, so that the benchmark runs on a model the macro audience already uses at scale. I treat only the household block and leave the New Keynesian aggregate machinery that closes the HANK economy to the literature cited above. Solver speed compounds in that setting: a general-equilibrium iteration can solve the household block thousands of times, each pass consuming the policies on the regular state grid, where the non-stochastic histogram method of Young (2010) carries the wealth distribution forward.

Nested methods conventionally handle the illiquid choice with a numerical first-order-condition search embedded in the backward step, exactly the numerical solve EGM${}^n$ is designed to remove; the decomposition below replaces it with a second closed-form inversion. Timing the two solvers measures what that substitution saves, net of the regrid and the regular-grid Euler inversion it adds.

3.1The Household Problem¶

A household enters the period with state $(b,a,e)$, where $b$ is incoming liquid wealth, $a$ is the incoming illiquid balance, and $e$ is idiosyncratic labor productivity following a discrete Markov chain with transition matrix $\mathrm{\Pi }(e^{\prime }|e)$. The letters follow the two-asset convention of the heterogeneous-agent literature. Following the Auclert et al. (2021) timing convention, $a$ and $b$ are post-decision balances: the choices this period, written $a^{\prime }$ and $b^{\prime }$, become next period’s state. Both balances pay out their gross returns, ${𝖱}_a=1+r_a$ on illiquid and ${𝖱}_b=1+r_b$ on liquid, labor income $y(e)$ is received, and market resources are $m={𝖱}_{b}b+{𝖱}_{a}a+y(e)$. The household chooses consumption $c\ge 0$ together with new balances $b^{\prime }\ge 0$ and $a^{\prime }\ge 0$ subject to the single within-period budget

The illiquid asset earns the higher return ($r_a>r_b$), but moving funds into or out of it costs $\mathrm{\Psi }$. Adjustment costs take the quadratic form of Auclert et al. (2021),

with offset ${\chi }_0>0$ and scale ${\chi }_1>0$. The cost vanishes when the household reinvests the illiquid interest ($a^{\prime }={𝖱}_{a}a$), is strictly convex in $a^{\prime }$, and scales with the post-return illiquid stock ${𝖱}_{a}a+{\chi }_0$, so that even small illiquid positions are costly to deplete. The recursive problem is

subject to the budget and the non-negativity constraints. Unlike the value-based curvilinear EGM of White (2015), the inversions below need only positive risk aversion, $\rho >0$, because each inverts a strictly positive marginal value rather than a utility level.

In one familiar limit, ${\chi }_1\to \infty$, the illiquid balance moves only at the gross return ($a^{\prime }\to {𝖱}_{a}a$) and the problem reduces to a one-asset consumption-saving model in $(b,e)$. A fixed cost paid on any adjustment ($a^{\prime }\ne {𝖱}_{a}a$) would generate a region of inaction and a non-convexity; that case leaves the smooth-EGM class and belongs to the G2EGM upper-envelope methods, the natural extension noted in the conclusion.

3.2Stage Decomposition and Backward Solution¶

The joint problem splits into three stages: an illiquid choice, a consumption choice, and a post-decision expectation. Under the stage convention of Section 2, $\mathrm{v}\equiv v^0$ is the first decision stage, $v^1$ is the value entering consumption, and $\mathrm{w}$ is the post-decision value, which carries its own letter here, following Powell (2011), because both inversions read off it. With two persistent states, derivatives are partial: they appear as subscripts on stage values ($v_{\ell }^1\equiv \partial v^1/\partial \ell$) and, because $\mathrm{w}$ carries no stage label, as superscripts on $\mathrm{w}$ (${\mathrm{w}}^{b^{\prime }}\equiv \partial \mathrm{w}/\partial b^{\prime }$). Productivity $e$ passes through every stage unchanged and resolves only in the expectation; every stage value therefore carries $e$ as a state. The expectation again appears only in the post-decision stage, so the equivalence argument of Section 2 applies unchanged. Where the consumption-labor stages shed a state variable apiece, here the state count stays at three throughout. Backward induction visits the stages in reverse, and both decision stages turn out to be EGM inversions, linked by envelope conditions.

The post-decision stage integrates over next period’s productivity:

Stage 1: The Consumption EGM with the Constrained Segment¶

Given the liquid resources $\ell$ left after the illiquid choice, the consumption stage chooses $c$:

The illiquid choice $a^{\prime }$ is already made when consumption is chosen, so at each fixed $a^{\prime }$ this stage is the consumption EGM of Section 2, run along the $b^{\prime }$ axis and resting on the same separable utility (Section 5). The stage reads the post-decision marginals ${\mathrm{w}}^{b^{\prime }},{\mathrm{w}}^{a^{\prime }}$ on the exogenous $(b^{\prime },a^{\prime })$ grid, inverts the Euler equation,

and returns $v_{\ell }^1,v_{a^{\prime }}^1$ on a regular grid of liquid resources $\ell$. The liquid borrowing constraint $b^{\prime }\ge 0$ is handled in the manner of Carroll (2006), by prepending the single point $(\ell =0,c=0,b^{\prime }=0)$ to the endogenous mesh. With $b^{\prime }$ starting at zero, the first bracket lies on the line $c=\ell$, and linear interpolation reproduces the constrained segment exactly. Through the envelope the illiquid marginal passes unchanged, $v_{a^{\prime }}^1={\mathrm{w}}^{a^{\prime }}$ evaluated at the Stage 1 optimum.

Stage 0: The Illiquid EGM with Closed-Form Inversion¶

The illiquid stage, given $(b,a,e)$, chooses $a^{\prime }$ and converts market resources into liquid resources:

EGM${}^n$ departs from NEGM at this stage, where the first-order condition for $a^{\prime }$ balances the marginal value of liquid against the marginal value of illiquid net of the adjustment-cost wedge:

The right-hand side is read directly from Stage 1’s regular-grid envelopes. A nested method searches over $a^{\prime }$, but EGM${}^n$ holds the gridded choice $a^{\prime }$ fixed and recovers the incoming stock $a$. Here the enabling structure is the invertible transition marginal of Section 5: the stage has no utility term of its own, and ${\mathrm{\Psi }}_1$ is strictly decreasing in the unknown $a$ because ${\mathrm{\Psi }}_{12}<0$. For the quadratic cost, ${\mathrm{\Psi }}_1(a^{\prime },a)={\chi }_{1}u$, where $u\equiv (a^{\prime }-{𝖱}_{a}a)/({𝖱}_{a}a+{\chi }_0)$ is the standardized adjustment wedge, and

where ${\mathrm{\Psi }}_{ij}$ denotes the corresponding second partial derivative.

Define the wedge target $\tau \equiv v_{a^{\prime }}^1/v_{\ell }^1-1$; the inversion is

a single multiplication and a ratio per grid point. The remaining endogenous quantities follow from the two accounting identities: endogenous market resources from the liquid-resources definition, $𝔪=\ell +a^{\prime }+\mathrm{\Psi }(a^{\prime },𝔞)$, and the endogenous liquid balance from the market-resources definition, $𝔟=(𝔪-{𝖱}_a𝔞-y(e))/{𝖱}_b$. There is no root-find and no iteration; the entire illiquid choice is closed form.^[11]

The inversion produces a curvilinear $(𝔟,𝔞)$ mesh while retaining the exogenous grid’s adjacency. Section 4 states the two geometric conditions ENGINE places on such a mesh, both of which I verify directly on every reported benchmark mesh. After ENGINE carries the illiquid choice to the regular state grid, the remaining policies follow by composition: the implied liquid resources pull the Stage 1 consumption function back to $(b,a,e)$, evaluated through the same Euler inversion, and the liquid choice is the budget residual (Online Appendix B). Figure 2 shows the result.^[12] The inversion bends the mesh most where the optimal adjustment is largest. Yet mesh lines of constant $\ell$ and constant $a^{\prime }$ never cross.

Boundaries: Two Economic Edges and One Grid-Design Rule¶

The inversions handle every interior solution; the boundaries come in two kinds. The economic boundaries are the non-negativity constraints, neither of which costs a subproblem. The liquid edge ($b^{\prime }=0$) is the consumption EGM’s own borrowing constraint, handled inside Stage 1 by the prepended constraint point. The illiquid edge ($a^{\prime }=0$) reduces the cell to a Stage 1 evaluation at the implied liquid resources $\ell =m-\mathrm{\Psi }(0,a)$, so the corner policy is read off Stage 1’s own $a^{\prime }=0$ slice, and each cell picks between the interior inversion and the corner by the sign of the unconstrained policy. With strictly concave utility and a convex adjustment cost, the stage policies are single-valued and continuous, so no discrete jumps arise and the solver never requires an upper envelope.

The grid boundary imposes a design rule with no analogue in one-block solvers. Sequential decomposition materializes the intermediate state $\ell$ on a regular grid between the stages, and that grid must cover whatever the second inversion can request. Inherited liquid wealth does not bound it: a household rich in illiquid assets finances liquidity by withdrawing, and because the marginal value of illiquid wealth is non-negative the first-order condition bounds the withdrawal wedge by $u^{\ast }\ge -1/{\chi }_1$, so $\ell$ reaches at most

in closed form, where $\overline{b}$ and $\overline{a}$ are the largest nodes of the liquid and illiquid grids and $\overline{y}$ the highest income realization.

Truncating the $\ell$-domain to the liquid state grid, the natural shortcut, strips mesh coverage at high illiquid wealth without warning and forces the realign step to extrapolate there. On this calibration the iteration then converges to a budget-infeasible pseudo-solution in the affected band, which a constraint-valid Euler diagnostic exposes at once (Section 4).^[13]

3.3Comparison to Nested EGM¶

I compare the two solvers on total time to a converged household solution, at matched grids, matched tolerance (10^-9), and matched linear interpolation. All timings are net of one-time compilation cost, best of three repetitions from a warm start.^[14] Table 1 reports both solvers’ times. The solvers share the consumption-stage EGM, the constrained segment, the post-decision integration, and the fixed-point iteration. NEGM solves the illiquid first-order condition directly on the regular state grid by a safeguarded bracketing iteration whose relative bracket tolerance must sit two orders below the outer norm, at 10^-11. With the bracket loosened to 10^-9, per-node root-find noise keeps the fixed-point iteration from meeting its own norm and the solve hits the iteration cap, so the comparator converges only with the tighter bracket (the archive’s inner-tolerance table records the check). EGM${}^n$ instead applies the closed-form inversion on its exogenous grid, then pays for a general ENGINE regrid and a regular-grid Euler inversion. The two backward steps are parallelized identically, so the comparison is matched at one thread and at sixteen.^[15]

NEGM and EGM${}^n$ converge in a comparable number of backward iterations (at $(\mathrm{40,60})$, $(\mathrm{80,120})$, and $(\mathrm{160,240})$, NEGM takes 473, 481, and 487 iterations and EGM${}^n$ takes 395, 410, and 424). The timing gap comes from the illiquid step, where EGM${}^n$ makes one closed-form evaluation and NEGM an iterated sequence of interpolated envelope evaluations. EGM${}^n$ still wins after paying for the general Stage-0 ENGINE regrid and regular-grid Euler inversion that NEGM does not need, since NEGM’s root-find returns the policy on the regular $(b,a)$ grid. Serially the factor is about thirteen at all three grids. At the largest, one serial NEGM solve takes just under three minutes against fourteen seconds for EGM${}^n$. At matched sixteen threads the factor grows from about 3.7 to 7.9 with the grid size. The root-find carries enough arithmetic per node to scale almost linearly with threads, whereas the EGM${}^n$ step is so cheap per node that at small grids its solve time is probably dominated by the fixed cost of coordinating the threads themselves. I expect the matched-thread ratio to keep climbing toward the serial factor at the grid sizes where computing time matters.

The two solvers also agree on the solved margins. Median relative distance between the converged policies stays below half a percent for both the liquid and the illiquid balance at every timing grid. At the $192\times 160$ policy-benchmark grid the two asset margins agree to 0.16 and 0.29 percent. On the model’s own liquid Euler equation, the comparison depends on where I take the residual. Off the grid, where both converged policies are interpolated and the residual puts the two methods on equal footing, EGM${}^n$’s median is the smaller at six of the seven reported grids (three timing, three diagnostic, and the policy benchmark), trailing only at the $(\mathrm{80,120})$ timing grid, by a quarter of a percent. At the three diagnostic grids the medians differ by less than fifteen percent, and EGM${}^n$’s advantage grows to factors of 1.4 and 1.6 at the two largest grids overall, the $(\mathrm{160,240})$ timing grid and the policy benchmark. On the regular grid, EGM${}^n$’s residual runs one to two orders of magnitude below NEGM’s, about $1.5\times {10}^{-5}$ against $3.9\times {10}^{-4}$ at $(\mathrm{40,60})$, because EGM${}^n$ delivers regular-grid consumption as an Euler-inversion output where NEGM reads it off an interpolated policy (Online Appendix B).^[16]

Each stage is itself a Bellman problem; every iteration of the backward step therefore also produces the marginal values $v_{\ell }^1$ and $v_{a^{\prime }}^1$ of the consumption stage on a regular grid. Their ratio minus one is the wedge $\tau \equiv v_{a^{\prime }}^1/v_{\ell }^1-1$, the adjustment-cost shadow price of converting illiquid into liquid wealth. A modeler can read it off to answer questions that turn on this object (how the adjustment friction varies over the wealth distribution, where it binds), with no auxiliary computation. NEGM shares this property (its root-find runs on the same stage envelopes), but solvers that treat the period as a single block with two controls, the G2EGM family among them, deliver policies without ever constructing the intermediate-stage values. Recovering the same objects on the whole stage grid, beyond the optimal choices, requires re-deriving the stage decomposition that EGM${}^n$ starts from.

4Multivariate Interpolation on Curvilinear Grids¶

EGM’s efficiency comes from working on exogenous grids on which calculations are straightforward, then using the inverted Euler equation to recover endogenous state variables. That recovery, however, gives up the grid’s rectilinear geometry, so standard multilinear interpolation, which assumes each dimension varies independently, cannot handle the curvilinear grid that results.

ENGINE relies on the index structure that endogenous grids inherit from exogenous grids despite geometric warping. It reduces multidimensional interpolation to a sequence of one-dimensional interpolations along these inherited indices, because row $k$ remains row $k$ after warping and neighbors remain neighbors, even when geometric regularity (uniform spacing, rectilinear alignment) is lost.

Consider a curvilinear grid indexed by $j\in \{1,\dots ,J\}$ and $k\in \{1,\dots ,K\}$, where each “row” (fixed $k$) contains $J$ points. The curvilinear interpolation of White (2015) (the curvilinear method below) locates the containing quadrilateral cell via sector-walking, then solves for normalized coordinates within that cell (a quadratic equation in 2D, Newton iteration in higher dimensions). ENGINE sidesteps both operations: for each of the $K$ rows it asks where a vertical line at $x^{\ast }$ intersects that row, a one-dimensional interpolation problem that finds where $x^{\ast }$ falls among the $J$ $x$-coordinates and interpolates the $y$-coordinate to produce an intermediate point. A final 1D interpolation across these $K$ intermediate points yields the answer. In implementation, a single query touches only the $O(\log K)$ rows its binary search probes, and the dense regrid reuses brackets across sorted targets.

4.1The ENGINE Algorithm¶

ENGINE works in index space, where bracketing is efficient, and can map back to physical coordinates because each endogenous point $(x_{jk},y_{jk})$ keeps the indices $(j,k)$ of the exogenous post-decision grid point that generated it through the EGM inversion. This index-to-mesh correspondence is a homeomorphism under conditions developed below, a result the appendix states and proves (Theorem 1). Figure 3 makes the correspondence concrete. Because the mapping preserves neighborhoods, a bracket search on the rectilinear index grid is a valid search of the warped physical mesh.

The full two-pass listing and boundary-extrapolation behavior are recorded in the Appendix. Figure 4 traces one query through both passes.

For a query point $(x^{\ast },y^{\ast })$, the interpolated value is

where ${\text{𝙴𝙽𝙶𝙸𝙽𝙴}}_k$ denotes Pass 2’s one-dimensional linear interpolation across the row index $k$, and ${\hat{y}}_k$ and ${\hat{f}}_k$ are the intermediate coordinates and values from Pass 1.

On a warped mesh the two passes do not commute. Walking the regular axis first moves the brackets less between rows.

For sorted targets, each search resumes from the previous target’s bracket and never moves backward, so total work is proportional to the number of targets plus their total bracket displacement.

4.2When ENGINE Applies¶

The pure consumption-saving problem illustrates the basic phenomenon. An exogenous grid of post-decision liquid assets $[\mathrm{a}]$ maps through the inverted Euler equation to an endogenous grid of market resources $[\mathrm{m}]$ with different spacing, the nonlinearity of marginal utility making the mapping non-uniform and binding constraints adding kinks. In one dimension, non-uniform linear interpolation handles the distortion. Higher-dimensional problems, however, inherit the warping in each dimension simultaneously. How much of the grid’s structure survives then determines the appropriate interpolation method.

Previous approaches to this problem include Delaunay triangulation Ludwig & Schön, 2018 and curvilinear interpolation White, 2015. Delaunay triangulation constructs a mesh of simplices over the scattered points and performs barycentric interpolation within each simplex. Although general and robust, triangulation is computationally expensive. Construction from unstructured point sets requires at least order $N\log N$ operations in the algebraic decision-tree model, a lower bound rather than an implementation artifact: sorting $N$ numbers reduces to triangulating $N$ points on a parabola Preparata & Shamos, 1985. The triangulation must also be recomputed when grids change during iterative solution procedures. The curvilinear method exploits the grid structure to map irregular quadrilateral sectors to the unit square for bilinear interpolation; sector location uses a visibility walk algorithm, and coordinate calculation in 2D uses a closed-form quadratic solution (3D and higher require Newton iteration).

ENGINE requires structural conditions on the endogenous grid. Analytic invertibility of a stage first-order condition (Proposition 1, Proposition 2) is not enough on its own, because the recovery map defined in the Introduction can still fold when continuation is nonconcave or when cross effects reverse orientation. Discrete choices and nonconcave regions are the leading exclusions.

ENGINE imposes two requirements, both geometric: cells cannot fold, and each coordinate must increase with its own index. A third property, a bound on bracket displacement between adjacent rows, governs how far ENGINE must search. The implementation verifies every windowed bracket and falls back to a full row search on a miss, so the search margin affects only speed.

Nonlinear EGM mappings can distort cells severely enough to invert them. The fold-free condition, the classical validity condition from structured mesh generation Knupp, 2001, ensures each cell maintains proper orientation.

Four-corner positivity rules out a fold; Lemma 1 shows that it delivers $𝒥(s,t)>0$ throughout the cell.^[17]

Fold-freeness ensures local cell validity, but ENGINE also requires that bracketing along rows and columns be well-defined. IPM ensures coordinates increase monotonically along each index direction.

Together, fold-free and IPM ensure the grid mapping is a homeomorphism, a continuous bijection with continuous inverse from index space to physical space, so the sequential brackets are well-defined within the interpolation region.^[19] Efficiency also depends on how far a query’s bracket can move between adjacent rows. The per-row maximum of that drift, the constant $\alpha$, determines the search window, and a smaller $\alpha$ reduces search work. Economically, $\alpha$ tends to grow where policies change rapidly relative to grid spacing, near occasionally binding constraints or under high value-function curvature, and when one axis’s spacing is much finer than the other’s. Refining one axis alone creates exactly that mismatch, so refinement should proceed at a common rate with comparable spacing across axes. The constant governs the single-query search window. By contrast, the cost of the dense regrid in the backward step depends on how far the brackets move in total across its sorted targets.

After each solve I compute the empirical estimate $\hat{\alpha }$, the largest bracket displacement observed between adjacent rows over a dense set of test points, as a mesh diagnostic. On the health-model meshes in Online Appendix A, $\hat{\alpha }$ ranges from 10 to 16, well above the default margin of two. Margins matched to those measurements avoid the fallback searches, and ENGINE completes every reported solve.^[20]

I verify directly whether an inverted stage’s recovery map delivers a qualifying mesh. On the two-asset block, exhaustive four-corner and row-and-column monotonicity tests confirm both conditions on every cell of the converged mesh at every reported benchmark refinement (Section 3). On the health model of Online Appendix A, every reported mesh of the two structured methods satisfies IPM, and the only folds on those meshes are sliver cells at the lowest health node. In two high-mortality calibrations of that appendix’s warp sweep, under half a percent of evaluation points fall in those cells, and in the other calibrations none do; the reported residuals include every final query inside one. The measured $\hat{\alpha }$ stays between two and three on every reported two-asset mesh, well below the health meshes’ ten to sixteen. Yet where it reaches three it exceeds the default margin of two, which the containment fallback absorbs. On the health meshes I set the margin from the measurement.

Occasionally binding constraints introduce kinks where the recovery map is not smooth. The implementations build the constrained and interior branches separately; on the two-asset block the checks above then verify the assembled mesh, kink-straddling cells included. In the health model of Online Appendix A the kink falls on a mesh line, so no cell straddles it, and the cells next to it are fold-free in every reported mesh of the two structured methods. Discrete choices instead create jumps that no continuous interpolant should cross. Stages with discrete choices therefore belong to upper-envelope methods such as G2EGM or DCEGM.

4.3Accuracy and Complexity¶

On rectangular grids the two passes reduce algebraically to standard bilinear interpolation, whose error for twice continuously differentiable functions is second order in the mesh width Boor, 2001. I measure accuracy on warped meshes directly. Under successive grid doublings the robust off-grid Euler residual falls by about a factor of five per doubling, and the converged policy’s distance to a high-resolution benchmark falls about threefold per doubling (Online Appendix B).

For the dense regrid onto a $P\times Q$ target grid in each backward step, each target costs at most the single-query bound, so the regrid is $O(PQ\log J\cdot \log K)$. Sorted target axes, however, let each bracket search resume where the previous one ended, so in the scaling benchmark the cost falls far below that bound. On meshes with monotone row shifts, whose measured $\hat{\alpha }$ nevertheless rises with mesh size, doubling $J=K$ from 32 to 1024 gives $N^{0.98}$ scaling in $N=JK$, about ten nanoseconds per output value at the largest mesh.

The two-dimensional implementation needs neither a construction step nor Newton iteration.^[21]

Choosing the Interpolation Order¶

Practitioners should walk the regular axis first, the one mapped to a strictly monotone endogenous coordinate with smaller bracket displacement. In the two-asset Stage 0 recovery domain, this is the liquid-resources axis $\ell$; adjustment costs warp the incoming illiquid coordinate.

One limiting case makes the second pass trivial. If the second coordinate is aligned across rows ($y_{jk}=y_k$) with the shared ordinates increasing ($y_{k+1}>y_k$), and every row satisfies $x_{j+1,k}>x_{jk}$, the aligned family supplies column monotonicity while the ordered rows supply row monotonicity and positive cell orientation.^[22]

On warped meshes, walking the regular axis first reduces bracket displacement and the effective constant $\alpha$. Reversal keeps IPM in the tested meshes but raises search work toward the $O(\log J\cdot \log K)$ worst case. It also changes the interpolant because sequential one-dimensional interpolation does not commute across warped coordinate lines. Re-solving the two-asset block in reverse order changes total solve time by less than ten percent at every grid, while the median relative consumption difference stays below 0.15 percent and robust Euler medians remain within two percent (Online Appendix B). In this application, order governs the worst-case search cost rather than converged accuracy.

Higher Dimensions¶

The sequential construction extends algorithmically to higher dimensions through a single one-dimensional pass per dimension. Its efficiency would require controlling bracket drift along each added dimension; the implementation and quantitative claims here are two-dimensional.

4.4Comparison with Alternative Interpolants¶

The three approaches to multidimensional EGM interpolation differ in how much of the endogenous grid’s structure they exploit. On a fold-free mesh satisfying IPM, the endogenous grid is the image of the rectilinear exogenous grid under a homeomorphism (Theorem 1), so the rectangular index lattice survives intact and every row and column remains sorted. Delaunay triangulation uses none of this structure; the curvilinear method uses adjacency (which cell neighbors which) but not order; ENGINE uses both:

Among these three methods, construction cost rises as the grid requirements fall, because a method that assumes less structure spends computation reconstructing it. ENGINE’s construction is free because EGM supplies the index lattice. The curvilinear method recomputes its per-cell sector data and the triangulation rebuilds order at each backward step. In the scaling benchmark the dense matched regrid runs at output size, while unstructured triangulation carries a construction cost of at least order $N\log N$ before its first query. A lattice-aware $O(N)$ triangulation would rely on the same index structure ENGINE exploits and would still be undefined outside the mesh’s convex hull.

Figure 5 contrasts the geometric work each method performs to reach the same query point.

The solve-time ratio understates the per-query advantage because interpolation is one of several cost components of a backward step, alongside the expectation step and the EGM inversions, so even a large per-query gain translates into a bounded reduction in total solve time. The asymptotic comparison in Table 2 covers only the interpolation component, where ENGINE can at best tie triangulation. Any comparison-based point location needs at least order $\log N$ comparisons to find the containing cell among $N$ candidates; Delaunay point location attains that floor, and ENGINE’s $O(\log J\cdot \log K)$ sits at most a row-search factor above it. The two passes multiply their logarithms rather than add them. Locating the query between rows takes $O(\log K)$ probes, each requiring its own $O(\log J)$ bracket search within a row.

Unlike the solver comparison of Section 3, the full-solve diagnostic of Online Appendix B holds the EGM${}^n$ solver fixed and swaps only its Stage-0 regrid between the two structured interpolants, scoring both by total solve time to convergence and by a robust Euler error that is valid at constraint-bound queries, since it compares consumption with the lesser of the Euler-implied and constraint-feasible levels.

On accuracy the diagnostic (Table 3) finds ENGINE and the curvilinear method interchangeable, agreeing on the robust Euler error within two percent at every grid and quantile. On speed, ENGINE completes the same solve roughly three to four and a half times faster, widening with the grid. The difference is pure interpolation-method cost, including the comparator’s per-step construction, since everything but the regrid step is identical across the solves. Measured directly at the converged meshes, one realign of the regular target grid costs ENGINE 0.02 to 0.26 milliseconds across the three grids against the curvilinear method’s 0.54 to 4.7, factors of 18 to 27. Both solves run this step on one thread, even in the sixteen-thread diagnostic, so these timings apply to it unchanged. Realigning alone accounts for more than four fifths of the per-iteration difference between the two solves. The asymptotic entry in Table 2 does not capture this per-call cost. At the converged fixed point, exactly one regular Stage-0 target lies outside the endogenous mesh at each grid. Online Appendix B reports the target totals and outside shares, the protocol, and boundary-sensitivity readings. Extending both axis maxima by half changes medium-grid consumption on the common interior set by 0.05 percent at the median query and by 4.3 percent at the worst.

The two-asset diagnostic includes only the two structured interpolants. Online Appendix A adds the third alternative, Delaunay triangulation, on the health investment model where it and the curvilinear method were originally proposed and tested. I use that model because it favors those two methods, so any advantage ENGINE shows there is conservative. There EGM${}^n$ with ENGINE solves 2.8 to 4.7 times faster than the simultaneous solver with the curvilinear interpolation of White (2015), and swapping only the interpolant on that solver reproduces the bulk of that speedup. Across nine calibrations chosen to bend the endogenous meshes hardest, the triangulation’s worst consumption residual lies between 0.25 and 0.66, while ENGINE’s never exceeds 0.20. Against the curvilinear method, however, the tails split: ENGINE’s 95th-percentile consumption residual is smaller at seven of nine calibrations, while its maximum is larger at six of nine.

On fold-free meshes that satisfy IPM, ENGINE matches the curvilinear method’s accuracy at lower cost; when IPM fails, ENGINE does not apply. Section 5 asks which economic problems produce the qualifying stages and recovery maps.

5When the Method Applies¶

When does EGM${}^n$ apply? The consumption-labor problem (Section 2) and two-asset problem (Section 3) illustrated two enabling structures: separable utility functions and explicit inversion of a transition marginal followed by state recovery. In economic terms, separability means the marginal utility of one choice does not depend on the level of another, while the second structure recovers a choice or incoming state from how that choice changes the states carried forward. Below, invertible transition marginal is shorthand for this full marginal-inversion and state-recovery condition.

5.1Splitting the Problem into Subproblems¶

Decomposition begins by counting the problem’s distinct controls. A problem with $n$ control variables might decompose into $n$ subproblems. Do not double-count: the consumption-saving choice ($c+a=m$) is one decision, not two. Similarly, labor-leisure is a single choice despite involving two variables.

Each control variable requires a structure that permits EGM inversion. Two broad sufficient structures are separable utility and an invertible transition marginal in the shorthand just defined. More generally, an interior stage qualifies when its first-order condition has an explicitly evaluable, single-valued inverse, the target lies in its range, the transition identities recover the remaining endogenous state, and feasibility and curvature confirm the candidate as the unique optimum. Where a stage fails these conditions, it uses standard optimization, and the gains accrue only from the stages that invert.^[23]

The ordering of subproblems should shed state variables early. When subproblems affect disjoint state dimensions, the ordering is immaterial. When they share one, however, poor sequencing propagates unnecessary state variables through later stages. In the consumption-labor problem, placing labor-leisure first resolves the wage rate before the consumption stage, keeping that subproblem one-dimensional. Choosing consumption first would force the labor decision to track both liquid assets and wages, doubling its dimensionality.^[24]

5.2Formal Conditions¶

The formal statements below use uppercase calligraphic functions ($\mathrm{V}$, $\mathrm{W}$, $\mathrm{U}$, $\mathrm{T}$, $\mathrm{G}$) for generic value, continuation, utility, within-stage transition, and next-period state functions, and $\mathrm{\Gamma }$ for the feasible-choice correspondence, the set of controls the constraints permit at a given state, with generic state variables $x$, $y$, $s$ in place of the model-specific variables: consumption $c$, assets $a$, and cash-on-hand $m$ of preceding sections; the notation reference of the Appendix collects the generic-stage symbols.

Separable Utility¶

Once the problem is split into subproblems, the Endogenous Grid Method applies to each qualifying subproblem during backward iteration from the terminal period.

Consider a generic scalar subproblem with differentiable utility, transition, and continuation-value functions:

where $\mathrm{W}(y)=\beta 𝔼[{\mathrm{V}}_{+1}(y)]$ is the continuation value.^[25] For an interior solution, the first-order condition is

For a concave scalar stage with an interval constraint, a corner solution is obtained by projecting the interior candidate onto the feasible interval, as Section 2 demonstrates for leisure.

The question is when this first-order condition can be inverted to yield the optimal control directly. Proposition 1 gives the standard case used in Section 2: separable utility with a transition linear in consumption or labor.

Each hypothesis does distinct work in the displayed formula. Because $x$ is recovered endogenously, $\varphi$ may not depend on $x$. Without strict concavity the first-order condition can hold at a local minimum or at one of several stationary points, so the formula could return a candidate that is not the optimum, and nothing in the step would flag it. Without a single-valued inverse the step must select among branches, a selection the one-shot formula does not provide. If $k$ depended on $x$, the right-hand side would contain the state the formula recovers, and the displayed inversion would not close.

Standard budget constraints in consumption-saving models satisfy the linearity assumption $\frac{\partial \mathrm{T}}{\partial a}=k$. For example, $a=m-c$ implies $\frac{\partial a}{\partial c}=-1$, so $k=-1$. Similarly, the labor income constraint $m=𝖱a+\theta \cdot n$ has $\frac{\partial m}{\partial n}=\theta$. The wage shock may enter $k$ because it is observed before the stage is solved; the unrecovered endogenous state may not.

Computationally, the inversion replaces a numerical root-find with a single function evaluation at each grid point.

To see when the method fails, consider non-separable CES preferences (Cobb-Douglas within the period, so the elasticity of substitution is one) $\mathrm{U}(c,z)=\frac{{(c^{\xi }{z}^{1-\xi })}^{1-\rho }}{1-\rho }$, where $\xi \in (\mathrm{0,1})$ is the weight on consumption. The first-order condition for consumption involves leisure: $\frac{\partial \mathrm{U}}{\partial c}=\xi c^{\xi (1-\rho )-1}{z}^{(1-\xi )(1-\rho )}$. Inverting for $c$ requires knowing $z$, and the symmetric condition for $z$ requires knowing $c$. The first-order conditions are coupled, so neither can be solved in isolation by EGM inversion, and the stages do not separate. Homotheticity, however, still allows a closed-form solution here: given the wage, the intratemporal condition makes leisure proportional to consumption, and substituting it back leaves a consumption marginal utility that inverts in closed form, the triangular case of Iskhakov (2015). Without such a substitution, coupled conditions require a nested numerical solve (NEGM) or upper envelope techniques (G2EGM).^[26]

Invertible Transition Marginals¶

Some subproblems have no within-period utility at all, only transitions affecting multiple post-decision states. The illiquid-choice stage of Section 3 is such a case, since moving funds between the accounts has no utility of its own and matters only through the two quantities it sets, the liquid resources left after the adjustment cost and the new illiquid balance. Proposition 2 shows that EGM still applies under the invertible-transition-marginal shorthand defined above, even without separable utility. Consider a problem with two endogenous state variables and two post-decision states, where the first transition may also carry the second state:

where the differentiable continuation value is

The first-order condition becomes

In case (i), additive separability keeps unrecovered states out of the control derivative. Case (ii) permits state-control coupling because the control itself sits on the exogenous grid, leaving only the unrecovered state in the first-order condition.

The literature offers two neighboring standards for when EGM applies. White (2015) characterizes when EGM applies to the period problem as a whole, requiring that the endogenous transition factor through a post-decision state that absorbs the constraints (whether a choice was feasible can be read off the post-decision state alone), with the controls entering that state through a full-rank map. Under those conditions the expectation is computed once per gridpoint. But the first-order system they leave behind may still demand a multidimensional root-find. That residual root-find disappears under the “triangular” transition structure of Iskhakov (2015). Case (i)'s additive separability is that structure, under which the system solves by substitution with no root-finding, though Triangular EGM requires the entire multidimensional first-order-condition system to be triangular. The stage conditions ask less than triangularity and more than the transition decomposition alone. Each stage’s system must collapse to a scalar closed-form inversion. Meeting the stage requirement does not imply meeting the system-level one. In the two-asset application of Section 3 the stage conditions hold, but the full period first-order-condition system is not triangular.

The health-investment stage of Online Appendix A is case (i); the illiquid stage of Section 3 is case (ii). ENGINE’s separate geometric requirements are fold-freeness and index monotonicity, which Section 4 tests cell by cell on every reported mesh. Neither EGM${}^n$ nor ENGINE applies under nonconcave continuation, which breaks the curvature check that confirms the optimum and can fold the mesh. Under a discrete choice the policy is discontinuous, so the single-valued inversion the stage conditions require fails. The grids such a choice produces need an upper envelope rather than structured interpolation.

5.3Classification of Solver Classes¶

The four solver classes differ along three problem features identified by the EGM literature: whether each stage has additively separable utility, whether its transition marginal is invertible in the sense of Proposition 2, and whether the stage has an occasionally binding constraint. Table 4 maps each (feature, solver) cell to the operation that solver applies: an inversion recovers the policy in closed form from the first-order condition, a nested root-find wraps an outer one-dimensional search around the next stage’s solution, and an upper envelope adjudicates across candidate constraint regimes.

In the two-asset example, EGM${}^n$ evaluates an explicit inverse where NEGM nests a root-find.

For $n$ independent scalar bounds under this sequential structure, a simultaneous upper envelope considers up to $2^n$ binding regimes, while the decomposition handles two branches in each of $n$ stages. At $n=2$ the counts coincide, so the Section 3 timing comparison does not measure this distinction. This is an architectural count, not a general complexity theorem: it requires that no cross-stage constraint regime be carried forward.

For a new model, the conditions work as a checklist: count the controls, order the stages to shed state variables early, and test each stage against the conditions stated at the start of this section, of which Proposition 1 and Proposition 2 are sufficient cases. A stage that passes inverts in closed form, while one that fails keeps standard optimization. ENGINE then interpolates any resulting mesh that is fold-free and index-monotone.

6Conclusion¶

Economically simultaneous decisions need not be computationally joint. Because both of the consumption-labor problem’s first-order conditions invert analytically, EGM${}^n$ needs neither grid search nor root-finding at either stage. The leisure margin receives the same closed-form EGM inversion that Carroll (2006) applied to consumption. The two-asset block replaces the nested method’s embedded numerical search with a second closed-form inversion. Section 3 documents the speedup at matched grids and tolerance: the serial gap is roughly a factor of thirteen, and sixteen matched threads still leave EGM${}^n$ 3.7 to 7.9 times faster. There the two solvers agree to a median of under half a percent on both asset margins; EGM${}^n$ carries the smaller regular-grid liquid Euler residual.

With EGM${}^n$ held fixed and only the Stage-0 regrid swapped, ENGINE cuts full-solve time by a factor of roughly three to four and a half with robust Euler errors within two percent. On the health model, EGM${}^n$ with ENGINE solves 2.8 to 4.7 times faster than the simultaneous solver with the curvilinear method.

EGM${}^n$ applies only when a stage has an explicitly evaluable, single-valued first-order-condition inverse and feasibility and curvature confirm the optimum. ENGINE separately requires a fold-free, index-monotone mesh; analytic inversion alone does not guarantee one. Both conditions are instead checked directly, cell by cell, on every reported mesh.

Stages that still require numerical optimization gain only the decomposition itself, which shrinks the grid search or root-find to the stage’s own lower dimension, the observation on which NEGM builds Druedahl, 2021. Discrete choices and severe non-convexities break the monotonicity both algorithms require and belong to G2EGM Druedahl & Jørgensen, 2017; in the paper’s consumption-labor and two-asset applications, the leisure and asset bounds are handled on separate branches that preserve the required ordering. Repeated estimation solves can accumulate the per-solve savings when their household blocks retain this stage structure. The two-asset consumption model of Kaplan & Violante (2014), the setting closest to this paper’s own application, turns on exactly this kind of interacting liquid and illiquid choice; life-cycle wealth accumulation Cagetti, 2003 and household balance-sheet dynamics Fagereng et al., 2021 pose the same structure.

Both algorithms are implemented in the Econ-ARK HARK toolkit Carroll et al., 2018, and the replication archive reproduces every figure and table in this paper.

Two extensions look tractable, one widening the problem class and one economizing the grid. A hybrid would invert the separable stages and hand only the stage carrying a discrete choice to G2EGM’s upper envelope, which would reach problems that mix the two structures without paying the envelope’s cost at every stage. Adaptive refinement would concentrate grid points where the policy bends, a location the warp sweep of Online Appendix A already measures. Can the index correspondence be rebuilt cheaply enough between refinements to keep ENGINE’s arithmetic advantage?

7Appendix¶

7.1The Grid Homeomorphism Theorem¶

The direct four-corner calculation verifies fold-freeness on every cell of each finite mesh reported in Section 3, kink-straddling cells included.

7.2The ENGINE Algorithm¶

Section 4 develops ENGINE and states its probe costs; this section records the single-query listing and the boundary behavior. The dense regrids in the benchmarks run the same two passes with brackets carried forward across sorted targets, as Section 4 describes.

Algorithm (ENGINE Two-Pass Interpolation). Inputs: grid $\{(x_{jk},y_{jk},f_{jk})\}$ for $j=1,\dots ,J$ and $k=1,\dots ,K$; query point $(x^{\ast },y^{\ast })$. Output: interpolated value $f(x^{\ast },y^{\ast })$. The listing uses subscripts $l$, $m$, and $h$ for the low, mid, and high indices. The routine $\mathrm{BracketSearch}$ verifies that the bracket it returns contains $x^{\ast }$ and falls back to a full row search on a miss.

Pass 1: binary search on rows.

1. Set $k_l\leftarrow 1$ and $k_h\leftarrow K$.

2. Locate the initial bracket at the bottom row, $j_l\leftarrow \mathrm{BinarySearch}(x^{\ast },{\{x_{j,k_l}\}}_j)$, and compute the intermediate coordinate ${\hat{y}}_l\leftarrow \mathrm{Interp1D}(x^{\ast };x_{j_l,k_l},x_{j_l+1,k_l},y_{j_l,k_l},y_{j_l+1,k_l})$.

3. While $k_h-k_l>1$:

   1. Set $k_m\leftarrow \lfloor (k_l+k_h)/2\rfloor$.

   2. Find $j_m\leftarrow \mathrm{BracketSearch}(x^{\ast },{\{x_{j,k_m}\}}_j,j_l)$, searching near $j_l$ within the configured bracket-displacement window.

   3. Compute ${\hat{y}}_m\leftarrow \mathrm{Interp1D}(x^{\ast };x_{j_m,k_m},x_{j_m+1,k_m},y_{j_m,k_m},y_{j_m+1,k_m})$.

   4. If the intermediate ordinate satisfies ${\hat{y}}_m\le y^{\ast }$, set $k_l\leftarrow k_m$, $j_l\leftarrow j_m$, and ${\hat{y}}_l\leftarrow {\hat{y}}_m$; otherwise set $k_h\leftarrow k_m$.

Pass 2: final interpolation.

1. Find $j_h\leftarrow \mathrm{BracketSearch}(x^{\ast },{\{x_{j,k_h}\}}_j,j_l)$.

2. Compute ${\hat{y}}_h\leftarrow \mathrm{Interp1D}(x^{\ast };x_{j_h,k_h},x_{j_h+1,k_h},y_{j_h,k_h},y_{j_h+1,k_h})$.

3. Compute ${\hat{f}}_l\leftarrow \mathrm{Interp1D}(x^{\ast };x_{j_l,k_l},x_{j_l+1,k_l},f_{j_l,k_l},f_{j_l+1,k_l})$.

4. Compute ${\hat{f}}_h\leftarrow \mathrm{Interp1D}(x^{\ast };x_{j_h,k_h},x_{j_h+1,k_h},f_{j_h,k_h},f_{j_h+1,k_h})$.

5. Set $t\leftarrow (y^{\ast }-{\hat{y}}_l)/({\hat{y}}_h-{\hat{y}}_l)$ and return $(1-t){\hat{f}}_l+t{\hat{f}}_h$.

EGM applications require an extrapolation rule: simulation occasionally produces state values slightly outside the solution grid because of shock realizations or numerical precision. Each 1D interpolation step handles out-of-bounds queries by extending the boundary slope, a continuous extension beyond the grid’s convex hull: for $x^{\ast }<x_1$, the value is $f_1+(x^{\ast }-x_1)(f_2-f_1)/(x_2-x_1)$, and symmetrically for $x^{\ast }>x_J$.

Linear extrapolation is unreliable near binding constraints, where policies have kinks and high curvature, and accurate in wealth-rich regions, where precautionary motives fade and policies become nearly linear. The asset grids in this paper extend well past their constraint regions, which confines any extrapolation to that smooth region.

7.3Notation Reference¶