Solving the Bewley-Huggett-Aiyagari Model
Heterogeneous Agent (HA) models feature idiosyncratic risk and incomplete markets.
Why is it hard? We must solve for both:
General Equilibrium: Prices (\(r\), \(w\)) depend on aggregate capital \(K\), which is the integral of individual asset holdings over the distribution \(\Gamma(a, z)\).
Solving an Aiyagari model in steady state typically involves these nested loops:
The individual state consists of:
To solve numerically, we restrict both to finite, discrete grids.
If income follows an AR(1) process: \[ \log z_{t} = \rho \log z_{t-1} + \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0, \sigma^2) \]
We approximate this with a discrete Markov Chain:
Standard Methods:
We define a grid of \(n_a\) points: \(A = \{a_1, a_2, \dots, a_{n_a}\}\) where \(a_1 = \underline{a}\).
\[ a_i = \underline{a} + (a_{max} - \underline{a}) \left( \frac{i-1}{n_a-1} \right)^\theta \] where \(\theta > 1\) (e.g., \(\theta = 2\) or \(3\)).
The standard Bellman equation approach: \[ V^{(k+1)}(a, z) = \max_{c, a'} \left\{ u(c) + \beta \sum_{z'} \Pi(z, z') V^{(k)}(a', z') \right\} \] Subject to the budget constraint.
Instead of the Value Function, iterate on the policy function using the Euler Equation.
Given a current guess of next period’s savings policy \(a'^{(k)}(a,z)\):
For each grid point \((a_i, z_j)\), define the residual function: \[R(x) = u'\!\left((1+r)a_i + w z_j - x\right) - \max\!\left\{ u'\!\left((1+r)a_i + wz_j - \underline{a}\right),\; \beta(1+r)\,\mathbb{E}_{z_j}\!\left[u'\!\left(c^{(k)}(x, z')\right)\right]\right\}\] where \(x = a'\) is the unknown savings choice.
Root-finding finds the \(x^* = a'^{(k+1)}(a_i, z_j)\) such that \(R(x^*) = 0\):
Repeat for all \((a_i, z_j)\) and iterate until \(\|a'^{(k+1)} - a'^{(k)}\|_{\infty} < \varepsilon\).
Faster than VFI, but requires a scalar root-finding call (e.g., bisection) at every grid point.
Carroll (2006): A massive speedup for solving HA models. Instead of fixing \(a\) and finding \(a'\) using root-finding, we heavily reduce the computational burden.
The trick: Fix tomorrow’s assets \(a'\), and solve for today’s endogenous assets \(a\).
Set up an exogenous grid for tomorrow’s assets: \(\vec{a}'\).
Use the Euler equation to directly compute consumption today corresponding to each \(a'\).
\[ c_i = \left( \beta (1+r) \mathbb{E}[u'(c(a'_i, z'))] \right)^{-1/\gamma} \]
Compute today’s endogenous assets using the budget constraint:
\[ a_i = \frac{c_i + a'_i - w \cdot z}{1+r} \]
Interpolate the policy function defined on the endogenous grid \((a_i, c_i)\) back onto the standard exogenous grid.
No root finding required!
Once the policy function \(a'(a, z)\) is converged, we need to know how many people are at each state. Let \(\Gamma(a, z)\) be the probability mass function.
We track the evolution of the distribution over the state space. \[ \Gamma_{t+1}(a', z') = \sum_{z} \Pi(z, z') \sum_{\{a | a'(a,z) = a'\}} \Gamma_t(a, z) \]
The target policy \(a'(a, z)\) typically does not land exactly on a grid point.
E.g., suppose optimal savings is \(a^*\), which falls between grid points \(a_k\) and \(a_{k+1}\).
Lottery Approach (Young 2010): Assign weights to the adjacent grid points.
The law of motion for the distribution can be written as a large linear system: \[ \vec{\Gamma}_{t+1} = M \vec{\Gamma}_t \] Where \(M\) is the transition matrix of size \((n_a \times n_z) \times (n_a \times n_z)\).
To find the steady-state distribution \(\Gamma^*\):
We have aggregated capital supply: \[ K^s(r) = \sum_{a, z} a \cdot \Gamma^*(a, z) \] (Assuming wage \(w\) is determined by \(r\) via the Firm’s FOC: \(w(r) = (1-\alpha) \left( \frac{\alpha}{r+\delta} \right)^{\frac{\alpha}{1-\alpha}}\)).
And we have capital demand from firms: \(K^d(r)\).
The Objective: Find \(r^*\) such that \(Excess Demand = K^d(r) - K^s(r) = 0\).
We treat \(ExcessDemand(r)\) as a black box function.
Algorithms:
scipy.optimize.brentq).Solving for transition dynamics involves these steps:
In any period \(t\), general equilibrium requires three objects to be mutually consistent:
| Object | Symbol | Determined by |
|---|---|---|
| Household policy | \(a'_t(a,z)\) | Euler equation given \((r_t, w_t)\) |
| Cross-sectional distribution | \(\Gamma_t(a,z)\) | Law of motion given \(a'_{t-1}\) |
| Prices | \((r_t, w_t)\) | Firm FOCs: \(r_t = F_K(K_t, Z_t)-\delta\), \(w_t = F_L(K_t, Z_t)\) |
Market clearing links them into a nonlinear system for the aggregate capital path \((K_t)_{t\geq 0}\): \[\underbrace{K_{t+1}}_{\text{supply}} = \mathcal{A}(K_t, Z_t, K_{t+1}, \ldots) \equiv \int a'_t(a,z; r_t, w_t)\, d\Gamma_t(a,z)\]
where \(\mathcal{A}\) encodes all the micro behavior. This is the equation whose solution is the equilibrium.
Let \(\bar{K}\) be the steady-state capital stock and \(\hat{K}_t = K_t - \bar{K}\) be the deviation. A first-order approximation yields:
\[\hat{K}_{t+1} = \phi_{KK}\, \hat{K}_t + \phi_{KZ}\, \hat{Z}_t\]
where the coefficients \(\phi_{KK}\), \(\phi_{KZ}\) are functions of two Jacobian matrices:
\[\mathbf{J}^{K} = \frac{\partial \mathcal{A}}{\partial r}\bigg|_{\bar{K}}, \qquad \mathbf{J}^{Z} = \frac{\partial \mathcal{A}}{\partial Z}\bigg|_{\bar{K}}\]
Each column of \(\mathbf{J}^K\) answers: “If \(r\) is 1% higher in period \(s\), how does aggregate capital supply in period \(t\) respond?”
For a finite horizon \(T\), stack all deviations into vectors \(\hat{\mathbf{K}} = (\hat{K}_0,\dots,\hat{K}_T)^\top\) and \(\hat{\mathbf{Z}} = (\hat{Z}_0,\dots,\hat{Z}_T)^\top\).
The full aggregate equilibrium system becomes:
\[\underbrace{(I - \mathbf{J}^K \mathbf{H}^{rK})}_{\text{$(T\!+\!1)\times(T\!+\!1)$ matrix}}\hat{\mathbf{K}} = \mathbf{J}^Z \hat{\mathbf{Z}}\]
where \(\mathbf{H}^{rK}\) maps capital to interest rates via firm FOCs. Key steps:
Libraries: Sequence Space Jacobians (Auclert et al.), Dolo.
Modern libraries (e.g., EconForge tools, Dolo, Sequence Space Jacobians) automate many of these steps, allowing us to focus on the economics!