# A deep learning method for solving stochastic optimal control problems
driven by fully-coupled FBSDEs

Shaolin Ji1 Shige Peng2 Ying Peng1 and Xichuan Zhang2

1Zhongtai Securities Institute for Financial Studies, Shandong University,
250100, China
E-mail: {jsl<EMAIL_ADDRESS>
2School of Mathematics, Shandong University, 250100, China
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we mainly focus on the numerical solution of high-dimensional
stochastic optimal control problem driven by fully-coupled forward-backward
stochastic differential equations (FBSDEs in short) through deep learning. We
first transform the problem into a stochastic Stackelberg differential
game(leader-follower problem), then a cross-optimization method (CO method) is
developed where the leader’s cost functional and the follower’s cost
functional are optimized alternatively via deep neural networks. As for the
numerical results, we compute two examples of the investment-consumption
problem solved through stochastic recursive utility models, and the results of
both examples demonstrate the effectiveness of our proposed algorithm.

Keywords stochastic optimal control, FBSDEs, deep learning, Stackelberg
differential game, recursive utility

##  1 Introduction

Bismut [Bis73] first introduced linear backward stochastic differential
equations (BSDEs in short) as the adjoint equation of the classical stochastic
optimal control problem. In 1990, Pardoux and Peng firstly proved the
existence and uniqueness of nonlinear BSDEs with Lipschitz condition [PP90].
Since then, the theory of BSDEs has been studied by many researchers and
applied in a wide range of areas, such as in stochastic optimal control and
mathematical finance [EKPQ97, Pen90]. When a BSDE is coupled with a (forward)
stochastic differential equation (SDE in short), the system is usually called
a forward-backward stochastic differential equation (FBSDE in short). We can
refer to the literatures in [PW99, Ant93, CM96, HP95, PT99, MPY94] which
studied the existence, uniqueness and the applications of coupled or fully-
coupled FBSDEs.

In 1993, Peng [Pen93] first established a local stochastic maximum principle
(SMP) for stochastic optimal control problems driven by FBSDEs. Then, the
local SMP for other related problems were studied by Dokuchaev and Zhou
[DZ99], Ji and Zhou [JZ06], Yong and Zhou [YZ99] ( see also the references
therein). When the control domain is non-convex, Yong [Yon10] studied a fully
coupled controlled FBSDE with mixed initial-terminal conditions. Wu [Wu13]
studied a stochastic recursive optimal control problem. Hu [Min17], Hu et al.
[HJX18] built the global stochastic maximum principle for stochastic optimal
control problems driven by BSDEs and FBSDEs respectively. On the other hand,
the dynamic programming principle (DPP) and related Hamilton-Jacobi-Bellman
(HJB) equations have been intensively studied by Li and Wei [LW14], Hu et al.
[HJX19] for this kind of stochastic optimal control problems. Furthermore, Hu
et al. [HJX20] revealed the relationship between the SMP and the DPP for a
stochastic optimal control problem where the system is governed by a fully
coupled FBSDE. However, few literature has studied the numerical method for
solving the stochastic optimal control problems driven by FBSDEs.

In this paper, we aim to solve the high-dimensional stochastic optimal control
problem where the state equation is governed by

$\left\\{\begin{array}[]{l}dx(t)=b(t,x(t),y(t),z(t),u(t))dt+\sigma(t,x(t),y(t),z(t),u(t))dB(t),\vspace{1ex}\\\
-dy(t)=l(t,x(t),y(t),z(t),u(t))dt-z(t)dB(t),\vspace{1ex}\\\ x(0)=a,\qquad
y(T)=g(x(T)).\end{array}\right.$ (1.1)

and the cost functional is defined by

$J(u(\cdot))=E\left[\int_{0}^{T}f(t,x(t),y(t),z(t),u(t))dt+h(x(T))+\gamma(y(0))\right].$
(1.2)

Recall that for the high-dimensional classical stochastic optimal control
problem where the state equation is a SDE, Han and E [HE16] solved it through
deep learning method. In more details, they developed a feed-forward neural
network to approximate the control $u(\cdot)$ and regard the cost functional
$J(u(\cdot))$ as the optimization objective (see also [CL]). Comparing with
the classical stochastic optimal control problem, the difficulty of our
problem lies in that we need to deal with the high-dimensional FBSDE (1.1)
firstly. And it is well-known that to obtain the numerical solution of (1.1)
itself is a difficult problem.

The traditional methods for solving the FBSDEs include the partial
differential equation (PDE in short) methods and the probabilistic methods,
such as [Tad12, BT04, JJZ08, BZ08, FZZ16, MT04, FZT16, HRO16]. However, most
of these methods can not deal with high-dimensional FBSDEs. In this paper, we
reformulate the fully-coupled FBSDE (1.1) as a stochastic optimal control
problem with the cost functional given by the error between the terminal
condition and the solution of the FBSDE,

$J(y_{0},z(\cdot))=E\left[|y(T)-g(x(T))|^{2}\right],$

and build two feed-forward neural networks to approximate the control
variables $(y_{0},z(\cdot))$. This idea for solving (1.1) is currently a
common method in solving high-dimensional BSDEs and FBSDEs (see e.g. [HJE18,
EHJ17, HL18, JPPZ20, HPW20, PWG21, BEJ19]).

Based on the above analysis, we put forward the following ideas to solve our
stochastic optimal control problem (1.1)-(1.2) by a novel deep learning
method. Firstly, we transform (1.1)-(1.2) into a stochastic Stackelberg
differential game(leader-follower problem). The goal of the follower is to
find a pair of optimal control $(y_{0},z(\cdot))$ which minimizes the
functional $J(y_{0},z(\cdot))$ under a given control $u(\cdot)$, while the
goal of the leader is to find an optimal control $u(\cdot)$ which minimizes
the cost functional $J(u(\cdot))$. Secondly, a cross-optimization method is
introduced to solve the high-dimensional stochastic Stackelberg differential
game.

Specifically, we construct a feed-forward neural network to approximate the
control $u(\cdot)$ and regard the cost functional $J(u(\cdot))$ as the
optimization objective for the leader. We also build two feed-forward neural
networks to approximate the control $(y_{0},z(\cdot))$ and regard
$J(y_{0},z(\cdot))$ as the optimization objective for the follower. As for the
parameter updates of the neural networks, we update the parameters
approximating $(y_{0},z(\cdot))$ through $J(y_{0},z(\cdot))$ and the
parameters approximating $u(\cdot)$ through $J(u(\cdot))$ alternatively. In
more details, in each training step, we first fix an approximated $u(\cdot)$
to perform $\kappa$ $(\kappa\geq 1,\kappa\in\mathbb{N}^{+})$ times updates for
the parameters approximating $(y_{0},z(\cdot))$, then based on the obtained
$(y_{0},z(\cdot))$, we perform one time update for the parameters
approximating $u(\cdot)$. The training will be repeated until a convergence
result is obtained. Note that under a given $u(\cdot)$, the follower should
choose the optimal $(y_{0},z(\cdot))$ such that $J(y_{0},z(\cdot))=0$. Here in
order to improve the computation efficiency, we use a relaxation method in the
parameter update process. That is, instead of forcing $J(y_{0},z(\cdot))=0$,
we control the number of updates by choosing an enough large coefficient
$\kappa$ so that $J(y_{0},z(\cdot))$ can obtain a sufficiently small value
which is close to $0$, but not strictly equal to $0$. We call this method the
cross-optimization method (CO method).

To show the feasibility of our proposed method, we compute the investment-
consumption problem for stochastic recursive utilities in a financial market.
The concept of stochastic recursive utility was first introduced by Duffie and
Epstein [DE92]. In fact, the stochastic recursive utility is associated with
the solution of a particular backward stochastic differential equation (BSDE).
From the BSDE point of view, El Karoui et al. [EKPQ97, EKPQ01] considered a
more general class of recursive utilities defined as the solutions of BSDEs.
In this paper, we compute two examples. In the first example, the stochastic
recursive utility is described by a linear BSDE. For this case, the
investment-consumption problem can be transformed to a classical stochastic
optimal control problem. In this way, we can compare our computation results
with those obtained by the numerical methods (such as the method in [HE16])
for solving classical stochastic optimal control problems. The numerical
results show that the value of the recursive utility functions obtained by the
two methods are very close. In the second example, we consider a more general
recursive utility problem whose generator contains the $z$ term (refer to
[EKPQ97, EKPQ01, CE02]). Chen and Epstein [CE02] studied the recursive utility
problem with $z$ term, in which $z$ is regarded as the volatility of utility.
The computation results show that our proposed method is also effective for
this general case.

This paper is organized as follows. In section 2, we describe the stochastic
optimal control problem driven by fully-coupled FBSDEs and reformulate it as a
stochastic Stackelberg differential game problem. In section 3, we present our
proposed CO method for solving the stochastic Stackelberg differential game
via deep learning. As the numerical results, we compute the investment-
consumption problem for stochastic recursive utilities in section 4. In
section 5, we make a brief conclusion.

##  2 Statement of the problem

###  2.1 The stochastic control problem driven by FBSDE

On a given complete probability space $(\Omega,\mathcal{F},\mathbb{P})$, let
$(B(t),t\geq 0)$ be a standard $d$-dimensional Brownian motion defined on a
finite interval $[0,T]$ for some given constant $T>0$. We denote

$\mathcal{F}_{t}=\sigma\\{B(s),0\leq s\leq t\\}$

the natural filtration of $\\{B(t)\\}_{0\leq t\leq T}$ and
$\mathbb{F}:=\\{\mathcal{F}_{t}\\}_{0\leq t\leq T}$. $\mathcal{F}_{0}$
contains all the $\mathbb{P}$-null set in $\mathcal{F}$ and
$\\{\mathcal{F}_{t}\\}_{t\geq 0}$ is right continuous.
$M^{2}(0,T;\mathbb{R}^{n})$ is denoted as the space of all mean square-
integrable $\mathcal{F}_{t}$-adapted and $\mathbb{R}^{n}$-valued processes. It
is a Hilbert space with the norm

$\|v(\cdot)\|_{M^{2}(0,T;\mathbb{R}^{n})}=\Big{(}\mathbb{E}\big{[}\int_{0}^{T}|v(t)|^{2}dt\big{]}\Big{)}^{1/2}.$

We also denote $L^{2}(\Omega,\mathcal{F}_{t},\mathbb{P})$ as

$L^{2}(\Omega,\mathcal{F}_{t},\mathbb{P})\triangleq\left\\{\xi|\xi\in\mathbb{R}^{n}\mbox{
is }\mathcal{F}_{t}\mbox{-measurable and
}\mathbb{E}\left[|\xi|^{2}\right]<\infty\right\\}.$

$|x|$ is denoted as the Euclidean norm of an element $x\in\mathbb{R}^{n}$ and
$\left\langle x,y\right\rangle$ is denoted as the Euclidean inner product of
elements $x,y\in\mathbb{R}^{n}$.

In this paper, we consider the following controlled fully-coupled FBSDE,

$\left\\{\begin{array}[]{l}dx(t)=b(t,x(t),y(t),z(t),u(t))dt+\sigma(t,x(t),y(t),z(t),u(t))dB(t),\vspace{1ex}\\\
-dy(t)=l(t,x(t),y(t),z(t),u(t))dt-z(t)dB(t),\vspace{1ex}\\\ x(0)=a,\qquad
y(T)=g(x(T)),\end{array}\right.$ (2.1)

where

$\displaystyle b:$
$\displaystyle\quad[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times
d}\times U\rightarrow\mathbb{R}^{n},$ $\displaystyle\sigma:$
$\displaystyle\quad[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times
d}\times U\rightarrow\mathbb{R}^{n\times d},$ $\displaystyle l:$
$\displaystyle\quad[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times
d}\times U\rightarrow\mathbb{R}^{m},$ $\displaystyle g:$
$\displaystyle\quad\mathbb{R}^{n}\rightarrow\mathbb{R}^{m},$

are given $C^{1}$ functions with respect to $(x,y,z,u)$. The process
$u(\cdot)$ taking value in a given nonempty convex set
$U\subseteq\mathbb{R}^{k}$ in the system (2.1) is called an admissible
control. For a given admissible control $u(\cdot)$, the corresponding solution
of the system (2.1) is denoted by

$(x(t),y(t),z(t))=(x^{u}(t),y^{u}(t),z^{u}(t)),\quad\mbox{ for all}\quad
t\in[0,T].$

The set of all admissible control $u(\cdot)$ is denoted as
$\mathcal{U}_{ad}[0,T]$, and the corresponding 4-tuple $(x(t),y(t),z(t),u(t))$
is called an admissible 4-tuple.

For a given control $u(\cdot)$, $b,\sigma,l,g$ are functions of $(t,x,y,z)$.
For these functions, we need to introduce some assumptions.

Given an $m\times n$ full-rank matrix $G$ and set

$w=\begin{pmatrix}x\\\ y\\\
z\end{pmatrix}\in\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times
d},\qquad A(t,w)=\begin{pmatrix}-G^{T}l\\\ Gb\\\ G\sigma\end{pmatrix}(t,w),$

where $G\sigma=(G\sigma_{1},G\sigma_{2},\cdots,G\sigma_{d})$, and

$\left\langle w^{1},w^{2}\right\rangle=\left\langle
x^{1},x^{2}\right\rangle+\left\langle y^{1},y^{2}\right\rangle+\left\langle
z^{1},z^{2}\right\rangle.$

######  Assumption 2.1.

  * (i)

$A(t,w)$ is uniformly Lipschitz with respect to $w$;

  * (ii)

for each $w\in\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times d}$,
$A(\cdot,w)\in M^{2}(0,T;\mathbb{R}^{n})$;

  * (iii)

$g(x)$ is uniformly Lipschitz with respect to $x\in\mathbb{R}^{n}$.

The following monotonic conditions firstly introduced in [PW99] are also
needed.

######  Assumption 2.2.

$\displaystyle\left\langle A(t,w)-A(t,w^{\prime}),w-w^{\prime}\right\rangle$
$\displaystyle\leq-\beta_{1}|G\hat{x}|^{2}-\beta_{2}(|G^{T}\hat{y}|^{2}+|G^{T}\hat{z}|^{2}),$
(2.2) $\displaystyle\left\langle
g(x)-g(x^{\prime}),G(x-x^{\prime})\right\rangle$
$\displaystyle\geq\mu_{1}|G\hat{x}|^{2},$ $\displaystyle\forall w=(x,y,z),\
w^{\prime}=(x^{\prime},y^{\prime},z^{\prime}),$ $\displaystyle\
\hat{x}=x-x^{\prime},\ \hat{y}=y-y^{\prime},\ \hat{z}=z-z^{\prime},$

or

$\displaystyle\left\langle A(t,w)-A(t,w^{\prime}),w-w^{\prime}\right\rangle$
$\displaystyle\geq\beta_{1}|G\hat{x}|^{2}+\beta_{2}(|G^{T}\hat{y}|^{2}+|G^{T}\hat{z}|^{2}),$
(2.3) $\displaystyle\left\langle
g(x)-g(x^{\prime}),G(x-x^{\prime})\right\rangle$
$\displaystyle\leq-\mu_{1}|G\hat{x}|^{2},$ $\displaystyle\forall w=(x,y,z),\
w^{\prime}=(x^{\prime},y^{\prime},z^{\prime}),$ $\displaystyle\
\hat{x}=x-x^{\prime},\ \hat{y}=y-y^{\prime},\ \hat{z}=z-z^{\prime},$

where $\beta_{1},\beta_{2}$ and $\mu_{1}$ are non-negative constants with
$\beta_{1}+\beta_{2}>0$, $\beta_{2}+\mu_{1}>0$. Moreover we have
$\beta_{1}>0,\mu_{1}>0$ (resp. $\beta_{2}>0$) when $m>n$ (resp. $m<n$).

We also introduce the following Lemma 1 without proof, and the proof of this
lemma can be found in [PW99, Wu98].

######  Lemma 1.

For any given admissible control $u(\cdot)$, let Assumption 2.1 and Assumption
2.2 hold. Then the FBSDE (2.1) has a unique adapted solution
$(x(\cdot),y(\cdot),z(\cdot))$.

In the following, we define the cost functional of the stochastic optimal
control problem by

$J(u(\cdot)):=E\left[\int_{0}^{T}f(t,x(t),y(t),z(t),u(t))dt+h(x(T))+\gamma(y(0))\right],$
(2.4)

where

$\displaystyle f:$
$\displaystyle\quad[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m\times
d}\times U\rightarrow\mathbb{R},$ $\displaystyle h:$
$\displaystyle\quad\mathbb{R}^{n}\rightarrow\mathbb{R},$
$\displaystyle\gamma:$
$\displaystyle\quad\mathbb{R}^{m}\rightarrow\mathbb{R},$

are given $C^{1}$ functions with respect to $(x,y,z,u)$. The 4-tuple process
$(x(t),y(t),z(t),u(t))$ satisfies equation (2.1). Then the stochastic control
problem driven by FBSDE can be described as following.

######  Problem 1.

The optimal control problem is to find an admissible control $u^{*}(\cdot)$
over $\mathcal{U}_{ad}[0,T]$ such that

$J(u^{*}(\cdot))=\inf_{u(\cdot)\in\mathcal{U}_{ad}[0,T]}J(u(\cdot)).$ (2.5)

If there exists a control $u^{*}(\cdot)$ which minimizes $J(u(\cdot))$ over
$\mathcal{U}_{ad}[0,T]$, then we call it an optimal control. (2.1) is called
the optimal state equation and its solution
$(x^{*}(\cdot),y^{*}(\cdot),z^{*}(\cdot))$ is called an optimal trajectory.

######  Remark.

In Problem 1 and Lemma 1, we study the stochastic control problem driven by
fully-coupled FBSDE (2.1), which needs some strong conditions (such as the
monotonic conditions Assumption 2.1 and Assumption 2.2). If the forward SDE in
(2.1) does not contain $y$ and $z$ terms, such strong conditions can be
relaxed [Pen93].

###  2.2 The Equivalent problem

As discussed in [KZ00, LZ01, Yon10, Wu13], $z$ can be regarded as a control
variable and $y(T)=g(x(T))$ in (2.1) can be seen as a terminal state
constraint. Then Problem 1 can be reformulated as the following stochastic
control problem to minimize

$J(u(\cdot),y_{0},z(\cdot))=E\left[\int_{0}^{T}f(t,x(t),y(t),z(t),u(t))dt+h(x(T))+\gamma(y(0))\right],$
(2.6)

over $u(\cdot)\in\mathcal{U}_{ad}[0,T],y_{0}\in\mathbb{R}^{m},z(\cdot)\in
M^{2}(0,T;\mathbb{R}^{m\times d})$, the state equation is given as

$\left\\{\begin{array}[]{l}dx(t)=b(t,x(t),y(t),z(t),u(t))dt+\sigma(t,x(t),y(t),z(t),u(t))dB(t),\vspace{1ex}\\\
-dy(t)=l(t,x(t),y(t),z(t),u(t))dt-z(t)dB(t),\vspace{1ex}\\\ x(0)=a,\qquad
y(0)=y_{0},\end{array}\right.$ (2.7)

with the terminal state constraint

$E\left[|y(T)-g(x(T))|^{2}\right]=0.$ (2.8)

In this paper, our idea is to relax the terminal state constraint as an
optimization objective which must be met firstly. Following this idea, the
above stochastic control problem (2.6) becomes a stochastic Stackelberg
differential game problem (leader-follower problem). In more details, we
assume that there is one follower and one leader, and their cost functionals
are given as follows

$\displaystyle J_{1}(u(\cdot),y_{0},z(\cdot))$
$\displaystyle:=E\left[|y(T)-g(x(T))|^{2}\right],$ (2.9) $\displaystyle
J_{2}(u(\cdot),y_{0},z(\cdot))$
$\displaystyle:=E\left[\int_{0}^{T}f(t,x(t),y(t),z(t),u(t))dt+h(x(T))+\gamma(y(0))\right].$

For any choice $u(\cdot)$ of the leader and a fixed initial state $a$, the
goal of the follower is to minimize the functional
$J_{1}(u(\cdot),y_{0},z(\cdot))$ over $(y_{0},z(\cdot))\in\mathcal{A}$ where
$\mathcal{A}$ is the set of all admissible $(y_{0},z(\cdot))$. Denote the
optimal $(y_{0},z(\cdot))$ of the follower as $(y_{0}^{u},z^{u}(\cdot))$ which
clearly depends on $u(\cdot)$ of the leader. Then, the goal of the leader is
to minimize the cost functional $J_{2}(u(\cdot),y_{0}^{u},z^{u}(\cdot))$ over
$u(\cdot)\in\mathcal{U}_{ad}[0,T]$.

In a more rigorous way, for $b,\sigma,l,g$ in (2.7) and $f,h,\gamma$ in (2.9),
we introduce the following additional assumptions.

######  Assumption 2.3.

  * (i)

$b,\sigma,l,g,f,h$ and $\gamma$ are continuously differential;

  * (ii)

the derivatives of $b,\sigma,l,g$ are bounded;

  * (iii)

the derivatives of $f$ are bounded by $C(1+|x|+|y|+|z|+|u|)$;

  * (iv)

the derivatives of $h$ and $\gamma$ are bounded by $C(1+|x|)$ and $C(1+|y|)$,
respectively.

Then the optimal control problem (2.5) (Problem 1) can be reformulated as the
following Stackelberg differential game.

######  Problem 2.

To find a map
$\hat{\alpha}(\cdot):\mathcal{U}_{ad}[0,T]\rightarrow\mathcal{A}$ and a
control $u^{*}(\cdot)$ such that

$\begin{cases}J_{1}(u(\cdot),\hat{\alpha}(u(\cdot))(\cdot))=\inf\limits_{(y_{0},z(\cdot))\in\mathcal{A}}J_{1}(u(\cdot),y_{0},z(\cdot))\qquad\forall
u(\cdot)\in\mathcal{U}_{ad}[0,T],\vspace{1ex}\\\
J_{2}(u^{*}(\cdot),\hat{\alpha}(u^{*}(\cdot))(\cdot))=\inf\limits_{u(\cdot)\in\mathcal{U}_{ad}[0,T]}J_{2}(u(\cdot),\hat{\alpha}(u(\cdot))(\cdot)).\end{cases}$
(2.10)

The state equation is given as (2.7).

Under Assumption 2.3, for any given $(u(\cdot),y_{0},z(\cdot))$ and the
initial state $a$, there exists a strong solution $(x(\cdot),y(\cdot))$ to
(2.7). If the above optimal pair $(\hat{\alpha}(\cdot),u^{*}(\cdot))$ in
Problem 2 exists, we have
$(y_{0}^{*},z^{*}(\cdot))=\hat{\alpha}(u^{*}(\cdot))(\cdot)$ and
$u^{*}(\cdot)$ is an optimal control of Problem 1. The advantage of Problem 2
is that it does not need the regularity/integrability on $z(\cdot)$, and the
state equation is a forward SDE instead of a backward SDE as in Problem 1. The
difficulty for solving Problem 2 is that it has to treat another optimization
goal $E\left[|y(T)-g(x(T))|^{2}\right]$.

For simplicity, we denote the leader’s cost functional
$J_{2}(u(\cdot),y_{0},z(\cdot))$ by $J(u(\cdot))$, and the follower’s cost
functional $J_{1}(u(\cdot),y_{0},z(\cdot))$ by $J(y_{0},z(\cdot))$ in the
following of the paper.

##  3 Deep neural network for solving the stochastic Stackelberg differential
game

In this section, we consider to solve Problem 2 with deep neural network. As
shown in section 2, Problem 2 is an optimal control problem with forward state
equations, it contains two optimization goals, and the follower’s cost
functional $J(y_{0},z(\cdot))$ should be optimized in priority to the leader’s
cost functional. In order to solve Problem 2, instead of using the penalty
method which multiply the follower’s cost functional $J(y_{0},z(\cdot))$ by a
penalty parameter and adding it to the leader’s cost functional $J(u(\cdot))$,
we propose a cross-optimization method (CO method in brief). We optimize these
two cost functionals alternatively and spend more computational cost on the
follower’s cost functional $J(y_{0},z(\cdot))$. In this way, on the one hand,
Problem 2 become unconstrained with respect to the control. On the other hand,
we do not have to choose the approximate penalty parameter, which is usually
hard to choose in the penalty method. The detailed process for optimizing
these two functionals can be found in subsection 3.2.

Before showing the approximation algorithm of Problem 2, we firstly need to
discretize the forward state equation (2.7). Let $N\in\mathbb{N}^{+}$ and
$0=t_{0}<t_{1}<t_{2}<\cdots<t_{N-1}<t_{N}=T$ so that

$\delta=\max_{0\leq i\leq N-1}|t_{i+1}-t_{i}|,$

is sufficiently small. Define $\Delta t_{i}=t_{i+1}-t_{i}$ and $\Delta
B_{t_{i}}=B(t_{i+1})-B(t_{i})$, where $B(t_{i})\sim\mathcal{N}(0,t_{i})$, for
$i=0,1,2,\cdots,N-1$. We apply the Euler-Maruyama scheme to the forward state
equation (2.7), then we have

$\left\\{\begin{array}[]{l}x_{t_{i+1}}=x_{t_{i}}+b(t_{i},x_{t_{i}},y_{t_{i}},z_{t_{i}},u_{t_{i}})\Delta
t_{i}+\sigma(t_{i},x_{t_{i}},y_{t_{i}},z_{t_{i}},u_{t_{i}})\Delta
B_{t_{i}},\vspace{1ex}\\\
y_{t_{i+1}}=y_{t_{i}}-l(t_{i},x_{t_{i}},y_{t_{i}},z_{t_{i}},u_{t_{i}})\Delta
t_{i}+z_{t_{i}}\Delta B_{t_{i}},\vspace{1ex}\\\ x_{t_{0}}=a,\qquad
y_{t_{0}}=y_{0}.\end{array}\right.$ (3.1)

For easier expression, we denote $x(t_{i})$ as $x_{t_{i}}$ here.

We use Monte Carlo sampling to approximate the expectations in the cost
functionals. Then, the leader’s cost functional (2.4) can be evaluated by the
following schemes:

$J(u(\cdot))=\dfrac{1}{M}\sum_{m=1}^{M}\left[\sum_{i=0}^{N-1}f(t_{i},x_{t_{i}}^{m},y_{t_{i}}^{m},z_{t_{i}}^{m},u_{t_{i}}^{m})\Delta
t_{i}+h(x_{t_{N}}^{m})+\gamma(y_{t_{0}}^{m})\right],$ (3.2)

where $M$ represents the number of Monte Carlo samples. Similarly, the
follower’s cost functional $J(y_{0},z(\cdot))$ can be approximated by

$J(y_{0},z(\cdot))=\dfrac{1}{M}\sum_{m=1}^{M}\left[|y_{t_{N}}^{m}-g(x_{t_{N}}^{m})|^{2}\right],$
(3.3)

Here the numbers of Monte Carlo samples $M$ in (3.2) and (3.3) may be
different. Without causing confusion, we use the same notation $M$ to
represent the Monte Carlo sampling number, and the same notation $J(u(\cdot))$
and $J(y_{0},z(\cdot))$ as in the continuous form to represent the
corresponding discrete cost functionals.

###  3.1 Neural network architecture

In this subsection, we present the neural network architecture for solving
Problem 2 approximately. Let $q\in\mathbb{N}^{+}$ be a sufficiently large
natural number, and
$\mathcal{NN}^{\theta_{u}}:[0,T]\times\mathbb{R}^{n}\mapsto U$ be a Borel
measurable function.

Now we let the function $\mathcal{NN}^{\theta_{u}}$ with suitable
$\theta_{u}\in\mathbb{R}^{q}$ be the approximation of the control $u$:

$\displaystyle u_{t_{i}}$
$\displaystyle=\mathcal{NN}^{\theta_{u}}(t_{i},x_{t_{i}};\theta_{u}),$ (3.4)

for $i\in\\{0,1,\cdots,N-1\\}$. $\theta_{u}$ is the trainable parameter in
this neural network. We represent the function $\mathcal{NN}^{\theta_{u}}$
with a multilayer feed-forward neural network of the form

$\mathcal{NN}^{\theta_{u}}=\psi\circ\mathcal{A}_{I}\circ\sigma_{I-1}\circ\mathcal{A}_{I-1}\circ\cdots\circ\sigma_{1}\circ\mathcal{A}_{1},$
(3.5)

where

  * •

$I$ is a positive integer specifying the depth of the neural network,

  * •

$\mathcal{A}_{1},\cdots,\mathcal{A}_{I}$ are functions of the form

$\mathcal{A}_{i}=w_{i}\alpha_{i}+b_{i}\in\mathbb{R}^{d_{i}},\qquad\text{for
}1\leq i\leq I,\qquad\alpha_{i}\in\mathbb{R}^{d_{i-1}},$

where the matrix weights $w_{i}$ and the bias vector $b_{i}$ are trainable
parameters such that $\theta_{u}=(w_{i},b_{i})_{1\leq i\leq I}$, $d_{i}$ is
the number of nodes in layer $i$, and $\alpha_{0}$ represents the inputs of
the neural network;

  * •

$\sigma_{I-1},\cdots,\sigma_{1}$ are the nonlinear activation functions, such
as sigmoid, ReLU, ELU, etc.;

  * •

$\psi:\mathbb{R}^{d_{I}}\mapsto U$ is a given function, in this paper, we set
$d_{I}=k$ and $U=\mathbb{R}^{k}$.

We treat the time $t$ as a part of input variables and the dimensions of the
input and output in the neural network are $(n+1)$ and $k$, respectively.
Figure 1 shows an example of a single neural network for $n=2$, $m=k=1$.

Figure 1: Representation of a single neural network with $I=4$,
$d_{1}=d_{2}=d_{3}=4$, $n=2$, $d_{4}=k=1$.

In this paper, as the time $t$ is treated as an input variable, we use a
common neural network, i.e. we share the parameters of the neural network
among all the time points. The neural network is constructed with 5 layers,
including 1 input layer with $(n+1)$ neurons, 3 hidden layers with $(n+10)$
neurons and 1 output layer with $k$ neurons. We adopt ReLU as the activation
functions through the network, and add a batch normalization layer for each
layer.

For approximating the controls $y_{0}$ and $z(\cdot)$, we also construct two
different feed-forward neural networks $\mathcal{NN}^{\theta_{y}}$ and
$\mathcal{NN}^{\theta_{z}}$, respectively.

$\left\\{\begin{array}[]{l}y_{0}=\mathcal{NN}^{\theta_{y}}(x_{0};\theta_{y})\vspace{1ex}\\\
z_{t_{i}}=\mathcal{NN}^{\theta_{z}}(t_{i},x_{t_{i}};\theta_{z}).\end{array}\right.$
(3.6)

The architectures of $\mathcal{NN}^{\theta_{y}}$ and
$\mathcal{NN}^{\theta_{z}}$ are similar to that of
$\mathcal{NN}^{\theta_{u}}$, and the trainable parameters of them are
$\theta_{y}$ and $\theta_{z}$, respectively. They have the same depth, the
same dimensions for the hidden layers, the same activation functions with the
neural network $\mathcal{NN}^{\theta_{u}}$, but different input and output
dimensions. The dimensions of the input and output for
$\mathcal{NN}^{\theta_{y}}$ are $n$ and $m$, while that of the input and
output for $\mathcal{NN}^{\theta_{z}}$ are $n+1$ and $m\times d$. We also add
a batch normalization layer for each layer in the two neural networks. The
whole neural network structure is shown in Figure 2.

Figure 2:  For easier notation, here we use $x_{i}$ to represent $x_{t_{i}}$,
and $y_{i}$, $z_{i}$, $u_{i}$ for $y_{t_{i}}$, $z_{t_{i}}$, $u_{t_{i}}$
respectively. In the whole neural network structure, $h_{y}$, $h_{z}$ and
$h_{u}$ represent the hidden layers of $\mathcal{NN}^{\theta_{y}}$,
$\mathcal{NN}^{\theta_{z}}$ and $\mathcal{NN}^{\theta_{u}}$, respectively. And
a common neural network is used for all the discrete time points.

###  3.2 Updating the parameters in the neural networks

From the formulation of Problem 2, two objective functionals $J(u(\cdot))$ and
$J(y_{0},z(\cdot))$ should be optimized, and the follower’s cost functional
$J(y_{0},z(\cdot))$ should be optimized in priority to the leader’s cost
functional. Which means that we should make much more effort on optimizing the
follower’s cost functional $J(y_{0},z(\cdot))$. A commonly used method for
solving this problem is the penalty method, the idea of which is to treat
$J(y_{0},z(\cdot))$ as a penalty and add it to the leader’s cost functional
$J(u(\cdot))$:

$\bar{J}(y_{0},z(\cdot),u(\cdot))=J(u(\cdot))+\mu J(y_{0},z(\cdot)),$ (3.7)

where $\mu$ is a sufficiently large penalty parameter. And the aim is to find
the optimal control $(y^{*}_{0},z^{*}(\cdot),u^{*}(\cdot))$ of (3.7), in this
way, Problem 2 is transformed to a classical stochastic optimal control
problem without state constraint.

In this paper, different from the penalty method (3.7), we develop a novel
deep learning method (the CO method) for solving this kind of stochastic
Stackelberg differential game. In the neural network architecture, we update
the network parameters through cross-optimization of the two cost functionals
(3.2) and (3.3). For a given network parameter $\theta_{u}$, we update the
network parameters $(\theta_{y},\theta_{z})$ through the discrete follower’s
cost functional (3.3). Then, under the approximation of $y_{0},z(\cdot)$ with
the updated $(\theta_{y},\theta_{z})$, we update the network parameters
$\theta_{u}$ through the discrete leader’s cost functional (3.2). In other
words, in each training step, the network parameters $\theta_{u}$ for
optimizing (3.2) and $(\theta_{y},\theta_{z})$ for optimizing (3.3) are
updated alternatively. As mentioned above, the minimum of the follower’s cost
functional $J(y_{0},z(\cdot))$ should be satisfied in priority to the leader’s
cost functional. Therefore in each training step, the update times of
$(\theta_{y},\theta_{z})$ will be much higher than that of $\theta_{u}$ in
practice.

Formally speaking, for a certain $\kappa\in\mathbb{N}^{+}$, a training step
contains two different sub-steps

$\left\\{\begin{array}[]{l}\theta^{l+1}_{u}=\Xi_{1}(\theta^{l}_{u},\nabla
J(u(\cdot))),\vspace{1ex}\\\
(\theta^{l+1}_{y},\theta^{l+1}_{z})=\Xi_{2}\circ\cdots\circ\Xi_{2}\circ\Xi_{2}(\theta^{l}_{y},\theta^{l}_{z},\nabla
J(y_{0},z(\cdot))),\quad\mbox{total $\kappa$ times,}\end{array}\right.$ (3.8)

where $\Xi_{1}(\theta^{l}_{u},\nabla J(u(\cdot)))$ and
$\Xi_{2}(\theta^{l}_{y},\theta^{l}_{z},\nabla J(y_{0},z(\cdot)))$ are given
functions of parameters and their gradients. They represent the parameter
update operations for the cost functionals (3.2) and (3.3), respectively. When
the parameter $\theta^{l}_{u}$ is updated through
$\Xi_{1}(\theta^{l}_{u},\nabla J(u(\cdot)))$, parameters
$(\theta^{l}_{y},\theta^{l}_{z})$ remain unchanged, and when parameters
$(\theta^{l}_{y},\theta^{l}_{z})$ are updated through
$\Xi_{2}^{1}(\theta^{l}_{y},\theta^{l}_{z},\nabla J(y_{0},z(\cdot)))$, the
parameter $\theta^{l}_{u}$ remains unchanged. Here we call $\kappa$ the
penalty updating coefficient, which represents the ratio of the number of
optimizations of $(\theta^{l}_{y},\theta^{l}_{z})$ to that of
$\theta^{l}_{u}$.

In the actual calculation process, since the updates of $\theta^{l}_{u}$ and
$(\theta^{l}_{y},\theta^{l}_{z})$ are carried out alternatively, it is
feasible to perform the first update for $\theta^{l}_{u}$ or
$(\theta^{l}_{y},\theta^{l}_{z})$. And we only need to ensure that $\kappa$ is
large. The pseudo code of the proposed algorithm is given in Algorithm 1.

Algorithm 1 The CO algorithm for solving stochastic control problem driven by
fully coupled FBSDE

1:The Brownian motion $\Delta B_{t_{i}}$, initial state $a$ and parameter
$(\theta_{y}^{0},\theta_{z}^{0},\theta_{u}^{0})$ and hyper-parameter
$\kappa^{\prime}$, etc.;

2:The precesses
$(x^{\ell}_{t_{i}},y^{\ell}_{t_{i}},z^{\ell}_{t_{i}},u^{\ell}_{t_{i}})$, etc.

3:for $\ell=0$ to $maxstep$ do

4: $x_{0}^{\ell,m}=a$, $loss=0$;

5:
$y_{0}^{\ell,m}=\mathcal{NN}^{\theta_{y}}(x_{0}^{\ell,m};\theta_{y}^{\ell})$;

6: $loss=loss+\gamma(y_{0}^{\ell,m})$

7: for $i=0$ to $N-1$ do

8:
$u^{\ell,m}_{t_{i}}=\mathcal{NN}^{\theta_{u}}(t_{i},x^{\ell,m}_{t_{i}};\theta_{u}^{\ell});$

9:
$z^{\ell,m}_{t_{i}}=\mathcal{NN}^{\theta_{z}}(t_{i},x^{\ell,m}_{t_{i}};\theta_{z}^{\ell});$

10:
$x^{\ell,m}_{t_{i+1}}=x^{\ell,m}_{t_{i}}+b(t_{i},x^{\ell,m}_{t_{i}},y^{\ell,m}_{t_{i}},z^{\ell,m}_{t_{i}},u^{\ell,m}_{t_{i}})\Delta
t_{i}+\sigma(t_{i},x^{\ell,m}_{t_{i}},y^{\ell,m}_{t_{i}},z^{\ell,m}_{t_{i}},u^{\ell,m}_{t_{i}})\Delta
B_{t_{i}};$

11:
$y^{\ell,m}_{t_{i+1}}=y^{\ell,m}_{t_{i}}-l(t_{i},x^{\ell,m}_{t_{i}},y^{\ell,m}_{t_{i}},z^{\ell,m}_{t_{i}},u^{\ell,m}_{t_{i}})\Delta
t_{i}+z^{\ell,m}_{t_{i}}\Delta B_{t_{i}};$

12:
$loss=loss+f(t_{i},x^{\ell,m}_{t_{i}},y^{\ell,m}_{t_{i}},z^{\ell,m}_{t_{i}},u^{\ell,m}_{t_{i}})\Delta
t_{i};$

13: end for

14:
$J(u(\cdot))=\dfrac{1}{M}\sum\limits_{m=1}^{M}\left[loss+h(x_{t_{N}}^{\ell,m})\right];$

15:
$J(y_{0},z(\cdot))=\dfrac{1}{M}\sum\limits_{m=1}^{M}\left[|y_{t_{N}}^{\ell,m}-g(x_{t_{N}}^{\ell,m})|^{2}\right];$

16: if $\ell\mbox{ mod }(\kappa+1)=0$ then

17: $\theta^{\ell}_{u}=\Xi_{1}(\theta^{\ell}_{u},\nabla J(u(\cdot)));$

18: else

19:
$(\theta^{\ell}_{y},\theta^{\ell}_{z})=\Xi_{2}(\theta^{\ell}_{y},\theta^{\ell}_{z},\nabla
J(y_{0},z(\cdot)));$

20: end if

21:
$(\theta_{y}^{\ell+1},\theta_{z}^{\ell+1},\theta_{u}^{\ell+1})=(\theta_{y}^{\ell},\theta_{z}^{\ell},\theta_{u}^{\ell});$

22:end for

23:$J(u(\cdot))=\dfrac{1}{M}\sum\limits_{m=1}^{M}\left[\sum\limits_{i=0}^{N-1}f(t_{i},x^{\ell,m}_{t_{i}},y^{\ell,m}_{t_{i}},z^{\ell,m}_{t_{i}},u^{\ell,m}_{t_{i}})\Delta
t_{i}+h(x_{t_{N}}^{\ell,m})+\gamma(y_{0}^{\ell,m})\right].$

In Algorithm 1, $maxstep$ represents the total number of parameter update
times. That is to say, we perform $maxstep/(\kappa+1)$ training steps. And in
each training step, there are $\kappa$ times updates for the parameters
$(\theta_{y},\theta_{z})$ and one time update for the parameter $\theta_{u}$.
In this CO method, the choices of $\Xi_{1}$ and $\Xi_{2}$ are independent, and
the training parameters, such as the learning rate, the optimizer, the neural
network architecture, can be chosen separately. And we can choose an
appropriate penalty updating coefficient $\kappa$ in $\Xi_{2}$ to make the
follower’s cost functional $\eqref{eq-dis-obj-z}$ small enough, instead of
simply choosing the approximate penalty coefficient $\mu$ in the penalty
method.

The above proposed algorithm can be regarded as a relaxation method for
solving Problem 2. In Problem 2, the optimization of the leader’s cost
functional $J(u(\cdot))$ should be based on that of the follower’s cost
functional $J(y_{0},z(\cdot))$, thus we should firstly guarantee that
$E\left[|y(T)-g(x(T))|^{2}\right]=0$. In our proposed CO method, we relax this
constraint to improve the computation efficiency. We make the value of
$\left[E|y(T)-g(x(T))|^{2}\right]$ sufficiently close to 0, but not strictly
equal to 0. Specifically, we make the coefficient $\kappa$ to be large enough
so that we can obtain an enough small value of the follower’s cost functional.

##  4 Numerical results

In this section, we show two optimal investment-consumption portfolio examples
solved through recursive utility with our proposed algorithm. If not specially
mentioned, we use a 5-layer fully connected neural network, 512 samples of
Brownian motion in the test set, and the number of time points is $N=25$. The
implementations are performed through TensorFlow on a Lenovo computer with a
2.40 Gigahertz (GHz) Inter Core i7 processor and 8 gigabytes (GB) random-
access memory (RAM).

Here we introduce a continuous and stochastic recursive utility problem which
was first introduced by Duffie and Epstein [DE92]. Suppose there are $n+1$
assets trading continuously in the market, one of which is a risk-free asset
whose price process is given as

$dP^{0}(t)=r(t)P^{0}(t)dt,\qquad P(0)^{0}=p^{0}>0,$ (4.1)

where $r(t)$ is the instantaneous rate of return. The other $n$ assets are
risky assets satisfying

$dP^{i}(t)=P^{i}(t)\left\\{\mu^{i}(t)dt+\left\langle\sigma^{i}(t),dB(t)\right\rangle\right\\},\qquad
P(0)^{i}=p_{0}^{i}>0,$ (4.2)

where $\mu^{i}(t):[0,T]\times\Omega\mapsto\mathbb{R}$ is the appreciation
rate, $\sigma^{i}(t):[0,T]\times\Omega\mapsto\mathbb{R}^{n}$ is the volatility
of the stocks and $B(t)$ is a standard Brownian motion valued in
$\mathbb{R}^{n}$. All the processes are assumed to be
$\mathcal{F}_{t}$-adapted.

An investor starts with a given initial wealth $x_{0}$ and his total wealth at
time $t$ is denoted as $x(t)$. The wealth process $\\{x(t)\\}_{0\leq t\leq T}$
is given by

$\left\\{\begin{array}[]{l}dx(t)=\left\\{x(t)r(t)+\sum\limits_{i=1}^{n}x(t)\pi^{i}(t)(\mu^{i}(t)-r(t))-c(t)\right\\}dt+\sum\limits_{i=1}^{n}x(t)\pi^{i}(t)\left\langle\sigma^{i}(t),dB(t)\right\rangle,\vspace{1ex}\\\
x(0)=x_{0}>0,\end{array}\right.$ (4.3)

where $\pi(t)=(\pi^{1}(t),\pi^{2}(t),\cdots,\pi^{n}(t))$ is called the
portfolio of the investor and $c(t)$ is called the consumption plan process.
$\pi^{i}(t)$ represents the proportion of total wealth $x(t)$ invested in the
$i$-th risky asset and is taking value in $[0,1]$, which means that the short-
selling is prohibited. The remaining fraction
$\pi^{0}(t)=1-\sum_{i=1}^{n}\pi^{i}(t)$ valued in $[0,1]$ is thus the
proportion of the left in the form of risk-free bond. In addition, we assume
that the consumption process is non-negative ($c(t)\geq 0)$.

The utility of the investor at time $t$ is denoted as a function of the
instantaneous consumption $c(t)$ and the future utility. The utility process
can be regarded as the solution of a BSDE given as

$-dy(t)=l(t,c(t),y(t),z(t))dt-\left\langle z(t),dB(t)\right\rangle,\qquad
y(T)=g(x(T)),$ (4.4)

where $g(x)$ is the terminal utility function. Under some regularity
assumptions (such as Assumption 2.1, Assumption 2.2), equation (4.4) has a
unique solution $(y(\cdot),z(\cdot))$. More precisely, the recursive utility
at time $t$ can be denoted by

$y(t)=E\left[\int_{t}^{T}l(s,c(s),y(s),z(s))ds+g(x(T))\Big{|}\mathcal{F}_{t}\right],$
(4.5)

where $\mathcal{F}_{t}$ represents the natural filtration associated with the
standard Brownian motion.

The goal of the investor is to choose the optimal investment and consumption
to maximize the utility at time zero. Mathematically, he aims to find the
optimal controls $\pi^{*}(\cdot)$ and $c^{*}(\cdot)$, such that

$y^{*}(0)=\sup_{(\pi(\cdot),c(\cdot))\in\mathcal{U}_{ad}[0,T]}y_{0},$ (4.6)

where $\mathcal{U}_{ad}[0,T]$ represents the set of all possible investment-
consumption strategies and $y_{0}$ is the recursive utility given as (4.5). We
rewrite $y_{0}$ as

$y_{0}=J(\pi(\cdot),c(\cdot))=E\left[\int_{0}^{T}l(t,c(t),y(t),z(t))dt+g(x(T))\right],$
(4.7)

and we suppose that $l$ is a generalized recursive utility function given as:

$\displaystyle l(t,c,y,z)=-\beta y-w(z)+u(c).$ (4.8)

Here $u(c)$ is the instantaneous utility function. The function $w(z)$ is a
differentiable function and can be regarded as a risk-aversion coefficient.

Actually, the recursive utility problem (4.6) is a special case of Problem 1,
in which $f(t,x,y,z,u)=0$, $h(x)=0$ and $\gamma(y)=y$, therefore we can
reformulate it to Problem 2, where the goal of the leader is to find the
optimal controls $(\pi^{*}(t),c^{*}(t))_{0\leq t\leq T}$ which maximize (4.7),
and the goal of the follower is to minimize

$J(y_{0},z(\cdot))=E\left[|y(T)-g(x(T))|^{2}\right].$ (4.9)

under the given controls $(\pi^{*}(t),c^{*}(t))$. Then we solve the problem
through our proposed CO algorithm presented in section 3.

In this example, in order to approximate the controls $\pi(t)$ and $c(t)$, we
need to construct two neural networks $\mathcal{NN}_{\pi}$, $\mathcal{NN}_{c}$
, and both of the controls are supposed as the feedback controls of the time
and the wealth process $(t,x(t))$. The neural network approximating $\pi(t)$
contains one $(n+1)$-dimensional input layer, three $100$-dimensional hidden
layers and one $(n+1)$-dimensional output layer. For the constraint of the
portfolio, $\pi^{i}(t)\in[0,1]$ and $\sum_{i=0}^{n}\pi^{i}(t)=1$ for
$i=0,\cdots,n$, we deal with it through a softmax function $\psi$. The
function $\psi$ is given as

$\psi(x)=\left(\dfrac{e^{x^{0}}}{\sum_{i=0}^{n}{e^{x^{i}}}},\dfrac{e^{x^{1}}}{\sum_{i=0}^{n}{e^{x^{i}}}},\cdots,\dfrac{e^{x^{n}}}{\sum_{i=0}^{n}{e^{x^{i}}}}\right),\qquad
x=(x^{0},\cdots,x^{n})\in\mathbb{R}^{n+1},$

which is often used in classification problems such as image classification.
Intuitively, it makes sense to choose the softmax function because choosing
the best investment asset based on the current state is essentially a
classification problem. Moreover, we use a function $\psi(x)=\max(x,0)$ to
deal with the constraint $c(t)\geq 0$. In addition, we construct the other two
neural networks to approximate $y_{0}$ and $z(t)$, respectively. In all, four
feed-forward neural networks are constructed in this example, the first for
approximating the consumption rate $c(t)$, the second for approximating the
wealth proportion $\pi(t)$, the third for approximating $y_{0}$, and the last
for approximating $z(t)$.

Note that in (4.9), $y_{0}$ is regarded as a control parameter which is
approximated through deep neural network. While at the same time, it is the
optimization objective in (4.7). In order to make difference, we call the
$y_{0}$ in (4.7) the integral form and that in the neural network the
parametric form. And we should emphasize that instead of optimizing the
parametric form $y_{0}$, we optimize the integral form
$y_{0}=J(\pi(\cdot),c(\cdot))$ with respect to $(\pi(\cdot),c(\cdot))$. The
reason is that the derivatives of the parametric form $y_{0}$ with respect to
$(\pi(\cdot),c(\cdot))$ do not exist. What’s more, we compare the integral
form with the parametric form of $y_{0}$. We measure the distance of the
values between the parametric form and the integral form and take it as a
criterion for the effectiveness of our proposed algorithm.

###  4.1 Case 1: a linear driver

In the first example, we consider a linear driver and set $w(z)\equiv 0$, then
the recursive utility functional we need to maximize is given as

$y_{0}=J(\pi(\cdot),c(\cdot))=E\left[\int_{0}^{T}(-\beta
y(t)+u(c(t)))dt+g(x(T))\right].$ (4.10)

In this case, (4.10) is equivalent to the stochastic functional

$y_{0}=J^{\prime}(\pi(\cdot),c(\cdot))=E\left[\int_{0}^{T}e^{-\beta
t}u(c(t))dt+e^{-\beta T}g(x(T))\right].$ (4.11)

Thus the stochastic recursive utility problem (4.3) and (4.7) can degenerate
to a classic stochastic optimal control problem as (4.3) and (4.11). There
have been some deep learning methods for solving the stochastic control
problem (4.3)-(4.11), such as the methods in [HE16, BHLP18, JPPZ21].

Here we suppose that $u(c)$ is a quadratic utility function

$u(c)=(c-c^{2}),$

and the terminal utility function $g(x)$ is

$g(x)=e^{-x}.$

We calculate the recursive utility problem with the previously mentioned
investment-consumption model. Let the wealth process and trading assets
processes satisfy (4.1)-(4.3). We set $r(t)=0.03$, $\mu^{i}(t)=0.05$,
$\sigma=(\sigma^{i},\cdots,\sigma^{n})=0.1I_{n}$ for $i=1,\cdots,n$, and the
initial wealth $x_{0}$ is equal to $100$. The parameters $\beta$ in
$l(t,c,y,z)$ is set to be $0.05$. As a comparison, we also calculate the
classic stochastic optimal control problem (4.3) and (4.11). Here we construct
a similar neural network with that in [HE16] and regard the discrete time
point $t$ as an input of the neural network.

The implementation results with different dimensions and different terminal
time are shown in Table 1. All the results are take from 5 independent runs.
The number of iterations for parameter update is 18,000 and $\kappa$ is set to
be 19. That is to say, after 19 iterations of optimization for the follower’s
cost functional (4.9), we perform one iteration of optimization for the
leader’s cost functional (4.7). And the total number of training steps is
$18,000/(19+1)=900$. We show the integral form and the parametric form of
$y_{0}$ (Inte. $y_{0}$ and Para. $y_{0}$ in the table, respectively), and
calculate their mean values and variances among 5 independent runs. We also
show the values calculated through the classic stochastic control problem
(4.3) and (4.11) (Clas. $y_{0}$ in the table).

Table 1: The implementation results of $y_{0}$ for case 1

| $T=0.25$ | $T=0.50$ | $T=0.75$ | $T=1.00$
---|---|---|---|---
| Mean | Var. | Mean | Var. | Mean | Var. | Mean | Var.
$n=10$ | Clas. $y_{0}$ | 0.06192 | 8.0998e-08 | 0.12323 | 7.2303e-08 | 0.18416 | 1.7761e-10 | 0.24314 | 3.9102e-07
Inte. $y_{0}$ | 0.06161 | 2.6347e-08 | 0.12269 | 9.8344e-08 | 0.18299 | 2.9634e-07 | 0.24216 | 7.3878e-07
Para. $y_{0}$ | 0.06146 | 3.1114e-07 | 0.12194 | 3.0065e-07 | 0.18217 | 3.3235e-06 | 0.24137 | 1.4374e-06
$n=20$ | Clas. $y_{0}$ | 0.06183 | 1.2111e-07 | 0.12301 | 3.9042e-07 | 0.18389 | 1.1229e-07 | 0.24359 | 3.3391e-07
Inte. $y_{0}$ | 0.06169 | 2.3650e-08 | 0.12268 | 2.9813e-07 | 0.18248 | 3.5300e-07 | 0.24202 | 4.7493e-07
Para. $y_{0}$ | 0.06169 | 2.9493e-07 | 0.12238 | 8.1997e-07 | 0.18149 | 1.0537e-06 | 0.24084 | 5.2556e-06
$n=30$ | Clas. $y_{0}$ | 0.06160 | 3.7891e-07 | 0.12282 | 6.4121e-07 | 0.18376 | 7.0530e-08 | 0.24315 | 8.0416e-07
Inte. $y_{0}$ | 0.06137 | 4.5255e-08 | 0.12203 | 1.9316e-07 | 0.18204 | 1.0402e-07 | 0.24268 | 6.7658e-07
Para. $y_{0}$ | 0.06036 | 3.8921e-07 | 0.12169 | 5.4307e-06 | 0.17921 | 4.2815e-07 | 0.24098 | 5.2413e-06
$n=40$ | Clas. $y_{0}$ | 0.06171 | 1.1584e-07 | 0.11055 | 1.3581e-03 | 0.18333 | 4.5646e-07 | 0.24321 | 8.7775e-07
Inte. $y_{0}$ | 0.06142 | 2.0930e-08 | 0.12226 | 2.3594e-07 | 0.18235 | 1.6639e-07 | 0.24228 | 5.0222e-07
Para. $y_{0}$ | 0.06129 | 4.6174e-07 | 0.12092 | 2.7358e-06 | 0.18113 | 3.5940e-06 | 0.23994 | 4.3413e-06
$n=50$ | Clas. $y_{0}$ | 0.06170 | 9.4358e-08 | 0.12303 | 2.7746e-07 | 0.18324 | 1.2237e-07 | 0.24126 | 1.9552e-06
Inte. $y_{0}$ | 0.06143 | 2.0046e-08 | 0.12222 | 6.9066e-08 | 0.18273 | 5.7814e-08 | 0.24134 | 3.6883e-07
Para. $y_{0}$ | 0.06077 | 2.7034e-06 | 0.12160 | 3.2163e-06 | 0.18115 | 1.3218e-06 | 0.23722 | 1.3412e-06

The results in Table 1 show that our proposed method is effective for the
recursive utility problem. The values of the utility $y_{0}$ between the
integral form and the parametric form are close, and both of the two forms are
close to the value of the classic form. We also exhibit the curves of the
utility values of different forms in Figure 3. The left figure shows the mean
and scope for the different forms of $y_{0}$ among 5 independent runs. We can
see that the values of the three different forms of $y_{0}$ are getting closer
with the increase of the number of iterations. At the same time, the value of
the integral form is getting larger which meets our objective to maximize the
recursive utility functional (4.7) (see the green curve and scope). The right
figure shows the mean and scope of the distance between the integral form and
the parametric form.

Figure 3: Case $n=50$ and $T=1.0$

In order to show the impact of the penalty updating coefficient $\kappa$ to
the convergence of the algorithm, we vary the value of $\kappa$ for $T=0.5$
and $n=10$. The results with different $\kappa$ are shown in Table 2 and
Figure 4. In Figure 4, the left figure shows that the distances between the
integral form and the parametric form are small when $\kappa>1$, which means
that the value of (4.9) can be small enough when we perform more iterations on
(4.9). The right figure shows the values of the utility functional are more
stable for the case $\kappa>1$ than the case $\kappa<1$. Here $\kappa<1$ means
that we perform more iterations of optimization on (4.7) than that on (4.9).

Figure 4: The distance of $y_{0}$ between the integral form and the parametric form and the values of $y_{0}$ with different $\kappa$ for $n=10$ and $T=0.5$. Table 2: The implementation results of $y_{0}$ with different $\kappa$ for case 1 | Inte. $y_{0}$ | Para. $y_{0}$ | Distance
---|---|---|---
| Mean | Var. | Mean | Var. | Mean | Var.
$\kappa=49$ | 0.12184 | 1.9438e-07 | 0.12056 | 1.9288e-06 | 0.00145 | 5.2404e-07
$\kappa=29$ | 0.12194 | 1.4726e-07 | 0.12045 | 2.4452e-06 | 0.00182 | 4.1077e-07
$\kappa=19$ | 0.12243 | 7.4394e-08 | 0.12182 | 6.7248e-07 | 0.00107 | 2.2293e-07
$\kappa=9$ | 0.12295 | 4.0809e-07 | 0.12347 | 1.2084e-08 | 0.00055 | 4.4437e-07
$\kappa=4$ | 0.12332 | 6.8988e-09 | 0.12326 | 2.2820e-09 | 0.00010 | 5.4356e-09
$\kappa=1$ | 0.12323 | 9.5936e-08 | 0.12218 | 5.6867e-08 | 0.00106 | 2.7134e-07
$\kappa=1/4$ | 0.12340 | 1.9741e-08 | 0.11848 | 1.3041e-07 | 0.00492 | 1.3164e-07
$\kappa=1/9$ | 0.12354 | 3.0263e-08 | 0.11218 | 6.6298e-07 | 0.01136 | 7.8541e-07
$\kappa=1/19$ | 0.12386 | 2.7079e-08 | 0.10002 | 3.1441e-07 | 0.02384 | 2.6049e-07
$\kappa=1/49$ | 0.12461 | 2.8228e-08 | 0.06496 | 9.7174e-08 | 0.05965 | 1.1475e-07

###  4.2 Case 2: a nonlinear driver

In this case, the function $w(z)$ is given as $w(z)=\frac{\nu}{2}|z|^{2}$,
$|z|^{2}=\left\langle z,z\right\rangle$, for $z\in\mathbb{R}^{n}$, where $\nu$
is a given constant $\nu=10$. The functions $u(c)$, $g(x)$ and the other
settings are the same as the linear driver.

The implementation results with different dimensions and different terminal
times for this example are shown in Table 3. We set $\kappa=19$ and the total
number of iterations to be 18,000. The computation results are similar to that
with the linear driver. The distances between the integral form and the
parametric form of $y_{0}$ are small at the end of the training. In Figure 5,
we show the curves for the values and the distances of $y_{0}$ between the
integral form and the parametric form with $n=50$ and $T=0.5$. We can see that
the distances between the integral form and the parametric form of $y_{0}$ are
getting closer to 0 when the number of iterations increases. Besides, the
value of the integral form is getting larger, which shows that we can find the
optimal investment-consumption strategy $(\pi^{*}(\cdot),c^{*}(\cdot))$ to
maximize the recursive utility functional (4.7).

Table 3: The implementation results of $y_{0}$ for case 2

| $T=0.25$ | $T=0.50$ | $T=0.75$ | $T=1.00$
---|---|---|---|---
| Mean | Var. | Mean | Var. | Mean | Var. | Mean | Var.
$n=10$ | Inte. $y_{0}$ | 0.06166 | 2.2275e-08 | 0.12293 | 7.4394e-08 | 0.18284 | 4.2424e-07 | 0.24182 | 1.0601e-06
Para. $y_{0}$ | 0.06110 | 4.1587e-07 | 0.12240 | 6.7248e-07 | 0.18236 | 6.0408e-07 | 0.23978 | 8.4018e-06
Distance | 0.00074 | 4.6313e-08 | 0.00064 | 2.2293e-07 | 0.00063 | 1.3188e-07 | 0.00236 | 2.5190e-06
$n=20$ | Inte. $y_{0}$ | 0.06156 | 1.3603e-07 | 0.12266 | 1.6195e-07 | 0.18284 | 9.5716e-08 | 0.24167 | 2.7468e-06
Para. $y_{0}$ | 0.06123 | 1.2727e-06 | 0.12302 | 3.9956e-07 | 0.18352 | 8.2581e-07 | 0.24052 | 1.1217e-05
Distance | 0.00066 | 3.2714e-07 | 0.00054 | 7.1423e-08 | 0.00117 | 9.8714e-08 | 0.00164 | 1.7542e-06
$n=30$ | Inte. $y_{0}$ | 0.06157 | 4.3009e-08 | 0.12270 | 8.5964e-08 | 0.18256 | 8.1253e-07 | 0.24257 | 7.8675e-07
Para. $y_{0}$ | 0.06126 | 4.5989e-07 | 0.12206 | 1.2343e-06 | 0.18123 | 4.9427e-06 | 0.24145 | 2.5426e-06
Distance | 0.00052 | 5.2667e-08 | 0.00109 | 3.3751e-07 | 0.00155 | 1.1978e-06 | 0.00112 | 5.7840e-07
$n=40$ | Inte. $y_{0}$ | 0.06145 | 2.8057e-08 | 0.12221 | 3.3040e-07 | 0.18304 | 6.4776e-07 | 0.24183 | 3.4355e-07
Para. $y_{0}$ | 0.06089 | 5.5206e-07 | 0.12065 | 4.3699e-06 | 0.18207 | 3.0924e-06 | 0.23964 | 1.2486e-06
Distance | 0.00071 | 1.4810e-07 | 0.00183 | 1.5946e-06 | 0.00097 | 9.1203e-07 | 0.00218 | 4.1711e-07
$n=50$ | Inte. $y_{0}$ | 0.06149 | 4.8995e-08 | 0.12251 | 1.1746e-07 | 0.18300 | 4.2698e-07 | 0.24273 | 9.6758e-08
Para. $y_{0}$ | 0.06109 | 7.2222e-07 | 0.12152 | 6.1477e-07 | 0.18203 | 1.6656e-06 | 0.24172 | 3.8806e-07
Distance | 0.00069 | 1.3734e-07 | 0.00099 | 2.3589e-07 | 0.00097 | 4.8985e-07 | 0.00101 | 1.3610e-07

Figure 5: Case $n=50$ and $T=0.5$. The left figure shows the mean and scope
for the values of $y_{0}$ among 5 independent runs. We can see that the two
different forms of $y_{0}$ are getting closer when the number of iteration
steps increases. Moreover, the value of the integral form of $y_{0}$ (see the
green curve and scope) is getting larger and then tends to be stable. The
right figure shows the mean and scope of the distances between the different
forms of $y_{0}$.

We also vary the value of $\kappa$ for $n=10$ and $T=0.5$. The computation
results are shown in Table 4 and Figure 6, which also exhibit that the
distances between the two different forms of $y_{0}$ are smaller and more
stable when $\kappa>1$.

Figure 6: The distances of $y_{0}$ between the two different forms (the integral form and the parametric form) and the values of $y_{0}$ with different $\kappa$ for $n=10$ and $T=0.50$. Table 4: The implementation results of $y_{0}$ with different $\kappa$ for case 2 | Inte. $y_{0}$ | Para. $y_{0}$ | Distance
---|---|---|---
| Mean | Var. | Mean | Var. | Mean | Var.
$\kappa=49$ | 0.12150 | 3.8599e-07 | 0.11983 | 2.0846e-06 | 0.00167 | 9.6376e-07
$\kappa=29$ | 0.12150 | 8.3208e-08 | 0.11965 | 7.7199e-07 | 0.00185 | 3.5862e-07
$\kappa=19$ | 0.12293 | 7.4394e-08 | 0.12240 | 6.7248e-07 | 0.00064 | 2.2293e-07
$\kappa=9$ | 0.12307 | 1.8615e-07 | 0.12350 | 6.5421e-10 | 0.00043 | 2.0767e-07
$\kappa=4$ | 0.12289 | 1.8920e-07 | 0.12322 | 1.3545e-09 | 0.00040 | 1.5793e-07
$\kappa=1$ | 0.12255 | 1.7449e-06 | 0.12219 | 1.3545e-09 | 0.00124 | 1.5793e-07
$\kappa=1/4$ | 0.12324 | 5.4207e-08 | 0.11862 | 2.6561e-08 | 0.00462 | 6.7300e-08
$\kappa=1/9$ | 0.12362 | 1.2905e-09 | 0.11254 | 2.7437e-08 | 0.01108 | 3.5992e-08
$\kappa=1/19$ | 0.12357 | 2.6459e-07 | 0.10039 | 9.0642e-09 | 0.02318 | 2.2544e-07
$\kappa=1/49$ | 0.12461 | 1.0065e-07 | 0.06392 | 7.9881e-10 | 0.06069 | 1.1117e-07

The above results demonstrate that our proposed method is effective for
solving the stochastic recursive utility problem. Besides, we should do more
optimizations on the follower’s cost functionals, as verified by the
implementation results.

##  5 Conclusion

In this paper, we propose a deep learning method for solving the stochastic
control problem driven by fully-coupled FBSDEs. We transform it into a
stochastic Stackelberg differential game and propose a CO method to solve this
new problem. The CO method has high flexibility and we can set the training
parameters of the neural networks seperately to optimize the leader’s cost
functional and the follower’s cost functional independently. In order to show
the performance of our algorithm, we give two optimal investment-consumption
portfolio examples solved through the stochastic recursive utility models. The
numerical results demonstrate remarkable performance.

## References

  * [Ant93] Fabio Antonelli.  Backward-forward stochastic differential equations.  Annals of Applied Probability, 3(3):777–793, 1993.
  * [BEJ19] Christian Beck, Weinan E, and Arnulf 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(4):1563–1619, 2019.
  * [BHLP18] Achref Bachouch, Come Hure, Nicolas Langrene, and Huyen Pham.  Deep neural networks algorithms for stochastic control problems on finite horizon, part 2: numerical applications.  arXiv: 1812.05916v1, 2018.
  * [Bis73] Jeanmichel Bismut.  Conjugate convex functions in optimal stochastic control.  Journal of Mathematical Analysis and Applications, 44(2):384–404, 1973.
  * [BT04] Bruno Bouchard and Nizar Touzi.  Discrete-time approximation and monte-carlo simulation of backward stochastic differential equations.  Stochastic Processes & Their Applications, 111(2):175–206, 2004\.
  * [BZ08] Christian Bender and Jianfeng Zhang.  Time discretization and markovian iteration for coupled fbsdes.  Annals of Applied Probability, 18(1):143–177, 2008.
  * [CE02] Zengjing Chen and Larry Epstein.  Ambiguity, risk, and asset returns in continuous time.  Econometrica, 70(4):1403–1443, 2002.
  * [CL] René Carmona and Mathieu 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.
  * [CM96] Jakša Cvitanić and Jin Ma.  Hedging options for a large investor and forward-backward sde’s.  The annals of applied probability, 6(2):370–398, 1996.
  * [DE92] Darrell Duffie and Larry G Epstein.  Stochastic differential utility.  Econometrica: Journal of the Econometric Society, pages 353–394, 1992.
  * [DZ99] Nikolai Dokuchaev and Xun Yu Zhou.  Stochastic controls with terminal contingent conditions.  Journal of Mathematical Analysis and Applications, 238(1):143–165, 1999.
  * [EHJ17] Weinan E, Jiequn Han, and Arnulf Jentzen.  Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations.  Communications in Mathematics & Statistics, 5(4):349–380, 2017\.
  * [EKPQ97] Nicole El Karoui, Shige Peng, and Marie Claire Quenez.  Backward stochastic differential equations in finance.  Mathematical finance, 7(1):1–71, 1997.
  * [EKPQ01] Nicole El Karoui, Shige Peng, and Marie Claire Quenez.  A dynamic maximum principle for the optimization of recursive utilities under constraints.  Annals of applied probability, pages 664–693, 2001.
  * [FZT16] Yu Fu, Weidong Zhao, and Zhou Tao.  Efficient spectral sparse grid approximations for solving multi-dimensional forward backward sdes.  Discrete and Continuous Dynamical Systems - Series B, 22(9), 2016\.
  * [FZZ16] Yu Fu, Weidong Zhao, and Tao Zhou.  Multistep schemes for forward backward stochastic differential equations with jumps.  Journal of Scientific Computing, 69(2):1–22, 2016.
  * [HE16] Jiequn Han and Weinan E.  Deep learning approximation for stochastic control problems.  NIPS Workshop on Deep Reinforcement Learning, 2016.
  * [HJE18] Jiequn Han, Arnulf Jentzen, and Weinan E.  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences of the United States of America, 115(34):8505–8510, 2018.
  * [HJX18] Mingshang Hu, Shaolin Ji, and Xiaole Xue.  A global stochastic maximum principle for fully coupled forward-backward stochastic systems.  SIAM Journal on Control and Optimization, 56(6):4309–4335, 2018\.
  * [HJX19] Mingshang Hu, Shaolin Ji, and Xiaole Xue.  The existence and uniqueness of viscosity solution to a kind of hamilton–jacobi–bellman equation.  SIAM Journal on Control and Optimization, 57(6):3911–3938, 2019\.
  * [HJX20] Hu, Mingshang, Ji, Shaolin, and Xue, Xiaole.  Stochastic maximum principle, dynamic programming principle, and their relationship for fully coupled forward-backward stochastic controlled systems.  ESAIM: COCV, 26:81, 2020.
  * [HL18] Jiequn Han and Jihao Long.  Convergence of the deep bsde method for coupled fbsdes.  arXiv:1811.01165, 2018.
  * [HP95] Ying Hu and Shige Peng.  Solution of forward-backward stochastic differential equations.  Probability Theory and Related Fields, 103(2):273–283, 1995.
  * [HPW20] Côme Huré, Huyên Pham, and Xavier Warin.  Deep backward schemes for high-dimensional nonlinear PDEs.  arXiv:1902.01599 [cs, math, stat], June 2020.
  * [HRO16] Thomas P Huijskens, Marjon Ruijter, and Cornelis W Oosterlee.  Efficient numerical fourier methods for coupled forward-backward sdes.  Journal of Computational and Applied Mathematics, 296:593–612, 2016\.
  * [JJZ08] M. A. Jin, Shen Jie, and Yanhong Zhao.  On numerical approximations of forward-backward stochastic differential equations.  Siam Journal on Numerical Analysis, 46(5):2636–2661, 2008.
  * [JPPZ20] Shaolin Ji, Shige Peng, Ying Peng, and Xichuan Zhang.  Three algorithms for solving high-dimensional fully coupled fbsdes through deep learning.  IEEE Intelligent Systems, 35(3):71–84, 2020.
  * [JPPZ21] Shaolin Ji, Shige Peng, Ying Peng, and Xichuan Zhang.  Solving stochastic optimal control problem via stochastic maximum principle with deep learning method.  2021\.
  * [JZ06] Shaolin Ji and Xun Yu Zhou.  A maximum principle for stochastic optimal control with terminal state constraints, and its applications.  Communications in Information and Systems, 6(4):321 – 338, 2006\.
  * [KZ00] Michael Kohlmann and Xunyu Zhou.  Relationship between backward stochastic differential equations and stochastic controls: a linear-quadratic approach.  SIAM Journal on Control and Optimization, 38(5):1392–1392, 2000\.
  * [LW14] Juan Li and Qingmeng Wei.  Optimal control problems of fully coupled fbsdes and viscosity solutions of hamilton–jacobi–bellman equations.  SIAM Journal on Control and Optimization, 52(3):1622–1662, 2014\.
  * [LZ01] Andrew E. B. Lim and Xunyu Zhou.  Linear-quadratic control of backward stochastic differential equations.  SIAM Journal on Control and Optimization, 2001.
  * [Min17] Hu Mingshang.  Stochastic global maximum principle for optimization with recursive utilities.  Probability, Uncertainty and Quantitative Risk, 2(1), 2017.
  * [MPY94] Jin Ma, Philip Protter, and Jiongmin Yong.  Solving forward-backward stochastic differential equations explicitly—a four step scheme.  Probability theory and related fields, 98(3):339–359, 1994.
  * [MT04] G. N. Milstein and M. V. Tretyakov.  Numerical algorithms for forward-backward stochastic differential equations connected with semilinear parabolic equations.  volume 28, pages 561–582, 2004.
  * [Pen90] Shige Peng.  A general stochastic maximum principle for optimal control problems.  Siam Journal on Control and Optimization, 28(4):966–979, 1990.
  * [Pen93] Shige Peng.  Backward stochastic differential equations and applications to optimal control.  Applied Mathematics and Optimization, 27(2):125–144, 1993.
  * [PP90] Etienne Pardoux and Shige Peng.  Adapted solution of a backward stochastic differential equation.  Systems & Control Letters, 14(1):55–61, 1990.
  * [PT99] Etienne Pardoux and Shanjian Tang.  Forward-backward stochastic differential equations and quasilinear parabolic pdes.  Probability Theory and Related Fields, 114(2):123–150, 1999.
  * [PW99] Shige Peng and Zhen Wu.  Fully coupled forward-backward stochastic differential equations and applications to optimal control.  Siam Journal on Control and Optimization, 37(3):825–843, 1999.
  * [PWG21] Huyên Pham, Xavier Warin, and Maximilien Germain.  Neural networks-based backward scheme for fully nonlinear PDEs.  SN Partial Differential Equations and Applications, 2:16, 2021.
  * [Tad12] Eitan Tadmor.  A review of numerical methods for nonlinear partial differential equations.  Bulletin of the American Mathematical Society, 49(4):507–554, 2012\.
  * [Wu98] Zhen Wu.  Maximum principle for optimal control problem of fully coupled forward-backward stochastic systems.  Systems Science and Mathematical Sciences, 11(3):249–259, 1998\.
  * [Wu13] Zhen Wu.  A general maximum principle for optimal control of forward backward stochastic systems.  Automatica, 49(5):1473–1480, 2013.
  * [Yon10] Jiongmin Yong.  Optimality variational principle for controlled forward-backward stochastic differential equations with mixed initial-terminal conditions.  SIAM Journal on Control and Optimization, 48(6):4119–4156, 2010\.
  * [YZ99] Jiongmin Yong and Xunyu Zhou.  Stochastic Controls-Hamiltonian System and HJB Equations.  Springer, 1999.