Deep Learning Methods for Mean Field Control Problems With Delay

We consider a general class of mean field control problems described by stochastic delayed differential equations of McKean–Vlasov type. Two numerical algorithms are provided based on deep learning techniques, one is to directly parameterize the optimal control using neural networks, the other is based on numerically solving the McKean–Vlasov forward anticipated backward stochastic differential equation (MV-FABSDE) system. In addition, we establish the necessary and sufficient stochastic maximum principle of this class of mean field control problems with delay based on the differential calculus on function of measures, and the existence and uniqueness results are proved for the associated MV-FABSDE system under suitable conditions.

Mathematical Subject Classification (2000): 93E20, 60G99, 68-04

1. Introduction

Stochastic games were introduced to study the optimal behaviors of agents interacting with each other. They are used to study the topic of systemic risk in the context of finance. For example, in Carmona et al. [1], the authors proposed a linear quadratic inter-bank borrowing and lending model, and solved explicitly for the Nash equilibrium with a finite number of players. Later, this model was extended in Carmona et al. [2] by considering delay in the control in the state dynamic to account for the debt repayment. The authors analyzed the problem via a probabilistic approach which relies on stochastic maximum principle, as well as via an analytic approach which is built on top of an infinite dimensional dynamic programming principle.

Both mean field control and mean field games are used to characterize the asymptotic behavior of a stochastic game as the number of players grows to infinity under the assumption that all the agents behave similarly, but with different notion of equilibrium. The mean field games consist of solving a standard control problem, where the flow of measures is fixed, and solving a fixed point problem such that this flow of measures matches the distribution of the dynamic of a representative agent. Whereas, the mean field control problem is a non-standard control problem in the sense that the law of state is present in the McKean–Vlasov dynamic, and optimization is performed while imposing the constraint of distribution of the state. More details can be found in Carmona and Delarue [3] and Bensoussan et al. [4].

In this paper, we considered a general class of mean field control problem with delay effect in the McKean–Vlasov dynamic. We derived the adjoint process associated with the McKean–Vlasov stochastic delayed differential equation, which is an anticipated backward stochastic differential equation of McKean–Vlasov type due to the fact that the conditional expectation of the future of adjoint process as well as the distribution of the state dynamic are involved. This type of anticipated backward stochastic differential equations (BSDE) was introduced in Peng and Yang [5], and for the general theory of BSDE, we refer Zhang [6]. The necessary and sufficient part of stochastic maximum principle for control problem with delay in state and control can be found in Chen and Wu [7]. Here, we also establish a necessary and sufficient stochastic maximum principle based on differential calculus on functions of measures as we consider the delay in the distribution. In the meantime, we also prove the existence and uniqueness of the system of McKean–Vlasov forward anticipated backward stochastic differential equations (MV-FABSDE) under some suitable conditions using the method of continuation, which can be found in Zhang [6], Peng and Wu [8], Bensoussan et al. [9], and Carmona et al. [2]. For a comprehensive study of FBSDE theory, we refer to Ma and Yong [10].

When there was no delay effect in the dynamic, Ma et al. [11] proved the relation between the solution to the FBSDE and quasi-linear partial differential equation (PDE) via “Four Step Scheme.” E et al. [12] and Raissi [13] explored the use of deep learning for solving these PDEs in high dimensions. However, the class of fully coupled MV-FABSDE considered in our paper has no explicit solution. Here, we present one algorithm to tackle the above problem by means of deep learning techniques. Due to the non-Markovian nature of the state dynamic, we apply the long short-term memory (LSTM) network, which is able to capture the arbitrary long-term dependencies in the data sequence. It also partially solves the vanishing gradient problem in vanilla recurrent neural networks (RNNs), as was shown in Hochreiter and Schmidhuber [14]. The idea of our algorithm is to approximate the solution to the adjoint process and the conditional expectation of the adjoint process. The optimal control is readily obtained after the MV-FABSDE being solved. We may also emphasize that the way that we present here for numerically computing conditional expectation may have a wide range of applications, and it is simple to implement. We also present another algorithm solving the mean field control problem by directly parameterizing the optimal control. Similar idea can be found in the policy gradient method in the regime of reinforcement learning [15] as well as in Han and E [16]. Numerically, the two algorithms that we propose in this paper yield the same results. Besides, our approaches are benchmarked to the case with no delay for which we have explicit solutions.

The paper is organized as follows. We start with an N-player game with delay, and let number of players goes to infinity to introduce a mean field control problem in section 2. Next, in section 3, we mathematically formulate the feedforward neural networks and LSTM networks, and we propose two algorithms to numerically solve the mean field control problem with delay using deep learning techniques. This is illustrated on a simple linear-quadratic toy model, however with delay in the control. One algorithm is based on directly parameterizing the control, and the other depends on numerically solving the MV-FABSDE system. In addition, we also provide an example of solving a linear quadratic mean field control problem with no delay both analytically, and numerically. The adjoint process associated with the delayed dynamic is derived, as well as the stochastic maximum principle is proved in section 4. Finally, the uniqueness and existence solution for this class of MV-FABSDE are proved under suitable assumptions via continuation method in section 5.

2. Formulation of the Problem

We consider an N-player game with delay in both state and control. The dynamic ${(X_t^i)}_{0\le t\le T}$ for player i ∈ {1, …, N} is given by a stochastic delayed differential equation (SDDE),

$\begin{array}{l}d{X}_t^i=b^i(t,{\text{X}}_t,{\text{X}}_{t-\tau },{\alpha }_t^i,{\alpha }_{t-\tau }^i)dt+{\sigma }^i(t,{\text{X}}_t,{\text{X}}_{t-\tau },{\alpha }_t^i,{\alpha }_{t-\tau }^i)d{W}_t^i,\\ t\in (0,T],\\ X_0^i=x_0^i,\\ X_t^i={\alpha }_t^i=0;t\in [-\tau ,0),& (2.1)\end{array}$

for T > 0, τ > 0 given constants, where ${\text{X}}_t=(X_t^1,\cdots ,X_t^N)$, and where ${({(W_t^i)}_{t\in \left[0,T\right]})}_{i=1,\cdots ,N}$ are N independent Brownian motions defined on the space $(\Omega ,{F},\mathbb{P})$, ${({{F}}_t)}_{0\le t\le T}$ being the natural filtration of Brownian motions.

$(b^i,{\sigma }^i):\left[0,T\right]\times \times {ℝ}^N\times {ℝ}^N\times A\times A\to ℝ\times ℝ,$

are progressively measurable functions with values in ℝ. We denote A a closed convex subset of ℝ, the set of actions that player i can take, and denote 𝔸 the set of admissible control processes. For each i ∈ {1, …, N}, A-valued measurable processes ${({\alpha }_t^i)}_{0\le t\le T}$ satisfy an integrability condition such that $𝔼\left[{\int }_{-\tau }^T|{\alpha }_t^i{|}^{2}dt\right]<+\infty$.

Given an initial condition ${\text{x}}_0=(x_0^1,\cdots ,x_0^N)\in {ℝ}^N$, each player would like to minimize his objective functional:

$\begin{array}{ll}{J}^i(\alpha )=𝔼\left[{\displaystyle {\int }_0^T}{f}^i(t,{\text{X}}_t,{\text{X}}_{t-\tau },{\alpha }_t^i)dt+g^i({\text{X}}_T)\right],& (2.2)\end{array}$

for some Borel measurable functions fⁱ:[0, T] × ℝ^N × ℝ^N × A → ℝ, and gⁱ:ℝ^N → ℝ.

In order to study the mean-field limit of (X_t)_t∈[0,T], we assume that the system (2.1) satisfy a symmetric property, that is to say, for each player i, the other players are indistinguishable. Therefore, drift bⁱ and volatility σⁱ in (2.1) take the form of

$\begin{array}{l}(b^i,{\sigma }^i)(t,{\text{X}}_t,{\text{X}}_{t-\tau },{\alpha }_t^i,{\alpha }_{t-\tau }^i)=(b^i,{\sigma }^i)(t,X_t^i,{\mu }_t^N,X_{t-\tau }^i,{\mu }_{t-\tau }^N,{\alpha }_t^i,{\alpha }_{t-\tau }^i),\end{array}$

and the running cost fⁱ and terminal cost gⁱ are of the form

$\begin{array}{l}{f}^i(t,{\text{X}}_t,{\text{X}}_{t-\tau },{\alpha }_t^i)=f^i(t,X_t^i,{\mu }_t^N,X_{t-\tau }^i,{\mu }_{t-\tau }^N,{\alpha }_t^i)\text{ and }{g}^i({\text{X}}_T)\\ =g^i(X_T^i,{\mu }_T^N),\end{array}$

where we use the notation ${\mu }_t^N$ for the empirical distribution of X = (X¹, ⋯ , X^N) at time t, which is defined as

$\begin{array}{l}{\mu }_t^N=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{\delta }_{X_t^j}.\end{array}$

Next, we let the number of players N go to +∞ before we perform the optimization. According to symmetry property and the theory of propagation of chaos, the joint distribution of the N dimensional process ${({\text{X}}_t)}_{0\le t\le T}={(X_t^1,\dots ,X_t^N)}_{0\le t\le T}$ converges to a product distribution, and the distribution of each single marginal process converges to the distribution of (X_t)_0≤t≤T of the following McKean–Vlasov stochastic delayed differential equation (MV-SDDE). For more detail on the argument without delay, we refer to Carmona and Delarue [3] and Carmona et al. [17].

$\begin{array}{l}d{X}_t=b(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t,{\alpha }_{t-\tau })dt\\ +\sigma (t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t,{\alpha }_{t-\tau })d{W}_t,t\in (0,T],\\ X_0=x_0,\\ X_t={\alpha }_t=0,t\in [-\tau ,0).& (2.3)\end{array}$

We then optimize after taking the limit. The objective for each player of (2.2) now becomes

$\begin{array}{ll}J(\alpha )=𝔼\left[{\displaystyle {\int }_0^T}f(X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t)dt+g(X_T,{\mu }_T)\right],& (2.4)\end{array}$

where we denote ${\mu }_t:={L}(X_t)$ the law of X_t.

3. Solving Mean-Field Control Problems Using Deep Learning Techniques

Due to the non-Markovian structure, the above mean-field optimal control problem (2.3 and 2.4) is difficult to solve either analytically or numerically. Here we propose two algorithms together with four approaches to tackle the above problem based on deep learning techniques. We would like to use two types of neural networks, one is called the feedforward neural network, and the other one is called Long Short-Term Memory (LSTM) network.

For a feedforward neural network, we first define the set of layers ${𝕄}_{d,h}^{\rho }$, for x ∈ ℝ^d, as

$\begin{array}{ll}{𝕄}_{d,h}^{\rho }:=\left\{M:{ℝ}^d\to {ℝ}^h|M(x)=\rho (Ax+b),A\in {ℝ}^{h\times d},b\in {ℝ}^h\right\}.& (3.1)\end{array}$

d is called input dimension, h is known as the number of hidden neurons, A ∈ ℝ^h×d is the weight matrix, b ∈ ℝ^h is the bias vector, and ρ is called the activation function. The following activation functions will be used in this paper, for some x ∈ ℝ,

$\begin{array}{c}{\rho }_{ReLU}(x):=x^{+}=max(0,x);{\rho }_s(x):=\frac{1}{1+e^{-x}};\\ {\rho }_{tanh}(x):=tanh(x);{\rho }_{Id}(x):=x.\end{array}$

Then feedforward neural network is defined as a composition of layers, so that the set of feedforward neural networks with l hidden layers we use in this paper is defined as

$\begin{array}{l}\mathbb{N}{\mathbb{N}}_{d_1,d_2}^l=\{\tilde{M}:{ℝ}^{d_1}\to {ℝ}^{d_2}|\tilde{M}=M_l◦\cdots ◦M_1◦M_0,\\ M_0\in {𝕄}_{d_1,h_1}^{{\rho }_{ReLU}},M_l\in {𝕄}_{h_l,d_2}^{{\rho }_{Id}},M_i\in {𝕄}_{h_i,h_{i+1}}^{{\rho }_{ReLU}},h_{·}\in {\mathbb{Z}}^{+},i=1,\dots ,l-1\}.& (3.2)\end{array}$

The LSTM network is one of RNN architectures, which are powerful for capturing long-range dependence of the data. It is proposed in Hochreiter and Schmidhuber [14], and it is designed to solve the shrinking gradient effects which basic RNN often suffers from. The LSTM network is a chain of cells. Each LSTM cell is composed of a cell state, which contains information, and three gates, which regulate the flow of information. Mathematically, the rule inside tth cell follows,

$\begin{array}{l}{\Gamma }_{f_t}={\rho }_s(A_f{x}_t+U_f{a}_{t-1}+b_f),\\ {\Gamma }_{i_t}={\rho }_s(A_i{x}_t+U_i{a}_{t-1}+b_i),\\ {\Gamma }_{o_t}={\rho }_s(A_o{x}_t+U_o{a}_{t-1}+b_o),\\ c_t={\Gamma }_{f_t}\odot c_{t-1}+{\Gamma }_{i_t}\odot {\rho }_{tanh}(A_c{x}_t+U_c{a}_{t-1}+b_c),\\ a_t={\Gamma }_{o_t}\odot {\rho }_{tanh}(c_t),& (3.3)\end{array}$

where the operator ⊙ denotes the Hadamard product. $({\Gamma }_{f_t},{\Gamma }_{i_t},{\Gamma }_{o_t})\in {ℝ}^h\times {ℝ}^h\times {ℝ}^h$ represents forget gate, input gate and output gate, respectively, h refers the number of hidden neurons. $x_t\in {ℝ}^d$ is the input vector with d features. $a_t\in {ℝ}^h$ is known as the output vector with initial value a₀ = 0, and $c_t\in {ℝ}^h$ is known as the cell state with initial value c₀ = 0. $A_{·}\in {ℝ}^{h\times d}$ are the weight matrices connecting input and hidden layers, $U_{·}\in {ℝ}^{h\times h}$ are the weight matrix connecting hidden and output layers, and b ∈ ℝ^h represents bias vector. The weight matrices and bias vectors are shared through all time steps, and are going to be learned during training process by back-propagation through time (BPTT), which can be implemented in Tensorflow platform. Here we define the set of LSTM network up to time t as

$\begin{array}{r}{𝕃𝕊𝕋𝕄}_{d,h,t}=\{M:{({ℝ}^d)}^t\times {ℝ}^h\times {ℝ}^h\to {ℝ}^h\times {ℝ}^h|M(x_0,\dots ,x_t,\\ a_0,c_0)=(a_t,c_t),c_t={\Gamma }_{f_t}\odot c_{t-1}+{\Gamma }_{i_t}\odot {\rho }_{tanh}(A_c{x}_t+U_c{a}_{t-1}+b_c),\\ a_t={\Gamma }_{o_t}\odot {\rho }_{tanh}(c_t),a_0=c_0=0\},& (3.4)\end{array}$

where Γ_f·, Γ_i·, Γ_o· are defined in (3.3).

In particular, we specify the model in a linear-quadratic form, which is inspired by Carmona et al. [2] and Fouque and Zhang [18]. The objective function is defined as

$\begin{array}{ll}J(\alpha )=𝔼\left[{\displaystyle {\int }_0^T}\left(\frac{1}{2}{\alpha }_t^2+\frac{c_f}{2}{(X_t-m_t)}^2\right)dt+\frac{c_t}{2}{(X_T-m_T)}^2\right],& (3.5)\end{array}$

subject to

$\begin{array}{l}d{X}_t=({\alpha }_t-{\alpha }_{t-\tau })dt+\sigma d{W}_t,t\in \left[0,T\right],\\ X_0=x_0,\\ X_t={\alpha }_t=0,t\in [-\tau ,0),& (3.6)\end{array}$

where σ, c_f, c_t > 0 are given constants, and $m_t:={\int }_{ℝ}xd{\mu }_t(x)$ denotes the mean of X at time t, and ${\mu }_t:={L}(X_t)$. In the following subsections, we solve the above problem numerically using two algorithms together with four approaches. The first two approaches are to directly approximate the control by either a LSTM network or a feedforward neural network, and minimize the objective (3.5) using stochastic gradient descent algorithm. The third and fourth approaches are to introduce the adjoint process associated with (3.6), and approximate the adjoint process and the conditional expectation of adjoint process using neural networks.

3.1. Approximating the Optimal Control Using Neural Networks

We first set Δt = T/N = τ/D for some positive integer N. The time discretization becomes

$-\tau =t_{-D}\le t_{-D+1}\le \cdots \le t_0=0=t_0\le t_1\le \cdots \le t_N\le T,$

for t_i − t_i−1 = Δt, i ∈ {−D + 1, ⋯ , 0, ⋯ , N − 1, N}. The discretized SDDE (3.6) according to Euler-Maruyama scheme now reads

$\begin{array}{ll}{X}_{t_{i+1}}=X_{t_i}+({\alpha }_{t_i}-{\alpha }_{t_{i-D}})\Delta t+\sigma \sqrt{\Delta t}\Delta W_{t_i},\text{ for }i\in \left\{0,\cdots ,N-1\right\},& (3.7)\end{array}$

where (ΔW_ti)_{0≤i≤N−1} are independent, normal distributed sequence of random variables with mean 0 and variance 1.

First, from the definition of open loop control, and due to non-Markovian feature of (3.6), the open-loop optimal control is a function of the path of the Brownian motions up to time t, i.e., α(t, (W_s)_0≤s≤t). We are able to describe this dependency by a LSTM network by parametrizing the control as a function of current time and the discretized increments of Brownian motion path, i.e.,

$\begin{array}{l}(a_{t_i},c_{t_i})={\phi }^1(t_i,{(\Delta W_s)}_{t_0\le s\le t_i}|{\Phi }^1)\text{ for }{\phi }^1\in {𝕃𝕊𝕋𝕄}_{2,h_1,t}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em} \hspace{0.17em}\hspace{1em}and }{\Phi }^1=(A_f,A_i,A_o,A_c,U_f,U_i,U_o,U_c,b_f,b_i,b_o,b_c),\\ \alpha (t_i,{(W_s)}_{t_0\le s\le t_i})\approx {\psi }^1(a_{t_i}|{\text{Ψ}}^1)\text{ for }{\psi }^1\in {𝕄}_{h_1,1}^{Id}\text{ and }{\text{Ψ}}^1=(A,b),& (3.8)\end{array}$

for some $h_1\in {\mathbb{Z}}^{+}$. We remark that the last dense layer is used to match the desired output dimension.

The second approach is again directly approximate the control but with a feedforward neural network. Due to the special structure of our model, where the mean of dynamic in (3.6) is constant, the mean field control problem coincides with the mean field game problem. In Fouque and Zhang [18], authors solved the associated mean field game problem using infinite dimensional PDE approach, and found that the optimal control is a function of current state and the past of control. Therefore, the feedforward neural network with l layers, which we use to approximate the optimal control, is defined as

$\begin{array}{l}{\alpha }_{t_i}(X_{t_i},{({\alpha }_s)}_{t_{i-D}\le s<t_i})\approx {\psi }^2(X_{t_i},{({\alpha }_s)}_{t_{i-D}\le s\le t_i}|{\text{Ψ}}^2)\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} for }{\psi }^2\in {\mathbb{N}\mathbb{N}}_{D+1,1}^l,{\text{Ψ}}^2=(A_0,b_0,\dots ,A_l,b_l).& (3.9)\end{array}$

From Monte Carlo algorithm, and trapezoidal rule, the objective function (3.5) now becomes

$\begin{array}{l}J=\frac{1}{M}{\displaystyle \sum _{j=1}^M}[(\frac{1}{2}{({\alpha }_{t_0}^{(j)})}^2+\frac{c_f}{2}{(X_{t_0}^{(j)}-{\bar{X}}_{t_0})}^2+\\ {\displaystyle \sum _{i=1}^{N-1}}({({\alpha }_{t_i}^{(j)})}^2+c_f{(X_{t_i}^{(j)}-{\bar{X}}_{t_i})}^2)\\ +\frac{1}{2}{({\alpha }_{t_N}^{(j)})}^2+\frac{c_f}{2}{(X_{t_N}^{(j)}-{\bar{X}}_{t_N})}^2)\frac{\Delta t}{2}+\frac{c_t}{2}{(X_{t_N}^{(j)}-{\bar{X}}_{t_N})}^2],& (3.10)\end{array}$

where M denotes the number of realizations and $\bar{X}:=\frac{1}{M}\sum _{j=1}^M{X}^{(j)}$ denotes the sample mean. After plugging in the neural network either given by (3.8) or (3.9), the optimization problem becomes to find the best set of parameters either (Φ¹, Ψ¹) or Ψ² such that the objective J(Φ¹, Ψ¹) or J(Ψ²) is minimized with respect to those parameters.

The algorithm works as follows:

ALGORITHM 1

Algorithm 1: Algorithms for solving mean field control problem with delay by directly approximating the optimal control using neural networks

In the following graphics, we choose x₀ = 0, c_f = c_t = 1, σ = 1, T = 10, τ = 4, Δt = 0.1, M = 4,000. For approach 1, the neural network φ ∈ ℕℕ, which is defined in (3.2), is composed of 3 hidden layers with d₁ = 42, h₁ = 64, h₂ = 128, h₃ = 64, d₂ = 1. For approach 2, the LSTM network φ ∈ 𝕃𝕊𝕋𝕄, which is defined in (3.4), consists of 128 hidden neurons. For a specific representative path, the underlying Brownian motion paths are approximately the same for different approaches. Figure 1 compares one representative optimal trajectory of the state dynamic and the control, and they coincide. Figure 2 plots the sample average of optimal trajectory of the state dynamic and the control, which are trajectories of approximately 0, and this is the same as the theoretical mean.

FIGURE 1

Figure 1. On the left, we compare one representative optimal trajectory of (X_ti)_t0≤ti≤tN. The plot on the right show the comparison of one representative optimal trajectory of (α_ti)_t0≤ti≤tN between approach 1 and approach 2.

FIGURE 2

Figure 2. On the left, we compare the sample mean of optimal trajectories of (X_ti)_t0≤ti≤tN. The plot on the right show the comparison of sample mean of trajectories of optimal control (α_ti)_t0≤ti≤tN between approach 1 and approach 2.

3.2. Approximating the Adjoint Process Using Neural Networks

The third and fourth approaches are based on numerically solving the MV-FABSDE system using LSTM network and feedforward neural networks. From section 4, we derive the adjoint process, and prove the sufficient and necessary parts of stochastic maximum principle. From (4.7), we are able to write the backward stochastic differential equation associated to (3.6) as,

$\begin{array}{ll}d{Y}_t=-c_f(X_t-m_t)dt+Z_{t}d{W}_t,t\in \left[0,T\right],& (3.11)\end{array}$

with terminal condition Y_T = c_t(X_T − m_T), and Y_s = 0 for s ∈ (T, T + τ]. The optimal control ${({\widehat{\alpha }}_t)}_{0\le t\le T}$ can be obtained in terms of the adjoint process Y_t from the maximum principle, and it is given by

${\widehat{\alpha }}_t=-Y_t+𝔼\left[Y_{t+\tau }|{{F}}_t\right].$

From the Euler-Maruyama scheme, the discretized version of (3.6) and (3.11) now reads, for i = {0, …, N − 1},

$\begin{array}{ll}{X}_{t_{i+1}}=X_{t_i}+({\widehat{\alpha }}_{t_i}-{\widehat{\alpha }}_{t_{i-D}})\Delta t+\sigma \sqrt{\Delta t}\Delta W_{t_i},\text{ where }\Delta W_{t_i}~N(0,1).& (3.12)\end{array}$

$\begin{array}{ll}{Y}_{t_{i+1}}=Y_{t_i}-c_f(X_{t_i}-{\bar{X}}_{t_i})\Delta t+Z_{t_i}\sqrt{\Delta t}\Delta W_{t_i},& (3.13)\end{array}$

where we use the sample average ${\bar{X}}_t=\frac{1}{M}\sum _{j=1}^M{X}_t^{(j)}$ to approximate the expectation of X_t. In order to solve the above MV-FABSDE system, we need to approximate ${(Y_{t_i},𝔼[Y_{t_{i+D}}|{{F}}_{t_i}],Z_{t_i})}_{0\le t_i\le t_N}$.

The third approach consists of approximating ${(Y_{t_i},𝔼[Y_{t_{i+D}}|{{F}}_{t_i}],Z_{t_i})}_{0\le t_i\le t_N}$ using three LSTM networks as functions of current time and the discretized path of Brownian motions, respectively, i.e.,

$\begin{array}{l}(a_{t_i}^Y,c_{t_i}^Y)={\phi }^Y(t_i,{(\Delta W_s)}_{t_0\le s\le t_i}|{\Phi }^Y)\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for }{\phi }^Y\in {𝕃𝕊𝕋𝕄}_{2,h_Y,t_i}\text{ and }\\ {\Phi }^Y=(A_f^Y,A_i^Y,A_o^Y,A_c^Y,U_f^Y,U_i^Y,U_o^Y,U_c^Y,b_f^Y,b_i^Y,b_o^Y,b_c^Y),\\ Y_{t_i}\approx {\psi }^Y(a_{t_i}^Y|{\text{Ψ}}^Y)\text{ for }{\psi }^Y\in {𝕄}_{h_Y,1}^{Id}\text{ and }{\text{Ψ}}^Y=(A^Y,b^Y),\end{array}$

$\begin{array}{l}(a_{t_i}^{𝔼Y},c_{t_i}^{𝔼Y})={\phi }^{𝔼Y}(t_i,{(\Delta W_s)}_{t_0\le s\le t_i}|{\Phi }^{𝔼Y})\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for }{\phi }^{𝔼Y}\in {𝕃𝕊𝕋𝕄}_{2,h_{𝔼Y},t_i}\text{ and }\\ {\Phi }^{𝔼Y}=(A_f^{𝔼Y},A_i^{𝔼Y},A_o^{𝔼Y},A_c^{𝔼Y},U_f^{𝔼Y},U_i^{𝔼Y},U_o^{𝔼Y},U_c^{𝔼Y},b_f^{𝔼Y},b_i^{𝔼Y},b_o^{𝔼Y},b_c^{𝔼Y}),\\ 𝔼\left[Y_{t_{i+D}}|{{F}}_{t_i}\right]\approx {\psi }^{𝔼Y}(a_{t_i}^{𝔼Y}|{\text{Ψ}}^{𝔼Y})\text{ for }{\psi }^{𝔼Y}\in {𝕄}_{h_{𝔼Y},1}^{Id}\text{ and }{\text{Ψ}}^{𝔼Y}=(A^{𝔼Y},b^{𝔼Y}),\end{array}$

$\begin{array}{l}(a_{t_i}^Z,c_{t_i}^Z)={\phi }^Z(t_i,{(\Delta W_s)}_{t_0\le s\le t_i}|{\Phi }^Z)\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for }{\phi }^Z\in {𝕃𝕊𝕋𝕄}_{2,h_Z,t_i}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}and }{\Phi }^Z=(A_f^Z,A_i^Z,A_o^Z,A_c^Z,U_f^Z,U_i^Z,U_o^Z,U_c^Z,b_f^Z,b_i^Z,b_o^Z,b_c^Z),\\ Z_{t_i}\approx {\psi }^Z(a_{t_i}^Z|{\text{Ψ}}^Z)\text{ for }{\psi }^Z\in {𝕄}_{h_Z,1}^{Id}\text{ and }{\text{Ψ}}^Z=(A^Z,b^Z),& (3.14)\end{array}$

for some $h_Y,h_{EY},h_Z\in {\mathbb{Z}}^{+}$. Again, the last dense layers are used to match the desired output dimension.

Since approach 3 consists of three neural networks with large number of parameters, which is hard to train in general, we would like to make the following simplification in approach 4 for approximating ${(Y_{t_i},𝔼[Y_{t_{i+D}}|{{F}}_{t_i}],Z_{t_i})}_{t_0\le t_i\le t_N}$ via combination of one LSTM network and three feedforward neural networks. Specifically,

$\begin{array}{l}(a_{t_i},c_{t_i})=\phi (t_i,{(\Delta W_s)}_{t_0\le s\le t_i}|\Phi )\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for }\phi \in {𝕃𝕊𝕋𝕄}_{2,h,t_i}\text{ and}\\ \Phi =(A_f,A_i,A_o,A_c,U_f,U_i,U_o,U_c,b_f,b_i,b_o,b_c),\\ Y_{t_i}\approx {\psi }^Y(a_{t_i}|{\text{Ψ}}^Y)\text{ for }{\psi }^Y\in {\mathbb{N}\mathbb{N}}_{h,1}^l\text{ and}\\ {\text{Ψ}}^Y=(A_0^Y,b_0^Y,\dots ,A_l^Y,b_l^Y),𝔼\left[Y_{t_{i+D}}|{{F}}_{t_i}\right]\\ \approx {\psi }^{𝔼Y}(a_{t_i}|{\text{Ψ}}^{𝔼Y})\text{ for }{\psi }^{𝔼Y}\in {\mathbb{N}\mathbb{N}}_{h,1}^l\text{ and }\\ {\text{Ψ}}^{𝔼Y}=(A_0^{𝔼Y},b_0^{𝔼Y},\dots ,A_l^{𝔼Y},b_l^{𝔼Y}),\\ Z_{t_i}\approx {\psi }^Z(a_{t_i}|{\text{Ψ}}^Z)\text{ for }{\psi }^Z\in {\mathbb{N}\mathbb{N}}_{h,1}^l\text{ and }\\ {\text{Ψ}}^Z=(A_0^Z,b_0^Z,\dots ,A_l^Z,b_l^Z).& (3.15)\end{array}$

In words, the algorithm works as follows. We first initialize the parameters (Θ^Y, Θ^EY, Θ^Z) = ((Φ^Y, Ψ^Y), (Φ^EY, Ψ^EY), (Φ^Z, Ψ^Z)) either in (3.14) or (Θ^Y, Θ^EY, Θ^Z) = ((Φ, Ψ^Y), Ψ^EY), Ψ^Z) in (3.15). At time 0, X₀ = x₀, $(Y_0,𝔼[Y_{t_D}|{{F}}_0],Z_0)\approx ({\phi }_0^Y({\Theta }^Y),{\phi }_0^{𝔼Y}({\Theta }^{𝔼Y}),{\phi }_0^Z({\Theta }^Z))$ for some network (φ^Y, φ^EY, φ^Z) given by either (3.14) or (3.15), and α₀ = −Y_t0+E[Y_tD|${{F}}_{t_0}$]. Next, we update X_ti+1 and Y_ti+1 according to (3.12), and the solution to the backward equation at t_i+1 is denoted by Ỹ_ti+1. In the meantime, Y_ti+1 is also approximated by a neural network. In such case, we refer to Ỹ_· as the label, and Y_· given by the neural network as the prediction. We would like to minimize the mean square error between these two. At time T, Y_tN is also supposed to match $c_t(X_{t_N}-{\bar{X}}_{t_N})$, from the terminal condition of (3.11). In addition, the conditional expectation $𝔼[Y_{t_{i+D}}|{{F}}_{t_i}]$ given by a neural network should be the best predictor of Ỹ_ti+D, which implies that we would like to find the set of parameters Θ^EY such that $𝔼[({Ỹ}_{t_{i+D}}-{\phi }_{t_i}^{EY}({\Theta }^{EY})){}^2]$ is minimized for all t_i ∈ {t₀, …, t_N−D}. Therefore, for M samples, we would like to minimize two objective functions L₁ and L₂ defined as

$\begin{array}{l}{L}_1({\Theta }^Y,{\Theta }^Z)=\frac{1}{M}[{\displaystyle \sum _{j=1}^M}{\displaystyle \sum _{i=0}^N}{({\phi }_{t_i}^{Y,(j)}-{Ỹ}_{t_i}^{(j)})}^2\\ +{\displaystyle \sum _{j=1}^M}{\left({\phi }_{t_N}^{Y,(j)}-c_t(X_{t_N}^{(j)}-{\bar{X}}_{t_N})\right)}^2],\\ L_2({\Theta }^{𝔼Y})=\frac{1}{M}{\displaystyle \sum _{j=1}^M}{\displaystyle \sum _{i=0}^{N-D}}{\left({\phi }_{t_i}^{𝔼Y,(j)}-{Ỹ}_{t_{i+D}}^{(j)}\right)}^2.& (3.16)\end{array}$

The algorithm works as follows:

ALGORITHM 2

Algorithm 2: Algorithms for solving mean field control problem with delay according to MV-FABSDE

Again, in the following graphics, we choose x₀ = 0, c_f = c_t = 1, σ = 1, T = 10, τ = 4, Δt = 0.1, M = 4, 000. In approach 3, each of the three LSTM networks approximating $Y_{t_i},𝔼[Y_{t_i+D}|{{F}}_{t_i}]$ and Z_ti consists of 128 hidden neurons, respectively. In approach 4, the LSTM consists of 128 hidden neurons, and each of the feedforward neural networks has parameters d₁ = 128, h₁ = 64, h₂ = 128, h₃ = 64, d₂ = 1. For a specific representative path, the underlying Brownian motion paths are the same for different approaches. Figure 3 compares one representative optimal trajectory of the state dynamic and the control via two approaches, and they coincide. Figure 4 plots the sample average of optimal trajectory of the dynamic and the control, which are trajectories of 0, which is the same as the theoretical mean. Comparing to Figures 1, 2, as well as based on numerous experiments, we find that given a path of Brownian motion, the two algorithms would yield similar optimal trajectory of state dynamic and similar path for the optimal control. From Figure 6, the loss L₁ as defined in (3.16) becomes approximately 0.02 in 1,000 epochs for both approach 3 and approach 4. This can also be observed from Figure 5, since the red dash line and the blue solid line coincide for both left and right graphs. In addition, from the righthand side of Figure 6, we observe the loss L₂ as defined in (3.16) converges to 50 after 400 epochs. This is due to the fact that the conditional expectation can be understood as an orthogonal projection. Figure 7 plots 64 sample paths of the process (Z_ti)_t0≤ti≤tN, which seems to be a deterministic function since σ is constant in this example. Finally, Figure 8 shows the convergence of the value function as number of epochs increases. Both algorithms arrive approximately at the same optimal value which is around 6 after 400 epochs. This confirms that the out control problem has a unique solution. In section 5, we show that the MV-FASBDE system is uniquely solvable. It is also observable that the first algorithm converges faster than the second one, since it directly paramerizes the control using one neural network, instead of solving the MV-FABSDE system, which uses three neural networks.

FIGURE 3

Figure 3. On the left, we compare one representative optimal trajectory of (X_ti)_t0≤ti≤tN. The plot on the right shows the comparison of one representative optimal trajectory of (α_ti)_t0≤ti≤tN between approach 3 and approach 4.

FIGURE 4

Figure 4. On the left, we compare the sample mean of optimal trajectories of (X_ti)_t0≤ti≤tN. The plot on the right shows the comparison of sample mean of trajectories of optimal control (α_ti)_t0≤ti≤tN between approach 3 and approach 4.

FIGURE 5

Figure 5. Plots of representative trajectories of [${(Y_{t_i})}_{t_0\le t_i\le t_N},{({Ỹ}_{t_i})}_{t_0\le t_i\le t_N},{(𝔼[Y_{t_{i+D}}|{{F}}_{t_i}])}_{t_0\le t_i\le t_N}$], from approach 3 (one the left) and from approach 4 (on the right).

FIGURE 6

Figure 6. Convergence of loss L₁ (on the left) and convergence of loss L₂ (on the right) as defined in (3.16) from approach 3 and approach 4.

FIGURE 7

Figure 7. 64 trajectories of (Z_ti)_t0≤ti≤tN based on approach 3 (on the left) and approach 4 (on the right).

FIGURE 8

Figure 8. Comparison convergence of objective values as in (3.10) among four approaches.

3.3. Numerically Solving the Optimal Control Problem With No Delay

Since the algorithms we proposed embrace the case with no delay, we illustrate the comparison between numerical results and the analytical results. By letting τ > T we obtain α_t−τ = 0 in (3.6), and we aim at solving the following linear-quadratic mean-field control problem by minimizing

$\begin{array}{ll}J(\alpha )=𝔼\left[{\displaystyle {\int }_0^T}\left(\frac{1}{2}{\alpha }_t^2+\frac{c_f}{2}{(X_t-m_t)}^2\right)dt+\frac{c_t}{2}{(X_T-m_T)}^2\right],& (3.17)\end{array}$

subject to

$\begin{array}{l}d{X}_t={\alpha }_{t}dt+\sigma d{W}_t,t\in \left[0,T\right],\\ X_0=x_0.& (3.18)\end{array}$

Again, from section 4, we find the optimal control

${\widehat{\alpha }}_t=-Y_t,$

where (Y_t, Z_t) is the solution of the following adjoint process,

$\begin{array}{ll}d{Y}_t=-c_f(X_t-m_t)dt+Z_{t}d{W}_t.& (3.19)\end{array}$

Next, we make the ansatz

$\begin{array}{ll}{Y}_t={\varphi }_t(X_t-m_t),& (3.20)\end{array}$

for some deterministic function ϕ_t, satisfying the terminal condition ϕ_T = c_t. Differentiating the ansatz, the backward equation should satisfy

$\begin{array}{ll}d{Y}_t=({\stackrel{.}{\varphi }}_t-{\varphi }_t^2)(X_t-m_t)dt+{\varphi }_t\sigma d{W}_t,& (3.21)\end{array}$

where ${\stackrel{\bullet }{\varphi }}_t$ denotes the time derivative of ϕ_t. Comparing with (3.19), and identifying the drift and volatility term, ϕ_t must satisfy the scalar Riccati equation,

$\begin{array}{ll}\{\begin{array}{l}{\stackrel{.}{\varphi }}_t={\varphi }_t^2-c_f,\\ {\varphi }_T=c_t,\end{array}& (3.22)\end{array}$

and the process Z_t should satisfy

$\begin{array}{ll}{Z}_t={\varphi }_t\sigma ,& (3.23)\end{array}$

which is deterministic. If we choose x₀ = 0, c_f = c_t = 1, T = 10, ϕ_t = 1 solves the Riccati equation (3.22), so that Z_t = 1, ∀t ∈ [0, T], and from (3.20), the optimal control satisfies

$\begin{array}{ll}{\widehat{\alpha }}_t=-(X_t-m_t).& (3.24)\end{array}$

Numerically, we apply the two deep learning algorithms proposed in the previous section. The first algorithm directly approximates the control. According to the open loop formulation, we set

$\begin{array}{l}(a_{t_i},c_{t_i})=\phi (t_i,{(\Delta W_s)}_{t_0\le s\le t_i}|\Phi )\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for }\phi \in {𝕃𝕊𝕋𝕄}_{2,h,t}\text{ and}\\ \Phi =(A_f,A_i,A_o,A_c,U_f,U_i,U_o,U_c,b_f,b_i,b_o,b_c),\\ \alpha (t_i,{(W_s)}_{t_0\le s\le t_i})\approx \psi (a_{t_i}|\text{Ψ})\text{ for }\psi \in {𝕄}_{h,1}^{Id}\text{ and Ψ}=(A,b),& (3.25)\end{array}$

for some h ∈ ℤ⁺. We remark that the last dense layer is used to match the desired output dimension. The second algorithm numerically solves the forward backward system as in (3.18) and (3.19). From the ansatz (3.20) and the Markovian feature, we approximate (Y_t, Z_t)_0≤t≤T using two feedforward neural networks, i.e.,

$\begin{array}{l}{Y}_{t_i}\approx {\psi }^1(t_i,X_{t_i}|{\text{Ψ}}^1)\text{ for }{\psi }^1\in {\mathbb{N}\mathbb{N}}_{2,1}^l,{\text{Ψ}}^1=(A_0^1,b_0^1,\dots ,A_l^1,b_l^1);\\ Z_{t_i}\approx {\psi }^2(t_i,X_{t_i}|{\text{Ψ}}^2)\text{ for }{\psi }^2\in {\mathbb{N}\mathbb{N}}_{2,1}^l,{\text{Ψ}}^2=(A_0^2,b_0^2,\dots ,A_l^2,b_l^2).& (3.26)\end{array}$

Figure 9 shows the representative optimal trajectory of ${(X_{t_i}-{\bar{X}}_{t_i})}_{t_0\le t_i\le t_N}$ and (α_ti)_t0≤ti≤tN from both algorithms, which are exactly the same. The symmetry feature can be seen from the computation (3.24). On the left of Figure 10 confirms the mean of the processes ${(X_{t_i}-{\bar{X}}_{t_i})}_{t_0\le t_i\le t_N}$ and (α_ti)_t0≤ti≤tN are 0 from both algorithms. The right picture of Figure 10 plots the data points of x against α, and we can observe that the optimal control α is linear in x as a result of (3.24), and the slope tends to be - 1, since ϕ = 1 solves the scalar Riccati equation (3.22). Finally, Figure 11 plots representative optimal trajectories of the solution to the adjoint equations (Y_t, Z_t). On the left, we observe that the adjoint process (Y_t)_0≤t≤T matches the terminal condition, and on the right, (Z_t)_0≤t≤T appears to be a deterministic process of value 1, and this matches the result we compute previously.

FIGURE 9

Figure 9. Representative optimal trajectory of ${(X_{t_i}-{\bar{X}}_{t_i})}_{t_0\le t_i\le t_N}$ and (α_ti)_t0≤ti≤tN from algorithm 1 (on the left) and from algorithm 2 (on the right).

FIGURE 10

Figure 10. The picture on the left plots the sample averages of (X_ti)_t0≤ti≤tN and (α_ti)_t0≤ti≤tN for both algorithm 1 and 2; The plot on the right shows the points of x against α for both algorithms.

FIGURE 11

Figure 11. The plot on the left shows representative trajectories of [(Y_ti)_t0≤ti≤tN, (Ỹ_ti)_t0≤ti≤tN]. The picture on the right plots 64 trajectories of (Z_ti)_t0≤ti≤tN.

4. Stochastic Maximum Principle for Optimality

In this section, we derive the adjoint equation associated to our mean field stochastic control problem (2.3) and (2.4). The necessary and sufficient parts of stochastic maximum principle have been proved for optimality. We assume

(H4.1) b, σ are differentiable with respect to (X_t, μ_t, X_t−τ, μ_t−τ, α_t, α_t−τ); f is differentiable with respect to (X_t, μ_t, X_t−τ, μ_t−τ, α); g is differentiable with respect to (X_T, μ_T). Their derivatives are bounded.

In order to simplify our notations, let θ_t = (X_t, μ_t, α_t). For 0 < ϵ < 1, we denote α^ϵ the admissible control defined by

${\alpha }_t^{ϵ}:={\alpha }_t+ϵ({\beta }_t-{\alpha }_t):={\alpha }_t+ϵ\Delta {\alpha }_t,$

for any (α)_0≤t≤T and (β)_0≤t≤T ∈ 𝔸. $X_t^{ϵ}:=X_t^{{\alpha }^{ϵ}}$ is the corresponding controlled process. We define

$\nabla X_t:={\displaystyle \underset{ϵ\to 0}{lim}}\frac{X_t^{ϵ}-X_t^{\alpha }}{ϵ}$

to be the variation process, which should follow the following dynamic for t ∈ (0, T],

$\begin{array}{l}d\nabla X_t=[{\partial }_{x}b(t,{\theta }_t,{\theta }_{t-\tau })\nabla X_t+{\partial }_{x_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })\nabla X_{t-\tau }\\ +\widetilde{𝔼}\left[{\partial }_{\mu }b(t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_t)\nabla {\tilde{X}}_t\right]+\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]\\ +{\partial }_{\alpha }b(t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_t+{\partial }_{{\alpha }_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{t-\tau }]dt\\ +[{\partial }_x\sigma (t,{\theta }_t,{\theta }_{t-\tau })\nabla X_t+{\partial }_{x_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\nabla X_{t-\tau }\\ +\widetilde{𝔼}\left[{\partial }_{\mu }\sigma (t,{\theta }_t,{\theta }_{t-\tau })(\widetilde{X_t})\nabla \widetilde{X_t}\right]+\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]\\ +{\partial }_{\alpha }\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_t+{\partial }_{{\alpha }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{t-\tau }]d{W}_t& (4.1)\end{array}$

with initial condition ∇X_t = Δα_t = 0, t ∈ [−τ, 0). $({\tilde{X}}_t,\nabla {\tilde{X}}_t)$ is a copy of (X_t, ∇X_t) defined on $(\tilde{\Omega },\tilde{{F}},\tilde{\mathbb{P}})$, where we apply differential calculus on functions of measure, see Carmona and Delarue [3] for detail. ∂_xb, ∂_xτb, ∂_μb, ∂_μτb, ∂_αb, ∂_ατb are derivatives of b with respect to (X_t, X_t−τ, μ_t, μ_t−τ, α_t, α_t−τ), respectively, and we use the same notation for ∂_·σ.

In the meantime, the Gateaux derivative of functional α → J(α) is given by

$\begin{array}{l}{\displaystyle \underset{ϵ\to 0}{lim}}\frac{J({\alpha }^{ϵ})-J(\alpha )}{ϵ}\\ =𝔼\left[{\partial }_{x}g(X_T,{\mu }_T)\nabla X_T+\widetilde{𝔼}\left[{\partial }_{\mu }g(X_T,{\mu }_T)({\tilde{X}}_T)\nabla {\tilde{X}}_T\right]\right]\\ +𝔼{\displaystyle {\int }_0^T}[{\partial }_{x}f({\theta }_t,X_{t-\tau },{\mu }_{t-\tau })\nabla X_t+\widetilde{𝔼}\left[{\partial }_{\mu }f({\theta }_t,X_{t-\tau },{\mu }_{t-\tau })({\tilde{X}}_t)\nabla {\tilde{X}}_t\right]\\ +{\partial }_{x_{\tau }}f({\theta }_t,X_{t-\tau },{\mu }_{t-\tau })\nabla X_{t-\tau }+\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}f({\theta }_t,X_{t-\tau },{\mu }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]\\ +{\partial }_{\alpha }f({\theta }_t,X_{t-\tau },{\mu }_{t-\tau })(\Delta {\alpha }_t)]dt& (4.2)\end{array}$

In order to determine the adjoint backward equation of (Y_t, Z_t)_0≤t≤T associated to (2.3), we assume it is of the following form:

$\begin{array}{l}d{Y}_t=-{\phi }_{t}dt+Z_{t}d{W}_t,t\in \left[0,T\right],\\ Y_T={\partial }_{x}g(X_T,{\mu }_T)+\widetilde{𝔼}\left[{\partial }_{\mu }g(X_T,{\mu }_T)({\tilde{X}}_T)\right],\\ Y_t=Z_t=0,t\in (T,T+\tau ]& (4.3)\end{array}$

Next, we apply integration by part to ∇X_t and Y_t. It yields

$\begin{array}{l}d(\nabla X_t{Y}_t)=Y_t[{\partial }_{x}b(t,{\theta }_t,{\theta }_{t-\tau })\nabla X_t+{\partial }_{x_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })\nabla X_{t-\tau }\\ +\widetilde{𝔼}\left[{\partial }_{\mu }b(t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_t)\nabla {\tilde{X}}_t\right]+\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]\\ +{\partial }_{\alpha }b(t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_t+{\partial }_{{\alpha }_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{t-\tau }]dt\\ +Y_t[{\partial }_x\sigma (t,{\theta }_t,{\theta }_{t-\tau })+{\partial }_{x_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\nabla X_{t-\tau }\\ +\widetilde{𝔼}\left[{\partial }_{\mu }\sigma (t,{\theta }_t,{\theta }_{t-\tau })(\widetilde{X_t})\nabla \widetilde{X_t}\right]+\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]\\ +{\partial }_{\alpha }\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_t+{\partial }_{{\alpha }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{t-\tau }]d{W}_t-{\phi }_t\nabla X_{t}dt\\ +\nabla X_t{Z}_{t}d{W}_t+Z_t[{\partial }_x\sigma (t,{\theta }_t,{\theta }_{t-\tau })\nabla X_t\\ +{\partial }_{x_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\nabla X_{t-\tau }+\widetilde{𝔼}\left[{\partial }_{\mu }\sigma (t,{\theta }_t,{\theta }_{t-\tau })(\widetilde{X_t})\nabla \widetilde{X_t}\right]\\ +\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]+{\partial }_{\alpha }\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_t\\ +{\partial }_{{\alpha }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{t-\tau }]dt\end{array}$

We integrate from 0 to T, and take expectation to get

$\begin{array}{l}𝔼\left[\nabla X_T{Y}_T\right]\\ =𝔼{\displaystyle {\int }_0^T}{Y}_t[{\partial }_{x}b(t,{\theta }_t,{\theta }_{t-\tau })\nabla X_t+{\partial }_{x_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })\nabla X_{t-\tau }\\ +\widetilde{𝔼}\left[{\partial }_{\mu }b(t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_t)\nabla {\tilde{X}}_t\right]+\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]\\ +{\partial }_{\alpha }b(t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{\tau }+{\partial }_{{\alpha }_{\tau }}b(t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{t-\tau }]dt\\ -𝔼{\displaystyle {\int }_0^T}{\phi }_t\nabla X_{t}dt+𝔼{\displaystyle {\int }_0^T}{Z}_t[{\partial }_x\sigma (t,{\theta }_t,{\theta }_{t-\tau })\nabla X_t\\ +{\partial }_{x_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\nabla X_{t-\tau }+\widetilde{𝔼}\left[{\partial }_{\mu }\sigma (t,{\theta }_t,{\theta }_{t-\tau })(\widetilde{X_t})\nabla \widetilde{X_t}\right]\\ +\widetilde{𝔼}\left[{\partial }_{{\mu }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })({\tilde{X}}_{t-\tau })\nabla {\tilde{X}}_{t-\tau }\right]+{\partial }_{\alpha }\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_t\\ +{\partial }_{{\alpha }_{\tau }}\sigma (t,{\theta }_t,{\theta }_{t-\tau })\Delta {\alpha }_{t-\tau }]dt& (4.4)\end{array}$

Using the fact that Y_t = Z_t = 0 for t ∈ (T, T + τ], we are able to make a change of time, and by Fubini's theorem, so that (4.4) becomes

$\begin{array}{l}𝔼\left[\nabla X_T{Y}_T\right]=𝔼{\displaystyle {\int }_0^T}(Y_t{\partial }_{x}b(t,{\theta }_t,{\theta }_{t-\tau })+Y_{t+\tau }{\partial }_{x_{\tau }}b(t+\tau ,{\theta }_{t+\tau },{\theta }_t)\\ +\widetilde{𝔼}\left[{\partial }_{\mu }b(t,{\tilde{\theta }}_t,{\tilde{\theta }}_{t-\tau })(X_t)\right]\\ +\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}b(t+\tau ,{\tilde{\theta }}_{t+\tau },{\tilde{\theta }}_t)(X_t))\nabla X_{t}dt\\ +𝔼{\displaystyle {\int }_0^T}({\partial }_{\alpha }b(t,{\theta }_t,{\theta }_{t-\tau })+{\partial }_{{\alpha }_{\tau }}b(t+\tau ,{\theta }_{t+\tau },{\theta }_t))\Delta {\alpha }_{t}dt\\ -𝔼[{\displaystyle {\int }_0^T}{\phi }_t\nabla X_t]dt+𝔼{\displaystyle {\int }_0^T}(Z_t{\partial }_x\sigma (t,{\theta }_t,{\theta }_{t-\tau })\\ +Z_{t+\tau }{\partial }_{x_{\tau }}\sigma (t+\tau ,{\theta }_{t+\tau },{\theta }_t)+\widetilde{𝔼}\left[{\partial }_{\mu }\sigma (t,{\tilde{\theta }}_t,{\tilde{\theta }}_{t-\tau })(X_t)\right]\\ +\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}\sigma (t+\tau ,{\tilde{\theta }}_{t+\tau },{\tilde{\theta }}_t)(X_t))\nabla X_{t}dt\\ +𝔼{\displaystyle {\int }_0^T}({\partial }_{\alpha }\sigma (t,{\theta }_t,{\theta }_{t-\tau })+{\partial }_{{\alpha }_{\tau }}\sigma (t+\tau ,{\theta }_{t+\tau },{\theta }_t))\Delta {\alpha }_{t}dt& (4.5)\end{array}$

Now we define the Hamiltonian H for $(t,x,\mu ,x_{\tau },{\mu }_{\tau },y,z,\alpha ,{\alpha }_{\tau })\in [0,T]\times ℝ\times {{P}}_2(ℝ)\times ℝ\times {{P}}_2(ℝ)\times ℝ\times ℝ\times A\times A$ as

$\begin{array}{l}H(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha ,{\alpha }_{\tau },y,z)\\ =b(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha ,{\alpha }_{\tau })y+\sigma (t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha ,{\alpha }_{\tau })z\\ +f(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha )& (4.6)\end{array}$

Using the terminal condition of Y_T, and plugging (4.5) into (4.2), and setting the integrand containing ∇X_t to zero, we are able to obtain the adjoint equation is of the following form

$\begin{array}{l}d{Y}_t=-\{{\partial }_{x}H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t,{\alpha }_{t-\tau },Y_t,Z_t)\\ +\widetilde{𝔼}\left[{\partial }_{\mu }H(t,{\tilde{X}}_t,{\mu }_t,{\tilde{X}}_{t-\tau },{\mu }_{t-\tau },{\tilde{\alpha }}_t,{\tilde{\alpha }}_{t-\tau },{Ỹ}_t,{\tilde{Z}}_t)(X_t)\right]\\ +𝔼\left[{\partial }_{x_{\tau }}H(t+\tau ,X_{t+\tau },{\mu }_{t+\tau },X_t,{\mu }_t,{\alpha }_{t+\tau },{\alpha }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t\right]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\tilde{X}}_{t+\tau },{\mu }_{t+\tau },{\tilde{X}}_t,{\mu }_t,{\tilde{\alpha }}_t,{\tilde{\alpha }}_{t-\tau },{Ỹ}_{t+\tau },{\tilde{Z}}_{t+\tau })\\ (X_{t+\tau })]|{{F}}_t]\}dt+Z_{t}d{W}_t\\ Y_T={\partial }_{x}g(X_T,{\mu }_T)+\widetilde{𝔼}\left[{\partial }_{\mu }g({\tilde{X}}_T,{\mu }_T)(X_T)\right].& (4.7)\end{array}$

Theorem 4.1. Let (α_t)_0≤t≤T ∈ A be optimal, (X_t)_0≤t≤T be the associated controlled state, and (_{Yt, Zt)0≤t≤T} be the associated adjoint processes defined in (4.7). For any β ∈ A, and t ∈ [0, T],

$\begin{array}{l}({\partial }_{\alpha }H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t,{\alpha }_{t-\tau },Y_t,Z_t)\\ +𝔼\left[{\partial }_{{\alpha }_{\tau }}H(t+\tau ,X_{t+\tau },{\mu }_{t+\tau },X_t,{\mu }_t,{\alpha }_{t+\tau },{\alpha }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t\right])\\ (\beta -{\alpha }_t)\ge 0a.e.& (4.8)\end{array}$

Proof. Given any (β_t)_0≤t≤T ∈ 𝔸, we perturbate α_t by ϵ(β_t−α_t) and we define ${\alpha }_t^{ϵ}:={\alpha }_t+ϵ({\beta }_t-{\alpha }_t)$ for 0 ≤ ϵ ≤ 1. Using the adjoint process (4.7), and apply integration by parts formula to (∇X_tY_t). Then plug the result into (4.2), and the Hamiltonian H is defined in (4.6). Also, since α is optimal, we have

$\begin{array}{l}0\le {\displaystyle \underset{ϵ\to 0}{lim}}\frac{J({\alpha }^{ϵ})-J(\alpha )}{ϵ}\\ =𝔼{\displaystyle {\int }_0^T}([{\partial }_{\alpha }H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)+\\ 𝔼\left[{\partial }_{{\alpha }_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t\right])({\beta }_t-{\alpha }_t)dt& (4.9)\end{array}$

Now, let $C\in {{F}}_t$ be an arbitrary progressively measurable set, and denote C′ the complement of C. We choose β_t to be ${\beta }_t:=\beta 1_C+{\alpha }_t{1}_{C^{\prime }}$ for any given β ∈ A. Then,

$\begin{array}{l}𝔼{\displaystyle {\int }_0^T}([{\partial }_{\alpha }H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)\\ +𝔼[{\partial }_{{\alpha }_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t])({\beta }_t-{\alpha }_t)1_{C}dt\ge 0,& (4.10)\end{array}$

which implies,

$\begin{array}{c}({\partial }_{\alpha }H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)\\ +𝔼[{\partial }_{{\alpha }_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t])(\beta -{\alpha }_t)\ge 0.\text{ a.e.}& (4.11)\end{array}$

□

Remark 4.2. When we further assume that H is convex in (α_t, α_t−τ), then for any β, β_τ ∈ A in Theorem 4.1, we have

$\begin{array}{r}H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t,{\alpha }_{t-\tau },Y_t,Z_t)\\ \le H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },\beta ,{\beta }_{\tau },Y_t,Z_t),\text{ a.e.}\end{array}$

as a direct consequence of (4.8).

Theorem 4.3. Let (α_t)_0≤t≤T ∈ 𝔸 be an admissible control. Let (X_t)_0≤t≤T be the controlled state, and (Y_t, Z_t)_0≤t≤T be the corresponding adjoint processes. We further assume that for each t, given Y_t and Z_t, the function (x, μ, x_τ, μ_τ, α, α_τ) → H(t, x, μ, x_τ, μ_τ, α, α_τ, Y_t, Z_t), and the function (x, μ) → g(x, μ) are convex. If

$\begin{array}{l}H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t,{\alpha }_{t-\tau },Y_t,Z_t)\\ ={\displaystyle \underset{{\alpha }^{\prime }\in 𝔸}{inf}}H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t^{\prime },{\alpha }_{t-\tau }^{\prime },Y_t,Z_t),& (4.12)\end{array}$

for all t, then (α_t)_0≤t≤T is an optimal control.

Proof. Let ${({\alpha }_t^{\prime })}_{0\le t\le T}\in 𝔸$ be a admissible control, and let ${(X_t^{\prime })}_{0\le t\le T}={(X_t^{{\alpha }^{\prime }})}_{{}_{0\le t\le T}}$ be the corresponding controlled state. From the definition of the objective function as in (2.4), we first use convexity of g, and the terminal condition of the adjoint process Y_t in (4.7), then use the fact that H is convex, and because of (4.12), we have the following

$\begin{array}{l}J(\alpha )-J({\alpha }^{\prime })=𝔼[g(X_T,{\mu }_T)-g(X_T^{\prime },{\mu }_T^{\prime })]\\ +𝔼{\displaystyle {\int }_0^T}[f(t,{\theta }_t,X_{t-\tau },{\alpha }_{t-\tau })-f(t,{\theta }_t^{\prime },X_{t-\tau }^{\prime },{\alpha }_{t-\tau }^{\prime })]dt\\ \le 𝔼[{\partial }_{x}g(X_T,{\mu }_T)(X_T-X_T^{\prime })+\widetilde{𝔼}[{\partial }_{\mu }g(X_T,{\mu }_T)({\tilde{X}}_T)({\tilde{X}}_T-{\tilde{X}}_T^{\prime })]]\\ +𝔼{\displaystyle {\int }_0^T}[f(t,{\theta }_t,X_{t-\tau },{\alpha }_{t-\tau })-f(t,{\theta }_t^{\prime },X_{t-\tau }^{\prime },{\alpha }_{t-\tau }^{\prime })]dt\\ =𝔼[({\partial }_{x}g(X_T,{\mu }_T)+Ẽ[{\partial }_{\mu }g({\tilde{X}}_T,{\mu }_T)(X_T)])(X_T-X_T^{\prime })]\\ +𝔼{\displaystyle {\int }_0^T}[f(t,{\theta }_t,X_{t-\tau },{\alpha }_{t-\tau })-f(t,{\theta }_t^{\prime },X_{t-\tau }^{\prime },{\alpha }_{t-\tau }^{\prime })]dt\\ =𝔼[Y_T(X_T-X_T^{\prime })]+𝔼{\displaystyle {\int }_0^T}[f(t,{\theta }_t,X_{t-\tau },{\alpha }_{t-\tau })\\ -f(t,{\theta }_t^{\prime },X_{t-\tau }^{\prime },{\alpha }_{t-\tau }^{\prime })]dt=𝔼{\displaystyle {\int }_0^T}[(b(t,{\theta }_t,{\theta }_{t-\tau })-b(t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime }))Y_t\\ +(\sigma (t,{\theta }_t,{\theta }_{t-\tau })-\sigma (t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime }))Z_t]dt-𝔼{\displaystyle {\int }_0^T}[({\partial }_{x}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)\\ +\widetilde{𝔼}[{\partial }_{\mu }H(t,{\tilde{\theta }}_t,{\tilde{\theta }}_{t-\tau },{Ỹ}_t,{\tilde{Z}}_t)(X_t)])(X_t-X_t^{\prime })]dt\\ -𝔼{\displaystyle {\int }_0^T}[(𝔼[{\partial }_{x_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\tilde{\theta }}_{t+\tau },{\tilde{\theta }}_t,{Ỹ}_{t+\tau },{\tilde{Z}}_{t+\tau })(X_t)]|{{F}}_t])(X_t-X_t^{\prime })]dt\\ +𝔼{\displaystyle {\int }_0^T}[H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)-H(t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime },Y_t,Z_t)]dt\\ +𝔼{\displaystyle {\int }_0^T}[(b(t,{\theta }_t,{\theta }_{t-\tau })-b(t,{\theta }_t,{\theta }_{t-\tau }))Y_t+(\sigma (t,{\theta }_t,{\theta }_{t-\tau })\\ -\sigma (t,{\theta }_t,{\theta }_{t-\tau }))Z_t]dt\le -𝔼{\displaystyle {\int }_0^T}[{\partial }_{x}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)(X_t-X_t^{\prime })\\ +\widetilde{𝔼}[{\partial }_{\mu }H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)({\tilde{X}}_t)({\tilde{X}}_t-{\widetilde{X_t}}^{\prime })]dt\\ -𝔼{\displaystyle {\int }_0^T}[{\partial }_{x_{\tau }}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)(X_{t-\tau }-X_{t-\tau }^{\prime })\\ +\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)({\tilde{X}}_{t-\tau })({\tilde{X}}_{t-\tau }-{\tilde{X}}_{t-\tau }^{\prime })]]dt\\ +𝔼{\displaystyle {\int }_0^T}[H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)-H(t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime },Y_t,Z_t)]dt\\ \le 𝔼{\displaystyle {\int }_0^T}[{\partial }_{\alpha }H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)({\alpha }_t-{\alpha }_t^{\prime })+{\partial }_{{\alpha }_{\tau }}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)\\ ({\alpha }_{t-\tau }-{\alpha }_{t-\tau }^{\prime })]dt\le 𝔼{\displaystyle {\int }_0^T}({\partial }_{\alpha }H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)\\ +𝔼[{\partial }_{{\alpha }_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t])({\alpha }_t-{\alpha }_t^{\prime })dt\le 0.& (4.13)\end{array}$

□

5. Existence and Uniqueness Result

Given the necessary and sufficient conditions proven in section 4, we use the optimal control ${({\widehat{\alpha }}_t)}_{0\le t\le T}$ defined by

$\begin{array}{l}\widehat{\alpha }(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },Y_t,Z_t,𝔼[Y_{t+\tau }|{{F}}_t],𝔼[Z_{t+\tau }|{{F}}_t])\\ =arg{\displaystyle \underset{\alpha \in A}{min}}H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\alpha }_t,{\alpha }_{t-\tau },Y_t,Z_t),& (5.1)\end{array}$

to establish the solvability result of the McKean–Vlasov FABSDE (2.3) and (4.7) for t ∈ [0, T]:

$\begin{array}{l}d{X}_t=b(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\widehat{\alpha }}_t,{\widehat{\alpha }}_{t-\tau })dt\\ +\sigma (t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\widehat{\alpha }}_t,{\widehat{\alpha }}_{t-\tau })d{W}_t,\\ d{Y}_t=-\{{\partial }_{x}H(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },{\widehat{\alpha }}_t,{\widehat{\alpha }}_{t-\tau },Y_t,Z_t)\\ +\widetilde{𝔼}[{\partial }_{\mu }H(t,{\tilde{X}}_t,{\mu }_t,{\tilde{X}}_{t-\tau },{\mu }_{t-\tau },{\tilde{\widehat{\alpha }}}_t,{\tilde{\widehat{\alpha }}}_{t-\tau },{Ỹ}_t,{\tilde{Z}}_t)(X_t)]\\ +𝔼[{\partial }_{x_{\tau }}H(t+\tau ,X_{t+\tau },{\mu }_{t+\tau },X_t,{\mu }_t,{\widehat{\alpha }}_{t+\tau },{\widehat{\alpha }}_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\tilde{X}}_{t+\tau },{\mu }_{t+\tau },{\tilde{X}}_t,{\mu }_t,{\tilde{\widehat{\alpha }}}_{t+\tau },{\tilde{\widehat{\alpha }}}_t,{Ỹ}_{t+\tau },{\tilde{Z}}_{t+\tau })\\ (X_t)]|{{F}}_t]\}dt+Z_{t}d{W}_t& (5.2)\end{array}$

with initial condition $X_0=x_0;X_t={\widehat{\alpha }}_t=0\mathrm{\text{for}}t\in [-\tau ,0)$ and terminal condition $Y_T={\partial }_{x}g(X_T,{\mu }_T)+\widetilde{𝔼}[{\partial }_{\mu }g({\tilde{X}}_T,{\mu }_T)(X_T)].$ In addition to assumption (H 4.1), we further assume

(H5.1) The drift and volatility functions b and σ are linear in x, μ, x_τ, μ_τ, α, α_τ. For all $(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha ,{\alpha }_{\tau })\in [0,T]\times ℝ\times {{P}}_2(ℝ)\times {{P}}_2(ℝ)\times A\times A$, we assume that

$\begin{array}{r}b(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha ,{\alpha }_{\tau })=b_0(t)+b_1(t)x+{\bar{b}}_1(t)m\\ +b_2(t)x_{\tau }+{\bar{b}}_2(t)m_{\tau }+b_3(t)\alpha +b_4(t){\alpha }_{\tau },\\ \sigma (t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha ,{\alpha }_{\tau })={\sigma }_0(t)+{\sigma }_1(t)x+{\bar{\sigma }}_1(t)m\\ +{\sigma }_2(t)x_{\tau }+{\bar{\sigma }}_2(t)m_{\tau }+{\sigma }_3(t)\alpha +{\sigma }_4(t){\alpha }_{\tau },& (5.3)\end{array}$

for some measurable deterministic functions $b_0,b_1,{\bar{b}}_1,b_2,{\bar{b}}_2,b_3,b_4,{\sigma }_0,{\sigma }_1,{\bar{\sigma }}_1,{\sigma }_2,{\bar{\sigma }}_2,{\sigma }_3,{\sigma }_4$ with values in ℝ bounded by R, and we have used the notation m = ∫xdμ(x) and m_τ = ∫xdμ_τ(x) for the mean of measures μ and μ_τ, respectively.

(H5.2) The derivatives of f and g with respect to (x, x_τ, μ, μ_τ, α) and (x, μ) are Lipschitz continuous with Lipschitz constant L.

(H5.3) The function f is strongly L-convex, which means that for any t ∈ [0, T], any $x,x^{\prime },x_{\tau },x_{\tau }^{\prime }\in ℝ$, any α, α′ ∈ A, any $\mu ,{\mu }^{\prime },{\mu }_{\tau },{\mu }_{\tau }^{\prime }\in {{P}}_2(ℝ)$, any random variables X and X′ having μ and μ′ as distribution, and any random variables X_τ and $X_{\tau }^{\prime }$ having μ_τ and ${\mu }_{\tau }^{\prime }$ as distribution, then

$\begin{array}{l}f(t,x^{\prime },{\mu }^{\prime },x_{\tau }^{\prime },{\mu }_{\tau }^{\prime },{\alpha }^{\prime })-f(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha )\\ -{\partial }_{x}f(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha )(x^{\prime }-x)\\ -{\partial }_{x_{\tau }}f(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha )(x_{\tau }^{\prime }-x_{\tau })\\ -𝔼[{\partial }_{\mu }f(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha )(X)·(X^{\prime }-X)]\\ -𝔼[{\partial }_{{\mu }_{\tau }}f(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha )(X_{\tau })·(X_{\tau }^{\prime }-X_{\tau })]\\ -{\partial }_{\alpha }f(t,x,\mu ,x_{\tau },{\mu }_{\tau },\alpha )({\alpha }^{\prime }-\alpha )\ge \kappa |{\alpha }^{\prime }-\alpha {|}^2.& (5.4)\end{array}$

The function g is also assumed to be L-convex in (x, μ).

Theorem 5.1. Under assumptions (H5.1-H5.3), the McKean–Vlasov FABSDE (5.2) is uniquely solvable.

The proof is based on continuation methods. Let λ ∈ [0, 1], consider the following class of McKean–Vlasov FABSDEs, denoted by MV-FABSDE(λ), for t ∈ [0, T]:

$\begin{array}{l}d{X}_t=(\lambda b(t,{\theta }_t,{\theta }_{t-\tau })+I_t^b)dt+(\lambda \sigma (t,{\theta }_t,{\theta }_{t-\tau })+I_t^{\sigma })d{W}_t,\\ d{Y}_t=-\{\lambda ({\partial }_{x}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)+\widetilde{𝔼}[{\partial }_{\mu }H(t,{\tilde{\theta }}_t,{\tilde{\theta }}_{t-\tau },{Ỹ}_t,{\tilde{Z}}_t)(X_t)]\\ +𝔼[{\partial }_{x_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\tilde{\theta }}_{t+\tau },{\tilde{\theta }}_t,{Ỹ}_{t+\tau },{\tilde{Z}}_{t+\tau })(X_t)]|{{F}}_t])\\ +I_t^f\}dt+Z_{t}d{W}_t,& (5.5)\end{array}$

where we denote θ_t = (X_t, μ_t, α_t), with optimality condition

$\begin{array}{r}{\alpha }_t=\widehat{\alpha }(t,X_t,{\mu }_t,X_{t-\tau },{\mu }_{t-\tau },Y_t,Z_t,𝔼[Y_{t+\tau }|{{F}}_t],\\ 𝔼[Z_{t+\tau }|{{F}}_t]),t\in [0,T],\end{array}$

and with initial condition X₀ = x₀; X_t = α_t = 0 for t ∈ [−τ, 0) and terminal condition

$\begin{array}{l}{Y}_T=\lambda \{{\partial }_{x}g(X_T,{\mu }_T)+\widetilde{𝔼}[{\partial }_{\mu }g({\tilde{X}}_T,{\mu }_T)(X_T)\}+I_T^g,\end{array}$

and Y_t = 0 for t ∈ (T, T + τ], where ${(I_t^b,I_t^{\sigma },I_t^f)}_{0\le t\le T}$ are some square-integrable progressively measurable processes with values in ℝ, and $I_T^g\in L^2(\Omega ,{{F}}_T,\mathbb{P})$ is a square integrable ${{F}}_T$-measurable random variable with value in ℝ.

Observe that when λ = 0, system (5.5) becomes decoupled standard SDE and BSDE, which has an unique solution. When setting $\lambda =1,I_t^b=I_t^{\sigma }=I_t^f=0$ for 0 ≤ t ≤ T, and $I_T^g=0$, we are able to recover the system of (5.2).

Lemma 5.2. Given λ₀ ∈ [0, 1), for any square-integrable progressively measurable processes ${(I_t^b,I_t^{\sigma },I_t^f)}_{0\le t\le T}$, and $I_T^g\in L^2(\Omega ,{{F}}_T,\mathbb{P})$, such that system FABSDE(λ₀) admits a unique solution, then there exists δ₀ ∈ (0, 1), which is independent on λ₀, such that the system MV-FABSDE(λ) admits a unique solution for any λ ∈ [λ₀, λ₀ + δ₀].

Proof. Assuming that $(\check{X},\check{Y},Ž,\check{\alpha })$ are given as an input, for any λ ∈ [λ₀, λ₀ + δ₀], where δ₀ > 0 to be determined, denoting δ: = λ − λ₀ ≤ δ₀, we take

$\begin{array}{l}{I}_t^b\leftarrow \delta [b(t,{\check{\theta }}_t,{\check{\theta }}_{t-\tau })]+I_t^b,\\ I_t^{\sigma }\leftarrow \delta [\sigma (t,{\check{\theta }}_t,{\check{\theta }}_{t-\tau })]+I_t^{\sigma },\\ I_t^f\leftarrow \delta [{\partial }_{x}H(t,{\check{\theta }}_t,{\check{\theta }}_{t-\tau },Y_t,Z_t)+\widetilde{𝔼}[{\partial }_{\mu }H(t,{\widetilde{\stackrel{ˇ}{\theta }}}_t,{\widetilde{\stackrel{ˇ}{\theta }}}_{t-\tau },{\widetilde{\stackrel{ˇ}{Y}}}_t,{\widetilde{Ž}}_t)(X_t)]\\ +𝔼[{\partial }_{x_{\tau }}H(t+\tau ,{\check{\theta }}_{t+\tau },{\check{\theta }}_t,{\check{Y}}_{t+\tau },{Ž}_{t+\tau })|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\widetilde{\stackrel{ˇ}{\theta }}}_{t+\tau },{\widetilde{\stackrel{ˇ}{\theta }}}_t,{\widetilde{\stackrel{ˇ}{Y}}}_{t+\tau },{\widetilde{Ž}}_{t+\tau })({\check{X}}_t)]|{{F}}_t]]+I_t^f,\\ I_T^g\leftarrow \delta [{\partial }_{x}g({\check{X}}_T,{\mu }_T)+\widetilde{𝔼}[{\partial }_{\mu }g({\widetilde{\stackrel{ˇ}{X}}}_T,{\mu }_T)({\check{X}}_T)]+I_T^g.& (5.6)\end{array}$

According to the assumption, let (X, Y, Z) be the solutions of MV-FABSDE(λ₀) corresponding to inputs $(\check{X},\check{Y},Ž)$, i.e., for t ∈ [0, T]

$\begin{array}{l}d{X}_t=({\lambda }_0{b}_t+\delta {\check{b}}_t+I_t^b)dt+({\lambda }_0{\sigma }_t+\delta {\check{\sigma }}_t+I_t^{\sigma })d{W}_t,\\ d{Y}_t=-\{{\lambda }_0({\partial }_x{H}_t+\widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]\\ +𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]+𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\tilde{H}}_{t+\tau }(X_t)]|{{F}}_t])\\ +\delta ({\partial }_x{Ȟ}_t+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]+𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t])+I_t^f\}dt+Z_{t}d{W}_t,& (5.7)\end{array}$

with initial condition, X₀ = x₀, X_s = α_s = 0 for s ∈ [−τ, 0), and terminal condition

$\begin{array}{ll}{Y}_T={\lambda }_0({\partial }_x{g}_T+\widetilde{𝔼}[{\partial }_{\mu }{\tilde{g}}_T(X_T)])+\delta ({\partial }_x{ǧ}_T+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])+I_T^g,& (5.8)\end{array}$

and Y_t = Z_t = 0 for t ∈ (T, T+τ], where we have used simplified notations,

$\begin{array}{l}{b}_t:=b(t,{\theta }_t,{\theta }_{t-\tau });{\check{b}}_t:=b(t,{\check{\theta }}_t,{\check{\theta }}_{t-\tau });\\ {\sigma }_t:=\sigma (t,{\theta }_t,{\theta }_{t-\tau });{\check{\sigma }}_t:=\sigma (t,{\check{\theta }}_t,{\check{\theta }}_{t-\tau });\\ {\partial }_x{H}_t:={\partial }_{x}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t);\\ \widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]:=\widetilde{𝔼}[{\partial }_{\mu }H(t,{\tilde{\theta }}_t,{\tilde{\theta }}_{t-\tau },{Ỹ}_t,{\tilde{Z}}_t)(X_t)]\\ 𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]:=𝔼[{\partial }_{x_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t];\\ 𝔼[\widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]|{{F}}_t]:=𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\tilde{\theta }}_{t+\tau },{\tilde{\theta }}_t,{Ỹ}_{t+\tau },{\tilde{Z}}_{t+\tau })\\ (X_t)]|{{F}}_t];\\ {\partial }_x{g}_T:={\partial }_{x}g(X_T,{\mu }_T);\widetilde{𝔼}[{\partial }_{\mu }{\tilde{g}}_T(X_T)]]\\ :=\widetilde{𝔼}[{\partial }_{\mu }g({\tilde{X}}_T,{\mu }_T)(X_T)]\\ \text{similar notation for }{\partial }_x{Ȟ}_t,\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)],𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t],\\ \text{and }𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t].& (5.9)\end{array}$

We would like to show that the map $\Phi :(\check{X},\check{Y},Ž,\check{\alpha })\to \Phi (\check{X},\check{Y},Ž,\check{\alpha })=(X,Y,Z,\alpha )$ is a contraction. Consider (ΔX, ΔY, ΔZ, Δα) = (X−X′, Y−Y′, Z−Z′, α−α′), where $(X^{\prime },Y^{\prime },Z^{\prime },{\alpha }^{\prime })=\Phi ({\check{X}}^{\prime },{\check{Y}}^{\prime },{Ž}^{\prime },{\check{\alpha }}^{\prime })$. In addition, for the following computation, we have used simplified notation:

$\begin{array}{l}\Delta b_t:=b(t,{\theta }_t,{\theta }_{t-\tau })-b(t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime });\\ \Delta {\check{b}}_t:=b(t,{\check{\theta }}_t,{\check{\theta }}_{t-\tau })-b(t,{\check{\theta }}_t^{\prime },{\check{\theta }}_{t-\tau }^{\prime });\\ \Delta {\sigma }_t:=\sigma (t,{\theta }_t,{\theta }_{t-\tau })-\sigma (t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime });\\ \Delta {\check{\sigma }}_t:=\sigma (t,{\check{\theta }}_t,{\check{\theta }}_{t-\tau })-\sigma (t,{\check{\theta }}_t^{\prime },{\check{\theta }}_{t-\tau }^{\prime })\\ {\partial }_x{g}_T:={\partial }_{x}g(X_T,{\mu }_T)-{\partial }_{x}g(X_T^{\prime },{\mu }_T);\\ \Delta \widetilde{𝔼}[{\partial }_{\mu }{\tilde{g}}_T(X_T)]]:=\widetilde{𝔼}[{\partial }_{\mu }g({\tilde{X}}_T,{\mu }_T)(X_T)]\\ -\widetilde{𝔼}[{\partial }_{\mu }g({\widetilde{\stackrel{ˇ}{X}}}_T^{\prime },{\mu }_T)(X_T^{\prime })]\\ \Delta {\partial }_x{H}_t:={\partial }_{x}H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)\\ -{\partial }_{x}H(t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime },Y_t,Z_t)\\ \Delta \widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]:=\widetilde{𝔼}[{\partial }_{\mu }H(t,{\tilde{\theta }}_t,{\tilde{\theta }}_{t-\tau },{Ỹ}_t,{\tilde{Z}}_t)(X_t)]\\ -\widetilde{𝔼}[{\partial }_{\mu }H(t,{\tilde{\theta }}_t^{\prime },{\tilde{\theta }}_{t-\tau }^{\prime },{Ỹ}_t,{\tilde{Z}}_t)(X_t^{\prime })]\\ \Delta 𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]:=𝔼[{\partial }_{x_{\tau }}H(t+\tau ,{\theta }_{t+\tau },{\theta }_t,Y_{t+\tau },Z_{t+\tau })|{{F}}_t]\\ -𝔼[{\partial }_{x_{\tau }}H(t+\tau ,{\theta }_{t+\tau }^{\prime },{\theta }_t^{\prime },Y_{t+\tau },Z_{t+\tau })|{{F}}_t]\\ \Delta 𝔼[\widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]|{{F}}_t]:=𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\tilde{\theta }}_{t+\tau },{\tilde{\theta }}_t,{Ỹ}_{t+\tau },\\ {\tilde{Z}}_{t+\tau })(X_t)]|{{F}}_t]-𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}H(t+\tau ,{\tilde{\theta }}_{t+\tau }^{\prime },{\tilde{\theta }}_t^{\prime },\\ {Ỹ}_{t+\tau },{\tilde{Z}}_{t+\tau })(X_t^{\prime })]|{{F}}_t]\\ \text{similar notation for}\Delta {\partial }_x{Ȟ}_t,\Delta \widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)],\Delta 𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t],\\ \text{and}\Delta 𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t].& (5.10)\end{array}$

Applying integration by parts to ΔX_tY_t, we have

$\begin{array}{l}d(\Delta X_t{Y}_t)\\ =Y_t\{[{\lambda }_0\Delta b_t+\delta \Delta {\check{b}}_t]dt+[{\lambda }_0\Delta {\sigma }_t+\delta \Delta {\check{\sigma }}_t]d{W}_t\}\\ -\Delta X_t\{{\lambda }_0({\partial }_x{H}_t+\widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]+𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\tilde{H}}_{t+\tau }(X_t)]|{{F}}_t])+\delta ({\partial }_x{Ȟ}_t+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]\\ +𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]+𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t])\}dt\\ +\Delta X_t{Z}_{t}d{W}_t+({\lambda }_0\Delta {\sigma }_t+\delta \Delta {\check{\sigma }}_t)Z_{t}dt.& (5.11)\end{array}$

After integrating from 0 to T, and taking expectation on both sides, we obtain

$\begin{array}{l}𝔼[\Delta X_T{Y}_T]\\ ={\lambda }_0𝔼{\displaystyle {\int }_0^T}(\Delta b_t{Y}_t+\Delta {\sigma }_t{Z}_t-\Delta X_t({\partial }_x{H}_t+\widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]\\ +𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]+𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\tilde{H}}_{t+\tau }(X_t)]|{{F}}_t]))dt\\ +\delta 𝔼{\displaystyle {\int }_0^T}(\Delta {\check{b}}_t{Y}_t+\Delta \check{\sigma }{Z}_t-\Delta X_t({\partial }_x{Ȟ}_t+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]\\ +𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]+𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t]))dt& (5.12)\end{array}$

In the meantime, from the terminal condition of Y_T given in (5.8), and since g is convex, we also have

$\begin{array}{l}𝔼[\Delta X_T{Y}_T]\\ =𝔼[\Delta X_T({\lambda }_0({\partial }_x{g}_T+\widetilde{𝔼}[{\partial }_{\mu }{\tilde{g}}_T(X_T)])+\delta ({\partial }_x{ǧ}_T\\ +\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])+I_T^g)]\ge {\lambda }_0𝔼[g(X_T,{\mu }_T)-g(X_T^{\prime }-{\mu }_T^{{}^{\prime }})]\\ +\delta \Delta X_T({\partial }_x{ǧ}_T+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])+\Delta X_T{I}_T^g& (5.13)\end{array}$

Following the proof of sufficient part of maximum principle and using (5.12), and (5.13), we find

$\begin{array}{l}{\lambda }_0(J(\alpha )-J({\alpha }^{\prime }))\\ ={\lambda }_0𝔼[g(X_T,{\mu }_T)-g(X_T^{\prime },{\mu }_T^{\prime })]\\ +{\lambda }_0𝔼{\displaystyle {\int }_0^T}[f(t,{\theta }_t,X_{t-\tau },{\mu }_{t-\tau })\\ -f(t,{\theta }_t^{\prime },X_{t-\tau }^{\prime },{\mu }_{t-\tau }^{\prime })]dt\\ \le 𝔼[\Delta X_T{Y}_T]-\delta \Delta X_T({\partial }_x{ǧ}_T+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])-\Delta X_T{I}_T^g\\ +{\lambda }_0𝔼{\displaystyle {\int }_0^T}[f(t,{\theta }_t,X_{t-\tau },{\mu }_{t-\tau })-f(t,{\theta }_t^{\prime },X_{t-\tau }^{\prime },{\mu }_{t-\tau }^{\prime })]dt\\ ={\lambda }_0𝔼{\displaystyle {\int }_0^T}[\Delta b_t{Y}_t+\Delta {\sigma }_t{Z}_t\\ -\Delta X_t({\partial }_x{H}_t+\widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]+𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\tilde{H}}_{t+\tau }(X_t)]|{{F}}_t])]dt\\ +\delta 𝔼{\displaystyle {\int }_0^T}[\Delta {\check{b}}_t{Y}_t+\Delta {\check{\sigma }}_t{Z}_t-\Delta X_t({\partial }_x{Ȟ}_t+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]\\ +𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]+𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }(X_t)]|{{F}}_t])]dt\\ +{\lambda }_0𝔼{\displaystyle {\int }_0^T}[H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)-H(t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime },Y_t,Z_t)]dt\\ -{\lambda }_0𝔼{\displaystyle {\int }_0^T}(\Delta b_t{Y}_t+{\Delta }_t\sigma Z_t)dt-\delta \Delta X_T({\partial }_x{ǧ}_T\\ +\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])-\Delta X_T{I}_T^g\\ ={\lambda }_0𝔼{\displaystyle {\int }_0^T}[H(t,{\theta }_t,{\theta }_{t-\tau },Y_t,Z_t)-H(t,{\theta }_t^{\prime },{\theta }_{t-\tau }^{\prime },Y_t,Z_t)\\ -\Delta X_t({\partial }_x{H}_t+\widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]+𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\tilde{H}}_{t+\tau }(X_t)]|{{F}}_t])]dt+\delta 𝔼{\displaystyle {\int }_0^T}[\Delta {\check{b}}_t{Y}_t\\ +\Delta \check{\sigma }{Z}_t-\Delta X_t({\partial }_x{Ȟ}_t+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]+𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t])]dt\\ -\delta \Delta X_T({\partial }_x{ǧ}_T+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])-\Delta X_T{I}_T^g\\ \le -𝔼{\displaystyle {\int }_0^T}{\lambda }_0\kappa |\Delta {\alpha }_t{|}^{2}dt+\delta 𝔼{\displaystyle {\int }_0^T}[\Delta {\check{b}}_t{Y}_t+\Delta {\check{\sigma }}_t{Z}_t\\ -\Delta X_t({\partial }_x{Ȟ}_t+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]+𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t])]dt-\delta \Delta X_T({\partial }_x{ǧ}_T\\ +\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])-\Delta X_T{I}_T^g& (5.14)\end{array}$

Reverse the role of α and α′, we also have

$\begin{array}{l}{\lambda }_0(J({\alpha }^{\prime })-J(\alpha ))\\ \le -𝔼{\displaystyle {\int }_0^T}{\lambda }_0\kappa |\Delta {\alpha }_t^{\prime }{|}^{2}dt+\delta 𝔼{\displaystyle {\int }_0^T}[\Delta {\check{b}}_t^{\prime }{Y}_t^{\prime }+\Delta {\check{\sigma }}_t^{\prime }{Z}_t^{\prime }\\ -\Delta X_t^{\prime }({\partial }_x{Ȟ}_t^{\prime }+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t^{\prime }({\check{X}}_t^{\prime })]+𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }^{\prime }|{{F}}_t]\\ +𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }^{\prime }({\check{X}}_t^{\prime })]|{{F}}_t])]dt-\delta \Delta X_T^{\prime }({\partial }_x{ǧ}_T\\ +\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T^{\prime })])-\Delta X_T^{\prime }{I}_T^g& (5.15)\end{array}$

Summing (5.14) and (5.15), using the fact that b and σ have the linear form, using change of time and Lipschitz assumption, it yields

$\begin{array}{l}2{\lambda }_0\kappa 𝔼{\displaystyle {\int }_0^T}|\Delta {\alpha }_t{|}^{2}dt\\ \le \delta 𝔼{\displaystyle {\int }_0^T}[\Delta {\check{b}}_t\Delta Y_t+\Delta \check{\sigma }\Delta Z_t-\Delta X_t(\Delta {\partial }_x{Ȟ}_t\\ +\Delta \widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t(X_t)]+\Delta 𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]\\ +\Delta 𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }(X_t)]|{{F}}_t])]dt+\delta \Delta X_T({\partial }_{x^{\prime }}{ǧ}_T^{\prime }\\ -{\partial }_x{ǧ}_T+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T^{\prime }({\check{X}}_T)]-\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)])\\ \le \frac{1}{2}𝔼{\displaystyle {\int }_0^T}[ϵ(|\Delta X_t{|}^2+|\Delta Y_t{|}^2+|\Delta Z_t{|}^2)\\ +\frac{1}{ϵ}{\delta }^2(|\Delta {\check{b}}_t{|}^2+|\Delta \check{\sigma }{|}^2+|\Delta {\partial }_x{Ȟ}_t+\Delta \widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t(X_t)]\\ +\Delta 𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]+\Delta 𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }(X_t)]|{{F}}_t]{|}^2)]dt\\ +\frac{1}{2}(ϵ|\Delta X_T{|}^2+\frac{1}{ϵ}{\delta }^2|{\partial }_{x^{\prime }}{ǧ}_T^{\prime }-{\partial }_x{ǧ}_T+\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T^{\prime }({\check{X}}_T)]\\ -\widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)]{|}^2)\le \frac{1}{2}ϵ𝔼[{\displaystyle {\int }_0^T}ϵ(|\Delta X_t{|}^2+|\Delta Y_t{|}^2\\ +|\Delta Z_t{|}^2+|\Delta {\alpha }_t{|}^2)dt+|\Delta X_T{|}^2]\\ +\frac{1}{2}\delta \frac{C}{ϵ}𝔼[{\displaystyle {\int }_0^T}(|\Delta {\check{X}}_t{|}^2+|\Delta {\check{Y}}_t{|}^2+|\Delta {Ž}_t{|}^2\\ +|\Delta {\check{\alpha }}_t{|}^2)]dt+|\Delta {\check{X}}_T{|}^2],& (5.16)\end{array}$

Next, we apply Ito's formula to $\Delta X_t^2$,

$\begin{array}{l}d\Delta X_t^2\\ =2\Delta X_{t}d{X}_t+d{\langle X,X\rangle }_t\\ =2\Delta X_t({\lambda }_0\Delta b_t+\delta \Delta {\check{b}}_t)dt+2\Delta X_t({\lambda }_0\Delta {\sigma }_t\\ +\delta \Delta {\check{\sigma }}_t)d{W}_t+{({\lambda }_0\Delta {\sigma }_t+\delta \Delta {\check{\sigma }}_t)}^{2}dt& (5.17)\end{array}$

Then integrate from 0 to T, and take expectation,

$\begin{array}{l}𝔼[|\Delta X_t{|}^2]\\ =2{\lambda }_0𝔼{\displaystyle {\int }_0^t}|\Delta X_s\Delta b_s|ds+2\delta 𝔼{\displaystyle {\int }_0^t}|\Delta X_s\Delta {\check{b}}_s|ds\\ +𝔼{\displaystyle {\int }_0^t}|{\lambda }_0\Delta {\sigma }_s+\delta \Delta {\check{\sigma }}_s{|}^{2}ds\\ \le {\lambda }_0𝔼{\displaystyle {\int }_0^t}(|\Delta X_s{|}^2+|\Delta b_s{|}^2)ds+𝔼{\displaystyle {\int }_0^t}(|\Delta X_s{|}^2+{\delta }^2|\Delta {\check{b}}_s{|}^2)ds\\ +𝔼{\displaystyle {\int }_0^t}(2{\lambda }_0^2|\Delta {\sigma }_s{|}^2+2{\delta }^2|\Delta {\check{\sigma }}_s{|}^2)ds\\ \le C𝔼{\displaystyle {\int }_0^{t+\tau }}(|\Delta X_s{|}^2+|\Delta {\alpha }_s{|}^2)ds+\delta C𝔼{\displaystyle {\int }_0^{t+\tau }}(|\Delta {\check{X}}_s{|}^2+|\Delta {\check{\alpha }}_s{|}^2)ds& (5.18)\end{array}$

From Gronwall's inequality, we can obtain

$\begin{array}{ll}{\displaystyle \underset{0\le t\le T}{sup}}𝔼[|X_t{|}^2]\le C𝔼{\displaystyle {\int }_0^T}|\Delta {\alpha }_t{|}^{2}dt+\delta C𝔼{\displaystyle {\int }_0^T}(|\Delta {\check{X}}_t{|}^2+|\Delta {\check{\alpha }}_t{|}^2)dt& (5.19)\end{array}$

Similarly, applying Ito's formula to $|\Delta Y_t{|}^2$, and taking expectation, we have

$\begin{array}{l}𝔼\left[|\Delta Y_t{|}^2+{\displaystyle {\int }_t^T}|\Delta Z_s{|}^{2}ds\right]\\ =2{\lambda }_0𝔼{\displaystyle {\int }_t^T}|\Delta Y_t(\Delta {\partial }_x{H}_t+\Delta \widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]\\ +\Delta 𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]+\Delta 𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\tilde{H}}_{t+\tau }(X_t)]|{{F}}_t])|\\ +2\delta 𝔼{\displaystyle {\int }_t^T}|\Delta Y_t(\Delta {\partial }_x{Ȟ}_t+\Delta \widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]+\Delta 𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]\\ +\Delta 𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t])|+𝔼|\Delta Y_T{|}^2\\ \le 𝔼{\displaystyle {\int }_t^T}(\frac{1}{ϵ}|\Delta Y_t{|}^2+ϵ{\lambda }_0^2|\Delta {\partial }_x{H}_t+\Delta \widetilde{𝔼}[{\partial }_{\mu }{\tilde{H}}_t(X_t)]\\ +\Delta 𝔼[{\partial }_{x_{\tau }}{H}_{t+\tau }|{{F}}_t]+\Delta 𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\tilde{H}}_{t+\tau }(X_t)]|{{F}}_t]{|}^2)dt\\ +𝔼{\displaystyle {\int }_t^T}(|\Delta Y_t{|}^2+{\delta }^2|\Delta {\partial }_x{Ȟ}_t+\Delta \widetilde{𝔼}[{\partial }_{\mu }{\widetilde{Ȟ}}_t({\check{X}}_t)]\\ +\Delta 𝔼[{\partial }_{x_{\tau }}{Ȟ}_{t+\tau }|{{F}}_t]+\Delta 𝔼[\widetilde{𝔼}[{\partial }_{{\mu }_{\tau }}{\widetilde{Ȟ}}_{t+\tau }({\check{X}}_t)]|{{F}}_t]{|}^2)dt\\ +𝔼|{\lambda }_0(\Delta {\partial }_x{g}_T+\Delta \widetilde{𝔼}[{\partial }_{\mu }{\tilde{g}}_T(X_T)])\\ +\delta (\Delta {\partial }_x{ǧ}_T+\Delta \widetilde{𝔼}[{\partial }_{\mu }{\widetilde{ǧ}}_T({\check{X}}_T)]){|}^2& (5.20)\end{array}$

Choose ϵ = 96 max{R², L}, and from assumption (H5.1 - H5.2) and Gronwall's inequality, we obtain a bound for ${sup}_{0\le t\le T}𝔼|\Delta Y_t{|}^2$; and then substitute the it back to the same inequality, we are able to obtain the bound for ${\int }_0^{T}E|Z_t{|}^{2}dt$. By combining these two bounds, we deduce that

$\begin{array}{l}𝔼\left[{\displaystyle \underset{0\le t\le T}{sup}}|Y_t{|}^2+{\displaystyle {\int }_0^T}|Z_t{|}^{2}dt\right]\\ \le C𝔼\left({\displaystyle \underset{0\le t\le T}{sup}}|\Delta X_t{|}^2+{\displaystyle {\int }_0^T}|\Delta {\alpha }_t{|}^{2}dt\right)\\ +\delta C𝔼[{\displaystyle \underset{0\le t\le T}{sup}}\left(|\Delta {\check{X}}_t{|}^2+|\Delta {\check{Y}}_t{|}^2\right)\\ +{\displaystyle {\int }_0^T}\left(|\Delta {Ž}_t{|}^2+|\Delta {\check{\alpha }}_t{|}^2\right)dt]& (5.21)\end{array}$

Finally, combining (5.19) and (5.21), and (5.16), we deduce

$\begin{array}{l}𝔼\left[{\displaystyle \underset{0\le t\le T}{sup}}|\Delta X_t{|}^2+{\displaystyle \underset{0\le t\le T}{sup}}|\Delta Y_t{|}^2+{\displaystyle {\int }_0^T}(|\Delta Z_t^2|+|\Delta {\alpha }_t{|}^2)dt\right]\\ \le \delta C𝔼\left[{\displaystyle \underset{0\le t\le T}{sup}}|\Delta {\check{X}}_t{|}^2+{\displaystyle \underset{0\le t\le T}{sup}}|\Delta {\check{Y}}_t{|}^2+{\displaystyle {\int }_0^T}\left(|\Delta {Ž}_t{|}^2+|\Delta {\check{\alpha }}_t{|}^2\right)dt\right]& (5.22)\end{array}$

Let ${\delta }_0=\frac{1}{2C}$, it is clear that the mapping Φ is a contraction for all δ ∈ (0, δ₀). It follows that there is a unique fixed point which is the solution of MV-FABSDE(λ) for λ = λ₀ + δ, δ ∈ (0, δ₀).□

Proof of Theorem 5.1: For λ = 0, FABSDE(0) has a unique solution. Using Lemma 5.2, there exists a δ₀ > 0 such that FBSDE(δ) has a unique solution for δ ∈ [0, δ₀], assuming (n − 1)δ₀ < 1 ≤ nδ₀. Following by a induction argument, we repeat Lemma 5.2 for n times, which gives us the existence of the unique solution of FABSDE(1).□

6. Conclusion

Overall, we presented a comprehensive study of a general class of mean-field control problems with delay effect. The state dynamics is described by a McKean–Vlasov stochastic delayed differential equation. We derive the adjoint process associated to the dynamics, which is in the form of an anticipated backward stochastic differential equation of McKean–Vlasov type. We also prove a version of stochastic maximum principle, which gives necessary and sufficient conditions for the optimal control. Furthermore, we prove the existence and uniqueness of this class of forward anticipated backward stochastic differential equations under suitable conditions.

However, due to the lack of explicit solutions, numerical methods are needed. The non-linear nature of the problem due to the McKean–Vlasov aspect combined with non-Markovianity due to delay rule out classical numerical methods. Our study show that deep learning methods can deal with these obstacles. we proposed two algorithms based on machine learning to numerically solve the control problem. One is to directly approximate the control using a neural network, while the loss is given by the objective function in the control problem. The other algorithm is based on the system of forward and backward stochastic differential equations. We approximate the adjoint processes (Y_·, Z_·) and the conditional expectation of the adjoint process $𝔼[Y_{·+\tau }|{{F}}_{·}]$ using neural networks. In this case, there are two loss functions that we need to minimize as shown in (3.16). The first loss is associated with the adjoint process Y_·, and the other one is related to $𝔼[Y_{·+\tau }|{{F}}_{·}]$. After minimizing the losses, the optimal control can be readily obtained from the adjoint processes. In addition, Figure 8 illustrates that the optimal value to the control problem, that is computed from different methods, converges. Moreover, our method also works when the control problem has no delay. As a sanity check, we provide an example of control problem without delay effect, and we show that the solution obtained from the algorithm matches the explicit solution.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This work supported by NSF grant DMS-1814091.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Carmona R, Fouque JP, Sun LH. Mean field games and systemic risk. Commun Math Sci. (2015) 13:911–33. doi: 10.4310/CMS.2015.v13.n4.a4

CrossRef Full Text | Google Scholar

2. Carmona R, Fouque JP, Mousavi SM, Sun LH. Systemic risk and stochastic games with delay. J Optim Appl. (2018) 179:366–99. doi: 10.1007/s10957-018-1267-8

CrossRef Full Text | Google Scholar

3. Carmona R, Delarue F. Probabilistic Theory of Mean Field Games with Applications I & II. Springer International Publishing (2018).

Google Scholar

4. Bensoussan A, Frehse J, Yam SCP. Mean Field Games and Mean Field Type Control Theory. New York, NY: Springer-Verlag (2013).

Google Scholar

5. Peng S, Yang Z. Anticipated backward stochastic differential equations. Ann Probab. (2009) 37:877–902. doi: 10.1214/08-AOP423

CrossRef Full Text | Google Scholar

6. Zhang J. Backward Stochastic Differential Equations. New York, NY: Springer-Verlag (2017).

Google Scholar

7. Chen L, Wu Z. Maximum principle for the stochastic optimal control problem with delay and application. Automatica. (2010) 46:1074–80. doi: 10.1016/j.automatica.2010.03.005

CrossRef Full Text | Google Scholar

8. Peng S, Wu Z. Fully coupled forward-backward stochastic differential equations and applications to optimal control. SIAM J Control Optim. (1999) 37:825–43. doi: 10.1137/S0363012996313549

CrossRef Full Text | Google Scholar

9. Bensoussan A, Yam SCP, Zhang Z. Well-posedness of mean-field type forward-backward stochastic differential equations. Stochastic Process Appl. (2015) 125:3327–54. doi: 10.1016/j.spa.2015.04.006

CrossRef Full Text | Google Scholar

10. Ma J, Yong J. Forward-Backward Stochastic Differential Equations and their Applications. Berlin; Heidelberg: Springer-Verlag (2007).

Google Scholar

11. Ma J, Protter P, Yong J. Solving forward-backward stochastic differential equations explicitly–a four step scheme. Probab Theory Relat Fields. (1994) 98:339–59.

Google Scholar

12. Weinan E, Han J, Jentzen A. Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations. Commun Math Stat. (2017) 5:349–80. doi: 10.1007/s40304-017-0117-6

CrossRef Full Text | Google Scholar

13. Raissi M. Forward-backward stochastic neural networks: deep learning of high-dimensional partial differential equations. (2018). arxiv[Preprint].arxiv:1804.07010.

Google Scholar

14. Hochreiter S, Schmidhuber J. Long short-term memory. Neural Comput. (1997) 9:1735–80.

PubMed Abstract | Google Scholar

15. Sutton RS, McAllester DA, Singh SP, Mansour Y. Policy gradient methods for reinforcement learning with function approximation. In: Advances in Neural Information Processing Systems (NIPS). (1999). p. 1057–63. Available online at: https://papers.nips.cc/paper/1713-policy-gradient-methods-for-reinforcement-learning-with-function-approximation

Google Scholar

16. Han J, Weinan E. Deep Learning Approximation for Stochastic Control Problems, Deep Reinforcement Learning Workshop. Conference on Neural Information Processing Systems (NIPS) (2016).

17. Carmona R, Delarue F, Lachapelle A. Control of McKean-Vlasov dynamics versus mean field games. Math Fin Econ. (2013) 7:131–66. doi: 10.1007/s11579-012-0089-y

CrossRef Full Text | Google Scholar

18. Fouque JP, Zhang Z. Mean field game with delay: a toy model. Risks. (2018) 6:1–17. doi: 10.3390/risks6030090

CrossRef Full Text | Google Scholar