1.12. Heterogeneous agents (HANK)
A heterogeneous-agent model (HANK, and more generally any
incomplete-markets / Aiyagari-type economy) is authored in the
ordinary .rs grammar with three additional constructs.
RISE reads the annotated file and expands it into a plain
Markov-switching DSGE; from that point the object behaves like any
other model – solve, filter, estimate, forecast,
simulate, irf and the decompositions all apply, and the
heterogeneity composes with regime switching.
A model file that uses none of these constructs is an ordinary representative-agent model and is parsed unchanged.
The three constructs are:
@heterogeneity_axis(...)– declare a dimension of the cross-section (a grid).@endogenous(individual)– declare the variables solved at the individual (per-grid-point) level.@Agg(...)– aggregate an individual expression across the cross-section.
The expanded model repeats the same few equations once per grid cell; template differentiation (on by default) differentiates each of those equations once instead of once per cell; on a two-asset model with 150 cells this cuts differentiation from about 11 minutes to 40 seconds and peak memory from 30 GB to 5 GB, with identical derivatives. See Template differentiation.
This chapter documents what each construct means and what it does. It is a parser-level interface: the model you write is the model that is solved. The solution itself is handled by RISE’s standard perturbation machinery, with the cross-section treated as an extension of Reiter’s method; that part is internal and not the subject of this chapter.
1.12.1. Heterogeneity axes
Each dimension of the cross-section is declared with one line:
@heterogeneity_axis(name, size, <option>)
name is the axis label, size is the number of grid points,
a positive integer: a literal (10), the name of a rise_flags
entry (@heterogeneity_axis(a, na, policy = ap) parsed with
dsge_model('m.rs', 'rise_flags', struct('na', 10))), or a
preparser expression such as @{2*na}. A name that is not among the
rise_flags is an error (hank_expand:axisSize). The same holds
for the value of span (see the bracket window section below). The third field
selects the kind of axis, in one of two ways: policy = <var>
makes it an endogenous axis, and evolution = <kind> selects
any other kind, where <kind> is one of transition (the
default), fixed or deterministic:
Declaration |
Meaning |
|---|---|
|
Endogenous (policy) axis. The grid (here |
|
Exogenous Markov axis. |
|
Fixed-type axis. A permanent characteristic that never changes (a household type); mass is conserved within each type. |
|
Deterministic axis. Each agent moves one grid point
forward per period; the |
For every axis, RISE automatically creates the grid-value
parameters name_1, name_2, ..., name_size. In addition:
an exogenous Markov axis creates the off-diagonal transition probabilities
name_htp_i_j(forinot equal toj); each diagonal probability is the residual1 - sum of its row's off-diagonals;a fixed-type axis creates the type masses
name_mass_i;a deterministic axis creates no transition parameters – its transition matrix is known (see Deterministic axes).
You only declare the axes; you supply their numeric content the same way as any other parameter – see Setting grid and transition values below.
Inside an equation, the axis name used bare (a, e) is
the coordinate of the current grid point – a fixed label, equal
to the grid value at that point. A coordinate is not a variable: it
carries no {t}, {t+1} or {t-1} and has no law of
motion.
1.12.2. Deterministic axes
A deterministic axis moves each agent one grid point forward per
period – s to s+1 – with no randomness. Its transition
matrix is therefore known, so RISE creates no name_htp_i_j
parameters (only the grid values name_i). This is the usual
representation of a lifecycle / age dimension.
The only freedom is what happens at the top of the grid, set by
boundary (default cyclic). For a four-point axis the
transition matrix \(T\) (row \(i\) = current point, column
\(j\) = next point; \(T_{ij}=1\) when an agent at \(i\)
moves to \(j\)) is, under each rule,
Every row has a single 1 – each agent has exactly one
destination – so mass is conserved under all three rules. They
differ only in the last row:
cyclicsends the top point back to the first (the matrix is a permutation): “die at the top, reborn at the bottom”. It is the default and the only rule that is non-degenerate on its own.absorbingholds agents at the top point (a terminal / catch-all bin); mass drifts up and the lower bins drain unless they are refilled by another (stochastic) axis.reflectingsends the top point back one step.
Because the destination is known, two things simplify relative to a
Markov axis: the distribution law of motion is sparser – each
cell receives from a single predecessor along the deterministic axis
– and an individual lead c{t+1} across a deterministic axis is
simply the value at the next point (there is no expectation to
take):
@heterogeneity_axis(a, 10, policy = ap) % assets
@heterogeneity_axis(s, 6, evolution = deterministic) % age, cyclic (default)
@heterogeneity_axis(s, 6, evolution = deterministic, boundary = absorbing)
1.12.3. Individual variables
The variables solved at each grid point are declared with the
individual qualifier:
@endogenous(individual)
c "consumption", ap "savings policy", mu "borrowing multiplier"
You declare each individual variable once. RISE replicates it
across the whole cross-section – one copy per grid cell. Optional
quoted descriptions are carried through to the per-cell variables,
as in the ordinary @endogenous block.
Aggregate variables, shocks and parameters are declared exactly as
usual with @endogenous, @exogenous and @parameters.
1.12.4. The individual block
You write each individual equilibrium condition once, in the
@model block, including its first-order conditions – RISE does
not take FOCs or reformulate anything. The single equation is
copied to every grid cell.
Two rules govern how symbols in an individual equation are read:
Individual variables carry time indices (
c{t},ap{t},c{t+1}); the bare axis coordinates do not (a,e).A lead on an individual variable –
c{t+1}– is read as the idiosyncratic conditional expectation: the expectation over the agent’s own movement across the cross-section. A lead on an aggregate variable –r{t+1}– keeps the ordinary \(E_t\). You simply write the lead; RISE forms the appropriate expectation.
For example, a consumption Euler equation:
c{t}^(-sigma) = beta*(1+r{t+1})*c{t+1}^(-sigma) + mu{t};
Here c{t+1} is individual (idiosyncratic expectation over the
agent’s next position), r{t+1} is aggregate (ordinary
\(E_t\)), and mu is the individual multiplier on the
borrowing limit.
A budget constraint using the bare coordinates:
c{t} = (1+r{t})*a + w{t}*e - ap{t};
a is the agent’s current asset position and e the current
income state – both fixed grid labels at this cell.
Occasionally-binding individual constraints (a borrowing limit) are written as your own kinky equation, exactly as the zero-lower-bound is written in the representative-agent case (see Occasionally-binding constraints):
max( -mu{t}, amin - ap{t} ) = 0;
The complementarity form is passed through untouched.
1.12.5. Aggregation: @Agg
@Agg(expr) is the cross-sectional aggregate of an
individual expression – its integral over the population
distribution. It is the only place the distribution enters the
model you write:
A{t} = @Agg(ap{t}); % aggregate savings = integral of ap
C{t} = @Agg(c{t}); % aggregate consumption
NE{t} = @Agg(e*n{t}); % effective labour supply
K{t} = @Agg(a); % aggregate assets = integral of the coordinate a
The argument is any individual expression – an individual variable,
a bare coordinate, or a product of them. @Agg is typically used
to form the aggregates that enter market-clearing conditions
(A{t} = B, K{t} = ...) and the aggregate block.
1.12.6. The aggregate block
Everything else is the ordinary representative-agent (“RANK”)
shell, written exactly as in a standard .rs file: firm
conditions, the monetary and fiscal rules, the New Keynesian
Phillips curve, market clearing, and the exogenous driving
processes. These equations reference the aggregates produced by
@Agg and the usual aggregate variables, with no special
syntax:
r{t} = alpha*exp(z{t})*(K{t}/L)^(alpha-1) - delta;
w{t} = (1-alpha)*exp(z{t})*(K{t}/L)^alpha;
A{t} = B; % bond-market clearing
z{t} = rho_z*z{t-1} + eps_z{t};
1.12.7. Setting grid and transition values
The grid values name_i, the Markov transitions name_htp_i_j
and the type masses name_mass_i are ordinary parameters. Set
them like any other parameter – in @parameterization, through
the macro rise_flags, or from the driver script with
set(m, 'parameters', ...). A typical driver computes a
discretised income process (for example a Rouwenhorst grid) for the
Markov axis and a chosen asset grid for the policy axis, then
assigns the resulting name_i and name_htp_i_j values:
m = dsge_model('my_hank.rs');
m = set(m, 'parameters', struct( ...
'a_1', 0, 'a_2', 0.26, ... % asset grid
'e_1', 0.54, 'e_2', 1.46, ... % income grid
'e_htp_1_2', 0.017, 'e_htp_2_1', 0.017)); % income transitions
The steady state is then solved (or imposed) as for any model, and
the object is ready for solve, irf, filter and the
rest.
1.12.8. Bin indicators: addressing a cell by position
The bare axis coordinate expands to the grid value at each cell, so two
bins that share a value cannot be told apart by it. The bin indicator
<axis>_at_<j> resolves at expansion time to the literal 1 at index
j on that axis and 0 elsewhere — it addresses a bin by position.
The canonical use is per-type income (Bayer–Lütticke’s entrepreneur, income
state 4, shares a worker’s productivity value):
c{t} + ap{t} = RB{t}*a + tauP*NW{t}*h*(1-h_at_4)
+ tauP*profitshare*Profits{t}*h*h_at_4;
Indicators cost no parameter and no equation (they are resolved to constants), and an out-of-range index is a parse error, never a silent zero.
1.12.9. The bracket window: span
Note
Most models should not use span. It is a declared reduction, and
its validity depends on the parameter values — which is what the rest of
this section is about guarding. The automatic, parameter-free
alternative is solve_prune_hank_constants (below), which discovers
the same constants and more without an assertion. Reach for span
only when the grid is so large that building the full symbolic model is
itself prohibitive.
By default every cell carries size-1 interpolation weights theta on
each endogenous axis, so the expansion grows quadratically in the grid.
The optional span narrows the live weights to the s brackets on
either side of the cell’s own index:
@heterogeneity_axis(a, 100, policy = ap, span = 3)
asserting that the policy from any cell lands within s bins of it.
Outside the window the clamp is a constant, carried exactly, so a covering
window (span = size-1) reproduces the default expansion byte-for-byte.
Measured growth: 7,692 endogenous at 60 asset points without a window,
1,420 with span = 3 — O(N²) to O(N).
The assertion is parameterization-dependent — move the parameters and a
window that was valid can stop being valid — so it is machine-guarded
rather than trusted. The het engine (next section) always solves on the
full grid and measures the landing distance d(i,j) — how many brackets
each cell’s policy moves it — at every solve, returning out.landing
(dmax, min_safe_span = dmax+1, and the per-cell map). Choose s
from min_safe_span; if a declared span is violated at the current
parameter values, the engine stops with a hard error
(het:engine:spanViolated) naming the minimal safe declaration. A
violated window would otherwise be dangerous precisely because the pruned
model stays well-formed: a solver can converge cleanly to a wrong steady
state. span must be a positive integer; span = 0 is rejected at
parse time (the top cell’s window would be empty).
1.12.10. Solving the steady state: the het engine
A heterogeneous-agent steady state is a nested problem — a household policy that is a contraction given prices, an invariant distribution that is linear given the policy, and an aggregate block closing prices — and a flat Newton over the expanded system does not solve it from any practical starting point. The in-core engine solves it in that nested form:
spec = struct();
spec.free = {'beta','vphi'}; % free parameters...
spec.targets = {'A', B ; 'NE', 1}; % ...against these aggregates
spec.prices = @(agg,P) my_closure(agg,P); % the aggregate-block closure
[ss, info] = rise.engine.dsge_tools.sstate.het.steady_state(m, spec);
m = solve(set(m,'sstate_imposed',true), 'sstate_bounds', ss);
The engine reads everything else from the parsed model: the individual
block is solved by Anderson-accelerated time iteration (all cells at
once, Fischer–Burmeister at each borrowing constraint), the cross-section by
the Tan (2020) factored forward iteration — it applies the endogenous
lottery and each exogenous chain as separate mode products, never forming the
fused transition, so it stays tractable with several exogenous axes (the
Young (2010) fused iteration remains available as a single-axis
reference). Aggregates come through the model’s @Agg definitions. The
per-cell Newton builds its Jacobian by complex-step differentiation —
exact, no step-size error — whenever the individual equations are holomorphic,
and by finite differences otherwise; the choice is automatic, a non-converging
complex-step solve retries once with finite differences (from where complex
step stopped), and out.jac reports which was used (opts.jac
overrides). The acceleration is guarded: where cells sit at a borrowing
corner the history of the mixing can turn nearly collinear, so a mixed step
much longer than the plain one is not taken, and a sweep whose step blows up
restarts from the best iterate so far. After repeated blow-ups the solve stops
with out.retcode 1 and the warning het:engine:diverged instead of
sweeping on; another start (opts.warm), damping (opts.damp) or grid
is then the remedy. The prices closure
encodes how prices respond to the aggregates: constant closures (normalized
NK) converge in one pass; feedback closures (Krusell–Smith, r/w from
K) iterate with damping (spec.price_damp). spec.free = {} handles
pure fixed-point models; targets may be function handles for ratios
(@(agg) agg.A/agg.C).
Validated against four external references — Auclert et al. one-asset (ten digits), McKay–Nakamura–Steinsson (1e-12), Krusell–Smith (5.6e-12), and Bayer–Lütticke (5,625 equations, certified 1e-12, ~30x the hand-rolled reference). The engine handles any number of endogenous (policy), transition (Markov), fixed (permanent-type) and deterministic (life-cycle) axes, with any number of lower-bound complementarity clauses, on a square individual block; a non-square block or a malformed axis declaration is the only hard error. Each axis kind is regression-tested against an independent reference: a deterministic axis against the equivalent permutation-matrix transition axis, a two-type fixed axis against mass-weighted per-type solves, a multi-shock model against its single joint-chain equivalent, and a dominated second asset against the one-asset solution.
1.12.11. Pruning the constants: solve_prune_hank_constants
The expansion creates one interpolation weight per bracket per cell because
parse time knows no parameters: it cannot know which weights the solved
policy will leave stuck at 0 or 1. At the steady state most of them
are stuck. A weight sitting on a rail is a constant, its deviation is
identically zero, and carrying it into the solution step is pure cost.
solve_prune_hank_constants removes them, automatically:
m = solve(m, 'sstate_bounds', ss, 'solve_prune_hank_constants', true);
There is nothing to choose and nothing to declare. The solver reads the
evaluated Jacobian and finds the variables whose equations pin their
first-order deviations at exactly zero — a saturated clamp differentiates
to a row whose only entry is its own diagonal, forcing dtheta = 0.
Those variables and their equations are removed, the smaller system goes
through the unchanged first-order machinery, and the solution is re-inflated
with exact zeros. The answer is identical with the option on or off;
only the size of the system handed to the solver changes.
The reduction is on the system, and happens before whichever solver runs — a direct decomposition for a constant-parameter model, the iterative MSRE solvers under regime switching. It is therefore tied to no particular solution method, which matters because the switching case is where it earns the most: see the mean-square-stability note below.
One qualification on “identical”. Where the solution comes from a direct
decomposition — the single-regime QZ — the two answers agree to machine
precision. Where it comes from an iterative scheme, as it does for
regime-switching models solved by the MSRE solvers, pruning changes the
size of the problem and therefore the iteration path, so the two runs stop
at two points that are each within the solver’s convergence tolerance of
the same fixed point. The residual difference is that tolerance, not an
error introduced by pruning: tighten solver{2}.TolFun and it follows
the tolerance down (measured on a two-regime loose-commitment model:
5.3e-8 at the default sqrt(eps), 2.5e-14 at 1e-14).
What it removes is wider than the weights: on the Auclert one-asset model it
takes out 182 of 292 variables — 161 saturated theta, 19 slack borrowing
multipliers (lambda == 0 on unconstrained cells), a policy pinned at the
borrowing limit, and a fixed-supply aggregate. No declared window can
express that set, because it is a property of the solution, not of the
model file.
Detection is a fixed point, not a single pass: killing one variable’s
columns can make another row trivial, so the scan repeats until nothing
more falls. A row counts as trivial only if it has a single nonzero, in the
same current-date column in every regime pair, with no shock or sigma
entry — that last condition is what protects a driven process whose
persistence happens to be zero, which otherwise has exactly the shape of a
pinned constant. Under regime switching a variable must additionally hold
the same steady-state level across regimes: one saturated at 1 in regime
A and 0 in regime B has zero deviations but a level jump, which the
misalignment machinery must keep carrying. If the reduction would not be
square the solver refuses, warns, and does the full solve.
Two consequences worth knowing. Pruned states’ columns in the transition matrix are set to zero by convention: they are unreachable on-model, since those variables never deviate, whereas an unpruned solve leaves small round-off there. And the option makes the mean-square-stability check feasible for switching heterogeneous models — that check builds a Kronecker product quadratic in the state count, and on a two-regime HANK it otherwise runs out of memory. That is the clearest illustration of the point above: nothing about it involves a constant-parameter decomposition, and it turns a check that could not run at all into one that finishes in seconds.
The option is binary, parameter-free and defaults to false. It applies
at first order; higher orders continue at full dimension. It is named for
where it is validated — the HANK replication suite — though the mechanism
itself is generic.
Note
Prefer this to span. The declared window is the fallback for the one
case pruning cannot reach: a grid so large that building the full
symbolic model is itself prohibitive. Pruning needs the model to exist
before it can read its Jacobian; span acts earlier, at construction.
This is distinct from the other two things called “pruning” in RISE:
simul_pruned discards explosive higher-order cross-terms when
simulating an already-solved model, and distribution/state reduction
(Reiter, Bayer–Lütticke) replaces the histogram by a low-dimensional
parameterization. Those change the answer by approximation; this one does
not change the answer at all.
The state reduction is available too, as a separate option:
solve_reduction writes the individual variables as their steady
state plus a basis times a few coordinates and solves and filters the small
system, at any order and with regime switching and occasionally-binding
constraints. See Model reduction.
1.12.12. Composing with regime switching
Heterogeneity and Markov switching are orthogonal. A switching
parameter block – @parameters(chain, N) name with its implicit
chain_tp_i_j transitions (see Model file language) –
is carried through the expansion unchanged, so the parameters of a
HANK model can switch across regimes exactly as in a
representative-agent model.
1.12.13. Worked example
A one-asset HANK (the replication target is the canonical sequence-jacobian one-asset example): an endogenous bond axis, an exogenous Rouwenhorst skill axis, individual consumption / savings / hours / borrowing multiplier, and a New Keynesian aggregate shell:
% ---- cross-section axes ----
@heterogeneity_axis(a, 10, policy = ap) % endogenous bond axis
@heterogeneity_axis(e, 2) % exogenous skill axis -> e_htp_i_j
% ---- individual variables ----
@endogenous(individual)
c "consumption", ap "bond policy", n "hours", lam "borrowing multiplier"
% ---- aggregate variables, shocks, parameters ----
@endogenous Y w pai r rstar Z Div Tax A NE C L
@exogenous eps_r eps_z
@parameters beta eis frisch vphi mkp kappa phi B rho_r rho_z amin
@model
% individual block (written ONCE)
"Euler" c{t}^(-1/eis) = beta*(1+r{t+1})*c{t+1}^(-1/eis) + lam{t};
"Labor" vphi*n{t}^(1/frisch) = w{t}*e*c{t}^(-1/eis);
"Budget" c{t} + ap{t} = (1+r{t})*a + w{t}*e*n{t} + (Div{t}-Tax{t})*e;
max( -lam{t}, amin - ap{t} ) = 0;
% aggregation (the only place the distribution appears)
A{t} = @Agg(ap{t});
C{t} = @Agg(c{t});
NE{t} = @Agg(e*n{t});
% aggregate New Keynesian block (ordinary RANK equations)
L{t} = Y{t}/Z{t};
r{t} = (1 + rstar{t-1} + phi*pai{t-1})/(1 + pai{t}) - 1;
Tax{t} = r{t}*B;
A{t} = B; % bond-market clearing
NE{t} = L{t}; % labour-market clearing
rstar{t} = rho_r*rstar{t-1} + (1-rho_r)*rstar{stst} + eps_r{t};
Z{t} = rho_z*Z{t-1} + (1-rho_z)*1 + eps_z{t};
@parameterization
eis = 0.5; frisch = 0.5; mkp = 1.2; kappa = 0.1; phi = 1.5;
B = 5.6; amin = 0.0; rho_r = 0.0; rho_z = 0.0;
% grid + transition values (set here or overwritten by the driver)
a_1 = 0; a_2 = 0.26; % ... a_10 = 150;
e_1 = 0.537870; e_2 = 1.462130;
e_htp_1_2 = 0.017; e_htp_2_1 = 0.017;
The individual block is written once; c{t+1} is the
idiosyncratic expectation; a and e are the bin coordinates;
@Agg forms the aggregates that enter market clearing. Once
parsed, the object is an ordinary RISE model.