Income and Wealth Distribution in Macroeconomics: A Continuous-Time Approach

Abstract

We recast the Aiyagari–Bewley–Huggett model of income and wealth distribution in continuous time. This workhorse model—as well as heterogeneous agent models more generally—then boils down to a system of partial differential equations, a fact we take advantage of to make two types of contributions. First, a number of new theoretical results: (1) an analytic characterization of the consumption and saving behaviour of the poor, particularly their marginal propensities to consume; (2) a closed-form solution for the wealth distribution in a special case with two income types; (3) a proof that there is a unique stationary equilibrium if the intertemporal elasticity of substitution is weakly greater than one. Second, we develop a simple, efficient and portable algorithm for numerically solving for equilibria in a wide class of heterogeneous agent models, including—but not limited to—the Aiyagari–Bewley–Huggett model.

1. Introduction

One of the key developments in macroeconomics research over the last three decades has been the incorporation of explicit heterogeneity into models of the macroeconomy. Fuelled by the increasing availability of high-quality micro data, the advent of more powerful computing methods as well as rising inequality in many advanced economies, such heterogeneous agent models have proliferated and are now ubiquitous. This is a welcome development for a number of reasons. First, it opens up the door to bringing micro data to the table in order to empirically discipline macro theories. Second, macroeconomists often want to analyse the welfare implications of particular shocks or policies. This is impossible without asking “who gains and who loses?”, that is, distributional considerations often cannot be ignored. Third, models with heterogeneity often deliver strikingly different aggregate implications than do representative agent models, for example, with respect to monetary and fiscal policies.¹

Despite the continuously increasing popularity of macroeconomic models with rich heterogeneity, the literature has suffered from a dearth of theoretical and analytical results. Little is known about the properties of consumption and saving behaviour in the presence of borrowing constraints, those of the resulting wealth distribution, and equilibrium uniqueness (or lack thereof). Instead, most studies rely on purely numerical analyses to characterize the implications of such theories. But even such computational approaches are often difficult and costly, particularly if the question at hand requires solving for the economy’s transition dynamics or if the model features non-differentiabilities or non-convexities.

In this article, we make some progress on these issues by recasting the standard incomplete market model of Aiyagari (1994), Bewley (1986), and Huggett (1993) in continuous time.² Our main contributions are twofold. First, we prove a number of new theoretical results about this workhorse model.³ Second, we develop a simple, efficient, and portable algorithm for numerically solving both stationary equilibria and transition dynamics of a wide class of heterogeneous agent models, including—but not limited to—the Aiyagari–Bewley–Huggett model.

Both types of contributions make use of an important property: when recast in continuous time, heterogeneous agent models boil down to systems of two coupled partial differential equations. The first of these is a Hamilton–Jacobi–Bellman (HJB) equation for the optimal choices of a single atomistic individual who takes the evolution of the distribution and hence prices as given. And the second is a Kolmogorov Forward (KF) equation characterizing the evolution of the distribution, given optimal choices of individuals.⁴ More generally, our approach is to cast heterogeneous agent models in terms of the mathematical theory of “Mean Field Games” (MFG) initiated by Lasry and Lions (2007).⁵ The system of coupled HJB and KF equations is known as the “backward–forward MFG system.”

In the context of the Aiyagari–Bewley–Huggett model, the HJB equation characterizes individuals’ optimal consumption and saving behaviour given a stochastic process for income; and the KF equation characterizes the evolution of the joint distribution of income and wealth. The two equations are coupled because optimal consumption and saving depend on the interest rate which is determined in equilibrium and hence depends on the wealth distribution. We start with a particularly parsimonious case: a Huggett (1993) economy in which idiosyncratic income risk takes the form of exogenous endowment shocks that follow a two-state Poisson process and in which individuals save in unproductive bonds that are in fixed supply. Later in the article, we extend many of our results to more general stochastic processes and to an Aiyagari (1994) economy in which individuals save in productive capital.

We prove three new theoretical results about the Aiyagari–Bewley–Huggett model. First, we provide an analytic characterization of the consumption and saving behaviour of the poor. We show that, under natural assumptions, an individual’s saving policy function behaves like

$$-\sqrt{2 \nu a}$$

in the vicinity of the borrowing constraint, where

$$a$$

is her wealth in deviations from this constraint and

$$\nu$$

is a constant that depends on parameters. Equivalently, her consumption function behaves like her total income plus

$$\sqrt{2 \nu a}$$

⁠. This characterization implies that (1) individuals necessarily hit the borrowing constraint in finite time after a long enough sequence of low-income shocks and (2) we have an intuitive characterization of the speed

$$\nu$$

at which an individual does so as well as her marginal propensity to consume (MPC) out of a windfall income gain. This MPC is higher the lower is the interest rate relative to the rate of time preference, the more willing to intertemporally substitute individuals are, or the higher is the likelihood of getting a high-income draw; it is non-monotone in the income received in low-income states (e.g. unemployment benefits). Understanding the theoretical determinants of MPCs is, of course, important for a large body of applied work.⁶ Second, we derive an analytic solution for the wealth distribution for a special case with two income types. This analytic solution provides a clean characterization of various properties of the wealth distribution, particularly the behaviour of its left and right tails. For example, a direct corollary of individuals hitting the borrowing constraint in finite time is that the wealth distribution features a Dirac point mass at this constraint. Third, we prove existence and uniqueness of a stationary equilibrium for general utility functions and income processes under the intuitive condition that the intertemporal elasticity of substitution (IES)

$$-u'(c)/(u''(c)c)$$

is weakly greater than one for all consumption levels

$$c$$

⁠.⁷ Without a uniqueness result the economy could, in principle, be subject to poverty traps and history dependence.⁸

In addition to these results, which are new also relative to the existing discrete-time literature, we extend some useful existing discrete-time results and concepts to continuous time. First, we adapt a number of results from Aiyagari (1994), e.g., that the wealth distribution has a finite upper bound and that a stationary equilibrium exists. Second, we characterize the saving behaviour of the wealthy and show that, with constant relative risk aversion (CRRA) utility, consumption and saving policy functions become linear for high wealth (Benhabib et al., 2015; Benhabib and Bisin, 2018). Third, we show how to define in continuous time marginal propensities to consume and save over discrete time intervals. This is not obvious and, at the same time, important for bringing the model to the data. Finally, a methodological contribution of our article is to show how to handle borrowing constraints in continuous time: conveniently, the borrowing constraint never binds in the interior of the state space and only shows up in a boundary condition. The consumption first-order condition always holds with equality, thereby sidestepping any complications due to “occasionally binding constraints.” Many of our proofs exploit this fact.

As already mentioned, our second main contribution is the development of a simple, efficient, and portable numerical algorithm for computing a wide class of heterogeneous agent models. The algorithm is based on a finite difference method and applies to the computation of both stationary and time-varying equilibria.⁹ We explain this algorithm in the context of the Aiyagari–Bewley–Huggett model. But the algorithm is, in fact, considerably more general and applies to any heterogeneous agent model with a continuum of atomistic agents (and without aggregate shocks). In Sections 6 and 7, we demonstrate the algorithm’s generality by applying it to other theories that feature non-convexities, a fat-tailed wealth distribution and multiple assets. Codes for these applications (and many more) are available at https://benjaminmoll.com/codes/ in Matlab as well as Python, Julia, and C++.

The first step of the algorithm is to solve the HJB equation for a given time path of prices. The second step is to solve the KF equation for the evolution of the joint distribution of income and wealth. Conveniently, after having solved the HJB equation, one obtains the time path of the distribution essentially “for free,” i.e., with very few lines of code. This is because the KF equation is the “transpose problem” of the HJB equation.¹⁰ The third step is to iterate and repeat the first two steps until an equilibrium time path of prices is found. For the first step, we make use of the theory of “viscosity solutions” to HJB equations (Crandall and Lions, 1983), and the corresponding theory for their numerical solution using finite difference methods (Barles and Souganidis, 1991). While much of our paper can be read without knowledge of viscosity solutions, we provide a brief introduction in Section 6.

Continuous time imparts a number of computational advantages relative to discrete time. As explained in more detail in Section 5.1, these relate to the handling of borrowing constraints, the numerical solution of first-order conditions and the fact that continuous-time problems with discretized state space are, by construction, very “sparse.” These computational advantages are reflected in the algorithm’s efficiency which we showcase in Section 5.6. At the same time, the algorithm is simple. Implementing it requires only some basic knowledge of matrix algebra and access to a software package that can solve sparse linear systems (e.g. Matlab). Finally, the algorithm is portable. For example, it applies without change to problems that involve non-differentiabilities and non-convexities. These are difficult to handle with standard discrete-time methods.¹¹ In contrast, viscosity solutions and finite difference methods are designed to handle non-differentiable and non-convex problems. To illustrate this, we use the same algorithm to compute equilibria of an economy in which the interplay of indivisible housing and mortgages with a down-payment constraint causes a non-convexity which can result in individual poverty traps and multiple stationary distributions, an idea going back to Galor and Zeira (1993) among others.

Besides hopefully being useful in their own right, our paper’s contributions are also the foundation for a number of generalizations that go beyond the setup that we consider in the present paper (as well as the extensions in Section 7). First, building on the finite difference method developed here, Ahn et al. (2017) and Fernández-Villaverde et al. (2019) develop computational methods for solving heterogeneous agent models with aggregate uncertainty in addition to idiosyncratic risk (as in Den Haan, 1997; Krusell and Smith, 1998). Second, Olivi (2018) leverages our continuous-time formulation to obtain powerful comparative statics and sufficient statistics in incomplete markets models—results we build on in our uniqueness proof. Third, Parra-Alvarez et al. (2017) discuss how to identify and estimate continuous-time Aiyagari-Bewley-Huggett models. Fourth, Nuño and Moll (2017) and Nuño and Thomas (2017) devise a method for computing social optima. Fifth, Shaker Akhtekhane (2017) and Moll (2018) extend our computational approach to economies with heterogeneous firms à la Hopenhayn (1992). Sixth, Ruttscheidt (2018), building on Ahn (2017), computes equilibria in economies with a large number of individual state variables (four or more) by marrying our finite difference method with a sparse grid approach (Bungartz and Griebel, 2004; Gerstner and Griebel, 2010). All of these generalizations build on the tools developed in the present article.

A large theoretical and quantitative literature studies environments in which heterogeneous households are subject to uninsurable idiosyncratic shocks. See Heathcote et al. (2009), Guvenen (2011), Quadrini and Ríos-Rull (2015), and Krueger et al. (2015) for recent surveys, and the textbook treatment in Ljungqvist and Sargent (2004). All of these are set in discrete time.

Much fewer papers have studied equilibrium models with heterogeneous households in continuous time.¹² All of these papers make “just the right assumptions” about the environment being studied so that equilibria can be solved explicitly (or at least characterized tightly).¹³ In contrast, our aim is to develop tools for solving and analysing models that do not permit closed-form solutions. Our methods apply as long as the model under consideration can be boiled down to an HJB equation and a KF equation, a feature shared by a wide class of heterogeneous agent models. These two approaches are clearly complementary: on the one hand, having explicit solutions is often extremely valuable for gaining intuition; on the other hand, restricting attention to environments for which these can be found may represent a sort of “analytic straitjacket” for some applications and the availability of more general methods may prove useful in such contexts.

One other paper by Bayer et al. (2019) also studies a continuous-time version of the standard Aiyagari–Bewley–Huggett model. The main differences between their paper and ours are: (1) they analyse a partial equilibrium framework whereas we consider a general equilibrium framework and (2) we develop a numerical algorithm for solving both stationary and time-varying equilibria.¹⁴ Barczyk and Kredler (2014a,b, 2018) study quantitative continuous-time Aiyagari–Bewley–Huggett models with imperfectly altruistic overlapping generations. We instead focus on the simpler workhorse version with infinitely lived individuals, allowing us to prove a number of new theoretical results. Finally, Rocheteau et al. (2015) propose an elegant alternative general equilibrium model with incomplete markets in continuous time. Market incompleteness in their framework stems from lumpy consumption expenditure shocks rather than idiosyncratic income risk. As a result their model features only one individual state variable and many results can be derived in closed form. The tradeoff is that their theory is further from the standard Aiyagari–Bewley–Huggett model that forms the backbone of much of modern macroeconomics.¹⁵

Section 2 lays out our continuous-time version of the workhorse macroeconomic model of income and wealth distribution in the parsimonious form due to Huggett (1993) while Section 3 briefly sketches the general approach of casting heterogeneous agent models as Mean Field Games. Section 4 contains our new theoretical results. Section 5 describes our computational algorithm for both stationary and time-varying equilibria and discusses computational advantages relative to existing discrete-time methods. Section 6 applies the algorithm to a problem with a non-convexity and provides a brief introduction to viscosity solutions. Section 7 discusses a number of generalizations and extensions and Section 8 concludes.

2. The Workhorse Model of Income and Wealth Distribution in Macroeconomics

To explain the logic of our approach in the simplest possible fashion, we present it in a context that should be very familiar to many economists: a general equilibrium model with incomplete markets and uninsured idiosyncratic labour income risk as in Aiyagari (1994), Bewley (1986), and Huggett (1993). We first do this in the context of an economy in which individuals save in unproductive bonds that are in fixed supply as in Huggett (1993). We later consider different ways of closing the model.

2.1. Setup

Individuals. There is a continuum of individuals that are heterogeneous in their wealth

$$a$$

and income

$$y$$

⁠. The state of the economy is the joint distribution of income and wealth. Individuals have standard preferences over utility flows from future consumption

$$c_t$$

discounted at rate

$$\rho \geq 0$$

⁠:

$$$$ \begin{equation}\label{eq:prefs}\mathbb{E}_0 \int_0^\infty e^{-\rho t} u(c_t)dt.\end{equation} $$$$

(1)

The function

$$u$$

is strictly increasing and strictly concave. An individual has an income

$$y_t$$

which is simply an endowment of the economy’s final good. His wealth takes the form of bonds and evolves according to

$$$$ \begin{equation}\label{eq:bc}\dot{a}_t = y_t + r_t a_t - c_t,\end{equation} $$$$

(2)

where

$$r_t$$

is the interest rate. Individuals also face a borrowing limit

$$$$ \begin{equation}\label{eq:borr}a_t \geq \underline{a},\end{equation} $$$$

(3)

where

$$-\infty<\underline{a} < 0$$

⁠.¹⁶ Finally, an individual’s income evolves stochastically over time. In particular, we assume that income follows a two-state Poisson process

$$y_t \in \{y_1,y_2\}$$

⁠, with

$$y_2>y_1$$

⁠. The process jumps from state

$$1$$

to state

$$2$$

with intensity

$$\lambda_1$$

and vice versa with intensity

$$\lambda_2$$

⁠. The two states can be interpreted as employment and unemployment so that

$$\lambda_1$$

is the job-finding rate and

$$\lambda_2$$

the job destruction rate. The two-state Poisson process is chosen for simplicity, and Section 7 extends the setup to more general income processes.

Individuals maximize (1) subject to (2), (3) and the process for

$$y_t$$

⁠, taking as given the evolution of the equilibrium interest rate

$$r_t$$

for

$$t\geq0$$

⁠. For future reference, we denote by

$$g_j(a,t),j=1,2$$

the density of the joint distribution of income

$$y_j$$

and wealth

$$a$$

⁠. As we will see, this distribution will typically feature a Dirac point mass at the borrowing constraint

$$\underline{a}$$

and we therefore write aggregates that integrate over the wealth distribution in terms of the corresponding cumulative distribution function (CDF)

$$G_j(a,t)$$

⁠.¹⁷

Equilibrium. The economy can be closed in a variety of ways. We here present the simplest possible way of doing this following Huggett (1993). We assume that the only price in this economy is the interest rate

$$r_t$$

which is determined by the requirement that, in equilibrium, bonds must be in fixed supply:

$$$$ \begin{equation}\label{eq:eq}\int_{\underline{a}}^\infty a dG_1(a,t) + \int_{\underline{a}}^\infty a dG_2(a,t) = B,\end{equation} $$$$

(4)

where

$$0 \leq B < \infty$$

⁠.

$$B=0$$

means that bonds are in zero net supply. Alternatively,

$$B$$

can be positive. For instance, a government could issue debt and sell it to individuals or there could be saving opportunities abroad. The economy can be closed in a number of alternative ways. For example, in Section 7, we let individuals save in productive capital and the interest rate equals the marginal product of capital of a representative firm as in Aiyagari (1994).

Useful utility functions. We have not imposed any assumptions on the utility function

$$u$$

besides it being strictly increasing and strictly concave. But in later parts of the paper, it will sometimes be instructive to specialize this utility function to either constant relative risk aversion (CRRA) utility

$$$$ \begin{equation}\label{eq:CRRApref} u(c) = \frac{c^{1-\gamma}}{1-\gamma}, \quad \gamma>0, \end{equation} $$$$

(5)

or to exponential utility

$$$$ \begin{equation}\label{eq:exppref} u(c) = -\frac{1}{\theta} e^{-\theta c}, \quad \theta>0. \end{equation} $$$$

(6)

2.2. Stationary equilibrium

Individuals’ consumption–saving decision and the evolution of the joint distribution of their income and wealth can be summarized with two differential equations: a HJB equation and a Kolmogorov Forward (or Fokker–Planck) equation. In a stationary equilibrium, these take the form:¹⁸

$$$$ \begin{align} \label{eq:HJBss}\rho v_j(a) =& \max_c \ u(c) + v_j'(a)(y_j + r a - c) + \lambda_j(v_{-j}(a) - v_j(a)),\\ \end{align} $$$$

(7)

$$$$ \begin{align} \label{eq:KFEss} 0 =& - \frac{d}{da}[s_j(a)g_j(a)] - \lambda_j g_j(a) + \lambda_{-j} g_{-j}(a), \end{align} $$$$

(8)

for

$$j=1,2$$

and where, throughout this article, we adopt the convention that

$$-j=2$$

when

$$j=1$$

and vice versa. The derivations of the HJB equation (7) and the KF equation (8) can be found in Supplementary Appendix B. The function

$$s_j$$

in (8) is the saving policy function, i.e. the optimally chosen drift of wealth

$$$$ \begin{equation}\label{eq:s}s_j(a) = y_j + r a - c_j(a), \quad \mbox{where} \quad c_j(a) = (u')^{-1}(v_j'(a)).\end{equation} $$$$

(9)

The domain of the differential equations (7) and (8) is

$$(\underline{a},\infty),$$

where

$$\underline{a}$$

is the borrowing limit.

The reader may wonder why the borrowing constraint (3) does not feature in the HJB equation (7). The reason is that, in our continuous-time formulation, the borrowing constraint never binds in the interior of the state space, i.e., for

$$a>\underline{a}$$

and as a result an undistorted first-order condition

$$u'(c_j(a))=v_j'(a)$$

holds everywhere.¹⁹ Intuitively, since wealth

$$a$$

is a continuously moving state variable, if it is strictly above the borrowing constraint today, it will still be strictly above the constraint an infinitesimal time interval later. Instead, the borrowing constraint gives rise to a state constraint boundary condition²⁰

$$$$ \begin{equation}\label{eq:stateconstraintss}v_j'(\underline{a}) \geq u'(y_j + r \underline{a}), \quad j=1,2. \end{equation} $$$$

(10)

To see why this is the appropriate boundary condition, note that the first-order condition

$$u'(c_j(a)) = v_j'(a)$$

still holds at

$$a=\underline{a}$$

⁠. The boundary condition (10) therefore implies

$$s_j(\underline{a}) = y_j + r\underline{a} - c_j(\underline{a}) \geq 0$$

⁠, i.e., it ensures that the borrowing constraint is never violated. Supplementary Appendix D derives the state constraint boundary condition (10) more rigorously and makes the connection to a somewhat more general strategy for imposing state constraints used in the mathematics literature, namely to look for a constrained viscosity solution of (7).²¹ The KF equation (8) requires no boundary condition at

$$\underline{a}$$

⁠: the state constraint is satisfied by virtue of

$$s_j$$

being the optimal saving policy function from the HJB equation (7).

Finally, the stationary interest rate

$$r$$

must satisfy the analogue of the market clearing condition (4)

$$$$ \begin{equation}\label{eq:eqss} S(r):= \int_{\underline{a}}^\infty a dG_1(a) + \int_{\underline{a}}^\infty a dG_2(a) = B. \end{equation} $$$$

(11)

The two ordinary differential equations (7) and (8) together with (9), (10), and the equilibrium relationship (11) fully characterize the stationary equilibrium of our economy. In the Mean Field Games (MFG) literature in mathematics this system of coupled HJB and KF equations is called a “backward–forward MFG system,” here in its stationary form.

2.3. Transition dynamics

Many interesting questions require studying transition dynamics, that is the evolution of the economy when the initial distribution of income and wealth does not equal the stationary distribution. The time-dependent analogue of the stationary system (7) to (11) is

$$$$ \begin{align} \label{eq:HJB}\rho v_j(a,t) =& \max_c \ u(c) + \partial_a v_j(a,t)(y_j + r(t) a - c) + \lambda_j(v_{-j}(a,t) - v_j(a,t)) + \partial_t v_j(a,t),\\ \end{align} $$$$

(12)

$$$$ \begin{align} \label{eq:KFE}\partial_t g_j(a,t) =& - \partial_a[s_j(a,t)g_j(a,t)] - \lambda_j g_j(a,t) + \lambda_{-j} g_{-j}(a,t),\\\ \end{align} $$$$

(13)

$$$$ \begin{align} \label{eq:st}s_j(a,t) &= y_j + r(t) a - c_j(a,t), \quad c_j(a,t) = (u')^{-1}(\partial_a v_j(a,t)), \end{align} $$$$

(14)

for

$$j=1,2$$

⁠, and together with the equilibrium condition (4). We here use the short-hand notation

$$\partial_a v = \partial v/\partial a$$

and so on, and as before

$$s_j$$

denotes the optimal saving policy function. The domain of the two partial differential equations (12) and (13) is

$$(\underline{a},\infty) \times \mathbb{R}^+$$

(though more on the time domain momentarily). The function

$$v_j$$

again satisfies a state constraint boundary condition similar to (10)

$$$$ \begin{equation}\label{eq:stateconstraint}\partial_a v_j(\underline{a},t) \geq u'(y_j + r(t) \underline{a}), \quad j=1,2. \end{equation} $$$$

(15)

The density

$$g_j$$

satisfies the initial condition

$$$$ \begin{equation}\label{eq:initial}g_j(a,0)=g_{j,0}(a). \end{equation} $$$$

(16)

The value function satisfies a terminal condition. In principle, the time domain is

$$\mathbb{R}^+$$

but in practice we work with

$$(0,T)$$

for

$$T$$

“large” and impose

$$$$ \begin{equation}\label{eq:terminal}v_j(a,T)=v_{j,\infty}(a),\end{equation} $$$$

(17)

where

$$v_{j,\infty}$$

is the stationary value function that solves the stationary problem (7) to (11).

The two partial differential equations (12) and (13) together with (14), the equilibrium relationship (4) and the boundary conditions (15) to (17) fully characterize the evolution of our economy. This is the time-dependent version of a “backward–forward MFG system.” It has two properties that are worth emphasizing. First, the two equations (12) and (13) are coupled: on one hand, an individual’s consumption-saving decision depends on the evolution of the interest rate which is in turn determined by the evolution of the distribution; on the other hand, the evolution of the distribution depends on individuals’ saving decisions. Second, the two equations run in opposite directions in time: the Kolmogorov Forward equation (13) runs forward (as indicated by its name) and looks backwards—it answers the question “given the wealth distribution today, savings decisions and the random evolution of income, what is the wealth distribution tomorrow?” In contrast, the HJB equation (12) runs backwards and looks forward—it answers the question “given an individual’s valuation of income and wealth tomorrow, how much will she save today and what is the corresponding value function today?”

3. General Heterogeneous Agent Models as MEAN FIELD GAMES

That a heterogeneous agent model boils down to a system of coupled HJB and KF equations is not special to the Aiyagari–Bewley–Huggett model. We here briefly sketch the general approach of casting such models in terms of the mathematical theory of Mean Field Games (MFGs) and provide some references. We also comment on an issue concerning HJB and KF equations, namely what the correct notion of a solution to these partial differential equations (PDEs)is.

3.1. General backward–forward MFG systems

Any heterogeneous agent model with a continuum of atomistic agents (and without aggregate shocks) can be written as a “backward–forward MFG system” of coupled HJB and KF equations. The system from the Aiyagari–Bewley–Huggett model generalizes in two obvious ways. First, the notion of an agent is abstract and also covers firms. Second, general heterogeneous agent models may feature

$$n$$

individual state variables rather than just two. The value function and distribution are then functions of

$$n$$

variables and the backward-forward MFG system is set in

$$n$$

dimensions. Supplementary Appendix C spells out the equations for the abstract,

$$n$$

-dimensional case. To make the mathematics literature accessible to economists, we there also explain the vector calculus notation that is typically used in this literature.

The mathematical theory of MFGs was initiated by Lasry and Lions (2007). Cardaliaguet (2013) and Ryzhik (2018) provide excellent and relatively accessible accounts of its current state. A natural question is whether this literature contains any “off-the-shelf” results on backward–forward MFG systems that apply to the economic models we want to study, say with regard to existence and uniqueness of solutions. As we explain in Supplementary Appendix C.5, the answer is “no” unfortunately: several of the typical features of heterogeneous agent models in economics mean that they are not special cases of the MFGs treated in mathematics.

3.2. Classical versus weak solutions of HJB and KF equations

A classical solution to a PDE or ordinary differential equation (ODE) is a solution that is differentiable as many times as needed to satisfy the corresponding equation. For example, classical solutions to the first-order HJB and KF equations (12) and (13) would need to be once differentiable. Similarly, classical solutions to second-order equations that arise for example if

$$y$$

follows a diffusion process would need to be twice differentiable. In general, we do not expect to find such classical solutions to either HJB or KF equations. For instance, the value function

$$v$$

may have kinks and the distribution

$$g$$

may feature Dirac point masses. Instead, we generally look for certain weak solutions of these equations, that is, solutions that may not be continuously differentiable or even continuous but still satisfy these equations in some sense. As we explain in Supplementary Appendices D and E, the correct notion for a weak solution of the HJB equation is a viscosity solution and that of the KF equation is a measure-valued solution.²² See Evans (2010, Section 1.3) and Tao (2008) for illuminating discussions on the role of weak solutions in the study of partial differential equations more generally.

Most of our article employs classical methods and the preceding paragraph equips the reader sufficiently well for those parts that do not. An exception is the model with a non-convexity due to indivisible housing in Section 6. We there discuss in more detail the usefulness of viscosity solutions.

4. Theoretical Results for the Aiyagari–Bewley–Huggett Model: Consumption, Saving, and Inequality

This Section presents theoretical results about our continuous-time version of the Aiyagari–Bewley–Huggett model, including the three new results emphasized in the introduction. Sections 4.1 to 4.5 analyse the HJB and KF equations (7) and (8) in partial equilibrium, i.e., taking as given a fixed interest rate

$$r$$

(assumed to be less than

$$\rho$$

which will be the equilibrium outcome). Section 4.6 then imposes market clearing (11) and considers the stationary equilibrium, particularly its existence and uniqueness. While we focus on the case with two income states for pedagogical reasons, all results except the closed-form wealth distribution in Proposition 3 generalize to alternative income processes (see Section 7).

4.1. An Euler equation

Our first few theoretical results concern the consumption and saving behaviour of individuals. Our characterization of individual behaviour uses the following Lemma.

Lemma 1

The consumption and saving policy functions
$$c_j(a)$$
and
$$s_j(a)$$
for
$$j=1,2$$
corresponding to the HJB equation (7) satisfy

$$$$ \begin{equation} \begin{split} \label{eq:Euler}(\rho -r) u'(c_j(a)) &= u''(c_j(a))c_j'(a)s_j(a) + \lambda_j(u'(c_{-j}(a)) - u'(c_j(a))),\\ s_j(a)&=y_j + ra - c_j(a). \end{split} \end{equation} $$$$

(18)

Proof

Differentiate the HJB equation (7) with respect to

$$a$$

(envelope condition) and use that

$$v_j'(a)=u'(c_j(a))$$

and hence

$$v_j''(a)=u''(c_j(a))c_j'(a)$$

⁠.

The differential equation (18) is an Euler equation. The right-hand side is simply the expected change of individual marginal utility of consumption

$$\mathbb{E}_t[d u'(c_j(a_t))]/dt$$

⁠.²³ Therefore (18) is equivalent to

$$$$ \frac{\mathbb{E}_t[d u'(c_j(a_t))]}{u'(c_j(a_t))} = (\rho - r)dt. $$$$

4.2. Consumption and saving behaviour of the poor

Our first main result is obtained by analysing the Euler equation (18) close to the borrowing constraint. The interesting case is when the behaviour at the constraint differs qualitatively from that of rich individuals. Whether this is the case depends crucially on two factors: the tightness of the borrowing constraint

$$\underline{a}$$

⁠, and the properties of the utility function at low levels of consumption. To focus on the interesting case, we make the following assumption.

Assumption 1

The coefficient of absolute risk aversion
$$R(c) := -u''(c)/u'(c)$$
remains finite at the borrowing limit

$$$$ \underline{R} := - \frac{u''(y_1 + r\underline{a})}{u'(y_1 + r\underline{a})}<\infty. $$$$

The next proposition shows that the borrowing constraint “matters” if this assumption holds. Standard utility functions satisfy

$$R(c)<\infty$$

for

$$c>0$$

but

$$\lim_{c \rightarrow 0} R(c)=\infty$$

⁠. For example, with CRRA utility (5),

$$R(c)=\gamma/c$$

⁠. Therefore, for such standard utility functions Assumption 1 is equivalent to

$$\underline{a}>-y_1/r$$

⁠, i.e., the borrowing constraint matters if it is tighter than the “natural borrowing constraint.” However, Assumption 1 is considerably weaker than this. In particular, if the utility function is such that absolute risk aversion remains bounded as consumption goes to zero,

$$R(0)<\infty$$

⁠, Assumption 1 holds and hence the constraint matters even with the natural borrowing constraint

$$\underline{a}=-y_1/r$$

⁠. An example is exponential utility (6) for which

$$\underline{R}=\theta<\infty$$

regardless of the tightness of the borrowing constraint. Summarizing, Assumption 1 says that either the borrowing constraint is tighter than the natural borrowing constraint or the coefficient of absolute risk aversion is bounded as consumption approaches zero (or both).²⁴ For completeness, the case in which Assumption 1 is violated is covered in Proposition 1’ in the Supplementary Appendix.

In what follows as well as elsewhere in the article, we use the following asymptotic notation: for any two functions

$$f$$

and

$$g$$

⁠, “

$$f(a)\sim g(a)$$

$$a \rightarrow \underline{a}$$

” is short-hand notation for

$$\lim_{a \rightarrow \underline{a}} f(a)/g(a)=1$$

⁠, i.e.,

$$f$$

“behaves like”

$$g$$

for

$$a$$

close to

$$\underline{a}$$

⁠.

Proposition 1 (MPCs and Saving at Borrowing Constraint)

Let

$$s_1(a)$$

and
$$c_1(a)$$
be the optimal saving and consumption policy functions of the low-income type. If
$$r<\rho$$
and Assumption 1 holds, then:

$$s_1(\underline{a})=0$$
but
$$s_1(a)<0$$
all
$$a>\underline{a}$$
⁠. That is, only individuals exactly at
$$\underline{a}$$
are constrained, whereas those with wealth
$$a>\underline{a}$$
are unconstrained and decumulate assets.
as
$$a\rightarrow \underline{a}$$
⁠, the saving and consumption policy functions of the low-income type and the corresponding instantaneous marginal propensity to consume satisfy
$$$$ \begin{align}\label{eq:approx} s_1(a) &\sim - \sqrt{2\nu_1} \sqrt{a - \underline{a}},\\ \notag c_1(a) &\sim y_1 + r a + \sqrt{2\nu_1}\sqrt{a - \underline{a}},\\ \end{align} $$$$
(19)

$$$$ \begin{align} \label{eq:MPC_inst} c_1'(a) &\sim r + \sqrt{\frac{\nu_1}{2(a - \underline{a})}},\\ \end{align} $$$$
(20)

$$$$ \begin{align} \label{eq:nu} \begin{split} \nu_1 :&= \frac{(\rho - r)u'(\underline{c}_{1}) + \lambda_1 (u'(\underline{c}_{1}) - u'(\underline{c}_{2}))}{-u''(\underline{c}_{1})}\\ &\approx (\rho-r)\text{IES}(\underline{c}_1)\underline{c}_1 + \lambda_1(\underline{c}_2 - \underline{c}_1), \end{split} \end{align} $$$$
(21)
where
$$\underline{c}_j = c_j(\underline{a}),j=1,2$$
and
$$\textrm{IES}(c):=-u'(c)/(u''(c)c)$$
⁠.²⁵ The derivatives of
$$c_1$$
and
$$s_1$$
are unbounded at the borrowing constraint,
$$c_1'(a) \rightarrow \infty$$
and
$$s_1'(a) \rightarrow -\infty$$
as
$$a\rightarrow \underline{a}$$
⁠.

The proof of the proposition, like those of all others, is in the Supplementary Appendix. The proof of the first part follows straight from the state constraint boundary condition (10). The second part of the proof follows from characterizing the limiting behaviour of the squared saving policy function

$$(s_1(a))^2$$

as wealth

$$a$$

approaches the borrowing constraint

$$\underline{a}$$

– hence the square root.

The consumption and saving behaviour in the Proposition is illustrated in Figure 1.

$Consumption and saving behaviour with $$r<\rho$$$

Figure 1

Consumption and saving behaviour with

$$r<\rho$$

Open in new tab Download slide

Importantly, the derivatives of type

$$1$$

’s consumption and saving policy functions become unbounded at the borrowing constraint. This unbounded derivative has an important implication, namely that individuals hit the borrowing constraint in finite time.

Corollary 1 (Hit Constraint in Finite Time)

If
$$r<\rho$$
and Assumption 1 holds, then the wealth of an individual with initial wealth
$$a_0$$
and successive low-income draws
$$y_1$$
hits the borrowing constraint at a finite time
$$T$$
and converges toward it at speed governed by
$$\nu_1$$
⁠:

$$$$ \begin{equation}\label{eq:a(t)} a(t) - \underline{a} \sim \frac{\nu_1}{2}\left(T-t\right)^2, \quad T:=\sqrt{\frac{2(a_0 - \underline{a})}{\nu_1}}, \quad 0 \leq t\leq T. \end{equation} $$$$

(22)

The result that the borrowing constraint is reached in finite time bears some similarity to optimal stopping time problems (see e.g. Stokey, 2009). Just like in stopping time problems, continuous time avoids a type of integer problem arising in discrete time: the borrowing constraint would be reached after a non-integer time period, but discrete time forces this to occur after an integer number of periods.²⁶

Proposition 1 features an intuitive formula (21) for the speed at which individuals hit the borrowing constraint,

$$\nu_1$$

⁠. In Section 4.4, we show that

$$\nu_1$$

is also the key quantity determining the marginal propensity to consume (MPC) out of a windfall income gain. We therefore postpone the discussion of formula (21) until that section.

Intuition for Proposition 1 and Corollary 1: two useful special cases. To understand the intuition for the square root in Proposition 1, the saving behaviour in Corollary 1 and the role of Assumption 1, we now consider two special cases for which analytic solutions are available. Both abstract from income uncertainty which is inessential to this point.²⁷

In the first special case, an individual has exponential utility (6), receives a deterministic income stream

$$y>0$$

⁠, faces a strict no-borrowing constraint

$$a\geq 0$$

and starts with some initial wealth

$$a_0>0$$

⁠. The corresponding Euler equation and budget constraint are

$$$$ \dot{c} = \frac{1}{\theta}(r-\rho), \qquad \dot{a} = y + ra - c. $$$$

Conjecture that at some time

$$T>0$$

⁠, individuals hit the borrowing constraint, i.e.

$$a(T)=0$$

and hence

$$c(T)=y$$

⁠. From the Euler equation, consumption for

$$t \leq T$$

$$$$ \begin{equation}\label{eq:c_xavier} c(t) = y + \nu(T-t), \quad \nu := \frac{\rho-r}{\theta} > 0 \end{equation} $$$$

(23)

The case

$$r=0$$

contains all the intuition.²⁸ Substituting into the budget constraint yields

$$\dot{a}(t) = -\nu(T-t)$$

with solution

$$$$ \begin{equation}\label{eq:a_xavier} a(t) = \frac{\nu}{2}(T-t)^2 \end{equation} $$$$

(24)

for

$$t \leq T$$

and where the constant of integration is zero because

$$a(T)=0$$

⁠. Since

$$a(0)=a_0$$

⁠, the hitting time is given by

$$T = \sqrt{2a_0/\nu}$$

⁠. These are the same expression as in Corollary 1. Figure 2 plots the time paths of consumption, saving and wealth and shows that consumption declines linearly toward

$$c(T)=y$$

while wealth declines quadratically towards

$$a(T)=0$$

⁠.

Figure 2

First special case in which borrowing constraint binds in finite time

Open in new tab Download slide

To understand the square root in the saving and consumption policy functions in Proposition 1, consider an individual at

$$t=0$$

with some initial wealth

$$a_0$$

⁠. From (23) and (24), we have

$$c(0) = y + \nu T$$

and

$$\dot{a}(0) = -\nu T$$

with

$$T = \sqrt{2a_0/\nu}$$

⁠. Writing consumption and saving in terms of the state variable

$$a$$

rather than calendar time, we have

$$c(a) = y + \sqrt{2 \nu a}$$

and

$$s(a) = - \sqrt{2 \nu a}$$

which are the square-root expressions from Proposition 1. This simple derivation also shows why the consumption policy function is concave in wealth

$$a$$

(Figure 1). As the individual approaches the borrowing constraint, both her consumption and wealth decline. If both consumption and wealth declined at the same speed, then consumption would be linear in wealth. Instead, wealth declines more rapidly than consumption – quadratically rather than linearly—see Figures 2(a) and (c)—and therefore consumption is strictly concave in wealth, and more so the higher is the speed

$$\nu$$

⁠.

In contrast, consider a second special case which is identical except that

$$y=0$$

and that the individual has CRRA utility (5). The Euler equation and budget constraint are then

$$$$ \frac{\dot{c}}{c} = \frac{1}{\gamma}(r-\rho), \qquad \dot{a} = ra - c. $$$$

It is easy to show that

$$\dot{a}(t) = -\eta a(t), c(t) = \left(r+\eta\right)a(t)$$

where

$$\eta:=\frac{\rho-r}{\gamma}$$

and

$$a(t) = a_0 e^{-\eta t}$$

⁠. The situation is depicted in Figure 3. As wealth decumulates towards the borrowing constraint, the rate of decumulation slows down more and more and individuals never actually hit the borrowing constraint. This is an immediate consequence of a linear saving policy function

$$s(a)=-\eta a$$

⁠. Turning this logic around, consumption is linear in wealth because both consumption and wealth decline toward the borrowing constraint at the same speed—see Figure 3(a) and (c)—rather than wealth declining faster as in the exponential case.

Figure 3

Second special case in which borrowing constraint never binds

Open in new tab Download slide

4.3. Consumption and saving behaviour of the wealthy

Proposition 1 characterizes consumption and saving behaviour close to the borrowing constraint. The following Proposition 2 characterizes this behaviour for large wealth levels. This will be useful below, when we characterize the upper tail of the wealth distribution.

Proposition 2 (Consumption and Saving Behaviour of the Wealthy)

Assume that

$$r<\rho$$

and that relative risk aversion
$$-cu''(c)/u'(c)$$
is bounded above for all
$$c$$
⁠.

Then there exists
$$a_{\max}<\infty$$
such that
$$s_j(a)<0$$
for all
$$a > a_{\max},j=1,2$$
⁠, and
$$s_2(a)\sim \zeta_2(a_{\max}-a)$$
as
$$a \rightarrow a_{\max}$$
for some constant
$$\zeta_2$$
.
In the special case of CRRA utility (5) individual policy functions are asymptotically linear in
$$a$$
⁠:
$$$$ \begin{equation} \label{eq:CRRAlim}s_j(a) \sim \frac{r- \rho}{\gamma} a, \quad c_j(a) \sim \frac{\rho - (1-\gamma) r}{\gamma} a \quad \mbox{as} \ a \rightarrow \infty \end{equation} $$$$
(25)

The first part of the Proposition is the analogue of Proposition 4 in Aiyagari (1993). The condition that

$$-cu''(c)/u'(c)$$

is bounded above for all

$$c$$

for example rules out exponential utility (6) in which case

$$\gamma(c)=\theta c$$

⁠.

The second part of the Proposition extends to continuous time a result by Benhabib et al. (2015) who have shown that, with CRRA utility, consumption, and saving policy functions are asymptotically linear for large wealth.²⁹ To understand the intuition for this linearity, consider a special case without income risk: individuals have CRRA utility, (5), deterministic labour income

$$y>0$$

⁠, and face the natural borrowing constraint

$$\underline{a}=-y/r$$

⁠. Consumption and saving policy functions then have a closed-form solution given by

$$$$ \begin{equation}\label{eq:CRRAdeterministic} s(a) = \frac{r-\rho}{\gamma}\left(a + \frac{y}{r}\right), \qquad c(a) = \frac{\rho - (1-\gamma) r}{\gamma} \left(a+\frac{y}{r}\right). \end{equation} $$$$

(26)

$$a\rightarrow \infty$$

⁠,

$$y/r$$

becomes small relative to

$$a$$

and the policy functions indeed satisfy (25).

The asymptotic linearity of consumption and saving policy functions with CRRA utility has played a key role in the literature. For instance, Krusell and Smith (1998) argue that this linearity explains their finding that the business cycle properties of a baseline heterogeneous agent model are virtually indistinguishable from its representative agent counterpart. Future studies may want to gauge the robustness of this result to relaxing the CRRA assumption.

4.4. Marginal propensities to consume and save

We now characterize marginal propensities to consume and save, defined as the changes in consumption and saving in response to a windfall increase in available funds

$$a$$

⁠. Propositions 1 and 2 characterize the slope of the consumption function

$$c_j'(a)$$

or, equivalently, the instantaneous MPC which captures the consumption gain (per time unit) after such a windfall over an infinitesimally small time interval. This is an interesting object but it does not correspond to what is measured in the data, namely the fraction of income consumed out of a windfall income gain over a discrete time interval. We here show how to characterize this more empirically relevant object.

Definition 1

The Marginal Propensity to Consume (MPC) over a period
$$\tau$$
is given by

$$$$ \begin{equation} \label{eq:MPC_tau}\textrm{MPC}_{j,\tau}(a) = C_{j,\tau}'(a), \quad \mbox{where} \quad C_{j,\tau}(a) = \mathbb{E}\left[\int_0^\tau c_{j}(a_t)dt \Big\vert a_0 = a,y_0 = y_j \right]\!. \end{equation} $$$$

(27)

Similarly, the Marginal Propensity to Save (MPS) over a period
$$\tau$$
is given by

$$$$ \begin{equation} \label{eq:MPS_tau} \textrm{MPS}_{j,\tau}(a) = S_{j,\tau}'(a), \quad \mbox{where} \quad S_{j,\tau}(a) = \mathbb{E}\left[a_\tau \Big\vert a_0 = a,y_0 = y_j \right]\!. \end{equation} $$$$

(28)

To get a feel for the behaviour of these objects and to see how they differ from their instantaneous counterparts

$$c_j'(a)$$

and

$$s_j'(a)$$

⁠, it is instructive to consider a time interval

$$\tau$$

that is sufficiently small so that individuals do not switch income state. Expected saving over a period

$$\tau$$

for the low-income type

$$S_{1,\tau}(a)$$

is then simply given by

$$a_\tau$$

starting from

$$a_0=a$$

⁠. But from Corollary 1 we know that

$$a_\tau - \underline{a} \sim \frac{\nu_1}{2}\left((T-\tau)^{+}\right)^2$$

with

$$T=\sqrt{2(a_0-\underline{a})/\nu_1}$$

⁠, i.e.

$$$$ \begin{equation}\label{eq:S_tau_det} S_{1,\tau}(a) \sim \frac{\nu_1}{2}\left(\left(\sqrt{ 2 (a - \underline{a})/\nu_1} - \tau\right)^{+}\right)^2 + \underline{a}. \end{equation} $$$$

(29)

Differentiating this expression and using the budget constraint, we get the following result.

Corollary 2

Assume
$$\tau$$
is sufficiently small that individuals with current income draw
$$y_1$$
do not receive the high-income draw
$$y_2$$
⁠,³⁰ that
$$r<\rho$$
⁠, and that Assumption 1 holds. Then

$$$$ \begin{align} \label{eq:MPC_tau_det}\textrm{MPC}_{1,\tau}(a) &\sim \min \left\{\tau \sqrt{\frac{\nu_1}{2(a-\underline{a})}},1\right\} + \tau r = \min\{\tau c_1'(a),1 + \tau r\}, \quad \mbox{as} \ a \rightarrow \underline{a}, \end{align} $$$$

(30)

where
$$c_1'(a)$$
is the instantaneous MPC from Proposition 1. Alternatively, (30) holds with equality in the special case with exponential utility, deterministic income and
$$\underline{a}=r=0$$
.

We make three observations. First, in contrast to the instantaneous MPC

$$c'_1(a)$$

which becomes unbounded as

$$a \rightarrow \underline{a}$$

⁠, the MPC over a time period

$$\tau$$

in (30) is bounded between zero and

$$1+\tau r$$

⁠. Second, for

$$a>\underline{a}$$

and

$$\tau$$

small enough, the marginal propensity to consume (30) is strictly decreasing in wealth

$$a$$

(i.e.

$$C_{1,\tau}(a)$$

⁠, is strictly concave in

$$a$$

⁠). Third, the key quantity determining the size of the MPC is

$$\nu_1$$

⁠, the speed at which individuals hit the borrowing constraint. We have already discussed the intuition for the last two properties: because wealth declines toward the borrowing constraint faster than consumption, the mapping from wealth to consumption is concave; the faster wealth declines, the more concave this mapping.

In contrast, consider the second special case of Section 4.2 in which individuals never hit the borrowing constraint. From

$$a(t) = a_0 e^{-\eta t}$$

with

$$\eta=:(\rho-r)/\gamma$$

⁠, saving over a period

$$\tau$$

is given by

$$S_\tau(a) = a e^{-\eta \tau}$$

⁠. Since both consumption and saving over a period

$$\tau$$

are linear in wealth

$$a$$

⁠, the marginal propensities to save and consume are independent of wealth:³¹

$$$$ \textrm{MPS}_\tau(a) = e^{-\eta \tau} \approx 1-\eta \tau, \quad \textrm{MPC}_\tau(a) = 1 - e^{-\eta \tau} + \tau r \approx \tau(\eta + r), \quad \eta:=\frac{\rho-r}{\gamma}. $$$$

Summarizing, when people hit the borrowing constraint in finite time, MPCs depend on wealth and, in particular, are higher for poorer people.

When individuals experience new income draws within the time interval

$$\tau$$

⁠, it is no longer possible to characterize the MPC and MPS as tightly as in Corollary 2 because we lack a characterization of

$$c_2(a)$$

in the vicinity of the borrowing constraint. However, the following Lemma shows how to easily compute the MPC numerically. The key idea is that

$$C_{j,\tau}(a)$$

defined in (27) is a conditional expectation that satisfies the Feynman-Kac formula which establishes a link between conditional expectations of stochastic processes and solutions to partial differential equations. Given

$$C_{j,\tau}(a)$$

⁠, we then compute

$$\textrm{MPC}_{j,\tau}(a)=C_{j,\tau}'(a)$$

⁠.

Lemma 2 (Computation of MPCs using Feynman–Kac formula)

The conditional expectation
$$C_{j,\tau}(a)$$
defined in (27) and therefore the MPC can be computed as
$$C_{j,\tau}(a) = \Gamma_j(a,0)$$
where
$$\Gamma_j(a,t)$$
satisfies the system of two PDEs

$$$$ 0 = c_j(a) + \partial_a \Gamma_j(a,t)s_j(a) + \lambda_j(\Gamma_{-j}(a,t) - \Gamma_j(a,t)) + \partial_t \Gamma_j(a,t), \quad j=1,2 $$$$

on
$$(\underline{a},\infty) \times (0,\tau)$$
⁠, with terminal condition
$$\Gamma_j(a,\tau)=0$$
for all
$$a$$
⁠.

Proof

This follows from a direct application of the Feynman–Kac formula.

Figure 4(a) plots the MPC computed in this way for the two income types and assuming that individuals have CRRA utility (5). For comparison, Figure 4(b) plots the “instantaneous MPC,” i.e. the slope of the consumption function.

Figure 4

MPCs across the wealth distribution

Open in new tab Download slide

As expected, the former is a smoother version of the latter and, in contrast to the latter, does not exceed

$$1+\tau r$$

⁠. As wealth

$$a \rightarrow \infty$$

the borrowing constraint becomes irrelevant and the slope of the consumption function converges to

$$\eta+r$$

with

$$\eta=(\rho-r)/\gamma$$

and therefore the MPC to

$$\tau(\eta+r)$$

⁠.

As an aside, in some applications a slightly altered version of the MPCs in Definition 1 is easier to map to the data. Empirical studies do not typically estimate MPCs out of an infinitesimal increase in resources. Instead, they estimate the increase in consumption in response to a discrete increase, say by

$$\$500$$

⁠. To this end, define

$$\textrm{MPC}_{j,\tau}^x(a) := (C_{j,\tau}(a+x) - C_{j,\tau}(a))/x.$$

This is the MPC out of

$$x$$

dollars over a period

$$\tau$$

⁠, i.e., a discrete counterpart to the MPC in (27). Kaplan et al. (2018) compute such discrete MPCs and compare them to various empirical studies such as Broda and Parker (2014), Misra and Surico (2014), Blundell et al. (2016), and Fagereng et al. (2016).

Using the analytic expression for

$$\nu_1$$

to better understand MPCs. As part of Proposition 1 we obtained an analytic expression (21) for

$$\nu_1$$

⁠, the speed of hitting the borrowing constraint. As just discussed,

$$\nu_1$$

is also the key quantity governing the size of MPCs. We can therefore use the formula (21) to examine how MPCs depend on various model parameters and to shed some light on various numerical results that may seem counterintuitive at first. For instance, consider the dependence of the low-income type’s

$$\textrm{MPC}_{1,\tau}(a)$$

on the low-income realization

$$y_1$$

⁠. This low-income realization may, for example, represent the size of unemployment benefits.

Figure 5(a) graphs this relationship for

$$y_1$$

ranging between

$$0$$

and

$$y_2=0.2$$

separately for various percentiles of the wealth distribution and assuming that individuals have CRRA utility (5). Perhaps surprisingly, the MPC is a hump-shaped function of the low-income realization

$$y_1$$

⁠. But formula (21) easily resolves the apparent mystery. With CRRA utility (5) so that the IES is constant

$$$$ \nu_1 \approx (\rho-r)\frac{\underline{c}_1}{\gamma} + \lambda_1(\underline{c}_2 - \underline{c}_1), $$$$

where

$$\underline{c}_1=y_1+r\underline{a}$$

⁠. An increase in

$$y_1$$

has two offsetting effects. The intuitive part is that as

$$y_1$$

increases, individuals are better insured against idiosyncratic income risk and therefore have a low MPC (as in models without risk and borrowing constraints). In the formula, as

$$y_1$$

increases toward

$$y_2$$

⁠,

$$\lambda_1(\underline{c}_2 - \underline{c}_1)$$

converges to zero and this results in a lower

$$\nu_1$$

⁠. But, there is an offsetting effect captured by the term

$$(\rho - r)\underline{c}_1/\gamma$$

⁠: if consumption conditional on hitting the constraint

$$\underline{c}_1$$

is high, individuals do not mind hitting the constraint as much. Hence they converge to it faster or, equivalently, have a higher MPC.

Figure 5

Dependence of MPCs on parameters

Open in new tab Download slide

Figure 5(b) instead graphs the dependence of the low-income type’s MPC on the realization of the high income

$$y_2$$

⁠. The MPC is increasing in

$$y_2$$

and the intuition can again be seen from our formula for

$$\nu_1$$

which shows that the MPC is higher the larger is the consumption gain from getting a high-income draw

$$\lambda_1(\underline{c}_2 - \underline{c}_1)$$

⁠. Other comparative statics are as follows: individuals have higher MPCs, the lower is the interest rate

$$r$$

relative to the rate of time preference

$$\rho$$

⁠, and the higher is the likelihood

$$\lambda_1$$

of getting a high-income draw (so that getting stuck at the constraint is less likely). Similarly, MPCs tend to be higher the higher is the IES and the tighter is the borrowing constraint, i.e. the closer to zero is

$$\underline{a}$$

⁠.

4.5. The stationary wealth distribution

We now present the paper’s second main theoretical result: an analytic solution to the Kolmogorov Forward equation characterizing the stationary distribution with two income types (8) for given individual saving policy functions. This analytic solution yields a number of insights about properties of the stationary wealth distribution, particularly at the borrowing constraint and in the right tail.

The derivation of this analytic solution is constructive and we therefore present it in the main text. Summing the KF equation (8) for the two income types, we have

$$\frac{d}{da}[s_1(a)g_1(a) + s_2(a)g_2(a)] = 0$$

for all

$$a$$

⁠, which implies that

$$s_1(a)g_1(a)+s_2(a)g_2(a)$$

equals a constant. Because any stationary distribution must be bounded, we must then have

$$s_1(a)g_1(a) + s_2(a)g_2(a) = 0$$

for all

$$a$$

⁠.³² Substituting into (8) and rearranging, we have

$$$$ \begin{align} \label{eq:gODE}g_j'(a) &= -\left(\frac{s_j'(a)}{s_j(a)} + \frac{\lambda_j}{s_j(a)} + \frac{\lambda_{-j}}{s_{-j}(a)}\right)g_j(a), \quad j=1,2. \end{align} $$$$

(31)

Importantly, (31) are two independent ODEs for

$$g_1$$

and

$$g_2$$

rather than the coupled system of two ODEs (8) we started out with. Together with two boundary conditions they can be solved separately. To obtain these boundary conditions, we simply impose that the densities integrate to the stationary mass of individuals with the respective income types:

$$$$ \begin{equation}\label{eq:integrate}m_1 + \int_{\underline{a}}^\infty g_1(a)da = \frac{\lambda_2}{\lambda_1+\lambda_2}, \quad \quad m_2 + \int_{\underline{a}}^\infty g_2(a)da = \frac{\lambda_1}{\lambda_1+\lambda_2},\end{equation} $$$$

(32)

where

$$m_1 = G_1(\underline{a})$$

and

$$m_2 = G_2(\underline{a})$$

are potential Dirac masses at the borrowing constraint. Solving the two ODEs (31) analytically and using our characterization of the optimal saving policy functions from Section 4.2 we obtain our second main theoretical result.

Proposition 3 (Stationary Wealth Distribution with Two Income Types)

If
$$r<\rho$$
⁠, relative risk aversion
$$-cu''(c)/u'(c)$$
is bounded above for all
$$c$$
⁠, and Assumption 1 holds, then there exists a unique stationary distribution given by

$$$$ \begin{equation}\label{eq:gSol} g_j(a) = \frac{\kappa_j}{s_j(a)}\exp\left(-\int_{\underline{a}}^a \left(\frac{\lambda_1}{s_1(x)} + \frac{\lambda_2}{s_2(x)}\right) dx\right),\quad j=1,2 \end{equation} $$$$

(33)

for some constants of integration
$$\kappa_1<0$$
and
$$\kappa_2>0$$
that satisfy
$$\kappa_1+\kappa_2=0$$
and are uniquely pinned down by (32). The stationary wealth distribution has the following properties:

(Close to the borrowing constraint) The stationary distribution of low-income types has a Dirac point mass at the borrowing constraint
$$\underline{a}$$
⁠, i.e., its CDF satisfies
$$G_1(\underline{a})=m_1>0$$
⁠. The Dirac point mass
$$m_1$$
can be found from the constants of integration
$$\kappa_1,\kappa_2$$
and is expressed in terms
$$\lambda_1,\lambda_2,s_1,s_2$$
in Supplementary Appendix equation (A.71). The CDF further satisfies
$$$$ \begin{equation}\label{eq:CDF1} G_1(a) \sim m_1 \exp\left(\lambda_1 \sqrt{2(a-\underline{a})/\nu_1} \right) \quad \mbox{as} \ a \downarrow \underline{a} \end{equation} $$$$
(34)
The stationary distribution of high-income types does not have a Dirac point mass at
$$\underline{a}$$
⁠, i.e., its CDF satisfies
$$G_2(\underline{a})= m_2 =0$$
⁠, and its density is in fact finite,
$$g_2(\underline{a})<\infty$$
⁠.
(In the right tail) The support of the stationary wealth distribution is bounded above at some
$$a_{\max}<\infty$$
defined in Proposition 2. It does not have a Dirac point mass at
$$a_{\max}$$
⁠.
(Smoothness) In contrast to the analogous discrete-time economy, the density of wealth is continuous and differentiable for all
$$a>\underline{a}$$
⁠, i.e., everywhere except at the borrowing constraint.

Corollary 3 in the Supplementary Appendix lists some additional but less central properties of the wealth distribution that follow from (33).

Part 1 of the Proposition follows immediately from Proposition 1 and that individuals hit the borrowing constraint in finite time if Assumption 1 holds. In this case, (1)

$$g_1$$

in (33) explodes as

$$a\rightarrow \underline{a}$$

and (2) there is a Dirac mass at

$$\underline{a}$$

⁠,

$$G_1(\underline{a})=m_1>0$$

⁠. This is illustrated in Figure 6(b). In particular, note the spike in the density

$$g_1(a)$$

$$a=\underline{a}$$

⁠. In contrast, if Assumption 1 is violated, then there is no Dirac mass. An alternative heuristic derivation of (34) provides some intuition. First, (8) implies that

$$G_1$$

satisfies

$$0 = -s_1(a)G_1'(a) - \lambda_1 G_1(a) + \lambda_2 G_2(a)$$

⁠. As

$$a \downarrow \underline{a}$$

⁠,

$$G_2(a)\rightarrow 0$$

and hence

$$G_1$$

satisfies

$$G_1'(a)/G_1(a) \sim -\lambda_1/s_1(a)$$

with solution

$$$$ G_1(a) \sim m_1 \exp\left(-\int_{\underline{a}}^a \frac{\lambda_1}{s_1(x)}dx \right), $$$$

for a constant of integration

$$m_1>0$$

⁠. Substituting in

$$s_1(a)\sim -\sqrt{2\nu_1}\sqrt{a-\underline{a}}$$

from Proposition 1 yields (34). In Supplementary Appendix A.6, we discuss in more detail the role of Assumption 1 using our two simple special cases from Section 4.2.

$Saving behaviour and stationary wealth distribution with $$r<\rho$$$

Figure 6

Saving behaviour and stationary wealth distribution with

$$r<\rho$$

Open in new tab Download slide

Part 2 of the Proposition states that the stationary wealth distribution in our economy is bounded above. Like discrete-time versions of Aiyagari–Bewley–Huggett economies with idiosyncratic labour income risk only, our model therefore has difficulties explaining the high observed wealth concentration in developed economies like the U.S. In particular, empirical wealth distributions seem to feature fat Pareto tails. Section 7 extends the model to feature such a fat-tailed stationary distribution by introducing a second, risky asset.

Part 3, which can also be seen in Figure 1, highlights an important difference between our continuous-time formulation and the traditional discrete-time one: except for the Dirac mass exactly at the borrowing constraint

$$\underline{a}$$

⁠, the wealth distribution is smooth for all

$$a>\underline{a}$$

⁠. This is true even though income follows a process with discrete states (a two-state Poisson process). In contrast, in discrete-time Aiyagari–Bewley–Huggett models with discrete-state income processes, this distribution features “spikes” on the interior of the state space.³³

4.6. Stationary equilibrium: existence and uniqueness

We construct stationary equilibria along the same lines as in Aiyagari (1994). That is, we fix an interest rate

$$r<\rho$$

⁠, solve the individual optimization problem (7), find the corresponding stationary distribution from (8), and then find the interest rate

$$r$$

that satisfies the market clearing condition (11), i.e.,

$$S(r)=B$$

⁠. While we continue to focus on the case of two income types for the sake of continuity, all results in this section generalize to any stationary Markovian process for income

$$y$$

⁠, e.g., continuous diffusion or jump-diffusion processes.

Figure 7 illustrates the typical effect of an increase in

$$r$$

on the solutions to the HJB equation (7) and the KF equation (8).

Figure 7

Effect of an increase in

$$r$$

on saving behaviour and stationary distribution.

Open in new tab Download slide

An increase in

$$r$$

from

$$r_L$$

$$r_H>r_L$$

leads to an increase in individual saving at most wealth levels and the stationary distribution shifts to the right. Aggregate saving

$$S(r)$$

as a function of the interest rate

$$r$$

typically looks like in Figure 8, i.e., it is increasing.

Figure 8

Equilibrium in the bond market

Open in new tab Download slide

A stationary equilibrium is then an interest rate

$$r$$

such that

$$S(r)=B$$

⁠. But, we have so far not proven that

$$S(r)$$

is increasing or that it intersects

$$B$$

and hence there may, in principle, be no or multiple equilibria. Existence of a stationary equilibrium can be proved with a graphical argument due to Aiyagari (1994) that is also the foundation for a number of existence results in the literature (e.g. Açıkgöz, 2018).

Proposition 4 (Existence of Stationary Equilibrium)

If relative risk aversion

$$-cu''(c)/u'(c)$$

is bounded above for all
$$c$$
and Assumption 1, then there exists a stationary equilibrium.

The logic behind the proof is simple. One can show that the function

$$S(r)$$

defined in (11) is continuous. To guarantee that there is at least one

$$r$$

such that

$$S(r)=B$$

⁠, it then suffices to show that

$$$$ \lim_{r \uparrow \rho} S(r) = \infty, \qquad \lim_{r \downarrow -\infty} S(r) = \underline{a}. $$$$

We next turn to our third main theoretical result: uniqueness of a stationary equilibrium.

Proposition 5 (Uniqueness of Stationary Equilibrium)

Assume that the intertemporal elasticity of substitution is weakly greater than one for all consumption levels

$$$$ \begin{equation}\label{eq:IES}\textrm{IES}(c):=-\frac{u'(c)}{u''(c)c} \geq 1 \quad \mbox{for all} \ c \geq 0,\end{equation} $$$$

(35)

and that the borrowing constraint takes the form of a strict no-borrowing limit
$$a \geq 0$$
⁠. Then:

Individual consumption
$$c_j(a;r)$$
is strictly decreasing in
$$r$$
for all
$$a > 0$$
and
$$j=1,2$$
.
Individual saving
$$s_j(a;r)$$
is strictly increasing in
$$r$$
for all
$$a > 0$$
and
$$j=1,2$$
.
An increase in the interest rate leads to a rightward shift in the stationary distribution in the sense of first-order stochastic dominance:
$$G(a;r) = G_1(a;r) + G_2(a;r)$$
is decreasing in
$$r$$
for all
$$a$$
in its support.
Aggregate saving
$$S(r)$$
is strictly increasing and hence our continuous-time version of Huggett’s economy has at most one stationary equilibrium.

We briefly sketch key steps in the proof. Part 1 makes use of an important result by Olivi (2018) who analyses the continuous-time income fluctuation problem put forth in the current paper and shows that the consumption response to a change in the interest rate can be decomposed into substitution and income effects as:

$$$$ \frac{\partial c_j(a)}{\partial r} = \underbrace{\frac{1}{u''(c_0)}\mathbb{E}_0\int_0^T e^{-\int_0^t \xi_s ds} u'(c_t)dt}_{\textrm{substitution effect} < 0} + \underbrace{\frac{1}{u''(c_0)}\mathbb{E}_0\int_0^T e^{-\int_0^t \xi_s ds} u''(c_t)(\partial_a c_t)a_t dt}_{\textrm{income effect} > 0} $$$$

with

$$\xi_t := \rho - r + \partial_a c_t>0$$

and where

$$T:= \inf\{t \geq 0 \vert a_t = 0\}$$

is the stopping time at which wealth reaches the borrowing constraint. Here the expectations are over sample paths

$$(a_t,y_t)$$

starting from

$$(a_0,y_0) = (a,y_j)$$

and

$$\partial_a c_t$$

is short-hand notation for the instantaneous MPC

$$c'_j(a_t)$$

⁠. Olivi (2018) further simplifies these substitution and income effects and expresses them in terms of potentially observable sufficient statistics. We instead pursue a different avenue: we show that a sufficient condition for the substitution effect to dominate the income effect and hence

$$\partial c_j(a)/\partial r < 0$$

for all

$$a>0$$

is that the IES is weakly greater than one.

Part 2 uses the budget constraint

$$s_j(a) = y_j + ra - c_j(a)$$

⁠. If consumption is decreasing in

$$r,$$

then saving is increasing in

$$r$$

for

$$a \geq 0$$

⁠. Because there is a positive mechanical effect of

$$r$$

on saving through interest income

$$ra$$

⁠, the assumption that the IES is greater than one is likely overly strong and saving may also be increasing in

$$r$$

with an IES less than one. Consistent with this, consider the simple deterministic example with CRRA utility in (26): whether consumption is increasing in

$$r$$

depends on the IES

$$1/\gamma$$

⁠; but saving

$$s(a)$$

is increasing in

$$a$$

independently of

$$1/\gamma$$

⁠.³⁴ Future work should try to prove uniqueness under weaker assumptions than the IES being greater than one. Either way,

$$\textrm{IES}(c)\geq 1$$

is an intuitive assumption and it includes the commonly used case of logarithmic utility

$$\textrm{IES}(c)= 1$$

⁠.

Part 3, first-order stochastic dominance, and Part 4, that aggregate saving is increasing in

$$r$$

⁠, both use a simple probabilistic argument: starting from a given initial wealth level, for any given income trajectory, a higher interest rate implies a higher saving flow and hence a higher wealth level. Uniqueness of the stationary equilibrium then follows immediately from the monotonicity of aggregate saving.³⁵

5. Computation

We now describe our algorithm for numerically computing equilibria of continuous-time heterogeneous agent models. We use a finite difference (FD) method based on work by Achdou and Capuzzo-Dolcetta (2010) and Achdou (2013) which is simple, efficient and easily extended to other environments. We explain our method in the context of the baseline heterogeneous agent model of Section 2. But the algorithm is, in fact, considerably more general and applies to any heterogeneous agent model with a continuum of atomistic agents (and without aggregate shocks). It is particularly well-suited for computing transition dynamics and solving problems with non-convexities, a fact we illustrate in Section 7 by computing equilibria of such economies. Codes for these applications (and many more) are available from https://benjaminmoll.com/codes/ in Matlab as well as Python, Julia and C++.

5.1. Computational advantages relative to discrete time

Before explaining our algorithm, we provide a brief overview of some of its computational advantages relative to traditional discrete-time methods. We here list four computational advantages that we consider crucial and that contribute notably to the efficiency gains over traditional methods. The first of these advantages is special to the solution of problems with borrowing constraints. The second to fourth advantages concern the solution of heterogeneous agent models more broadly (e.g. models with heterogeneous firms).

To appreciate the first two advantages, contrast the first-order condition of the continuous-time income fluctuation problem (7) with that in the analogous discrete-time problem. For concreteness also assume CRRA utility (5) so that the two conditions are

$$$$ \begin{align} \label{eq:FOCcont}c^{-\gamma} &= v_j'(a), \quad j=1,2 \qquad \mbox{in continuous time and}\\ \end{align} $$$$

(36)

$$$$ \begin{align} \label{eq:FOCdisc}c^{-\gamma} &\geq \beta \sum_{k=1}^2 \pi_{jk} v'_k(a'), \quad a' = y_j + (1+r)a - c, \quad j=1,2 \end{align} $$$$

(37)

in discrete time, where

$$0 < \beta < 1$$

is a discount factor and

$$\pi_{jk}=\Pr(y'=y_k|y=y_j)$$

are the entries of the Markov transition matrix for the analogous discrete-time income process. The first advantage of our continuous-time approach is that, as explained in Section 2.2, the borrowing constraint (3) only shows up in the boundary condition (10) and therefore the first-order condition (36) holds with equality everywhere in the interior of the state space. In contrast, the discrete-time first-order condition (37) holds with complementary slackness and therefore is an inequality. This is because the borrowing constraint may bind one time period ahead. Continuous time allows us to completely sidestep any technical difficulties arising due to such occasionally binding constraints.

Second and related, the first-order condition in (36) is “static” in the sense that it only involves contemporaneous variables. Given (a guess for) the value function

$$v_j(a)$$

it can be solved by hand:

$$c_j(a) = (v_j'(a))^{-1/\gamma},j=1,2$$

⁠. In contrast, the discrete-time condition (37) defines the optimal choice only implicitly. Typical solution methods therefore employ costly root-finding operations. Our continuous-time approach again sidesteps this difficulty.³⁶ Intuitively, in discrete-time dynamic programming, the first-order condition is concerned with tradeoffs between “today” and “tomorrow” while all relevant information from tomorrow onwards is encoded in the continuation value function; in contrast, in continuous-time dynamic programming, the value function encodes all relevant information from today onwards and the first-order condition is therefore static.³⁷

The third advantage of continuous time is a form of “sparsity.” To solve the HJB and KF equations (7) and (8), we discretize these so that their solution boils down to solving systems of linear equations. The resulting matrices are typically extremely sparse, namely “tridiagonal” or at least “block-tridiagonal.” This sparsity generates considerable efficiency gains because there are well-developed routines for solving sparse linear systems, either implemented as part of commercial software packages like Matlab or open-source libraries like SuiteSparse. The reason that tridiagonal matrices arise is that a discretized continuous-time process either stays at the current grid point, takes one step to the left or one step to the right. But, it never jumps.³⁸

Fourth, in all heterogeneous agent models, there is a tight link between solving the HJB and KF equations. One can typically “kill two birds with one stone” in the sense that, having computed the solution to the HJB equation one gets the solution to the KF equation “for free”: the matrix in the discretized version of the latter is the transpose of the matrix in that of the former. The underlying mathematical reason is that the KF equation is the “transpose problem” of the HJB equation or, more precisely, that the differential operator in the KF equation is the adjoint of the operator in the HJB equation.³⁹

5.2. Bird’s eye view of algorithm for stationary equilibria

Our aim is to calculate stationary equilibria—functions

$$v_1,v_2$$

and

$$g_1,g_2$$

and a scalar

$$r$$

satisfying (7), (8), and (11)—given a specified function

$$u$$

⁠, and values for the parameters

$$\rho,\lambda_1,\lambda_2$$

and

$$\underline{a}$$

⁠. Transition dynamics are the subject of Section 5.5. Before, we describe the algorithm in detail, we provide a bird’s eye view of the algorithm’s general structure. We focus on two distinct challenges. First, the HJB and KF equations describing a stationary equilibrium are coupled and one therefore has to iterate on them somehow. Second, solving these differential equations requires approximating the value function and distribution.

Iterating on the equilibrium system. From a bird’s eye perspective our algorithm for solving the stationary equilibrium shares many similarities with algorithms typically used to solve discrete-time heterogeneous agent models. In the context of our Huggett economy, we use a bisection algorithm on the stationary interest rate. We begin an iteration with an initial guess

$$r^0$$

⁠. Then for

$$\ell=0,1,2,...$$

we follow

Given
$$r^\ell$$
⁠, solve the HJB equation (7) using a FD method and calculate
$$s_{j}^{\ell}(a)$$
⁠.
Given
$$s_{j}^{\ell}(a)$$
⁠, solve the KF equation (8) for
$$g_{j}^{\ell}(a)$$
using a FD method.
Given
$$g_{j}^{\ell}(a)$$
⁠, compute the net supply of bonds
$$S(r_\ell) = \int_{\underline{a}}^\infty a (g_{1}^{\ell}(a) + g_{2}^{\ell}(a))da$$
and update the interest rate: if
$$S(r^{\ell})>B$$
⁠, decrease it to
$$r^{\ell+1}<r^{\ell}$$
and vice versa.

When

$$r^{\ell+1}$$

is close enough to

$$r^{\ell}$$

⁠, we call

$$(r^{\ell},v_{1}^{\ell},v_{2}^{\ell},g_{1}^{\ell},g_{2}^{\ell})$$

a stationary equilibrium. As already noted, this continuous-time algorithm is extremely close to typical discrete-time algorithms. It instead differs in the solutions of the dynamic programming equation and the equation for the distribution, and it is this difference that is the source of the resulting efficiency gains.

Discretization of the equilibrium system. In order to solve the differential equations (7) and (8), the value function and distribution need to be approximated in some fashion. We explain our approach—a FD method—in more detail in the next two subsections. But a brief sketch is as follows. In a nutshell, the key idea is that this FD method transforms our system of differential equations into a system of sparse matrix equations. With this goal in mind, we approximate both

$$v_1,v_2$$

and

$$g_1,g_2$$

$$I$$

discrete points in the space dimension,

$$a_i, i=1,...,I$$

⁠. Denote the value function and distribution along this discrete grid using the vectors

$$\mathbf{v}=\left(v_1(a_{1}),...,v_1(a_{I}),v_2(a_{1}),...,v_2(a_{I})\right)^{\rm T}$$

and

$$\mathbf{g}=\left(g_1(a_{1}),...,g_1(a_{I}),g_2(a_{1}),...,g_2(a_{I})\right)^{\rm T}$$

⁠; both

$$\mathbf{v}$$

and

$$\mathbf{g}$$

are of dimension

$$2I$$

⁠, the total number of grid points in the individual state space. The end product of our discretization method will be the following system of matrix equations:

$$$$ \begin{align} \label{eq:HJB_discrete}\rho \mathbf{v} &= \mathbf{u}(\mathbf{v}) + \mathbf{A}(\mathbf{v};r)\mathbf{v},\\ \end{align} $$$$

(38)

$$$$ \begin{align} \label{eq:KFE_discrete}\mathbf{0} &= \mathbf{A}(\mathbf{v};r)^{\rm T}\mathbf{g},\\ \end{align} $$$$

(39)

$$$$ \begin{align} \label{eq:eq_discrete}B &= S(\mathbf{g};r). \end{align} $$$$

(40)

The first equation is the discretized HJB equation (7), the second equation is the discretized KF equation (8) and the third equation is the discretized market clearing condition (11). The

$$2I \times 2I$$

matrix

$$\mathbf{A}(\mathbf{v};r)$$

has the interpretation of a transition matrix that captures the evolution of the idiosyncratic state variables in the discretized state space. It turns out to be extremely sparse.

$$\mathbf{A}(\mathbf{v};r)^{\rm T}$$

in the second equation denotes the transpose of that same matrix, i.e. the discretized KF equation is the “transpose problem” of the discretized HJB equation. As already noted, (38)–(40) is a system of sparse matrix equations that is easy to solve on a computer by following (the analogues of) Steps 1 to 3 in the previous paragraph.

5.3. Step 1: Solving the HJB equation

For Step 1, we solve the HJB equation (7) using a FD method. We now explain this approach. The Supplementary Appendix contains a more detailed explanation.

Theory for numerical solution of HJB equations (Barles–Souganidis). Before we explain our approach, we note that there is a well-developed theory concerning the numerical solution of HJB equations using FD schemes in the same way as there is a well-developed theory concerning the numerical solution of discrete-time Bellman equations. The key result is due to Barles and Souganidis (1991) who have proven that, under certain conditions, the solution to a FD scheme converges to the (unique viscosity) solution of the HJB equation. The interested reader should consult Barles and Souganidis’ original (and relatively accessible) paper or the introduction by Tourin (2013). In short, for their result to hold, the FD scheme needs to satisfy three conditions: (1) “monotonicity”, (2) “stability” and (3) “consistency.” Here, it suffices to note that (2) and (3) are typically easy to satisfy and, in practice, the main difficulty is to design a FD scheme that is “monotone.”

Finite difference method. We here explain the FD method for solving the stationary HJB equation for the special case with no income uncertainty

$$y_1=y_2=y$$

⁠:

$$$$ \begin{equation}\label{eq:HJBsimple}\rho v(a) = \max_c \ u(c) + v'(a)(y + ra - c). \end{equation} $$$$

(41)

The generalization to income risk is straightforward. As already mentioned, the FD method approximates the function

$$v$$

$$I$$

discrete points in the space dimension,

$$a_i,i=1,...,I$$

⁠. We use an equispaced grid with distance

$$\Delta a$$

between points and denote

$$v_{i} := v(a_i)$$

⁠. We approximate the derivative

$$v_{i}'=v'(a_i)$$

with either a forward or a backward difference approximation

$$$$ \begin{equation*} \begin{split} v'(a_i) \approx \frac{v_{i+1} - v_{i}}{\Delta a}=: v_{i,F}' \qquad \mbox{or} \qquad v'(a_i) \approx \frac{v_{i} - v_{i-1}}{\Delta a}=: v_{i,B}'. \end{split} \end{equation*} $$$$

The FD approximation to (41) is then

$$$$ \begin{equation}\label{eq:vvv}\rho v_i = u(c_i) + v_i' s_i, \quad s_i:=y + ra_i - c_i, \quad c_i = (u')^{-1}(v_i'), \quad i=1,...,I \end{equation} $$$$

(42)

where

$$v_i'$$

is either the forward or backward difference approximation. Which of the two approximations is used where in the state space is extremely important because this choice determines whether Barles and Souganidis’ monotonicity condition is satisfied.

Upwinding. As just mentioned, it is important whether and when a forward or a backward difference approximation is used. The ideal solution 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, saving

$$s_{i} = y + ra_{i} - c_{i}$$

⁠) is positive and to use a backward difference whenever it is negative. This is intuitive: if saving is positive, what matters is how the value function changes when wealth increases by a small amount; and vice versa when saving is negative. The right thing to do is therefore to approximate the derivative in the direction of the movement of the state. To this end use the notation

$$s_i^+ = \max\{s_i,0\}$$

⁠, i.e.

$$s_i^+$$

is “the positive part of

$$s_i$$

” and analogously

$$s_i^{-} = \min\{s_i,0\}$$

⁠. The upwind version of (42) is

$$$$ \begin{equation}\label{eq:vFDi}\rho v_i = u(c_i) + \frac{v_{i+1} - v_i}{\Delta a} s_{i}^+ + \frac{v_{i} -v_{i-1}}{\Delta a} s_{i}^-, \qquad i=1,...,I,\end{equation} $$$$

(43)

where

$$c_i = (u')^{-1}(v_i')$$

uses an FD approximation

$$v_i'$$

that depends on the sign of

$$s_i$$

in the same way. This simplified exposition ignores two important and related issues. First, that the HJB equation (41) is highly non-linear due to the presence of the max operator, and therefore so is its FD approximation (43). It therefore has to be solved using an iterative scheme and one faces a choice between using so-called “explicit” and “implicit” schemes. Related, from the first-order condition

$$c_i = (u')^{-1}(v_i')$$

⁠, saving

$$s_i$$

and consumption

$$c_i$$

themselves depend on whether the forward or backward approximation is used so (43) has a circular element to it. The solution to both these issues is described in detail in the Supplementary Appendix.

The upwind FD scheme for the HJB equation (43) can be conveniently written in matrix notation. Denoting by

$$\mathbf{v}=(v_1,...,v_I)^{\rm T}$$

the vector collecting the value function at different grid points, we have the matrix equation (38). The matrix

$$\mathbf{A}(\mathbf{v};r)$$

has a special structure: first, it is sparse; more precisely, it is tridiagonal: all entries are zero except for those on the main diagonal, the first diagonal below this, and the first diagonal above the main diagonal. Second, all diagonal entries are negative and given by

$$\frac{s_i^{-}}{\Delta a} - \frac{s_i^{+}}{\Delta a} \leq 0$$

and all off-diagonal entries are positive and given by

$$-\frac{s_i^{-}}{\Delta a} \geq 0$$

and

$$\frac{s_i^{+}}{\Delta a} \geq 0$$

⁠. Third, all rows of

$$\mathbf{A}(\mathbf{v};r)$$

sum to zero. All these properties are extremely intuitive. In effect, the FD method approximates the continuous state variable

$$a$$

with a discrete-state Poisson process on the grid

$$a_i,i=1,...,I$$

and the matrix

$$\mathbf{A}(\mathbf{v};r)$$

summarizes the corresponding Poisson intensities. The properties noted above are precisely the properties that a Poisson transition matrix needs to satisfy. For these reasons we will sometimes refer to

$$\mathbf{A}(\mathbf{v};r)$$

as “transition matrix.”

Boundary conditions and handling the borrowing constraint. Besides guaranteeing that the Barles–Souganidis monotonicity condition holds, an upwind scheme like (43) has an additional advantage: the handling of boundary conditions. First, consider the upper end of the state space

$$a_I$$

⁠. If it is large enough, saving is negative

$$s_I<0$$

so that

$$s_I^+=0$$

⁠.⁴⁰ Therefore from (43) the forward difference is never used at the upper end of the state space. As a result no boundary condition is needed. Next, consider the lower end of the state space and how to impose the state constraint boundary condition (10) which holds with equality only when the constraint binds. This is where the special structure of the upwind scheme comes in: we set

$$v_{1,B}' = u'(y + r a_1)$$

but only for the backward difference approximation and not for the forward difference approximation

$$v_{1,F}'$$

which is instead computed as

$$(v_2-v_1)/\Delta a$$

⁠. Then, we let the upwind scheme select by itself whether this boundary condition is used. From (43), we see that the boundary condition is only imposed if it would be the case that

$$s_1<0$$

⁠; but it is not used if

$$s_1>0$$

⁠. This ensures that the borrowing constraint is never violated.

5.4. Step 2: Solving the Kolmogorov forward equation

For Step 2, consider the stationary KF equation (8). We again discretize the equation using a FD scheme. In contrast to the HJB equation which is non-linear in

$$v$$

⁠, the KF equation is linear in

$$g$$

⁠. Its discretized counterpart can therefore be solved in one iteration. There are a number of admissible FD schemes, but one is particularly convenient and well-founded: the discretization (39) which involves the transpose of the transition matrix

$$\mathbf{A}(\mathbf{v},r)$$

⁠.

The deep underlying reason for this choice of discretization is that the KF equation actually is the “transpose” problem of the HJB equation. More precisely, the differential operator in the KF equation (8) is the adjoint of the operator in the HJB equation (7), the “infinitesimal generator.”⁴¹ Our transpose discretization of the KF equation (39) is not only well-founded mathematically; it is also extremely convenient: having solved the HJB equation, the solution of the Kolmogorov Forward equation is essentially “for free.”

The same numerical method—building the matrix

$$\mathbf{A}$$

and then working with its transpose – can also be used when solving problems that involve only the KF equation, e.g., because the optimal decision rules can be solved for analytically. This approach is, for example, pursued in Jones and Kim (2018) and Gabaix et al. (2016).

Does the presence of a Dirac point mass in the stationary wealth distribution

$$g_1$$

cause problems for our FD method? Supplementary Appendix F.2 explains why the answer is “no.” First, we show theoretically that the only implication of the Dirac mass is that some care is required when interpreting the numerical output, in particular the first element of the vector

$$\mathbf{g}$$

(corresponding to the density of income type

$$y_1$$

at the point

$$a=\underline{a}$$

⁠). Second, we use the analytic solution for the wealth distribution in Proposition 3 as a test case for our numerical algorithm and show that it performs extremely well in practice unless the wealth grid is very coarse.

5.5. Computing transition dynamics

The algorithm to calculate time-varying equilibria—functions

$$v_1,v_2,g_1,g_2,$$

and

$$r$$

satisfying (4), (12), and (13) given an initial condition (16) and a terminal condition (17)—is the natural generalization of that used to compute stationary equilibria. The main complication is that we now need to iterate on an entire function

$$r(t)$$

rather than just a scalar

$$r$$

⁠. We begin an iteration with an initial guess

$$r^0(t),t \in (0,T)$$

⁠. Then for

$$\ell=0,1,2,...,$$

we follow

Given
$$r^\ell(t)$$
and (17), solve the HJB equation (12), marching backward in time. Calculate the saving policy function
$$s_j^\ell(a,t)$$
⁠.
Given
$$s_j^\ell(a,t)$$
and (16), solve the KF equation, marching forward in time, for
$$g_{j}^{\ell}(a,t)$$
⁠.
Given
$$g_{j}^\ell(a,t)$$
⁠, compute the net supply of bonds
$$S^\ell(t) = \int_{\underline{a}}^\infty a (g_{1}^{\ell}(a,t) + g_{1}^{\ell}(a,t))da$$
and update the interest rate as
$$r^{\ell+1}(t) = r^{\ell}(t) - \xi(t) \frac{d S^\ell(t)}{dt}$$
where
$$\xi(t)>0$$
⁠.⁴²

When

$$r^{\ell+1}(t)$$

is close enough to

$$r^{\ell}(t)$$

for all

$$t$$

⁠, we call

$$(r^{\ell},v^\ell_{1},v_{2}^{\ell},g_{1}^{\ell},g_{2}^{\ell})$$

an equilibrium.

The FD method used for computing the time-dependent HJB and KF equations for a given time path

$$r^\ell(t)$$

is similar to that for their stationary counterparts. In addition to discretizing wealth

$$a$$

⁠, we now also discretize time

$$t$$

on a grid

$$t^n,n=1,...,N$$

⁠, here with equal-sized time steps of length

$$\Delta t$$

⁠. Denoting by

$$\mathbf{v}^n$$

and

$$\mathbf{g}^n$$

the stacked, discretized value function and distribution at time

$$t^n$$

⁠, the time-dependent counterpart to (38)–(40) is

$$$$ \begin{equation}\label{eq:time-discretized} \begin{split} \rho \mathbf{v}^n &= \mathbf{u}(\mathbf{v}^{n+1}) + \mathbf{A}(\mathbf{v}^{n+1};r^n)\mathbf{v}^n + \frac{\mathbf{v}^{n+1} - \mathbf{v}^n}{\Delta t} ,\\ \frac{\mathbf{g}^{n+1} - \mathbf{g}^n}{\Delta t} &= \mathbf{A}(\mathbf{v}^n;r^n)^{\rm T}\mathbf{g}^{n+1},\\ B &= S(\mathbf{g}^n), \end{split} \end{equation} $$$$

(44)

for time steps

$$n=1,...,N$$

⁠, with terminal condition

$$\mathbf{v}^N = \mathbf{v}$$

where

$$\mathbf{v}$$

is the steady state solution to (38) and with initial condition

$$\mathbf{g}^1 = \mathbf{g}_0$$

⁠.

An alternative to Step 3 in the algorithm above is to view Steps 1 and 2 as defining excess supply as a function of the entire time path

$$r(t),t\geq 0$$

and to solve it by means of a root-finding method like Newton’s method or a quasi-Newton method like Broyden’s method. This is easiest to understand in the context of the discretized system (44): given an

$$N$$

-dimensional vector of interest rates

$$\mathbf{r}=(r^1,...,r^N)$$

⁠, the system (44) defines an

$$N$$

-dimensional excess supply function

$$\mathbf{f}:\mathbb{R}^N \rightarrow \mathbb{R}^N$$

⁠, with the

$$n$$

th element equal to excess supply at the

$$n$$

th time step,

$$S(\mathbf{g}^n)-B$$

⁠. One then uses a root-finding method to find the vector

$$\mathbf{r}^*$$

such that

$$\mathbf{f}(\mathbf{r}^*)=\mathbf{0}$$

⁠. All this is explained in more detail in the Supplementary Appendix.

Our computational method for transition dynamics can, of course, also be used to compute (non-linear) impulse responses to unanticipated aggregate shocks (“MIT shocks”). It should be straightforward to use a linearized counterpart to compute linear impulse responses to small MIT shocks along the lines of Boppart et al. (2018) and Auclert et al. (2019) to obtain further speed gains.

5.6. Performance and comparison to a discrete-time method

Section 5.1 listed a number of computational advantages of continuous time relative to discrete time. We here substantiate these claims and compare the computational performance of our method to that of a state-of-the-art discrete-time method. We also explain how one may assess the accuracy of our solution method and discuss the relation to traditional discrete-time accuracy metrics like Euler equation errors.

A test problem. As a basis for these comparisons, we use a standard partial-equilibrium income fluctuation problem with a fixed interest rate. We focus on a partial-equilibrium problem because our strategy for iterating on equilibrium prices and computing transition dynamics is identical to that of standard discrete-time methods, i.e., any speed gain due to continuous time will necessarily occur in partial equilibrium. However, we now consider a version with a richer income process. We consider both a continuous-time and a discrete-time version and specify these to be as comparable as possible. For example, the discrete-time version features an AR(1) process for the logarithm of income and the continuous-time version features an Ornstein–Uhlenbeck process, a diffusion process that is the natural continuous-time analogue of an AR(1) process, both with comparable persistence and standard deviation and with a stationary mean of one. Other details about the specification and parameterization are in Supplementary Appendix F.1. We solve the continuous-time version using our FD method and the discrete-time version using the endogenous grid method (Carroll, 2006).

Assessing accuracy. When comparing different computational methods, we want to understand which method is faster while keeping the accuracy of the numerical solution constant. We therefore require a metric for assessing this accuracy. A challenge in this regard is that, standard discrete-time accuracy metrics such as Euler equations errors are not applicable in continuous time. To see this, consider the analogous discrete-time problem in Section 5.1 with Euler equation (37). As explained by Santos (2000), the rationale for examining the residuals in this Euler equation is that it is the first-order condition of the maximization problem in the Bellman equation. And by bounding the error in this first-order condition, one can bound the error in the policy function and more importantly the value function, i.e. the welfare loss from suboptimal behaviour due to numerical error. But for HJB equations like (41) and the associated finite-difference approximation (42), there is no error in the first-order condition (36) because it can be solved by hand. Instead, any error in the numerical solution of this PDE stems only from the finite-difference approximation of its derivatives. This is explained in more detail in Supplementary Appendix F.1 where we also briefly discuss other candidate accuracy metrics from the mathematics literature on HJB equations.

Given this, we use the following pragmatic approach. We first solve the two income fluctuation problems using an extremely fine grid with

$$10,000$$

wealth grid points. We then treat this solution as the “true solution” and assess the accuracy of numerical solutions with coarser grids in terms of the discrepancy from this solution. We present results for two accuracy metrics: (1) the mean percentage error in the policy function averaged over all wealth and income states and (2) the percentage error in stationary aggregate consumption.

Speed-accuracy tradeoff. Figure 9(a) plots speed-accuracy tradeoffs for the continuous- and discrete-time solution methods of our test problem, using the policy function error as our accuracy metric. Each blue circle corresponds to a continuous-time computation but with a different number of grid points ranging from a very coarse discretization with 25 grid points to an extremely fine one with 10,000 grid points (the “true” solution). Similarly, each red cross corresponds to a discrete-time computation. For each computation, the figure plots the time until the algorithm converged measured in seconds (vertical axis) against the policy function error (horizontal axis). Starting with a coarse grid in the lower right and increasing the number of grid points, the computations become slower but more accurate and we move toward the upper left. By varying the number of grid points, the blue line with circles therefore traces out a continuous-time speed-accuracy tradeoff. Similarly, the red line with crosses traces out the analogous discrete-time tradeoff.

Figure 9

Computational speed and accuracy: continuous versus discrete time

Notes: The figure reports speed and accuracy measures for the numerical solution of an income fluctuation problem in both continuous and discrete time. See Supplementary Appendix F.1 for a detailed description of the exercise. Each blue circle correspond to a continuous-time computation with a different number of grid points ranging from 25 to 10,000 grid points (the “true” solution). Similarly, each red cross corresponds to a discrete-time computation. (a) uses the policy function error as the accuracy metric and (b) uses the percentage error in aggregate consumption. The code is available at https://benjaminmoll.com/comparison/.

Open in new tab Download slide

The key takeaway from the figure is that the continuous-time speed-accuracy tradeoff strictly dominates its discrete-time counterpart: for any given policy function error, the continuous-time method is always faster; conversely, for any given computational speed, the continuous-time method is always more accurate. Supplementary Appendix Figure A.1(a) reports the ratio of the computational times for different accuracy levels. It shows that the continuous-time method is at least twice as fast as the discrete-time method but can be more than 30 times faster if high accuracy (low error) is required.

Figure 9(b) repeats this exercise but instead using the percentage error in aggregate consumption as the accuracy metric. The computations are now somewhat more time-intensive because they require computing stationary distributions in addition to policy functions. The difference in computational performance is even more striking: for a given level of accuracy the continuous-time method is between 10 and 500 times faster—see Supplementary Appendix Figure A.1(a) which plots the relative speed.

General equilibrium and transition dynamics. As already noted, our strategy for iterating on equilibrium prices and computing transition dynamics is the same as in standard discrete time problems and we therefore do not conduct a comparison between the two. We nevertheless briefly comment on our method’s performance for computing these. All of Figures 1, 6, and 8 for the Huggett economy earlier in the paper were computed using a Matlab implementation of the algorithm in Section 5.2. Even though we work with a fine wealth grid with

$$I=1000$$

grid points, solving for a stationary equilibrium takes about 0.25 seconds on a MacBook Pro laptop computer. Next, consider the corresponding transition dynamics. With

$$I=1000$$

wealth grid points,

$$N=400$$

time steps and the same hardware, computing (12) and (13) for a fixed time path

$$r(t)$$

takes about 2 s. Iterating on

$$r(t)$$

until an equilibrium transition is found takes about 4 minutes (even though market clearing conditions like (4) that implicitly define prices are difficult to impose during transitions).⁴³

5.7. Finite difference methods in economics and alternatives

Candler (1999) has previously used a FD method to solve HJB equations arising in economics and also discusses upwinding. Our numerical method adds to his in three dimensions. First and most obviously, we consider coupled HJB and KF equations rather than just the HJB equation in isolation: the system (7), (8), and (11) rather than just (7). Second, even when considered in isolation, our HJB equation differs from Candler’s because it features a borrowing constraint—a ubiquitous feature of heterogeneous agent models—and we design an upwind method that respects this constraint. Third, we show that our solution method has well-developed theoretical underpinnings by making the connection to the Barles and Souganidis (1991) theory.

Another method that is closely related to FD methods and has been previously used in economics is the Markov-chain approximation method (MCAM) of Kushner and Dupuis (2013). See, e.g., Golosov and Lucas (2007) and Barczyk and Kredler (2014a,b, 2018). One way of viewing our FD method is as a simple special case of MCAM. As we explained in Section 5.3, our FD method effectively approximates the law of motion for continuous state variables with a discrete-state Poisson process; that is, we use the FD method to build an approximating Markov chain. MCAM is a more general approach to building such approximating Markov chains, e.g., via trinomial trees.

Besides the FD and Markov-chain approximation (MCA) methods, there are many alternative methods for solving partial differential equations in general and HJB and KF equations in particular. Examples include finite-element, finite-volume, and semi-Lagrangean methods as well as approximation via orthogonal (e.g. Chebyshev) polynomials. In principle, these other methods can also be used to solve heterogeneous agent models of the type discussed here; in particular by following the same Steps 1 to 3 laid out in Section 5.2 but simply exchanging the solution method used within Steps 1 and 2. There is no sense in which our FD method dominates these other methods, some of which may even be more accurate for coarse discretizations. We nevertheless prefer the FD method for two reasons. First, it is transparent and easy to implement: in case the algorithm spits out junk, it is usually easy to track down the problem. Second, it delivers a useful symmetry for the HJB and KF equations: the transpose property discussed in Section 5.4 which is typically not shared by other methods. Finally, because the FD method is fast, choosing fine grids is usually sufficiently cheap and the potentially higher accuracy of alternative methods for coarse discretizations is therefore not relevant.

6. Non-convexities and the Power of Viscosity Solutions

Many important economic problems involve non-convexities. These are difficult to handle with standard discrete-time methods. In contrast, viscosity solutions and finite difference methods are designed to handle such problems. To illustrate this, we now present an economy in which the interplay of indivisible housing and mortgages with a down-payment constraint causes a non-convexity. We show that the same algorithm that we used in the standard Aiyagari–Bewley–Huggett model can be used without change. We also use this problem to provide a brief introduction to viscosity solutions and to comment on their usefulness.

6.1. Non-convexities: indivisible housing, mortgages, poverty traps

We here provide a parsimonious example of a theory that features a non-convexity: individuals can take out a mortgage to buy houses subject to a down-payment constraint and housing is indivisible. The purpose of this subsection is not to propose a quantitatively realistic model of housing; rather it is to showcase what kind of models can be solved with our computational algorithm.

Setup. Individuals have preferences over consumption

$$c_t$$

and housing services

$$h_t$$

⁠:

$$$$ \mathbb{E}_0 \int_0^\infty e^{-\rho t} U(c_t,h_t)dt. $$$$

They can borrow and save in a riskless bond

$$b_t$$

and buy housing at price

$$p$$

⁠. The key restriction is that there are no houses below some threshold size

$$h_{\min}>0$$

⁠. That is, an individual can either not own a house

$$h_t=0$$

or one that is larger than

$$ h_{\min}$$

⁠; compactly:

$$$$ h_t \in \{0,[h_{\min},\infty)\}. $$$$

An individual’s budget constraint is

$$\dot{b}_t + p \dot{h}_t = y_t + rb_t - c_t$$

⁠. As before

$$y_t\in \{y_1,y_2\}$$

follows a two-state Poisson process. When buying a house, the individual can take out a mortgage and borrow up to a fraction

$$\theta \in [0,1]$$

of the value of the house:

$$$$ -b_t \leq \theta p h_t. $$$$

Equivalently, the down-payment needs to be at least a fraction

$$1-\theta$$

of the house’s value. The interest rate

$$r$$

and house price

$$p$$

are determined in equilibrium. Housing is in fixed supply normalized to unity and bonds are in zero net supply.

Stationary equilibrium. It is convenient to work with net worth

$$a_t := b_t + ph_t$$

which follows

$$\dot{a}_t = y_t + r(a_t- ph_t) - c_t$$

⁠. The borrowing constraint then becomes

$$p h_t \leq \phi a_t$$

with

$$\phi :=\frac{1}{1-\theta}$$

⁠. Denote the set of admissible housing choices by

$$\mathcal{H}(a) = \{h: ph \leq \phi a \} \cap \{0,[h_{\min},\infty)\}$$

⁠. A stationary equilibrium is fully characterized by the following system of equations

$$$$ \begin{align} \label{eq:HJBss_housing}\rho v_j(a) =& \max_{c,h \in \mathcal{H}(a)} \ U(c,h) + v'_j(a)(y_j + r(a- ph) - c) + \lambda_j(v_{-j}(a) - v_j(a)),\\ \notag 0 =& - \frac{d}{da}[s_j(a)g_j(a)] - \lambda_j g_j(a) + \lambda_{-j} g_{-j}(a),\\ \notag 1 &= \int_0^\infty (h_1(a)dG_1(a)+h_2(a)dG_2(a)), \quad 0 = \int_0^\infty (b_1(a)dG_1(a)+b_2(a)dG_2(a)), \end{align} $$$$

(45)

where

$$c_j(a),h_j(a),s_j(a) = y_j + r(a - ph_j(a)) - c_j(a)$$

and

$$b_j(a)=a-ph_j(a)$$

are the optimal consumption, housing, saving and bond holding policy functions.

In what follows, we solve this equilibrium system under the additional assumption that utility is quasi-linear

$$U(c,h)=u(c+f(h))$$

⁠. This assumption is convenient because the optimal housing choice separates from the consumption-saving problem and allows for a simple connection to theories with non-convex production technologies (Skiba, 1978). That being said, the model can easily be solved numerically for general preferences

$$u(c,h)$$

⁠. Assuming quasi-linear utility and defining

$$q=c+f(h)$$

⁠, (45) becomes

$$$$ \begin{align} \label{eq:HJBsimplified_housing}\rho v_j(a) &= \max_{q} u(q) + v'_j(a)(y_j + F(a) - q) + \lambda_j(v_{-j}(a) - v_j(a)),\\ \end{align} $$$$

(46)

$$$$ \begin{align} \label{eq:fdef}F(a) &= \max_{h \in \mathcal{H}(a)} \{f(h) + r(a-ph)\}. \end{align} $$$$

(47)

The function

$$F$$

is the pecuniary benefit from optimally allocating wealth

$$a$$

between housing and bonds. Figure 10(a) and (b) plot the solution to the optimal housing choice problem as a function of wealth

$$a$$

⁠: panel (a) plots the housing policy function

$$h(a)$$

and panel (b) plots the pecuniary benefit

$$F(a)$$

⁠.

Figure 10

Model with indivisible housing: policy and value functions and multiple stationary distributions.

Open in new tab Download slide

The vertical line in the two graphs is at

$$a^* := p h_{\min}/\phi$$

⁠, the down-payment necessary to buy the smallest available house

$$h_{\min}$$

⁠. An individual with wealth

$$a_t$$

below this threshold simply cannot buy a house at time

$$t$$

⁠, and hence

$$h(a)=0$$

for

$$a \leq a^*$$

⁠. As her wealth increases above

$$a^*$$

⁠, the individual is first up against the constraint

$$p h \leq \phi a$$

so that the size of her house increases linearly with wealth; when her wealth is large enough, she chooses the unconstrained house size given by the solution to

$$f'(h) = rp$$

⁠. Importantly the function

$$F$$

in panel (b) is convex–concave as a function of wealth

$$a$$

as in some theories of economic growth (Skiba, 1978).

Figure 10(c) plots the resulting saving policy function. The black, dashed horizontal line is at zero, i.e., saving is positive above that line and negative below. Optimal saving has the typical feature of problems with non-convexities: for each income type, there is a threshold wealth level (the “Skiba point”) below which individuals decumulate assets and above which they accumulate assets. In Figure 10(c), these points are where the saving policy functions intersect zero while sloping upward. As usual, each “Skiba point” is strictly below the point of the non-convexity

$$a^*$$

(the dashed vertical line). Figure 10(d) plots the corresponding value functions: importantly, they feature convex kinks at the “Skiba point”. The next section discusses how the theory of viscosity solutions deals with such non-differentiability.

Since the theory features classic poverty trap dynamics, there can be multiple stationary wealth distributions. Figures 10(e) and (f) confirm this possibility: they plot two possible stationary wealth distributions. In fact, there is a continuum of stationary wealth distributions. In results not shown here due to space constraints, we have also computed the model’s transition dynamics. Not surprisingly given the discussion thus far, the economy features history dependence in the sense that initial conditions determine where the economy ends up in the long run. As already noted, the point of this subsection is not to argue quantitatively that the presence of indivisible housing and down-payment constraints creates history dependence. Rather it is to showcase the possibilities of our computational algorithm.

6.2. Viscosity Solutions

Because the value function features kinks the HJB equation does not have a classical solution (Section 3.2). Instead, the solution is a particular type of weak solution: a viscosity solution (Crandall and Lions, 1983). Supplementary Appendix D defines this solution concept and provides an overview of the corresponding theory. We here provide a brief summary.

To this end, consider a simplified version of (46) without labour income

$$y_1=y_2=0$$

⁠:

$$$$ \begin{equation}\label{eq:viscosity_example} \rho v(a) = \max_{q} \left\{u(q) + v'(a)(F(a)- q)\right\}, \end{equation} $$$$

(48)

with

$$F$$

defined in (47). As already noted,

$$F$$

is convex–concave (Figure 10(b)). Note the similarity to theories of economic growth with convex–concave production functions (Skiba, 1978) and to theories of entrepreneurship with financial constraints and non-convexities in production, either due to fixed costs in production or to an occupational choice (see e.g. Buera, 2009; Buera et al., 2011; Buera and Shin, 2013; Buera et al., 2015). Our theoretical and computational approaches also carry over to such theories.

Since, we expect the value function

$$v$$

to feature a kink, we require a solution concept to the HJB equation (48) that allows for such points of non-differentiability. In particular, we need to answer the question: what do we mean by saying “

$$v$$

satisfies (48)” if the derivative

$$v'$$

does not exist at some point? The basic idea of viscosity solutions is as follows. We replace the derivative

$$v'$$

at a point where it does not exist with the derivative

$$\phi'$$

of a smooth function

$$\phi$$

(a “test function”) that “touches

$$v$$

⁠.” We then define a viscosity solution as a function

$$v$$

that satisfies an alternative equation that features

$$\phi'$$

instead of

$$v'$$

⁠. Definition 1 in Appendix D spells this out formally. To give the reader a flavour, consider a convex kink in

$$v$$

as in Figure 10(d) and denote the kink point by

$$a^*$$

⁠. Informally, a viscosity solution to (48) is defined as a function

$$v$$

that satisfies

$$$$ \begin{equation}\label{eq:viscosity_intuition} \rho v(a^*) \geq \max_{q}\left\{ u(q) + \phi'(a^*)(F(a^*)- q)\right\} \end{equation} $$$$

(49)

for any smooth function

$$\phi$$

that “touches

$$v$$

from below” and that satisfies another symmetric condition for potential concave kinks. The rough intuition for (49) is as follows. The derivative

$$v'$$

in (48) is the marginal continuation value. At points where this derivative does not exist, we replace the non-differentiable continuation value

$$v$$

with a differentiable one

$$\phi$$

and make use of a “monotonicity” condition familiar from discrete-time dynamic programming: a higher continuation value implies that also today’s value must be weakly higher. In particular, the inequality in (49) is due to

$$\phi$$

representing a worse continuation value than

$$v$$

⁠. The viscosity solution is then defined as the function

$$v$$

such that this monotonicity condition holds for any such function

$$\phi$$

⁠. Importantly, one can also prove that this solution is unique and that it equals the value function of the corresponding sequence problem. Summarizing, the viscosity solution is the natural solution from the point of view of optimal control theory while being able to handle kinks. All of this is discussed in more detail in Supplementary Appendix D.

As discussed in Section 5.1 a computational advantage of continuous time is that first-order conditions are static and can often be solved by hand. This is still true in the presence of non-convexities. For instance, the first-order condition in (48) is simply

$$u'(q) = v'(a)$$

and the optimal policy

$$q^*(a) = (u')^{-1}(v'(a))$$

⁠. In contrast, the analogous discrete-time first-order condition

$$u'(F(a)-a') = \beta v'(a')$$

is no longer sufficient and typically has multiple solutions (because

$$v'$$

jumps at the kink). Intuitively, the discrete-time first-order condition is concerned with tradeoffs between “today” and “tomorrow” while all relevant information from tomorrow onwards is encoded in the continuation value function. But this means that one still has to worry about “crossing the non-convexity” between today and tomorrow. In continuous time, instead, the value function encodes all relevant information about the future from today onwards, thereby eliminating this problem.

Finally, the Barles–Souganidis theory still applies in the presence of kinks and therefore our computational algorithm can be applied without change.

7. Generalizations and Extensions

We here outline a number of generalizations and extensions of the baseline Huggett model in Sections 2 to 5. Details are in Supplementary Appendix G. The purpose of this brief section is to showcase the generality of our methods and in particular the portability of our computational algorithm.

More general income processes. Our baseline model assumed that income

$$y_t$$

takes one of two values, high and low. Supplementary Appendix G.1 extends many of our results to an environment with a continuum of productivity types.⁴⁴ In particular, the computational algorithm laid out in Section 5 carries over without change. This is true even though the system of equations describing an equilibrium will be a system of PDEs rather than a system of ODEs.

Aiyagari model. Supplementary Appendix G.2 extends our results to the case where individuals save in productive capital and income takes the form of labour income as in Aiyagari (1994).

Soft borrowing constraints. Supplementary Appendix G.3 considers a wedge between borrowing and saving rates and characterize the implication for saving behaviour and the wealth distribution. This form of “soft constraint” can explain the empirical observation that wealth distributions typically have a spike at zero net worth and mass both to the left and the right of zero.

Fat tails. As shown in Proposition 3, the stationary wealth distribution in the Huggett economy with a bounded income process is bounded above. More generally, any income process with a thin-tailed stationary distribution generates a thin-tailed wealth distribution. This property of the model is problematic vis-à-vis empirical wealth distributions which typically feature fat upper tails. In Supplementary Appendix G.4, we extend the Huggett model of Section 2 to feature a fat-tailed stationary wealth distribution by introducing idiosyncratic investment risk (Nardi and Fella, 2017; Benhabib and Bisin, 2018).

Multiple assets with adjustment costs. The model in Sections 6 featured two assets: bonds and housing. But portfolio adjustment between the two assets was costless so that they collapsed into one state variable, net worth. Such costless portfolio adjustment is often a bad assumption, in particular when modeling illiquid assets such as housing. On our computational website, we therefore extend our algorithm to multiple assets with kinked (but convex) adjustment costs. Kaplan et al. (2018) argue that such a two-asset structure is important for understanding the monetary transmission mechanism.

Stopping time problems. Problems with multiple assets may also feature non-convex adjustment costs like fixed costs. Individuals then solve stopping time problems (Stokey, 2009). Instead of solving a standard HJB equation, the value function solves a “HJB Variational Inequality” (HJBVI, Øksendal, 1995; Tourin, 2013).⁴⁵ On our computational website, we also generalize our algorithm to such stopping time problems.⁴⁶ McKay and Wieland (2019), Guerrieri et al. (2020), and Laibson et al. (2020) use this algorithm to study durables and housing investment as well as mortgage refinancing decisions. The method also promises to be useful in other applications, e.g. problems involving default by individuals (see e.g. Livshits et al., 2007) or by sovereign states (see e.g. Aguiar et al., 2013; Bornstein, 2020).

8. Conclusion

This article makes two types of contributions. First, we prove a number of new theoretical results about the Aiyagari–Bewley–Huggett model, the workhorse theory of income and wealth distribution in macroeconomics: (1) an analytic characterization of the consumption and saving behaviour of the poor, particularly their marginal propensities to consume; (2) a closed-form solution for the wealth distribution in a special case with two income types; (3) a proof that there is a unique stationary equilibrium if the intertemporal elasticity of substitution is weakly greater than one. Second, we develop a simple, efficient and portable algorithm for numerically solving both stationary equilibria and transition dynamics of a wide class of heterogeneous agent models, including—but not limited to—this model. Both types of contributions were made possible by recasting the Aiyagari–Bewley–Huggett model in continuous time, thereby transforming it into a system of partial differential equations.

It is our hope that the methods developed in this article, particularly the numerical algorithm, will also prove useful in other applications. One potential application is to spatial theories of trade and development as in Rossi-Hansberg (2005) and Allen and Arkolakis (2014). These theories typically feature a continuum of producers and households distributed over a continuum of locations. In dynamic versions, space would simply be an additional variable in the HJB and KF equations. A challenge would be how to solve for equilibrium prices which are typically functions of space rather than a small number of (potentially time-varying) scalars. Related, a second avenue for future research is to explore richer interactions between individuals. In the class of theories, we have considered here, individuals interact only through prices. But for many questions of interest, richer interactions may be important: for instance there may be more “local” interactions in the form of knowledge spillovers or diffusion (see e.g. Perla and Tonetti, 2014; Lucas and Moll, 2014; Burger et al., 2016; Benhabib et al., 2017; Papanicolaou et al., 2020). In principle, the apparatus put forward in this article—the backward–forward MFG system of coupled HJB and KF equations—is general enough to encompass such richer models.

The editor in charge of this paper was Christian Hellwig.

Acknowledgments

This version supersedes an earlier version of the paper entitled “Heterogeneous Agent Models in Continuous Time.” It is supplemented by two online appendices https://benjaminmoll.com/HACT_appendix/ and https://benjaminmoll.com/HACT_Numerical_Appendix/ as well as a website with codes https://benjaminmoll.com/codes/. We are grateful to Fernando Alvarez, Adrien Auclert, Dave Backus, Daniel Barczyk, Roland Bénabou, Jess Benhabib, Jocelyn Boussard, Paco Buera, Lorenzo Caliendo, Dan Cao, Gabe Chodorow-Reich, Wouter Den Haan, Xavier Gabaix, Fatih Guvenen, Mark Huggett, Mariacristina De Nardi, Greg Kaplan, Tatiana Kirsanova, Nobu Kiyotaki, Matthias Kredler, Ellen McGrattan, Giuseppe Moscarini, Galo Nuño, Ezra Oberfield, Alan Olivi, Jesse Perla, Matt Rognlie, Tony Smith, Ivan Werning, Wei Xiong, Stan Zin, and seminar participants at various institutions for useful comments. We also thank Déborah Sanchez for stimulating discussions in early stages of this project and SeHyoun Ahn, Riccardo Cioffi, Adrien Couturier, Xiaochen Feng, Soroush Sabet, Rui Sousa, and Max Vogler for outstanding research assistance.

Supplementary Data

Supplementary data are available at Review of Economic Studies online. And the replication packages are available at https://doi.org/10.5281/zenodo.4357052.

Data Availability Statement

The data underlying this article are available in Zenodo, at https://doi.org/10.5281/zenodo.4357052. A dataset used for the appendix was also derived from sources in the public domain: 2007 Survey of Consumer Finances, accessible at https://www.federalreserve.gov/econres/scf_2007.htm.

Footnotes

Deaton (2016) succinctly summarizes the second and third reasons: “Aggregation needs to be seen, not as a nuisance, but as a hallmark of seriousness [...] While we often must focus on aggregates for macroeconomic policy, it is impossible to think coherently about national well-being while ignoring inequality and poverty, neither of which is visible in aggregate data. Indeed, and except in exceptional cases, macroeconomic aggregates themselves depend on distribution.”

Another important early reference is Imrohoroğlu (1989) who studies a model with both idiosyncratic and aggregate risk. The presence of a storage technology means that the model is set in partial equilibrium (the interest rate is exogenous and equals zero). Because of this difference, we refer to the Aiyagari–Bewley–Huggett model throughout the paper, but without wanting to diminish Imrohoroğlu’s important contribution.

Of course and as is well-known, the unadorned Aiyagari–Bewley–Huggett model is not sufficiently rich to be an empirically realistic theory of income and wealth distribution. Understanding its theoretical properties is nevertheless important, simply because it forms the backbone of much of modern macroeconomics.

The “Kolmogorov Forward equation” is also often called “Fokker–Planck equation.” Because the term “Kolmogorov Forward equation” seems to be somewhat more widely used in economics, we will use this convention throughout the article. But these are really two different names for the same equation.

The theory of “Mean Field Games” is a general and rigorous framework for the analysis of dynamic, stochastic games with a continuum of players. The name comes from an analogy to the continuum limit taken in “Mean Field theory” which approximates large systems of interacting particles by assuming that these interact only with the statistical mean of other particles. In general, MFGs can be written in terms of a so-called “Master equation” which reduces to the “backward–forward MFG system” in the case without aggregate uncertainty. For more on MFGs, see e.g. Guéant et al. (2011) and Cardaliaguet (2013).

The distribution of MPCs determines, for example, the efficacy of fiscal stimulus (e.g. Kaplan and Violante, 2014; Hagedorn et al., 2017), the transition mechanism of monetary policy (e.g. Auclert, 2019; Kaplan et al., 2018), the effect of a credit crunch or house price movements on consumer spending (e.g. Guerrieri and Lorenzoni, 2017; Berger et al., 2018), and the extent to which inequality affects aggregate demand (e.g. Auclert and Rognlie, 2018; Auclert and Rognlie, 2017).

The uniqueness result additionally assumes that individuals cannot borrow. A key step in our proof is an important result by Olivi (2018). Contemporaneous work by Light (2018) derives a uniqueness result under a similar condition on the IES in a more restrictive discrete-time setting: an Aiyagari economy with CRRA utility and Cobb–Douglas production (as well as no borrowing).

For the same reasons, one of the first results that every graduate student learns is that the neoclassical growth model—the representative-agent counterpart to the Aiyagari model—features a unique steady state.

Our numerical method is based on Achdou and Capuzzo-Dolcetta (2010) and Achdou (2013) but modified to handle the particular features of heterogeneous agent models, in particular borrowing constraints.

10.

More precisely, the differential operator in the KF equation is the adjoint of the differential operator in the HJB equation. The adjoint of an operator is the infinite-dimensional analogue of a matrix transpose.

11.

Because first-order conditions are no longer sufficient and standard envelope theorems do not apply.

12.

See, for example, Jovanovic (1979), Moscarini (2005), Alvarez and Shimer (2011), Moll (2014), Stokey (2014), Vindigni et al. (2015), Jones and Kim (2018), Jones (2015), Toda and Walsh (2015), Benhabib et al. (2016), Cao and Luo (2017), and Kasa and Lei (2018). Miao (2005), Luttmer (2007, 2011, 2015, and Benhabib et al. (2017) analyse theories with heterogeneous producers. These papers, like ours, all study economies with a continuum of heterogeneous agents yielding a system of coupled HJB and KF equations. In contrast, other papers study environments with a finite number of heterogeneous agents (typically equal to two). For example, see Scheinkman and Weiss (1986) and applications of their framework by Conze et al. (1993) and Lippi et al. (2013).

13.

Similarly, there are also several discrete-time approaches for retaining tractability in environments with heterogeneous households (e.g. Bénabou, 2002; Krebs, 2003; Heathcote et al., 2014).

14.

Additionally, Bayer, Rendall, and Wälde assume a “natural borrowing constraint” implying that individuals never actually hit that constraint. Another difference is that they characterize individuals’ saving behaviour in terms of a differential equation for its consumption policy function whereas we work with the HJB equation. Lise (2013) is another paper studying a continuous-time partial-equilibrium setting.

15.

Wang (2007) proposes an elegant continuous-time Aiyagari–Bewley–Huggett model that can be solved analytically but at the cost of making two non-standard assumptions on preferences: CARA utility and discount rates that are increasing in past consumption (in the absence of the second assumption, constant absolute risk aversion (CARA) utility implies exploding wealth inequality, and non-existence of a stationary distribution).

16.

As discussed in detail in Aiyagari (1994), if the borrowing limit

$$\underline{a}$$

is less tight than the so-called “natural borrowing limit”, the constraint

$$a_t \geq \underline{a}$$

will never bind and the “natural borrowing limit” will be the effective borrowing limit. In a stationary equilibrium with

$$r>0$$

⁠, the “natural borrowing limit” is

$$a_t \geq -y_1/r$$

where

$$y_1$$

is the lowest income. In an equilibrium with a time-varying interest rate

$$r_t$$

the natural borrowing constraint is

$$a_t \geq -y_1 \int_t^\infty \exp\left(-\int_t^s r_\tau d\tau \right)ds$$

⁠. The natural borrowing constraint ensures that

$$a_t$$

never becomes so negative that the individual cannot repay her debt even if she chooses zero consumption thereafter.

17.

Suppressing dependence on time

$$t$$

⁠, the CDF will decompose as

$$G_j(a) = \tilde{G}_j(a) + m_j \delta_{a = \underline{a}}$$

where

$$d\tilde{G}_j(a) = g_j(a)da$$

⁠,

$$m_j$$

is the mass at

$$\underline{a}$$

⁠, and

$$\delta$$

is the Dirac delta function. Equivalently, the integral of any function

$$\varphi$$

$$\int_{\underline{a}}^\infty \varphi(a) dG_j(a) = \int_{\underline{a}}^\infty \varphi(a) g_j(a)da + \varphi(\underline{a}) m_j$$

⁠. Also note that

$$\int_{\underline{a}}^\infty dG_1(a) + \int_{\underline{a}}^\infty dG_2(a) = 1$$

⁠, that is,

$$G_j(a)$$

is the unconditional distribution of wealth for a given productivity type

$$j=1,2$$

⁠.

18.

The system can be written more compactly as two non-linear partial differential equations in

$$v_j$$

and

$$g_j$$

only, that do not involve a max operator: define the Hamiltonian

$$H(p) = \max_c u(c) - p c$$

⁠, write the saving policy function as

$$s_j(a) = y_j + r a + H'(v_j'(a))$$

and write (7) and (8) as

$$$$ \begin{equation*} \begin{split} \rho v_j(a) &= H(v_j'(a)) + v_j'(a)(y_j + r a) + \lambda_j(v_{-j}(a)-v_j(a)),\\ 0 &= - \frac{d}{da}[(y_j + r a + H'(v_j'(a))) g_j(a)] - \lambda_j g_j(a) + \lambda_{-j} g_{-j}(a). \end{split} \end{equation*} $$$$

19.

This is in contrast to discrete-time formulations where there is a set

$$[\underline{a},a^*)$$

with

$$a^*>\underline{a}$$

such that type

$$1$$

’s borrowing constraint binds for all

$$a\in [\underline{a},a^*)$$

and hence the first-order condition is distorted.

20.

Note that this inequality has very little to do with the inequality in discrete-time first-order conditions due to occasionally binding borrowing constraints—see e.g. equation (37) later in the article. In fact, the two inequalities go in opposite directions. Even though both inequalities result from the presence of borrowing constraints, the logic behind them is completely different.

21.

See Soner (1986a, b) and Capuzzo-Dolcetta and Lions (1990).

22.

The standard notion of a measure-valued solution is only defined on the interior of the state space and therefore cannot be used to deal with a Dirac point mass at the boundary, a feature that arises in our application. We show in Supplementary Appendix E how to extend the standard notion to take into account this possibility.

23.

This uses the extension of Ito’s formula to Poisson processes:

$$\mathbb{E}_t[d u'(c_j(a_t))] = [u''(c_j(a_t))c_j'(a_t)s_j(a_t) + \lambda_j(u'(c_{-j}(a_t)) - u'(c_j(a_t)))]dt$$

⁠.

24.

Assumption 1 is also what determines consumption and saving behaviour at the constraint for a wide class of less standard utility functions, say, with subsistence concerns. For example, in the Stone–Geary case

$$u'(c)=(c-\underline{c})^{-\gamma}$$

for

$$c \geq \underline{c}$$

⁠, the constraint “matters” if

$$\underline{a}>-(y_1-\underline{c})/r$$

⁠. We have also been asked why Assumption 1 is an assumption on risk aversion rather than marginal utility: intuition suggests that a sufficient condition for the constraint to “matter” should be that

$$u'(y_1+r \underline{a})<\infty$$

⁠. This is not true. A counterexample is

$$\underline{a}=-y_1/r$$

and Guiso and Paiella (2008) utility

$$u'(c) = \exp(-\theta c^{\varepsilon}/\varepsilon)$$

with

$$0<\varepsilon<1$$

⁠. Marginal utility at the constraint is then indeed bounded,

$$u'(y_1+r\underline{a})=1<\infty$$

⁠, but risk aversion is not,

$$\underline{R}= - \lim_{a \rightarrow \underline{a}}\theta (y_1+ra)^{\varepsilon-1} = \infty$$

⁠, violating Assumption 1. As Proposition 1’ in the Supplementary Appendix shows, this implies that individuals never hit the constraint in finite time.

25.

Type

$$1$$

’s consumption at the borrowing constraint is
$$\underline{c}_1 = y_1 + r\underline{a}$$
and type
$$2$$
’s consumption
$$\underline{c}_2>\underline{c}_1$$
is a more complicated object determined by the HJB equation (7) for
$$j=2$$
⁠.

26.

A discrete-time analogue of Corollary 1 is derived in Huggett (1997) who proves that for a sufficiently long sequence of “bad” shocks, individuals hit the borrowing constraint (the argument is part of the proof of his Lemma 1). Unlike our formula, Huggett’s result does not include a characterization of the speed at which this happens, i.e., there is no analogue of our

$$\nu_1$$

⁠.

27.

We are indebted to Xavier Gabaix for suggesting the first special case. Also see Holm (2018) who characterizes consumption behaviour with deterministic income, a borrowing constraint and hyperbolic absolute risk aversion (HARA) utility.

28.

In the case

$$r\neq 0$$

⁠, (24) generalizes to

$$a(t) = \frac{\nu}{r^2} \left(r(T-t) - 1 + e^{-r (T-t)}\right)$$

which satisfies

$$a(t) \sim \frac{\nu}{2}(T-t)^2$$

$$t\rightarrow T$$

(apply l’Hôpital’s rule twice to see that

$$\lim_{t \rightarrow T}\frac{a(t)}{(T-t)^2} = \frac{\nu}{2}$$

⁠. Also all other properties of consumption and saving behaviour emphasized here are unchanged for

$$t$$

close to

$$T$$

or, equivalently, for

$$a_0$$

close to

$$\underline{a}$$

⁠.

29.

The proof makes use of a simple homogeneity property: for all

$$\xi>0$$

⁠, the value function

$$v$$

expressed as a function of wealth

$$a$$

and income

$$y$$

satisfies

$$v(\xi a,y) = \xi^{1-\gamma} v(a,y/\xi)$$

⁠. That is, doubling wealth

$$a$$

⁠, besides scaling everything by a factor

$$2^{1-\gamma}$$

⁠, effectively halves income

$$y$$

⁠. Therefore, as wealth becomes large, it is as if the individual had no labour income. And it is well-known that the consumption-saving problem with CRRA and without labour income has an analytic solution with linear policy functions given by (25).

30.

Alternatively, we can take

$$\lambda_1 \rightarrow 0$$

so that the low-income state is close to being absorbing.

31.

Similar expressions can be obtained in the special case from Section 4.3 in which saving and consumption policy functions are linear throughout the wealth distribution and given by (26).

32.

Under Assumption 1, we expect

$$s_1(\underline{a})=0,s_2(\underline{a})>0,g_2(\underline{a})>0$$

so how can

$$s_1(\underline{a})g_1(\underline{a}) + s_2(\underline{a})g_2(\underline{a}) = 0$$

? The answer is that, as

$$a \rightarrow \underline{a}$$

⁠,

$$g_1$$

explodes and

$$s_1g_1$$

converges to a negative constant. Proposition 3 confirms this: see (33) which implies

$$s_j(\underline{a})g_j(\underline{a}) = \kappa_j,j=1,2$$

with

$$\kappa_1<0,\kappa_2>0$$

and

$$\kappa_1+\kappa_2=0$$

⁠.

33.

See for example Figure 1 in Imrohoroğlu (1989) and Figure 17.7.1 in Ljungqvist and Sargent (2004). To see why this must happen, consider a discrete-time Huggett economy with two income states. All individuals with wealth

$$a=\underline{a}$$

who get the high-income draw choose the same wealth level

$$a'>\underline{a}$$

⁠. So if there is a Dirac mass at

$$\underline{a}$$

⁠, there must also be a Dirac mass at

$$a'$$

⁠. But all individuals with wealth

$$a'$$

who get the high-income draw, also choose the same wealth level

$$a''$$

⁠. So there must also be a Dirac mass at

$$a''>a'$$

⁠. And so on. Through this mechanism, the Dirac mass at the borrowing constraint “spreads into the rest of the state space.” In continuous time individuals instead leave the borrowing constraint in a smooth fashion (here according to a Poisson process, i.e., a process with a continuously distributed arrival time).

34.

Differentiating

$$s(a)$$

in (26) yields

$$\frac{\partial s(a)}{\partial r} = \frac{1}{\gamma} \left(a + \frac{y}{r}\frac{\rho}{r}\right)$$

which is positive for all

$$a \geq -y/r$$

as long as

$$r<\rho$$

⁠.

35.

Açıkgöz (2018) provides a numerical example of multiplicity in an Aiyagari model with an IES of

$$1/6.5$$

⁠, i.e., considerably below one.

36.

In this regard, it shares some similarities with the “endogenous grid method” of Carroll (2006). The difference is that in continuous time this also works with “exogenous grids.”

37.

The first-order condition (36) also does not involve an expectation operator as in (37), i.e., a summation or costly numerical integral over future income states. Instead, the HJB equation (7) captures the stochastic evolution of income with an additive terms

$$\lambda_j (v_{-j}(a) - v_j(a))$$

that does not affect the first-order condition.

38.

Except of course if the process is a Poisson process, i.e., if jumps are “built in.” That being said, the sparsity property survives as long as there is at least one continuously moving state variable (like wealth), i.e. not all individual state variables follow discrete-state Poisson processes.

39.

In principle, one can use an analogous approach in discrete time: form a Markov transition matrix over all states and use it to both iterate backward over value functions and forward over distributions. This method is less popular than it should be and researchers often solve for distributions by Monte-Carlo simulation.

40.

In fact, in the special case without uncertainty analysed in the present section and under the assumption

$$r<\rho$$

⁠, saving is negative everywhere in the state space. The condition that

$$a_I$$

needs to be large enough really only matters for the case with uncertainty. That the case without uncertainty is special also applies to the subsequent discussion: without uncertainty the state constraint (10) always holds with equality.

41.

The “infinitesimal generator” is the continuous-time analogue of a discrete-time transition matrix, and the adjoint of an operator is the infinite-dimensional analogue of a matrix transpose. In our context, the infinitesimal generator captures the evolution of the process in

$$(a,y_j)$$

-space. This operator – let us denote it by

$$\mathcal{A}$$

– is defined as follows: for any vector of functions

$$[f_1(a),f_2(a)]^{\rm T}$$

$$$$ \mathcal{A}\left[\begin{matrix}f_1(a) \\ f_2(a) \end{matrix}\right] = \left[\begin{matrix} f_1'(a) s_1(a) + \lambda_1(f_2(a) - f_1(a)) \\ f_2'(a) s_2(a) + \lambda_2(f_1(a) - f_2(a)) \end{matrix}\right]\!. $$$$

Next, one can show that the operator in the KF equation (8) is the adjoint of this operator: denoting by

$$\mathcal{A}^{*}$$

the adjoint of

$$\mathcal{A}$$

⁠, (8) is

$$0 = \mathcal{A}^*\left[\begin{matrix}g_1(a) \\ g_2(a) \end{matrix}\right]$$

⁠. Equation (39) is the discretized version of this problem.

42.

A good initial guess satisfies

$$r^0(T) = r^*$$

where

$$r^*$$

corresponds to the new stationary equilibrium. With such an initial guess it is convenient to choose

$$\xi(t)$$

with

$$\xi(T)=0$$

⁠, e.g.,

$$\xi(t)= \xi_0\left(e^{-\xi_1 t} - e^{-\xi_1 T} \right),\xi_0,\xi_1>0$$

⁠.

43.

In contrast, computing transitions for the Aiyagari model in Supplementary Appendix G.2, where prices are explicit functions of the aggregate capital stock, takes only 1 minute and 40 seconds. The Matlab code for the stationary equilibrium and transition dynamics of the Huggett model are available at https://benjaminmoll.com/huggett_equilibrium_iterate/ and http://benjaminmoll.com/huggett_transition/. The code for the Aiyagari model is at http://benjaminmoll.com/aiyagari_poisson_MITshock/.

44.

Among the theoretical results, we extend Propositions 1, 2, 4, and 5. That is, all propositions from Section 4 except Propositions 3 (the analytic solution for the stationary distribution with two income types).

45.

Economists typically tackle such problems by imposing a “smooth pasting” condition on the boundary between an inaction region and an adjustment region. While this approach is convenient in one dimension when this boundary is a simple threshold, it becomes impractical in multiple dimensions. The HJBVI approach instead avoids imposing a smooth pasting condition (which becomes a result rather than an imposed axiom) and multi-dimensional problems pose no conceptual problem over one-dimensional one.

46.

See the applications labelled “Stopping Time Problem” at https://benjaminmoll.com/codes/.

References

AÇIKGÖZ,

O. T.

(

2018

), “

On the Existence and Uniqueness of Stationary Equilibrium in Bewley Economies with Production

”,

Journal of Economic Theory

173

–

Month:	Total Views:
February 2021	4
March 2021	64
April 2021	45
May 2021	31
June 2021	56
July 2021	87
August 2021	96
September 2021	232
October 2021	419
November 2021	367
December 2021	321
January 2022	2,526
February 2022	1,014
March 2022	750
April 2022	687
May 2022	659
June 2022	510
July 2022	509
August 2022	496
September 2022	607
October 2022	638
November 2022	654
December 2022	599
January 2023	554
February 2023	588
March 2023	690
April 2023	733
May 2023	725
June 2023	650
July 2023	561
August 2023	532
September 2023	642
October 2023	776
November 2023	642
December 2023	503
January 2024	781
February 2024	744
March 2024	674
April 2024	779
May 2024	593

Article Contents

Income and Wealth Distribution in Macroeconomics: A Continuous-Time Approach

Abstract

1. Introduction

2. The Workhorse Model of Income and Wealth Distribution in Macroeconomics

2.1. Setup

2.2. Stationary equilibrium

2.3. Transition dynamics

3. General Heterogeneous Agent Models as MEAN FIELD GAMES

3.1. General backward–forward MFG systems

3.2. Classical versus weak solutions of HJB and KF equations

4. Theoretical Results for the Aiyagari–Bewley–Huggett Model: Consumption, Saving, and Inequality

4.1. An Euler equation

4.2. Consumption and saving behaviour of the poor

4.3. Consumption and saving behaviour of the wealthy

4.4. Marginal propensities to consume and save

4.5. The stationary wealth distribution

4.6. Stationary equilibrium: existence and uniqueness

5. Computation

5.1. Computational advantages relative to discrete time

5.2. Bird’s eye view of algorithm for stationary equilibria

5.3. Step 1: Solving the HJB equation

5.4. Step 2: Solving the Kolmogorov forward equation

5.5. Computing transition dynamics

5.6. Performance and comparison to a discrete-time method

5.7. Finite difference methods in economics and alternatives

6. Non-convexities and the Power of Viscosity Solutions

6.1. Non-convexities: indivisible housing, mortgages, poverty traps

6.2. Viscosity Solutions

7. Generalizations and Extensions

8. Conclusion

Acknowledgments

Supplementary Data

Data Availability Statement

Footnotes

References

Supplementary data

Citations

Views

Altmetric

Email alerts

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only