Achdou et al., Online Appendix: Numerical Methods for 'Income and Wealth Distribution in Macroeconomics'
Summary
An online appendix by Yves Achdou, Jiequn Han, Jean-Michel Lasry, Pierre-Louis-Lions and Benjamin Moll on numerical methods for “Income and Wealth Distribution in Macroeconomics: A Continuous-Time Approach,” cited as Achdou et al. (2020). It explains how to solve the continuous time heterogeneous agent models of that paper with a finite difference method, starting with a continuous time version of Huggett (1993). The appendix sets out explicit and implicit methods for the Hamilton-Jacobi-Bellman (HJB) equation and states that the implicit method is the preferred approach because it is both more efficient and more stable. Its closing pages show results for a two-asset model with a Pareto tail of the wealth distribution and list references.
Summary drafted by a model from the document's text below and checked by script against that text before publication. It is a navigation aid, not a reading of what the document proves. Where AI is used
- Issuer
- Source documents
- Document type
- Countercyclical Unemployment Benefits A General Equilibrium Analysis Of Transition Dynamic
- Case
- Countercyclical Unemployment Benefits A General Equilibrium Analysis Of Transition Dynamics
Cite this document
“Achdou et al., Online Appendix: Numerical Methods for 'Income and Wealth Distribution in Macroeconomics'.”
Original at doi.org. Archived at Pandemic Darlings, https://pandemicdarlings.com/source-documents/countercyclical-unemployment-benefits-a-general-equilibrium-analysis-of-transition-dynamics/
, SHA-256 214d6347a8f8ff1883dc722358038c5f4dc2453089cc0087fa9b46987e2f5d08.Full text
Online Appendix: Numerical Methods for
“Income and Wealth Distribution in Macroeconomics:
A Continuous-Time Approach”1
Yves Achdou, Jiequn Han, Jean-Michel Lasry, Pierre-Louis-Lions, Benjamin Moll
This is an online Appendix to Achdou et al. (2020). It is concerned with the numerical solution
using a finite difference method of the continuous time heterogeneous agent models presented
in that paper. Also see the discussion in Section 4 of the paper, particulary the material on
the conditions of Barles and Souganidis (1991).
Huggett Economy
We start by solving a continuous time version of Huggett (1993) which is arguably the simplest
heterogeneous agent model that captures many of the features of richer models. The economy
can be represented by the following system of equations which we aim to solve numerically:
ρv1 (a) = max u(c) + v10 (a)(z1 + ra − c) + λ1 (v2 (a) − v1 (a)) (1)
c
ρv2 (a) = max u(c) + v20 (a)(z2 + ra − c) + λ2 (v1 (a) − v2 (a)) (2)
c
d
0=− [s1 (a)g1 (a)] − λ1 g1 (a) + λ2 g2 (a) (3)
da
d
0 = − [s2 (a)g2 (a)] − λ2 g2 (a) + λ1 g1 (a) (4)
Z da ∞ Z ∞
1= g1 (a)da + g2 (a)da (5)
a a
Z ∞ Z ∞
0= ag1 (a)da + ag2 (a)da ≡ S(r) (6)
a a
where z1 < z2 and sj (a) = zj + ra − cj (a) and cj (a) = (u0 )−1 (vj (a)) are optimal savings and
consumption. Finally, there is a state constraint a ≥ a. The first order condition u0 (cj (a)) =
vj0 (a) still holds at the borrowing constraint. However, in order to respect the constraint we
need sj (a) = zj + ra − cj (a) ≥ 0. Combining this with the FOC, the state constraint motivates
a boundary condition
vj0 (a) ≥ u0 (zj + ra), j = 1, 2 (7)
1
We thank SeHyoun Ahn for fantastic research assistance. Matthieu Gomez provided extremely useful
comments and in particular showed us how to clearly think about non-uniform grids in Appendix Section 7.
1
We use a finite difference method. A useful reference is Candler (1999). We first explain
how to solve the Hamilton-Jacobi-Bellman (HJB) equation (1) and (2), and then turn to the
Kolmogorov Forward (Fokker-Planck) equation (3) and (4).
Section 4 explains how the setup and the solution method can be generalized to an en-
vironment where productivity z is continuous and follows a diffusion rather than a two-state
Poisson process. Finally, a useful “warm-up problem” is to solve the HJB equation with no un-
certainty, λj = 0. See http://www.princeton.edu/~moll/HACTproject/HACT_Additional_
Codes.pdf. All algorithms are available as Matlab codes from https://benjaminmoll.com/
codes/. We are especially indebted to SeHyoun Ahn for showing us how to use matlab’s sparse
matrix routines to increase speed by an order of magnitude.
1 HJB Equation
We use a finite difference method and approximate the functions (v1 , v2 ) at I discrete points
in the space dimension, ai , i = 1, ..., I. We use equispaced grids, denote by ∆a the distance
between grid points, and use the short-hand notation vi,j ≡ vj (ai ) and so on. The derivative
0
vi,j = vj0 (ai ) is approximated with either a forward or a backward difference approximation
vi+1,j − vi,j 0
vj0 (ai ) ≈ ≡ vi,j,F
∆a (8)
vi,j − vi−1,j
vj0 (ai ) ≈ 0
≡ vi,j,B
∆a
The finite difference approximation to (1) and (2) is
0
ρvi,j = u(ci,j ) + vi,j (zj + rai − ci,j ) + λj (vi,−j − vi,j ), j = 1, 2
0 −1 0
(9)
ci,j = (u ) (vi,j )
0
where vi,j is either the forward or the backward difference approximation. There are two
complications. The first question is when to use a forward and when a backward difference
approximation. It turns out that this is actually quite important for the stability properties
of the scheme. The second is that the HJB equations are highly non-linear, and therefore so
is the system of equations (9). It therefore has to be solved using an iterative scheme (rather
than simply inverting a matrix).
There are two options that differ in how the value function is updated: a so-called “explicit”
method and an “implicit” method. As a general rule, the implicit method is the preferred
approach because it is both more efficient and more stable/reliable. However, the explicit
method is easier to explain so we turn to it first.
2
1.1 Explicit Method
See matlab program HJB_stateconstraint_explicit.m. One starts with an initial guess
0
0
vj0 = (v1,j , ..., vI,j ), j = 1, 2 and then updates vjn , n = 1, ... according to
n+1 n
vi,j − vi,j n 0
n
+ ρvi,j = u(cni,j ) + (vi,j ) (zj + rai − cni,j ) + λj (vi,−j
n n
− vi,j ) (10)
∆
where cni,j = (u0 )−1 [(vi,j
n 0
) ]. The parameter ∆ is the step size of the explicit method. It can
be shown that the explicit method only converges if ∆ is not too large (it has to satisfy the
so-called “CFL condition”, see e.g. p.181 in Candler (1999)). An advantage of implicit methods
discussed in section 1.2 is that the step size ∆ can be arbitrarily large.
Upwind Scheme. As already mentioned, it is important whether and when a forward or a
backward difference approximation is used. The correct way of doing this is to use a so-called
“upwind scheme.” The rough idea is to use a forward difference approximation whenever the
drift of the state variable (here, savings sni,j = zj + rai − cni,j ) is positive and to use a backwards
difference whenever it is negative.2 In practice, this is done as follows: first compute savings
0 0
according to both the backwards and forward difference approximations vi,j,F and vi,j,B
si,j,F = zj + rai − (u0 )−1 (vi,j,F
0
), si,j,B = zj + rai − (u0 )−1 (vi,j,B
0
)
where we suppress n superscripts for notational simplicity. Then use the following approxima-
0
tion for vi,j :
0 0 0 0
vi,j = vi,j,F 1{si,j,F >0} + vi,j,B 1{si,j,B <0} + v̄i,j 1{si,j,F ≤0≤si,j,B } (11)
where 1{·} denotes the indicator function. The meaning of the last term is as follows. First
0 0
note that since v is concave in a, we have vi,j,F < vi,j,B and so si,j,F < si,j,B . Therefore, for some
grid points i, si,j,F ≤ 0 ≤ si,j,B . At these grid points, we set savings equal to zero and hence
0
set the derivative of the value function equal to v̄i,j = u0 (zj + rai ). The fact that v is concave
also means that we do not have to worry about the case where both si,j,F > 0 and si,j,B < 0:
because concavity implies si,j,F < si,j,B , this cannot happen.
Under some circumstances (and in more general applications e.g. problems with non-
convexities) it can happen that the value function is not concave. The question then arises
what to do when both si,j,F > 0 and si,j,B < 0. The following upwind scheme works well in
practice in this case. Define an indicator for the problematic case in which both si,j,F > 0
2
Note that treatments of finite difference methods concerned with solving PDEs forward in time usually define
an “upwind scheme” in the opposite way (forward difference when the drift is negative, backwards difference
whenever it is positive). See e.g. https://en.wikipedia.org/wiki/Upwind_scheme. The difference is that
solving HJB equations (even stationary ones) amounts to solving PDEs backwards in time given a terminal
condition (rather than forward in time given an initial condition). The two seemingly different definitions of
the term “upwind scheme” are the same when taking this difference into account.
3
unique
and si,j,B < 0, 1both i,j := 1{si,j,B ≤0≤si,j,F } , and an indicator for the unproblematic case 1i,j =
1{si,j,F <0 and si,j,B >0} + 1{si,j,F <0 and si,j,B <0} . Next, define the forward and backward Hamilto-
0
nians Hi,j,F := u(ci,j,F ) + vi,j,F si,j,F and similarly for Hi,j,B . Finally, use the upwind scheme
0
vi,j 0
= vi,j,F (1{si,j,F >0} 1unique
i,j + 1{Hi,j,F ≥Hi,j,B } 1both
i,j )
0
+ vi,j,B (1{si,j,B <0} 1unique
i,j + 1{Hi,j,F <Hi,j,B } 1both
i,j ) (12)
0
+ v̄i,j 1{si,j,F ≤0≤si,j,B }
Intuitively, in the problematic case when both si,j,F > 0 and si,j,B < 0, this upwind scheme uses
as the “tie breaker” the rule to use the derivative in the direction in which the gain according
to the Hamiltonians Hi,j,B and Hi,j,F is larger.
State Constraint. The state constraint (7) is enforced by setting
0
v1,j,B = u0 (zj + ra1 ), j = 1, 2
From (11) it can then be seen that the state constraint is imposed whenever the forward
difference approximation would result in negative savings s1,j,F ≤ 0. Otherwise if s1,j,F > 0
0
the forward difference approximation v1,j,F is used at the boundary, implying that the value
function “never sees the state constraint.” At the upper end of the state space, the upwind
method should make sure that a backward-difference approximation is used. In practice, it can
sometimes help stability of the algorithm to simply impose a state constraint a ≤ amax where
amax is the upper end of the bounded state space used for computations (this can be achieved
0
by setting vI,j,F = u0 (zj + raI )).
Initial Guess. A natural initial guess is the value function of “staying put”
0 u(zj + rai )
vi,j = .
ρ
Summary of Algorithm. Summarizing, the algorithm for finding a solution to the HJB
0
equation (1) and (2) is as follows. Guess vi,j , i = 1, ..., I, j = 1, 2 and for n = 0, 1, 2, ... follow
n 0
1. Compute (vi,j ) using (8) and (11).
2. Compute cn from cni,j = (u0 )−1 [(vi,j
n 0
)]
3. Find v n+1 from (10).
4. If v n+1 is close enough to v n : stop. Otherwise, go to step 1.
4
One can show that, for a small enough ∆, this algorithm satisfies the three conditions of Barles
and Souganidis (1991) (monotonicity, consistency, stability). See the discussion in Section 4 of
Achdou et al. (2020).
1.2 Implicit Method
See matlab program HJB_stateconstraint_implicit.m. Relative to the explicit scheme in
(10), an implicit differs in how v n is updated. In particular, v n+1 is now implicitly defined by
the equation
n+1 n
vi,j − vi,j n+1 n+1 0 n+1 n+1
+ ρvi,j = u(cni,j ) + (vi,j ) (zj + rai − cni,j ) + λj (vi,−j − vi,j )
∆
Note the n + 1 superscripts on the right-hand side of the equation.3 The main advantage of
the implicit scheme is that the step size ∆ can be arbitrarily large.
Upwind Scheme. As was the case for the explicit method, we need to use an “upwind
scheme.” As above, the idea is still to use the forward difference approximation whenever the
drift of the state variable is positive and the backward difference approximation whenever it is
negative. We use the following finite difference approximation to (1) and (2).
n+1 n
vi,j − vi,j n+1 0 n+1 0
n+1
+ ρvi,j =u(cni,j ) + (vi,j,F ) [zj + rai − cni,j,F ]+ + (vi,j,B ) [zj + rai − cni,j,B ]−
∆ (13)
n+1 n+1
+ λj [vi,−j − vi,j ]
n 0 n 0
where cni,j = (u0 )−1 [(vi,j ) ] and (vi,j ) is given by (11).4 For any number x, the notation x+
means “the positive part of x”, i.e. x+ = max{x, 0} and analogously x− = min{x, 0}, i.e.
[zj + rai − cni,j,F ]+ = max{zj + rai − cni,j,F , 0} and [zj + rai − cni,j,B ]− = min{zj + rai − cni,j,B , 0}.
Equation (13) constitutes a system of 2 × I linear equations, and it can be written in
matrix notation using the following steps. Substituting the definition of the derivatives (8),
and defining sni,j,F = zj + rai − cni,j,F and similarly for sni,j,B , (13) is
n+1 n n+1 n+1 n+1 n+1
vi,j − vi,j vi+1,j − vi,j vi,j − vi−1,j
n+1
+ ρvi,j = u(cni,j ) + (sni,j,F )+ + (sni,j,B )− + λj [vi,−j
n+1 n+1
− vi,j ]
∆ ∆a ∆a
3
Strictly speaking, the present method is a “semi-implicit method.” A fully implicit method would feature
n + 1 superscripts also on ci,j . Such a fully implicit scheme can be solved using a Newton method, which ends
up looking very similar to the iterative scheme outlined here.
4
As noted in the discussion of the explicit scheme in Section 1.1, this works well when the value function is
concave. When the value function is not concave, we can again use the scheme (12) instead of (11).
5
Collecting terms with the same subscripts on the right-hand side:
n+1 n
vi,j − vi,j n+1 n+1 n+1 n+1 n+1
+ ρvi,j = u(cni,j ) + vi−1,j xi,j + vi,j yi,j + vi+1,j zi,j + vi,−j λj where
∆
(sni,j,B )−
xi,j = − ,
∆a (14)
(sni,j,F )+ (sni,j,B )−
yi,j = − + − λj ,
∆a ∆a
(sni,j,F )+
zi,j =
∆a
n+1 n+1
Note that importantly x1,j = zI,j = 0, j = 1, 2 so v0,j and vI+1,j are never used. Equation (14)
is a system of 2 × I linear equations which can be written in matrix notation as:
1 n+1
(v − v n ) + ρv n+1 = un + An v n+1 (15)
∆
where
y1,1 z1,1 0 · · · 0 λ1 0 0 ··· 0
x
2,1 y2,1 z2,1 0 · · · 0 λ1 0 0 ··· u(cn1,1 )
..
0 x3,1 y3,1 z3,1 0 · · · 0 λ1 0 0
.
. ... ... ... ... ... ... ... ... .. ..
.. .
.
0 ... ... n
u(cI,1 )
xI,1 yI,1 0 0 0 0 λ1
An = , un =
u(cn )
λ2 0 0 0 0 y1,2 z1,2 0 0 0
1,2
.
..
0 λ2 0 0 0 x2,2 y2,2 z2,2 0 0
0 0 λ2 0 0 0 x3,2 y3,2 z3,2 0 ..
.
.. .. .. .. .. .. .. ..
0 0 . . . . . . . . u(cnI,2 )
0 ··· ··· 0 λ2 0 ··· 0 xI,2 yI,2
This system can in turn be written as
n n+1 n n 1 1 n
B v =b , B = + ρ I − An , bn = un + v (16)
∆ ∆
Equation (16) can be solved very efficiently in matlab using sparse matrix routines. To check
that one has constructed the intensity matrix A correctly, the matlab function spy is a conve-
nient tool. Figure 1 plots an example of spy’s output with 30 × 2 grid points (we usually use
many more).
Finally, it is instructive to consider the case with an infinite updating step size ∆1 = 0 and
to write the linear system (15) as
ρv n+1 = un + An v n+1 (17)
6
0
10
20
30
40
50
60
0 10 20 30 40 50 60
nz = 177
Figure 1: Visualization of the matrix A using Matlab’s spy function
It can be seen that (17) is just a way of writing the discretized version of the HJB equations
(1) and (2) in matrix form. In particular the matrix An encodes the evolution of the stochastic
process (at , zt ). The finite difference method basically approximates this process with a discrete
Poisson process with a transition matrix An summarizing the corresponding Poisson intensi-
ties. Note that An satisfies all the properties a Poisson transition matrix needs to satisfy. In
particular, all rows sum to zero, diagonal elements are non-positive and off-diagonal elements
are non-negative (all entries in a row being zero would mean that the state remains fixed over
time). We will therefore sometimes refer to An as “Poisson transition matrix” or “intensity
matrix.” All this will be useful in section 2 below when we solve the Kolmogorov Forward
equation (3) and (4).
Summary of Algorithm. The algorithm is exactly the same as above, except that the
0
updating step uses (13) or equivalently (16). Guess vi,j , i = 1, ..., I, j = 1, 2 and for n = 0, 1, 2, ...
follow
n 0
1. Compute (vi,j ) using (8) and (11).
2. Compute cn from cni,j = (u0 )−1 [(vi,j
n 0
)]
3. Find v n+1 from (16).
4. If v n+1 is close enough to v n : stop. Otherwise, go to step 1.
7
One can show that this algorithm satisfies the three conditions of Barles and Souganidis (1991)
(monotonicity, consistency, stability) regardless of the size of ∆. See the discussion in Section
4 of Achdou et al. (2020). This is the big advantage of an implicit scheme over an explicit
scheme.
2 Kolmogorov Forward (Fokker-Planck) Equation
See matlab code huggett_partialeq.m We now turn to the solution of (3) and (4), which also
have to satisfy (5). The rough idea is to discretize these as
0 = −[si,j gi,j ]0 − λj gi,j + λ−j gi,−j (18)
I
X I
X
1= gi,1 ∆a + gi,2 ∆a (19)
i=1 i=1
(Instead of (19), one could also use a slightly more accurate trapezoidal rule, but results are
virtually identical given the fine grid size.) Because (3) and (4) are linear in g1 and g2 so is the
finite difference approximation. As a result, no iterative procedure like the one for the HJB
equation is needed and the equation can be solved in one step.
Upwind Scheme. There is again a question when to use a forward and a backward approxi-
mation for the derivative [si,j gi,j ]0 . It turns out that the most convenient/correct approximation
is as follows:
(sni,j,F )+ gi,j − gi−1,j (sni−1,j,F )+ gi+1,j (sni+1,j,B )− − gi,j (sni,j,B )−
− − − gi,j λj + gi,−j λ−j = 0
∆a ∆a
Note that because g0,j and gI+1,j are outside the state space, the density at these points is zero
and so (s0,j,F )+ and (sI+1,j,B )− are never used. The reason why the approximation above is
desirable is as follows. Collecting terms, we can write
gi−1,j zi−1,j + gi,j yi,j + gi+1,j xi+1,j + gi,−j λ−j = 0
(sni,j+1,B )−
xi+1,j = −
∆a
(si,j,F )+ (sni,j,B )−
n
yi,j = − + − λj
∆a ∆a
(sni,j−1,F )+
zi−1,j =
∆a
8
The reason this is the preferred approximation is that it can be written in matrix form in a
way that is closely related to the approximation used for the HJB equation
AT g = 0 (20)
where AT is the transpose of the intensity matrix A from the HJB equation (17) (An from
the final HJB iteration). This makes sense: the operation is exactly the same as that used
for finding the stationary distribution of a discrete Poisson process (continuous-time Markov
chain). The matrix A captures the evolution of the stochastic process and to find the stationary
distribution, one solves the eigenvalue problem AT g = 0. There is therefore a deep reason why
one wants to use the transpose of the intensity matrix A. For interested readers, this can be
made more precise using some tools from the theory of differential operators: one can write
the HJB equations (1) and (2) in terms of a differential operator A, the so-called “infinitesimal
generator” of the process. Similarly, the Kolmogorov Forward equations (3) and (4) can be
written in terms of an operator A∗ . An operator is the infinite-dimensional analogue of a matrix.
And the analogue of a matrix transpose is the so-called “adjoint” of an operator. It turns out
that the operator in the Kolmogorov Forward equation A∗ is the “adjoint” of the operator in the
HJB equation A. Putting things together, A is simply the discretized infinitesimal generator
whereas AT is the discretized version of its adjoint, the “Kolmogorov Forward operator.”
Besides making sense, this approximation is also convenient: once one has constructed the
matrix A for solving the HJB equation using an implicit method, almost no extra work is
needed.
To solve the eigenvalue problem (20) while imposing (19), the simplest procedure is as
follows. Fix gi,j = 0.1 (any other number will do as well) for an arbitrary (i, j), to then solve
the system for some g̃ and then to renormalize gi,j = g̃i,j /( Ii=1 g̃i,1 ∆a + Ii=1 g̃i,2 ∆a). Fixing
P P
gi,j = 0.1 is achieved by replacing the corresponding entry of the zero vector in (20) by 0.1,
and the corresponding row of AT by a row of zeros everywhere except for one on the diagonal.
Without this “dirty fix,” the matrix AT is singular and so cannot be inverted.
Alternatively, the eigenvalue problem (20) can be solved using a pre-built routine for nu-
merical eigenvalue problems. For example, MATLAB’s eigs function is well suited.5
As shown in Achdou et al. (2020), the wealth distribution of the low-income type g1 features
a Dirac mass at the borrowing constraint a (the left boundary of the state space). When
discretizing the distribution using a finite difference method, there is technically a Dirac mass
at every point in the state space. The algorithm therefore simply ignores the Dirac mass at the
boundary and treats it like any other point. Nevertheless, as we will see below, the numerical
solution has a clearly visible spike at a.
5
The particular command [g,val]=eigs(A’,1,’lr’) seems to work well. It returns the eigenvalue v and
eigenvector val of AT with the largest real part (which is v= 0). Finally, for very large problems (e.g. with
three or four state variables) iterative methods such as bicgstab may be preferable.
9
2.1 Results
Figure 2 (a) plots the functions s1 (a) (solid blue line), s2 (a) (solid green line). Note that
s01 (a) → −∞ as a → a as expected. Figure 2 (b) plots the associated densities g1 (a) and g2 (a).
s1 (a) g1 (a)
s2 (a) g2 (a)
Saving, sj (a) Densities, gj (a)
a Wealth, a a Wealth, a
Figure 2: Savings Policy Function and Implied Wealth Distribution
3 Equilibrium
3.1 Asset Supply
See matlab code huggett_asset_supply.m. After having solved (1) to (5), the asset supply
function S(r) defined in (6) can be easily computed. We approximate it as
I
X I
X
S(r) ≈ ai gi,1 ∆a + ai gi,2 ∆a
i=1 i=1
Figure 3 plots asset supply as a function of the interest rate. It looks as expected: in particular,
supply is bounded below by the borrowing constraint and S(r) → ∞ as r → ρ.
3.2 Finding the Equilibrium Interest Rate
See matlab code huggett_equilibrium_iterate.m. The equilibrium interest rate can easily
be found using a bisection method: the obvious idea is to increase r whenever S(r) < 0 and
decrease r whenever S(r) > 0. See the code for details.
10
r=ρ
a=a
r
S(r)
B
Figure 3: Asset Supply S(r)
4 Transition Dynamics and “MIT Shocks”
See matlab codes huggett_transition.m which needs input from huggett_initial.m and
huggett_terminal.m (the initial and terminal conditions). Besides solving transition dynamics
from an arbitrary initial condition, the same algorithm can be used to study the economy’s
impulse response after an “MIT shock,” i.e. an unanticipated (zero probability) shock followed
by a deterministic transition. We provide an example in Section 6.3.
The system to be solved is:
Z ∞ Z ∞
0 =S(r(t)) = ag1 (a, t)da + ag2 (a, t)da (21)
a a
ρvj (a, t) = max u(c) + ∂a vj (a, t)[zj + r(t)a − c] + λj [v−j (a, t) − vj (a, t)] + ∂t vj (a, t), (22)
c
∂t gj (a, t) = − ∂a [sj (a, t)gj (a, t)] − λj gj (a, t) + λ−j g−j (a, t) (23)
sj (a, t) =zj + r(t)a − cj (a, t), cj (a, t) = (u0 )−1 (∂a vj (a, t)) (24)
The bond market clearing condition can be written as
Z ∞ Z ∞
0= s1 (a, t)g1 (a, t)da + s2 (a, t)g2 (a, t)da
a a
We solve this system using the following algorithm. Guess a function r0 (t) and then for ` =
0, 1, 2, ... follow
1. Given r` (t), solve the HJB equation (22) with terminal condition vj` (a, T ) = vj (a) back-
ward in time to compute the time path of vj` (a, t). Also compute the implied saving policy
function s`j (a, t)
2. Given s`j (a, t), solve the Kolmogorov Forward equation (23) with initial condition gj (a, 0) =
11
gj0 (a) forward in time to calculate the time path for gj` (a, t).
3. Given s`j (a, t) and gj` (a, t) calculate
Z ∞ Z ∞
`
S (t) = ag1` (a, t)da + ag2` (a, t)da
a a
`
4. Update r`+1 (t) = r` (t) − ξ dSdt(t) , where ξ > 0.
5. Stop when r`+1 is sufficiently close to r` (t).
4.1 Solving the Time-Dependent HJB Equation (Step 1)
Approximate the value function at I discrete points in the wealth dimension and N discrete
n
points in the time dimension, and use the shorthand notation vi,j = vj (ai , tn ). The discrete
approximation to the time-dependent HJB (22) is
n+1 n
n 0
vi,j − vi,j
n
ρvi,j = u(cn+1
i,j ) + (vi,j ) [zj + r
n+1
ai − cn+1 n n
i,j ] + λj [vi,−j − vi,j ] + (25)
∆t
N
with terminal condition vi,j = vj (ai ). Given v n+1 , this system can be solved for v n exactly as
in Section 1.2. In particular, one can write this in matrix notation as
1 n+1
ρv n = un+1 + An+1 v n + (v − vn) (26)
∆t
where An+1 is defined in an analogous fashion to Section 1.2 and still has the interpretation
of the transition matrix of the discretized stochastic process for (at , zt ). Now each n has the
interpretation of a time step instead of an iteration on the stationary value function. The reason
for this similarity in the algorithm is that intuitively a stationary value function can be found
by solving a time-dependent problem and going far enough back in time, i.e. as t → −∞.
4.2 Solving the Time-Dependent Kolmogorov Forward Eq. (Step 2)
Analogously to the value function, we approximate the density at J discrete points in the
wealth dimension and N discrete points in the time dimension, and use the shorthand notation
n
gi,j = gj (ai , tn ). Similarly to Section 2 one can directly make use of the transition matrix An
defined when solving the time-dependent HJB equation (Section 4.1). Given an initial condition
0
gi,j = gj,0 (ai ), the Kolmogorov Forward equation (23) is then easily solved. One here has the
option of using either an explicit method
g n+1 − g n
= (An )T g n ⇒ g n+1 = ∆t(An )T g n + g n (27)
∆t
12
or an implicit method
g n+1 − g n
= (An )T g n+1 ⇒ g n+1 = (I − ∆t(An )T )−1 g n .
∆t
Note that these schemes preserve mass: starting from any initial distribution g 0 that sums to
one, all future g n ’s also sum to one. This follows from the fact that the rows of the intensity
matrices An sum to zero. The implicit scheme is also guaranteed to preserve the positivity of
g for arbitrary time steps ∆t.
4.3 Results: Transition Dynamics
Figure 4 plots the time path for the equilibrium interest rate in response to a permanent increase
in “unemployment risk”, λ2 . Figure 5 plots the densities g1 (a, t) and g2 (a, t) at various points in
time during the transition. For comparison, the Figure also plots densities in the initial steady
state (dashed lines).
Equilibrium Interest Rate, r(t)
0.03
0.02
0.01
0
−0.01
−0.02
0 2 4 6 8 10
Year
Figure 4: Time Path of Equilibrium Interest Rate
13
t = 0.1 t=2
3 3
g (a) g (a)
1 1
2.5 g (a) 2.5 g (a)
2 2
D e ns it ie s , g i ( a) D e ns it ie s , g i ( a)
2 2
1.5 1.5
1 1
0.5 0.5
0 0
−0.1 0 0.1 0.2 0.3 0.4 0.5 −0.1 0 0.1 0.2 0.3 0.4 0.5
We alt h, a We alt h, a
t=5 t=∞
3 3
g (a) g (a)
1 1
2.5 g2(a) 2.5 g2(a)
D e ns it ie s , g i ( a) D e ns it ie s , g i ( a)
2 2
1.5 1.5
1 1
0.5 0.5
0 0
−0.1 0 0.1 0.2 0.3 0.4 0.5 −0.1 0 0.1 0.2 0.3 0.4 0.5
We alt h, a We alt h, a
Figure 5: Dynamics of Wealth Distribution
5 Generalization to Diffusion Process
The system to be solved is:
σ 2 (z)
ρv(a, z) = max u(c) + ∂a v(a, z)[z + ra − c] + µ(z)∂z v(a, z) + ∂zz v(a, z) (28)
c 2
1
0 = −∂a [s(a, z)g(a, z)] − ∂z [µ(z)g(a, z)] + ∂zz [σ 2 (z)g(a, z)] (29)
Z ∞Z ∞ 2
1= g(a, z)dadz (30)
0 a
Z ∞Z ∞
0= ag(a, z)dadz ≡ S(r) (31)
0 a
14
We assume that the z-process gets reflected at some z and z̄. One can show that this gives rise
to the following boundary conditions for v:6
0 = ∂z v(a, z) = ∂z v(a, z̄)
Finally, we have the state constraint boundary condition:
∂a v(a, z) ≥ u0 (z + ra), all z. (32)
We again use a finite difference method and use the short-hand notation v(ai , zj ) = vi,j . Note
that we changed notation slightly and i now indexes wealth and j indexes productivity.
5.1 HJB Equation
See matlab program HJB_diffusion_implicit.m. With a diffusion process, an explicit method
becomes extremely inefficient so we here only explain the solution of the HJB with an implicit
method. The derivative in the a dimension is again approximated using an upwind method,
i.e. using either a forward or a backward difference approximation depending on the sign of
the drift:
vi,j − vi−1,j
∂a,B vi,j =
∆a (33)
vi+1,j − vi,j
∂a,F vi,j =
∆a
Similarly, we also use an upwind method in the z-direction. For the second-order derivative,
we use a central difference approximation. Hence:
vi,j − vi,j−1
∂z,B vi,j =
∆z
vi,j+1 − vi,j
∂z,F vi,j =
∆z
vi,j+1 − 2vi,j + vi,j−1
∂zz vi,j =
(∆z)2
Analogously to the model with Poisson shocks, v n+1 is now implicitly defined by the equation
n+1 n
vi,j − vi,j n+1 n n+1 n n+1
σj2 n+1
+ ρvi,j = u(ci,j ) + ∂a vi,j [zj + rai − ci,j ] + µj ∂z vi,j + ∂zz vi,j
∆ 2
Note the n + 1 superscripts on the right-hand side of the equation. The main advantage of the
implicit scheme is that the step size ∆ can be arbitrarily large.
6
See e.g. Section 3.5 in Dixit (1993).
15
Upwind Scheme. We again need to use an “upwind scheme.” As above, the idea is still to
use the forward difference approximation whenever the drift of the state variable is positive
and the backward difference approximation whenever it is negative. We use the following finite
difference approximation to (28).
n+1 n
vi,j − vi,j n+1
+ ρvi,j n+1
=u(cni,j ) + ∂a,F vi,j n+1
[zj + rai − cni,j,F ]+ + ∂a,B vi,j [zj + rai − cni,j,B ]−
∆ (34)
n+1 + n+1 −
σj2 n+1
+ ∂z,F vi,j µj + ∂z,B vi,j µj + ∂zz vi,j
2
Equation (34) constitutes a system of I × J linear equations, and it can be written in matrix
notation using the following steps. Substituting the definition of the derivatives (33), and
defining sni,j,F = zj + rai − cni,j,F and similarly for sni,j,B , (34) is
n+1 n n+1 n+1 n+1 n+1
vi,j − vi,j vi+1,j − vi,j vi,j − vi−1,j
n+1
+ ρvi,j =u(cni,j ) + (sni,j,F )+ + (sni,j,B )−
∆ ∆a ∆a
n+1 n+1 n+1 n+1 n+1 n+1 n+1
vi,j+1 − vi,j vi,j − vi,j−1 − σj2 vi,j+1 − 2vi,j + vi,j−1
+ µ+
j + µ j +
∆z ∆z 2 (∆z)2
Collecting terms with the same subscripts on the right-hand side
n+1 n
vi,j − vi,j n+1 n+1 n+1 n+1 n+1 n+1
+ ρvi,j = u(cni,j ) + vi−1,j xi,j + vi,j (yi,j + υj ) + vi+1,j zi,j + vi,j−1 χj + vi,j+1 ζj
∆
(sni,j,B )−
xi,j = − ,
∆a
(sni,j,F )+ (sni,j,B )−
yi,j = − + ,
∆a ∆a
(sni,j,F )+
zi,j = (35)
∆a
− 2
µj σj
χj = − +
∆z 2(∆z)2
µ−j µ+j σj2
υj = − −
∆z ∆z (∆z)2
µ+ j σj2
ζj = +
∆z 2(∆z)2
n+1 n+1
Note that importantly x1,j = zI,j = 0 for all j so v0,j and vI+1,j are never used. At the
boundaries in the j dimension, the equations become
n+1 n
vi,1 − vi,1 n+1 n+1 n+1 n+1 n+1
+ ρvi,1 = u(cni,1 ) + vi−1,1 xi,1 + vi,1 (yi,1 + υ1 + χ1 ) + vi+1,1 zi,1 + vi,2 ζ1
∆
n+1 n
vi,J − vi,J n+1 n+1 n+1 n+1 n+1
+ ρvi,J = u(cni,J ) + vi−1,J xi,J + vi,J (yi,J + υJ + ζJ ) + vi+1,J zi,J + vi,J−1 χJ
∆
16
v −v
where, in the first equation, we have used that ∂z,B vi,1 = i,1∆z i,0 = 0 and hence vi,0 = vi,1 .
v −vi,J
Similarly, in the second equation, ∂z,F vi,J = i,J+1
∆z
= 0 and hence vi,J+1 = vi,J . Note that
we here defined the boundary conditions relative to the points j = 0 and j = J + 1 and used
the values vi,0 and vi,J+1 . These points are sometimes called “ghost nodes”. Equation (35) is a
system of I × J linear equations which can be written in matrix notation as:
1 n+1
(v − v n ) + ρv n+1 = un + An v n+1 (36)
∆
where v n is a vector of length I × J with entries (v1,1 , ..., vI,1 , v1,2 , ..., vI,2 , ..., vI,J ) and An =
e n + C where the (I × J) × (I × J) matrices A
A e n and C are
y1,1 z1,1 0 ··· ··· ··· ··· ··· ··· ··· ··· ··· 0
... ... ... ... ... ... ... ... ..
x2,1 y2,1 z2,1 0 .
.. .. .. .. .. .. .. .. .. .. .. ..
0 . . . . . . . . . . . .
.
.. .. .. .. .. .. .. .. ..
. xI,1 yI,1 0 0 . . . . . . .
.
.. .. .. .. .. .. .. ..
. 0 0 y1,2 z1,2 0 . . . . . .
.
.. .. .. .. .. .. .. ..
. . 0 x2,2 y2,2 z2,2 0 . . . . .
e n = ...
A .. .. .. .. .. .. .. .. .. .. .. ..
. . . . . . . . . . .
.
.. .. .. .. .. .. .. .. ..
. . . . . 0 xI,2 yI,2 0 . . . .
.
.. .. .. .. .. .. .. .. .. .. .. .. ..
. . . . . . . . . . . .
.
.. .. .. .. .. .. .. .. ..
. . . . . . . 0 y1,J z1,J 0 .
. . . . . . . .
.. .. .. .. .. .. .. . . 0 x2,J y2,J z2,J
0
.. .. .. .. .. .. .. .. .. .. .. .. ..
. . . . . . . . . . . . .
0 ... ... ... ... ... ... ... ... ... 0 xI,J yI,J
17
υ1 + χ1 0 ··· ··· 0 ζ1 0 · · · · · · ··· ··· ··· 0
... ... ... ... ... ... ..
0 υ1 + χ 1 0 0 ζ1 0 .
.. ... ... ... ... ... ... ... ... ... ... ... ..
. .
.. .. .. .. .. ..
0
. 0 υ1 + χ 1 0 . . 0 ζ1 0 . . .
. .. .. .. .. ..
χ
2 0 0 υ2 0 . . 0 ζ2 0 . .
. .. .. .. ..
0
χ2 0 0 υ2 0 . . 0 ζ2 0 .
C = ... ... ..
.
.. .. .. ..
. . . .
..
.
..
.
..
.
..
0 0 .
.. .. .. .. .. ..
. . 0 χ2 0 . 0 υ2 0 . . 0 .
. ..
.. ... ... .. .. .. .. .. .. .. .. ..
. . . . . . . . . .
.
.. . .. . .. . .. .. .. .. ..
0 χJ 0 . . υJ + ζJ 0 . .
.
.. .. .. .. .. .. ..
. . . . 0 χJ 0 . 0 υJ + ζJ 0 .
.. . .. . .. . .. .. .. .. .. .. .. .. ..
. . . . . . . .
. 0
0 ... ... ... . . . . . . . . . 0 χJ 0 ... 0 υJ + ζJ
Equation (36) can again be solved very efficiently in matlab. Figure 6 again plots a visualization
of the intensity matrix A in practice (using Matlab’s spy function).
0
500
1000
1500
2000
2500
3000
3500
4000
0 500 1000 1500 2000 2500 3000 3500 4000
nz = 15749
Figure 6: Visualization of the matrix A using spy in model with diffusion process
18
5.2 Kolmogorov Forward Equation and Equilibrium
See matlab program huggett_diffusion_partialeq.m. The Kolmogorov Forward equation is
solved exactly as in section 2, and the equilibrium is found as in section 3.
5.3 Results
Figures 7 and 8 plot the functions s(a, z) and g(a, z).
0.6
0.4
0.2
Savings s(a,z)
0
-0.2
-0.4
-0.6
-0.8
0 1.5
5
10
1
15
20
Wealth, a 25 0.5 Productivity, z
Figure 7: Savings Policy Function in Huggett Model with Diffusion
6 Aiyagari Model
We now briefly explain how to solve the Aiyagari model in Section 6 of Achdou et al. (2020).
As in the paper we focus on the case where productivity follows a diffusion process so as
to explain how to handle that case. Of course, it is also straightforward to solve an Aiya-
gari model in which income follows a two-state Poisson process. Codes for this case are
also available: aiyagari_poisson_steadystate.m, aiyagari_poisson_asset_supply.m and
aiyagari_poisson_MITshock.m.
19
0.25
0.2
Density g(a,z)
0.15
0.1
0.05
0
0 1.5
2
4 1
6
Wealth, a 8 0.5 Productivity, z
Figure 8: Wealth Distribution in Huggett Model with Diffusion
6.1 Steady State
See matlab program aiyagari_diffusion_equilibrium.m, Julia program aiyagari_diffusion_
equilibrium.jl and C++ program aiyagari_diffusion_equilibrium.cpp. A steady state
or stationary equilibrium can be represented by the following system of equations which we
aim to solve numerically:
1
ρv(a, z) = max u(c) + ∂a v(a, z)(wz + ra − c) + ∂z v(a, z)µ(z) + ∂zz v(a, z)σ 2 (z) (37)
c 2
1 2
0 =∂a (s(a, z)g(a, z)) − ∂z (µ(z)g(a, z)) + ∂zz (σ (z)g(a, z)) (38)
2
r = ∂K F (K, 1) − δ, w = ∂L F (K, 1), (39)
Z z̄ Z ∞
K= ag(a, z)dadz (40)
z a
on (a, ∞)×(z, z̄), where s(a, z) = wz+ra−c(a, z), c(a, z) = (u0 )−1 (∂a v(a, z)) and with boundary
conditions
u0 (wz + ra) ≥ ∂a v(a, z), all z
∂z v(a, z) = 0, ∂z v(a, z̄) = 0, all a
20
The algorithm for solving the HJB and KF equations is the same as in Sections 5.1 and 5.2.
To find the equilibrium wage and interest rate w and r, we use a fixed point algorithm on the
scalar K. Alternatively, one can express w as a function of r using (39) and (40) and use a
bisection method to solve for the equilibrium r.
6.2 Transition Dynamics
See matlab program aiyagari_diffusion_transition.m. The algorithm for solving the HJB
and KF equations is the natural generalization to the time-dependent case of that outlined in
Sections 5.1 and 5.2. To solve for the equilibrium time paths of the wage and interest rate, we
use a fixed point algorithm on the function K(t).
6.3 “MIT Shocks”
See matlab program aiyagari_poisson_MITshock.m. We consider the version of the Aiyagari
model with Poisson income shocks and compute the impulse response to a negative aggregate
productivity shock that mean reverts over time. This shock is modeled as an “MIT shock,” i.e.
an unexpected (zero probability) shock followed by a deterministic transition. More precisely,
we assume that the aggregate production function is Yt = Ft (K, L) = At K α L1−α and aggregate
productivity follows a deterministic version of an Ornstein-Uhlenbeck process (the continuous-
time analogue of an AR(1) process):
dAt = ν(Ā − At )dt.
The parameter ν governs the speed of mean reversion (one can show that Corr(At , At+s ) =
e−νs ). Figure 9 plots the impulse response to this productivity MIT-shock. In particular note
the measures of income and wealth inequality in the last two panels.
6.4 Visualizing Evolution of Wealth Distribution as Movie
After running aiyagari_diffusion_transition.m or aiyagari_poisson_MITshock.m, you
can run make_movie.m to make the movie of the transition of the distribution.
After you run make_movie.m, you will have distribution.avi file. One can use ffmpeg
(available at https://www.ffmpeg.org/) to convert from avi file to mp4 file. In particular
type
ffmpeg -i distribution.avi -c:v libx264 -g 30 -pix_fmt yuv420p distribution.mp4
and this will create a mp4 file that is accepted by pretty much anything (and it works on
website with html5).7
7
”-pix_fmtyuv420p” is only necessary while browsers and media players are being updated. It is using old
format because ffmpeg was updated recently, and most media players have not been.
21
Aggregate Productivity Aggregate Capital Stock Wage
0.1005 0.3 0.039
0.1 0.299 0.0388
0.0995 0.298 0.0386
0.099 0.297 0.0384
0.0985 0.296 0.0382
0.098 0.295 0.038
0.0975 0.294 0.0378
0.097 0.293 0.0376
0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100
Year Year Year
Interest Rate Wealth Gini Income Gini
0.049 0.896 0.242
0.0485
0.894
0.241
0.048
0.892
0.0475 0.24
0.047 0.89
0.239
0.0465 0.888
0.046 0.238
0.886
0.0455
0.237
0.884
0.045
0.0445 0.882 0.236
0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100
Year Year Year
Figure 9: Impulse Response to Negative Productivity MIT-shock in Aiyagari Model
If you do not want to go through the pain of compiling your own ffmpeg binary and have a
Mac, you can download a pre-compiled binary at http://ffmpegmac.net/. Copy the binary
ffmpeg into the same directory as your avi file and run (don’t forget the “./” in the beginning)
./ffmpeg -i distribution.avi -c:v libx264 -g 30 -pix_fmt yuv420p distribution.mp4
The resulting movie is at http://www.princeton.edu/~moll/HACTproject/distribution.
mp4.
7 Non-uniform Grids
In all previous sections, we worked with uniformly spaced grids. But for many applications, we
may want to economize on grid points. This can be achieved by working with non-uniform grids
and “putting grid points” at points in the state space where the value and density functions
have the most curvature.
22
7.1 HJB equation with non-uniform grid
Extending our algorithm for the HJB equation to the case of non-uniform grids is straightfor-
ward. Denoting by ∆ai,+ = ai+1 − ai and ∆ai,− = ai − ai−1 , the forward and backward distance
between two grid points, we simply change (8) to
vi+1,j − vi,j vi+1,j − vi,j
vj0 (ai ) ≈ = 0
≡ vi,j,F
ai+1 − ai ∆ai,+
(41)
vi,j − v i−1,j vi,j − vi−1,j
vj0 (ai ) ≈ = 0
≡ vi,j,B
ai − ai−1 ∆ai,−
Following the same steps as above, we end up with
n+1 n
vi,j − vi,j n+1 n+1 n+1 n+1 n+1
+ ρvi,j = u(cni,j ) + vi−1,j xi,j + vi,j yi,j + vi+1,j zi,j + vi,−j λj where
∆
(42)
(sni,j,B )− (sni,j,F )+ (sni,j,B )− (sni,j,F )+
xi,j = − , yi,j = − + − λj , zi,j =
∆ai,− ∆ai,+ ∆ai,− ∆ai,+
This can again be written in matrix form (15) where the intensity matrix An has the same
structure as above but with the entries in (42). The rest of the algorithm is unchanged.
If an approximation to the second derivative of vj is needed as in Section 8, a good candidate
approximation is
∆ai,− vi+1,j − (∆ai,− + ∆ai,+ )vi,j + ∆ai,+ vi−1,j
vj00 (ai ) ≈ 1 (43)
2
(∆ai,+ + ∆ai,− ) ∆ai,− ∆ai,+
This approximation can be derived from a Taylor approximation to vj 8 and it can be seen that
with ∆ai,− = ∆ai,+ = ∆a, it reduces to the standard second-derivative approximation in the
case with uniform grids:
vi+1,j − 2vi,j + vi−1,j
vj00 (ai ) ≈
(∆a)2
8
Consider a second-order Taylor approximation to vj around ai :
1
vi+1,j ≈ vi,j + ∆ai,+ vj0 (ai ) + (∆ai,+ )2 vj00 (ai )
2
1
vi−1,j ≈ vi,j − ∆ai,− vj0 (ai ) + (∆ai,− )2 vj00 (ai )
2
Multiply the first equation by ∆ai,− and the second equation by ∆ai,+ and add the two equations
1
∆ai,− vi+1,j + ∆ai,+ vi−1,j ≈ (∆ai,− + ∆ai,+ )vi,j + (∆ai,+ + ∆ai,− ) ∆ai,− ∆ai,+ vj00 (ai )
2
Rearranging yields (43).
23
7.2 Kolmogorov Forward equation with non-uniform grid
Extending the Kolmogorov Forward equation to the case of a non-uniform grid requires more
work. Section 2 suggests working with the transpose of the intensity matrix A and to solve
g n+1 − g n
= AT g n
∆t
or the implicit analogue. The problem with this scheme is that it is not guaranteed to preserve
mass: starting with an initial distribution g that integrates/sums to one, the total mass may
converge to a different number over time (including zero or infinity).
To see this consider a simplified example with one income type only z1 = z2 in which case
the intensity matrix is of size I × I with entries Ai,i0 given by
(sni,B )− (sni,F )+ (sni,B )− (sni,F )+
Ai,i−1 = xi = − , Ai,i = yi = − + , Ai,i+1 = zi = (44)
∆ai,− ∆ai,+ ∆ai,− ∆ai,+
P
From (44), we can see that the rows of A still sum to zero, i0 Ai,i0 = 0. Therefore
I I I
!
X g n+1 − g n
i i
X X
= Ai0 ,i gi0 = 0
i=1
∆t i0 =1 i=1
P gn+1 −gn
In the case of a uniform grid, this also implies mass preservation Ii=1 i ∆t i ∆a = 0. However,
in the case of a non-uniform grid what we want instead is that an appropriate approximation
P gn+1 −gn
of the integral equals zero, i.e. something like Ii=1 i ∆t i ∆ai = 0 (where the ∆ai ’s depend
on the precise integral approximation method).
The following solution to this problem seems to work well in practice, in particular it
preserves mass and the positivity of g.9 First, approximate the integral of g with the trapezoidal
rule 10
Z amax I−1 I
n 1X X
˜ i
g(a, t )da ≈ (ai+1 − ai )(gi+1 + gi ) = gi ∆a
a 2 i=1 i=1
1
2 ∆ai,+ , i=1
˜ i = 1 (∆ai,+ + ∆ai,− ),
∆a i = 2, ..., I − 1
2
1
2
∆ai,− , i=I
The key idea is now that, rather than working with the vector g with elements gi , we work
˜ i and which must therefore satisfy
directly with a vector g̃ whose elements are g̃i = gi ∆a
9
We have not yet been able to derive a rigorous theoretical justification for the proposed scheme.
10
See http://en.wikipedia.org/wiki/Trapezoidal_rule#Non-uniform_grid
24
PI 1
i=1 g̃i = 1. Given an initial condition g̃ we simply solve the following analogue of (27):
g̃ n+1 − g̃ n
= AT g̃ n .
∆t
Following the same logic as above, the condition Ii0 =1 Ai,i0 = 0 guarantees mass preservation
P
PI PI ˜ n
i=1 g̃i = i=1 gi ∆ai = 1. We can then always back out the true distribution g from gi =
˜ i ). In matrix form, we have g̃ = Dg where D is a diagonal matrix with elements
g̃i /(∆a
∆ãi , i = 1, ..., I. Therefore underlying distribution g can also be found from g = D−1 g̃.11
It is also straightforward to show that the same approach – working with the rescaled density
g̃ = Dg – also applies to the case with multiple income states (or a continuum).
8 Aiyagari Model with Fat-tailed Wealth Distribution
See matlab code fat_tail_partialeq.m. This section shows how to extend our computational
methods to the model with two assets and a fat-tailed wealth distribution from Section 6 of
Achdou et al. (2020).
8.1 Model Setup
The system of equations to be solved is
1
ρvj (a) = max u(c) + vj0 (a)(zj + ra + (R − r)k − c) + vj00 (a)σ 2 k 2 + λj (v−j (a) − vj (a)) (45)
c,k≤a+φ 2
2
d 1 d
0 = − [sj (a)gj (a)] + 2
[σ 2 kj (a)2 gj (a)] − λj gj (a) + λ−j g−j (a) (46)
da 2 da
Z ∞ Z ∞ Z ∞ Z ∞
k1 (a)g1 (a)da + k2 (a)g2 (a)da = ag1 (a)da + ag2 (a)da (47)
a a a a
11
A previous version of this Appendix advocated working with a rescaled version of the intensity matrix A
rather than a rescaled version of the density g: replace the intensity matrix A in (20) or (27) with an alternative
PI ˜ i0 = 0 and then solve
intensity matrix à that satisfies i0 =1 Ãi,i0 ∆a
g n+1 − g n
= ÃT g n .
∆t
The rescaled intensity matrix was given by à = DAD−1 where again D is a diagonal matrix with elements
∆ãi , i = 1, ..., I. The two approaches are equivalent given that g̃ = Dg. However, we prefer the approach in
the text because it is much more transparent and easier to explain. We are grateful to Matthieu Gomez for
suggesting the approach in the current version.
25
Optimal consumption and choice of risky assets are
cj(a) = vila) (48)
k;(a) = min { a “ a+ a} . (49)
Boundary Conditions In theory, the HJB equation (45) is defined on (a, co) but in practice
it has to be solved on a bounded interval (@,dmax). A non-trivial issue concerns the question
what boundary condition to impose at Gpax. We use the asymptotic behavior of the value
function in Lemma 2 in the paper to motivate boundary conditions as follows.’* For large a,
we have
~ ~ _
v;(a) = %oj tt j;a 7
for unknown constants Up; and 01,;. Hence, we impose the following boundary condition
N
Uj (Gmax) = —7U} (max) /Gmax- (50)
To solve (45), what we really need is a boundary condition for the term xu" (a)k(a)*. From
(49) and (50)
R-r
k; (Gmax) = “yo2 mex (51)
a Qn a 24
7 Ki (dina) U5 (imax) = ~ > ki (dinax) VU; (Gmax) / max
2 (52)
oy _ (R 7 r)
— VU; (Gmax)§, g —_ Dy? max
We will condition (52) below when solving (45) using a finite difference method. Finally, it
sometimes helps numerical stability to impose a state constraint @ < dmax. This is equivalent
to Cj(@max) > 2% +1TQmax + (R—1)kj(Gmax) or using (51)
(R _ r)? ~yY
U;(Qmax) < (:, + T Amax + yee max :
8.2 Finite Difference Method for HJB Equation
The steps follow closely the solution method of the one-asset model. We therefore only outline
the main differences. As before we use an implicit upwind method. In contrast to before,
the HJB equation (45) now involves the second derivative of the value function. For sake of
transparency, we here explain our finite difference method for a uniform grid. The code in
fat_tail_partialeq.m instead uses a non-uniform grid as explained in Section 7 so as to be
!2We thank Matthieu Gomez for suggesting this boundary condition.
26
able to put more points in the region of the state space where policy functions have a lot of
curvature, i.e. close to the borrowing constraint.
Defining sni,j,F = zj + (R − r)ki,j + rai − cni,j,F and similarly for sni,j,B , the discretization of
(45) is
n+1 n n+1 n+1 n+1 n+1
vi,j − vi,j vi+1,j − vi,j vi,j − vi−1,j
n+1 n
+ ρvi,j = u(ci,j ) + n +
(si,j,F ) + (sni,j,B )−
∆ ∆a ∆a
n+1 n+1 n+1
n+1 n+1 σ 2 2 vi+1,j − 2vi,j + vi−1,j
+ λj [vi,−j − vi,j ] + ki,j
2 (∆a)2
Collecting terms with the same subscripts on the right-hand side:
n+1 n
vi,j − vi,j n+1 n+1 n+1 n+1 n+1
+ ρvi,j = u(cni,j ) + vi−1,j xi,j + vi,j yi,j + vi+1,j zi,j + vi,−j λj where
∆
(sni,j,B )− σ 2 ki,j2
xi,j = − + ,
∆a 2 (∆a)2
(53)
(sni,j,F )+ (sni,j,B )− 2
2
ki,j
yi,j = − + −σ − λj ,
∆a ∆a (∆a)2
(sni,j,F )+ σ 2 ki,j 2
zi,j = +
∆a 2 (∆a)2
At the upper boundary a = amax = aI , we make use of (52) and write the approximation as
n+1 n n+1 n+1 n+1 n+1
vI,j − vI,j vI+1,j − vI,j vI,j − vI−1,j
n+1 n
+ ρvI,j = u(cI,j ) + n +
(sI,j,F ) + (snI,j,B )−
∆ ∆a ∆a
n+1 n+1
n+1 n+1
vI,j − vI−1,j
+ λj [vI,−j − vI,j ] + ξ
∆a
so that the corresponding entries of (53) become.
(snI,j,B )− ξ
xI,j = − − ,
∆a ∆a
(snI,j,F )+ (snI,j,B )− ξ
yI,j = − + + − λj , (54)
∆a ∆a ∆a
(snI,j,F )+
zI,j =
∆a
Equations (53) and (54) is a system of 2 × I linear equations which can be written in matrix
notation like equation (15)
1 n+1
(v − v n ) + ρv n+1 = un + An v n+1
∆
and that can be solved efficiently.
27
8.3 Finite Difference for Kolmogorov Forward Equation
The solution of the Kolmogorov Forward equation (46) is exactly as in Section 2, with one
difference: because of the second-order term one has to decide what to do at the upper end
of the state space amax . The cleanest solution is to impose an artificial reflecting barrier. To
this end, consider the “intensity matrix” A. Rather than using (54) at the upper end of
the state space, construct the transition matrix according to (53). But then move all entries
corresponding to the (non-existent) grid point I + 1 to the entry corresponding to I:
(snI,j,B )− σ 2 kI,j 2
x̃I,j = xI,j = − + ,
∆a 2 (∆a)2
(snI,j,B )− σ 2 kI,j2
(55)
ỹI,j = yI,j + zI,j = − − λj ,
∆a 2 (∆a)2
z̃I,j = 0.
The interpretation is that whenever the process would leave the state space according to the
discretized law of motion (if it would go to point I + 1), it is “reflected” back in (back down to
point I).13
8.4 Results
See Figure 10.
9 Accuracy of Finite Difference Method
See Appendix F.1 of Achdou et al. (2020) available at https://benjaminmoll.com/HACT_
appendix/ for various accuracy checks.
13
This condition can be derived more rigorously as follows. For simplicity consider the case without produc-
tivity shocks z1 = z2 so that the process for wealth is dat = s(at )dt + σk(at )dWt . Impose a reflecting barrier at
a = amax . Then the infinitesimal generator corresponding to this process is given by
σ2
(Af )(a) = s(a)f 0 (a) + k(a)2 f 00 (a)
2
with the boundary condition corresponding to a reflecting barrier: f 0 (amax ) = 0. Its discrete version is:
fi+1 − fi fi − fi−1 σ 2 2 fi+1 − 2fi + fi−1
(Af )i = s+ + s− + k
i
∆a i
∆a 2 i (∆a)2
= xi fi−1 + yi fi + zi fi−1
where xi , yi , zi are analogous to (53). The discretized boundary condition is f 0 (aI ) ≈ (fI+1 − fI )/(∆a) = 0 or
fI+1 = fI . Therefore (Af )I = xI fI−1 + (yI + zI )fI = x̃I fI + ỹI fI with ỹI = yI + zI i.e. just like in (55).
28
s1 (a) k1 (a)
s2 (a) k2 (a)
Risky Assets, ki (a)
Savings, si (a)
a Wealth, a a Wealth, a
(a) Saving Policy Function (b) Investment in Risky Assets
g1 (a) log(f 1 (x))
g2 (a) log(f 2 (x))
Slope = -ζ
Densities, gi (a) log (fi (x))
a Wealth, a Log Wealth, x = log (a)
(c) Stationary Densities (d) Pareto Tail of Stationary Distribution
Figure 10: Optimal Choices and Pareto Tail of Wealth Distribution in Two-Asset Model
References
Achdou, Yves, Jiequn Han, Jean-Michel Lasry, Pierre-Louis Lions, and Benjamin
Moll. 2020. “Income and Wealth Distribution in Macroeconomics: A Continuous-Time Ap-
proach.” London School of Economics Working Paper.
Barles, G., and P. E. Souganidis. 1991. “Convergence of approximation schemes for fully
nonlinear second order equations.” Asymptotic Analysis, 4: 271–283.
Candler, Graham V. 1999. “Finite-Difference Methods for Dynamic Programming Prob-
lems.” In Computational Methods for the Study of Dynamic Economies. Cambridge, Eng-
land:Cambridge University Press.
29
Dixit, Avinash. 1993. The Art of Smooth Pasting. Fundamentals of Pure and Applied Eco-
nomics 55, The Routledge.
Huggett, Mark. 1993. “The risk-free rate in heterogeneous-agent incomplete-insurance
economies.” Journal of Economic Dynamics and Control, 17(5-6): 953–969.
30
File and source
- File
- countercyclical-unemployment-benefits-a-general-equilibrium-analysis-of-transition-dynamics.pdf
- Size
- 1,407,523 bytes
- SHA-256
- 214d6347a8f8ff1883dc722358038c5f4dc2453089cc0087fa9b46987e2f5d08
- Our copy
- countercyclical-unemployment-benefits-a-general-equilibrium-analysis-of-transition-dynamics.pdf
- Original
- doi.org