# Deep Learning for Principal-Agent Mean Field Games

Steven Campbell  Department of Statistical Sciences, University of Toronto.
<EMAIL_ADDRESS>,  Yichao Chen  Department of Mathematics,
University of Toronto<EMAIL_ADDRESS>,  Arvind Shrivats
Department of Operations Research & Financial Engineering, Princeton
University<EMAIL_ADDRESS>and  Sebastian Jaimungal  Department of
Statistical Sciences, University of Toronto<EMAIL_ADDRESS>http://sebastian.statistics.utoronto.ca

###### Abstract.

Here, we develop a deep learning algorithm for solving Principal-Agent (PA)
mean field games with market-clearing conditions – a class of problems that
have thus far not been studied and one that poses difficulties for standard
numerical methods. We use an actor-critic approach to optimization, where the
agents form a Nash equilibria according to the principal’s penalty function,
and the principal evaluates the resulting equilibria. The inner problem’s Nash
equilibria is obtained using a variant of the deep backward stochastic
differential equation (BSDE) method modified for McKean-Vlasov forward-
backward SDEs that includes dependence on the distribution over both the
forward and backward processes. The outer problem’s loss is further
approximated by a neural net by sampling over the space of penalty functions.
We apply our approach to a stylized PA problem arising in Renewable Energy
Certificate (REC) markets, where agents may rent clean energy production
capacity, trade RECs, and expand their long-term capacity to navigate the
market at maximum profit. Our numerical results illustrate the efficacy of the
algorithm and lead to interesting insights into the nature of optimal PA
interactions in the mean-field limit of these markets.

SC would like to acknowledge support from the Natural Sciences and Engineering
Research Council of Canada (Alexander Graham Bell Canada Graduate Scholarship
Application No. CGSD3-535625-2019).

SJ would like to acknowledge support from the Natural Sciences and Engineering
Research Council of Canada (grants RGPIN-2018-05705, RGPAS-2018-522715).

##  1\. Introduction

The study of equilibria in systems of interacting agents is a field that has
commanded considerable attention in recent years. Understanding how firms
optimize their strategic decision making while interacting with other
optimizing firms has broad applicability to areas in nature, social sciences,
and finance. In this work, we apply advances in deep learning (DL) to propose
an algorithm that provides a solution to a broad class of complex problems
regarding the equilibrium of a large number of small competing interactive
agents, commonly known as mean field games (MFGs). Moreover, the class of
problems we study incorporates additional layers of complexity, as we consider
equilibria within a principal-agent (PA) formulation, otherwise known as a
principal-agent mean field games (PA-MFGs). Furthermore, the problems that we
study may include market clearing conditions which introduces an additional
interaction into the PA-MFG problems that did not exist before.

Specifically, we consider a PA-MFG problem where the principal induces an
infinite collection of agents to work on their behalf through an incentivized
payment structure. Given the principal’s incentivized payment structure, the
agents form a Nash equilibrium to optimize their goals (the inner problem).
The principal then aims to find a payment structure that best incentivizes the
agents for the principal’s needs (the outer problem). This is a highly non-
trivial problem, both in finding the Nash equilibrium and in finding the
principal’s optimal incentive structure.

Our methodology utilizes deep neural networks (DNN) to solve a McKean-Vlasov
Forward Backward Stochastic Differential Equation (MV-FBSDE) that is
associated with the solution of the inner problem’s MFG. Simultaneously, we
dynamically estimate the principal’s loss surface and implement a gradient
descent algorithm to solve the outer problem. The DNN parametrization of the
MV-FBSDE’s solution allows us to solve complex and high dimensional problems
that classical numerical approaches, such as Least Squares Monte Carlo (LSMC)
or iterative semi-analytic schemes, fail to provide efficient solutions to.
Moreover, the dynamic estimation and concurrent updates for the inner MFG
problem and outer PA problem expedites convergence to a solution.

An application of this particular problem is readily found in the world of
environmental regulation, and particularly emissions markets. Examples of
these include carbon cap-and-trade or Renewable Energy Certificate (REC)
markets. Indeed, regulators may be viewed as principals who choose a market
design to induce some combination of revenue generation and environmental
compliance, with the regulated power generators being agents who interact with
one another through the market and attempt to navigate the regulatory
environment to maximize profit. Thus, the latter half of this paper discusses
the PA-MFG problem in a REC market and presents the results of our proposed
algorithm for a synthetic principal and set of agents.

[56, 10, 9] study optimal behavior within REC (and C&T) markets at the
individual firm level, and [54] from the perspective of a social planner.
These market designs have been discussed in some of the aforementioned works
([10, 9] in particular) as well as in [19, 45]. To the best of the authors’
knowledge, however, this is the first work to explicitly consider a PA-MFG
approach to efficient emission market design.

PA games in general are well-studied. Canonically, they consider a single
principal and agent, with [41] being an early rigorous contribution to this
field. Recent years have seen renewed efforts, thanks in large part to [53]
who investigates a class of infinite-horizon PA games where the the principal
rewards the agent continuously. [53] leverages the dynamic nature of the
agent’s value function arising from the dynamic programming principle. This
allows the principal’s problem to be formulated as a tractable optimal control
problem, and this approach was made rigorous and expanded upon in [21] and
[20], using second order BSDEs.

These works focus on the case where a single principal contracts a single
agent. Many real-world situations (such as the aforementioned REC markets),
however, involve a principal influencing a number of agents. With large
populations of agents, finding equilibria is generally intractable, which
motivates the idea of applying MFG theory to the PA problem. This is exactly
the motivation taken by [23] and [15], the former in a continuous state space,
and the latter in a discrete state space.

Naturally, [23, 15], as well as our own, rely heavily on MFG theory. MFG
theory started with the seminal works of [42, 43] and [47, 48, 49]. The
methodology developed in these concurrent works provides approximate Nash
equilibria in stochastic differential games with symmetric interactions and a
large number of players. Any individual player is impacted by the others
through the empirical distribution of their states, a fact which makes these
problems generally intractable. The key methodological insight of MFG theory
is that the problem is simplified by taking the number of players to infinity,
as the empirical distribution can be replaced by a mean-field distribution.
This allows for a Nash equilibrium to be found in the infinite player limit,
which crucially can be proven to provide an approximate Nash equilibrium for
the finite player game.

Over the past 15 years, MFGs have exploded in popularity, with numerous
applications to engineering ([2, 46, 22]), economics ([35, 33]), and
mathematical finance ([11, 44, 25, 26, 28, 5, 50, 16, 17, 30, 24]). There have
also been efforts to apply MFGs to environmental regulation, as in [55], as
well as commodities ([1, 36, 18, 4, 51]).

The intersection of MFGs and DL has also grown recently. The probabilistic
approach to solving for Nash equilibria in MFGs leads to an MV-FBSDE system.
To solve FBSDE systems, the deep BSDE method [39, 57, 3] has been remarkably
successful. A rigorous theoretical foundation for this method is developed in
[40]. In a related setting, DL techniques have also been used to solve
$N$-player stochastic differential games [37, 38]. More recently, [12, 13]
directly treat MFGs and mean-field type control (MFC) problems in an DL
setting, propose efficient numerical algorithms for arriving at a solution,
and carry out a detailed analysis of their convergence. Similarly, [29] uses
DL techniques to tackle MFC problems with delay. [14] provides an excellent
overview of techniques in this area and their applications to finance.

The contributions of this work are threefold. Firstly, we further generalize
the application of the deep BSDE method seen in [39, 57, 3, 40, 12, 13, 14] to
MV-FBSDEs of the type that have dependence on the law of the adjoint process.
This arises, for example, in settings where there is market clearing in
addition to equilibrium (see, e.g., [55]). Secondly, we develop an alternating
optimization method to solve the principal’s problem concurrently with the
solution to the MV-FBSDEs. In the process, we introduce an equivalent
formulation of the PA problem that lends itself to unconstrained optimization.
To the best of our knowledge, this is the first application of DL techniques
to the PA-MFG setting. Lastly, we implement our algorithm on a non-trivial PA-
MFG problem of independent interest in REC markets. The result yields several
interesting insights into the optimal interactions that arise in the mean-
field regime.

The remainder of this paper is organized as follows. Section 2 provides the
mathematical framework in which we operate, including the formal definition of
the general principal agent problem we wish to solve. Following that, we
present our algorithm in Section 3. We conclude by discussing an application
of our approach to environmental markets in Section 4 and present the
algorithm’s solution to this problem in Section 5.

##  2\. The General Principal Agent mean-field Game Problem

In a typical PA game, there are two classes of decision makers: (i) the agent,
and (ii) the principal. The agent has the ability to influence some output
process, denoted by $X$, through the modification of some control process,
denoted by $\alpha$. The principal delegates the management of this output
process to the agent and is assumed to derive some benefit from it. They
compensate the agent at some terminal time $T$ for the output they produce.
This provides an incentive for the agent to manage $X$ properly. Principals
cannot observe the control of the agent. Rather, they only observe the agent’s
output process. This asymmetry of information leads to moral hazard.

In this work, we assume the principal compensates the agent only through
Markovian functions of the output process. This compensation may be negative,
in which case funds flow from agent to principal. The problem the principal
faces is: What is the optimal compensation structure to balance the output
process that the agent manages and the amount the agent is compensated (or
charged) for their efforts?

The agent wishes to maximize profit, treating the reward structure offered by
the principal as exogenous. Typically, PA problems include a mechanism
referred to as a reservation cost which ensures the agent receives a minimum
reward for their effort. The reservation cost reflects agents preferences to
work on behalf of a principal only if they can expect a minimum payment.

A natural extension to the PA problem described above is to consider what
happens when multiple agents are working on behalf of the principal – such as
those that arise in regulation and oversight. In this case, the principal
(regulator) sets rules to oversee a set of agents (regulated firms) to induce
desired agents’ behavior. Environmental regulation is a clear example of such
a situation, and in Section 4 we provide a specific example of it. Here, the
reservation cost may represent the maximum burden the principal can place on
the set of regulated agents for the regulations to be politically feasible.

In the multiple agent setting, the agents may also interact with each other,
adding a complication to the optimization that an individual agent must solve
to behave optimally. In essence, each agent must solve a multi-player
stochastic game, which is generally intractable. Using an MFG approximation of
the many-player game helps to partially resolve the intractability.

MFGs derive from the idea that stochastic games with large but finite amounts
of agents may be approximated with an infinite number of agents. The MFG
solution then provides strategies that result in near Nash equilibria for the
finite player game. In this work, we do not concern ourselves with the
translation of the infinite player solution back to a solution of the finite
player problem. Rather, we consider only the infinite player problem directly.

Specifically, we consider a PA game with an infinite number of agents. Agents
belong to one of a finite set of sub-populations or agent type. Within sub-
populations agents are assumed to have identical preferences, while across
sub-populations, they are assumed to (possibly) vary through the parameters to
their cost functions.

The advantage of considering the infinite player setting is that the MFG is
significantly simpler to solve than the finite player game. This is due to the
assumption that an individual agent is insignificant as the number of agents
tends to infinity. Thus, we may model the agents’ interactions through a mean-
field distribution representing the population distribution of output
processes. In contrast, finite player games require the agents to interact
directly, which is significantly more complicated.

Therefore, in this PA-MFG setting, our goals are as follows. For a given
reward structure of the principal, we seek the collection of controls that
results in a Nash equilibrium across agents. Similarly, we seek the reward
structure which benefits the principal the most, taking into account the
equilibrium that ensues for any particular admissible reward structure. We
begin by formulating the agent’s problem.

###  2.1. The agents’ problem in a generic PA-MFG

We consider the agents’ problem with $N\rightarrow\infty$ agents in total,
each belonging to one of $K$ different sub-populations. The relevant time
horizon is defined to be ${\mathfrak{T}}:=[0,T]$. We work on the filtered
probability space
$(\Omega,\mathcal{F},(\mathcal{F})_{t\in{\mathfrak{T}}},\mathbb{P})$. All
processes are assumed to be $\mathcal{F}$-adapted unless otherwise stated. We
define the filtration specifically later in this section.

We denote the set of agent types by ${\mathcal{K}}:=\\{1,\cdots,K\\}$, the
subset of agents belonging to sub-population $k\in{\mathcal{K}}$ by
${\mathfrak{N}}_{k}$, and the set of all agents by ${\mathfrak{N}}$. From the
perspective of the agents, the principal exogenously specifies the payment.

Agent $i$, belonging to sub-population $k$, seeks to modulate their controls,
denoted by $(\alpha_{t}^{i})_{t\in{\mathfrak{T}}}$. $\alpha_{t}^{i}$ is
assumed to take values in a set $A\subset\mathbb{R}^{j}$. Through modulating
these controls, the agent modifies the output process which satisfies the
stochastic differential equation (SDE)

(2.1)
$dX_{t}^{i}=b^{k}(t,X_{t}^{i},\bm{\mu_{t}},\alpha_{t}^{i})dt+\sigma^{k}(t)dW_{t}^{i};\qquad
X_{0}^{i}=\xi^{i}\sim F(\Theta^{k})$

Here, $W=\\{W^{i}=(W_{t}^{i})_{t\in{\mathfrak{T}}},i\in{\mathfrak{N}}\\}$ is a
set of $N$ $m$-dimensional independent Wiener processes, and $W^{i}$ is
progressively measurable with respect to the filtration
$\mathcal{F}^{W}:=(\mathcal{F}_{t}^{W})_{t\in{\mathfrak{T}}}\subset\mathcal{F}$.
The function
$b:{\mathfrak{T}}\times\mathbb{R}^{d}\times\mathcal{P}(\mathbb{R}^{d})\times
A:\rightarrow\mathbb{R}^{d}$ is a deterministic and measurable function which
specifies the drift of the output process, and the function
$\sigma:{\mathfrak{T}}\rightarrow\mathbb{R}^{d\times m}$ is a deterministic
and measurable function which specifies the diffusion coefficient of the
output process. $\Theta^{k}$ is a vector valued set of parameters which may
vary across sub-populations.

We make one brief notational remark. Where appropriate, we use a superscript
$\cdot^{(k)}$ rather than a superscript $\cdot^{i}$ to denote quantities
relating to a representative (unspecified index) agent belonging to sub-
population $k$. For example $X_{t}^{(k)}$ represents the output process at
time $t$ for an representative (generic) agent from sub-population $k$. We
also define $\bm{X_{t}}:=(X_{t}^{(k)})_{k\in{\mathcal{K}}}$ as the collection
of $K$ output processes representing the joint output processes from
representative agents from each of the sub-populations.

Next, we make the following definition of the mean-field distribution of
states $\bm{\mu}$.

######  Definition 2.1 (Mean-Field Distribution of States).

In the infinite-player setting, we denote the mean-field distribution of the
output process for agents in sub-population $k\in{\mathcal{K}}$ by
$\mu^{(k)}$. Specifically, we introduce the flow of measures

(2.2)
$\mu^{(k)}=(\mu_{t}^{(k)})_{t\in{\mathfrak{T}}},\qquad\mu_{t}^{(k)}\in{\mathcal{P}}(\mathbb{R}^{d}),\;\forall
t\in{\mathfrak{T}}$

for all $k\in{\mathcal{K}}$, where $\mathcal{P}(\mathbb{R}^{d})$ represents
the space of probability measures on $\mathbb{R}^{d}$, such that
$\mu_{t}^{(k)}(A)$ is the probability that a representative agent from sub-
population $k$ has an output process belonging to the set
$A\in\mathcal{B}(\mathbb{R}^{d})$, at time $t$. Furthermore, we define

(2.3)
$\bm{\mu}=(\\{\mu_{t}^{(k)}\\}_{k\in{\mathcal{K}}})_{t\in{\mathfrak{T}}}$

to be the flow of the collection of all mean-field measures.

We also make an assumption about the initial condition of the output process.

######  Assumption 2.2.

The initial states $\\{\xi^{i}\\}_{i\in{\mathfrak{N}}_{k}}$ are identically
distributed, mutually independent, and independent of $\mathcal{F}^{W}$.
Moreover, there exists a $c>0$ such that $\mathbb{E}[\|\xi^{i}\|^{2}]\leq
c<\infty$, for all $i\in{\mathfrak{N}}_{k}$ and $k\in{\mathcal{K}}$.

Next, we denote by

(2.4) $\mathcal{G}_{t}^{i}:=\sigma\left((X_{u}^{i})_{u\in[0,t]}\right),$

the filtration that the $i$-th agent must adapt their strategy to, which is
the sigma field generated by the $i$-th firm’s output process. Note that we
assume that all firms have knowledge of the initial distribution (but not the
actual value of) other firms’ output processes. The full filtration
$\mathcal{F}=(\mathcal{F}_{t})_{t\geq 0}$ is defined as
$\mathcal{F}_{t}=\bigvee_{i\in{\mathfrak{N}}}\mathcal{G}_{t}^{i}$. The set of
square integrable controls is defined as

(2.5)
$\mathbb{H}_{t}^{2}:=\left\\{\alpha:\Omega\times{\mathfrak{T}}\rightarrow
A\subset\mathbb{R}^{j}\;\bigg{|}\;\mathbb{E}\left[{\int_{0}^{T}}(\alpha_{t})^{2}\,dt\right]<\infty\right\\}.$

The admissible set of controls for the agent are defined in the assumption
below.

######  Assumption 2.3.

The set of admissible controls for firm $i\in{\mathfrak{N}}$ is

(2.6)
$\displaystyle{\mathcal{A}}^{i}:=\left\\{\alpha\in\mathbb{H}_{t}^{2}\,\,\text{s.t.
}\alpha\text{ is $\mathcal{G}^{i}$-adapted}\,\,\right\\}.$

Each agent chooses their controls with the goal of minimizing the cost they
incur through their acceptance of the principal’s contract. Specifically, for
a fixed $\bm{\mu}$, the $i$-th firm (belonging to sub-population $k$) aims to
minimize the cost functional $J^{A,i}:{\mathcal{A}}^{i}\rightarrow\mathbb{R}$,
where

(2.7)
$J^{A,i}(\alpha;\bm{\mu})=\mathbb{E}\left[\int_{0}^{T}f^{k}(t,X_{t}^{i},\bm{\mu_{t}},\alpha_{t}^{i})dt+g(X_{T}^{i})\;\bigg{\rvert}\;\mathcal{G}_{0}^{i}\right].$

The functions $f^{k}$ and $g$ are commonly referred to as the running costs
and terminal costs, respectively. The terminal cost can be viewed as the
control of the principal, but in the context of the agents’ problem, it is
exogenous and fixed. This terminal cost may also be a function of the terminal
mean-field distribution, but we focus on the case when it does not for
simplicity.

Agents within a sub-population have the same cost parameters, and
consequently, those agents act, in equilibrium, in a similar manner. Each
individual agent’s strategy, however, is adapted to their own inventory and as
such, agents’ strategies are not identical, even within the same sub-
population.

As previously stated, all agents seek to minimize their own costs, and we seek
the optimal control for all agents simultaneously. Specifically, we seek a
Nash equilibrium. This is a collection of controls
$\\{(\alpha^{i,\star})\in{\mathcal{A}}^{i}:i\in{\mathfrak{N}}\\}$ and mean-
field flow $\bm{\mu^{\star}}$ (induced by the collection of controls) such
that for any other admissible control $\alpha^{i}$ we have:

(2.8) $\displaystyle J^{A,i}(\alpha^{i,\star},\bm{\mu^{\star}})$
$\displaystyle\leq J^{A,i}(\alpha^{i},\bm{\mu^{\star}})\qquad\forall
i\in{\mathfrak{N}}_{k},\,\forall k\in{\mathcal{K}},$ (2.9)
$\displaystyle\mathbb{P}_{X_{t}^{(k)}}$
$\displaystyle=\mu_{t}^{(k)}\qquad\forall k\in{\mathcal{K}}.$

Equation (2.8) describes the optimality of the collection of controls, and
indicates that no agent can improve their position by unilaterally deviating
from it. Equation (2.9) ensures that the collection of controls implies a
controlled output process for the representative agent (recall agents are
assumed to have identical preferences within sub-populations) from each sub-
population that is consistent with the mean-field distribution used in the
value function.

For an exogenous $g$, let
$V^{A,k}(g)=\inf_{\alpha\in{\mathcal{A}}^{\cdot}}J^{A,\cdot}(\alpha;\bm{\mu^{\star}})$
represent the optimally controlled value function for a representative agent
in sub-population $k$, $\forall\;k\in{\mathcal{K}}$.

To solve (2.8) and (2.9), we ultimately express the problem in terms of an MV-
FBSDE. There are many possible approaches to obtaining such an MV-FBSDE. Two
common ones are applying the so-called probabilistic method of MFGs (see [8]
for a detailed pedagogical treatment) or tools from variational analysis (see
e.g., [55, 16, 27]).

In this work, we assume that the agents’ MFG problem has already been cast
into its MV-FBSDE form.

######  Assumption 2.4.

We assume the Nash equilibria characterized by (2.8) and (2.9) may be
expressed through the solution of an MV-FBSDE. That is, for each
$i\in{\mathfrak{N}}_{k}$,

(2.10) $\displaystyle dX_{t}^{i}$
$\displaystyle=\varphi^{k}(t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},Y_{t}^{i},\mathbb{P}_{\bm{Y_{t}}})dt+\sigma^{k}(t)dW_{t}^{i};\qquad
X_{0}^{i}=\xi^{i}\sim F(\Theta^{k})$ (2.11) $\displaystyle dY_{t}^{i}$
$\displaystyle=\rho^{k}(t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},Y_{t}^{i},\mathbb{P}_{\bm{Y_{t}}})dt+Z_{t}^{i}dW_{t}^{i};\qquad
Y_{T}^{i}=\partial_{x}g(X_{T}^{i}).$

Specifically, we assume solving (2.10)-(2.11) results in a triple
$(X_{t}^{i,\star},Y_{t}^{i,\star},Z_{t}^{i,\star})$ such that the optimal
$\alpha^{i,\star}$ can be recovered from the solution and $\alpha^{i,\star}$
solves (2.8) with $\bm{\mu^{\star}}:=\mathbb{P}_{\bm{X_{t}^{\star}}}$ (which
naturally solves (2.9)). Here,
$\bm{X_{t}^{\star}}:=(X_{t}^{\star,(k)})_{k\in{\mathcal{K}}}$, is the
collection of $K$ output processes from representative agents across the sub-
populations which solve the MV-FBSDE above.

######  Remark 2.5.

When there is no market clearing,
$\varphi^{k}(t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},Y_{t}^{i},\mathbb{P}_{\bm{Y_{t}}})=b^{k}((t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},\hat{\alpha}^{k}(t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},Y_{t}^{i}))$
and
$\rho^{k}(t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},Y_{t}^{i},\mathbb{P}_{\bm{Y_{t}}})=-\partial_{x}\mathcal{H}^{k}(t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},Y_{t}^{i},\hat{\alpha}^{k}(t,X_{t}^{i},\mathbb{P}_{\bm{X_{t}}},Y_{t}^{i}))$
where $\mathcal{H}$ represents the Hamiltonian of the problem and there is no
need for dependence on the law $\mathbb{P}_{\bm{Y_{t}}}$. Specifically,
$\mathcal{H}^{k}(t,x,\mu,y,\alpha)=\langle
b^{k}(t,x,\mu,\alpha),y\rangle+f^{k}(t,x,\mu,\alpha)$ and the functions
$\hat{\alpha}^{k}$ are the optimizers of the $\mathcal{H}^{k}$.

This is a rich class of problems for which existence and uniqueness results
are well-developed. See [6] for a set of assumptions that results in such an
MV-FBSDE, for a solution exists and is unique. Furthermore, this class
includes MFGs that involve endogenous price formation, such as the MFG problem
discussed in [55], [31], and [34]. In Section 4, we discuss an MFG model with
endogenous price formation pertaining to REC markets. It is these cases of
endogenous price formation that can give rise to a representation where there
is dependence of the drift coefficients on the law of the adjoint process.

###  2.2. Principal’s problem in a generic PA-MFG

The principal aims to determine the terminal payment structure $g(X_{T})$ that
optimizes their own criterion. The principal derives a benefit from the output
process and compensates the agent for said output at the terminal time $T$.

We restrict ourselves to payment structures such that the set
$\mathcal{M}(g):=\\{(\alpha,\bm{\mu}):(\alpha,\bm{\mu})\text{ solves
}\eqref{eq:MFG_optimal}-\eqref{eq:MFG_consistency}\\}$ is non-empty. That is,
the principal chooses over a set of functions such that the agents’ MFG is
well-posed. Furthermore, the principal ensures the agents’ reservation cost is
met (in expectation), across all sub-populations.

We denote the set $\Xi$ as the set of admissible terminal payment structures
for the principal to choose among. Let $\mathbf{V}^{\bm{\mu}}(g):=(\int
V^{A,k}(g)d\mu_{0}^{(k)})_{k\in{\mathcal{K}}}$, $\Pi\geq 0$ be a full rank
$s\times|{\mathcal{K}}|$ matrix (where the inequality is understood
componentwise) and $\mathbf{R}_{0}$ be a $s$-dimensional vector of reservation
costs. In particular, we say

(2.12)
$\Xi:=\left\\{g:\mathbb{R}^{d}\rightarrow\mathbb{R}\,\bigg{\lvert}\,\mathbb{E}[g(X^{(k)}_{T})^{2}]<\infty,\,g\text{
convex},\,\mathcal{M}(g)\neq\emptyset,\,\Pi\mathbf{V}^{\bm{\mu}^{\star}}(g)\leq\mathbf{R}_{0}\,\forall
k\in{\mathcal{K}}\right\\}.$

If $\Pi=\mathrm{Id}$, then the constraint corresponds to a classic upper bound
on the reservation costs for all agents. If
$\Pi=\bm{\pi}\in\Delta_{|{\mathcal{K}}|}$ (where $\Delta_{\cdot}$ is the unit
simplex), then the constraint corresponds to a restriction on the (weighted)
average cost across sub-populations.

If $\mathcal{M}(g)$ is not a singleton, we assume the agent chooses the one
that is best for the principal. As a result, the principal’s problem may be
written as the solution to

(2.13)
$V^{P}:=\inf_{g\in\Xi}\mathbb{E}\left[\sum_{k\in{\mathcal{K}}}\pi_{k}\left(-U\left(X_{T}^{(k)}\right)-\lambda
g(X_{T}^{(k)})\right)\right],$

where $U$ is a concave and non-decreasing utility function, and $\lambda>0$ is
an exogenous principal parameter indicating their willingness to trade off
between the output process and the contract they pay. Here, $\pi_{k}$
represents the fraction of the total population attributable to sub-population
$k$ in the mean-field limit. The objective is to find $V^{P}$, the
corresponding function $g$, and use them to characterize the resulting Nash
equilibrium among the agents.

In general, the principal’s problem is challenging to solve analytically, even
with the simplifications of uncontrolled volatility and restricting to Markov
contract structures. The general approach to solving variations of these
problems analytically is outlined by [23], which was inspired by the work of
[21] and [53]. Broadly speaking, the idea is to equate the solution of the MFG
with the solution of a BSDE, and use this BSDE to define a set of coupled
McKean-Vlasov SDEs, which describe (among other things) the dynamics of the
agents’ optimal value process. From here, the principal’s problem can be
turned into a MFC problem, with the agents’ MFG solution embedded into the
problem definition. This can either be solved using PDEs (as in [23]) or using
the probabilistic approach (as in [7]). In both cases, however, analytic
solutions are impossible in all but the simplest of situations. Classical
numerical approaches exist, but as always, one faces the Curse of
Dimensionality. Instead, we take a DL approach and present an optimization
algorithm that can more flexibly solve these type of problems.

##  3\. Numerical Algorithm

###  3.1. Approximating the Principal’s Objective

We begin with some simplifying assumptions on $g$. In particular, we restrict
to functions that are continuous, monotone, convex, and take a real valued
argument. This is a classic setting for many MFGs. Convexity is often required
to derive the form of the MV-FBSDE in Assumption 2.4 and monotonicity is
natural as it relates to reward/cost structures that are monotone in the
terminal state. It is known that such functions $g$ can be approximated
uniformly on compact sets by piecewise linear functions (see for instance
[52]). Moreover, a convex piecewise linear function can be constructed through
a positive linear combination of call and put payoffs. We use this insight to
formulate the following.

######  Assumption 3.1.

The function $g:\mathbb{R}\to\mathbb{R}$ is continuous and can be written as

(3.1)
$g(x)=\phi_{0}+\sum_{j=1}^{\bar{N}}w_{j}(x-R_{j})^{+},\quad\text{or}\quad
g(x)=\phi_{0}+\sum_{j=1}^{\bar{N}}w_{j}(R_{j}-x)^{+},$

for weights $\phi_{0}\in\mathbb{R},w_{j}\in\mathbb{R}_{+}$ and knot points
$R_{j}\in\mathbb{R}$, $j=1,...,{\bar{N}}$.

######  Remark 3.2.

For uniqueness and existence results, we often also require $g$ to be
differentiable with a Lipschitz continuous derivative. However, since this
holds in the formulation above, except at a finite number of points, and we
use a simulation approach the lack of differentiability at these points does
not affect the implementation. Nevertheless, while we do not treat it
explicitly here, it is straightforward to regularize the piecewise linear
approximations to produce smooth versions of them.

Without loss of generality we treat the first case in (3.1). For a fixed set
of knot points the principal’s objective function now becomes

$\mathbb{E}\left[\sum_{k\in{\mathcal{K}}}\pi_{k}\left(-U(X_{T}^{\star,(k)})-\lambda\phi_{0}-\lambda\sum_{j=1}^{\bar{N}}w_{j}(R_{j}-X_{T}^{\star,(k)})^{+}\right)\right],$

where the optimally controlled state processes are associated with the mean-
field equilibrium induced by the principal’s input weights $\mathbf{w}$. The
following simple lemma will prove useful as it allows us to eliminate one of
the problem constraints.

######  Lemma 3.3.

Under Assumption 3.1 the principal’s problem can be equivalently expressed as

(3.2)
$\inf_{\mathbf{w}\in\hat{\Xi}}\left(\lambda\max_{i\in[s]}\frac{\Pi^{\top}_{i}\mathbf{V}^{\bm{\mu}^{\star}}(\hat{g})-R_{0,i}}{\Pi_{i}^{\top}\mathbf{1}}+\mathbb{E}\left[\sum_{k\in{\mathcal{K}}}\pi_{k}\left(-U(X_{T}^{\star,(k)})-\lambda\hat{g}(X_{T}^{\star,(k)})\right)\right]\right)$

where $X_{T}^{\star,(k)}$ is the optimally controlled process of a
representative agent from population $k\in{\mathcal{K}}$ in equilibrium,

(3.3) $\hat{g}(x):=\sum_{j=1}^{\bar{N}}w_{j}(R_{j}-x)^{+},$

and

(3.4)
$\hat{\Xi}:=\left\\{\mathbf{w}\in\mathbb{R}^{\bar{N}}\bigg{\lvert}\mathbb{E}[\hat{g}(X_{T}^{\star,(k)})^{2}]<\infty,\,\mathcal{M}(\hat{g})\neq\emptyset,\
k\in{\mathcal{K}}\right\\}.$

Here $\Pi_{i}$ is the vector corresponding to the $i$th row of the full rank
matrix $\Pi\in\mathbb{R}_{+}^{s\times|\mathcal{K}|}$, $[s]:=\\{1,...,s\\}$,
$R_{0,i}$ is the $i$th component of $\mathbf{R}_{0}$, and the agents’ optimal
controls for $\hat{g}$ and $g$ coincide.

###### Proof.

Consider in the definition of $\Xi$ the constraint

$\mathbf{V}^{\bm{\mu}^{\star}}(g)\leq\mathbf{R}_{0}.$

Note that the agents’ optimal control is invariant to the choice of
$\phi_{0}$. Using this observation we may pull $\phi_{0}$ out from the agents’
objective. Re-arranging yields that

$\phi_{0}\Pi\mathbf{1}\leq\mathbf{R}_{0}-\Pi\mathbf{V}^{\bm{\mu}^{\star}}(\hat{g}).$

From this and the positivity of $\Pi$ we can conclude

$\phi_{0}\Pi\mathbf{1}\leq\mathbf{R}_{0}-\Pi\mathbf{V}^{\bm{\mu}^{\star}}(\hat{g})\
\ \ \iff\ \ \
\phi_{0}\leq\min_{i\in[s]}\frac{R_{0,i}-\Pi^{\top}_{k}\mathbf{V}^{\bm{\mu}^{\star}}(\hat{g})}{\Pi_{k}^{\top}\mathbf{1}}.$

We can divide by $\Pi_{k}^{\top}\mathbf{1}$ since $\Pi$ is full rank and thus
has no rows that are identically $0$. Furthermore, if $\phi_{0}$ is optimal
then equality holds in the above. If not, the principal can always improve
their objective by increasing $\phi_{0}$ to this upper bound since it has no
impact on the agents’ controls - a contradiction. On the other hand, if we
restrict to the above equality, minimizing the principal’s objective over
$\mathbf{w}\in\hat{\Xi}$ determines the minimizing value of $\phi_{0}$. The
claim in the lemma then follows by substitution. ∎

We also enforce the following standard assumptions that simplifies the
optimization.

######  Assumption 3.4.

The following conditions hold:

  1. (1)

For any $\alpha\in\mathcal{A}^{k}$, $\phi_{0}\in\mathbb{R}$, and
$\mathbf{w}=(w_{1},\dots,w_{\bar{N}})^{\top}\in\mathbb{R}^{\bar{N}}_{+}$, the
drift coefficient $b^{k}(t,x,p,\alpha)$ and diffusion coefficient
$\sigma^{k}(t)$ are such that
$\mathbb{E}\left[g(X^{(k)}_{T})^{2}\right]<\infty$, for all
$k\in{\mathcal{K}}$.

  2. (2)

There exists a $E\subseteq\mathbb{R}^{\bar{N}}$ such that $\mathbf{w}\in E$ if
and only if $\mathcal{M}(g)\not=\emptyset$.

  3. (3)

There exists an almost-everywhere differentiable mapping $\psi$ such that
$\psi[\mathbb{R}^{\bar{N}}]=E$.

This assumption allows us to replace the constraint $\mathbf{w}\in\hat{\Xi}$
with $\mathbf{w}\in E$ which is straightforward to implement numerically. In
particular, if Assumption 3.4(3) holds we can optimize directly on
$\mathbb{R}^{\bar{N}}$ by defining $\mathbf{w}:=\psi(\mathbf{u})$ for
$\mathbf{u}\in\mathbb{R}^{\bar{N}}$.

###  3.2. Approximating the MV-FBSDE

To solve the PA problem we first need a mechanism for solving the agent
problem through the MV-FBSDE (2.10)-(2.11) for different weights $\mathbf{w}$
specified by the principal. We tackle this by discretizing said MV-FBSDE, and
parameterizing the co-adjoint process and initial condition by neural nets
akin to the deep BSDE/FBSDE techniques employed in [40, 39, 32, 13]. We note
that a new addition in our setting is the dependence on the law of the adjoint
process. We briefly introduce a new notation, whereby $[\cdot]^{l,(k)}$
represents the $l$-th sampled path of $[\cdot]$ for a representative agent
from sub-population $k$ in our discretization. We use $l$ to avoid any
possible confusion with agent indexing as used in prior sections, and do not
mix agent indexing and representative indexing to maintain clarity.

Letting $\mathbb{T}=\\{t_{0},...,t_{M}\\}$ be a discrete set of points with
$t_{0}=0$ and $t_{M}=T$ we define an ensemble of neural nets
$Y^{\theta_{0}^{(k)}}$, $(Z^{\theta_{m}^{(k)}})_{m=1}^{M}$ and introduce the
time discretized version of (2.10)-(2.11):

(3.5) $\displaystyle X_{t_{m}}^{l,(k)}$
$\displaystyle=X_{t_{m-1}}^{l,(k)}+\varphi^{k}(t,X_{t_{m-1}}^{l,(k)},\hat{\mathbb{P}}_{\mathbf{X_{t_{m-1}}}},Y_{t_{m-1}}^{l,(k)},\hat{\mathbb{P}}_{\mathbf{Y_{t_{m-1}}}})\Delta
t+\sigma^{k}(t_{m-1})\Delta W^{i}_{t_{m}}$ (3.6) $\displaystyle X_{0}^{l,(k)}$
$\displaystyle=\xi^{l}$ (3.7) $\displaystyle Y_{t_{m}}^{l,(k)}$
$\displaystyle=Y_{t_{m-1}}^{l,(k)}+\rho^{k}(t,X_{t_{m-1}}^{l,(k)},\hat{\mathbb{P}}_{\mathbf{X_{t_{m-1}}}},Y_{t_{m-1}}^{l,(k)},\hat{\mathbb{P}}_{\mathbf{Y_{t_{m-1}}}})\Delta
t+Z^{\theta_{m}^{(k)}}(X_{t_{m-1}}^{l,(k)},Y_{t_{m-1}}^{l,(k)})\Delta
W^{l}_{t_{m}}$ (3.8) $\displaystyle Y_{0}^{l,(k)}$
$\displaystyle=Y^{\theta_{0}^{(k)}}(X_{0}^{l,(k)})$

for $\Delta t=T/M$, $\Delta W_{t_{m}}^{i}=W^{l}_{t_{m}}-W^{l}_{t_{m-1}}$ and
$m=1,...,M$. Here $\bm{\theta}:=(\theta_{m}^{(k)})_{m\in[N],k\in\mathcal{K}}$
are the network parameters, and $\hat{\mathbb{P}}_{\cdot}$ denotes the
empirical distribution over the samples. To solve the original MV-FBSDE we
have also sampled finitely many paths ($\mathcal{N}_{k}$ paths) for each sub-
population $k\in\mathcal{K}$. For compactness of notation we write
$\mathbf{X}^{(k)}_{t}:=(X^{1,(k)}_{t},...,X_{t}^{\mathcal{N}_{k},(k)})^{\top},\mathbf{Y}^{(k)}_{t}:=(Y^{1,(k)}_{t},...,Y_{t}^{\mathcal{N}_{k},(k)})^{\top}$
for $k\in\mathcal{K}$ and $t\in\mathbb{T}$ to represent the samples of each
population. When the superscript $k$ is omitted, we refer to the samples
across all sub-populations.

The intuition of the algorithm is that we optimize over the parameters
$\bm{\theta}$ in order to match the terminal condition for $Y_{T}^{l,(k)}$. To
this end, we define the MV-FBSDE loss function on a set of sample paths as

(3.9)
$\mathcal{L}_{F}(\bm{\theta})=\frac{1}{|\mathcal{K}|}\sum_{k\in\mathcal{K}}\frac{1}{\mathcal{N}_{k}}||\mathbf{Y}^{(k)}_{t_{M}}-\partial_{x}g(\mathbf{X}^{(k)}_{t_{M}})||_{2}^{2}.$

We say that we have an approximate solution to the MV-FBSDE when the loss
function has converged within a specified tolerance to $0$. This loss function
implicitly depends on the weights $\mathbf{w}$ through the definition of $g$.

###  3.3. Learning the Principal Objective Function

To tackle the principal problem we discretize the running costs using a
Trapezoidal rule to optimize the objective with respect to the weights
$\mathbf{w}$. Here, for a fixed set of weights, we simulate $\mathcal{N}_{k}$
sample paths for each sub-population $k\in{\mathcal{K}}$ to approximate the
expectation. We also assume a uniform grid $\mathbb{T}$ so that we may write

$\displaystyle\mathrm{Running\ Costs}_{l}^{(k)}$ $\displaystyle\approx\Delta
t\sum_{m=1}^{M-1}f^{k}(t,X^{l,(k)}_{t_{m}},\hat{\mathbb{P}}_{\mathbf{X}_{t_{m}}},\hat{\alpha}_{t_{m}}^{l})+\frac{\Delta
t}{2}\sum_{m\in\\{0,M\\}}f^{k}(t,X^{l,(k)}_{t_{m}},\hat{\mathbb{P}}_{\mathbf{X}_{t_{m}}},\hat{\alpha}_{t_{m}}^{l}).$

Our principal loss on a sample is given by

(3.10)
$\hat{\mathcal{L}}_{P}(\mathbf{w})=\lambda\max_{i\in[s]}\frac{\Pi^{\top}_{i}\hat{\mathbf{V}}^{\hat{\mathbb{P}}}(\hat{g})-R_{0,i}}{\Pi_{i}^{\top}\mathbf{1}}+\sum_{k\in{\mathcal{K}}}\frac{\pi_{k}}{\mathcal{N}_{k}}\sum_{l=1}^{\mathcal{N}_{k}}\left(-U(X_{T}^{l,(k)})-\lambda\hat{g}(X_{T}^{l,(k)})\right),$

where the $k^{th}$ component of $\hat{\mathbf{V}}^{\hat{\mathbb{P}}}(\hat{g})$
is given by:

$\hat{\mathbf{V}}_{k}^{\hat{\mathbb{P}}}(\hat{g}):=\frac{1}{\mathcal{N}_{k}}\sum_{l=1}^{\mathcal{N}_{k}}\left(\mathrm{Running\
Costs}_{l}^{(k)}+\hat{g}(X^{l,(k)}_{t_{M}})\right).$

We use this to estimate the true principal loss at $\mathbf{w}$.

1 Initialize forward network parameters $\mathbf{\theta}$. Initialize
principal parameters $\Upsilon$. Initialize memory buffer $\mathfrak{M}$.
Initialize $\mathbf{u}^{(0)}$

2 $N_{O}:=$ max # of outer steps, $N_{S}\geq\bar{N}:=$ local sample size,
$N_{F}:=$ # MV-FBSDE steps, $N_{A}:=$ # loss approximation steps, $N_{P}:=$ #
principal loss steps, $(\mathcal{N}_{k})_{k\in{\mathcal{K}}}:=$ MV-FBSDE
sample size, $N_{B}:=$ # loss approximation batch size;

3 Define weight error tolerance $TOL$ and MV-FBSDE error tolerance $TOL_{F}$;

4 Define update rule for sampling radius $\epsilon$, network parameters
$\bm{\theta},\Upsilon$ and values $\mathbf{u}^{(\cdot)}$;

5 foreach  _$n$ $\mathrm{in}$ $[N_{O}]$_ do

6 Initialize local samples
$\mathcal{U}:=(\mathbf{u}_{i}^{(j)})_{i\in[N_{S}]}\subset
B_{\epsilon}(\mathbf{u}^{(j)})$;

7 foreach  _$\mathbf{u}^{\prime}\in\mathcal{U}$_ do

8 foreach  _$k\ \mathrm{in}\ [N_{F}]$_ do

9 Sample paths of the discretized MV-FBSDE (3.5)-(3.8) by generating initial
conditions, Brownian noise, and passing through the ensemble network
$\bm{\theta}$ using the sample sizes $(\mathcal{N}_{k})_{k\in{\mathcal{K}}}$;

10 Compute the MV-FBSDE loss $\mathcal{L}_{F}(\bm{\theta})$ (see (3.9));

11 Update $\bm{\theta}$;

12 if  _$\mathcal{L}_{F}(\bm{\theta}) <TOL_{F}$_ then

13 break;

14

15 end if

16

17 end foreach

18 Estimate the principal loss via
$\hat{\mathcal{L}}_{P}(\psi(\mathbf{u}^{\prime}))$ (see (3.10));

19 Store
$(\mathbf{u}^{\prime},\hat{\mathcal{L}}_{P}(\psi(\mathbf{u}^{\prime})))$ in
$\mathfrak{M}$;

20

21 end foreach

22 foreach  _$k\ \mathrm{in}\ [N_{A}]$_ do

23 Sample from $\mathfrak{M}$ a batch
$\mathbf{U}:=(\mathbf{u}_{m})_{m\in[N_{B}]}$ and associated principal losses
$\hat{\mathcal{L}}_{P}(\psi(\mathbf{U})):=(\hat{\mathcal{L}}_{P}(\psi(\mathbf{u}_{m})))_{m\in[N_{B}]}$;

24 Compute the MSE loss
$||\mathcal{L}^{\Upsilon}_{P}(\psi(\mathbf{U}))-\hat{\mathcal{L}}_{P}(\psi(\mathbf{U}))||_{2}^{2}$;

25 Update $\Upsilon$;

26

27 end foreach

28 foreach  _$k$ $\mathrm{in}$ $[N_{P}]$_ do

29 Compute $\mathcal{L}^{\Upsilon}_{P}(\psi(\mathbf{u}^{(j)}))$;

30 Update $\mathbf{u}^{(j)}\mapsto\mathbf{u}^{(j+1)}$;

31 Increment $j\mapsto j+1$;

32

33 end foreach

34 Update $\epsilon$;

35 if  _$||\mathbf{u}^{(j)}-\mathbf{u}^{(j-N_{P})}||_{2} <TOL$_ then

36 break;

37

38 end if

39

40 end foreach

41

Algorithm 1 Principal-Agent Optimization

Given this estimate of the principal’s loss we aim to optimize it over the
weights, which requires gradients. To compute an approximate gradient, we
estimate the principal loss function over the weights locally. In particular,
we sample $N$ points uniformly in an $\epsilon$-ball about the current
estimate of the weights. The principal loss on each of these samples may be
estimated by solving the MV-FBSDE. With these estimates of the loss within an
$\epsilon$-ball, we use a second neural net,
$\mathcal{L}_{P}^{\Upsilon}(\mathbf{w})$ with parameters $\Upsilon$ that fits
the target $\hat{\mathcal{L}}_{P}(\mathbf{w})$. From this estimate, we may
obtain estimates of the gradients and perform an approximate gradient update.
For simplicity, the optimization over weights proceeds on
$\mathbb{R}^{\bar{N}}$ by using $\psi$ from Assumption 3.4. More specifically,
we take gradient of $\mathcal{L}_{P}^{\Upsilon}(\psi(\mathbf{u}))$ with
respect to the input $\mathbf{u}$ and perform a gradient step to determine the
next value of $\mathbf{u}$. This, in turn, defines updated values of
$\mathbf{w}$.

###  3.4. Algorithm

We now have all the components necessary to describe the algorithm for solving
the PA problem fully. Our approach may be viewed as having two main
components: an inner optimization (for the agents) and an outer optimization
(for the principal). The approach works at a high level by (i) initializing
$\mathbf{u}^{(0)}\in\mathbb{R}^{\bar{N}}$, (ii) performing $N_{F}$
optimization steps for the inner agent problem, (iii) fixing the agents’
controls, estimate the principal’s loss function, and (iv) performing $N_{P}$
gradient updates of the outer problem. This is repeated until the Principal’s
weights and the MV-FBSDE loss are within a given tolerance. The details are
provided in Algorithm 1.

##  4\. Example

As a motivating example of a real-world problem where our methodology could be
applied, we consider Renewable Energy Certificate (REC) markets. The optimal
behavior of agents in REC markets is covered in detail from a MFG perspective
in [55], however, we extend their model in two significant ways: first by
extending the problem to include the regulator as a principal, and second by
expanding the state and control space of the model for a more realistic
formulation.

###  4.1. The Basics of REC Markets

REC markets belong to the class of so-called market-based emissions regulation
policies, the most well known of which are carbon cap-and-trade (C&T) markets.

REC markets are a closely related alternative to C&T markets. In REC markets,
a regulator sets a floor on the amount of energy generated from renewable
sources for each firm (based on a percentage of their total energy
generation), and provides certificates for each MWh of energy produced via
these means. This is also known as a Renewable Portfolio Standard (RPS), and
typically apply to private Load Serving Entities (LSEs), also known as
electricity suppliers. To ensure compliance, each firm must surrender
certificates totaling the floor at the end of a compliance period, with a
monetary penalty paid for each lacking certificate. The certificates are
traded assets, allowing regulated LSEs to make a choice about whether to
produce electricity from renewable means themselves, or purchase the
certificates on the market (or a mix of both). In either case, whether
purchasing or producing clean, RECs induce the production of clean energy, as
all purchased RECs must have been generated through clean generation means.

In practice, these systems regulate multiple consecutive and disjoint
compliance periods, which are linked together through a mechanism called
banking, where unused allowances in period-$n$ can be carried over to
period-$(n+1)$. However, we consider a simpler single-period framework in this
example.

###  4.2. Modeling of REC markets

We discuss the single-period framework for REC markets with the following
rules. The market governs a compliance period from $[0,T]$ denoted as
${\mathfrak{T}}$. A firm obtains RECs in each period, with their terminal RECs
denoted by $X_{T}$. At time $T$ a firm must pay (or possibly is paid)
$C(X_{T})$, which represents a generic penalty chosen by the regulator, which
is made known to the firms. We assume there are no costs after time $T$.

Firms receive RECs by generating electricity through a particular energy
source (depending on the market, this could be solar, nuclear, etc.). One REC
typically corresponds to one MWh of electricity produced via the target energy
source. A firm may also purchase or sell RECs on the market. After $T$, all
firms forfeit any remaining RECs. $T$ can be thought of as ‘the end of the
world’ – there are no costs associated with any time after this.

####  4.2.1. The Agent Problem in REC Markets

Agents belonging to sub-population $k\in{\mathcal{K}}$ are assumed to have a
baseline REC generation rate $h_{t}^{k}$ at which they generate with zero
marginal cost. As discussed in Section 2.1, agents preferences are assumed to
be identical within sub-populations, but are distinct across them. All agents
are minor agents; that is, they have no market impact individually.

Agents have the ability to further increase $h_{t}^{k}$ through their rate of
expansion activities. In essence, this models building more REC generation
capacity. In the context of a Solar REC market, this can be thought of as
installing more solar panels. This control is denoted by
$(\alpha_{t}^{i})_{t\in{\mathfrak{T}}}$, with the units being REC capacity
added per unit time.

The agents also have the ability to rent short-term REC generation capacity
for an assumed quadratic cost, denoted by the control
$(g_{t}^{i})_{t\in{\mathfrak{T}}}$. Once again, this represents a rate, this
time of generation capacity rented per unit time.

Finally, agents may engage in the REC market and purchase/sell RECs as needed,
with their purchasing rate in the market denoted by
$(\Gamma_{t}^{i})_{t\in{\mathfrak{T}}}$. Generation expansion and rental must
be positive, while trading rate may be positive or negative.

We denote the collection of expansion rates, generation rental rates, and
trading rates by $\bm{\alpha_{t}}:=(\alpha_{t}^{i})_{i\in{\mathfrak{N}}}$,
$\bm{g_{t}}:=(g_{t}^{i})_{i\in{\mathfrak{N}}}$ and
$\bm{\Gamma_{t}}:=(\Gamma_{t}^{i})_{i\in{\mathfrak{N}}}$, respectively. The
addition of expansion generation as a control differentiates this example from
[55]. This model could be made even more realistic by incorporating a delay
between the decision to expand and the increase to the baseline generation
rate.

The agents therefore must keep track of two state variables: their REC
inventory, and their REC capacity. The relevant state variables for agent $i$
in sub-population $k$ are denoted $X^{i}$ and $C^{i}$ and they satisfy the
SDEs

(4.1) $\displaystyle dX_{t}^{i}$
$\displaystyle=(h_{t}^{k}+C_{t}^{i}+\Gamma_{t}^{i})dt+\sigma^{k}dW_{t}^{i},\qquad$
$\displaystyle X_{0}^{i}=\xi^{i}\sim F(\Theta^{k})$ (4.2) $\displaystyle
dC_{t}^{i}$ $\displaystyle=\alpha_{t}^{i}dt+0\,dB_{t}^{i},\qquad$
$\displaystyle C_{0}^{i}=0.$

We include the diffusive term in $C_{t}^{i}$ even though its coefficient is 0
to make clear that in principle it may be stochastic as well – e.g., a firm
pays for and plans expansion capacity but the capacity is not delivered.

As in Section 2.1, let $\mu_{t}^{(k)}$ represent the mean-field distribution
of states. That is for all $t\in{\mathfrak{T}}$,
$\mu_{t}^{(k)}\in\mathcal{P}(\mathbb{R}^{2})$. Similarly,
$\bm{\mu}=(\\{\mu_{t}^{(k)}\\}_{k\in{\mathcal{K}}})_{t\in{\mathfrak{T}}}$.
Moreover, $\mu_{t}^{(k),X}$ and $\mu_{t}^{(k),C}$ denotes the marginal mean-
field distributions of the denoted state. 2.2 and 2.3 are assumed to hold.

The agents seek to profit maximize, and in particular, they aim to optimize
the following cost functional:

(4.3) $\displaystyle J^{A,i}(\alpha,g,\Gamma;\bm{\mu})=\mathbb{E}\biggl{[}$
$\displaystyle\int_{0}^{T}\frac{\zeta^{k}}{2}(g_{u}^{i})^{2}+\frac{\gamma^{k}}{2}(\Gamma^{i}_{u})^{2}+\frac{\beta^{k}}{2}(\alpha_{u}^{i})^{2}+S^{\bm{\mu}}_{u}\Gamma^{i}_{u}du+\phi_{0}+\sum_{j=1}^{N}w_{j}(R_{j}-X_{T}^{i})^{+}\;\rvert\;\mathcal{G}_{0}^{i}\biggr{]},$

with the form of the terminal condition given in 3.1. $S^{\bm{\mu}}$
represents the equilibrium REC price. This cost functional can be thought of
as the sum of various costs. The controls all have quadratic costs associated
with them, with the latter two terms representing the cost (or profits) from
trading and the non-compliance penalty imposed by the principal, respectively.
$\zeta^{k},\gamma^{k},\beta^{k}$ are scalar cost parameters which are
identical for firms within a sub-population, and possibly different across
them.

The REC price is obtained endogenously through a market clearing condition, as
in [55, 31, 34]. The clearing conditions imply that

(4.4)
$\lim_{N\rightarrow\infty}\tfrac{1}{N}\sum_{i\in{\mathfrak{N}}}\Gamma_{t}^{i}=0,$

and this condition induces an expression for the REC Price $S^{\bm{\mu}}$
provided in (4.12). The clearing condition amounts to the average trading rate
(across agents) vanishing at all times in the limit of large number of
players. Note that in the individual agents’ cost functional (4.3), the mean-
field distribution arises only through the REC price.

As in Section Section 2.1, the agents seek to find a Nash equilibrium as
defined in (2.8) and (2.9). One can show (following the steps of [55], or the
probabilistic approach espoused in [6]) that the solution for agent $i$ in
sub-population $k$ for all $i\in{\mathfrak{N}}_{k}$, $k\in{\mathcal{K}}$, can
be found through the solution to the following MV-FBSDE

(4.5) $\displaystyle dX_{t}^{i}$
$\displaystyle=(h_{t}^{k}-\left(\tfrac{1}{\zeta^{k}}+\tfrac{1}{\gamma^{k}}\right)Y_{t}^{i,X}-\tfrac{1}{\gamma^{k}}S_{t}^{\bm{\mu}}+C_{t}^{i})dt+\sigma^{k}dW_{t}^{i},\;$
$\displaystyle X_{0}^{i}=\xi^{i},$ (4.6) $\displaystyle dC_{t}^{i}$
$\displaystyle=-\tfrac{1}{\beta^{k}}Y_{t}^{i,C}dt+0dB_{t}^{i},\;$
$\displaystyle C_{0}^{i}=0,$ (4.7) $\displaystyle dY_{t}^{i,X}$
$\displaystyle=Z_{t}^{i,X}dW_{t}^{i},\;$ $\displaystyle
Y_{T}^{i,X}=\partial_{x}g(X_{T}^{i}),$ (4.8) $\displaystyle dY_{t}^{i,C}$
$\displaystyle=-Y_{t}^{i,X}dt+Z_{t}^{i,C}dB_{t}^{i},\qquad$ $\displaystyle
Y_{T}^{i,C}=0,$

with the optimal controls given by

(4.9) $\displaystyle\alpha_{t}^{i,\star}$
$\displaystyle=-\tfrac{Y_{t}^{i,C}}{\beta^{k}},$ (4.10) $\displaystyle
g_{t}^{i,\star}$ $\displaystyle=-\tfrac{Y_{t}^{i,X}}{\zeta^{k}},$ (4.11)
$\displaystyle\Gamma_{t}^{i,\star}$
$\displaystyle=-\tfrac{Y_{t}^{i,X}-S_{t}^{\bm{\mu}}}{\gamma^{k}},\quad\text{and}$
(4.12) $\displaystyle S_{t}^{\bm{\mu}}$
$\displaystyle=\frac{-1}{\sum_{k\in{\mathcal{K}}}\frac{\pi_{k}}{\gamma_{k}}}\sum_{k\in\mathcal{K}}\frac{\pi_{k}}{\gamma^{k}}\mathbb{E}\left[Y_{t}^{(k),X}\right].$

Therefore, the agents’ problem fits into the formulation described in Section
2.1. We now discuss the principal’s problem.

####  4.2.2. The Principal’s Problem in REC Markets

The principal aims to solve the optimization problem:

$\displaystyle\inf_{\phi_{0}\in\mathbb{R},\mathbf{w}\in\mathbb{R}_{+}^{\bar{N}}}J^{P}\left(\phi_{0},\mathbf{w}\right)=\inf_{\phi_{0}\in\mathbb{R},\mathbf{w}\in\mathbb{R}_{+}^{\bar{N}}}\mathbb{E}\left[\sum_{k\in\mathcal{K}}\pi_{k}\left(-\lambda\left(\phi_{0}+\sum_{j=1}^{m}w_{j}(R_{j}-X_{T}^{(k),\star})^{+}\right)-X_{T}^{(k),\star}\right)\right],$

where $X_{t}^{(k),\star}$ represents the optimal controlled REC inventory for
a representative agent in sub-population $k$, and
$\mathbf{w}=(w_{1},\dots,w_{\bar{N}})^{\top}$ subject to the constraint

$\bm{\pi}^{\top}\mathbf{V}^{\bm{\mu}^{\star}}(g)\leq R_{0},$

where $\bm{\pi}\in\Delta_{|{\mathcal{K}}|}$ is the weight vector whose entries
correspond to the fraction of the market in the mean-field that is
attributable the sub-populations $k\in{\mathcal{K}}$. When expanded this may
be written

$\sum_{k=1}^{K}\pi_{k}\int
J^{A,(k)}(\alpha^{(k),\star},g^{(k),\star},\Gamma^{(k),\star};\bm{\mu^{\star}})d\mu_{0}^{(k),\star}\leq
R_{0},$

where $(\alpha^{(k),\star},g^{(k),\star},\Gamma^{(k),\star})$ represent the
collection of optimal controls across representative agents from each group
$k\in\mathcal{K}$ and $J^{A,(k)}$ denotes their cost functional.
$\bm{\mu}^{\star}$ represents the induced mean-field flow from the optimal
controls. The constant $R_{0}$ represents the reservation cost.

It is possible (via Lemma 3.3) to see that the principal’s problem is
equivalent to the modified problem:

$\inf_{\mathbf{w}\in\mathbb{R}_{+}^{N}}\mathbb{E}\left[\sum_{k\in\mathcal{K}}\pi_{k}\left(-\lambda\left(R_{0}-\mathrm{Running\
Costs}^{k,\star}\right)-X_{T}^{k,\star}\right)\cdot\right]$

where $\mathrm{Running\
Costs}^{k,\star}=\int_{0}^{T}\frac{\zeta^{k}}{2}(g_{u}^{(k),\star})^{2}+\frac{\gamma^{k}}{2}(\Gamma^{(k),\star}_{u})^{2}+\frac{\beta^{k}}{2}(\alpha_{u}^{(k),\star})^{2}+S^{\bm{\mu}^{\star}}_{u}\Gamma^{(k),\star}_{u}du$.
In this representation we have

$\phi_{0}=R_{0}-\sum_{k\in\mathcal{K}}\pi_{j}\mathbb{E}\left[\mathrm{Running\
Costs}^{k,\star}+\sum_{j=1}^{N}w_{j}(R_{j}-X_{T}^{k,\star})^{+}\right].$

Note that $X_{T}^{k,\star}$ and the running costs depend implicitly on
$\mathbf{w}$ through the optimal controls.

Both the agent and principal’s problems fit into the overall framework
discussed in Section 2.1 and Section 2.2, more specifically, Algorithm 1 may
be applied to them. Moreover, this problem is much harder to solve using the
methods described in [55], with the increased state space posing a challenge
to the semi-analytic iterative scheme proposed therein. Thus, the algorithm
unlocks a problem that is otherwise infeasibly difficult with standard
numerical methods.

##  5\. Results

We next proceed to implementing the algorithm on the PA REC problem and
discuss the results. We first consider a relatively simpler model of the REC
market example discussed above, as a method to establish some evidence of the
appropriateness of the solution output by Algorithm 1. Specifically, we first
consider the scenario where the terminal cost $g(X_{T})$ as in (3.1) has only
a single knot point (i.e. $\bar{N}$ = 1). Upon validating our solution in this
setting, we proceed to implement our algorithm on a similar synthetic example
where $\bar{N}=10$. This allows for a richer class of terminal penalty
functions to be searched over from the principal’s perspective, and results in
a more interesting optimization problem. In this section, we use the
parameters given in Table 1 and Table 2 throughout. The distribution of
initial inventory $X_{0}^{(k)}\sim\xi^{(k)}=\mathcal{N}(v^{k},\eta^{k})$, for
all agents in sub-population $k$, for all $k\in{\mathcal{K}}$.

$\Delta_{t}$ | $T$ | $K$ | $\lambda$ | $R_{0}$
---|---|---|---|---
$1/52$ | $1$ | $2$ | $6$ | $0$
Table 1. Compliance Parameters Sub-population | $\pi_{k}$ | $h^{k}$ | $\sigma^{k}$ | $\zeta^{k}$ | $\gamma^{k}$ | $v^{k}$ | $\eta^{k}$ | $\beta^{k}$
---|---|---|---|---|---|---|---|---
$k=1$ | $0.25$ | $0.2$ | $0.1$ | $1.75$ | $1.25$ | $0.6$ | $0.1$ | 1.0
$k=2$ | $0.75$ | $0.5$ | $0.15$ | $1.25$ | $1.75$ | $0.2$ | $0.1$ | 1.0
Table 2. Model Parameters

These parameters are not calibrated to data, but chosen to provide stylized
differences between the sub-populations that will be reflected in their
optimal controls, as in [55] and [56]. In particular, we note that a firm from
population 1 is more likely to have a much greater initial allocation of RECs
than one from population 2. Furthermore, they have greater rental costs, but a
lower cost of trading RECs, as well as a significantly lower baseline REC
generation capacity. These differences in parameters results in optimal
controls that delineate sub-populations for their particular needs in this
system.

###  5.1. Single knot non-compliance penalty

For the single knot model, we fix $g$ such that it has a single knot point at
$R_{1}=0.9$. Assuming $R_{0}=0$, the terminal condition $\partial_{x}g$ of the
MV-FBSDE becomes a step function

$\displaystyle\begin{cases}\partial_{x}g(x)=-w&x\leq 0.9,\\\
\partial_{x}g(x)=0&x>0.9.\end{cases}$

We apply Algorithm 1 in this case, and in Figure Figure 1 plot (left) the
trajectory of $w$ across the optimization, and (right) the estimated principal
loss across a set of candidate $w$’s. Specifically, we conducted a grid search
for $w\in[0.05,0.40]$ by applying our algorithm and saving the trained
ensemble model for each $w$, which we denote by $\bm{\theta}_{w}$. For each of
these models, we sampled $100$ new batches of Brownian noise and initial
conditions, collecting the resulting means and standard errors of the
principal losses. From these plots, we see that the optimal weight suggested
by the algorithm quickly entered a band roughly between $0.18$ and $0.24$ and
then eventually settles near $0.205$. Meanwhile, the curve of principal loss
appears to be convex in $w$, with the minimum occurring right around this
number, albeit with a relatively flat area surrounding it.

Figure 1. Trajectory of $w$ during optimization (left) and estimated principal
loss as a function of the weight $w$ (right) within the single knot test case
($R_{1}=0.9)$. Parameters as in Table 1 and Table 2.

###  5.2. Multiple knot non-compliance penalty

The single knot example discussed above provides some level of evidence that
our algorithm can properly identify the optimal non-compliance penalty of the
principal. We now progress to a more involved example. In this subsection we
fix $10$ nodes from $0.8$ to $1.16$ with an equal spacing of $0.04$ (see Table
3 for the precise knot points). With the multi-knot setting, we are able to
optimize over a richer class of non-increasing convex terminal costs, and one
that in fact encompasses the simple single knot setting.

In Figure 2, we plot the optimal weight trajectories throughout the
optimization (left), as well as the principal loss trajectory (right). The
principal loss converges quickly to approximately $-0.81$, while it takes
longer for the weights $w_{i}$ themselves to converge. This indicates that the
minimizer of the principal loss is located in a relatively flat region. It is
worth noting that the converged $w$ has many of its coordinates close to zero.
In fact, only the weights for the knots at $0.8,0.84,1.04$ and $1.08$ RECs are
non-negligible (see Table 3.

Figure 2. Trajectory of optimal weights (left) and principal loss (right) within the multiple knot test case. Parameters as in Table 1 and Table 2. Knots as in Table 3. Node | $0.8$ | $0.84$ | $0.88$ | $0.92$ | $0.96$ |  $1$. | $1.04$ | $1.08$ | $1.12$ | $1.16$
---|---|---|---|---|---|---|---|---|---|---
Weight | $0.0193$ | $0.0483$ | $0.0050$ | $0.0025$ | $0.0014$ | $0.0018$ | $0.1097$ | $0.0101$ | $0.0077$ | $0.0063$
Table 3. Knots and corresponding optimal weight in non-compliance penalty for
multiple knot test case. Parameters as in Table 1 and Table 2.

We use the output of the optimization to plot the optimal non-compliance
penalty of the principal, and the distribution of induced non-compliance costs
across the samples of representative agents from the differing sub-populations
in Figure 3. Note that the samples do not represent a finite set of agents
competing in the same market and rather represent samples from the
(approximate) mean-field distribution.

We observe that the non-compliance penalty $g$ is much steeper for $x<1.04$
than that for $x\geq 1.04$ , indicating a higher marginal penalty for
insufficient inventory below $1.04$. Turning our attention to the right panel
in Figure 3, we see that a representative firm from population $2$ (blue) pays
more on average in non-compliance penalties than one from population $1$
(red).

Figure 3. Optimal terminal condition (left) and distribution of non-compliance
costs paid by the agents (right). Population 1 indicated in red, population 2
indicated in blue. Parameters as in Table 1 and Table 2. Knots as in Table 3.
Parameters as in Table 1 and Table 2.

There are many other quantities of interest with respect to the two sub-
populations, such as their controls and state processes. In particular, we are
interested in how they evolve across the period, as well as their
distributions at $T$.

In Figure 4, we plot the sampled trajectories of the inventories for
representative agents in each sub-population, as well as a histogram of their
terminal inventory, with the non-negligible knot points indicated by dashed
lines. In both sub-populations, a considerable number of firms pay the highest
marginal penalty for insufficient inventory, indicated by the notable
probability mass to the left of the vertical line $x=0.8$ in the right panel.
We also include various terminal inventory percentiles for each sub-population
in Table 4.

Figure 4. Inventory paths (left) and terminal inventory (right). Population 1 indicated in red, population 2 indicated in blue. The lines denote the non-negligible nodes in the terminal penalty $g$. Sub-population |  $10$-th percentile |  $25$-th percentile |  $50$-th percentile |  $75$-th percentile |  $90$-th percentile
---|---|---|---|---|---
$k=1$ | $0.79$ | $0.85$ | $0.92$ | $0.99$ | $1.05$
$k=2$ | $0.7$ | $0.78$ | $0.88$ | $0.97$ | $1.05$
Table 4. Percentiles of the terminal inventories across sub-populations.
Parameters as in Table 1 and Table 2.

Interestingly, we observe from Figure 4 that representative agents from
population 1 had a significantly higher average initial inventory than
population 2, but their terminal inventories are far closer. This occurs for a
variety of reasons. For one, representative agents from population 2 are
induced more to acquire RECs due to the nature of the compliance penalty. They
also have a much greater initial baseline of REC generation ($h_{t}^{k})$, and
a lower rental cost ($\zeta^{k})$, allowing them to acquire RECs more quickly.

By plotting the measure flow for each sub-population across time, as in Figure
5, we can see this behavior more clearly. We can also see that the mean of the
inventory distribution for population $1$ shifts to the right as time
progressed, with the standard deviation of the distribution appearing fairly
constant. On the other hand, the mean of the inventory distribution for
population $2$ shifts to the right more aggressively as time goes to $1$, with
the standard deviation of the distribution increasing. This suggests greater
variance in the REC acquisition rate of a representative firm from population
$2$.

Figure 5. Inventories for population $1$ (left) and population $2$ (right)
throughout time. Parameters as in Table 1 and Table 2.

We now next investigate the agents’ optimal controls. Firms may obtain RECs in
one of three ways: rental generation, expansion generation, and trading. We
wish to understand which of these primarily drove REC acquisitions, and to
that end, we plot the acquisition rates and totals for each sub-population
across each of these methods in Figure 6.

Figure 6. Expansion rate and total expansion (left column), rental generation
rate and cumulative rental generation (middle column), trading rate and net
trading position (right column) across sampled paths of representative agent
from each sub-population. Population 1 indicated by red, population 2
indicated by blue. Parameters as in Table 1 and Table 2.

From the left panels of Figure 6, we observe expansion rates for each sub-
population taper off to 0 as time progresses, as expected. Firms accrue less
benefit from expanding as the terminal time approaches, as they have less time
to reap the fruits of their labor. We find the sampled representative firms
from sub-population 2 generally expand more than those from sub-population 1,
despite the fact that they have the same cost parameters associated with
expansion ($\beta^{k}$). Again, this is expected due to their generally lower
initial inventory.

From the middle panels, we see similar behaviour with rental generation. The
rental generation rates for population $2$ were greater than the generation
rates for population $1$ throughout the time horizon, resulting in greater
cumulative rental generation. This is expected due to the lower cost parameter
firms from population 2 have for rental generation, and their lower initial
inventory incentivizing them to acquire RECs quickly.

From the right panels, we observe the net trading volume for population $2$
was bigger than that of population $1$. This means representative firms from
population $1$ were more likely to be net buyers of RECs than representative
firms from population $0$. Taken together, these plots provide reasoning for
and internal consistency with the previously exemplified behavior of a
representative firm from sub-population $2$ generally acquiring more RECs than
one from sub-population $1$.

Figure 7. Total generation(top left), rental generation(top right), total
trading volume (bottom left) and net trading volume (bottom right) for
Population $1$ in red and Population $2$ in blue.

Figure Figure 7 plots the distribution across samples of representative firms
from each sub-population of total generation, rental generation, total trading
volume, and rental trading volume. The increased variance in total generation
for population $2$ is apparent, which is consistent with the variance of the
inventory of a representative firm from said population increasing in time, as
seen in Figure 5.

Finally, with the same argument as in [55], one can show that the equilibrium
price process is constant. This fact is observed in Figure 8, and can be
intuited by noting in (4.12) that the equilibrium price is in terms of the
expectation of the adjoint processes $(Y_{t}^{(k),X})_{t\in\mathfrak{T}}$
which are, in fact, martingales.

Figure 8. Approximate equilibrium price path.

##  6\. Conclusion

In this paper we introduced a multi-step DL algorithm for solving Principal-
Agent mean field Games (PA-MFGs) and illustrated its use on a previously
intractable problem in REC markets that involves market clearing. This
application demonstrates the flexibility of our method and leads to several
insights into the nature of equilibria in REC markets which are of independent
interest. This study also suggests several future directions which we describe
below.

  1. (i)

This work leaves open the question of convergence guarantees for the numerical
algorithm presented in Section 3. As in the analysis conducted in [12, 13], it
is of interest to understand how well, and under what assumptions, our
numerical scheme approximates the solution of the original PA-MFG problem.
[37] lays a theoretical foundation for the deep BSDE method in the general
case of coupled FBSDEs, and show that the error converges to zero given the
universal approximation capability of NNs. This could be a good starting point
for convergence results of our proposed algorithm, as our algorithm for PA-MFG
problem is an alternating optimization whose inner loop involves solving
coupled FBSDEs with a similar deep neural net approach.

  2. (ii)

In Section 4 we introduce a stylistic model for SREC markets. In the setting
we provide, we only consider a single compliance period. In practice, firms
and regulators in these markets will interact over consecutive periods.
Allowing for interaction across these periods through the banking of REC
inventories is an important extension that will bridge the gap between the
existing single period theory, and the multi-period setting that arises in
practice. One clear change that we anticipate is that the rate at which agents
expand their capacity will no longer be curtailed at the end of the first
period as the new capacity can now be carried through time. This allows agents
to “plan ahead” and anticipate their obligations in future periods.

Nonetheless, the problem formulation, algorithm, and example presented in this
work represent a clear step forward within the exciting field of PA-MFGs and
their solutions, as well as their tantalizing applicability to the world of
environmental regulation and beyond.

## References

  * [1] R. Aïd, M. Basei, and H. Pham, A McKean-Vlasov approach to distributed electricity generation development, arXiv preprint arXiv:1705.01302v3, (2019).
  * [2] M. Aziz and P. E. Caines, A mean field game computational methodology for decentralized cellular network optimization, IEEE transactions on control systems technology, 25 (2016), pp. 563–576.
  * [3] C. Beck, E. Weinan, and A. Jentzen, Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations, Journal of Nonlinear Science, 29 (2019), pp. 1563–1619.
  * [4] I. Brown, J. Funk, and R. Sircar, Oil prices & dynamic games under stochastic demand, Available at SSRN 3047390, (2017).
  * [5] P. Cardaliaguet and C.-A. Lehalle, Mean field game of controls and an application to trade crowding, Mathematics and Financial Economics, 12 (2018), pp. 335–363.
  * [6] R. Carmona and F. Delarue, Probabilistic analysis of mean-field games, SIAM Journal on Control and Optimization, 51 (2013), pp. 2705–2734.
  * [7] , Forward–backward stochastic differential equations and controlled mckean–vlasov dynamics, The Annals of Probability, 43 (2015), pp. 2647–2700.
  * [8] , Probabilistic Theory of Mean Field Games with Applications II: Mean Field Games with Common Noise and Master Equations, vol. 84, Springer, 2018.
  * [9] R. Carmona, M. Fehr, and J. Hinz, Optimal stochastic control and carbon price formation, SIAM Journal on Control and Optimization, 48 (2009), pp. 2168–2190.
  * [10] R. Carmona, M. Fehr, J. Hinz, and A. Porchet, Market design for emission trading schemes, Siam Review, 52 (2010), pp. 403–452.
  * [11] R. Carmona and D. Lacker, A probabilistic weak formulation of mean field games and applications, The Annals of Applied Probability, 25 (2015), pp. 1189–1231.
  * [12] R. Carmona and M. Laurière, Convergence analysis of machine learning algorithms for the numerical solution of mean field control and games: Ii–the finite horizon case, arXiv preprint arXiv:1908.01613, (2019).
  * [13] , 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 (2021), pp. 1455–1485.
  * [14] , Deep learning for mean field games and mean field control with applications to finance, arXiv preprint arXiv:2107.04568, (2021).
  * [15] R. Carmona and P. Wang, Finite-state contract theory with a principal and a field of agents, Management Science, 67 (2021), pp. 4725–4741.
  * [16] P. Casgrain and S. Jaimungal, Mean field games with partial information for algorithmic trading, arXiv preprint arXiv:1803.04094, (2018).
  * [17] , Mean-field games with differing beliefs for algorithmic trading, Mathematical Finance, 30 (2020), pp. 995–1034.
  * [18] P. Chan and R. Sircar, Fracking, renewables, and mean field games, SIAM Review, 59 (2017), pp. 588–615.
  * [19] M. Coulon, J. Khazaei, and W. B. Powell, Smart-srec: A stochastic model of the new jersey solar renewable energy certificate market, Journal of Environmental Economics and Management, 73 (2015), pp. 13–31.
  * [20] J. Cvitanić, D. Possamaï, and N. Touzi, Moral hazard in dynamic risk management, Management Science, 63 (2017), pp. 3328–3346.
  * [21] , Dynamic programming approach to principal–agent problems, Finance and Stochastics, 22 (2018), pp. 1–37.
  * [22] B. Djehiche, A. Tcheukam, and H. Tembine, Mean-field-type games in engineering, arXiv preprint arXiv:1605.03281, (2017).
  * [23] R. Elie, T. Mastrolia, and D. Possamaï, A tale of a principal and many, many agents, Mathematics of Operations Research, 44 (2019), pp. 440–467.
  * [24] D. Evangelista and Y. Thamsten, Finite population games of optimal execution, arXiv, (2020).
  * [25] D. Firoozi, Mean field games and optimal execution problems: hybrid and partially observed major minor systems, PhD Thesis, McGill University, 2019\.
  * [26] D. Firoozi and P. E. Caines, The execution problem in finance with major and minor traders: A mean field game formulation, in Annals of the International Society of Dynamic Games (ISDG): Advances in Dynamic and Mean Field Games, vol. 15, Birkhäuser Basel, 2017, pp. 107–130.
  * [27] D. Firoozi, S. Jaimungal, and P. E. Caines, Convex analysis for lqg systems with applications to major–minor lqg mean–field game systems, Systems & Control Letters, 142 (2020), p. 104734.
  * [28] D. Firoozi, A. Pakniyat, and P. E. Caines, A mean field game - hybrid systems approach to optimal execution problems in finance with stopping times, in Proceedings of the 56th IEEE Conference on Decision and Control (CDC), Melbourne, Australia, Dec. 2017, pp. 3144–3151.
  * [29] J.-P. Fouque and Z. Zhang, Deep learning methods for mean field control problems with delay, Frontiers in Applied Mathematics and Statistics, 6 (2020), p. 11.
  * [30] G. Fu, P. Graewe, U. Horst, and A. Popier, A mean field game of optimal portfolio liquidation, arXiv preprint arXiv:1804.04911, (2018).
  * [31] M. Fujii and A. Takahashi, A mean field game approach to equilibrium pricing with market clearing condition, arXiv preprint arXiv:2003.03035, (2020).
  * [32] M. Germain, J. Mikael, and X. Warin, Numerical resolution of mckean-vlasov fbsdes using neural networks, arXiv preprint arXiv:1909.12678, (2019).
  * [33] D. Gomes, L. Lafleche, and L. Nurbekyan, A mean-field game economic growth model, in 2016 American Control Conference (ACC), 2016, pp. 4693–4698.
  * [34] D. Gomes and J. Saúde, A mean-field game approach to price formation in electricity markets, arXiv preprint arXiv:1807.07088, (2018).
  * [35] D. Gomes, R. M. Velho, and M.-T. Wolfram, Socio-economic applications of finite state mean field games, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 372 (2014), pp. 1–17.
  * [36] P. J. Graber and C. Mouzouni, On mean field games models for exhaustible commodities trade, (To appear) ESAIM: COCV, (2019).
  * [37] J. Han and R. Hu, Deep fictitious play for finding markovian nash equilibrium in multi-agent games, in Mathematical and Scientific Machine Learning, PMLR, 2020, pp. 221–245.
  * [38] J. Han, R. Hu, and J. Long, Convergence of deep fictitious play for stochastic differential games, arXiv preprint arXiv:2008.05519, (2020).
  * [39] J. Han, A. Jentzen, and E. Weinan, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences, 115 (2018), pp. 8505–8510.
  * [40] J. Han and J. Long, Convergence of the deep bsde method for coupled fbsdes, Probability, Uncertainty and Quantitative Risk, 5 (2020), pp. 1–33.
  * [41] B. Holmstrom and P. Milgrom, Aggregation and linearity in the provision of intertemporal incentives, Econometrica: Journal of the Econometric Society, (1987), pp. 303–328.
  * [42] M. Huang, P. E. Caines, and R. P. Malhamé, Large-population cost-coupled LQG problems with nonuniform agents: individual-mass behavior and decentralized $\epsilon$-Nash equilibria, IEEE transactions on automatic control, 52 (2007), pp. 1560–1571.
  * [43] M. Huang, R. P. Malhamé, P. E. Caines, et al., Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle, Communications in Information & Systems, 6 (2006), pp. 221–252.
  * [44] X. Huang, S. Jaimungal, and M. Nourian, Mean-field game strategies for optimal execution, Applied Mathematical Finance, 26 (2019), pp. 153–185.
  * [45] J. Khazaei, M. Coulon, and W. B. Powell, Adapt: A price-stabilizing compliance policy for renewable energy certificates: The case of srec markets, Operations Research, 65 (2017), pp. 1429–1445.
  * [46] A. C. Kizilkale, R. Salhab, and R. P. Malhamé, An integral control formulation of mean field game based large scale coordination of loads in smart grids, Automatica, 100 (2019), pp. 312 – 322.
  * [47] J. M. Lasry and P. L. Lions, Jeux à champ moyen. i - le cas stationnaire, Comptes Rendus de l’Académie des Sciences, 343 (2006), pp. 619–625.
  * [48] , Jeux à champ moyen. ii - horizon fini et contrôle optimal, Comptes Rendus de l’Académie des Sciences, 343 (2006), pp. 679–684.
  * [49] J.-M. Lasry and P.-L. Lions, Mean field games, Japanese journal of mathematics, 2 (2007), pp. 229–260.
  * [50] C.-A. Lehalle and C. Mouzouni, A mean field game of portfolio trading and its consequences on perceived correlations, arXiv preprint arXiv:1902.09606, (2019).
  * [51] M. Ludkovski and X. Yang, Mean field game approach to production and exploration of exhaustible commodities, arXiv preprint arXiv:1710.05131, (2017).
  * [52] M. J. Marsden, On uniform spline approximation, Journal of Approximation Theory, 6 (1972), pp. 249–253.
  * [53] Y. Sannikov, A continuous-time version of the principal-agent problem, The Review of Economic Studies, 75 (2008), pp. 957–984.
  * [54] J. Seifert, M. Uhrig-Homburg, and M. Wagner, Dynamic behavior of co2 spot prices, Journal of Environmental Economics and Management, 56 (2008), pp. 180–194.
  * [55] A. Shrivats, D. Firoozi, and S. Jaimungal, A mean-field game approach to equilibrium pricing, optimal generation, and trading in solar renewable energy certificate (srec) markets, arXiv preprint arXiv:2003.04938, (2020).
  * [56] A. Shrivats and S. Jaimungal, Optimal generation and trading in solar renewable energy certificate (srec) markets, Applied Mathematical Finance, 27 (2020), pp. 99–131.
  * [57] E. Weinan, J. Han, and A. Jentzen, Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Communications in Mathematics and Statistics, 5 (2017), pp. 349–380.