Olivier Guéant
email: olivier.gueant@ann.jussieu.fr

Charles-Albert Lehalle
email: clehalle@cheuvreux.com

Joaquin Fernandez Tapia

Bruno Bouchard

Optimal Portfolio Liquidation with Limit Orders

, we link the optimal trade-schedule to the price of the limit orders that have to be sent to the limit order book to optimally liquidate a portfolio. Most practitioners address these two issues separately: they compute an optimal trading curve and they then send orders to the markets to try to follow it. The results obtained here solve simultaneously the two problems. As in a previous paper that solved the "intra-day market making problem" [19], the interactions of limit orders with the market are modeled via a Poisson process pegged to a diffusive "fair price" and a Hamilton-Jacobi-Bellman equation is used to solve the trade-off between execution risk and price risk. Backtests are finally carried out to exemplify the use of our results.

Introduction

Optimal scheduling of large orders in order to control the overall trading costs with a tradeoff between market impact (demanding to trade slow) and market risk (urging to trade fast) has been proposed in the litterature in the late nineties mainly by Bertsimas and Lo [START_REF] Bertsimas | Optimal control of execution costs[END_REF] and Almgren and Chriss [START_REF] Almgren | Optimal execution of portfolio transactions[END_REF]. If the original approach has been recently generalized in several directions (see for instance [START_REF] Almgren | Optimal Trading in a Dynamic Market[END_REF][START_REF] Almgren | Optimal execution with nonlinear impact functions and tradingenhanced risk[END_REF][START_REF] Dang | Optimal trading with transient price impact. A comparison of discrete and continuous approach[END_REF][START_REF] Forsyth | Optimal Trade Execution: A Mean-Quadratic-Variation Approach[END_REF][START_REF] Peter | A Hamilton Jacobi Bellman approach to optimal trade execution[END_REF][START_REF] Gatheral | Transient Linear Price Impact and Fredholm Integral Equations[END_REF][START_REF] Lehalle | Rigorous Strategic Trading: Balanced Portfolio and Mean-Reversion[END_REF][START_REF] Lorenz | Mean-variance optimal adaptive execution[END_REF][START_REF] Schied | Risk aversion and the dynamics of optimal liquidation strategies in illiquid markets[END_REF]), only few attempts have been made to drill down the model at the level of the interactions with the order books. The more noticeable proposal is the one by Obizhaeva and Wang [START_REF] Obizhaeva | Optimal Trading Strategy and Supply/Demand Dynamics[END_REF], followed and sophisticated by Alfonsi, Fruth and Schied [START_REF] Alfonsi | Optimal execution strategies in limit order books with general shape functions[END_REF] or Predoiu, Shaikhet and Shreve [START_REF] Predoiu | Optimal Execution of a General One-Sided Limit-Order Book[END_REF]. This branch of the optimal trading literature 1 focuses on the dynamics initiated by aggressive orders hitting a martingale and resilient order-book 2 , ignoring then trading by passive orders.

Just recall that during the continuous auction processes implemented by most electronic trading pools, market participants send their open interests (i.e. buy or sell orders) to a queuing system where a "first in first out" queue stands at each possible price. If a buy (respectively sell) order reaches a queue of sell (resp. buy) orders, a transaction occurs (see for instance [START_REF] Lehalle | High Frequency Simulations of an Order Book: a Two-Scales Approach[END_REF] for more explanations and modelings). Orders generating trades are said to be aggressive or liquidity-consuming ones; orders filling queues are said to be passive or liquidity-providing ones. In practice, most trading algorithms are as passive as possible (a typical balance for a scheduling algorithm is around 60% of passively obtained trades -see [START_REF] Lehalle | The Established Liquidity Fragmentation Affects all Investors[END_REF]). The economic literature first explored and studied these interactions between orders sent to a continuous auction system by different actors from a global efficiency viewpoint 3 . With the fragmentation of equity markets in the US and in Europe, the issue of linking optimal posting prices to the optimal liquidation of a portfolio is more and more important. A trading algorithm has to find an optimal scheduling or rhythm for its trading, but also to choose a sequel of prices and quantities of orders to send to the markets to follow this optimal rhythm as much as possible.

This paper answers this "optimal scheduling and posting" problem as a whole. Indeed, and in contrast with most of the preceding literature 4 , we propose a new approach which is liquidity-providing oriented: liquidation strategies involve limit orders and not only market orders. As a consequence, the classical trade-off of the literature between market impact, or execution costs, and price risk disappears in our setting since no execution cost is incurred. However, since the broker does not know when his orders are going to be executed -if at all -, a new risk is borne: (non-)execution risk. If a limit ask order is inserted in the order book, probability of execution and eventually time of execution will depend on the price of the order. In our framework, the flow of trades "hitting" a passive order at a distance δ from the "fair price" S t (modeled by a Brownian Motion) follows an adapted Poisson process of intensity A exp(-kδ). It means that the further away from the "fair price" an order is posted, the less transactions it will obtain. If the limit order price is far above the best ask price, the trading gain may be high but execution is not guaranteed and the broker is exposed to the risk of a price decrease. On the contrary, if the limit order price is near the best ask price, or even reduces the market bid-ask spread, gains will be small but the probability of execution will be higher, resulting in faster trading and less price risk.

Our modeling paradigm for optimal liquidation is in fact rooted to recent works in other areas of algorithmic trading. Market making models developed by Ho and Stoll [START_REF] Ho | Optimal dealer pricing under transactions and return uncertainty[END_REF] or more recently by Avellaneda and Stoikov [START_REF] Avellaneda | High-frequency trading in a limit order book[END_REF] for "high frequency market making in an order-book" are instances of the use of limit orders in the financial literature. Our approach is more precisely inspired from the resolution of the "market making" problem in a companion paper [START_REF] Guéant | Dealing with the inventory risk[END_REF].

The main limitations of our models are twofold. First, our framework deals with one trading venue only. Hence we do not model a "smart order routing" (SOR) mechanism but it can be seen as the consolidation of all available trading venues, without taking directly into account the potential specificity of each of them5 . Second, market impact is not modeled in our framework, be it permanent or transient. In previous works on optimal trade scheduling, market impact was typically rendered by an explicit model because trading was done with liquidity-consuming orders. The same kind of assumptions could have been made on executed orders but, the introduction of limit orders in the literature on optimal trading being quite recent, there is still no model for the market impact of liquidity-providing orders.

To the authors' knowledge, it is indeed the first proposal to drill down to passive orders modeling to solve the optimal trade scheduling for large orders 6 . It is organized as follows: in the first part, we present the setting of the model and the main hypotheses on execution. The second part is devoted to the resolution of the partial differential equations arising from the control problem. Part 3, deals with three special cases, namely the time-asymptotic case, the "no-volatility" benchmark and a limiting case in which the trader has a large incentive to liquidate before the end. These special cases provide closed-form formulae allowing us to better understand the forces at stake. Then, in part 4, we carry out comparative statics and discuss the way optimal strategies depend on the model parameters. Finally, in part 5, we show how our approach can be used in practice for optimal liquidation.

Setup of the model

We consider a trader who has to liquidate a portfolio containing a large quantity q 0 of a given stock. We suppose that the reference price of the stock (that can be considered the mid-price or the best bid quote for example) moves as a brownian motion with a drift:

dS t = µdt + σdW t
The trader under consideration will continuously propose an ask quote denoted S a t = S t + δ a t and will hence sell shares according to the rate of arrival of aggressive orders at the prices he quotes. His inventory q, that is the quantity he holds, is given by q t = q 0 -N a t where N a is the jump process giving the number of shares he sold 7 . This jump process is supposed to be a Poisson process and to simplify the exposition we consider that jumps are of unitary size 8 . Arrival rates obviously depend on the price S a t quoted by the trader and we assume, in accordance with most datasets, that intensity λ a associated to N a is of the following form:

λ a (δ a ) = A exp(-kδ a ) = A exp(-k(s a -s))
This means that the cheaper the order price, the faster it will be executed.

As a consequence of his trades, the trader has an amount of cash whose dynamics is given by:

dX t = (S t + δ a t )dN a t
The trader has a time horizon T to liquidate the portfolio and his goal is to optimize the expected utility of his P&L at time T . We will focus on CARA utility functions and we suppose that the trader optimizes: sup

(δ a t )t∈A E [-exp (-γX T + q T (S T -b))]
where A is the set of predictable processes on [0,T], bounded from below, where γ is the absolute risk aversion characterizing the trader, where X T is the amount of cash at time T , where q T is the remaining quantity of shares in the inventory at time T and where b is a cost one has to incur to liquidate the eventual remaining quantity.

Resolution

Hamilton-Jacobi-Bellman equation

The optimization problem set up in the preceding section can be solved using classical Bellman tools. To this purpose, we introduce the Hamilton-Jacobi-Bellman equation associated to the optimization problem:

(HJB) ∂ t u(t, x, q, s) + µ∂ s u(t, x, q, s) + 1 2 σ 2 ∂ 2 ss u(t, x, q, s)

+ sup δ a λ a (δ a ) [u(t, x + s + δ a , q -1, s) -u(t, x, q, s)] = 0
with the final condition:

u(T, x, q, s) = -exp (-γ(x + q(s -b)))
and the boundary condition:

u(t, x, 0, s) = -exp (-γx)
To solve the Hamilton-Jacobi-Belmann equation, we will use a change of variables that transforms the PDEs in a system of linear ODEs. Proposition 1 (A system of linear ODEs). Let's consider a family of functions (w q ) q∈N solution of the linear system of ODEs (S) that follows: ∀q ∈ N, ẇq (t) = (αq 2 -βq)w q (t) -ηw q-1 (t) with w q (T ) = e -kqb and w 0 = 1, where α = k 2 γσ 2 , β = kµ and η = A(1

+ γ k ) -(1+ k γ ) .
Then u(t, x, q, s) = -exp(-γ(x + qs))w q (t) -γ k is solution of (HJB). Now, using this system of ODEs, we can find the optimal quotes through a verification theorem:

Theorem 1 (Verification theorem and optimal quotes). Let's consider the solution w of the system (S) of Proposition 1. Then, u(t, x, q, s) = -exp(-γ(x + qs))w q (t) -γ k is the value function of the control problem and the optimal ask quote can be expressed as:

δ a * (t, q) = ( 1 k ln ( w q (t) w q-1 (t) ) + 1 γ ln ( 1 + γ k ) )

Numerical example

Proposition 1 and Theorem 1 provide a way to solve the Hamilton-Jacobi-Bellman equation and to derive the optimal quotes for a trader willing to liquidate a portfolio. To exemplify these results, considering a very large value 9 for b, we compute the optimal quotes when a quantity of shares up to 6 times the ATS has to be sold within 5 minutes (Figure 1). Optimal strategy δ a * (t, q) (in Ticks) for an agent willing to sell a quantity of shares up to 6 times the ATS within 5 minutes (µ = 0 (Tick.s -1 ), σ = 0.3 (Tick.s -1 2 ), A = 0.1 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

We clearly see that the optimal quotes depend on time and inventory in a monotonic way. Indeed, a trader with a lot of shares to get rid of wants to trade fast and will therefore propose a low price. On the contrary a trader with only a few shares in his portfolio may be willing to benefit from a trading opportunity and will send an order with a higher price because the risk he bears allows him to trade more slowly. Now, coming to the time-dependence of the quotes, a trader with a given number of shares will, ceteris paribus, lower his quotes as the time horizon gets closer. At the limit when t is close to the time horizon T , the optimal quotes tends to very low values as we considered very large b. In practice very negative values for the ask quotes have to be understood as market orders on the bid side. Also, we see that optimal quotes have an asymptotic behavior as time horizon increases. The associated limiting case will be studied in depth in the next section.

Finally, the trading curve associated to the above optimal quotes can be obtained by Monte-Carlo simulations as exemplified on Figure 2.

-1 ), σ = 0.3 (Tick.s -1 2 ), A = 0.1 (s -1 ), k = 0.3 (Tick -1
) and γ = 0.05 (Tick -1 )). Dotted line: Trading curve calculated using the Almgren-Chriss framework with the same parameters and the impact function

f (v) = 20v 1+ 3 4 .

Special cases

The above equations can be solved explicitly for w and hence for the optimal quotes using the above verification theorem. However, the resulting closed-form expressions are not really tractable and do not provide any intuition on the behavior of the optimal quotes. Three special cases are now considered for which simpler closed-form formulae can be derived. We start first with the limiting behavior of the quotes, when the time horizon T tends to infinity. We then consider a benchmark case in which there is no volatility (σ = 0) as an approximation for low-volatility cases and we finally consider, by analogy with the classical literature, the behavior of the solution as the liquidation cost b increases.

Asymptotic behavior as T → +∞

We have seen on Figure 1 that the optimal quotes seem to exhibit an asymptotic behavior.

We are in fact going to prove that δ a * (0, q) tends to a limit as the time horizon T tends to infinity.

Proposition 2 (Asymptotic behavior of the optimal quotes). Let's suppose that 10 µ γσ 2 < 1 2 . Let's consider the solution w of the system (S) of Proposition 1. Then:

lim T →+∞ w q (0) = η q q! q ∏ j=1 1 αj -β
The resulting asymptotic behavior for the optimal ask quote of Theorem 1 is:

lim T →+∞ δ a * (0, q) = 1 k ln ( A 1 + γ k 1 αq 2 -βq
) 10 This condition is the same as α > β.

The "no-volatility" benchmark

Though unrealistic, we now concentrate on a benchmark case in which there is no volatility (σ = 0). In that benchmark case, the trader bears no price risk and the only risk he faces is linked to the non-execution of his orders.

In this framework, we can derive tractable formulae for w and for the optimal quotes. Proposition 3 (The "no-volatility" benchmark). Assume that σ = 0 and that there is no drift (µ = 0). Let's define:

w q (t) = q ∑ j=0 η j j! e -kb(q-j) (T -t) j
Then w defines a solution of the system (S) and the optimal quote is:

δ a * (t, q) =   1 k ln   ∑ q j=0 η j j! e -kb(q-j) (T -t) j ∑ q-1 j=0 η j j! e -kb(q-1-j) (T -t) j   + 1 γ ln ( 1 + γ k )  
It is noteworthy that, in the b = 0 case corresponding to a questionable absence of liquidation cost, we can find the results of [START_REF] Bayraktar | Liquidation in Limit Order Books with Controlled Intensity[END_REF] as a special case of ours when γ tends to 0.

Limiting behavior as b → +∞

Let us now consider the limiting case b → +∞. Sending b to infinity corresponds to a situation in which more and more incentives are given to the trader for complete liquidation before time T . If we look at the Almgren-Chriss-like literature on optimal execution, the authors are often assuming that q T = 0 11 . Hence, if one writes the value functions associated to most liquidity-consuming optimal strategies, it turns out that they are equal to -∞ at the horizon time T except when the inventory is equal to nought (hence b = +∞, in this paper). However, here, due to the uncertainty on execution, we cannot write a well-defined control problem when b is equal to +∞. Rather, we are interested in the limiting behavior when b → +∞ and the limit we obtain provide approximations for when the incentive to liquidate before time T is large. By analogy with the initial literature on optimal liquidation [START_REF] Almgren | Optimal execution of portfolio transactions[END_REF], we can also have some limiting results on the so-called trading curve. Hereafter we denote w b,q (t) the solution of the system (S) for a given liquidation cost b, δ a * b (t, q) the associated optimal quote and q b,t the resulting process modeling the number of stocks in the portfolio.

Proposition 4 (Form of the solutions and trading intensity). The limiting solution lim b→+∞ w b,q (t) is of the form A q v q (t) where v does not depend on A. The limit of the trading intensity lim b→+∞ Ae -kδ a * b does not depend on A. If we define the trading curve as the average12 evolution of the inventory, i.e. V b (t) := E[q b,t ], then the limit of the trading curve V (t) = lim b→+∞ V b (t) does not depend on A.

More results can be obtained in the no-volatility case: Proposition 5 (The "no-volatility" benchmark, b → +∞). Assume that σ = 0 and consider first the no-drift case (µ = 0).

lim b→+∞ w b,q (t) = η q q! (T -t) q
The limit of the optimal quote is given by:

lim b→+∞ δ a * b (t, q) = ( 1 k ln ( η q (T -t) ) + 1 γ ln ( 1 + γ k ) )
that can also be written:

lim b→+∞ δ a * b (t, q) = 1 k ln ( A 1 + γ k 1 q (T -t)
)

The limit of the associated trading curve is

V (t) = q 0 ( 1 -t T ) 1+ γ k .
If µ ̸ = 0, similar results can be obtained:

lim b→+∞ w b,q (t) = η q q! ( e β(T -t) -1 β ) q
The limit of the optimal quote is:

lim b→+∞ δ a * b (t, q) = ( 1 k ln ( η q e β(T -t) -1 β ) + 1 γ ln ( 1 + γ k ) )
that can also be written:

lim b→+∞ δ a * b (t, q) = 1 k ln ( A 1 + γ k 1 q e β(T -t) -1 β )
The limit of the associated trading curve is

V (t) = q 0 ( 1-e -β(T -t) 1-e -βT ) 1+ γ k .
4 Comparative statics

Intuition from the above cases

Before going to comparative statics numerically, we focus on the particular cases for which tractable closed form formulae have just been derived. These closed-form formulae provide an approximation for the optimal quotes and trading curves respectively when t is far from T , when volatility is small or when liquidation cost is large. Hence, we can start to discuss the role played by the different parameters in these contexts.

Focusing first on the optimal quotes, we see both in the asymptotic case and in the "novolatility" benchmark case (when b → ∞) that, quite naturally, the optimal ask quote is an increasing function of A. When trading intensity is high, many trades occur and a limit order inserted far above the reference price has indeed a large probability to be executed. Risk aversion also plays an important role. Both in the asymptotic case and in the "novolatility" benchmark case (also this last one misses part of the picture, ignoring price risk), a very risk-adverse agent will indeed be willing to reduce execution risk and will submit orders at low price. Now, as far as k is concerned, the dependence of the optimal quote on k is ambiguous because the interpretation of k depends on the optimal quote itself. An increase in k corresponds indeed to a decrease in the probability to be executed for prices above the reference price. However, due to the exponential form for execution intensity, the exact opposite is true for quotes below reference price, an increase in k implying an increase in the probability to be executed for quotes below reference price. Now, turning to the trading curve in the "no-volatility" benchmark case when b → ∞, we see that the agent trades fast when risk aversion γ is large or when k is small. This is clearly instanced by Figure 3. Also, the agent trades faster when the price is expected to decrease to avoid selling at low price. This is clearly instanced by Figure 4.

General case

Now, let's turn to the general case to take account of volatility and thus of the price risk dimension of the problem.

As far as the drift is concerned, quotes are naturally increasing with µ. If indeed the trader expects the price to move down, he is going to send orders at low prices to be executed fast and thus reduce the risk to suffer from the decrease in price. This is well exemplified in Figure 5. 2 ), A = 0.1 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

Now, coming to volatility, the optimal quotes depend on σ in a monotonic way. If there is an increase in volatility, the price risk increases. In order to reduce this additional price risk the trader will send orders at cheaper price. This is what we observe numerically on Figure 6. 2 ) (dots), σ = 0.3 (Tick.s -1 2 ) (solid line), σ = 0.6 (Tick.s -1

2 ) (dotted line), A = 0.1 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

Now, coming to A, we observe numerically that the optimal quote is an increasing function of A (see Figure 7). 2 ), A = 0.05 (s -1 ) (dots), A = 0.10 (s -1 ) (solid line), A = 0.15 (s -1 ) (dotted line), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

Then as far as the dependence of the optimal quotes on k is concerned, the above discussion on the interpretation of k is still valid. Since high volatility may more often induce optimal quotes below the reference price, the dependence on k is not the same at all times and for all values of the inventory (Figure 8 and Figure 9). 2 ), A = 0.10 (s -1 ), k = 0.2 (Tick -1 ) (dots), k = 0.3 (Tick -1 ) (solid line), k = 0.4 (Tick -1 ) (dotted line) and γ = 0.05 (Tick -1 )) 2 ), A = 0.10 (s -1 ), k = 0.2 (Tick -1 ) (dots), k = 0.3 (Tick -1 ) (solid line), k = 0.4 (Tick -1 ) (dotted line) and γ = 0.05 (Tick -1 ))

Finally, turning to the risk aversion γ, two effects are at stake that go in the same direction. The risk aversion is indeed common for both price risk and execution risk. Hence if risk aversion increases, the trader will try to reduce both price risk and execution risk, thus selling at cheaper price. We indeed see on Figure 10 that optimal quotes are decreasing in γ. 2 ), A = 0.10 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.01 (Tick -1 ) (dots), γ = 0.05 (Tick -1 ) (solid line), γ = 0.1 (Tick -1 ) (dotted line))

Applications

Before using the above model in reality, we need to discuss some features of the model that need to be adapted before any backtest is possible. First of all, the model is continuous in both time and space while the real control problem under scrutiny is intrinsically discrete in space, because of the tick size, and in time, because orders have a certain priority and changing position too often reduces the actual chance to be reached by a market order. Hence, the model has to be reinterpreted in a discrete way. In terms of prices, quotes must not be between two ticks and we decided to ceil or floor the optimal quotes with probabilities that depend on the respective proximity to the neighboring quotes. In terms of time, an order is sent to the market and is not canceled nor modified for a given period ∆t, unless a trade occurs and, though perhaps partially, fills the order. Now, when a trade occurs and changes the inventory or when an order stayed in the order book for longer than ∆t, then the optimal quote is updated and, if necessary, a new order is inserted. Now, concerning the parameters, σ, A and k can be calibrated easily on trade-by-trade limit order book data while γ has to be chosen. However, it is well known by practitioners that A and k have to depend at least on the actual market bid-ask spread. Since we do not explicitly take into account the underlying market, there is no market bid-ask spread in the model. Thus, we simply chose to calibrate A and k as functions of the market bid-ask spread, making then an off-model hypothesis.

Turning to the backtests, they were carried out with trade-by-trade data and we assumed that our orders were entirely filled when a trade occurred above the ask price quoted by the agent. Our goal here is just to provide examples in various situations and, to exemplify the practical use of this model, we carried out several backtests13 on the French stock AXA, either on very short periods (5 to 10 minutes) or on slightly longer periods of a few hours. 11 and12) consist in getting rid of a quantity of shares equal to 3 times the ATS 14 . The periods have been chosen to capture the behavior in both bullish and bearish markets.

The first two examples (Figures

On Figure 11, we see that the first order is executed after 50 seconds. Then, since the trader has only 2 times the ATS left in his inventory, he sends an order at a higher price. Since the market price moves up, the second order is executed in the next 30 seconds, in advance on the average schedule. This is the reason why the trader places a new order far above the best ask. Since this order is not executed within the time window ∆t, it is canceled and new orders are successively inserted with lower prices. The last trade happens less than 1 minute before the end of the period. Overall, on this example, the strategy works far better than a market order (even ignoring execution costs).

On Figure 12, we see the use of the strategy in a bearish period. The first order is executed rapidly and since the market price goes down, the trader's last orders are only executed at the end of the period when prices of orders are lowered substantially since selling becomes of utmost importance. Practically, this obviously raises the question of linking a trend detector to these optimal liquidation algorithms. -1 ) to get rid of a quantity of shares equal to 3 times the ATS within 5 minutes. Top: quotes of the trader (bold line), market best bid and ask quotes (thin lines). Trades are represented by dots. Bottom left: evolution of the inventory. Bottom right: P&L. Now, another example is shown on Figure 13 where a trader wants to trade a quantity of shares equal to 6 times the ATS within 10 minutes. The analysis is similar to the preceding ones, the market being even more bullish than in the first example. Finally, the model can also be used on longer periods and we exhibit the use of the algorithm on a period of two hours, to sell a quantity of shares equal to 20 times the ATS (Figure 14).

Conclusion

As claimed in the introduction, this paper is, to authors' knowledge, the first proposal to optimize the trade scheduling of large orders with small passive orders. Thanks to an innovative change of variables, it provides an explicit solution in two steps: (1) solve an ODEs, (2) deduce the optimal price of the order to be sent to the market.

The choices of modeling made here have been more extensively discussed in [START_REF] Guéant | Dealing with the inventory risk[END_REF]. It is nevertheless worthy to underline that no explicit model of what could be called "passive market impact" (i.e. the perturbations of the price formation process by liquidity provision) is used here. Just note that up to now, no quantitative model of this type of impact has been proposed in the literature. Thanks to very promising and recent studies of the multidimensional point processes governing the arrival of orders (see for instance the link between the imbalance in the order flow and the moves of the price studied in [START_REF] Cont | Price Dynamics in a Markovian Limit Order Book Market[END_REF] or [START_REF] Cont | The Price Impact of Order Book Events[END_REF], or interesting properties of Hawkes-like models in [START_REF] Bacry | Modeling microstructure noise with mutually exciting point processes[END_REF]), we can hope for obtaining models of this kind in the near future. The authors will try to embed them into the HJB framework used here. In between, an on-going work dedicated to model dependencies between the Brownian motion supporting S t and the Poisson process N a t is under consideration.

Appendix

Proof of Proposition 1 and Theorem 1:

First, let us remark that a solution (w q ) q of (S) exists and is unique and that, by immediate induction, its components are strictly positive for all times. Then, let's introduce u(t, x, q, s) = -exp (-γ(x + qs)) w q (t) -γ k .

We have:

∂ t u + µ∂ s u + 1 2 σ 2 ∂ 2 ss u = - γ k ẇq (t) w q (t) u -γqµu + γ 2 σ 2 2 q 2 u
Now, concerning the non-local part of the equation, we have:

sup δ a λ a (δ a ) [u(t, x + s + δ a , q -1, s) -u(t, x, q, s)] = sup δ a
Ae -kδ a u(t, x, q, s)

[ exp (-γδ a ) ( w q-1 (t) w q (t) ) -γ k -1 ]
The first order condition of this problem corresponds to a maximum and writes:

(k + γ) exp (-γδ a ) ( w q-1 (t) w q (t)

) -γ k = k
Hence we introduce the candidate δ a * for the optimal control:

δ a * = 1 k ln ( w q (t) w q-1 (t) ) + 1 γ ln ( 1 + γ k
) and sup δ a λ a (δ a ) [u(t, x + s + δ a , q -1, s) -u(t, x, q, s)]

= - γ k + γ
A exp(-kδ a * )u(t, x, q, s)

= -A γ k + γ ( 1 + γ k ) -k γ w q-1 (t) w q (t) u(t, x, q, s)
Hence, putting the three terms together we get:

∂ t u(t, x, q, s) + µ∂ s u(t, x, q, s) + 1 2 σ 2 ∂ 2 ss u(t, x, q, s) + sup δ a λ a (δ a ) [u(t, x + s + δ a , q -1, s) -u(t, x, q, s)] = - γ k ẇq (t) w q (t) u -γµqu + γ 2 σ 2 2 q 2 u -A γ k + γ ( 1 + γ k ) -k γ w q-1 (t) w q (t) u = - γ k u w q (t) [ ẇq (t) + kµqw q (t) - kγσ 2 2 q 2 w q (t) + A ( 1 + γ k ) - ( 1+ k γ ) w q-1 (t) ] = 0
Now, noticing that the boundary and terminal conditions for w q are consistent with the conditions on u, we get that u verifies (HJB). Now, we need to verify that u is indeed the value function associated to the problem and to prove that our candidate (δ a * t ) t is indeed the optimal control. To that purpose, let's consider a control ν ∈ A and let's consider the following processes for τ ∈ [t, T ]:

dS t,s τ = µdτ + σdW τ , S t,s t = s dX t,x,ν τ = (S τ + ν τ )dN a τ , X t,x,ν t = x dq t,q,ν τ = -dN a τ , q t,q,ν t = q
where the Poisson process has stochastic intensity (λ τ ) τ with λ τ = Ae -kντ 1 q τ -≥1 15 .

Now, let us write Itô's formula for u since u is smooth:

u(T, X t,x,ν T -, q t,q,ν T -, S t,s T ) = u(t, x, q, s) + ∫ T t ( ∂ τ u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) + µ∂ s u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) + σ 2 2 ∂ 2 ss u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) ) dτ + ∫ T t ( u(τ, X t,x,ν τ -+ S t,s τ + ν τ , q t,q,ν τ --1, S t,s τ ) -u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) ) λ τ dτ + ∫ T t σ∂ s u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ )dW τ + ∫ T t ( u(τ, X t,x,ν τ -+ S t,s τ + ν τ , q t,q,ν τ --1, S t,s τ ) -u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) ) dM a τ
where M a is the compensated Poisson process associated to N a for the intensity process (λ τ ) τ . Now, we have to ensure that the last two integrals consist of martingales so that their mean is 0. To that purpose, let us notice that ∂ s u = -γu and hence we just have to prove that:

E [∫ T t u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) 2 dτ ] < +∞ E [∫ T t u(τ, X t,x,ν τ -+ S t,s τ + ν τ , q t,q,ν τ --1, S t,s τ ) λ τ dτ ] < +∞ and E [∫ T t u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) λ τ dτ ] < +∞
Now, remember that the process q t,q,ν takes values between 0 and q and that t ∈ [0, T ]. Hence, ∃ε > 0, w q (t) > ε for the values of t and q under scrutiny and:

u(τ, X t,x,ν τ , q t,q,ν τ , S t,s τ ) 2 ≤ ε -2γ k exp ( -2γ(X t,x,ν τ + q t,q,ν τ S t,s τ ) ) ≤ ε -2γ k exp ( -2γ(x -q∥ν -∥ ∞ + 2q inf τ ∈[t,T ] S t,s τ 1 inf τ ∈[t,T ] S t,s τ <0 ) ) ≤ ε -2γ k exp (-2γ(x -q∥ν -∥ ∞ )) ( 1 + exp ( -2γq inf τ ∈[t,T ] S t,s τ ))
Hence:

E [∫ T t u(τ, X t,x,ν τ , q t,q,ν τ , S t,s τ ) 2 dτ ] = E [∫ T t u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) 2 dτ ] ≤ ε -2γ k exp (-2γ(x -q∥ν -∥ ∞ ))(T -t) ( 1 + E [ exp ( -2γq inf τ ∈[t,T ] S t,s τ )]) ≤ ε -2γ k exp (-2γ(x -q∥ν -∥ ∞ ))(T -t) ( 1 + E [ exp ( -2γq inf τ ∈[t,T ] S t,s τ )]) ≤ ε -2γ k exp (-2γ(x -q∥ν -∥ ∞ ))(T -t) ( 1 + e -2γqs E [ exp ( 2γqσ √ T -t|Y | )]) < +∞
where the last inequalities come from the reflection principle with Y ∼ N (0, 1) and the fact that E [ e C|Y | ] < +∞ for any C ∈ R. Now, the same argument works for the second and third integrals, noticing that ν is bounded from below and that λ is bounded.

Hence, since we have, by construction 16∂ τ u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) + µ∂ s u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) + σ 2 2 ∂ 2 ss u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ ) + ( u(τ, X t,x,ν τ -+ S t,s τ + ν t , q t,q,ν τ --1, S t,s τ ) -u(τ, X t,x,ν τ -, q t,q,ν τ -, S t,s τ )

) λ τ ≤ 0 we obtain that E [ u(T, X t,x,ν T , q t,q,ν T , S t,s T )

] = E
[ u(T, X t,x,ν T -, q t,q,ν T -, S t,s T )

] ≤ u(t, x, q, s) and this is true for all ν ∈ A. Since for ν = δ a * we have an equality in the above inequality we obtain that:

sup ν∈A E [ u(T, X t,x,ν T , q t,q,ν T , S t,s T )

] ≤ u(t, x, q, s) = E [ u(T, X t,x,δ a * T , q t,q,δ a * T , S t,s T )

]

This proves that u is the value function and that δ a * is optimal.

Proof of Proposition 2:

We have that ∀q ∈ N, ẇq (t) = (αq 2 -βq)w q (t) -ηw q-1 (t)

Hence if we consider for a given Q ∈ N the vector w(t) =      w 0 (t) w 1 (t) . . .

w Q (t)     
we have that w ′ (t) = M w(t) where:

M =             0 0 • • • • • • • • • 0 -η α -
(Q -1) 2 -β(Q -1) 0 0 • • • • • • 0 -η αQ 2 -βQ             with w(T ) =      1 e -kb
. . .

e -kbQ     
. Hence we know that, if we consider a basis (f 0 , . . . , f Q ) of eigenvectors (f j being associated to the eigenvalue αj 2 -βj), there exists (c 0 , . . . , c Q ) ∈ R Q+1 independent of T such that:

w(t) = Q ∑ j=0
c j e -(αj 2 -βj)(T -t) f j Consequently, since we assumed that α > β, we have that w ∞ := lim T →+∞ w(0) = c 0 f 0 . Now, w ∞ is characterized by: (αq 2 -βq)w ∞ q = ηw ∞ q-1 , q > 0 w ∞ 0 = 1

As a consequence we have:

w ∞ q = η q q! q ∏ j=1 1 αj -β
The resulting asymptotic behavior for the optimal ask quote is: lim T →+∞ δ a * (0, q) = 1 k ln

( A 1 + γ k 1 αq 2 -βq )

Proof of Proposition 3:

The result of Proposition 3 is obtained by induction. For q = 0 the result is obvious. Now, if the result is true for some q we have that:

ẇq+1 (t) = - q ∑ j=0 η j+1
j! e -kb(q-j) (T -t) j

Hence:

w q+1 (t) = e -kb(q+1) + q ∑ j=0 η j+1 (j + 1)! e -kb(q-j) (T -t) j+1 w q+1 (t) = e -kb(q+1) + q+1 ∑ j=1 η j j! e -kb(q-j+1) (T -t) j w q+1 (t) = q+1 ∑ j=0 η j j! e -kb(q+1-j) (T -t) j

This proves the results for w and then the result follows for the optimal quote.

Proof of Proposition 4:

Because the solutions depend continuously on b, we can directly get interested in the limiting equation: ∀q ∈ N, ẇq (t) = (αq 2 -βq)w q (t) -ηw q-1 (t) with w q (T ) = 1 q=0 and w 0 = 1.

Then, if we define v q (t) = lim b→+∞ w b,q (t) A q , v solves: ∀q ∈ N, vq (t) = (αq 2 -βq)v q (t) -ηv q-1 (t) with v q (T ) = 1 q=0 and v 0 = 1, where η = η A is independent of A.

Figure 1 :

 1 Figure1: Optimal strategy δ a * (t, q) (in Ticks) for an agent willing to sell a quantity of shares up to 6 times the ATS within 5 minutes (µ = 0 (Tick.s -1 ), σ = 0.3 (Tick.s -1 2 ), A = 0.1 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

Figure 2 :

 2 Figure 2: Solid line: Trading curve for an agent willing to sell a quantity of shares equal to 6 times the ATS within 5 minutes (µ = 0 (Tick.s -1 ), σ = 0.3 (Tick.s -12 ), A = 0.1 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 )). Dotted line: Trading curve calculated using the Almgren-Chriss framework with the same parameters and the impact function f (v) = 20v 1+ 3 4 .

Figure 3 :

 3 Figure 3: Trading curves of an agent willing to sell 6 times the ATS for different sets of parameters (in the no-drift case). Dots: risk neutral case γ = 0. Dotted line: k = 0.3 (Tick -1 ), γ = 0.05 (Tick -1 ). Solid line: k = 0.2 (Tick -1 ), γ = 0.1 (Tick -1 )

Figure 4 :

 4 Figure 4: Trading curves of an agent willing to sell 6 times the ATS for different values of the drift. µ = 1/150 (Tick.s -1 ) (dots), µ = 0 (Tick.s -1 ) (dotted line), µ = -1/150 (Tick.s -1 ) (solid line) k = 0.2 (Tick -1 ), γ = 0.1 (Tick -1 )

Figure 5 :

 5 Figure 5: Dependence on µ of δ a * for an agent willing to sell a quantity of shares up to 6 times the ATS in the next 5 minutes. (µ = -1/150 (Tick.s -1 ) (dots), µ = 0 (Tick.s -1 ) (dotted line), µ = 1/150 (Tick.s -1 ) (solid line), σ = 0.3 (Tick.s -12 ), A = 0.1 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

Figure 6 :

 6 Figure 6: Dependence on σ of δ a * for an agent willing to sell a quantity of shares up to 6 times the ATS in the next 5 minutes. (µ = 0 (Tick.s -1 ), σ = 0 (Tick.s -12 ) (dots), σ = 0.3 (Tick.s -1 2 ) (solid line), σ = 0.6 (Tick.s -12 ) (dotted line), A = 0.1 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

Figure 7 :

 7 Figure 7: Dependence on A of δ a * for an agent willing to sell a quantity of shares up to 6 times the ATS in the next 5 minutes. (µ = 0 (Tick.s -1 ), σ = 0.3 (Tick.s -12 ), A = 0.05 (s -1 ) (dots), A = 0.10 (s -1 ) (solid line), A = 0.15 (s -1 ) (dotted line), k = 0.3 (Tick -1 ) and γ = 0.05 (Tick -1 ))

Figure 8 :

 8 Figure 8: Dependence on k of δ a * for an agent willing to sell a quantity of shares up to 6 times the ATS in the next 5 minutes. (µ = 0 (Tick.s -1 ), σ = 0.3 (Tick.s -12 ), A = 0.10 (s -1 ), k = 0.2 (Tick -1 ) (dots), k = 0.3 (Tick -1 ) (solid line), k = 0.4 (Tick -1 ) (dotted line) and γ = 0.05 (Tick -1 ))

Figure 9 :

 9 Figure 9: Dependence on k of δ a * for an agent willing to sell a quantity of shares up to 6 times the ATS in the next 5 minutes. (µ = 0 (Tick.s -1 ), σ = 3 (Tick.s -12 ), A = 0.10 (s -1 ), k = 0.2 (Tick -1 ) (dots), k = 0.3 (Tick -1 ) (solid line), k = 0.4 (Tick -1 ) (dotted line) and γ = 0.05 (Tick -1 ))

Figure 10 :

 10 Figure 10: Dependence on γ of δ a * for an agent willing to sell a quantity of shares up to 6 times the ATS in the next 5 minutes. (µ = 0 (Tick.s -1 ), σ = 0.3 (Tick.s -12 ), A = 0.10 (s -1 ), k = 0.3 (Tick -1 ) and γ = 0.01 (Tick -1 ) (dots), γ = 0.05 (Tick -1 ) (solid line), γ = 0.1 (Tick -1 ) (dotted line))

Figure 11 :Figure 12 :

 1112 Figure 11: Backtest example on AXA (November 5 th 2010). The strategy is used with γ = 10 (euro -1 ) to get rid of a quantity of shares equal to 3 times the ATS within 5 minutes. Top: quotes of the trader (bold line), market best bid and ask quotes (thin lines). Trades are represented by dots. Bottom left: evolution of the inventory. Bottom right: P&L.

Figure 13 :

 13 Figure 13: Backtest example on AXA (November 5 th 2010). The strategy is used with γ = 10 (euro -1 ) to get rid of a quantity of shares equal to 6 times the ATS within 10 minutes. Top: quotes of the trader (bold line), market best bid and ask quotes (thin lines). Trades are represented by dots. Bottom left: evolution of the inventory. Bottom right: P&L.

Figure 14 :

 14 Figure 14: Backtest example on AXA (November 8 th 2010). The strategy is used with γ = 1 (euro -1 ) to get rid of a quantity of shares equal to 20 times the ATS within 2 hours. Top: quotes of the trader (bold line), market best bid and ask quotes (thin lines). Trades are represented by dots. Bottom left: evolution of the inventory. Bottom right: P&L.

see[START_REF] Foucault | Competition for Order Flow and Smart Order Routing Systems[END_REF] for a study of the effect of "smart order routing" on competitive trading venues, or[START_REF] Bloomfield | The "Make or Take" Decision in an Electronic Market: Evidence on the Evolution of Liquidity[END_REF] for a study of the efficiency of order matching mechanisms.

Kratz and Schöneborn [25] proposed an approach with both market orders and access to dark pools.

The actual models studying the SOR problem across several venues are more focused on routing across Dark Pools specifically (see[START_REF] Nevmyvaka | Electronic Trading in Order-Driven Markets: Efficient Execution[END_REF] for a statistical approach or[START_REF] Pagès | Optimal split of orders across liquidity pools: a stochatic algorithm approach[END_REF] for a probabilistic one).

It recently came to our knowledge that a similar approach for liquidation with limit orders is being developed independently by E. Bayraktar and M. Ludkovski[START_REF] Bayraktar | Liquidation in Limit Order Books with Controlled Intensity[END_REF]. Their approach uses the same framework as ours to describe the price process and the execution mechanism. However, they only consider risk-neutral traders and consequently they do not take account of risk, be it price risk, or execution risk.

Once the whole portfolio is liquidated, we assume that the trader remains inactive.

It is important to notice that 1 share may be understood as 1 bunch of shares, each bunch being of the same size, typically the average trade size (hereafter ATS) or a fraction of it. If one wants to replace orders of size 1 by orders of size δq in the model, this can be done easily. However, the framework of the model imposes to trade with orders of constant size, an hypothesis that is an approximation of reality since orders may in practice be partially filled.

In all computations, we considered very large values for b to force the trader to use every endeavor to liquidate the portfolio before time T .

The authors most often consider target problems in which the target can always be attained.

In our case indeed, since the execution process is not deterministic, the trading curve associated to the optimal strategy can only be defined on average.

No drift in prices is assumed in the strategy used for backtesting.

In the backtests we do not deal with quantity and priority issues in the order books and supposed that our orders were always entirely filled.

This intensity being bounded since ν is bounded from below.

This inequality is also true when the portfolio is empty because of the boundary conditions.

Hence v q (t) is independent of A. Now, for the trading intensity we have:

and this does not depend on A.

Proof of Proposition 5:

Using the preceding proposition, we can now reason in terms of v and look for a solution of the form v q (t) = h(t) q q! . Then,

From Theorem 1, we obtain the limit of the optimal quote:

Using the expression for η, this can also be written:

Now, the for the trading intensity we get:

) q β e β(T -t) -1 Hence, because the limit of the intensity is proportional to q, the limit V (t) of the trading curve is characterized by the following ODE:

Solving this equation, we get:

) 1+ γ k When β = 0 (i.e. µ = 0) we proceed in the same way or by a continuity argument.


Olivier Guéant
email: olivier.gueant@ann.jussieu.fr

Mean field games equations with quadratic Hamiltonian: a specific approach *

. We use a change of variables to transform the mean field games (MFG) equations into a system of simpler coupled partial differential equations, in the case of a quadratic Hamiltonian. This system is then used to exhibit a monotonic scheme to build solutions of the MFG equations. Effective numerical methods based on this constructive scheme are presented and numerical experiments are carried out.

Introduction

Mean field games equations have been introduced by J.-M. Lasry and P.-L. Lions [START_REF] Lasry | Jeux à champ moyen i. le cas stationnaire[END_REF][START_REF] Lasry | Jeux à champ moyen ii. horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF] to describe the dynamic equilibrium of stochastic differential games involving a continuum of players.

In the time-dependent case, these equations write:

(HJB)

∂ t u + σ 2 2 ∆u + H(∇u) = -f (x, m) (K) ∂ t m + ∇ • (mH (∇u)) = σ 2 2 ∆m
with prescribed initial condition m(0, •) = m 0 (•) ≥ 0 and terminal condition u(T, •) = u T (•), where u and m are scalar functions defined on [0, T ] × Ω, Ω typically being (0, 1) d .

In this paper, we focus on the particular case of quadratic hamiltonian H(p) = p 2 2 . In this special case, a change of variables have been introduced by O. Guéant, J.-M. Lasry and P.-L. Lions in [START_REF] Guéant | Mean field games and applications[END_REF] to write the mean field games equations as two coupled heat equations with similar source terms. If indeed we introduce φ = exp u σ 2 and ψ = m exp -u σ 2 then the system reduces to:

∂ t φ + σ 2 2 ∆φ = - 1 σ 2 f (x, φψ)φ ∂ t ψ - σ 2 2 ∆ψ = 1 σ 2 f (x, φψ)ψ with φ(T, •) = exp u T (•) σ 2
and ψ(0, •) = m 0 (•) φ(0,•) .

We use this system to exhibit a constructive scheme for solutions to the mean field games equations. This scheme and the method used to prove it indeed converges have been proposed by J.-M. Lasry and P.-L. Lions. It starts with ψ 0 = 0 and builds recursively two sequences (φ n+ 1 2 ) n and (ψ n+1 ) n using the following equations:

∂ t φ n+ 1 2 + σ 2 2 ∆φ n+ 1 2 = - 1 σ 2 f (x, φ n+ 1 2 ψ n )φ n+ 1 2 ∂ t ψ n+1 - σ 2 2 ∆ψ n+1 = 1 σ 2 f (x, φ n+ 1 2 ψ n+1 )ψ n+1 with φ n+ 1 2 (T, •) = exp u T (•) σ 2
and ψ n+1 (0, •) = m 0 (•)

φ n+ 1 2 (0,•)
.

Then, φ and ψ are obtained as the monotonic limit of the two sequences (φ n+ 1 2 ) n and (ψ n ) n under usual assumptions of f . Moreover, this scheme provides a new numerical method to solve mean field games that is completely different from the methods already proposed in the literature. Finite different schemes on u and m have been proposed, along with methods based on an equivalent formulation of the mean field game partial differential equations in the form of an optimization problem (see [START_REF] Achdou | Mean field games: numerical methods for the planning problem[END_REF], [START_REF] Achdou | Mean field games: Numerical methods[END_REF], [START_REF] Lachapelle | Computation of mean field equilibria in economics[END_REF], [START_REF] Guéant | A reference case for mean field games models[END_REF] for the different proposals). Here, we build a discrete counterpart to the above constructive scheme and the resulting algorithm to approximate (φ, ψ) consists in a sequence of totally implicit finite difference schemes that can be tackled with Newton methods. The main difference with the preceding literature arises from the monotonicity properties of the scheme and from the peculiarity that the equilibrium m is approached by a sequence that does not verify, except at the limit, the mass conservation principle.

In section 1, we recall the change of variables and derive the associated system of coupled parabolic equations. Then, section 2 is devoted to the introduction of the functional framework and we prove the main monotonicity properties of the system. Section 3 presents the constructive scheme and proves, as explained in P.-L. Lions's course (see [START_REF] Lions | Cours au collège de france : Théorie des jeux à champs moyens[END_REF]), that we can have two monotonic sequences converging towards φ and ψ. Section 4 uses the same ideas as in the preceding sections, but in a discrete setting, to provide numerical schemes to approximate φ and ψ numerically. Stability and convergence of the scheme are proved. Finally, section 5 presents the numerical experiments carried out and discusses its properties.

1 From mean field games equations to a forward-backward system of heat equations with source terms

We consider the mean field games equations introduced in [START_REF] Lasry | Jeux à champ moyen i. le cas stationnaire[END_REF][START_REF] Lasry | Jeux à champ moyen ii. horizon fini et contrôle optimal[END_REF][START_REF] Lasry | Mean field games[END_REF] in the case of a quadratic Hamiltonian. These partial differential equations, hereafter denoted (MFG) are considered on the domain [0, T ] × Ω, Ω standing for (0, 1) d , and consist in the following equations:

(HJB)

∂ t u + σ 2 2 ∆u + 1 2 |∇u| 2 = -f (x, m) (K) ∂ t m + ∇ • (m∇u) = σ 2 2 ∆m
with:

• Boundary conditions: ∂u ∂ n = ∂m ∂ n = 0 on (0, T ) × ∂Ω • Terminal condition: u(T, •) = u T (•) a given payoff whose regularity is to be precised • Initial condition: m(0, •) = m 0 (•) ≥ 0 a given positive function in L 1 (Ω), typically a probability distribution function.

The change of variables introduced in [START_REF] Guéant | Mean field games and applications[END_REF] is recalled in the following proposition:

Proposition 1. Let us consider a smooth solution (φ, ψ) of the following system (S) with φ > 0:

∂ t φ + σ 2 2 ∆φ = - 1 σ 2 f (x, φψ)φ (E φ ) ∂ t ψ - σ 2 2 ∆ψ = 1 σ 2 f (x, φψ)ψ (E ψ )
with:

• Boundary conditions: ∂φ ∂ n = ∂ψ ∂ n = 0 on (0, T ) × ∂Ω • Terminal condition: φ(T, •) = exp u T (•) σ 2 . • Initial condition: ψ(0, •) = m 0 (•) φ(0,•)
Then (u, m) = (σ 2 ln(φ), φψ) defines a solution of (MFG).

Proof:

Let us start with (HJB).

∂ t u = σ 2 ∂ t φ φ , ∇u = σ 2 ∇φ φ ∆u = σ 2 ∆φ φ -σ 2 |∇φ| 2 φ 2 Hence ∂ t u + σ 2 2 ∆u + 1 2 |∇u| 2 = σ 2 ∂ t φ φ + σ 2 2 ∆φ φ = σ 2 φ - 1 σ 2 f (x, φψ)φ = -f (x, m)
Now, for the equation (K):

∂ t m = ∂ t φψ + φ∂ t ψ ∇ • (∇um) = σ 2 ∇ • (∇φψ) = σ 2 [∆φψ + ∇φ • ∇ψ] ∆m = ∆φψ + 2∇φ • ∇ψ + φ∆ψ Hence ∂ t m + ∇ • (∇um) = ∂ t φψ + φ∂ t ψ + σ 2 [∆φψ + ∇φ • ∇ψ] = ψ ∂ t φ + σ 2 ∆φ + φ∂ t ψ + σ 2 ∇φ • ∇ψ = ψ σ 2 2 ∆φ - 1 σ 2 f (x, φψ)φ + φ σ 2 2 ∆ψ + 1 σ 2 f (x, φψ)ψ + σ 2 ∇φ • ∇ψ = σ 2 2 ∆φψ + σ 2 ∇φ • ∇ψ + σ 2 2 φ∆ψ = σ 2 2 ∆m
This proves the result since the boundary conditions and the initial and terminal conditions are coherent. Now, we will focus our attention on the study of the above system of equations (S) and use it to design a constructive scheme for the couple (φ, ψ) and thus for the couple (u, m) under regularity assumptions.

Properties of (S)

To study the system (S) we introduce several hypotheses on f : we suppose that it is a decreasing function of its second variable, continuous in that variable and uniformly bounded. Moreover, to simplify the exposition1 we suppose that f ≤ 0. It's noteworthy that the monotony hypothesis is to be linked to the usual proof of uniqueness for the mean field games equations (see [START_REF] Lasry | Mean field games[END_REF]). Now let's introduce the functional framework we are working in. Let us note P ⊂ C([0, T ], L 2 (Ω)) the natural set for parabolic equations:

g ∈ P ⇐⇒ g ∈ L 2 (0, T, H 1 (Ω)) and ∂ t g ∈ L 2 (0, T, H -1 (Ω))
and let's introduce P = {g ∈ P, g ≥ }.

Proposition 2. Suppose that u T ∈ L ∞ (Ω). ∀ψ ∈ P 0 , there is a unique weak solution φ to the following equation (E φ ):

∂ t φ + σ 2 2 ∆φ = - 1 σ 2 f (x, φψ)φ (E φ ) with ∂φ ∂ n = 0 on (0, T ) × ∂Ω and φ(T, •) = exp u T (•) σ 2 .
Hence Φ : ψ ∈ P 0 → φ ∈ P is well defined.

Moreover, ∀ψ ∈ P 0 , φ = Φ(ψ)

∈ P for = exp -1 σ 2 ( u T ∞ + f ∞ T ) Proof:
Let us consider ψ ∈ P 0 .

Existence of a weak solution φ:

Let us introduce F ψ : ϕ ∈ L 2 (0, T, L 2 (Ω)) → φ weak solution of the following linear parabolic equation:

∂ t φ + σ 2 2 ∆φ = - 1 σ 2 f (x, ϕψ)φ with ∂φ ∂ n = 0 on (0, T ) × ∂Ω and φ(T, •) = exp u T (•) σ 2
.

By classical linear parabolic equations theory, φ is in P ⊂ L 2 (0, T, L 2 (Ω)).

Our goal is to use Schauder's fixed-point theorem on F ψ .

Compactness:

Usual energy estimates (see [START_REF] Evans | Partial Differential Equations[END_REF]) give that there exists a constant C that only depends on u T ∞ , σ and f ∞ such that:

∀(ψ, ϕ) ∈ P 0 × L 2 (0, T, L 2 (Ω)), F ψ (ϕ) L 2 (0,T,H 1 (Ω)) + ∂ t F ψ (ϕ) L 2 (0,T,H -1 (Ω)) ≤ C Hence F ψ maps the closed ball B L 2 (0,T,L 2 (Ω)) (0, C) to a compact subset of B L 2 (0,T,L 2 (Ω)) (0, C).

Continuity:

Let us now prove that F ψ is a continuous function.

Let us consider a sequence (ϕ n ) n of L 2 (0, T, L 2 (Ω)) with ϕ n → n→∞ ϕ in the L 2 (0, T, L 2 (Ω)) sense. Let us write φ n = F ψ (ϕ n ). We know from the above compactness result that we can extract from (φ n ) n a new sequence denoted (φ n ) n that converges in the L 2 (0, T, L 2 (Ω)) sense toward a function φ. To prove that F ψ is continuous, we then need to show that φ cannot be different from F ψ (ϕ). Now, because of the energy estimates, we know that φ is in P and that we can extract another subsequence (still denoted (φ n ) n ), such that:

• φ n → φ in the L 2 (0, T, L 2 (Ω)) sense. • ∇φ n ∇φ weakly in L 2 (0, T, L 2 (Ω)). • ∂ t φ n ∂ t φ weakly in L 2 (0, T, H -1 (Ω)).
and

• ϕ n → ϕ almost everywhere.
By definition we have that ∀w ∈ L 2 (0, T, H 1 (Ω)):

T 0 ∂ t φ n (t, •), w(t, •) H -1 (Ω),H 1 (Ω) dt - σ 2 2 T 0 Ω ∇φ n (t, x) • ∇w(t, x)dxdt = - 1 σ 2 T 0 Ω f (x, ϕ n (t, x)ψ(t, x))φ n (t, x)w(t, x)dxdt
By weak convergence, the left hand side of the above equality converges toward

T 0 ∂ t φ(t, •), w(t, •) H -1 (Ω),H 1 (Ω) dt - σ 2 2 T 0 Ω ∇φ(t, x) • ∇w(t, x)dxdt
Since f is a bounded continuous function and using dominated convergence theorem, the right hand side converges toward

- 1 σ 2 T 0 Ω f (x, ϕ(t, x)ψ(t, x))φ(t, x)w(t, x)dxdt Hence φ = F ψ (ϕ).
Schauder's Theorem: By Schauder's theorem, we then know that there exists a fixed-point φ to F ψ and hence a weak solution to the nonlinear parabolic equation (E φ ).

Positiveness of φ:

Let us consider a solution φ as above. If I(t) = 1 2 Ω (φ(t, x) -) 2 dx, then:

I (t) = - Ω ∂ t φ(t, x)φ(t, x) -dx = - Ω ∇φ(t, x) • ∇(φ(t, x) -) - 1 σ 2 f (x, φ(t, x)ψ(t, x))φ(t, x)φ(t, x) -dx = - Ω -|∇φ(t, x)| 2 1 φ(t,x)≤0 + 1 σ 2 f (x, φ(t, x)ψ(t, x))(φ(t, x) -) 2 dx = Ω |∇φ(t, x)| 2 1 φ(t,x)≤0 - Ω 1 σ 2 f (x, φ(t, x)ψ(t, x))(φ(t, x) -) 2 dx ≥ 0
Since I(T ) = 0 and I ≥ 0, we know that I = 0. Hence, φ is positive.

Uniqueness:

Let us consider two weak solutions φ 1 and φ 2 to the equation (E φ ).

Let us introduce J(t) = 1 2 Ω (φ 2 (t, x) -φ 1 (t, x)) 2 dx. We have:

J (t) = Ω (∂ t φ 2 (t, x) -∂ t φ 1 (t, x))(φ 2 (t, x) -φ 1 (t, x))dx = - Ω 1 σ 2 (f (x, φ 2 (t, x)ψ(t, x))φ 2 (t, x) -f (x, φ 1 (t, x)ψ(t, x))φ 1 (t, x)) (φ 2 (t, x) -φ 1 (t, x))dx + Ω |∇φ 2 (t, x) -∇φ 1 (t, x)| 2 dx
Because of our assumptions on f , the function ξ ∈ R + → 1 σ 2 f (x, ψξ)ξ is a decreasing function.

Hence, since φ 1 and φ 2 are positive, J (t) ≥ 0. Since J(T ) = 0 and J ≥ 0, we know that J = 0. Hence, φ 1 = φ 2 .

Lower bound to φ:

We can get a lower bound to φ through a subsolution taken as the solution of the following ordinary differential equation:

φ (t) = 1 σ 2 f ∞ φ(t) φ(T ) = exp - u T ∞ σ 2
Let us indeed consider K(t) = 1 2 Ω ((φ(t) -φ(t, x)) + ) 2 dx. We have:

K (t) = Ω (φ (t) -∂ t φ(t, x))(φ(t) -φ(t, x)) + dx = Ω 1 σ 2 f ∞ φ(t)(φ(t) -φ(t, x)) + + |∇φ(t, x)| 2 1 φ(t)-φ(t,x)≥0 + 1 σ 2 f (x, φ(t, x)ψ(t, x))φ(t, x)(φ(t) -φ(t, x)) + dx ≥ 1 σ 2 Ω f ∞ φ(t) + f (x, φ(t, x)ψ(t, x))φ(t, x) (φ(t) -φ(t, x)) + dx ≥ 1 σ 2 Ω ( f ∞ + f (x, φ(t, x)ψ(t, x))) φ(t, x)(φ(t) -φ(t, x)) + dx ≥ 0
Since K(T ) = 0 and K ≥ 0, we know that K = 0.

Hence, φ(t, x) ≥ φ(t) = e -u T ∞ σ 2 exp -1 σ 2 f ∞ (T -t)
≥ and the result follows. Now, we turn to a monotonicity result regarding Φ.

Proposition 3. ∀ψ 1 ≤ ψ 2 ∈ P 0 , Φ(ψ 1 ) ≥ Φ(ψ 2 )
Proof:

Let us introduce φ 1 = Φ(ψ 1 ) and φ 2 = Φ(ψ 2 ).

Let us introduce I(t) = 1 2 Ω ((φ 2 (t, x) -φ 1 (t, x)) + ) 2 dx. We have:

I (t) = Ω (∂ t φ 2 (t, x) -∂ t φ 1 (t, x))(φ 2 (t, x) -φ 1 (t, x)) + dx = Ω |∇φ 2 (t, x) -∇φ 1 (t, x)| 2 1 φ 2 (t,x)-φ 1 (t,x)≥0 - 1 σ 2 (f (x, φ 2 (t, x)ψ 2 (t, x))φ 2 (t, x) -f (x, φ 1 (t, x)ψ 1 (t, x))φ 1 (t, x)) (φ 2 (t, x) -φ 1 (t, x)) + dx ≥ 1 σ 2 Ω (f (x, φ 1 (t, x)ψ 1 (t, x))φ 1 (t, x) -f (x, φ 2 (t, x)ψ 2 (t, x))φ 2 (t, x)) (φ 2 (t, x) -φ 1 (t, x)) + dx ≥ 1 σ 2 Ω (f (x, φ 1 (t, x)ψ 2 (t, x))φ 1 (t, x) -f (x, φ 2 (t, x)ψ 2 (t, x))φ 2 (t, x)) (φ 2 (t, x) -φ 1 (t, x)) + dx ≥ 0
Hence, since I(T ) = 0 and I ≥ 0, we know that I = 0. Consequently, φ 1 ≥ φ 2 .

We now turn to the second equation (E ψ ) of the system (S).

Proposition 4. Let us fix > 0 and suppose that m 0 ∈ L 2 (Ω). ∀φ ∈ P , there is a unique weak solution ψ to the following equation (E ψ ):

∂ t ψ - σ 2 2 ∆ψ = 1 σ 2 f (x, φψ)ψ (E ψ )
with ∂ψ ∂ n = 0 on (0, T ) × ∂Ω and ψ(0, •) = m 0 (•) φ(0,•) .

Hence Ψ : φ ∈ P → ψ ∈ P is well defined.

Moreover, ∀φ ∈ P , ψ = Ψ(φ) ∈ P 0 .

Proof:

The proof of existence and uniqueness of a weak solution ψ ∈ P is the same as in Proposition 2. The only thing to notice is that the initial condition ψ(0, •) is in L 2 (Ω) because m 0 ∈ L 2 (Ω) and φ bounded from below by > 0. Now, to prove that ψ ≥ 0, let's introduce I(t) = 1 2 Ω (ψ(t, x) -) 2 dx, then:

I (t) = - Ω ∂ t ψ(t, x)ψ(t, x) -dx = - Ω -∇ψ(t, x) • ∇(ψ(t, x) -) + 1 σ 2 f (x, φ(t, x)ψ(t, x))ψ(t, x)ψ(t, x) -dx = - Ω |∇ψ(t, x)| 2 1 ψ(t,x)≤0 - 1 σ 2 f (x, φ(t, x)ψ(t, x))(ψ(t, x) -) 2 dx ≤ 0
Since I(0) = 0 and I ≥ 0, we know that I = 0. Hence, ψ is positive. Now, we turn to a monotonicity result regarding Ψ.

Proposition 5.

∀φ 1 ≤ φ 2 ∈ P , Ψ(φ 1 ) ≥ Ψ(φ 2 ) Proof: Let us introduce ψ 1 = Ψ(φ 1 ) and ψ 2 = Ψ(φ 2 ). Let us introduce I(t) = 1 2 Ω ((ψ 2 (t, x) -ψ 1 (t, x)) + ) 2 dx.
We have:

I (t) = Ω (∂ t ψ 2 (t, x) -∂ t ψ 1 (t, x))(ψ 2 (t, x) -ψ 1 (t, x)) + dx = 1 σ 2 (f (x, φ 2 (t, x)ψ 2 (t, x))ψ 2 (t, x) -f (x, φ 1 (t, x)ψ 1 (t, x))ψ 1 (t, x)) (ψ 2 (t, x) -ψ 1 (t, x)) + dx - Ω |∇ψ 2 (t, x) -∇ψ 1 (t, x)| 2 1 ψ 2 (t,x)-ψ 1 (t,x)≥0 dx ≤ 1 σ 2 Ω (f (x, φ 2 (t, x)ψ 2 (t, x))ψ 2 (t, x) -f (x, φ 1 (t, x)ψ 1 (t, x))ψ 1 (t, x)) (ψ 2 (t, x) -ψ 1 (t, x)) + dx ≤ 1 σ 2 Ω (f (x, φ 1 (t, x)ψ 2 (t, x))ψ 2 (t, x) -f (x, φ 1 (t, x)ψ 1 (t, x))ψ 1 (t, x)) (ψ 2 (t, x) -ψ 1 (t, x)) + dx ≤ 0 Now, I(0) = 1 2 Ω m 0 (x) 1 φ 2 (0,x) -1 φ 1 (0,x) + 2 dx = 0. Hence since I ≥ 0, we know that I = 0. Consequently, ψ 1 ≥ ψ 2 .
We will use these properties to design a constructive scheme for the couple (φ, ψ).

3 A constructive scheme to solve the system (S)

The scheme we consider involves two sequences (φ n+ 1 2 ) n and (ψ n ) n that are built using the following recursive equations:

ψ 0 = 0 ∀n ∈ N, φ n+ 1 2 = Φ(ψ n ) ∀n ∈ N, ψ n+1 = Ψ(φ n+ 1 2 ) Theorem 1. Suppose that u T ∈ L ∞ (Ω) and that m 0 ∈ L 2 (Ω).
Then, the above scheme has the following properties:

• (φ n+ 1 2
) n is a decreasing sequence of P where is as in Proposition 2. • (ψ n ) n is an increasing sequence of P 0 , bounded from above by Ψ( )

• (φ n+ 1 2 , ψ n ) n converges for almost every (t, x) ∈ (0, T ) × Ω, and in L 2 (0, T, L 2 (Ω)), towards a couple (φ, ψ).
• (φ, ψ) ∈ P × P 0 is a weak solution of (S).

Proof:

By immediate induction we get from Propositions 2 and 4 that the two sequences are well defined and in the appropriate spaces. Now, as far as monotonicity is concerned we have that ψ 1 = Ψ(φ 1 2 ) ≥ 0 = ψ 0 . Hence, if for a given n ∈ N we have ψ n+1 ≥ ψ n , then Proposition 3 gives:

φ n+ 3 2 = Φ(ψ n+1 ) ≤ Φ(ψ n ) = φ n+ 1 2
Applying now the function Ψ we get:

ψ n+2 = Ψ(φ n+ 3 2 ) ≥ Ψ(φ n+ 1 2 ) = ψ n+1
By induction we then have that (φ n+ 1 2 ) n is decreasing and (ψ n ) n is increasing.

Moreover, since φ n+ 1 2 ≥ , we have that

ψ n+1 = Ψ(φ n+ 1 2 ) ≤ Ψ( ).
Now this monotonic behavior allows to define two limit functions φ and ψ in L 2 (0, T, L 2 (Ω)) and the convergence is almost everywhere and in L 2 (0, T, L 2 (Ω)). Now, we want to show that (φ, ψ) is a weak solution of (S) and to this purpose we use the energy estimates of the parabolic equations.

We know that there exists a constant C > 0, that only depends on u T ∞ , σ and f ∞ such that:

∀n ∈ N, φ n+ 1 2 L 2 (0,T,H 1 (Ω)) + ∂ t φ n+ 1 2 L 2 (0,T,H -1 (Ω)) ≤ C
Hence, φ ∈ P and we can extract a subsequence (φ n + 1 2 ) n such that:

• φ n + 1 2 → φ in the L 2 (0, T, L 2 (Ω)
) sense and almost everywhere.

• ∇φ n + 1 2 ∇φ weakly in L 2 (0, T, L 2 (Ω)). • ∂ t φ n + 1 2 ∂ t φ weakly in L 2 (0, T, H -1 (Ω)). Now, for w ∈ L 2 (0, T, H 1 (Ω)) T 0 ∂ t φ n + 1 2 (t, •), w(t, •) H -1 (Ω),H 1 (Ω) dt - σ 2 2 T 0 Ω ∇φ n + 1 2 (t, x) • ∇w(t, x)dxdt = - 1 σ 2 T 0 Ω f (x, φ n + 1 2 (t, x)ψ n (t, x))φ n + 1 2 (t, x)w(t, x)dxdt
Using the weak convergence stated above, the continuity hypothesis on f and the dominated convergence theorem for the last term we get that ∀w ∈ L 2 (0, T, H 1 (Ω))

T 0 ∂ t φ(t, •), w(t, •) H -1 (Ω),H 1 (Ω) dt - σ 2 2 T 0 Ω ∇φ(t, x) • ∇w(t, x)dxdt = - 1 σ 2 T 0 Ω f (x, φ(t, x)ψ(t, x))φ(t, x)w(t, x)dxdt
Hence, for almost every t ∈ (0, T ) and ∀v ∈ H 1 (Ω),

∂ t φ(t, •), v H -1 (Ω),H 1 (Ω) - σ 2 2 Ω ∇φ(t, x) • ∇v(x)dx = - 1 σ 2 Ω f (x, φ(t, x)ψ(t, x))φ(t, x)v(x)dx
and the terminal condition is obviously satisfied.

Obviously, we also have φ ≥ . Now, turning to the second equation, the proof works the same. The only additional thing to notice is that

ψ(0, •) = lim n→∞ ψ n+1 (0, •) = lim n→∞ m 0 φ n+ 1 2 (0, •) = m 0 φ(0, •)
where the limits are in the L 2 (Ω) sense.

Hence (φ, ψ) is indeed a weak solution of (S).

We exhibited a way to build a solution to the system (S). However, it should be noticed that the two functions φ and ψ have been introduced to transform the initial (MFG) system into two simpler partial differential equations. The functions we are interested in, as far as the mean field game is concerned, are rather u and m.

To well understand the nature of the change of variables and of the constructive scheme, let us introduce the sequence (m n+1 ) n where m n+1 = φ n+ 1 2 ψ n+1 . From Theorem 1, we know that (m n+1 ) n converges almost everywhere and in L 1 toward the function m = φψ for which we have conservation of mass along the trajectory. However, this property is not true for m n+1 as it is stated in the following proposition: Proposition 6. Let us consider n ∈ N and let's denote M n+1 (t) = Ω m n+1 (t, x)dx the total mass of m n+1 at date t. Then, there may be a loss of mass along the trajectory in the sense that:

d dt M n+1 (t) = Ω ψ n+1 (t, x)φ n+ 1 2 (t, x) f x, ψ n+1 (t, x)φ n+ 1 2 (t, x) -f x, ψ n (t, x)φ n+ 1 2 (t, x) ≤ 0 Proof:
From the regularity obtained above, we can write:

d dt M n+1 (t) = ∂ t φ n+ 1 2 (t, •), ψ n+1 (t, •) H -1 (Ω),H 1 (Ω) + ∂ t ψ n+1 (t, •), φ n+ 1 2 (t, •) H -1 (Ω),H 1 (Ω) = σ 2 2 Ω ∇φ n+ 1 2 (t, x) • ∇ψ n+1 (t, x)dx - 1 σ 2 Ω φ n+ 1 2 (t, x)ψ n+1 (t, x)f x, φ n+ 1 2 (t, x)ψ n (t, x) dx - σ 2 2 Ω ∇φ n+ 1 2 (t, x) • ∇ψ n+1 (t, x)dx + 1 σ 2 Ω φ n+ 1 2 (t, x)ψ n+1 (t, x)f x, φ n+ 1 2 (t, x)ψ n+1 (t, x) dx = Ω ψ n+1 (t, x)φ n+ 1 2 (t, x) f x, ψ n+1 (t, x)φ n+ 1 2 (t, x) -f x, ψ n (t, x)φ n+ 1 2 (t, x) ≤ 0
This property shows that the constructive scheme is rather original since it basically consists in building a probability distribution function using a sequence of functions in L 1 that only has the right total mass asymptotically. For the same reason, although there will be no systematic loss of mass, the numerical scheme we present in the next section consists in approximating probability distribution functions without taking care to mass conservation.

A numerical scheme

The constructive scheme of Theorem 1 allows to design a monotonic numerical scheme to solve the system (S) and in turn to approximate u and m.

For the sake of simplicity, we present the scheme with Ω = (0, 1) but the results would be similar in higher dimensions. We also assume that u T and m 0 are bounded functions.

Let us consider a uniform subdivision (t 0 , . . . , t I ) of [0, T ] where t i = i∆t for i ∈ I = {0, . . . , I} and a uniform subdivision (x 0 , . . . , x J ) of the interval [0, 1] where x j = j∆x for j ∈ J = {0, . . . , J} .

We consider a recursive sequence of finite difference schemes in which ψ n i,j and φ

n+ 1 2 i,j
stand for the approximations of ψ n and φ n+ 1 2 at point (t i , x j ). For convenience, we also define

ψ n i,-1 = ψ n i,0
and ψ n i,J+1 = ψ n i,J to take account of Neumann conditions at the border.

The numerical scheme is as follows:

ψ 0 i,j = 0, (i, j) ∈ I × J and for n ≥ 0:

φ n+ 1 2 i+1,j -φ n+ 1 2 i,j ∆t + σ 2 2 φ n+ 1 2 i,j+1 -2 φ n+ 1 2 i,j + φ n+ 1 2 i,j-1 (∆x) 2 = - 1 σ 2 f (x j , φ n+ 1 2 i,j ψ n i,j ) φ n+ 1 2 i,j , i ∈ I -{I}, j ∈ J φ n+ 1 2 I,j = exp u T (x j ) σ 2 , j ∈ J ψ n+1 i+1,j -ψ n+1 i,j ∆t - σ 2 2 ψ n+1 i+1,j+1 -2 ψ n+1 i+1,j + ψ n+1 i+1,j-1 (∆x) 2 = 1 σ 2 f (x j , φ n+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j , i ∈ I -{I}, j ∈ J ψ n+1 0,j = m 0 (x j ) φ n+ 1 2 0,j , j ∈ J
To study this scheme, let's introduce M = M I+1,J+1 (R) and

M = {(m i,j ) i∈I,j∈J ∈ M, ∀(i, j) ∈ I × J , m i,j ≥ }
Let us start with the discrete counterpart of Proposition 2.

Proposition 7. ∀ ψ ∈ M 0 , there is a unique solution φ ∈ M to the following equation:

φ i+1,j -φ i,j ∆t + σ 2 2 φ i,j+1 -2 φ i,j + φ i,j-1 (∆x) 2 = - 1 σ 2 f (x j , φ i,j ψ i,j ) φ i,j , i ∈ I -{I}, j ∈ J with φ I,j = exp u T (x j ) σ 2
, ∀j ∈ J and the conventions φ i,-1 = φ i,0 , φ i,J+1 = φ i,J .

Hence Φ :

ψ ∈ M 0 → φ ∈ M is well defined. Moreover, ∀ ψ ∈ M 0 , φ = Φ( ψ) ∈ M for = exp -1 σ 2 ( u T ∞ + f ∞ T ) Proof:
Let us consider ψ ∈ M 0 .

Existence and positiveness of a solution φ:

Let us introduce F ψ : ϕ ∈ M → φ solution of the following equations:

φ i+1,j -φ i,j ∆t + σ 2 2 φ i,j+1 -2 φ i,j + φ i,j-1 (∆x) 2 = - 1 σ 2 f (x j , ϕ i,j ψ i,j ) φ i,j , i ∈ I -{I}, j ∈ J with φ I,j = exp u T (x j ) σ 2
, ∀j ∈ J and the conventions φ i,-1 = φ i,0 , φ i,J+1 = φ i,J .

Such a matrix φ exists by straightforward backward induction. We indeed know φ I,j , ∀j ∈ J , and if we know φ i+1,j , ∀j ∈ J for some i ∈ I -{I}, then we have:

A    φ i,0
. . .

φ i,J    =    φ i+1,0 . . . φ i+1,J   
where:

A =           1 + ∆tβ i,0 -σ 2 ∆t 2(∆x) 2 0 • • • 0 -σ 2 ∆t 2(∆x) 2 1 + ∆tβ i,1 -σ 2 ∆t 2(∆x) 2 . . . . . . 0 . . . . . . . . . 0 . . . . . . -σ 2 ∆t 2(∆x) 2 1 + ∆tβ i,J-1 -σ 2 ∆t 2(∆x) 2 0 • • • 0 -σ 2 ∆t 2(∆x) 2 1 + ∆tβ i,J          
where

β i,j = σ 2 (∆x) 2 -1 σ 2 f (x j , ϕ i,j ψ i,j ), ∀i ∈ I -{I}, ∀j ∈ J -{0, J}, β i,0 = σ 2 2(∆x) 2 -1 σ 2 f (x 0 , ϕ i,0 ψ i,0 ) and β i,J = σ 2 2(∆x) 2 -1 σ 2 f (x J , ϕ i,J ψ i,J ), ∀i ∈ I -{I}.
Since A is an M -matrix, φ i,j is well defined ∀j ∈ J , and by immediate induction φ is well defined. Now, F ψ is well-defined and continuous since f is continuous.

To prove that there is a fixed point to F ψ , we will use Brouwer's theorem. The only thing that needs to be proved is that ∃C > 0, ∀ ϕ ∈ M, max (i,j)∈I×J F ψ ( ϕ) i,j ≤ C. More precisely, we are going to prove that if φ = F ψ ( ϕ) then:

∀(i, j) ∈ I × J , 0 ≤ φ i,j ≤ max j∈J φ I,j ≤ exp u T σ 2 ∞
Let us suppose first that φ = min (i,j)∈I×J φ i,j < 0.

Then, let's consider (i, j) ∈ (i , j ) ∈ I × J | φ i,j = φ, ∀(i , j ) ∈ I × J , φ i ,j = φ ⇒ i ≤ i . We know that i = I and hence:

0 < φ i+1,j -φ i,j ∆t + σ 2 2 φ i,j+1 -2 φ i,j + φ i,j-1 (∆x) 2 = - 1 σ 2 f (x j , ϕ i,j ψ i,j ) φ i,j
Since f ≤ 0, this gives φ i,j = φ ≥ 0 contrary to our hypothesis. Hence, ∀(i, j) ∈ I × J , φ i,j ≥ 0.

Similarly let's suppose that φ = max (i,j)∈I×J φ i,j > max j∈J φ I,j . Then, let's consider a couple

(i, j) ∈ (i , j ) ∈ I × J | φ i,j = φ, ∀(i , j ) ∈ I × J , φ i ,j = φ ⇒ i ≤ i .
We have that i = I and hence:

0 > φ i+1,j -φ i,j ∆t + σ 2 2 φ i,j+1 -2 φ i,j + φ i,j-1 (∆x) 2 = - 1 σ 2 f (x j , ϕ i,j ψ i,j ) φ i,j
But then, φ i,j < 0 and this is not possible. Hence, we indeed have that the maximum is attained for i = I, namely max (i,j)∈I×J φ i,j = max j∈J φ I,j . Now, from Brouwer's theorem, we know that there exists φ ∈ M 0 , such that F ψ ( φ) = φ.

Uniqueness of φ:

Let us suppose that there are two fixed points φ and φ to F ψ and imagine for instance that ∃(i, j) ∈ I × J such that φ i,j > φ i,j . Then, let's consider δ = φ -φ , δ = max (i,j)∈I×J δ i,j and (i, j) ∈ (i , j ) ∈ I × J | δ i,j = δ, ∀(i , j ) ∈ I × J , φ i ,j = δ ⇒ i ≤ i . We know that i = I and hence:

0 > δ i+1,j -δ i,j ∆t + σ 2 2 δ i,j+1 -2 δ i,j + δ i,j-1 (∆x) 2 = - 1 σ 2 f (x j , φ i,j ψ i,j ) φ i,j -f (x j , φ i,j ψ i,j ) φ i,j
But since φ i,j > φ i,j ≥ 0 as above, the right hand side of the above equation must be greater than 0, in contradiction with the above inequality. Hence, φ ≤ φ and, by symmetry, φ = φ . As a consequence Φ : ψ ∈ M 0 → φ ∈ M is well defined.

Lower bound for φ:

Let us consider φ i,j = exp -u T ∞ σ 2 1 1+∆t f ∞ σ 2
I-i solution of the following problem:

φ i+1,j -φ i,j ∆t + σ 2 2 φ i,j+1 -2 φ i,j + φ i,j-1 (∆x) 2 = 1 σ 2 f ∞ φ i,j , i ∈ I -{I}, j ∈ J with φ I,j = exp -u T ∞ σ 2
, ∀j ∈ J and the conventions φ i,-1 = φ i,0 , φ i,J+1 = φ i,J . Now, let's write δ = φ -φ and let's suppose that δ = max (i,j)∈I×J δ i,j > 0. Then, let's consider

(i, j) ∈ (i , j ) ∈ I × J | δ i,j = δ, ∀(i , j ) ∈ I × J , δ i ,j = δ ⇒ i ≤ i .
By construction, we know that i = I and hence:

0 > δ i+1,j -δ i,j ∆t + σ 2 2 δ i,j+1 -2 δ i,j + δ i,j-1 (∆x) 2 = 1 σ 2 f ∞ φ i,j + f (x j , φ i,j ψ i,j ) φ i,j
But then, φ i,j < φ i,j in contradiction with the hypothesis.

In conclusion, we have ∀(i, j) ∈ I × J :

φ i,j ≥ φ i,j ≥ e- u T ∞ σ 2 1 1 + ∆t f ∞ σ 2 I-i ≥ e -u T ∞ σ 2 1 1 + ∆t f ∞ σ 2 T ∆t

≥

Let us now turn to a monotonicity property of the function Φ.

Proposition 8.

∀ ψ 1 ≤ ψ 2 ∈ M 0 , Φ( ψ 1 ) ≥ Φ( ψ 2 )
Proof:

Let us consider ψ 1 ≤ ψ 2 ∈ M 0 and let's denote φ 1 = Φ( ψ 1 ) and φ 2 = Φ( ψ 2 ). Now, let's write δ = φ 2 -φ 1 and let's suppose that δ = max (i,j)∈I×J δ i,j > 0. Then, let's consider (i, j) ∈ (i , j ) ∈ I × J | δ i,j = δ, ∀(i , j ) ∈ I × J , δ i ,j = δ ⇒ i ≤ i . We know that i = I and hence:

0 > δ i+1,j -δ i,j ∆t + σ 2 2 δ i,j+1 -2 δ i,j + δ i,j-1 (∆x) 2 = - 1 σ 2 f (x j , φ 2,i,j ψ 2,i,j ) φ 2,i,j -f (x j , φ 1,i,j ψ 1,i,j ) φ 1,i,j
But then, f (x j , φ 2,i,j ψ 2,i,j ) φ 2,i,j -f (x j , φ 1,i,j ψ 1,i,j ) φ 1,i,j > 0 in contradiction with the hypothesis for the chosen (i, j).

Hence, φ 1 ≥ φ 2 .
Let us now turn to the second equation. We have exactly the same result as above with the same proof: Proposition 9. ∀ φ ∈ M , there is a unique solution ψ ∈ M to the following equation:

ψ i+1,j -ψ i,j ∆t - σ 2 2 ψ i+1,j+1 -2 ψ i+1,j + ψ i+1,j-1 (∆x) 2 = 1 σ 2 f (x j , φ i+1,j ψ i+1,j ) ψ i+1,j , 0 ≤ i ≤ I-1, 0 ≤ j ≤ J with ψ 0,j = m 0 (x j ) φ 0,j
, ∀j ∈ J and the conventions ψ i,-1 = ψ i,0 , ψ i,J+1 = ψ i,J .

Hence

Ψ : φ ∈ M → ψ ∈ M is well defined. Moreover, ∀ φ ∈ M , ψ = Ψ( φ) ∈ M 0 .
Similarly, the same monotonicity result can be proved:

Proposition 10. ∀ φ 1 ≤ φ 2 ∈ M , Ψ( φ 1 ) ≥ Ψ( φ 2 )
Now, using the functions introduced in the above propositions, we can write the numerical scheme in a more compact fashion:

ψ 0 = 0 ∀n ∈ N, φ n+ 1 2 = Φ ψ n ∀n ∈ N, ψ n+1 = Ψ φ n+ 1 2
Using the monotonicity properties, we get that the sequences ψ n n and φ n+ 1 2 n converge monotonically. More precisely we have that: Proposition 11. The numerical scheme verifies the following properties:

• φ n+ 1 2 n
is a decreasing sequence of M where is as in Proposition 7.

• ψ n n is an increasing sequence of M 0 , bounded from above by m 0 ∞ .

• φ n+ 1 2 , ψ n n converges towards a couple ( φ, ψ) ∈ M × M 0 that verifies:

φ i+1,j -φ i,j ∆t + σ 2 2 φ i,j+1 -2 φ i,j + φ i,j-1 (∆x) 2 = - 1 σ 2 f (x j , φ i,j ψ i,j ) φ i,j , i ∈ I -{I}, j ∈ J φ I,j = exp u T (x j ) σ 2 , ∀j ∈ J ψ i+1,j -ψ i,j ∆t - σ 2 2 ψ i+1,j+1 -2 ψ i+1,j + ψ i+1,j-1 (∆x) 2 = 1 σ 2 f (x j , φ i+1,j ψ i+1,j ) ψ i+1,j , i ∈ I-{I}, j ∈ J ψ 0,j = m 0 (x j ) φ 0,j , ∀j ∈ J Proof:
Let us prove by induction that

∀n ∈ N, 0 ≤ ψ n ≤ ψ n+1 ≤ Ψ( ), ≤ φ n+ 3 2 ≤ φ n+ 1 2
First, for n = 0 we have that:

ψ 1 = Ψ( Φ( ψ 0 )) ≥ 0 = ψ 0 Consequently, φ 1 2 = Φ( ψ 0 ) ≥ Φ( ψ 1 ) = φ 3 2
and because of Proposition 7, we know that φ

1 2 ≥ φ 3 2 ≥ . Hence, ψ 1 = Ψ( ψ 1 2 ) ≤ Ψ( ).
Now, let's suppose that the result is true for some n ∈ N, then:

ψ n ≤ ψ n+1 ⇒ φ n+ 1 2 = Φ( ψ n ) ≥ Φ( ψ n+1 ) = φ n+ 3 2 ⇒ ψ n+1 = Ψ( φ n+ 1 2 ) ≤ Ψ( φ n+ 3 2 ) = ψ n+2
Now, by Proposition 7, we know that φ n+ 3 2 ≥ and thus ψ n+2 = Ψ( ψ n+ 3 2 ) ≤ Ψ( ) and with the same discrete maximum principle as above, we have that Ψ( ) is bounded by m 0 ∞ . Consequently, ∀(i, j) ∈ I × J , ( φ n+ 1 2 i,j ) n is a decreasing sequence bounded from below by that thus converges towards a value denoted φ i,j . Similarly, ∀(i, j) ∈ I × J , ( ψ n i,j ) n is an increasing sequence bounded from above that thus converges towards a value denoted ψ i,j . The resulting couple ( φ, ψ) ∈ M × M 0 straightforwardly verifies the equations of the proposition. Now, to study convergence, let's introduce a norm • on the sets M by:

∀m = (m i,j ) i,j ∈ M, m 2 = sup 0≤i≤I 1 J + 1 J j=0 m 2 i,j
Let us also define the consistency errors ηn+ 1 2 and ηn+1 associated to the equations that define respectively φ n+ 1 2 and ψ n+1 :

ηn+ 1 2 i,j = φn+ 1 2 i+1,j - φn+ 1 2 i,j ∆t + σ 2 2 φn+ 1 2 i,j+1 -2 φn+ 1 2 i,j + φn+ 1 2 i,j-1 (∆x) 2 + 1 σ 2 f (x j , φn+ 1 2 i,j ψn i,j ) φn+ 1 2 i,j ηn+1 i,j = ψn+1 i+1,j -ψn+1 i,j ∆t - σ 2 2 ψn+1 i+1,j+1 -2 ψn+1 i+1,j + ψn+1 i+1,j-1 (∆x) 2 - 1 σ 2 f (x j , φn+ 1 2 i+1,j ψn+1 i+1,j ) ψn+1 i+1,j
where φn+ 1 2 i,j = φ n+ 1 2 (t i , x j ) and ψn+1 i,j = ψ n+1 (t i , x j ).

We are now ready to enounce a theorem that gives recursive stability bounds for the scheme.

Theorem 2. Let us suppose, in addition to the hypotheses made above, that f is uniformly Lipschitz with respect to its second variable, i.e.

∃K, ∀x ∈ (0, 1),

∀ξ 1 , ξ 2 ∈ R + , |f (x, ξ 2 ) -f (x, ξ 1 )| ≤ K|ξ 2 -ξ 1 | Then, ∀ν > 0, ∃C > 0, ∀I, J ∈ N such that 1 ∆t > 1 + K σ 2 max e u T σ 2 2 ∞ , m 0 2 ∞ 2 + ν, we have ∀n ∈ N: φ n+ 1 2 -φn+ 1 2 2 ≤ C ψ n -ψn 2 + C ηn+ 1 2 2 ψ n+1 -ψn+1 2 ≤ C φ n+ 1 2 -φn+ 1 2 2 + C ηn+1 2 Proof:
Let us denote for n ∈ N, δ

n+ 1 2 i,j = φ n+ 1 2 i,j - φn+ 1 2 i,j
and δ n i,j = ψ n i,j -ψn i,j .

By definition we have for i ∈ I -{I}, j ∈ J :

δ n+ 1 2 i+1,j -δ n+ 1 2 i,j ∆t + σ 2 2 δ n+ 1 2 i,j+1 -2δ n+ 1 2 i,j + δ n+ 1 2 i,j-1 (∆x) 2 = - 1 σ 2 f (x j , φ n+ 1 2 i,j ψ n i,j ) φ n+ 1 2 i,j -f (x j , φn+ 1 2 i,j ψn i,j ) φn+ 1 2 i,j -η n+ 1 2 i,j
Hence:

δ n+ 1 2 i,j 2 = δ n+ 1 2 i+1,j δ n+ 1 2 i,j + σ 2 ∆t 2(∆x) 2 δ n+ 1 2 i,j δ n+ 1 2 i,j+1 -2δ n+ 1 2 i,j + δ n+ 1 2 i,j-1 + ∆t σ 2 δ n+ 1 2 i,j f (x j , φ n+ 1 2 i,j ψ n i,j ) φ n+ 1 2 i,j -f (x j , φn+ 1 2 i,j ψn i,j ) φn+ 1 2 i,j + δ n+ 1 2 i,j ηn+ 1 2 i,j ∆t ⇒ 1 2 δ n+ 1 2 i,j 2 ≤ 1 2 δ n+ 1 2 i+1,j 2 + σ 2 ∆t 2(∆x) 2 δ n+ 1 2 i,j δ n+ 1 2 i,j+1 -2δ n+ 1 2 i,j + δ n+ 1 2 i,j-1 + ∆t σ 2 δ n+ 1 2 i,j f (x j , φ n+ 1 2 i,j ψ n i,j ) φ n+ 1 2 i,j -f (x j , φn+ 1 2 i,j ψn i,j ) φn+ 1 2 i,j + ∆t 2 δ n+ 1 2 i,j 2 + ∆t 2 ηn+ 1 2 i,j 2 Now, δ n+ 1 2 i,j f (x j , φ n+ 1 2 i,j ψ n i,j ) φ n+ 1 2 i,j -f (x j , φn+ 1 2 i,j ψn i,j ) φn+ 1 2 i,j = δ n+ 1 2 i,j f (x j , φ n+ 1 2 i,j ψ n i,j ) φ n+ 1 2 i,j -f (x j , φ n+ 1 2 i,j ψn i,j ) φ n+ 1 2 i,j + φ n+ 1 2 i,j - φn+ 1 2 i,j f (x j , φ n+ 1 2 i,j ψn i,j ) φ n+ 1 2 i,j -f (x j , φn+ 1 2 i,j ψn i,j ) φn+ 1 2 i,j ≤ δ n+ 1 2 i,j f (x j , φ n+ 1 2 i,j ψ n i,j ) φ n+ 1 2 i,j -f (x j , φ n+ 1 2 i,j ψn i,j ) φ n+ 1 2 i,j ≤ K δ n+ 1 2 i,j φ n+ 1 2 i,j φ n+ 1 2 i,j ψ n i,j -φ n+ 1 2 i,j ψn i,j ≤ K δ n+ 1 2 i,j φ n+ 1 2 i,j 2 δ n i,j ≤ K 2 φ n+ 1 2 2 ∞ δ n+ 1 2 i,j 2 + δ n i,j 2 ≤ K 2 φ 1 2 2 ∞ δ n+ 1 2 i,j 2 + δ n i,j 2
Thus:

J j=0 δ n+ 1 2 i,j 2
1 -∆t - ∆t σ 2 K φ 1 2 2 ∞ ≤ J j=0 δ n+ 1 2 i+1,j 2 + σ 2 ∆t (∆x) 2 J j=0 δ n+ 1 2 i,j δ n+ 1 2 i,j+1 -2δ n+ 1 2 i,j + δ n+ 1 2 i,j-1 + ∆t σ 2 K φ 1 2 2 ∞ J j=0 δ n i,j 2 + ∆t J j=0 ηn+ 1 2 i,j 2 ≤ J j=0 δ n+ 1 2 i+1,j 2 - σ 2 ∆t (∆x) 2 J-1 j=0 δ n+ 1 2 i,j+1 -δ n+ 1 2 i,j 2 + ∆t σ 2 K φ 1 2 2 ∞ (J + 1) δ n 2 + ∆t(J + 1) ηn+ 1 2 2 ≤ J j=0 δ n+ 1 2 i+1,j 2 + ∆t σ 2 K φ 1 2 2 ∞ (J + 1) δ n 2 + ∆t(J + 1) ηn+ 1 2 2
We know that φ

1 2 ∞ ≤ e u T σ 2 ∞ . Hence, α(∆t) = 1 -∆t -∆t σ 2 K φ 1 2 2 ∞
lies in (0, 1) by hypothesis.

This gives:

J j=0 δ n+ 1 2 i,j 2 ≤ 1 α(∆t) J j=0 δ n+ 1 2 i+1,j 2 + ∆t α(∆t) (J + 1) 1 σ 2 K φ 1 2 2 ∞ δ n 2 + ηn+ 1 2 2
By immediate induction we have:

J j=0 δ n+ 1 2 i,j 2 ≤ 1 α(∆t) I-i J j=0   δ n+ 1 2 I,j =0    2 + ∆t I-1 l=i 1 α(∆t) l-i (J + 1) 1 σ 2 K φ 1 2 2 ∞ δ n 2 + ηn+ 1 2 2 ≤ 1 α(∆t) I T (J + 1) 1 σ 2 K φ 1 2 2 ∞ δ n 2 + ηn+ 1 2 2 ≤ 1 α(∆t) T ∆t T (J + 1) 1 σ 2 K φ 1 2 2 ∞ δ n 2 + ηn+ 1 2 2
Thus:

δ n+ 1 2 2 ≤ 1 α(∆t) T ∆t T 1 σ 2 K φ 1 2 2 ∞ δ n 2 + ηn+ 1 2 2 Now, lim ∆t→0 1 α(∆t) T ∆t = exp 1 + 1 σ 2 K φ 1 2 2 ∞
T and hence there exists a constant C 1 such that:

δ n+ 1 2 2 ≤ C 1 δ n 2 + C 1 ηn+ 1 2 2
This is the first part of the Theorem. For the second part, we need to control the difference arising from the initial condition. By definition we have for i ∈ I -{I}, j ∈ J :

δ n+1 i+1,j -δ n+1 i,j ∆t - σ 2 2 δ n+1 i+1,j+1 -2δ n+1 i+1,j + δ n+1 i+1,j-1 (∆x) 2 = 1 σ 2 f (x j , φ n+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j -f (x j , φn+ 1 2 i+1,j ψn+1 i+1,j ) ψn+1 i+1,j -ηn+1 i+1,j
Hence:

δ n+1 i+1,j 2 = δ n+1 i+1,j δ n+1 i,j + σ 2 ∆t 2(∆x) 2 δ n+1 i+1,j δ n+1 i+1,j+1 -2δ n+1 i+1,j + δ n+1 i+1,j-1 + ∆t σ 2 δ n+1 i+1,j f (x j , φ n+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j -f (x j , φn+ 1 2 i+1,j ψn+1 i+1,j ) ψn+1 i+1,j -δ n+1 i+1,j ηn+1 i+1,j ∆t ⇒ 1 2 δ n+1 i+1,j 2 ≤ 1 2 δ n+1 i+1,j 2 + σ 2 ∆t 2(∆x) 2 δ n+1 i+1,j δ n+1 i+1,j+1 -2δ n+1 i+1,j + δ n+1 i+1,j-1 + ∆t σ 2 δ n+1 i+1,j f (x j , φ n+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j -f (x j , φn+ 1 2 i+1,j ψn+1 i+1,j ) ψn+1 i+1,j + ∆t 2 δ n+1 i+1,j 2 + ∆t 2 ηn+1 i+1,j 2
Now,

δ n+1 i+1,j f (x j , φ n+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j -f (x j , φn+ 1 2 i+1,j ψn+1 i+1,j ) ψn+1 i+1,j = δ n+1 i+1,j f (x j , φ n+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j -f (x j , φn+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j + ψ n+1 i,j -ψn+1 i,j f (x j , φn+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j -f (x j , φn+ 1 2 i+1,j ψn+1 i+1,j ) ψn+1 i+1,j ≤ δ n+1 i+1,j f (x j , φ n+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j -f (x j , φn+ 1 2 i+1,j ψ n+1 i+1,j ) ψ n+1 i+1,j ≤ K δ n+1 i+1,j ψ n+1 i+1,j φ n+ 1 2 i+1,j ψ n+1 i+1,j - φn+ 1 2 i+1,j ψ n+1 i+1,j ≤ K δ n+1 i+1,j ψ n+1 i+1,j 2 δ n+ 1 2 i+1,j ≤ K 2 ψ n+1 2 ∞ δ n+ 1 2 i+1,j 2 + δ n+1 i+1,j 2 ≤ K 2 ψ 2 ∞ δ n+ 1 2 i+1,j 2 + δ n+1 i+1,j 2
Thus:

J j=0 δ n+1 i+1,j 2 1 -∆t - ∆t σ 2 K ψ 2 ∞ ≤ J j=0 δ n+1 i,j 2 + σ 2 ∆t (∆x) 2 J j=0 δ n+1 i+1,j δ n+1 i+1,j+1 -2δ n+1 i+1,j + δ n+1 i+1,j-1 + ∆t σ 2 K ψ 2 ∞ J j=0 δ n+ 1 2 i+1,j 2 + ∆t J j=0 ηn+1 i+1,j 2 ≤ J j=0 δ n+1 i,j 2 - σ 2 ∆t (∆x) 2 J-1 j=0 δ n+1 i+1,j+1 -δ n+1 i+1,j 2 + ∆t σ 2 K ψ 2 ∞ (J + 1) δ n+ 1 2 2 + ∆t(J + 1) ηn+1 2 ≤ J j=0 δ n+1 i,j 2 + ∆t σ 2 K ψ 2 ∞ (J + 1) δ n+ 1 2 2 + ∆t(J + 1) ηn+1 2 Let us introduce β(∆t) = 1 -∆t -∆t σ 2 K ψ 2 ∞
. By hypothesis, β lies in (0, 1).

This gives:

β(∆t) J j=0 δ n+1 i+1,j 2 ≤ J j=0 δ n+1 i,j 2 + ∆t(J + 1) 1 σ 2 K ψ 2 ∞ δ n+ 1 2 2 + ηn+1 2
By immediate induction we have:

(β(∆t)) i+1 J j=0 δ n+1 i+1,j 2 ≤ J j=0 δ n+1 0,j 2 + ∆t i l=0 (β(∆t)) l (J + 1) 1 σ 2 K ψ 2 ∞ δ n+ 1 2 2 + ηn+1 2
Hence:

J j=0 δ n+1 i+1,j 2 ≤ 1 β(∆t) I J j=0 δ n+1 0,j 2 + 1 β(∆t) I T (J + 1) 1 σ 2 K ψ 2 ∞ δ n+ 1 2 2 + ηn+1 2
Now,

δ n+1 0,j = m 0 (x j ) φ n+ 1 2 0,j - m 0 (x j ) φn+ 1 2 0,j
Consequently:

|δ n+1 0,j | ≤ m 0 ∞ 1 φ n+ 1 2 0,j - 1 φn+ 1 2 0,j ≤ m 0 ∞ 2 δ n+ 1 2 0,j
We end up with:

J j=0 δ n+1 i+1,j 2 ≤ 1 β(∆t) I m 0 ∞ 2 2 (J + 1) δ n+ 1 2 2 + 1 β(∆t) I T (J + 1) 1 σ 2 K ψ 2 ∞ δ n+ 1 2 2 + ηn+1 2 ⇒ δ n+1 2 ≤ 1 β(∆t) I m 0 ∞ 2 2 (J + 1) δ n+ 1 2 2 + 1 β(∆t) I T (J + 1) 1 σ 2 K ψ 2 ∞ δ n+ 1 2 2 + ηn+1 2
Now, as above, lim ∆t→0 1 β(∆t)

T ∆t = exp 1 + 1 σ 2 K ψ 2 ∞
T and hence there exists a constant C 2 such that:

δ n+1 2 ≤ C 2 δ n+ 1 2 2 + C 2 ηn+1 2
Eventually, C = max(C 1 , C 2 ) provides the bounds given in the theorem

These bounds allow to prove, under regularity assumptions, that the scheme is convergent. Let us start with the convergence of the schemes defining respectively φ n+ 1 2 and ψ n .

Theorem 3. Let us suppose, in addition to the hypotheses of Theorem 2 that u T , m 0 and f are so that ∀n ∈ N,

φ n+ 1 2 , ψ n ∈ C 1,2 ([0, T ] × [0, 1]). Then ∀n ∈ N: lim ∆t,∆x→0 φ n+ 1 2 -φn+ 1 2 = 0 lim ∆t,∆x→0 ψ n+1 -ψn+1 = 0 Proof:
Let us fix n ≥ 0. We have by immediate induction from Theorem 2 that, as soon as ∆t is small enough, that there exists a constant C independent of ∆t and ∆x so that:

φ n+ 1 2 -φn+ 1 2 2 ≤ n k=0 C 2k+1 ηn+ 1 2 -k 2 + C 2k+2 ηn-k 2 and ψ n+1 -ψn+1 2 ≤ n k=0 C 2k+1 ηn+1-k 2 + C 2k+2 ηn+ 1 2 -k 2
Hence, because our assumptions imply that the consistency errors tend to 0 as ∆t and ∆x tend to 0, we have that:

lim ∆t,∆x→0 φ n+ 1 2 -φn+ 1 2 = 0
and lim ∆t,∆x→0

ψ n+1 -φn+1 = 0
For the global convergence of the numerical scheme we have the following theorem: We have seen that (φ n+ 1 2 ) n converges monotonically towards φ. Hence, we know from Dini's theorem that the convergence is in fact uniform. Consequently:

∀ε > 0, ∃n 0 , ∀n ≥ n 0 , ∀∆t, ∆x > 0, φn+ 1 2 -φ ≤ φn+ 1 2 -φ ∞ ≤ ε
Hence, ∀n ≥ n 0 , we know from the preceding result that: lim sup ∆t,∆x→0

φ n+ 1 2 -φ ≤ lim sup ∆t,∆x→0 φ n+ 1 2 -φn+ 1 2 + φn+ 1 2 -φ ≤ ε
Thus, we get that:

lim n→∞ lim sup ∆t,∆x→0 φ n+ 1 2 -φ = 0
The same proof holds regarding ψ.

We have proved that the scheme was convergent and we provided a bound on ∆t that guarantees stability of the scheme. Although this bound is certainly not the best one, a remark must be done regarding how it depends on parameters and especially on σ. It's indeed important to notice that the bound in Theorem 2 is

1 + K σ 2 max e u T σ 2 2 ∞ , m 0 2 ∞ 2 = 1 + K σ 2 max e u T σ 2 2 ∞ , m 0 2 ∞ exp 2 σ 2 ( u T ∞ + f ∞ T )
and this bound tends very rapidly to +∞ as σ tends to 0. As a consequence, the numerical scheme may not be useful for small values of σ.

We will empirically tackle the question of the influence of σ in the next section but before turning to the numerical experiments, let's highlight another important property of this numerical scheme. In most cases in which a probability distribution function is involved, the numerical scheme is built so that mass is preserved along the trajectory. Here, as we have seen above, since m is only reconstructed at the end from ψ and φ, there is nothing like mass conservation. However, a difference arises with Proposition 6 since there is theoretically no systematic loss of mass for a fixed n as t goes from 0 to T . The evolution of total mass is given by the following proposition: Proposition 12. Let us introduce m n+1 = φ n+ 1 2 ψ n+1 . The difference in total mass between two successive times at step n + 1 is:

1 J + 1 J j=0 m n+1 i+1,j - 1 J + 1 J j=0 m n+1 i,j = ∆t σ 2 1 J + 1 J j=0 ψ n+1 i+1,j φ n+ 1 2 i,j f ψ n+1 i+1,j φ n+ 1 2 i+1,j -f ψ n i,j φ n+ 1 2 i,j
Then, under the hypotheses of Theorem 4, we have that:

∀i ∈ I, lim n→∞ lim sup ∆t,∆x→0 1 J + 1 J j=0 m n+1 i,j = Ω m 0 (x)dx
Proof:

1 J + 1 J j=0 m n+1 i+1,j - 1 J + 1 J j=0 m n+1 i,j = 1 J + 1 J j=0 ψ n+1 i+1,j φ n+ 1 2 i+1,j - 1 J + 1 J j=0 ψ n+1 i,j φ n+ 1 2 i,j = 1 J + 1 J j=0 ψ n+1 i+1,j φ n+ 1 2 i+1,j -φ n+ 1 2 i,j + 1 J + 1 J j=0 ψ n+1 i+1,j -ψ n+1 i,j φ n+ 1 2 i,j = - ∆t J + 1   σ 2 2(∆x) 2 J j=0 ψ n+1 i+1,j φ n+ 1 2 i,j+1 -2 φ n+ 1 2 i,j + φ n+ 1 2 i,j-1 + 1 σ 2 J j=0 ψ n+1 i+1,j φ n+ 1 2 i,j f φ n+ 1 2 i,j ψ n i,j   + ∆t J + 1   σ 2 2(∆x) 2 J j=0 φ n+ 1 2 i,j ψ n+1 i+1,j+1 -2 ψ n+1 i+1,j + ψ n+1 i+1,j-1 + 1 σ 2 J j=0 ψ n+1 i+1,j φ n+ 1 2 i,j f φ n+ 1 2 i+1,j ψ n+1 i+1,j   = - ∆t J + 1   - σ 2 2(∆x) 2 J-1 j=0 ψ n+1 i+1,j+1 -ψ n+1 i+1,j φ n+ 1 2 i,j+1 -φ n+ 1 2 i,j + 1 σ 2 J j=0 ψ n+1 i+1,j φ n+ 1 2 i,j f φ n+ 1 2 i,j ψ n i,j   + ∆t J + 1   - σ 2 2(∆x) 2 J-1 j=0 ψ n+1 i+1,j+1 -ψ n+1 i+1,j φ n+ 1 2 i,j+1 -φ n+ 1 2 i,j + 1 σ 2 J j=0 ψ n+1 i+1,j φ n+ 1 2 i,j f φ n+ 1 2 i+1,j ψ n+1 i+1,j   = ∆t σ 2 1 J + 1 J j=0 ψ n+1 i+1,j φ n+ 1 2 i,j f ψ n+1 i+1,j φ n+ 1 2 i+1,j -f ψ n i,j φ n+ 1 2 i,j
Hence, if the assumptions are the same as in Theorem 4, we have that:

1 J + 1 J j=0 m n+1 i+1,j - 1 J + 1 J j=0 m n+1 i,j ≤ ∆t σ 2 1 J + 1 K ψ n+1 ∞ φ n+ 1 2 ∞ J j=0 ψ n+1 i+1,j φ n+ 1 2 i+1,j -ψ n i,j φ n+ 1 2 i,j ≤ ∆t σ 2 1 J + 1 K e u T σ 2 ∞ m 0 ∞ J j=0 ψ n+1 i+1,j φ n+ 1 2 i+1,j -φ n+ 1 2 i,j + φ n+ 1 2 i,j ψ n+1 i+1,j -ψ n i,j ≤ ∆t σ 2 1 J + 1 K e u T σ 2 ∞ m 0 ∞ max e u T σ 2 ∞ , m 0 ∞ J j=0 φ n+ 1 2 i+1,j -φ n+ 1 2 i,j + ψ n+1 i+1,j -ψ n i,j Now, 1 J + 1 J j=0 φ n+ 1 2 i+1,j -φ n+ 1 2 i,j ≤ 1 J + 1 J j=0 φ n+ 1 2 i+1,j - φn+ 1 2 i+1,j + φn+ 1 2 i+1,j - φn+ 1 2 i,j + φn+ 1 2 i,j -φ n+ 1 2 i,j ≤ 2 φ n+ 1 2 -φn+ 1 2 + 1 J + 1 J j=0 φn+ 1 2 i+1,j - φn+ 1 2 i,j and 1 J + 1 J j=0 ψ n+1 i+1,j -ψ n i,j ≤ ψ n+1 -ψn+1 + 1 J + 1 J j=0 ψn+1 i+1,j -ψn i,j + ψ n -ψn
As a consequence, ∀i ∈ I :

1 J + 1 J j=0 m n+1 i,j - 1 J + 1 J j=0 m n+1 0,j ≤ T σ 2 K e u T σ 2 ∞ m 0 ∞ max e u T σ 2 ∞ , m 0 ∞ ×   2 φ n+ 1 2 -φn+ 1 2 + ψ n+1 -ψn+1 + ψ n -ψn + 1 J + 1 sup l∈I J j=0 φn+ 1 2 l+1,j - φn+ 1 2 l,j + ψn+1 l+1,j -ψn l,j  
Using Theorem 4, we obtain that lim n→∞ lim sup ∆t,∆x→0

1 J + 1 J j=0 m n+1 i,j - 1 J + 1 J j=0 m n+1 0,j = 0 But, 1 J+1 J j=0 m n+1 0,j = 1 J+1 J j=0 m 0 (x j
) and this expression tends to Ω m 0 (x)dx when ∆x tends to 0 since m 0 is continuous from the hypotheses of Theorem 4.

Numerical examples

Let us consider, as an example, the following problem on Ω = (0, 1). Agents want to live near the center of the domain but prefer to live far from the others. This is modeled by the following2 function f :

f (x, ξ) = -16(x -1/2) 2 -0.1 max(0, min(5, ξ))
At the end of the period, when t = T = 0.5, we suppose that agents do not have any incentive, namely u T (x) = 0. We suppose that volatility is σ = 1. We assume that the population is initially distributed as:

m 0 (x) = µ(x) 1 0 µ(x )dx where µ(x) = 1 + 0.2 cos π 2x - 3 2 2
In this case, we solved the numerical scheme with a Newton method for ∆t = 0.01 and ∆x = 0.02 and we obtained the following results for3 (φ, ψ, u, m, ∇u). The iterations in n were stopped once the difference φ n+ 1 2 ψ n+1 -φ n-1 2 ψ n ∞ was below 10 -7 (here 5 steps).

We see that quite rapidly the population as a whole gets close to what would be the stationary equilibrium. Then, just before the time horizon, the agents have no more incentive to be around 0 and do not pay the cost associated to the stationary equilibrium. Hence the distribution spreads just before the end. Now, we can study, using the preceding choice of u T and f , the complexity of the scheme and its convergence speed.

We first consider the time spent by the algorithm for φ n+ 1 2 ψ n+1 -φ n-1 2 ψ n ∞ to be below 10 -7 for the same parameters as above except ∆t and ∆x that vary: We see that the computing time is linear in the number of points in space if the number of points in time is fixed and similarly that the computing time is linear in the number of points in time if the number of points in space is fixed. Now, if ∆t and ∆x are proportional, we obtain that the computing time is quadratic: Empirically, we therefore have that the computing time is of order O (∆x∆t) -1 .

Then, we can also consider a reference solution calculated on a grid with 301 points in both space and time and find the convergence speed with respect to this reference solution. We consider first the computation of a solution -the algorithm being stopped as soon as φ n+ 1 2 ψ n+1 -φ n-1 2 ψ n ∞ is below 10 -7 -for ∆x = 1 150 and different values of ∆t and calculate the error in discrete uniform norm: We consider then the computation of a solution for ∆t = 1 300 and different values of ∆x and calculate the error in discrete uniform norm4 : We see that the scheme seems to be of order 1 in both time and space. Now, we can look at the evolution of the total mass of m n+1 , as n grows. Because of Proposition 6 and Proposition 12 we overall expect to have a loss of mass as t grows, although this may not be true when ψ n and ψ n+1 are really close, namely for large n. This is what we observe empirically: We see that the smaller the σ, the slower the algorithm. This is due to the increasing number of steps in n to obtain φ n+ 1 2 ψ n+1 -φ n-1 2 ψ n ∞ < 10 -7 and is also due to the increasing number of steps in the Newton scheme inside each step. The rationale behind this is that a small change in σ does not change u nor m by a lot. However, φ = e u σ 2 and ψ = me -u σ 2 are largely influenced by σ: decreasing σ increases φ and decreases ψ by a lot.

Theorem 4 .

 4 Let us suppose, in addition to the hypotheses made in Theorem 2 that u T , m 0 and f are so that ∀n ∈ N,φ n+ 1 2 , ψ n ∈ C 1,2 ([0, T ] × [0, 1]) and φ, ψ ∈ C 1,2 ([0, T ] × [0, 1]) Then: lim n→∞ lim sup ∆t,∆x→0 φ n+ 1 2 -φ = 0 lim n→∞ lim sup ∆t,∆x→0ψ n+1 -ψ = 0Proof:

Figure 1 :Figure 2 :

 12 Figure 1: Solution for φ

Figure 3 :Figure 4 :

 34 Figure 3: Solution for u

Figure 5 :

 5 Figure 5: Solution for α = ∇u

Figure 6 :

 6 Figure 6: Left: Computing time for ∆t = 0.005 and different values of ∆x. Right: Computing time for ∆x = 0.01 and different values of ∆t.

  time) log(number of points in space)

Figure 7 :

 7 Figure 7: Computing time for ∆x = 2∆t for different values of ∆t.

Figure 8 :

 8 Figure 8: Error in discrete uniform norm for φ (left) and ψ (right) with ∆x = 1 150 and different values of ∆t. The reference solution has been calculated with ∆x = 2∆t = 1 300 .

Figure 9 :

 9 Figure 9: Error in discrete uniform norm for φ (left) and ψ (right) with ∆t = 1 300 and different values of ∆x. The reference solution has been calculated with ∆x = 2∆t = 1 300 .

Figure 10 :

 10 Figure 10: Evolution of the total mass of m n+1 . Top left: n = 0, Top right; n = 1, Bottom left: n = 2, Bottom right: n = 3.

Finally, to comeFigure 11 :

 11 Figure 11: Computing time for ∆x = 0.02 and ∆t = 0.005 for different values of σ.

In terms of the initial mean field game problem, the optimal control ∇u and the subsequent distribution m are not changed if we subtract f ∞ to f .

The bounds on ξ guarantee that f is bounded. In practice, these bounds are not binding.

α = ∇u is the optimal control in this context.

The fact that the convergence speeds up for small ∆x in the second set of graphs may be due to the fact that we did not consider a closed-form reference solution but, rather, an approximation of the solution calculated on a grid.