# Global Solutions to Master Equations for Continuous Time Heterogeneous Agent
Macroeconomic Models

Zhouzhou Gu111Princeton, Department of Economics. Email<EMAIL_ADDRESS>Mathieu Laurière222Shanghai Frontiers Science Center of Artificial
Intelligence and Deep Learning; NYU-ECNU Institute of Mathematical Sciences at
NYU Shanghai. Email<EMAIL_ADDRESS>Sebastian Merkel333University
of Exeter, Department of Economics. Email<EMAIL_ADDRESS>Jonathan
Payne444Princeton, Department of Economics. Email<EMAIL_ADDRESS>555We
would like to thank Adam Rebei for outstanding research assistance throughout
this project. We are also very grateful to the comments and discussion from
Fernando Alvarez, Mark Aguiar, Yacine Ait-Sahalia, Marlon Azinovic, Johannes
Brumm, Markus Brunnermeier, Rene Carmona, Wouter Den Haan, Mikhail Golosov,
Goutham Gopalakrishna, Lars Hansen, Ji Huang, Felix Kubler, Greg Kaplan,
Moritz Lenel, Ben Moll, Galo Nuño, Thomas J. Sargent, Simon Scheidegger, Jesús
Fernández-Villaverde, Yucheng Yang, Jan Žemlička, and Shenghao Zhu. We are
also thankful for the comments from participants in the DSE 2023 Conference on
Deep Learning for Solving and Estimating Dynamic Models, the CEF 2023
Conference, the IMSI 2023 Conference, the 11th Western Conference on
Mathematical Finance, the Macro-finance Reading Group at USC Marshall, the BFI
Workshop at Chicago, and the Princeton Civitas Seminar. We thank Princeton
Research Computing Cluster for their resources.

(September 3, 2024)

###### Abstract

We propose and compare new global solution algorithms for continuous time
heterogeneous agent economies with aggregate shocks. First, we approximate the
agent distribution so that equilibrium in the economy can be characterized by
a high, but finite, dimensional non-linear partial differential equation. We
consider different approximations: discretizing the number of agents,
discretizing the agent state variables, and projecting the distribution onto a
finite set of basis functions. Second, we represent the value function using a
neural network and train it to solve the differential equation using deep
learning tools. We refer to the solution as an Economic Model Informed Neural
Network (EMINN). The main advantage of this technique is that it allows us to
find global solutions to high dimensional, non-linear problems. We demonstrate
our algorithm by solving important models in the macroeconomics and spatial
literatures (e.g. Krusell and Smith, (1998), Khan and Thomas, (2007), Bilal,
(2023)).

_JEL:_ C63, C68, E60, E27.
_Keywords:_ Heterogeneous agents, computational methods, deep learning,
inequality, mean field games, continuous time methods, aggregate shocks,
global solution.

##  1 Introduction

Macroeconomists have great interest in studying models with heterogeneous
agents and aggregate shocks. Recent work has characterized equilibria in these
models recursively in continuous time using “master equations” from the mean
field game theory literature. A major challenge with these recursive
formulations is that the agent distribution of agents states becomes an
aggregate state variable and so the state space becomes infinite dimensional.
There have been attempts to resolve this issue by using perturbations in the
distribution space (e.g. Bilal, (2023), Alvarez et al., (2023), Bhandari et
al., (2023)). In this paper we use tools from the deep learning literature to
characterize global solutions to the master equation. This necessitates
constructing a finite dimensional approximation to the cross-sectional
distribution of agent states. Most existing deep learning approaches are in
discrete time and replace the agent continuum by a finite collection of
agents. We develop and compare the three main approaches for approximating the
distribution in continuous time: imposing a finite number of agents,
discretizing the state space, and projecting onto a collection of finite basis
functions. For each approximation, we show how to characterize general
equilibrium as a high but finite dimensional differential master equation and
how to customize deep learning techniques to compute global numerical
solutions to the differential equation.

We develop solution techniques for a class of continuous time dynamic,
stochastic, general equilibrium economic models with the following features.
There is a large collection of price-taking agents who face uninsurable
idiosyncratic and aggregate shocks. Given their belief about the evolution of
aggregate state variables, agents choose control processes to solve dynamic
optimization problems. When making their decisions, agents face financial
frictions that constrain their behaviour, which potentially breaks
“aggregation” results and makes the distribution of agent states an aggregate
state variable. Solving for the rational expectations equilibrium reduces to
solving a “master” partial differential equation (PDE) that summarizes both
the agent optimization behaviour (from the Hamilton-Jacobi-Bellman equation
(HJBE)) and the evolution of the distribution (from the Kolmogorov Forward
Equation (KFE)). A canonical example of this type of environment is the
continuous time version of the Krusell and Smith, (1998) model, which we (and
many others) use to our solution.

Our solution approach approximates the infinite dimensional master equation by
a finite, but high, dimensional PDE and then uses deep learning to solve the
high dimensional equation. We consider three different approaches for reducing
the dimension of the master equation: the _finite-agent_ , _discrete-state_ ,
and _projection_ methods. The finite agent method approximates the
distribution by a large, finite number of agents. The discrete state method
approximates the distribution by discretizing the agent state space so the
density becomes a collection of masses at grid points. The projection method
approximates the distribution by a linear combination of finitely many basis
functions. Most deep learning macroeconomic papers are in discrete time and
have focused on the finite agent approximation. One exception is Huang, 2023b
, which uses a discrete-state approximation. All of our distribution
approximations preserve the full non-linearity of the model, which is a key
difference to perturbation based approaches.

We solve the finite dimensional approximation to the master equation by
adapting the Deep Galerkin Method (DGM) and the Physics Informed Neural
Networks (PINN) method developed in the applied mathematics and physics
literatures. This involves approximating the value function by a neural
network and then using stochastic gradient descent to train the neural network
to minimise a loss function that summarizes the average error in the master
equation. We calculate average errors by randomly sampling over points in the
state space. We refer to our approach as training Economic Model Informed
Neural Networks (EMINNs).

Although our deep learning approach is relatively simple to describe at an
abstract level, there are many intricate implementation details. A
particularly important detail is choosing how to sample the states on which we
evaluate the master differential equation error. This choice is less relevant
for discrete time approaches where the economy is simulated to calculate
expectations and also for other continuous time deep learning papers that
don’t have high dimensional distributions as states. We consider three
approaches for sampling the distribution. The first is _moment sampling_ that
draws selected moments of the distribution and then samples agent
distributions that satisfy the drawn moments. The second is _mixed steady
state sampling_ that samples random mixtures of steady state distributions at
different fixed aggregate states. The third is _ergodic sampling_ that samples
by regularly simulating the model economy based on the current estimate of the
model solution. We find that moment sampling is most effective for the finite-
agent approximation, a combination of mixed steady state sampling and ergodic
sampling is most effective for the discrete state space approximation, and a
combination of moment sampling and ergodic sampling is most effective for the
projection approach. Deep learning techniques are sometimes referred to as
“breaking the curse of dimensionality” but this is not really true in the
sense that we cannot train a neural network by showing it a dense set of
distributions. Instead, a key “art” to training a neural network is to sample
from an intelligently chosen subspace that helps the neural network to learn
the equilibrium functional relationships in an economically interesting part
of the state space.

In Section 5, we illustrate our techniques by solving a continuous time
version of Krusell and Smith, (1998) using all three of our methods. We take
this model as a “test case” because there are well established approximate
solutions using traditional techniques to which we can compare. We show that
our methods generate solutions with a low master equation error that produce
very similar simulations to contemporary approaches such as Fernández-
Villaverde et al., (2023). In addition, for the version of the model without
aggregate shocks (the Aiyagari, (1994) model), we show that we can match the
finite difference solution to high accuracy and, for the finite-agent
approximation, solve for transition dynamics without a shooting algorithm.

Although all our methods are ultimately effective for solving the Krusell and
Smith, (1998) model, we find they have different strengths and weaknesses. One
point of difference is how much information about the model solution is
required for sampling. We find that only the finite agent method can be solved
without ergodic sampling, which makes its solution more robust away from the
ergodic mean. A second point of difference is how large and complicated the
neural network needs to be. The discrete state space requires a fine grid
(approximately 200 points) while the finite agent method requires a moderate
number of agents (approximately 40) and the projection method only requires 5
basis functions. A third point of difference is the complication of the law
motion for the parameters approximating the distribution. For the finite agent
method, this is simply the law of motion for agent idiosyncratic states.
However, for the discrete state method it is the law of motion for the mass
between grid points and for the projection method it is the law of motion for
projection coefficients, both of which are hard to compute. How to balance
these trade-offs depends upon how the distribution interacts with agent
decision making and the complexity of the KFE. Ultimately, our conclusions
are:

  1. 1.

When only the mean of the distribution enters into the pricing equations, then
the finite agent method with moment sampling is robust and powerful.

  2. 2.

When many aspects of the distribution are important for prices and the KFE
does not involve complicated derivatives, then the discrete state method with
mixed steady state and ergodic sampling is most successful. However, with
complicated KFEs it is a difficult method.

  3. 3.

When it is clear how to customize the projection approach to the economic
problem, then it can effectively solve the model with low dimensions.

In Section 6, we solve two additional models that have more complicated master
equations and illustrate the power of different methods: a model with
heterogeneous firms (similar to a continuous time version of Khan and Thomas,
(2008)) and a spatial model (taken from Bilal, (2023)). We solve the firm
model using a finite agent method and show that this technique copes well with
the additional non-linearity that comes with pricing firm equity. We solve the
spatial model using a discrete set of locations (since finite agent and
projection approximations infeasible) and show this technique is effective
when different parts of the distribution affect different prices.

_Literature Review_ : Our paper is part of the literature attempting to solve
heterogeneous agent continuous time (HACT) general equilibrium economic models
with aggregate shocks. Many researchers have shown that it can be useful to
cast economic models in continuous time (e.g. Achdou et al., 2022a , Ahn et
al., (2018), Kaplan et al., (2018)). Recent papers (e.g. Bilal, (2023)) have
made progress by characterizing equilibrium using the “master equation”
approach developed in the mathematics literature (see Lions, (2011),
Cardaliaguet et al., (2019)). However, the continuous time literature has
lacked a global solution technique for these models, particularly when written
recursively in master equation form. Instead, many of the existing techniques
in the literature focus on perturbation in the distribution and other
aggregate state variables. Our main contribution is to offer a global solution
to the master equation for HACT models with aggregate shocks.

Technically, our approach is part of a literature attempting to use deep
learning to numerically characterize solutions to dynamic general equilibrium
models with heterogeneous agents and aggregate shocks, which in the
mathematics literature are mean-field games with common noise. There have been
many papers attempting to solve discrete time versions of these models (e.g.
Azinovic et al., (2022), Han et al., (2021), Maliar et al., (2021), Kahou et
al., (2021), Bretscher et al., (2022), Germain et al., 2022a , Azinovic and
Žemlička, (2023)). We are focused on deep learning for continuous time models,
which has been the focus for the mathematics literature but less so for the
economics deep learning literature. For continuous time systems, there are two
broad approaches for using deep learning: (i) a “probabilistic” approach that
attempts to solve a system of stochastic differential equations and (ii) an
“analytic” approach that attempts solve a system of partial differential
equations. Papers using the former approach follow Han et al., (2018) and
train the neural network to learn optimal controls across simulations of a
discrete time approximation to the system of stochastic differential equations
(e.g Fouque and Zhang, (2020), Min and Hu, (2021), Carmona and Laurière,
(2022), Germain et al., 2022b , Huang, 2023b , Huang, 2023a , Huang and Yu,
(2024)). This approach has many similarities to the discrete time deep
learning methods in economics, which also use simulations to train the optimal
control problem. Our paper follows the latter (i.e., analytic) approach, which
involves minimizing the loss in the differential equation across randomly
sampled points. We view this as the “true” continuous time approach in the
sense that it does not involve discretizing the time dimension in order to
train the neural network. To implement the analytic solution approach, we
build on the DGM and PINN methods developed by Sirignano and Spiliopoulos,
(2018), Raissi et al., (2017), and Li et al., (2022) to solve problems in the
physics and applied mathematics literatures. The key difficulty we resolve is
that none of the existing deep learning papers are directly applicable to
solving our master equations. Relative to the discrete-time and probabilistic
approaches, we need to minimize loss in a differential equation, which
requires developing new sampling approaches and resolving how to evaluate
derivatives (rather than simulating to approximate expectations). Relative to
the DGM and PINN literatures, we have to work out how handle models with
forward looking optimizing agents and market clearing conditions. Some papers
in the mean-field literature have attempted to address some these issues (e.g.
Al-Aradi et al., (2022); Carmona and Laurière, (2021)) but they do not handle
aggregate shocks (see e.g. Hu and Lauriere, (2022) for a recent survey).

The major challenge with using the analytic approach to solve master equations
with a general state space is that we need to develop and sample from a finite
dimensional representation of the distribution. There are some existing papers
in the mean field literature that work with low dimensional discrete state
spaces (e.g. Perrin et al., (2022), Cohen et al., (2024)) and some papers in
the economics literature that solve continuous time models without complicated
distributions or with one-dimensional distribution approximations (e.g. Duarte
et al., (2024), Gopalakrishna, (2021), Fernández-Villaverde et al., (2023),
Sauzet, (2021), Achdou et al., 2022b , Barnett et al., (2023)). By contrast,
we develop and compare a range of finite dimensional approximations: reducing
to finite agents, discretizing the state space, and projection onto a finite
set of basis functions. An important contribution of our paper is to
systematically explore the full set of distributional approximations and
develop sampling approaches that are customized to those distributional
approximations.

More generally, our technique (and other deep learning techniques) expand the
traditional techniques in economic for solving heterogeneous agent models with
aggregate shocks. One traditional approach has been to fit a statistical
approximation to the law of motion for the key aggregate state variables (e.g.
Krusell and Smith, (1998), Den Haan, (1997), Fernández-Villaverde et al.,
(2023)). As has been extensively discussed in the literature, this approach
works well when the law of motion for the key state variables can be
efficiently approximated as a function of key moments of the distribution (and
so the economy is very close to permitting “aggregation”). By contrast, our
approach can handle economies without near-aggregation results. A second
approach is to take a type of perturbation in the aggregate state and then
solve the resulting problem with matrix algebra (e.g. Reiter, (2002), Reiter,
(2008), Reiter, (2010), Winberry, (2018), Ahn et al., (2018), Auclert et al.,
(2021), Bilal, (2023), Bhandari et al., (2023)). By contrast, we solve the
model globally and so can handle partial differential equations with global
non-linearity. A final approach is to take a low dimensional projection of the
distribution (e.g. Prohl, (2017), Schaab, (2020)). Our work is complementary
to these papers in that it allows for more general, higher dimensional
projections.

This paper is organised as follows. Section 2 describes the general economic
environment that we will be studying and derives the master equation. Section
3 describes the different finite dimensional approximations to the master
equation. Section 4 describes the solution approach. Section 5 applies our
algorithm to a continuous time version of Krusell and Smith, (1998). Section 6
provides additional examples. Section 7 concludes with practical lessons.

##  2 Economic Model

In this section, we outline the class of economic models for which our
techniques are appropriate. At a high level, in economics terminology, we are
solving continuous time, general equilibrium models with a distribution of
optimizing agents who face idiosyncratic and aggregate shocks. In maths
terminology, we are solving mean field games with common noise.

###  2.1 Environment

_Setting:_ The model is in continuous time with infinite horizon. There is an
exogenous one-dimensional666 For ease of exposition, we restrict attention to
models with one-dimensional aggregate shocks. The method is no different for
the case with multi-dimensional aggregate shocks.  aggregate state variable,
$z_{t}$, which evolves according to:

$\displaystyle dz_{t}={}$
$\displaystyle\mu_{z}(z_{t})dt+\sigma_{z}(z_{t})dB_{t}^{0},\quad z_{0}\hbox{
given},$ (2.1)

where $B_{t}^{0}$ denotes a common Brownian motion process. We let
$\mathcal{F}_{t}^{0}$ denote the filtration generated by $B_{t}^{0}$.

_Agent Problem:_ The economy is populated by a continuum of agents, indexed by
$i\in I=[0,1]$. Each agent $i$ has an idiosyncratic state vector,
$x^{i}\in\mathcal{X}\subset\mathbb{R}^{N_{x}}$, that follows:

$\displaystyle dx_{t}^{i}={}$
$\displaystyle\mu_{x}(c_{t}^{i},x_{t}^{i},z_{t},q_{t})dt+\sigma_{x}(c_{t}^{i},x_{t}^{i},z_{t},q_{t})dB_{t}^{i}+\varsigma_{x}(c_{t}^{i},x_{t}^{i},z_{t},q_{t})dJ_{t}^{i},\quad
x^{i}_{0}\hbox{ given},$ (2.2)

where $c_{t}^{i}$ is a one-dimensional control variable chosen by the agent,
$q_{t}\in\mathcal{Q}$ is a collection of aggregate prices in the economy that
will be determined endogenously in equilibrium, $B_{t}^{i}$ denotes an
$N_{B}$-dimensional idiosyncratic Brownian motion process, $J_{t}^{i}$ denotes
a 1-dimensional idiosyncratic Poisson jump process, $\mu_{x}(\cdot)$ is
$N_{x}$–dimensional, $\sigma_{x}(\cdot)$ is an $N_{x}\times N_{B}$–dimensional
matrix, and $\varsigma_{x}(\cdot)$ is $N_{x}$-dimensional.777 We consider a
one-dimensional control and $J_{t}^{i}$ process for notational simplicity. We
do not consider the case where $dB_{t}^{0}$ directly impacts the evolution of
idiosyncratic states.  We let $\lambda(x,z)$ denote the rate at which Poisson
jump shocks arrive given idiosyncratic state $x$ and aggregate state $z$. We
let $\mathcal{F}_{t}^{i}$ denote the filtration generated by $B_{t}^{i}$,
$J_{t}^{i}$, and $B_{t}^{0}$.

Each agent $i$ has a belief about the stochastic price process
$\tilde{\mathsf{q}}=\\{\tilde{q}_{t}:t\geq 0\\}$ adapted to
$\mathcal{F}_{t}^{0}$. Given their belief, agent $i$ chooses their control
process, $\mathsf{c}^{i}=\\{c_{t}^{i}:t\geq 0\\}$ adapted to
$\mathcal{F}_{t}^{i}$, to solve:

$\displaystyle V(x_{0}^{i},z_{0})={}$
$\displaystyle\max_{\mathsf{c}^{i}\in\mathcal{C}}\mathbb{E}_{0}\left[\int_{0}^{\infty}e^{-\rho
t}u(c_{t}^{i},z_{t},q_{t})dt\right]\ s.t.\ \eqref{eq:dstate_agg},\
\eqref{eq:dstate_i},$ (2.3)

where $\rho>0$ is a discount parameter, $u(c_{t}^{i},z_{t},g_{t})$ is the flow
benefit the agent gets and
$\mathcal{C}=\\{c_{t}^{i}\in\mathcal{C}(x^{i}_{t},z_{t},q_{t}):t\geq 0\\}$ is
the set of admissible controls, where $\mathcal{C}(x,z,q)$ denotes the set of
possible actions for a player whose current state is $x$, when the aggregate
state is $z$ and the prices are $q$. This constraint set incorporates any
“financial frictions” that restrict agent choices. A classic example in
economics is that $x_{t}^{i}$ represents agent wealth and the control must
keep $x_{t}^{i}\geq\underline{x}$. We assume that $u$ is increasing and
concave in its first argument.

_Distributions and Markets:_ Let $\mathcal{D}(x_{t}^{i}|\mathcal{F}_{t}^{0})$
be the population distribution across $x_{t}^{i}$ at time $t$, for a given
history $\mathcal{F}_{t}^{0}$. We assume that
$\mathcal{D}(x_{t}^{i}|\mathcal{F}_{t}^{0})$ has a density $g_{t}$, where for
continuous components of $x$ the density is with respect to the Lebesgue
measure and for discrete components of $x$ the density is with respect to the
counting measure. This means we restrict attention to distributions without
mass points in the continuous variables, which potentially requires us to
“smooth” constraints in the model, as we discuss in Section 2.2. We further
assume that $g_{t}\in\mathcal{G}$, where $\mathcal{G}$ is a subspace of the
set of square integrable functions defined on $\mathcal{X}$.

We assume that the economy contains a collection of markets with price vector
$q_{t}$ and market clearing conditions that allow us to solve for $q_{t}$
explicitly in terms of $z_{t}$ and $g_{t}$:

$\displaystyle q_{t}=Q(z_{t},g_{t}),\quad\forall t\geq 0,$ (2.4)

We assume that markets are incomplete in the sense that agents cannot trade
claims directly (or indirectly) on their idiosyncratic shocks $dB_{t}^{i}$ and
$dJ_{t}^{i}$. This means that the idiosyncratic shocks generate a non-
degenerate cross sectional distribution of agent states.

_Equilibrium:_ Given an initial density $g_{0}$, an equilibrium for this
economy consists of a collection of $\mathcal{F}_{t}^{0}$–adapted stochastic
processes, $\\{\mathsf{g},\mathsf{q},\mathsf{z}\\}$, and
$\mathcal{F}_{t}^{i}$–adapted stochastic processes, $\mathsf{c}^{i}$, for each
$i\in I$, that satisfy the following conditions: (i) each agent’s control
process $\mathsf{c}^{i}$ solves problem (2.3) given their belief that the
price process is $\tilde{\mathsf{q}}$, (ii) the equilibrium prices
$\mathsf{q}$ satisfy market clearing condition (2.4), and (iii) agent beliefs
about the price process are consistent with the optimal behaviour of other
agents in the sense that $\tilde{\mathsf{q}}=\mathsf{q}$.

###  2.2 Recursive Representation of Equilibrium

We assume that there exists an equilibrium that is recursive in the aggregate
state variables: $(z,g)\in\mathcal{Z}\times\mathcal{G}$. In this case, beliefs
about the price process can be characterized by beliefs about the evolution of
the distribution $dg_{t}(x)=\tilde{\mu}_{g}(x,z_{t},g_{t})dt$ since the price
vector $q$ can be expressed explicitly in terms of $(z,g)$. We assume that in
this equilibrium, the value function of the agent takes the form
$V:\mathcal{X}\times\mathcal{Z}\times\mathcal{G}\rightarrow\mathbb{R}$,
$(x,z,g)\mapsto V(x,z,g)$, where $V$ is twice differentiable with respect to
$x$ and $z$ and Frechet differentiable with respect to $g$.

_Hamilton Jacobi Bellman Equation (HJBE):_ In principle, the constraint
$c_{t}^{i}\in\mathcal{C}(x^{i}_{t},z_{t},q_{t})$ restricts $x$ to a subspace
$\mathcal{X}^{c}\subset\mathcal{X}$ and imposes the constraint on the
recursive formulation:

$\displaystyle\Psi(c,x,Q(z,g);\ V,D_{x}V,D_{x}^{2}V)\geq 0,\qquad
x\in\mathcal{X}^{c}.$ (2.5)

However, in order to help the neural network train the value function,
throughout we replace the hard constraint on the set of admissible controls by
a flow utility penalty $\psi(c,x,Q(z,g))$ that is larger the “more” that (2.5)
is violated (and where we have left the dependence on the value function
implicit). We provide explicit examples of $\psi$ in our applications in
Section 5 and Section 6. Given a belief about the evolution of the
distribution, $\tilde{\mu}_{g}$, the agent’s value function, $V$, and optimal
choice of control, $c$, jointly solve the HJBE:

$\displaystyle 0=\max_{c\in\mathcal{C}(x,z,g)}\Big{\\{}{}$ $\displaystyle-\rho
V(x,z,g)+u(c,z,Q(z,g))+\psi(c,x,Q(z,g))$ (2.6)
$\displaystyle+(\mathcal{L}_{x}^{c;Q}+\mathcal{L}_{z}+\mathcal{L}_{g}^{\tilde{\mu}_{g}})V(x,z,g)\Big{\\}}$
(2.7)

where the operators are defined as:

$\displaystyle\mathcal{L}_{x}^{c;Q}V(x,z,g)={}$
$\displaystyle\sum_{j=1}^{N_{x}}\partial_{x_{j}}V(x,z,g)\mu_{x,j}(c,x,z,Q(z,g))$
(2.8)
$\displaystyle+\frac{1}{2}\sum_{j,k=1}^{N_{x}}\partial_{x_{j},x_{k}}V(x,z,g)\Sigma_{x,jk}(c,x,z,Q(z,g))$
(2.9)
$\displaystyle+\lambda(x,z)\left(V(x+\varsigma_{x}(c,x,z,Q(z,g)),z,g)-V(x,z,g)\right)$
(2.10) $\displaystyle\mathcal{L}_{z}V(x,z,g)={}$
$\displaystyle\partial_{z}V(x,z,g)\mu_{z}(z)+\frac{1}{2}\partial_{zz}V(x,z,g)\left(\sigma_{z}(z)\right)^{2}$
(2.11) $\displaystyle\mathcal{L}_{g}^{\tilde{\mu}_{g}}V(x,z,g)={}$
$\displaystyle\int_{\mathcal{X}}D_{g}V(x,z,g)(y)\tilde{\mu}_{g}(y,z,g)dy$
(2.12)

where $D_{g}V(x,z,g):\mathcal{X}\rightarrow\mathbb{R}$, $y\mapsto
D_{g}V(x,z,g)(y)$ is the kernel of the Riesz representation of the Frechet
derivative of $V$ with respect to the distribution at $(x,z,g)$,
$\Sigma_{x}(\cdot):=\sigma_{x}(\cdot)\left(\sigma_{x}(\cdot)\right)^{\top}$,
$a_{j}$ denotes element $j$ of vector $a$ and $A_{jk}$ denotes element $(j,k)$
of matrix $A$. The recursive optimal control, $c^{\ast}(x,z,g)$, is
characterised by the first order condition:

$0=\partial_{c}\left(u(c,z,Q(z,g))+\psi(c,x,Q(z,g))+\mathcal{L}^{c;Q}_{x}V(x,z,g)\right)\big{|}_{c=c^{*}(x,z,g)}$

_Kolmogorov Forward Equation (KFE):_ Denote the recursive equilibrium optimal
control of the individual agents by $c^{*}(x,z,g;\tilde{\mu}_{g})$. Compared
with the above notation, we add $\tilde{\mu}_{g}$ to stress the fact that the
belief of the agent may differ from the true $\mu_{g}$. In Appendix
LABEL:asec:genericModel, we show that, for a given $z$ path, the evolution of
the distribution under the optimal control can be characterized by the
Kolmogorov Forward Equation (KFE):888 Observe that there is no noise term in
the KFE because $dB_{t}^{0}$ does not directly impact the evolution of
idiosyncratic states. ,999 For notational convenience, we assume in equation
(2.13) that all variables are continuous. The same formula extends naturally
to discrete variables if the integrals in the discrete dimensions are
interpreted as sums and all derivatives in that dimension are set to zero.

$\displaystyle dg_{t}(x)={}$
$\displaystyle\mu_{g}(x,z_{t},g_{t};c^{*},\tilde{\mu}_{g})dt,\quad\text{where}$
(2.13) $\displaystyle\mu_{g}(x,z,g;c^{*},\tilde{\mu}_{g}):={}$
$\displaystyle-\sum_{j=1}^{N_{x}}\partial_{x_{j}}\left[\mu_{x,j}(c^{*}(x,z,g;\tilde{\mu}_{g}),x,z,Q(z,g))g(x)\right]$
(2.14)
$\displaystyle+\frac{1}{2}\sum_{j,k=1}^{N_{x}}\partial_{x_{j},x_{k}}\left[\Sigma_{x,jk}(c^{*}(x,z,g;\tilde{\mu}^{g}),x,z,Q(z,g))g(x)\right]$
(2.15)
$\displaystyle+\lambda(g(x-\breve{\varsigma}(x,z,g;\tilde{\mu}_{g}))|I-D_{x}\breve{\varsigma}(x,z,g;\tilde{\mu}_{g})|-g(x))$
(2.16)

where $\breve{\varsigma}$ is defined so that
$X=x-\breve{\varsigma}(x,z,g;\tilde{\mu}_{g})$ is equivalent to
$X+\varsigma(c^{*}(X,z,g;\tilde{\mu}_{g}),X,z,Q(z,g))=x$. Under this recursive
characterization, the belief consistency condition becomes that
$\mu_{g}=\tilde{\mu}_{g}$.

_Master Equation:_ We follow the approach of Lions, (2011) and characterize
the equilibrium in one PDE, which is often referred to as the “master
equation” of the “mean field game”. This formulation is particularly
convenient when the evolution of the economy is subject to aggregate shocks
and the evolution of the aggregate state variables cannot be determined
deterministically. Conceptually, the master equation is related to the HJBE
(2.7) but it imposes the belief consistency by putting the equilibrium KFE
into the HJBE. Formally, the equilibrium value function $V(x,z,g)$ is the
solution to the following “master equation”:

$\displaystyle 0=\mathcal{L}V(x,z,g)={}$ $\displaystyle-\rho
V(x,z,g)+u(c^{*}(x,z,g),z,Q(z,g))+\psi(c^{*}(x,z,g),x,Q(z,g))$ (2.17)
$\displaystyle+(\mathcal{L}_{x}^{c^{*}(x,z,g);Q}+\mathcal{L}_{z}+\mathcal{L}_{g}^{\mu_{g}})V(x,z,g)$
(2.18)

where all operators are as defined above. The mathematical challenge of
working with the master equation is that it contains an infinite dimensional
variable, $g$, and a derivative with respect to this variable. This poses a
collection of mathematical difficulties for the mean field game theory
literature, which attempts to find conditions under which the infinite
dimensional master equation is well defined and has a solution. We refer to
Cardaliaguet et al., (2019); Bensoussan et al., (2015) for more details. By
contrast, we are focused on numerical approximation, which means that we need
to find a finite dimensional approximation to the distribution, convert the
master equation into a high, but finite, dimensional PDE, and develop
techniques for solving the PDE. The goal of our paper is use deep learning
techniques to characterize such a finite dimensional numerical approximation
without using a perturbation approach.

###  2.3 Model Generality

Our model set up nests many canonical models in macroeconomics such as
continuous time versions of Krusell and Smith, (1998) and Khan and Thomas,
(2008). However, we have also made strong assumptions, which we discuss below:

  1. (i).

We have assumed that the prices, $q$, can be expressed explicitly in closed
form as functions of the aggregate state variables $(z,g)$, in equation (2.4).
In models with complicated asset pricing and market clearing conditions we
need to solve an auxiliary fixed point problem in order to solve for the
recursive representation of the price. In Gopalakrishna et al., (2024), some
of the authors to this paper extend the methodology to resolve the challenges
of introducing portfolio choice and long-term assets with complicated pricing.

  2. (ii).

We have assumed that the distribution only enters the master equation through
the pricing function. If the model had search and matching frictions, then the
distribution would impact which agents are likely to be matched and so enter
the master equation in a more complicated way. In Payne et al., (2024), some
of the authors to this paper extend the methodology to resolve these
challenges.

  3. (iii).

We have assumed that the aggregate Brownian motion $B^{0}$ does not directly
shock the evolution of idiosyncratic states. Instead it changes prices and so
indirectly changes agent controls. This means that the KFE (2.13) does not
have an aggregate noise term $dB_{t}^{0}$. The solution approach easily
generalizes to handle KFEs with aggregate noise if it does not lead to
significantly more complicated market clearing conditions (a common extension
studied is Fernández-Villaverde et al., (2023)). In Gopalakrishna et al.,
(2024), some of the authors to this paper consider such cases (and more
complicated extensions with long-term asset pricing).

Ultimately, the goal of this paper is not to offer a toolbox to solve every
macroeconomic model because deep learning techniques need to be customized to
particular master equations. Instead, we explain and compare how to use
different types of distribution approximations to solve continuous time models
with heterogeneous agents. This offers a foundational set of techniques that
are being extended in other papers.

##  3 Finite Dimensional Master Equation

In this section, we study finite dimensional approximations to the population
distribution $g$. Let $\hat{\varphi}\in\hat{\Phi}\subseteq\mathbb{R}^{N}$ be a
finite dimensional parameter vector, which could represent the collection of
agents, the bins in a histogram, or the coefficients in a projection. Let
$\hat{G}$ be a mapping from parameters to distributions:

$\displaystyle\hat{G}:\hat{\Phi}\to\mathcal{G},\quad\hat{\varphi}\mapsto\hat{G}(\hat{\varphi})=\hat{g}$
(3.1)

We look for an approximation to the value function of the form
$\hat{V}:\mathcal{X}\times\mathcal{Z}\times\hat{\Phi}\rightarrow\mathbb{R}$,
$(x,z,\hat{\varphi})\mapsto\hat{V}(x,z,\hat{\varphi})$ that satisfies an
approximate master equation:

$\displaystyle 0=\hat{\mathcal{L}}\hat{V}(x,z,\hat{\varphi}):={}$
$\displaystyle-\rho\hat{V}(x,z,\hat{\varphi})+u(\hat{c}^{*}(x,z,\hat{\varphi}),z,\hat{Q}(z,\hat{\varphi}))+\psi(\hat{c}^{*}(x,z,\hat{\varphi}),x,\hat{Q}(z,\hat{\varphi}))$
(3.2)
$\displaystyle+(\mathcal{L}_{x}^{\hat{c}^{*}(x,z,\hat{\varphi});\hat{Q}}+\mathcal{L}_{z}+\hat{\mathcal{L}}_{g})\hat{V}(x,z,\hat{\varphi}),$
(3.3)

where $\hat{Q}(z,\hat{\varphi}):=Q(z,\hat{G}(\hat{\varphi}))$ and
$\hat{c}^{*}$ is defined to satisfy:

$\displaystyle
0=\partial_{c}\left(u(c,z,\hat{Q}(z,\hat{\varphi}))+\psi(c,x,\hat{Q}(z,\hat{\varphi}))+\mathcal{L}^{c;\hat{Q}}_{x}V(x,z,\hat{\varphi})\right)\big{|}_{c=\hat{c}^{*}(x,z,\hat{\varphi})}$
(3.4)

The operators $\mathcal{L}_{x}^{c;Q}$ and $\mathcal{L}_{z}$ are the same as
before. However, the operator $\hat{\mathcal{L}}_{g}$, is different because it
depends upon the way that the distribution is approximated.

We characterize three distribution approximation approaches. The first
approach approximates the distribution with a finite collection of agents. The
second approach approximates the continuous state space by a finite collection
of grid points. The third approach projects the distribution onto a finite set
of basis functions. In each case, we show how $\hat{G}$ and
$\hat{\mathcal{L}}_{g}$ are defined. We then discuss the relationship to
perturbation approaches.

###  3.1 Finite Agent Approximation

_Distribution approximation:_ In this approach, we restrict the model so that
the economy contains a large but finite number of agents $N<\infty$. In this
case, the finite dimensional parameter vector is the vector of agent
positions: $\hat{\varphi}_{t}:=\left(x_{t}^{i}\right)_{i\leq N}$. The mapping
$\hat{G}$ takes the agents positions and computes the empirical measure:
$\hat{G}(\hat{\varphi})=\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{t}^{i}}$ where
$\delta_{x_{t}^{i}}$ denotes a Dirac mass at $x_{t}^{i}$. The evolution of
$\hat{\varphi}_{t}$ is simply the law of motion for each agents $x_{t}^{i}$,
as described by equation (2.2).

_The Operator $\hat{\mathcal{L}}_{g}$:_ The market clearing condition now
becomes $q_{t}=Q(z_{t},\hat{\varphi}_{t})=Q(z,\hat{G}(\hat{\varphi}_{t}))$, as
described in the introduction.101010 In the examples we consider, it is
straightforward to extend $Q$ to an empirical measure because $Q$ only depends
upon moments of the distribution. For more complicated cases, we could fit a
kernel to the empirical measure to approximate smooth density.  However, to
maintain the price taking assumption in the finite agent model, we impose that
agent $i$ behaves as if their individual actions do not influence prices.
Formally, this means that agent $i$ perceives the pricing function to be
$q_{t}=Q(z_{t},\hat{\varphi}_{t}^{-i})$, where
$\hat{\varphi}_{t}^{-i}=\\{x_{t}^{j}\in N^{-i}\\}$ is the position of the
other agents $N^{-i}:=\\{j\leq N:j\neq i\\}$. Ultimately, this will ensure
that the neural network trains the policy rules as if the agents believe that
their assets do not influence the market prices. The approximate distribution
impact operator $\hat{\mathcal{L}}_{g}$ is given by:

$\displaystyle\begin{aligned}
(\hat{\mathcal{L}}_{g}\hat{V})(x^{i},z,\hat{\varphi}):={}&\sum_{j\neq
i}D_{\hat{\varphi}^{j}}\hat{V}(x^{i},z,\hat{\varphi})\mu_{x}(\hat{c}^{*}(x^{j},z,\hat{\varphi}),x^{j},z,\hat{Q}(z,\hat{\varphi}^{-j}))\\\
{}&+\sum_{j\neq
i}\frac{1}{2}\text{tr}\left\\{\Sigma_{x}(x^{j},z,\hat{Q}(z,\hat{\varphi}^{-j}))D^{2}_{\hat{\varphi}^{j}}\hat{V}(x^{i},z,\hat{\varphi})\right\\}\\\
{}&+\sum_{j\neq
i}\lambda(x^{j},z)\left(\hat{V}(x^{i},z,\hat{\varphi}^{-i}+\Delta_{j})-\hat{V}(x^{i},z,\hat{\varphi}^{-i})\right),\end{aligned}$
(3.5)

where for simplicity we have used the vector and trace notation (rather than
the summation notation in Section 2) and where
$D_{\hat{\varphi}^{j}}\hat{V}(x^{i},z,\hat{\varphi})$ is the partial gradient
of $\hat{V}$ with respect to the $j$-th point in $\hat{\varphi}$ (i.e.,
$x^{j}_{t}$), and $\Delta_{k}=0$ for $k\neq j$ and
$\Delta_{k}=\varsigma_{x}(x^{k},z,\hat{Q}(z,\hat{\varphi}^{-k}))$ for $k=j$.

_Convergence properties:_ The solution $V(x^{i},z,\hat{\varphi})$ is expected
to converge to the solution of the master equation (2.18) as the number of
agents, $N$, grows to infinity so long $V$ is smooth. Intuitively, this relies
on the idiosyncratic noise in the population distribution averaging out as the
population becomes large, see e.g. Sznitman, (1991). Such results have been
extended to systems with equilibrium conditions by the mean field games
literature; see e.g. Cardaliaguet et al., (2019); Lacker, (2020); Delarue et
al., (2020).

###  3.2 Discrete State Space Approximation

We consider $N$ points in the state space, denoted by
$\xi_{1},\dots,\xi_{N}\in\mathcal{X}$. We approximate $g$ by a vector
$\hat{\varphi}=(\hat{\varphi}_{1},\dots,\hat{\varphi}_{N})\in\mathbb{R}^{N}$,
whose values represent the masses at $\xi_{1},\dots,\xi_{N}$. The mapping
$\hat{G}$ then takes the form:
$\hat{G}(\hat{\varphi})=\sum_{n=1}^{N}{\hat{\varphi}_{n}\delta_{\xi_{n}}}$,
where again $\delta_{\xi_{n}}$ denotes a Dirac mass at $\xi_{n}$.
Conceptually, the finite agent approximation fixes the mass associated to each
$x_{t}^{i}$ and allows the $x_{t}^{i}$ values to move whereas the discrete
state space approximation fixes the grid points $\xi_{n}$ and allows the
masses at each grid point to change.

We denote by $\hat{\varphi}_{t}$ the vector at time $t$. Its evolution is
given by an ordinary differential equation in dimension $N$ of the form:

$d\hat{\varphi}_{t}=\mu_{\hat{\varphi}}(z_{t},\hat{\varphi}_{t})dt$ (3.6)

describing the evolution of mass at $\xi_{1},\dots,\xi_{N}$. The right-hand
side needs to be obtained using information from the KFE (2.13). In our
numerical examples we use a finite difference approximation to the KFE to
derive $\mu_{\hat{\varphi}}$, analogous to the approximation described in
Achdou et al., 2022a . However, the technique can be applied to other types of
approximations, like the finite volume method used by Huang, 2023b . In
Appendix LABEL:asec:master_equations, we provide an explicit example for the
Krusell and Smith, (1998) model using a finite difference scheme. Here, we
describe a generic approximation of the KFE equation (2.13) using a finite
difference scheme:

$\displaystyle\mu_{\hat{\varphi}}(z,\hat{\varphi}):={}$
$\displaystyle-\sum_{j=1}^{N_{x}}\hat{\partial}_{x_{j}}\left[\bm{\mu}_{x,j}(z,\hat{\varphi})\odot\hat{\varphi}\right]+\frac{1}{2}\sum_{j,k=1}^{N_{x}}\hat{\partial}_{x_{j},x_{k}}\left[\bm{\Sigma}_{x,jk}(z,\hat{\varphi})\odot\hat{\varphi}\right]$
$\displaystyle+\lambda(\Delta_{\breve{\bm{\varsigma}}(z,\hat{\varphi})}[\hat{\varphi}]\odot\bm{a}(z,\hat{\varphi})-\hat{\varphi})$

where $\odot$ denotes the element-wise product of vectors (Hadamard product),
and $\bm{\mu}_{x,j}(z,\hat{\varphi})$, $\bm{\Sigma}_{x,jk}(z,\hat{\varphi})$,
$\breve{\bm{\varsigma}}_{j}(z,\hat{\varphi})$,
$\bm{a}(z,\hat{\varphi})\in\mathbb{R}^{N}$ and
$\breve{\bm{\varsigma}}(z,\hat{\varphi})\in(\mathbb{R}^{N_{x}})^{N}$ are the
vectors representing the values on the grid points:

$\displaystyle\bm{\mu}_{x,j}(z,\hat{\varphi})$
$\displaystyle=\left(\mu_{x,j}(\hat{c}^{*}(\xi_{n},z,\hat{\varphi}),\xi_{n},z,\hat{Q}(z,\hat{\varphi})):n=1,...,N\right)$
$\displaystyle\bm{\Sigma}_{x,jk}(z,\hat{\varphi})$
$\displaystyle=\left(\Sigma_{x,jk}(\hat{c}^{*}(\xi_{n},z,\hat{\varphi}),\xi_{n},z,\hat{Q}(z,\hat{\varphi})):n=1,...,N\right)$
$\displaystyle\breve{\bm{\varsigma}}_{j}(z,\hat{\varphi})$
$\displaystyle=(\breve{\varsigma}_{j}(\xi_{n},z,\hat{\varphi}):n=1,...,N)$
$\displaystyle\breve{\bm{\varsigma}}(z,\hat{\varphi})$
$\displaystyle=(\breve{\varsigma}(\xi_{n},z,\hat{\varphi})):n=1,...,N)$
$\displaystyle\bm{a}(z,\hat{\varphi})$
$\displaystyle=(|I-((\hat{\partial}_{x_{k}}[\breve{\bm{\varsigma}}_{j}(z,\hat{\varphi})])_{n})_{j,k=1,...,N_{x}}|:n=1,...,N).$

Informally, $\hat{\partial}_{x_{j}}$ and $\hat{\partial}_{x_{j},x_{k}}$ takes
the function values on the grid points and return finite difference
approximations of the first derivative w.r.t. $x_{j}$ and the second
derivative w.r.t. $(x_{j},x_{k})$. Informally,
$\Delta_{\breve{\bm{\varsigma}}}$ takes the function values on the grid points
and returns interpolated values at shifted inputs
$\xi_{n}-\breve{\bm{\varsigma}}_{n}$, which may not be on the grid. More
formally, $\hat{\partial}_{x_{j}}$, $\hat{\partial}_{x_{j},x_{k}}$, and
$\Delta_{\breve{\bm{\varsigma}}}$ are operators
$\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ that map vectors
$\bm{f}\in\mathbb{R}^{N}$ of function values of a function $f$ on the grid
(i.e., $\bm{f}_{n}=f(\xi_{n})$) into vectors of the same size, such that,
$(\hat{\partial}_{x_{j}}[\bm{f}])_{n}\approx\partial_{x_{j}}f(\xi_{n})$,
$(\hat{\partial}_{x_{j},x_{k}}[\bm{f}])_{n}\approx\partial_{x_{j},x_{k}}f(\xi_{n})$,
and $(\Delta_{\breve{\bm{\varsigma}}}[\bm{f}])_{n}\approx
f(\xi_{n}-\breve{\bm{\varsigma}}_{n})$.111111 The precise forms of
$\hat{\partial}_{x_{j}}$, $\hat{\partial}_{x_{j},x_{k}}$, and
$\Delta_{\breve{\bm{\varsigma}}}$ are implementation-specific. In the most
common arrangement, grid points are positioned on a rectangular grid, the
operators $\hat{\partial}_{x_{j}}$ and $\hat{\partial}_{x_{j},x_{k}}$ compute
approximate derivatives at each point based on differences in function values
to neighbor grid points, and the operator $\Delta_{\breve{\bm{\varsigma}}}$
implements linear interpolation on the regular grid.

_The Operator $\hat{\mathcal{L}}_{g}$:_ The approximate distribution impact
operator $\hat{\mathcal{L}}_{g}$ is defined by:

$(\hat{\mathcal{L}}_{g}\hat{V})(x,z,\hat{\varphi})=\sum_{n=1}^{N}\mu_{\hat{\varphi},n}(z,\hat{\varphi})\,\partial_{\hat{\varphi}_{n}}\hat{V}(x,z,\hat{\varphi}),$
(3.7)

where $\mu_{\hat{\varphi},n}(z,\hat{\varphi})$ denotes the $n$-th coordinate
of the $N$-dimensional vector $\mu_{\hat{\varphi}}(z,\hat{\varphi})$.

_Convergence properties:_ Once again, we expect the solution of the
approximate master equation $V(x,z,\hat{g})$ to converge towards the solution
of the true master equation (2.18) so long as $V$ is smooth. Convergence of
mean field games with discrete state spaces to mean field games with
continuous state spaces has been proved by Bayraktar et al., (2018);
Hadikhanloo and Silva, (2019) without noise or with idiosyncratic noise, and
by Bertucci and Cecchin, (2022) with common noise.

###  3.3 Projection onto Basis

_Distribution approximation:_ In this approach, we represent the distribution
$g_{t}$ by functions of the form:

$\hat{g}_{t}(x)=\hat{G}(\hat{\varphi}_{t})(x):=b_{0}(x)+\sum_{n=1}^{N}{\hat{\varphi}_{n,t}b_{n}(x)},$

where
$\hat{\varphi}_{t}:=(\hat{\varphi}_{1,t},...,\hat{\varphi}_{N,t})\in\hat{\Phi}:=\mathbb{R}^{N}$
is a vector of real-valued coefficients and $b_{0},b_{1},...,b_{N}$ is a
collection of linearly independent real-valued functions on $\mathcal{X}$ that
satisfy

$\int_{\mathcal{X}}b_{n}(x)dx=\begin{cases}1,&n=0\\\ 0,&n\geq 1\end{cases}$
(3.8)

and we refer to as the basis of the projection. To complete this approximation
description we need to specify the the law of motion for $\hat{\varphi}_{t}$
and the choice of basis.

As in the discrete state space approach, the evolution of $\hat{\varphi}_{t}$
is given by an ODE of the form (3.6) with drift $\mu_{\hat{\varphi}}$. To
derive $\mu_{\hat{\varphi}}$, we would like to use the KFE (2.13) directly.
However, if we substitute our projection $\hat{g}_{t}$ into the KFE, then we a
get a linear equation system:

$\displaystyle\sum_{n=1}^{N}{\mu_{\hat{\varphi},n,t}b_{n}(x)}=\mu^{g}(\hat{c}_{t}^{*},x,z_{t},\hat{g}_{t})$
(3.9)

for $\mu_{\hat{\varphi}}$ that has generically no solution for $N<\infty$
because we typically have an infinite number of equations (for each
$x\in\mathcal{X}$) but only $N$ degrees of freedom to solve these equations.
To make progress, we need to generate a finite set of equations to solve. One
option would be to discretize the $x$-dimension, as in the previous technique.
However, we instead find it helpful to work with a more general discretization
to a set of $M\geq N$ “test functions”
$\phi:\mathcal{X}\rightarrow\mathbb{R}$. This is useful because it ultimately
allows us to choose test functions that focus the approximation accuracy on
pre-selected statistics of the distribution that are likely to be economically
important. For example, if prices $q_{t}=Q(z_{t},g_{t})$ depend only on
certain moments of the distribution such as the mean, it makes sense to choose
test functions that represent these moments. To make this idea precise, we
start with the integral form of the KFE (2.13) before integration by parts:

$\displaystyle\frac{d\int_{\mathcal{X}}{\phi(x)g_{t}(x)dx}}{dt}$
$\displaystyle=\int_{\mathcal{X}}{\mathcal{L}_{x}^{c^{\ast}(x,z_{t},g_{t});Q}\phi(x,z_{t},g_{t})g_{t}(x)dx}=:\mu_{\phi}(z_{t},g_{t};c^{\ast})$

that has to hold for all the $M$ test functions. This variant of the KFE
describes the time evolution of the statistics
$\int_{\mathcal{X}}{\phi(x)g_{t}(x)dx}$ of the distribution.121212We recover
the original KFE (2.13) for the Dirac delta distributions as test functions,
i.e. $\phi=\delta_{x}$ for $x\in\mathcal{X}$. We define the coefficient drifts
$\mu_{\hat{\varphi},n}(z,\hat{\varphi})$ (conditional on the aggregate state
$(z,\hat{\varphi})$) as the regression coefficients that minimize, in the
least-squares sense, the $M$ linear regression residuals131313In practice, the
integral terms appearing in the residual formula have to be computed either
analytically (if possible) or by numerical quadrature.

$\varepsilon_{m}(z,\hat{\varphi}):=\mu_{\phi_{m}}(z,\hat{G}(\hat{\varphi});\hat{c}^{\ast})-\sum_{n=1}^{N}{\mu_{\hat{\varphi},n}(z,\hat{\varphi})\int_{X}{\phi_{m}(x)b_{n}(x)dx}},\qquad
m=1,...,M.$

Ultimately, this means that our approximation is most accurate the on chosen
statistics encoded in the test functions.

The approximation presented so far works, in principle, for any choice of
basis. Here we propose a basis that approximately tracks the persistent
dimensions of $g_{t}$ while neglecting those dimensions that mean-revert fast.
These persistent dimensions of the distribution are related to certain
eigenfunctions of the differential operator characterizing the KFE (2.13).
Because this differential operator is generally time-dependent and stochastic,
we first replace it by a time-invariant operator
$\overline{\mathcal{L}}^{KF}$. For example, this could be the steady-state
operator $\mathcal{L}^{KF,ss}$ that is defined as the KFE operator at the
steady state, $g_{t}=g^{ss}$ in a simplified model without common noise.141414
Our chosen basis approximately tracks the persistent dimensions of $g_{t}$ if
$\tilde{\mathcal{L}}^{KF}$ is similar to full stochastic operator
$\mathcal{L}^{KF}_{t}$. This is plausible for
$\tilde{\mathcal{L}}^{KF}=\mathcal{L}^{KF,ss}$ because both operators share
many similar features.  Let $\\{b_{i}\,:\,i\geq 0\\}$ be the set of
eigenfunctions of $\overline{\mathcal{L}}^{KF}$ with corresponding eigenvalues
$\\{\lambda_{i}\in\mathbb{C}\,:\,i\geq 0\\}$. If the dynamics prescribed by
the KFE are locally stable around the steady state $g^{ss}$, there is one
eigenvalue $\lambda_{0}=0$ with eigenfunction $b_{0}=g^{ss}$ and all remaining
eigenvalues have negative real part, $\Re\lambda_{i}<0$. We pick as our basis
$b_{0},b_{1},...,b_{N}$ the $N+1$ eigenfunctions corresponding to eigenvalues
with real parts closest to zero as these represent the most persistent
components. We provide further details on this basis choice in Appendix
LABEL:sec:eigenfunction-basis.

_The Operator $\hat{\mathcal{L}}_{g}$:_ The approximate distribution impact
operator $\hat{\mathcal{L}}_{g}$ is defined by:

$(\hat{\mathcal{L}}_{g}\hat{V})(x,z,\hat{\varphi})=\sum_{n=1}^{N}\mu_{\hat{\varphi},n}(z,\hat{\varphi})\,\partial_{\hat{\varphi}_{n}}\hat{V}(x,z,\hat{\varphi}).$

_Convergence:_ There are fewer results about convergence for projections in
the continuous time mean-field-game literature. However, in discrete time,
Prohl, (2017) has proven convergence results for particular projections.

###  3.4 Relationship to Continuous Time Perturbation Approaches

All approaches offer a global characterization of the aggregate dynamics. This
opens up the possibility of studying the full nonlinear dynamics but requires
non-trivial approximations to the KFE to preserve the nonlinearity. Here, we
briefly contrast our approaches to continuous time perturbation approaches
that replace the nonlinear evolution of aggregate states by a simpler linear
or quadratic one. We compare to state space perturbations rather than sequence
space perturbations since they are closer to our work.

Ahn et al., (2018) adapts and extends the perturbation approach of Reiter,
(2002, 2009) to continuous-time settings. They first discretize the state
space and then linearize the model with respect to aggregate dynamics. This
allows the authors to solve for equilibrium using conventional matrix algebra
techniques. In our setup, this is analogous to linearizing the finite-
dimensional master equation that results from our discrete state space
approximation and imposing their dimension reduction.

The perturbation method of Bilal, (2023) works directly with the master
equation formulation. Relative to traditional techniques, they take linear or
quadratic perturbations of the master equation with respect to the aggregate
states $(z,g)$. The approximation by a linear or quadratic functional form
reduces the complexity of the problem because one only needs to solve for the
perturbation coefficients not for the value function on the infinite-
dimensional state space of the full nonlinear problem.151515 Note that, when
solving the remaining equation numerically, Bilal, (2023) chooses a finite-
difference method, which still requires a discretization of the state space
akin to our discrete state space approach.  By contrast we obtain a finite-
dimensional master equation by means of approximating the distribution, which
allows us to preserve the full nonlinearity of the problem.

Alvarez and Lippi, (2022) use Fourier methods to analytically characterize
impulse responses functions to shocks to the cross-sectional distribution
$g_{t}$ in a model that features a linear KFE and a constant
$z_{t}$-state.161616Alvarez et al., (2023) employ similar methods in a
nonlinear model but linearize dynamics first. That approach is closely related
to our projection method based on an eigenfunction basis. Specifically, if the
KFE is truly linear, then the following proposition holds:

######  Proposition 1.

If the KFE (2.13) is of the form $dg_{t}=\mathcal{L}^{KF}g_{t}dt$ with a
(constant) linear operator $\mathcal{L}^{KF}$, $b_{0}=g^{ss}$, $b_{1}$, …,
$b_{N}$ are eigenfunctions of $\mathcal{L}^{KF}$ with eigenvalues
$\lambda_{0}=0$, $\lambda_{1}$, …, $\lambda_{N}$, and
$g_{0}-g^{ss}=\sum_{n=1}^{N}\varphi_{i,0}b_{i}\in\operatorname{span}\\{b_{1},...,b_{N}\\}$,
then $g_{t}$ satisfies equation (2.13) for all $t\geq 0$ if and only if, for
all $t\geq 0$,

$g_{t}=g_{t}-g^{ss}=\sum_{i=1}^{N}{\varphi_{i,0}e^{\lambda_{i}t}b_{i}}.$
(3.10)

###### Proof.

See Appendix LABEL:sec:eigenfunction-basis. ∎

This says that, if the KFE is linear and we use
$\overline{\mathcal{L}}^{KF}=\mathcal{L}^{KF}$ to form the eigenfunction
basis, then we are in the non-generic case in which we can approximate the KFE
perfectly because the equation system (3.9) has a solution even though there
are only finite degrees of freedom. The connection to the impulse response
formula provided by Alvarez and Lippi, (2022) is as follows. If we consider an
initial shock that moves $g_{0}$ away from $g^{ss}$ in a way such that
$g_{0}-g^{ss}$ is spanned by $b_{1},...,b_{N}$, then equation (3.10) tells us
exactly how the distribution evolves over time. We can further integrate this
equation with respect to any function $f:\mathcal{X}\rightarrow\mathbb{R}$ to
obtain a formula for the dynamics of the statistic $\int{f(x)g_{t}(x)dx}$ over
time that resembles the impulse response representation of Alvarez and Lippi,
(2022).171717A difference to Alvarez and Lippi, (2022) is that, in the context
of their model, they know closed-form solutions for the eigenfunctions and
they consider $N=\infty$, so that they can consider any initial distribution
displacement $g_{0}-g^{ss}$.

Our computational approach is different from Proposition 1 in that we do not
presume a linear KFE and so allow for non-linear distributional dynamics. In
this case, the KFE can no longer be approximated perfectly on a finite basis
and so there are no closed-form expressions in the spirit of equation (3.10).
Instead, we need to make additional choices about how to approximate the KFE,
which leads to the more complicated procedure outlined in Section 3.3.

###  3.5 Comparison Between Our Approximation Approaches

Table 1 summarizes the key differences between the distribution
approximations: a finite population, a discrete state space, and a projection
onto a finite set of basis functions.

|  Finite Population |  Discrete State |  Projection
---|---|---|---
Distribution approx. |  $\frac{1}{N}\sum_{i=1}^{N}\delta_{x_{t}^{i}}$ |  $\sum_{n=1}^{N}{\hat{\varphi}_{n,t}\delta_{\xi_{n}}}$ |  $\sum_{n=0}^{N}{\hat{\varphi}_{n,t}b_{n}(x)}$
KFE approx. |  Evolution of other agents’ states |  Evolution of mass between discrete states |  Evolution of projection coefficients
Table 1: Comparison of Distribution Approximations

We discuss how these approximations compare with regard to computational
difficulties:

  1. 1.

Dimensionality ($N$): The approximation dimension needs to be large enough to
capture sufficient shape in the distribution. The projection method can
potentially have the lowest dimension if the choice of basis is efficient
(e.g. $N=5$ in our Krusell and Smith, (1998) example). The finite population
needs to be large enough to average out idiosyncratic noise (e.g. $N=40$). The
discrete state space needs to be sufficiently fine to approximate the
derivatives in the KFE, which means it needs a high dimension (e.g. $N=200$).

  2. 2.

Customization: For the finite agent approach, we just choose the number of
agents, $N$. For the discrete state space approach, we choose a grid and a
method for approximating the KFE on the grid points. For the projection
method, we choose a set of basis function and a set of statistics on which to
minimize the error. In this sense, the projection is potentially lower
dimensional because we need to make more intricate choices in the setup.

  3. 3.

Computational “bottlenecks”: Each method has its own computational
bottlenecks. For the finite agent method, we need to switch agent positions in
the neural network approximation to $V$ when we calculate the derivatives of
$V$ with respect to the other agent positions in equation (3.5). For the
discrete state space method, the dimensionality of the approximation is the
main computational problem. For the projection method, determining the drift
$\mu_{\hat{\varphi}}$ of the distribution approximation is computationally
involved as we need to solve a linear regression problem and compute several
integrals by quadrature for every single evaluation of the master equation.

##  4 Solution Method

All approaches in Section 3 lead to finite approximations to the density,
$\hat{g}$, and the master equation operator $\hat{\mathcal{L}}$. However, the
resulting master equations are high dimensional and so cannot be solved by
traditional numerical techniques. Instead, we represent the solution to the
approximate master equation by a neural network and deploy tools from the
“deep learning” literature to “train” the neural network to solve the
approximate master equation.

###  4.1 Neural Network Approximations

A neural network is a type of parametric functional approximation that is
built by composing affine and non-linear functions in a chain or “network”
structure (see Goodfellow et al., (2016)). We let
$\hat{X}:=\\{x,z,\hat{\varphi}\\}$ denote the collection of inputs into the
neural network representation of the value function. We denote the neural
network approximation to the value function by
$\hat{V}(\hat{X})\approx\mathbb{V}(\hat{X};\Theta)$, where $\Theta$ are the
neural network parameters that depend upon the architecture. There are many
types of neural networks. The simplest form is a “feedforward” neural network
which is defined by:

$\displaystyle\begin{aligned}
h^{(1)}={}&\phi^{(1)}(W^{(1)}\hat{X}+b^{(1)})&&&\ldots\text{Hidden layer 1}\\\
\vdots\\\ h^{(H)}={}&\phi^{(H)}(W^{(H)}h^{(H-1)}+b^{(H)})&&&\ldots\text{Hidden
layer $H$}\\\ o={}&W^{(H+1)}h^{(H)}+b^{(H+1)}&&&\ldots\text{Output layer}\\\
\mathbb{V}={}&\phi^{(H+1)}(o)&&&\ldots\text{Output}\end{aligned}$ (4.1)

where the $\\{h^{(i)}\\}_{i\leq H}$ are vectors referred to as “hidden layers”
in the neural network, $\\{W^{(i)}\\}_{i\leq(H+1)}$ are matrices referred to
as the “weights” in each layer, $\\{b^{(i)}\\}_{i\leq(H+1)}$ are vectors
referred to as the “biases” in each layer, $\\{\phi^{(i)}\\}_{i\leq(H+1)}$ are
non-linear functions applied element-wise to each affine transformation and
referred to as “activation functions” for each layer. The length of hidden
layer, $h^{(i)}$, is defined as the number of neurons in the layer, which we
refer to as $\\#h^{(i)}$. The total collection of parameters is denoted by
$\Theta=\\{W^{(i)},b^{(i)}\\}_{i\leq(H+1)}$. The goal of deep learning is to
train the parameters, $\Theta$, to make $\mathbb{V}(\cdot;\Theta)$ a close
approximation to $\hat{V}$.

The neural network defined in (4.1) is called a “feedforward” network because
hidden layer $i$ cannot depend on hidden layers $j>i$. This is in contrast to
a “recurrent” neural network where any hidden layer can be a function of any
other hidden layer. It is called “fully connected” if all the entries in the
weight matrices can be non-zero so each layer can use all the entries in the
previous layer. In this paper, we will consider a fully connected
“feedforward” network to be the default network. This is because these
networks are the quickest to train and so we typically start by trying out
this approach. However, there are applications where we find that more
complicated neural network architectures are useful. In particular, we find
that the type of recurrent neural network suggested by the Deep Galerkin
Method in Sirignano and Spiliopoulos, (2018) is helpful for discrete state and
projection approximations.

###  4.2 Solution Algorithm

We train the neural network to learn parameters $\Theta$ that minimize the
error in the master equation and boundary conditions. We describe the key
steps in Algorithm 1.181818 The generic pseudo-code given in Algorithm 1 can
be modified in practice. For example, instead of fixing a precision threshold,
one can fix a number of iterations, and instead of using a fixed sequence of
learning rates, one can use an adaptive method, such as Adam. Essentially, the
algorithm samples random points in the discretized state space
$\\{x,z,\hat{\varphi}\\}$, calculates the master equation error on those
points under the current neural network approximation, and then updates the
neural network parameters to decrease the master equation error. In practice,
our loss consists of two terms: $\mathcal{E}^{e}$ for the average mean squared
master equation error and $\mathcal{E}^{s}$, which is used to incorporate
information about the shape of the solution (e.g., monotonicity or concavity).
Specific examples of $\mathcal{E}^{s}$ will be discussed in the following
sections. In the deep learning literature, this approach is sometimes referred
to as “unsupervised” learning (e.g. Azinovic et al., (2022)) because we do not
have direct observations of the value function, $V(x,z,\hat{g})$, and instead
have to learn it indirectly via the master equation.

Input : Initial neural network parameters $\Theta^{0}$, number of sample
points $M$, positive weights $\kappa^{e}$ and $\kappa^{s}$ on the master
equation errors; sequence of learning rates $\\{\alpha_{n}:n\geq 0\\}$,
precision threshold $\epsilon$

Output : A neural network approximation
$(x,z,\hat{\varphi})\mapsto\mathbb{V}(x,z,\hat{\varphi};\Theta)$ of the value
function.

1: Initialize neural network object $\mathbb{V}(x,z,\hat{\varphi};\Theta^{0})$
with parameters $\Theta^{0}$.

2: while Loss $>$ $\epsilon$, in iteration $n$ do

3: Generate $M$ new sample points,
$S^{n}=\\{(x_{m},z_{m},\hat{\varphi}_{m})\\}_{m\leq M}$.

4: Calculate the weighted average error:
$\displaystyle\mathcal{E}(\Theta^{n},S^{n})=\kappa^{e}\mathcal{E}^{e}(\Theta^{n},S^{n})+\kappa^{s}\mathcal{E}^{s}(\Theta^{n},S^{n})$
(4.2) where $\mathcal{E}^{e}$ is the average mean square master equation
error:  $\displaystyle\mathcal{E}^{e}(\Theta^{n},S^{n}):={}$
$\displaystyle\frac{1}{|S^{n}|}\sum_{(x,z,\hat{\varphi})\in
S^{n}}|(\hat{\mathcal{L}}\mathbb{V}(x,z,\hat{\varphi};\Theta^{n}))|^{2}$ (4.3)
in which the derivatives in the operator $\hat{\mathcal{L}}$ are calculated
using automatic differentiation and $\mathcal{E}^{s}$ is a penalty for a
“wrong” shape that depends upon the specific problem.

5: Update the parameters using stochastic gradient descent:
$\displaystyle\Theta^{n+1}=\Theta^{n}-\alpha_{n}D_{\Theta}\mathcal{E}(\Theta^{n},S^{n})$
where $D_{\Theta}\mathcal{E}$ is the gradient (i.e., vector differential)
operator.

6: end while

Algorithm 1 Pseudo Code for Generic Solution Algorithm

Although the algorithm is straightforward to describe at the high level,
implementing the deep learning training scheme successfully involves a lot of
complicated decisions. We discuss these decisions in detail for solving the
Krusell and Smith, (1998) in Sections 5.2 and
LABEL:asec:implementation_details. Here we discuss sampling, which an
implementation detail that is particular challenging, and shape constraints.

####  4.2.1 Sampling

A very important aspect of our solution algorithm is the approach used to
sample the set of training points $S^{n}$ in each iteration. Sampling
ultimately has to be tailored to the specific application at hand. This is a
difference between deep learning for continuous time and discrete time
techniques. Discrete time models need to calculate expectations and so
typically need to use simulation to approximate the expectation operator.
Continuous time models replace the expectation term in the Bellman equation by
the derivative terms in the HJBE and then sample points on which to evaluate
the HJBE. This gives continuous time techniques more flexibility in how to
sample but can also make the sampling task harder. The following general
considerations are relevant when deciding on how to sample.

We can sample the idiosyncratic state $x$, the aggregate exogenous state $z$,
and the distribution state $\hat{\varphi}$ independently with different
approaches. The separation between $x$ and $\hat{\varphi}$ deserves particular
emphasis in the context of the finite agent approximation, for which also
$\hat{\varphi}$ contains a sample of $x$-points (for the other agents). This
is because we can focus the sampling on regions of the idiosyncratic state
space with high curvature without having to increase the number of agents so
that simulations have many agents in those regions.

Sampling $x$ and $z$ is less complicated because their dimension is usually
relatively low. For example, in macro models a typical dimension of $x$ is 2-3
and a typical dimension of $z$ is 1-5. For these variables, we typically
sample from a pre-specified statistical distribution such as the uniform or
normal. The sampling can be refined by using “active sampling” (see, e.g.,
Gopalakrishna, (2021) and Lu et al., (2021)), which adapts the sampling during
training to actively learn in regions where the algorithm is having trouble
minimizing the loss function. This is achieved by regularly inspecting the
losses during training and adding training points to regions with the largest
losses.

Sampling the distribution approximation $\hat{\varphi}$ is significantly more
complicated. This is because it is typically high-dimensional and so we can
only train the neural network on a very small subset of the total possible
distributions. In this sense, deep learning does not break the “curse of
dimensionality”. Instead, it gives flexibility to train on a useful subspace
that gives enough information to the value function for economically relevant
distributions. This means that choosing the right subspace on which to sample
is very important for the algorithm to converge in reasonable time (or at
all). Ultimately, this requires us to use some information about the model
solution in the sampling. We have focused on the following three sampling
schemes:

  1. (i)

_Moment sampling:_ We first draw samples for selected moments of the
distribution that are important for calculating prices
$\hat{Q}(z,\hat{\varphi})$. We then sample $\hat{\varphi}$ from a distribution
that satisfies the moments drawn in the first step. E.g., in many
macroeconomic models, the distribution mean is particularly important for
calculating prices. In this case, we would sample from the mean and then draw
$\hat{\varphi}$ from a distribution with that mean.191919 This final step
could involve sampling from a pre-specified distribution (e.g. uniform) or
from an ergodic distribution, as described in (iii).

  2. (ii)

_Mixed steady state sampling:_ We first solve for the steady state for a
collection of fixed aggregate states $z$. This needs to be done only once
before training begins. We then draw in each training step random mixtures of
this collection of steady state distributions and then perturb with additional
noise.202020Without these perturbations, the random mixtures remain strictly
confined to a subspace whose dimension (the number of steady states in the
collection) is typically much smaller than $\dim{\hat{\Phi}}$.

  3. (iii)

_Ergodic sampling:_ We adapt the training sample dynamically by regularly
simulating the model economy based on the candidate solution for the value
function from a previous iteration.

Two additional issues arise in sampling schemes that adapt the training sample
dynamically, such as active sampling (for $x$ and $z$) and ergodic sampling
(for $\hat{\varphi}$). First, these schemes only adapt the sampling
distribution in a meaningful way if the current guess for the value function
is sufficiently good. It is therefore advisable to start with a pre-specified
sampling distribution in early training and switch to a dynamic sampling
scheme later on. Second, dynamically adapting the sample might lead to
instability of training due to feedback effects between the training sample
and the trained solution. We have found this issue to be particularly relevant
for ergodic sampling. To mitigate the issue, we have found the following work
well: (i) to use ergodic sampling only for a fraction of the training sample
and (ii) to update training points by simulating over a small time interval
starting from the end the last simulation rather than repeatedly taking very
long simulations to approximate the ergodic distribution.212121 In this sense,
despite the name “ergodic sampling”, we do not insist on drawing training
points from the ergodic distribution implied by the current value function
candidate. Instead, the ergodic distribution is reached only gradually once
the value function is close to convergence.

At a high level, we find the following sampling strategies are useful for the
different types of distribution approximations. For the finite agent
approximation, we find moment sampling to be simple and effective because the
$\hat{\varphi}$ variables have a natural interpretation as the idiosyncratic
($x$) states of the other agents in the population and so it is
straightforward to determine a region of “typical” values to sample from. For
the discrete state space approximation, we find the most stable approach is to
start with mixed steady state sampling and move to ergodic sampling once the
neural network starts to converge. For the projection approximation, we find
that a combination of moment sampling and ergodic sampling is effective. To
put moment sampling to work with projections requires a rotation of the basis
functions so as to isolate the components that correspond to the selected
moments, see Appendix LABEL:asec:master_equations for details. We discuss all
three sampling strategies in detail for our Krusell and Smith, (1998) example
in Appendix LABEL:asec:implementation_details:sampling.

####  4.2.2 Constraining the Value Function Shape

We find that the neural network can converge to “cheat solutions” that
approximately solve the differential equation by setting derivatives to zero.
One way of helping with this is to include terms in the loss function the
penalize undesirable curvature of the value function to enforce the correct
shape (e.g. penalizing non-monotonicity or non-concavity). There are other,
complementary potential ways of addressing this issue: (i) initializing the
neural network to match an initial guess satisfying known properties of the
value function so that when the neural network is trained to minimize the PDE
loss, it converges to a (local) minimizer which has the same desired
properties, (ii) sampling from a sufficiently large part of the state space
that the neural network realizes there is curvature in all dimensions, and
(iii) choosing an architecture which satisfies the constraints (or at least
makes it easier for the neural network to satisfy these constraints). We find
a combination of these approaches to be helpful.

####  4.2.3 Extension: Boundary Conditions

Algorithm 1 is suitable for solving a problem that does not include boundary
conditions. This is sufficient for our purposes because we “soften” any hard
constraints in the problem by replacing them with a utility flow penalty (see
Section 2.2). In principle, the solution algorithm could be extended to a
problem that does require boundary conditions. In this case, we would need to
sample a second sample $S^{b}$ of training points on the boundary in step 3
and, in step 4, add an additional error term $\mathcal{E}^{b}$ (with
corresponding weight $\kappa^{b}$) for the mean-squared error of boundary
condition residuals evaluated at the training set $S^{b}$, as described by
Sirignano and Spiliopoulos, (2018). All other aspects of the algorithm would
remain unchanged. That being said, while conceptually straightforward, we have
found in the context of our examples that the inequality boundary conditions
arising from hard constraints represent a significant difficulty for the
robustness of our solution algorithm. Specifically, we observe that the neural
network only learns an accurate solution if the weights on equation residuals
($\kappa^{e}$) and the boundary condition ($\kappa^{b}$) in the loss function
are well-calibrated, which we find difficult to achieve for inequality
constraints. Replacing the hard constraints by a soft penalty makes the method
much more robust.

####  4.2.4 Extension: Staggered updating

In Algorithm 1, when evaluating the master equation residuals to determine
$\mathcal{E}^{e}$, we compute the policy $\hat{c}^{\ast}$ as a function of
$\mathbb{V}(\cdot;\Theta)$ according to equation (3.4). In principle, there is
no need to parameterize the policy $\hat{c}^{\ast}$ with a separate neural
network. However, in practice there are reasons for doing so, as suggested by
Duarte et al., (2024) and Al-Aradi et al., (2022). First, if equation (3.4)
does not have a closed form solution, the computational cost of solving this
equation numerically at each training point may be prohibitively high. In this
case, we could include a separate neural network for $\hat{c}^{\ast}$ with
separate parameters $\Theta^{c}$ and use this network in place of the true
solution to equation (3.4) to evaluate the master equation errors. This would
necessitate adding steps into Algorithm 1 that involve updating the
$\hat{c}^{\ast}$-network for a given $\mathbb{V}$-network. Second, even if
equation (3.4) can be solved, adding a separate neural network for
$\hat{c}^{\ast}$ allows us to slow down the updating of the policy function
akin to a “Howard improvement algorithm”. To do so, we can make the updating
of $\Theta^{c}$ infrequent, effectively fixing the same policy rule
$\hat{c}^{\ast}$ for several iterations in the training of $\Theta$. In some
implementations for our numerical example we use this variant of our baseline
algorithm. We have found that this can help with stability, particularly when
starting from a poor (e.g. random) initial guess.

##  5 Example: Uninsurable Income Risk and TFP Shocks

A canonical macroeconomic model with heterogeneous agents and aggregate risk
is Krusell and Smith, (1998), which we refer to as the KS model. In this
section, we illustrate how our three solution approaches can be deployed to
solve the continuous time version of the model described in Achdou et al.,
2022a  and Ahn et al., (2018).

###  5.1 Model Specification

_Setting:_ There is a perishable consumption good and a durable capital stock
with depreciation rate $\delta$. The economy consists of a unit continuum
$I=[0,1]$ of households and a representative firm. The representative firm
controls the production technology, which produces consumption goods according
to the production function $Y_{t}=e^{z_{t}}K_{t}^{\alpha}L_{t}^{1-\alpha}$,
where $K_{t}$ is the capital rented at time $t$, $L_{t}$ is the labor hired at
time $t$, and $z_{t}$ is the aggregate productivity, which follows an
Ornstein-Uhlenbeck (OU) process $dz_{t}=\eta(\overline{z}-z_{t})dt+\sigma
dB_{t}^{0}$.

_Heterogeneous households:_ Each household $i\in[0,1]$ has discount rate
$\rho$ and gets flow utility $u(c_{t}^{i})=(c_{t}^{i})^{1-\gamma}/(1-\gamma)$
from consuming $c_{t}^{i}$ consumption goods at time $t$. Each household has
two idiosyncratic states $x_{t}^{i}=[a_{t}^{i},l_{t}^{i}]$, where $a_{t}^{i}$
is the household’s net wealth and $l_{t}^{i}\in\\{l_{1},l_{2}\\}$ is the
household’s labor endowment, where $l_{1}<l_{2}$ so $l_{1}$ is interpreted as
unemployment and $l_{2}$ is interpreted as employment. Labor endowments switch
idiosyncratically between $l_{1}$ and $l_{2}$ at Poisson rate
$\lambda(l_{t}^{i})$. Households choose consumption $c_{t}^{i}$ and their
idiosyncratic state evolves according to:

$\displaystyle dx_{t}^{i}={}$ $\displaystyle d\begin{bmatrix}a_{t}^{i}\\\
l_{t}^{i}\end{bmatrix}=\begin{bmatrix}s(a_{t}^{i},l_{t}^{i},c_{t}^{i},r_{t},w_{t})\\\
0\end{bmatrix}dt+\begin{bmatrix}0\\\
\check{l}_{t}^{i}-l_{t}^{i}\end{bmatrix}dJ_{t}^{i}$ (5.1)

where $r_{t}$ is the return on household wealth, $w_{t}$ is the wage rate,
$\check{l}_{t}^{i}$ is the complement of $l_{t}^{i}$, $J_{t}^{i}$ is an
idiosyncratic Poisson process with arrival rate $\lambda(l_{t}^{i})$, and the
agent’s saving function is given by $s(a,l,c,r,w)=wl+ra-c$. So, $g_{t}$
denotes the population density across $\\{a_{t}^{i},l_{t}^{i}\\}$ at time $t$,
given a filtration $\mathcal{F}_{t}^{0}$ generated by the sequence of
aggregate productivity shocks.

_Assets, markets, and financial frictions:_ Each period, there are competitive
markets for goods, capital rental, and labor. We use goods as the numeraire.
We let $r_{t}$ denote the rental rate on capital, $w_{t}$ denote the wage rate
on labor, and $q_{t}=[r_{t},w_{t}]$ denote the price vector. Given $g_{t}$ and
$z_{t}$, firm optimization and market clearing imply that $r_{t}$ and $w_{t}$
solve:

$\displaystyle r_{t}{}$ $\displaystyle=\alpha
e^{z_{t}}K_{t}^{\alpha-1}L^{1-\alpha}-\delta,$ $\displaystyle w_{t}{}$
$\displaystyle=(1-\alpha)e^{z_{t}}K_{t}^{\alpha}L^{-\alpha},$ (5.2)
$\displaystyle K_{t}{}$
$\displaystyle=\sum_{j\in\\{1,2\\}}\int_{\mathbb{R}}ag_{t}(a,l_{j})da$
$\displaystyle L{}$
$\displaystyle=\sum_{j\in\\{1,2\\}}\int_{\mathbb{R}}l_{j}g_{t}(a,l_{j})da.$
(5.3)

where we assume that economy starts with the steady state labor distribution.
So, in the terminology of Section 2, we can write the prices explicitly as
functions of $(g_{t},z_{t})$:

$\displaystyle q_{t}=\begin{bmatrix}r_{t}\\\
w_{t}\end{bmatrix}=\begin{bmatrix}\alpha
e^{z_{t}}\left(\sum_{j\in\\{1,2\\}}\int_{\mathbb{R}}ag_{t}(a,l_{j})da\right)^{\alpha-1}L^{1-\alpha}-\delta\\\
(1-\alpha)e^{z_{t}}\left(\sum_{j\in\\{1,2\\}}\int_{\mathbb{R}}ag_{t}(a,l_{j})da\right)^{\alpha}L^{-\alpha}\end{bmatrix}=:Q(g_{t},z_{t})$
(5.4)

Asset markets are incomplete so households cannot insure their idiosyncratic
labor shocks. Instead, households can trade claims to the aggregate capital
stock in a competitive asset market. The original Krusell and Smith, (1998)
model imposes the “borrowing constraint” that each agent’s net asset position,
$a_{t}^{i}$, must satisfy $a_{t}^{i}\geq\underline{a}$, where $\underline{a}$
is an exogenous “borrowing limit”. This generates an inequality boundary
constraint and a mass point, as discussed in Achdou et al., 2022a . However,
as discussed in Section 2.2, to help the neural network train more smoothly,
we instead follow Brzoza-Brzezina et al., (2015) and introduce a penalty
function $\psi$ at the left boundary, replacing the agent flow utility by
$U(a_{t},c_{t})=u(c_{t})+\textbf{1}_{a_{t}\leq a_{lb}}\psi(a_{t})$. The
penalty function we use here is the quadratic function:
$\psi(a)=-\frac{1}{2}\kappa(a-a_{lb})^{2}$ where $\kappa$ is a positive
constant and $a_{lb}>\underline{a}$.

_Master equation:_ Let $c^{*}((a,l),z,g)$ denote the equilibrium optimal
household control. Then, for the ABH model, the master equation (2.18)
becomes:

$\displaystyle 0=\mathcal{L}V((a,l),z,g)={}$ $\displaystyle-\rho
V((a,l),z,g)+u(c^{*}((a,l),z,g))+\textbf{1}_{a\leq a_{lb}}\psi(a)$ (5.5)
$\displaystyle+(\mathcal{L}_{x}+\mathcal{L}_{z}+\mathcal{L}_{g})V((a,l),z,g)$
(5.6)

where the operators become:

$\displaystyle\mathcal{L}_{x}V((a,l),z,g)={}$
$\displaystyle\partial_{a}V((a,l),z,g)s((a,l),c^{\ast}((a,l),z,g),r(z,g),w(z,g))$
(5.7) $\displaystyle+\lambda(l)(V((a,\check{l}),z,g)-V((a,l),z,g))$ (5.8)
$\displaystyle\mathcal{L}_{z}V((a,l),z,g)={}$
$\displaystyle\partial_{z}V((a,l),z,g)\eta(\overline{z}-z)+\frac{1}{2}\sigma^{2}\partial_{zz}V((a,l),z,g)$
(5.9) $\displaystyle\mathcal{L}_{g}V((a,l),z,g)={}$
$\displaystyle\sum_{j\in\\{1,2\\}}\int_{\mathbb{R}}D_{g_{j}}V((a,l),z,g)(b)\mu_{g}((b,l_{j})z,g)db$
(5.10)

where $D_{g_{j}}$ is the Frechet derivative with respect to the marginal
density $g(\cdot,l_{j})$ and the KFE is:

$\displaystyle\mu_{g}((a,l),z,g)={}$
$\displaystyle-\partial_{a}\left[s\left((a,l),c^{\ast}((a,l),z,g),r(z,g),w(z,g)\right)g(a,l)\right]$
(5.11) $\displaystyle+\lambda(\check{l})g(a,\check{l})-\lambda(l)g(a,l),$
(5.12)

where $\check{l}$ denotes the complement of $l$, where $r(z,g)$ and $w(z,g)$
solve the system of equations (5.4), and the optimal consumption policy
satisfies the first order optimality condition:

$\displaystyle\partial_{a}V((a,l),z,g)=u^{\prime}(c^{\ast}((a,l),z,g)).$
(5.13)

In the next sections, we solve this master equation numerically using
Algorithm 1. Because the optimal control is a function of the
$\partial_{a}V((a,l),z,g)$, it will turn out to be more convenient to solve
the master equation for the partial derivative, which we denote by
$W((a,l),z,g):=\partial_{a}V((a,l),z,g)$. The parameters that we use in
numerical experiments are in Appendix
LABEL:asec:implementation_details:parameters.

###  5.2 Implementation Details

The implementation details for all methods are summarized in Table 2. Here we
discuss key features about: the neural network structure, the sampling, the
loss function, and the training. We provide a more detailed description of the
implementation in Appendix LABEL:asec:implementation_details.

|  Finite Population |  Discrete State |  Projection
---|---|---|---
Neural Network |  |  |
(i) Structure |  Fully connected feed-forward |  Recurrent with embedding |  Recurrent with embedding
(ii) Activation, $(\phi^{(i)})_{i\leq H}$ |  tanh |  tanh |  tanh
(ii) Output, $\phi^{(H+1)}$ |  soft-plus |  elu activation and factor $(a_{0}+a)^{-\tilde{\eta}}$ |  elu activation and factor $(a_{0}+a)^{-\tilde{\eta}}$
(iii) Layers, $H$ |  5 |  recurrent: 3
embedding: 2 |  recurrent: 3
embedding: 2
(iv) Neurons, $\\#|h|$ |  64 |  recurrent: 100
embedding: 128 |  recurrent: 100
embedding: 64
(v) Initialization |  $W(a,\cdot)=e^{-a}$ |  random |  random
(vi) Auxiliary networks |  none |  consumption |  consumption
Sampling |  |  |
(i) $(a,l)$ |  Active sampling $[a_{min},a_{max}]\times\\{y_{1},y_{2}\\}$ |  Uniform sampling $[a_{min},a_{max}]\times\\{y_{1},y_{2}\\}$ |  Uniform sampling $[a_{min},a_{max}]\times\\{y_{1},y_{2}\\}$
(ii) $(\hat{\varphi}_{i})_{i\leq N}$ |  $\bullet$ Moment sampling: sample $r$ then random distribution of agents to generate $r$ |  $\bullet$ Mixed steady-state sampling
$\bullet$ Ergodic sampling |  $\bullet$ Moment sampling: sample $K$ then other coefficients uniformly
$\bullet$ Ergodic sampling
(iii) $z$ |  $U[z_{min},z_{max}]$ |  $U[z_{min},z_{max}]$ |  $U[z_{min},z_{max}]$
Loss Function |  |  |
(i) Master equation |  $\mathcal{E}^{e}$ |  $\mathcal{E}^{e}$ |  $\mathcal{E}^{e}$
(ii) Constraints |  $\partial_{a}W(a,\cdot)<0$ and $\partial_{z}W(a,\cdot)<0$ |  $\partial_{a}W(a,\cdot)<0$ and $\partial_{z}W(a,\cdot)<0$ |  $\partial_{a}W(a,\cdot)<0$ and $\partial_{z}W(a,\cdot)<0$
(iii) Weights |  $\kappa^{e}=100$, $\kappa^{s}=1$ |  $\kappa^{e}=\kappa^{s}=1$ |  $\kappa^{e}=\kappa^{s}=1$
Training |  |  |
(i) Learning rate |  $10^{-4}$ |  Decaying from
$3\times 10^{-4}$ to $10^{-6}$ |  Decaying from
$3\times 10^{-4}$ to $10^{-6}$
(ii) Optimizer |  ADAM |  ADAM |  ADAM
Table 2: Key Implementation Details

First, consider the neural network structure. For the finite agent
approximation, we use a “plain vanilla” fully connected feed-forward neural
network with 5 layers and 64 neurons per layer. We choose a $\tanh$ activation
function in all hidden layers and a softplus activation in the output layer to
ensure that the output is always positive. For the discrete state and
projection methods, we use an architecture that feeds a fully connected feed-
forward network (referred to as “embedding”) into the one proposed by
Sirignano and Spiliopoulos, (2018) (referred to as “recurrent”). The former
network is used to preprocess the distribution state $\hat{\varphi}$ before
the result is passed to the main (recurrent) network. The embedding portion
has an output dimension of 10, 2 layers, and 128 neurons per layer for the
discrete state space method and 64 neurons per layer for the projection
method. The recurrent portion of the network has 3 layers and 100 neurons per
layer. We use a tanh activation function in all hidden layers and apply an elu
activation in the last layer of the recurrent network to ensure positivity of
the output. In addition, we multiply the neural network output by the factor
$(a_{0}+a)^{-\tilde{\eta}}$, where $a_{0},\tilde{\eta}\geq 0$ are non-
trainable shape parameters. The additional factor aids training by reducing
the amount of curvature in the marginal value function that must be captured
by the neural network. Finally, in both the discrete state and the projection
method, we use the variant of our solution algorithm with staggered updating
of the consumption function and parameterize the latter by a separate
auxiliary neural network that has the same structure as the one for the
marginal value function with the exception that we do not multiply the network
output by $(a_{0}+a)^{-\tilde{\eta}}$ in the output layer.

Second, consider the sampling. We find moment sampling is sufficient for the
finite agent method but not for the other techniques. For the discrete state
method, our main approach is ergodic sampling. As discussed in Section 4.2.1,
this sampling scheme is only meaningful once a sufficiently good guess of the
value function exists. We therefore start with mixed steady state sampling in
early training and gradually increase the fraction of the training set taken
from the ergodic sample. For the projection method, we use a variant of moment
sampling to sample capital stock $K$ from a uniform distribution. To be able
to do this, we rotate the basis, such that the first basis vector points in
the direction of increasing first moments (in the $a$-dimension) and all other
basis vectors are orthogonal to the first. This rotation leaves the space of
approximate distributions unaffected but leads to a one-to-one relationship
between the aggregate capital stock and the first component of the
distribution state $\hat{\varphi}$. We explain this in detail in Appendix
LABEL:asec:master_equations. For the remaining components of $\hat{\varphi}$
we use a combination of uniform and ergodic sampling and gradually increase
the fraction of the training set taken from the ergodic sample.

Third, consider the loss function. In all cases, we impose concavity
constraints on $V$, which are equivalent to monotonicity constraints on $W$.
In addition, we find it is very helpful to impose a sign constraint on
$\partial_{z}W$. This is because the KS model has limited curvature with
respect to $z$ and so the neural network tends to find approximate solutions
where the $z$ component is ignored.

Finally, consider the neural network training. In all cases, we choose an ADAM
optimizer to perform the gradient descent steps. We fix the learning rate at
$10^{-4}$ throughout for the finite agent approach, whereas we use a learning
rate schedule with decaying learning rate (from $3\times 10^{-4}$ to
$10^{-6}$) for the remaining approaches.

###  5.3 Results

We solve the Krusell-Smith model using all three methods. We first discuss the
accuracy of the results for the full model. We then discuss the model without
aggregate shocks and possible ways to use neural network training for
calibration.

####  5.3.1 Mean Reverting TFP Process

For each approach, the error in the neural network approximation to the master
equation is shown in Table 3 below. Evidently all approaches have master
equation losses to the approximate order of $10^{-5}$, which we interpret as
convergence. In Appendix LABEL:sec:robustness we show that all three methods
are robust across different training runs.

| Master equation training loss
---|---
Finite Agent NN | $3.037\times 10^{-5}$
Discrete State Space NN | $9.639\times 10^{-5}$
Projection NN | $8.506\times 10^{-6}$
Table 3:  Neural Networks’ final losses for solving the KS master equation.

In order to sense check our solution, we compare output from our neural
networks to output from traditional approaches. Unfortunately, we do not have
a clear benchmark solution because there is no traditional technique that
provides an arbitrarily precise solution to the model with aggregate shocks.
We choose to compare to the recent approach suggested by Fernández-Villaverde
et al., (2023), which uses a neural network to approximate a statistical law
of motion but solves the Master equation using a finite difference scheme. It
is understood that this technique is good approximation for the KS model. We
make our comparison by computing sample paths from all of our solution
approaches. For the discrete state space and projection methods, we can
generate sample paths by iterating the approximate equilibrium KFE. For the
finite agent method, computing the sample path is more complicated because we
need to average over the transition matrices generated by many draws from the
finite population. We describe this in detail in Algorithm
LABEL:alg:transition_KS in the Appendix LABEL:subsec:KS_simulation.

Figure 1 offers a visual inspection of the difference between our neural
network solutions and Fernández-Villaverde et al., (2023) for a random path of
productivity shocks. The upper-left panel shows the draw from the Ornstein-
Uhlenbeck process: $dz_{t}=\eta(\overline{z}-z)dt+\sigma dB^{0}_{t}$. The
upper-right compares the evolution of capital stock. The second row compares
the evolution of prices. The third and fourth rows compare the population
distribution at various times in the simulation. Evidently we get a similar
path for all variables across all the methods.

In Figure 2, we generate multiple random TFP paths, $z_{t}$, and show the
evolution of our neural network solutions and the Fernández-Villaverde et al.,
(2023) solution in a “fan chart” that displays percentiles for the evolution
of the population. In particular, we generate 1000 TFP paths starting from
$z_{0}=0$ and calculate the corresponding aggregate capital evolution paths.
At each time $t$, for each solution, we compute the pth-percentile of the
capital stock across simulation paths. We then plot the time series of the
percentiles for p equal to 10%, 30%, 50%, 70%, and 90%. Evidently, the finite
agent and projection methods are a close match to the Fernández-Villaverde et
al., (2023). The discrete state space method does a good job at central
percentiles but has more difficulty at the extremes. This reflects our general
experience that the discrete state space approximation is the most difficult
to work with for the KS model and relies most heavily on ergodic sampling,
which decrease it accuracy away from the ergodic mean.

Figure 1: Simulations for the Krusell-Smith Model. The top left plot is the
TFP shock path, the top right panel is the aggregate relative capital change.
The second row left plot shows the relative change in the capital return and
the second row right plot shows the relative change in the wage rate. The
plots on rows three and four show the distribution at different times in the
simulation. The labels “NN, FA”, “NN, DS”, and “NN, P” refer to solutions from
the finite agent, discrete state, and projection neural networks respectively.
_“FV”_ refers to the solution from Fernández-Villaverde et al., (2023).
Subscript sss refers to the stochastic steady state at $z=0$. Figure 2:
Forecasted aggregate capital dynamics starting from the stochastic steady
state (sss) for the Krusell-Smith Model. The left plot is the fan chart for
the TFP shock path, generated from the OU process with initial condition
$z_{0}=0$. The right panel is the time series plot for relative change in
aggregate capital at percentiles 10%, 30%, 50%, 70%, 90% (from the lowest to
the highest). The labels “NN, FA”, “NN, DS”, and “NN, P” refer to solutions
from the finite agent, discrete state, and projection neural networks
respectively. _“FV”_ refers to the solution from Fernández-Villaverde et al.,
(2023).

####  5.3.2 Calibration

As has been suggested by a number of papers (e.g. Duarte et al., (2024)), a
potential benefit of deep learning algorithms is that we can include the
parameters, $\zeta$, as additional inputs into the neural network,
$\hat{V}(\hat{X},\zeta)\approx\mathbb{V}(\hat{X},\zeta;\Theta)$, and then
train the neural network using sampling from both $\hat{X}$ and $\zeta$. In
principle, this means that we can train the model once to get a solution
across the parameter space and then use $\mathbb{V}$ to calculate the moments
for different parameters to calibrate the model. In practice, this requires
that the distribution approximation and sampling do not have high dependence
on the solution. This means that, in our example, the finite agent method is
the natural candidate for testing this approach because it only requires
moment sampling in the training. In Appendix LABEL:subsec:calibration, we
provide a simple example that calibrates the borrowing constraint to match a
particular capital-to-labor ratio. For the other methods this would be much
more difficult because their training requires ergodic sampling.

####  5.3.3 Fixed Aggregate Productivity

For fixed aggregate productivity, $z_{t}=z$, our example model is the
continuous time version of the Aiyagari-Bewley-Huggett (ABH) model discussed
in Achdou et al., 2022a . This model has a precise finite difference solution
and so acts as a more detailed “check” for the accuracy of our solution
technique. Table 4 plots the mean squared error (MSE) between our steady state
solution and the finite difference solution. Figure 3 plots the steady state
consumption policy rule, value function derivative, probability density
function (pdf), and cumulative distribution function (cdf) for the solutions
from the finite agent, discrete state space, and finite difference methods. It
does not make sense to compare to the steady state for the projection method
because the basis has been constructed from steady state solution. Evidently,
the neural network solutions align very closely to the finite difference
solution.

| Master equation loss | MSE(NN, FD)
---|---|---
Finite Agent NN | $3.135\times 10^{-5}$ | $4.758\times 10^{-5}$
Discrete State Space NN | $9.303\times 10^{-6}$ | $6.591\times 10^{-5}$
Table 4: Neural Networks’ final losses for solving the ABH master equation.
Master equation loss is the mean squared error of residuals. MSE(NN,FD) is the
mean squared difference between steady state consumption from the neural
network and finite difference solutions on a wealth grid.  Figure 3:
Comparison between the neural network and finite difference solutions for the
Aiyagari model. The top left plot shows the consumption policy. The top right
shows the derivative of the value function, the bottom left shows the pdf, and
the bottom right shows the cdf. The labels _“NN,FA”_ , _“NN,FS”_ , and _“FD”_
refers to solutions from the finite agent neural network, the discrete state
space neural network, and finite difference respectively.

We can also consider the transition path following an unexpected shock to
aggregate productivity (a so-called “MIT” shock). This makes little sense for
the discrete state space approximation and the projection method because both
rely on types of ergodic sampling and so have difficulty with unanticipated
shocks. However, we were able to train our finite agent approximation using
moment sampling and so it makes sense to consider how successfully the neural
network approximation can handle transition paths. We do the comparison in
Appendix LABEL:subsubsec:transition_dynamics_ABH and find it is remarkably
successful even for large shocks.

###  5.4 Comparison of Techniques

Although the different distribution approximation approaches can all be
effective, we find they have different strengths and weaknesses for solving
the KS model. The finite-agent method is very robust in a number of ways: the
neural network can be trained with the moment sampling procedure, the
algorithm only requires a moderate number of agents (approximately 40), and we
can successfully add parameters as auxiliary states so the model can be solved
across the state and parameter space at the same time.

By contrast, we find the discrete-state method to be difficult to work with.
We believe many of the issues come from having to approximate the derivatives
in the KFE on the discrete state space. One challenge is that this requires a
fine grid for agent wealth (approximately 200 grid points in our example). A
related challenge is that the training samples need to come from relatively
smooth densities and so ergodic sampling is very important. Ultimately, this
makes training the KS model using the discrete state method slow and
complicated. In Section 6.2, we consider a spatial model with locations that
are ex-ante different and agents that choose where to locate instead of having
a consumption saving decision. This means the model is effectively impossible
to solve using the finite agent method but ends up being straightforward to
solve when there are a discrete number of locations. Follow up work in Payne
et al., (2024) extends our approach to search and matching models and also
finds that the discrete state approximation is very effective for that class
of models. This offers suggestive evidence that the discrete state method is
most helpful when the KFE does not contain complicated derivatives.

Finally, the projection method brings a different set of trade-offs. Both the
finite agent and discrete state methods feed relatively large state spaces
into the neural network and then let the neural network work out how to
structure the approximate value function. By contrast, the projection method
requires more choices ex-ante. This allows us to work with a much lower
dimensional approximation (approximately 5 basis functions) and to choose
which statistics of the distribution we want to match most closely in our
approximation. For these reasons, we believe the projection method offers new
advantages for the macroeconomics deep learning literature.

##  6 Additional Examples

We close the paper by considering two additional examples: a model of firms
with capital adjustment costs (similar to a continuous time version Khan and
Thomas, (2007)) and a dynamic spatial model (similar to the model solved using
perturbation in Bilal, (2023)). We have chosen these examples to help
understand the power of the different techniques. The first model is well
suited to a finite agent approximation while the second model can only be
solved by using a discrete set of locations.

###  6.1 Heterogeneous Firms With Adjustment Costs

_Setting._ There is a perishable consumption good and durable capital stock.
The economy consists of a representative household and a continuum of firms
$i\in[0,1]$ that invest in capital and hire labor. There are competitive
markets for goods, labor, firm equity and risk free bonds but capital is not
tradable. We use goods as the numeraire and denote the labor wage by $w_{t}$,
the bond interest rate by $r_{t}$, and the price of equity in firm $i$ by
$p_{t}^{i}$. The representative household chooses consumption, $C_{t}$, labour
supply, $L_{t}$, and the wealth invested in firm $i$ equity, $E_{t}^{i}$, to
solve:

$\displaystyle\max_{C_{t},L_{t},\\{e^{i}\\}}\mathbb{E}\int_{0}^{\infty}e^{-\rho
t}U(C_{t},L_{t})dt,\quad\text{where}\quad
U(C_{t},L_{t})=\frac{1}{1-\gamma}\left(C_{t}-\chi\frac{L_{t}^{1+\varphi}}{1+\varphi}\right)^{1-\gamma}$
(6.1)

and where $\rho$ is the continuous time discount factor, $\gamma$ is the
coefficient of relative risk aversion, $\chi$ determines the disutility of
labor, $\varphi$ is the Frisch elasticity of labor supply, and the household
is subject to the budget constraint:

$\displaystyle{dA_{t}}=-C_{t}dt+w_{t}L_{t}dt+\left[\int_{i}E_{t}^{i}\left(\frac{\pi_{t}^{i}dt+dp_{t}^{i}}{p_{t}^{i}}\right)di\right],\quad\text{where
}\int_{i}E_{t}^{i}di=A_{t}.$ (6.2)

where $\pi_{t}^{i}$ is firm profit. So, following standard analysis, their
stochastic discount factor (or “state-price”) and labor supply are given by:

$\displaystyle\Lambda_{t}={}$
$\displaystyle\partial_{C}U(C_{t},L_{t})=\left(C_{t}-\chi\frac{L_{t}^{1+\varphi}}{1+\varphi}\right)^{-\gamma},$
$\displaystyle L_{t}={}$
$\displaystyle\left(\frac{w}{\chi}\right)^{1/\varphi}$ (6.3)

_Heterogeneous Firms:_ Each firm $i\in[0,1]$ has a production function
$e^{z_{t}}\epsilon_{t}^{i}(k_{t}^{i})^{\theta}(l_{t}^{i})^{\nu}$, where
$z_{t}$ is aggregate TFP following a mean reverting process
$dz_{t}=\eta(\overline{z}-z_{t})dt+\sigma dB_{t}^{0}$ and $\epsilon_{t}^{i}$
is the idiosyncratic shock taking values from
$\\{\epsilon_{L},\epsilon_{H}\\}$, with switching rate
$\lambda_{L},\lambda_{H}$ respectively. The firm can pay dividends or invest
to accumulate capital but faces the adjustment cost
$\psi(n,k)=\frac{\chi_{1}}{2}\frac{n^{2}}{k}$, where $n$ is investment. So,
each firm has idiosyncratic state $x_{t}^{i}=[k_{t}^{i},\epsilon_{t}^{i}]$,
which evolves according to:

$\displaystyle dx_{t}^{i}={}$ $\displaystyle d\begin{bmatrix}k_{t}^{i}\\\
\epsilon_{t}^{i}\end{bmatrix}=\begin{bmatrix}n_{t}^{i}-\delta k_{t}^{i}\\\
0\end{bmatrix}dt+\begin{bmatrix}0\\\
\check{\epsilon}_{t}^{i}-\epsilon_{t}^{i}\end{bmatrix}dJ_{t}^{i}.$ (6.4)

where $J_{t}^{i}$ is the idiosyncratic Poisson process. Taking the wage rate
and the households SDF, $\Lambda_{t}$, as given, firm $i$ chooses labor,
$l_{t}^{i}$, and investment, $n_{t}^{i}$, to maximize their price of equity:

$\displaystyle\max_{l^{i},n^{i}}\left\\{p_{0}^{i}=\mathbb{E}_{0}\int_{0}^{\infty}\frac{\Lambda_{t}}{\Lambda_{0}}e^{-\rho
t}\pi_{t}^{i}dt\right\\}\quad s.t.\quad\eqref{eq:firm:lomx}$ (6.5)

where
$\pi_{t}^{i}:=e^{z_{t}}\epsilon_{t}^{i}(k_{t}^{i})^{\theta}(l_{t}^{i})^{\nu}-w_{t}l_{t}^{i}-\left(n_{t}^{i}+\psi(n_{t}^{i},k_{t}^{i})\right)$
is firm profit. So, for this model $g_{t}$ denotes the population density
across $[k_{t}^{i},\epsilon_{t}^{i}]$ at time t, given a filtration
$\mathcal{F}_{t}^{0}$ generated by the sequence of aggregate productivity
shocks.

_Market Clearing._ Market clearing for goods, labor, and firm equity implies:

$\displaystyle C_{t}={}$
$\displaystyle\sum_{\epsilon=\epsilon_{L},\epsilon_{H}}\int\left[e^{z_{t}}\epsilon_{t}(k_{t})^{\theta}(l_{t})^{\nu}-\left(n_{t}+\psi(n_{t},k_{t})\right)\right]g_{t}(k,\epsilon)dk$
(6.6) $\displaystyle L_{t}={}$
$\displaystyle\sum_{\epsilon=\epsilon_{L},\epsilon_{H}}\int
l_{t}(k,\epsilon)g_{t}(k,\epsilon)dk,\quad E_{t}^{i}=p_{t}^{i}$ (6.7)

So, referring back to our general notation, the price vector that can be
expressed explicitly in terms of $(z,g)$ is $q_{t}=[w_{t},\Lambda_{t}]$, which
satisfies:

$\displaystyle w_{t}={}$
$\displaystyle\left(\sum_{\epsilon=\epsilon_{L},\epsilon_{H}}\int\left(\frac{1}{\nu
z_{t}\epsilon
k^{\theta}}\right)^{\frac{\varphi}{\nu-1}}g_{t}(k,\epsilon)dk\right)^{\frac{1-\nu}{1+\varphi-\nu}},$
(6.8) $\displaystyle\Lambda_{t}={}$
$\displaystyle\left(\sum_{\epsilon=\epsilon_{L},\epsilon_{H}}\int\left[e^{z}\epsilon(k)^{\theta}(l_{t}(k,\epsilon))^{\nu}-n_{t}(k,\epsilon)-\psi(n_{t}(k,\epsilon),k)\right]g_{t}(k,\epsilon)dk-\frac{\chi\left(\frac{w_{t}}{\chi}\right)^{\frac{1+\varphi}{\varphi}}}{1+\varphi}\right)^{-\gamma}$
(6.9)

The price of equity in firm $i$ is the value function of firm $i$,
$p_{t}^{i}=V_{t}^{i}(k_{t},\epsilon_{t})$, which we need to solve for
numerically.

_Master equation._ The aggregate states are $(z_{t},g_{t})$. We let
$V_{t}(k,\epsilon)=V(k,\epsilon,z_{t},g_{t})$ denote the value function in
recursive form.

######  Proposition 2 (Master Equation).

Let $l^{*}((k,\epsilon),z,g)$ and $n^{*}((k,\epsilon),z,g)$ denote the optimal
firm policy functions. Then the master equation (2.18) is given by:

$\displaystyle 0={}$ $\displaystyle-\rho
V((k,\epsilon),z,g)+\Lambda(z,g)\pi_{i}(l^{*}((k,\epsilon),z,g),n^{*}((k,\epsilon),z,g),k,\epsilon,w(z,g))$
(6.10)
$\displaystyle+(\mathcal{L}_{x}+\mathcal{L}_{z}+\mathcal{L}_{g})V((k,\epsilon),z,g)$
(6.11)

where the operators are defined as:

$\displaystyle\mathcal{L}_{x}V((k,\epsilon),z,g)={}$
$\displaystyle\partial_{k}V((k,\epsilon),z,g)(n^{*}((k,\epsilon),z,g)-\delta
k)$ (6.12)
$\displaystyle+\lambda(\epsilon)(V((k,\check{\epsilon}),z,g)-V((k,\epsilon),z,g))$
(6.13) $\displaystyle\mathcal{L}_{z}V((k,\epsilon),z,g)={}$
$\displaystyle\partial_{z}V((k,\epsilon),z,g)\eta(\bar{z}-z)+\frac{1}{2}\partial_{zz}V((k,\epsilon),z,g)\sigma_{Z}^{2}$
(6.14) $\displaystyle\mathcal{L}_{g}V((k,\epsilon),z,g)={}$
$\displaystyle\sum_{j\in\\{1,2\\}}\int_{\mathbb{R}}D_{g_{j}}V((k,\epsilon),z,g)(b)\mu_{g}((b,\epsilon_{j})z,g)db$
(6.15)

where $D_{g_{j}}$ is the Frechet derivative with respect to the marginal
$g(\cdot,l_{j})$ and the KFE is:

$\displaystyle\mu_{g}((k,\epsilon),z,g)={}$
$\displaystyle-\partial_{k}\left[\left(n^{*}((k,\epsilon),z,g)-\delta
k\right)g(k,\epsilon)\right]+\lambda(\check{\epsilon})g(k,\check{\epsilon})-\lambda(\epsilon)g(k,\epsilon),$
(6.16)

The optimal firm policies satisfy:

$\displaystyle n^{*}((k,\epsilon),z,g)={}$
$\displaystyle\frac{k}{\chi_{1}}(\partial_{k}V((k,\epsilon),z,g)-1),$
$\displaystyle l^{*}((k,\epsilon),z,g)={}$ $\displaystyle\left(\frac{w}{\nu
z\epsilon k^{\theta}}\right)^{\frac{1}{\nu-1}}$ (6.17)

###### Proof.

See Appendix LABEL:asubsec:firms. ∎

For computation, it is convenient to define the scaled value function
$\breve{V}:=V/\Lambda$. Then, the new Master equation that we take to the
computer is:

$\displaystyle\breve{\mathcal{L}}^{h}\breve{V}(k,\epsilon,z,g)=\mathcal{L}^{h}\breve{V}(k,\epsilon,z,g)+\frac{\mu_{\Lambda}(z,g)}{\Lambda(z,g)}\breve{V}(k,\epsilon,z,g)+\partial_{z}\breve{V}(k,\epsilon,z,g)\partial_{z}\Lambda(z,g)\sigma_{z}^{2}$
(6.18)

The main technical difference compared to the KS master equation is that the
drift of the household SDF appears in the effective discount rate. This
introduces additional feedback which makes the master equation harder train.
As a result, we start the training without $\mu^{\Lambda}$ and then introduce
the term once the value function has started to converge. We discuss the
details in Appendix LABEL:asubsec:firms. We train the model and get master
equation loss of $7.408\times 10^{-6}$. In Figure 4, we show two plots from
the solution: the investment policy rule and the ergodic distribution for
firms in the low and high state.

Figure 4: Illustration of model solution for the firm dynamics model. The left
panel depicts the investment policy in the stochastic steady state and the
depreciation line (dashed) as a function of firm’s own capital level. The
right panel depicts the (marginal) ergodic density for high and low
productivity firms.

###  6.2 Dynamic Spatial Model

_Setting._ 222222 We describe here a streamlined version of the model
specified in Bilal, (2023) that does not include housing. These changes are
without loss of generality in the sense that, up to a change of parameters,
our resulting master equation is isomorphic to Bilal, (2023). The economy
consists of a finite set of locations $j\in\\{1,...,J\\}$, a continuum of
workers $i\in[0,1]$, and a representative competitive firm at each location
$j$ with an owner who consumes their profits. At any given date $t$, each
worker $i$ resides at a specific location $j$. Workers receive infrequent
moving opportunities. There is a perishable consumption good, tradable across
locations. The representative firm at location $j$ produces consumption goods
according to the production function
$Y_{j,t}=\exp(\beta_{j}+\chi_{j}z_{t})L_{j,t}^{1-\alpha}$, where $L_{j,t}$ is
the labor hired from workers residing at location $j$. As in the KS model,
$z_{t}$ is the aggregate productivity, which follows the same process as in
Section 5.1. $\beta_{j}>0$ is a time-invariant location-specific productivity
shifter and $\chi_{j}>0$ is a parameter that governs the sensitivity of
production in location $j$ to variation in aggregate productivity $z_{t}$.

_Heterogeneous workers:_ Each worker $i\in[0,1]$ has discount rate $\rho$ and
gets flow utility $u(c_{t}^{i})=(c_{t}^{i})^{1-\gamma}/(1-\gamma)$ from
consuming $c_{t}^{i}$ consumption goods at time $t$. The worker is endowed
with one unit of labor, which they supply inelastically in the labor market at
location $j_{t}^{i}$ where they currently reside. The worker has no access to
financial markets and therefore consumes wage income each period,
$c_{t}^{i}=w_{j_{t}^{i},t}$, where $w_{j,t}$ denotes the wage at location $j$
and time $t$. The worker’s idiosycratic state $x_{t}^{i}=j_{t}^{i}$ consists
solely of the location of residence. With arrival rate $\mu>0$, the worker
receives an idiosyncratic opportunity to move location. Upon receiving such an
opportunity, the worker draws idiosyncratic i.i.d. additive preference shocks
for each potential destination $j^{\prime}$ according to a Gumbel distribution
with mean $0$ and inverse scale parameter $\nu$ and then chooses one location
$j^{\prime}$ as their new residence. When moving from $j$ to $j^{\prime}$, the
worker also incurs a moving disutility of $\tau_{j,j^{\prime}}\geq 0$. The
structure of this location choice problem implies that, after optimization,
$j_{t}^{i}$ follows a, possibly time-inhomogeneous, continuous-time Markov
chain with transition rate $\mu\pi_{j,j^{\prime},t}$ from $j$ to $j^{\prime}$
where

$\pi_{j,j^{\prime},t}=\pi_{j,j^{\prime}}(V_{t}):=\frac{e^{\nu(V_{t}(j^{\prime})-\tau_{j,j^{\prime}})}}{\sum_{k=1}^{J}e^{\nu(V_{t}(k)-\tau_{j,k})}}.$

Here, $V_{t}(j)$ denotes the value function of a worker with idiosyncratic
state $j$. The distribution $g_{t}$ is the population density across
$\\{j_{t}^{i}\\}$ at time $t$, given a filtration $\mathcal{F}_{t}^{0}$
generated by the sequence of aggregate productivity shocks.

_Assets, markets, and financial frictions:_ Each period, there are competitive
markets for goods and, in each location, for labor. We use goods as the
numeraire. The price vector $q_{t}=[w_{j,t}:j=1,...,J]$ is given by the
collection of all local wages. Given $g_{t}$ and $z_{t}$, firm optimization
and market clearing imply that the wages $w_{j,t}$ solve:

$\displaystyle w_{j,t}{}$
$\displaystyle=w_{j}(z_{t},g_{t}):=(1-\alpha)\exp(\beta_{j}+\chi_{j}z_{t})(g_{t}(j))^{-\alpha}.$
(6.19)

_Master equation._ The aggregate states are $(z_{t},g_{t})$. We let
$V_{t}(j)=V(j,z_{t},g_{t})$ denote the value function in recursive form.

######  Proposition 3 (Master Equation).

The master equation (2.18) is given by:

$\displaystyle 0=-\rho
V(j,z,g)+u(w_{j}(z,g))+(\mathcal{L}_{x}+\mathcal{L}_{z}+\mathcal{L}_{g})V(j,z,g)$
(6.20)

where the operators are defined as:

$\displaystyle\mathcal{L}_{x}V(j,z,g)$
$\displaystyle=\mu\left(\frac{1}{\nu}\log\left(\sum_{j^{\prime}=1}^{J}e^{\nu(V(j^{\prime},z,g)-\tau_{j,j^{\prime}})}\right)-V(j,z,g)\right)$
$\displaystyle\mathcal{L}_{z}V(j,z,g)$
$\displaystyle=\partial_{z}V(j,z,g)\eta(\overline{z}-z)+\frac{1}{2}\sigma^{2}\partial_{zz}V(j,z,g)$
$\displaystyle\mathcal{L}_{g}V(j,z,g)$
$\displaystyle=\sum_{j^{\prime}=1}^{J}\partial_{g(j^{\prime})}V(j,z,g)\mu_{g}(j^{\prime},z,g),$

the KFE is:

$\displaystyle\mu_{g}(j,z,g)=\mu\left(\sum_{k=1}^{J}\pi_{k,j}(V(\cdot,z,g))g(k)-g(j)\right),$
(6.21)

and the conditional moving probabilities $\pi_{j^{\prime}\prime,j^{\prime}}$
are as defined previously. Up to a transformation of model parameters, this
master equation is equivalent to the one in Bilal, (2023).

###### Proof.

See Appendix LABEL:asubsec:spatial:model_solution. ∎

_Model solution._ Because the model’s idiosyncratic state space is discrete
and all individual density values $g_{t}(j)$ matter directly for prices, the
discrete state space method is the most appropriate for this model. We solve
the model with $J=50$ locations. We provide further details on the
parameterization and numerical implementation in Appendices
LABEL:asubsec:spatial:parameters and LABEL:asubsec:spatial:implementation,
respectively. Here, we only note that, with regard to moving costs
$\tau_{jj^{\prime}}$, we create a simple example based on a “cluster
structure” with a central cluster and several periphery clusters that makes is
more costly to move from/to some locations (the periphery) than others. The
terminal loss after training is $8.403\times 10^{-05}$, similar in magnitude
to the losses for our KS example. Further details on the loss decay are
provided in Appendix LABEL:asubsec:spatial:loss_decay.

Figure 5: Illustration of model solution for the dynamic spatial model. The
left panels depict the wage distribution in the stochastic steady state,
$g=g^{sss}$, $z=0$, as a function of location $j$ (top) and location-specific
productivity $\beta_{j}$ (bottom). The right panels depict the impact effects
of a hypothetical “recession shock” that moves the aggregate state to
$(z,g)=(-0.02,g^{sss})$. The top right panel shows the relative change in
wages in a recession and the bottom right panel the resulting net migration
flows as a proportion of the population ($\mu_{g,j,t}/g_{t}(j)$). In all
panels but the bottom left, locations $j$ are sorted by their sensitivity to
aggregate productivity ($\chi_{j}$) in ascending order. Blue bars/dots depict
periphery locations and red bars/dots depict central locations.

Figure 5 depicts some aspects of the computed model solution. The figure
illustrates wages by location in the stochastic steady state (top left panel)
and how wages depend on the local productivity (bottom left panel) as well as
how a recession (reduction in $z_{t}$) affects wages and population flows.
Except for the scatter plot, locations in the figure are sorted by their
sensitivity to the aggregate productivity shock $\chi_{j}$. Note that there
are several dimensions of ex-ante heterogeneity across locations, so that the
non-monotonic variation in the plots does not indicate training errors. The
overall picture that emerges appears to be economically meaningful. Wages in
harder-hit locations fall more in a recession (top right panel). As a result,
migration flows in the recession are mostly directed from high-$\chi_{j}$ to
low-$\chi_{j}$ regions (bottom right). The exception are primarily very
productive regions in which wages are high to begin with. The cluster
structure of moving costs matters mainly for the baseline level of wages which
tends to be lower in high-productivity regions (bottom left): these regions
attract more workers due to a higher option value of moving.

##  7 Conclusion

This paper proposes a new algorithm that uses deep learning to globally
characterize numerical solutions to continuous time heterogeneous agent
economies with aggregate shocks. We demonstrate our algorithm by solving
canonical models in the macroeconomics literature. Although deep learning
algorithms are straightforward to describe, we find that implementing them
successfully requires careful consideration of the training details. We close
by collecting some practical lessons from our experiences:

  1. 1.

Working out the correct sampling approach is very important.

  2. 2.

Neural networks find it easier to work with smooth constraints.

  3. 3.

Enforcing shape constraints helps with speed and stability.

  4. 4.

Starting with a simple, solvable model is helpful for tuning hyperparameters.

Ultimately, we believe this is only the beginning what these techniques can
achieve. We have demonstrated that we can now globally solve continuous time
heterogeneous agent economies with aggregate shocks, which opens up new and
exciting areas of research.

## References

  * (1) Achdou, Y., Han, J., Lasry, J.-M., Lions, P.-L., and Moll, B. (2022a).  Income and wealth distribution in macroeconomics: A continuous-time approach.  The Review of Economic Studies, 89(1):45–86.
  * (2) Achdou, Y., Lasry, J.-M., and Lions, P. L. (2022b).  Simulating numerically the Krusell-Smith model with neural networks.  arXiv preprint arXiv:2211.07698.
  * Ahn et al., (2018) Ahn, S., Kaplan, G., Moll, B., Winberry, T., and Wolf, C. (2018).  When inequality matters for macro and macro matters for inequality.  NBER macroeconomics annual, 32(1):1–75.
  * Aiyagari, (1994) Aiyagari, S. R. (1994).  Uninsured idiosyncratic risk and aggregate saving.  The Quarterly Journal of Economics, 109(3):659–684.
  * Al-Aradi et al., (2022) Al-Aradi, A., Correia, A., Jardim, G., de Freitas Naiff, D., and Saporito, Y. (2022).  Extensions of the deep Galerkin method.  Applied Mathematics and Computation, 430:127287.
  * Alvarez and Lippi, (2022) Alvarez, F. and Lippi, F. (2022).  The analytic theory of a monetary shock.  Econometrica, 90(4):1655–1680.
  * Alvarez et al., (2023) Alvarez, F., Lippi, F., and Souganidis, P. (2023).  Price setting with strategic complementarities as a mean field game.  Econometrica, 91(6):2005–2039.
  * Auclert et al., (2021) Auclert, A., Bardóczy, B., Rognlie, M., and Straub, L. (2021).  Using the sequence-space Jacobian to solve and estimate heterogeneous-agent models.  Econometrica, 89(5):2375–2408.
  * Azinovic et al., (2022) Azinovic, M., Gaegauf, L., and Scheidegger, S. (2022).  Deep equilibrium nets.  International Economic Review, 63(4):1471–1525.
  * Azinovic and Žemlička, (2023) Azinovic, M. and Žemlička, J. (2023).  Economics-inspired neural networks with stabilizing homotopies.  arXiv preprint arXiv:2303.14802.
  * Barnett et al., (2023) Barnett, M., Brock, W., Hansen, L. P., Hu, R., and Huang, J. (2023).  A deep learning analysis of climate change, innovation, and uncertainty.  arXiv preprint arXiv:2310.13200.
  * Bayraktar et al., (2018) Bayraktar, E., Budhiraja, A., and Cohen, A. (2018).  A numerical scheme for a mean field game in some queueing systems based on Markov chain approximation method.  SIAM Journal on Control and Optimization, 56(6):4017–4044.
  * Bensoussan et al., (2015) Bensoussan, A., Frehse, J., and Yam, S. C. P. (2015).  The master equation in mean field theory.  Journal de Mathématiques Pures et Appliquées, 103(6):1441–1474.
  * Bertucci and Cecchin, (2022) Bertucci, C. and Cecchin, A. (2022).  Mean field games master equations: from discrete to continuous state space.  arXiv preprint arXiv:2207.03191.
  * Bhandari et al., (2023) Bhandari, A., Bourany, T., Evans, D., and Golosov, M. (2023).  A perturbational approach for approximating heterogeneous-agent models.
  * Bilal, (2023) Bilal, A. (2023).  Solving heterogeneous agent models with the master equation.  Technical report, National Bureau of Economic Research.
  * Bilal and Rossi-Hansberg, (2023) Bilal, A. and Rossi-Hansberg, E. (2023).  Anticipating climate change across the united states.  Technical report, National Bureau of Economic Research.
  * Bretscher et al., (2022) Bretscher, L., Fernández-Villaverde, J., and Scheidegger, S. (2022).  Ricardian business cycles.  Available at SSRN.
  * Brzoza-Brzezina et al., (2015) Brzoza-Brzezina, M., Kolasa, M., and Makarski, K. (2015).  A penalty function approach to occasionally binding credit constraints.  Economic Modelling, 51:315–327.
  * Cardaliaguet et al., (2019) Cardaliaguet, P., Delarue, F., Lasry, J.-M., and Lions, P.-L. (2019).  The master equation and the convergence problem in mean field games:(ams-201).  Princeton University Press.
  * Carmona and Laurière, (2021) Carmona, R. and Laurière, M. (2021).  Convergence analysis of machine learning algorithms for the numerical solution of mean field control and games I: the ergodic case.  SIAM Journal on Numerical Analysis, 59(3):1455–1485.
  * Carmona and Laurière, (2022) Carmona, R. and Laurière, M. (2022).  Convergence analysis of machine learning algorithms for the numerical solution of mean field control and games: II—the finite horizon case.  The Annals of Applied Probability, 32(6):4065–4105.
  * Cohen et al., (2024) Cohen, A., Laurière, M., and Zell, E. (2024).  Deep backward and galerkin methods for the finite state master equation.  arXiv preprint arXiv:2403.04975.
  * Delarue et al., (2020) Delarue, F., Lacker, D., and Ramanan, K. (2020).  From the master equation to mean field game limit theory: Large deviations and concentration of measure.  Annals of Probability, 48(1):211–263.
  * Den Haan, (1997) Den Haan, W. (1997).  Solving Dynamic Models with Aggregrate Shocks and Heterogeneous Agents.  Macroeconomic Dynamics, 1(2):355–386.
  * Duarte et al., (2024) Duarte, V., Duarte, D., and Silva, D. (2024).  Machine learning for continuous-time finance.
  * Fernández-Villaverde et al., (2018) Fernández-Villaverde, J., Hurtado, S., and Nuño, G. (2018).  Financial Frictions and the Wealth Distribution.  Working Paper, pages 1–51.
  * Fernández-Villaverde et al., (2023) Fernández-Villaverde, J., Hurtado, S., and Nuno, G. (2023).  Financial frictions and the wealth distribution.  Econometrica, 91(3):869–901.
  * Fouque and Zhang, (2020) Fouque, J.-P. and Zhang, Z. (2020).  Deep learning methods for mean field control problems with delay.  Frontiers in Applied Mathematics and Statistics, 6:11.
  * (30) Germain, M., Laurière, M., Pham, H., and Warin, X. (2022a).  DeepSets and their derivative networks for solving symmetric PDEs.  Journal of Scientific Computing, 91(2):63.
  * (31) Germain, M., Mikael, J., and Warin, X. (2022b).  Numerical resolution of mckean-vlasov fbsdes using neural networks.  Methodology and Computing in Applied Probability, pages 1–30.
  * Glorot and Bengio, (2010) Glorot, X. and Bengio, Y. (2010).  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings.
  * Goodfellow et al., (2016) Goodfellow, I., Bengio, Y., and Courville, A. (2016).  Deep learning.  MIT press.
  * Gopalakrishna, (2021) Gopalakrishna, G. (2021).  Aliens and continuous time economies.  Swiss Finance Institute Research Paper, (21-34).
  * Gopalakrishna et al., (2024) Gopalakrishna, G., Gu, Z., and Payne, J. (2024).  Asset pricing, participation constraints, and inequality.  Princeton Working Paper.
  * Hadikhanloo and Silva, (2019) Hadikhanloo, S. and Silva, F. J. (2019).  Finite mean field games: fictitious play and convergence to a first order continuous mean field game.  Journal de Mathématiques Pures et Appliquées, 132:369–397.
  * Han et al., (2018) Han, J., Jentzen, A., and E, W. (2018).  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510.
  * Han et al., (2021) Han, J., Yang, Y., and E, W. (2021).  DeepHAM: A global solution method for heterogeneous agent models with aggregate shocks.  arXiv preprint arXiv:2112.14377.
  * Hu and Lauriere, (2022) Hu, R. and Lauriere, M. (2022).  Recent developments in machine learning methods for stochastic control and games.  ssrn.4096569.
  * (40) Huang, J. (2023a).  Breaking the curse of dimensionality in heterogeneous-agent models: A deep learning-based probabilistic approach.  Available at SSRN 4649043.
  * (41) Huang, J. (2023b).  A probabilistic solution to high-dimensional continuous-time macro and finance models.
  * Huang and Yu, (2024) Huang, J. and Yu, J. (2024).  Applications of deep learning-based probabilistic approach to “combinatorial” problems in economics.
  * Kahou et al., (2021) Kahou, M. E., Fernández-Villaverde, J., Perla, J., and Sood, A. (2021).  Exploiting symmetry in high-dimensional dynamic programming.  Technical report, National Bureau of Economic Research.
  * Kaplan et al., (2018) Kaplan, G., Moll, B., and Violante, G. L. (2018).  Monetary policy according to hank.  American Economic Review, 108(3):697–743.
  * Khan and Thomas, (2007) Khan, A. and Thomas, J. K. (2007).  Inventories and the business cycle: An equilibrium analysis of (s, s) policies.  American Economic Review, 97(4):1165–1188.
  * Khan and Thomas, (2008) Khan, A. and Thomas, J. K. (2008).  Idiosyncratic shocks and the role of nonconvexities in plant and aggregate investment dynamics.  Econometrica, 76(2):395–436.
  * Krusell and Smith, (1998) Krusell, P. and Smith, A. A. (1998).  Income and Wealth Heterogeneity in the Macroeconomy.  Journal of Political Economy, 106(5):867–896.
  * Lacker, (2020) Lacker, D. (2020).  On the convergence of closed-loop Nash equilibria to the mean field game limit.  Annals of applied probability: an official journal of the Institute of Mathematical Statistics, 30(4):1693–1761.
  * Li et al., (2022) Li, J., Yue, J., Zhang, W., and Duan, W. (2022).  The deep learning Galerkin method for the general stokes equations.  Journal of Scientific Computing, 93(1):1–20.
  * Lions, (2011) Lions, P.-L. (2007-2011).  Lectures at College de France.
  * Lu et al., (2021) Lu, L., Meng, X., Mao, Z., and Karniadakis, G. E. (2021).  DeepXDE: A deep learning library for solving differential equations.  SIAM review, 63(1):208–228.
  * Maliar et al., (2021) Maliar, L., Maliar, S., and Winant, P. (2021).  Deep learning for solving dynamic economic models.  Journal of Monetary Economics, 122:76–101.
  * Min and Hu, (2021) Min, M. and Hu, R. (2021).  Signatured deep fictitious play for mean field games with common noise.  In International Conference on Machine Learning, pages 7736–7747. PMLR.
  * Payne et al., (2024) Payne, J., Rebei, A., and Yang, Y. (2024).  Deep learning for search and matching models.  Princeton Working Paper.
  * Perrin et al., (2022) Perrin, S., Laurière, M., Pérolat, J., Élie, R., Geist, M., and Pietquin, O. (2022).  Generalization in mean field games by learning master policies.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 9413–9421.
  * Prohl, (2017) Prohl, E. (2017).  Discetizing the Infinite-Dimensional Space of Distributions to Approximate Markov Equilibria with Ex-Post Heterogeneity and Aggregate Risk.
  * Raissi et al., (2017) Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2017).  Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations.  arXiv preprint arXiv:1711.10561.
  * Reiter, (2002) Reiter, M. (2002).  Recursive computation of heterogeneous agent models.  manuscript, Universitat Pompeu Fabra, (July):25–27.
  * Reiter, (2008) Reiter, M. (2008).  Solving heterogeneous-agent models by projection and perturbation.  Journal of Economic Dynamics and Control, 33:649–665 Contents.
  * Reiter, (2009) Reiter, M. (2009).  Solving heterogeneous-agent models by projection and perturbation.  Journal of Economic Dynamics and Control, 33(3):649–665.
  * Reiter, (2010) Reiter, M. (2010).  Approximate and Almost-Exact Aggregation in Dynamic Stochastic Heterogeneous-Agent Models.  Technical report, Vienna Institute for Advanced Studies.
  * Sauzet, (2021) Sauzet, M. (2021).  Projection methods via neural networks for continuous-time models.  Available at SSRN 3981838.
  * Schaab, (2020) Schaab, A. (2020).  Micro and macro uncertainty.  Available at SSRN 4099000.