Numerical Method for FBSDEs of McKean-Vlasov Type
NNUMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE
JEAN-FRANÇOIS CHASSAGNEUX ˚ , DAN CRISAN : AND FRANÇOIS DELARUE ; Abstract.
This paper is dedicated to the presentation and the analysis of a numericalscheme for forward-backward SDEs of the McKean-Vlasov type, or equivalently forsolutions to PDEs on the Wasserstein space. Because of the mean field structure ofthe equation, earlier methods for classical forward-backward systems fail. The schemeis based on a variation of the method of continuation. The principle is to implementrecursively local Picard iterations on small time intervals.We establish a bound for the rate of convergence under the assumption that thedecoupling field of the forward-bakward SDE (or equivalently the solution of the PDE)satisfies mild regularity conditions. We also provide numerical illustrations. Introduction
In this paper, we investigate a probabilistic numerical method to approximate thesolution of the following non-local PDE B t U p t, x, µ q ` b p x, U p t, x, µ q , ν q ¨ B x U p t, x, µ q`
12 Tr rB xx U p t, x, µ q a p x, µ qs ` f ` x, U p t, x, µ q , B x U p t, x, µ q σ p x, µ q , ν ˘ ` ż R d B µ U p t, x, µ qp υ q ¨ b p υ, U p t, υ, ν q , ν q d µ p υ q ` ż R d
12 Tr rB x B µ U p t, x, µ qp υ q a p υ, µ qs d µ p υ q “ , (1)for p t, x, µ q P r , T q ˆ R d ˆ P p R d q with the terminal condition U p T, ¨q “ g p¨q , where ν is a notation for the image of the probability measure µ by the mapping R d Q x ÞÑp x, U p t, x, µ qq P R d . Above, a p x, µ q “ r σσ : sp x, µ q . The set P p R d q is the set of probabil-ity measures with a finite second-order moment, endowed with the Wasserstein distancei.e. W p µ, µ q : “ inf π ˆż R d ˆ R d | x ´ x | d π p x, x q ˙ , for p µ, µ q P P p R d q ˆ P p R d q , the infimum being taken over the probability distributions π on R d ˆ R d whose marginals on R d are respectively µ and µ .Whilst the first two lines in (1) form a classical non-linear parabolic equations, the lasttwo terms are non-standard. Not only are they non-local, in the sense that the solution orits derivatives are computed at points υ different from x , but also they involve derivativesin the argument µ , which lives in a space of probability measures. In this regard, thenotation B µ U p t, x, µ qp υ q denotes the so-called Wasserstein derivative of the function U in the direction of the measure, computed at point p t, x, µ q and taken at the continuous ∗ Laboratoire de Probabilités et Modèles aléatoires, Université Paris Diderot. [email protected] † Department of Mathematics, Imperial College London. [email protected] ‡ Laboratoire Jean-Alexandre Dieudonné, Université de Nice Sophia-Antipolis. [email protected] a r X i v : . [ m a t h . P R ] M a r JEAN-FRANÇOIS CHASSAGNEUX, DAN CRISAN AND FRANÇOIS DELARUE coordinate υ . We provide below a short reminder of the construction of this derivative,as introduced by Lions, see [12] or [17, Chap. 5].These PDEs arise in the study of large population stochastic control problems, either ofmean field game type, see for instance [12, 13, 20, 30] or [18, Chap. 12] and the referencestherein, or of mean field control type, see for instance [9, 10, 20, 33]. In both cases, U plays the role of a value function or, when the above equation is replaced by a systemof equations of the same form, the gradient of the value function. Generally speaking,these types of equations are known as “master equations”. We refer to the aforementionedpapers and monographes for a complete overview of the subject, in which existence anduniqueness of classical or viscosity solutions have been studied. In particular, in ourprevious paper [20], we tackled classical solutions by connecting U with a system of fullycoupled Forward-Backward Stochastic Differential Equations of the McKean-Vlasov type(MKV FBSDE), for which U plays the role of a decoupling field. We also refer to [18,Chap. 12] for a similar approach.In the current paper, we build on this link to design our numerical method.The connection between U and FBSDEs may be stated as follows. Basically, U may bewritten as U p t, x, µ q “ Y t,x,µt for all p t, x, µ q P r , T s ˆ R d ˆ P p R d q , where Y t,x,µ togetherwith p X t,x,µ , Z t,x,µ q solves the following standard FBSDE: X t,x,µs “ x ` ż st b p X t,x,µr , Y t,x,µr , r X t,ξr , Y t,ξr sq d r ` ż st σ p X t,x,µr , r X t,x,µr sq d W r (2) Y t,x,µs “ g p X t,x,µT , r X t,ξT sq ` ż Ts f p X t,x,µr , Y t,x,µr , Z t,x,µr , r X t,ξr , Y t,ξr sq d r ´ ż Ts Z t,x,µr ¨ d W r , (3)which is parametrized by the law of the following MKV FBSDE: X t,ξs “ ξ ` ż st b p X t,ξr , Y t,ξr , r X t,ξr , Y t,ξr sq d r ` ż st σ p X t,ξr , r X t,ξr sq d W r (4) Y t,ξs “ g p X t,ξT , r X t,ξT sq ` ż Ts f p X t,ξr , Y t,ξr , Z t,ξr , r X t,ξr , Y t,ξr sq d r ´ ż Ts Z t,ξr ¨ d W r , (5)where p W t q ď t ď T is a Brownian motion and ξ has µ as distribution. In the previousequations and in the sequel, we use the notation r θ s for the law of a random variable θ . In particular, in the above, we have that r ξ s “ µ . So, to obtain an approximationof U p t, x, µ q given by the initial value of (3), our strategy is to approximate the system(4)-(5) as its solution appears in the coefficients of (2)-(3). In this regard, our approachis probabilistic.Actually, our paper is not the first one to address the numerical approximation ofequations of the type (1) by means of a probabilistic approach. In its PhD disserta-tion, Alanko [4] develops a numerical method for mean field games based upon a Picarditeration: Given the proxy for the equilibrium distribution of the population (which isrepresented by the mean field component in the above FBSDE), one solves for the valuefunction by approximating the solution of the (standard) BSDE associated with the con-trol problem; given the solution of the BSDE, we then get a new proxy for the equilibriumdistribution and so on... Up to a Girsanov transformation, the BSDE associated with thecontrol problem coincides with the backward equation in the above FBSDEs. In [4], theGirsanov transformation is indeed used to decouple the forward and backward equationsand it is the keystone of the paper to address the numerical impact of the change ofmeasure onto the mean field component. Under our setting, this method would more or UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 3 less consist in solving for the backward equation given a proxy for the forward equationand then in iterating, which is what we call the
Picard method for the FBSDE system.Unfortunately, convergence of the Picard iterations is a difficult issue, as the convergenceis known in small time only, see the numerical examples in Section 4 below. It is indeedwell-known that Picard theorem only applies in small time for fully coupled problems.In this regard, it must be stressed that our system (4)-(5) is somehow doubly coupled,once in the variable x and once in the variable µ , which explains why a change measuredoes not permit to decouple it entirely.The goal of our paper is precisely to go further and to propose an algorithm whoseconvergence is known on any interval of a given length (observe that the convergenceis not studied in [4]). In the classical case, this question has been addressed by severalauthors, among which [21, 22] and [7], but all these methods rely on the Markov structureof the problem. Here, the Markov property is true but at the price of regarding the entire R d ˆ P p R d q as state space: The fact that the second component is infinite dimensionalmakes intractable the complexity of these approaches. To avoid any similar problem, weuse a pathwise approach for the forward component; it consists in iterating successivelythe Picard method on small intervals, all the Picard iterations being implemented with atree approximation of the Brownian motion. This strategy is inspired from the methodof continuation, the parameter in the continuation argument being the time length T itself. The advantage for working on a tree is twofold: as we said, we completely bypassany Markov argument; also, we get, not only, an approximation of the system (4)-(5)but also, for free, an approximation of the system (2)-(3), which “lives” on a subtreeobtained by conditioning on the initial root. We prove that the method is convergentand provide a rate of convergence for it. Numerical examples are given in Section 4. Ofcourse, the complexity remains pretty high in comparison with the methods developedin the classical non McKean-Vlasov case. This should not come as a surprise since, aswe already emphasized, the problem is somehow infinite dimensional.We refer the interested reader to the following papers for various numerical methods,based upon finite differences or variational approaches, for mean field games: [1, 2, 3]and [6, 26, 25]. Recently, a Markov chain approximation method was also suggested in[5].The paper is organized as follows. The method for the system (4)-(5) is exposed inSection 2. The convergence is addressed in Section 3. In Section 4, we explain how tocompute in practice U p t, x, µ q (and thus approximate (2)-(3)) from the approximation ofthe sole (4)-(5) and we present some numerical results validating empirically the conver-gence results obtained in Section 3. We collect in the appendix some key results for theconvergence analysis.2. A new algorithm for coupled forward backward systems
As announced right above, we will focus on the approximation of the following typeof McKean-Vlasov forward-backward stochastic differential equation: dX t “ b ` X t , Y t , r X t , Y t s ˘ dt ` σ ` X t , r X t s ˘ dW t ,dY t “ ´ f ` X t , Y t , Z t , r X t , Y t s ˘ d t ` Z t ¨ dW t , t P r , T s ,Y T “ g ` X T , r X T s ˘ and X “ ξ , (6)for some time horizon T ą . Throughout the analysis, the equation is regarded on a com-plete filtered probability space p Ω , F , F , P q , equipped with a d -dimensional F -Brownian JEAN-FRANÇOIS CHASSAGNEUX, DAN CRISAN AND FRANÇOIS DELARUE motion p W t q ď t ď T . To simplify, we assume that the state process p X t q ď t ď T is of thesame dimension. The process p Y t q ď t ď T is -dimensional. As a result, p Z t q ď t ď T is d -dimensional.In (6), the three processes p X t q ď t ď T , p Y t q ď t ď T and p Z t q ď t ď T are required to be F -progressively measurable. Both p X t q ď t ď T and p Y t q ď t ď T have continuous trajectories.Generally speaking, the initial condition X is assumed to be square-integrable, butat some point, we will assume that X belongs to L p p Ω , F , P ; R d q , for some p ą .Accordingly, p X t q ď t ď T , p Y t q ď t ď T and p Z t q ď t ď T must satisfy: ~p X, Y, Z q~ r ,T s : “ E „ sup ď t ď T ` | X t | ` | Y t | ˘ ` ż T | Z t | dt { ă 8 . The domains and codomains of the coefficients are defined accordingly. The assump-tion that σ is assumed to be independent of the variable y is consistent with the globalsolvability results that exist in the literature for equations like (6). For instance, it coverscases coming from optimization theory for large mean field interacting particle systems.We refer to our previous paper [20] for a complete overview on the subject, togetherwith the references [8, 12, 17, 18, 19]. In light of the examples tackled in [20], the factthat b is independent of z may actually seem more restrictive, as it excludes cases whenthe forward-backward system of the McKean-Vlasov type is used to represent the valuefunction of the underlying optimization problem. It is indeed a well-known fact that,with or without McKean-Vlasov interaction, the value function of a standard optimiza-tion problem may be represented as the backward component of a standard FBSDEwith a drift term depending upon the z variable. This says that, in order to tackle theaforementioned optimization problems of the mean field type by means of the numericalmethod investigated in this paper, one must apply the algorithm exposed below to thePontryagin system. The latter one is indeed of the form (6), provided that Y is allowedto be multi-dimensional. (Below, we just focus on the one-dimensional case, but theadaptation is straightforward.)In fact, our choice for assuming b to be independent of z should not come as a surprise.The same assumption appears in the papers [21, 22] dedicated to the numerical analysisof standard FBSDEs, which will serve us as a benchmark throughout the text. Seehowever Remark 4.Finally, the fact that the coefficients are time-homogeneous is for convenience only.As a key ingredient in our analysis, we use the following representation result given ine.g. Proposition 2.2 in [20], Y ξt : “ U p t, X ξt , r X ξt sq , (7)where U : r , T s ˆ R d ˆ P p R d q Ñ R is assumed to be the classical solution, in the senseof [20, Definition 2.6], to (1). In this regard, the derivative with respect to the measureargument is defined according to Lions’ approach to the Wasserstein derivative. In short,the lifting ˆ U of U to L p Ω , F , P ; R d q , which we define by ˆ U p t, x, ξ q “ U p t, x, r ξ sq , t P r , T s , x P R d , ξ P L p Ω , F , P ; R d q , is assumed to be Fréchet differentiable. Of course, this makes sense as long as the space p Ω , F , P q is rich enough so that, for any µ P P p R d q , there exists a random variable ξ P L p Ω , F , P ; R d q such that ξ „ µ . So, in the sequel, p Ω , F , P q is assumed to be atom-less, which makes it rich enough. A crucial point with Lions’ approach to Wassersteindifferential calculus is that the Fréchet derivative of ˆ U in the third variable, which can be UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 5 identified with a square-integrable random variable, may be represented at point p t, x, ξ q as B µ U p t, x, r ξ sqp ξ q for a mapping B µ U p t, x, µ qp¨q : R d Q v ÞÑ B µ U p t, x, µ qp v q P R d . Thislatter function plays the role of Wasserstein derivative of U in the measure argument. Todefine a classical solution, it is then required that R d Q v ÞÑ B µ U p t, x, µ qp v q is differen-tiable, both B µ U and B v B µ U being required to be continuous at any point p t, x, µ, v q suchthat v is in the support of µ . Assumptions.
Our analysis requires some minimal regularity assumptions on the co-efficients b , σ , f and the function U . As for the coefficients functions, we assume thatthere exists a constant Λ ě such that:- p H q : The functions b , σ , f and g are Λ -Lipschitz continuous in all the variables, thespace P p R d q being equipped with the Wasserstein distance W . Moreover, the function σ is bounded by Λ .We now state the main assumptions on U , see Remark 1 for comments.- p H q : for any t P r , T s and ξ P L p Ω , F t , P ; R d q , the McKean-Vlasov forward-backwardsystem (6) set on r t, T s instead of r , T s with X t “ ξ as initial condition at time t has aunique solution p X t,ξs , Y t,ξs , Z t,ξs q t ď s ď T ; in parallel, U is the classical solution, in the senseof [20, Definition 2.6], to (1); and U and its derivatives satisfy | U p t, x, µ q ´ U p t, x, µ q| ` |B x U p t, x, µ q ´ B x U p t, x, µ q| ď Λ W p µ, µ q , (8) |B x U p t, x, µ q| ` }B µ U p t, x, r ξ sqp ξ q} ď Λ , (9) |B xx U p t, x, µ q| ` }B υ B µ U p t, x, r ξ sqp ξ q} ď Λ , (10)and |B xx U p t, x, µ q ´ B xx U p t, x , µ q| ď Λ | x ´ x | , (11)for p t, x, x , ξ q P r , T s ˆ R d ˆ R d ˆ L p Ω , F , P ; R d q and µ, µ P P p R d q . Also, we requirethat | U p t ` h, x, r ξ sq ´ U p t, x, r ξ sq| ` |B x U p t ` h, x, r ξ sq ´ B x U p t, x, r ξ sq| ď Λ h ` ` | x | ` } ξ } ˘ , (12)and for all h P r , T q , p t, x q P r , T ´ h s ˆ R d , ξ P L p Ω , F , P ; R d q and v, v P R d , |B υ B µ U p t, x, r ξ sqp υ q ´ B υ B µ U p t, x, r ξ sqp υ q| ď Λ t ` | υ | α ` | υ | α ` } ξ } α u | υ ´ υ | , (13)for some α ą . Remark 1. In [20] , it is shown that, under some conditions on the coefficients b , f and σ , the PDE (1) has indeed a unique classical solution which satisfies the assumption p H q .(1) Estimate (13) is obtained by combining Definition 2.6 and Proposition 3.9 in [20] .A major difficulty in the analysis provided below is the fact that α may be largerthan 1, in which case the Lipschitz bound for the second order derivative is super-linear. This problem is proper to the McKean-Vlasov structure of the equationand does not manifest in the classical setting, compare for instance with [21, 22] .Below, we tackle two cases: the case when α ď , which has been investigatedin [13] and [18, Chap. 12] under stronger conditions on the coefficients, and thecase when α ą but U is bounded.(2) Estimates (8) - (12) are required to control the convergence error when the coeffi-cients ( b or f ) depend on Z . JEAN-FRANÇOIS CHASSAGNEUX, DAN CRISAN AND FRANÇOIS DELARUE (a) The estimate (8) can be retrieved from the computations made in [20] . Seethe comments at the bottom of page 60, near equation p . q .(b) The estimate (12) comes from the theory of FBSDEs (without McKean-Vlasov interaction). Indeed, using the Lipschitz property of U and B x U inthe variable µ , it suffices to prove | U p t ` h, x, r X t,ξt ` h sq ´ U p t, x, r X t,ξt sq| ` |B x U p t ` h, x, r X t,ξt ` h sq ´ B x U p t, x, r X t,ξt sq|ď Λ h ` ` | x | ` } ξ } ˘ . As stated in Proposition 2.2 in [20] , for ξ „ µ , U p s, x, r X t,ξs sq “ u t,µ p s, x q where u t,µ is solution to a quasi-linear PDE. Then the estimate (12) followsfrom standard results on non-linear PDEs, see e.g. Theorem 2.1 in [21] .In comparison with the assumption used in [21] , the condition p H q is more demanding.In [21] , there is no need for assuming the second-order derivative to be Lipschitz in space.This follows from the fact that, here, we approximate the Brownian increments by randomvariables taking a small number of values, whilst in [21] , the Brownian increments areapproximated by a quantization grid with a larger number of points. In this regard, ourapproach is closer to the strategy implemented in [22] . Description.
The goal of the numerical method exposed in the paper is to approx-imate U . The starting point is the formula (6) and, quite naturally, the strategy is toapproximate the process p X ξ , Y ξ , Z ξ q : “ p X ,ξ , Y ,ξ , Z ,ξ q .Generally speaking, this approach raises a major difficulty, as it requires to handle thestrongly coupled forward-backward structure of (6). Indeed, theoretical solutions to (6)may be constructed by means of basic Picard iterations but in small time only, whichcomes in contrast with similar results for decoupled forward or backward equations forwhich Picard iterations converge on any finite time horizon. In the papers [21, 22] –whichdeal with the non McKean-Vlasov case–, this difficulty is bypassed by approximating thedecoupling field U at the nodes of a time-space grid. Obviously, this strategy is hopelessin the McKean-Vlasov setting as the state variable is infinite dimensional; discretizing iton a grid would be of a non-tractable complexity. This observation is the main rationalefor the approach exposed below.Our method is a variation of the so-called method of continuation . In full generality,it consists in increasing step by step the coupling parameter between the forward andbackward equations. Of course, the intuition is that, for a given time length T , thePicard scheme should converge for very small values of the coupling parameter. The goalis then to insert the approximation computed for a small coupling parameter into thescheme used to compute a numerical solution for a higher value of the coupling parameter.Below, we adapt this idea, but we directly regard T itself as a coupling parameter. So weincrease T step and by step and, on each step, we make use of a Picard iteration basedon the approximations obtained at the previous steps.This naturally motivates the introduction of an equidistant grid (cid:60) “ t r “ , . . . , r N “ T u of the time interval r , T s , with r k “ kδ and δ “ TN for N ě . In the following weshall consider that δ is “small enough” and state more precisely what it means in themain results, see Theorem 5 and Theorem 7. UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 7
For ď k ď N ´ , we consider intervals I k “ r r k , T s and on each interval, thefollowing FBSDE, for ξ P L p F r k q (which is a shorter notation for L p Ω , F r k , P ; R d q ): X t “ ξ ` ż tr k b ` X s , Y s , r X s , Y s s ˘ d s ` ż tr k σ p X s , r X s sq d W s , (14) Y t “ g ` X T , r X T s ˘ ` ż Tt f ` X s , Y s , Z s , r X s , Y s s ˘ d s ´ ż Tt Z s ¨ d W s . (15)Picard iterations. We need to compute backwards the value of U p r k , ξ, r ξ sq for some ξ P L p F r k q , ď k ď N ´ . We are then going to solve the FBSDE (14)-(15) on theinterval I k . As explained above, the difficulty is the arbitrariness of T : When k is large, I k is of a small length, but this becomes false as k decreases. Fortunately, we can rewritethe forward-backward system on a smaller interval at the price of changing the termi-nal boundary condition. Indeed, from p H q , we know that p X r k ,ξs , Y r k ,ξs , Z r k ,ξs q r k ď s ď r k ` solves: " X t “ ξ ` ş tr k b ` X s , Y s , r X s , Y s s ˘ d s ` ş tr k σ ` X s , r X s s ˘ d W s ,Y t “ U ` r k ` , X r k ` , r X r k ` s ˘ ` ş r k ` t f ` X s , Y s , Z s , r X s , Y s s ˘ d s ´ ş r k ` t Z s ¨ d W s , for t P r r k , r k ` s .If δ is small enough, a natural approach is to introduce a Picard iteration schemeto approximate the solution of the above equation. To do so, one can implement thefollowing recursion (with respect to the index j ): X jt “ ξ ` ş tr k b ` X js , Y js , r X js , Y js s ˘ d s ` ş tr k σ ` X js , r X js s ˘ d W s ,Y jt “ U ` r k ` , X j ´ r k ` , r X j ´ r k ` s ˘ ` ş r k ` t f ` X j ´ s , Y js , Z js , r X j ´ s , Y js s ˘ d s ´ ş r k ` t Z js ¨ d W s . with p X s “ ξ ` ş tr k b ` X s , , r X s , s ˘ d s ` ş tr k σ ` X s , r X s s ˘ d W s q r k ď s ď r k ` and p Y s “ q r k ď s ď r k ` .It is known that, for δ small enough, p X j , Y j , Z j q Ñ j Ñ8 p X, Y, Z q , in the sense that ~p X j ´ X, Y j ´ Y, Z j ´ Z q~ r r k ,r k ` s Ñ j Ñ8 .But in practice we will encounter three main difficulties.(1) The procedure has to be stopped after a given number of iterations J .(2) The above Picard iteration assumes the perfect knowledge of the map U at time r k , but U is exactly what we want to compute...(3) The solution has to be discretized in time and space. Ideal recursion.
We first discuss 1) and 2) above. The main idea is to use a recursivealgorithm (with a new recursion, but on the time parameter).Namely, for k ď N ´ , we assume that we are given a solver which computes solver[ k ` ]( ξ ) “ U p r k ` , ξ, r ξ sq ` (cid:15) k ` p ξ q , (16)where (cid:15) is an error made, for any ξ P L p F r k ` q . We shall sometimes refer to solver[ k ` ]( ¨ ) as “the solver at level k ` ”.Taking these observations into account, we first define an ideal solver, which assumesthat each Picard iteration in the approximation of the solution of the forward-backwardsystem can be perfectly computed. We denote it by picard[]() . Accordingly, we identify(for the time being) solver[ k ` ]() with picard[ k ` ]() . Given picard[ k ` ]() , JEAN-FRANÇOIS CHASSAGNEUX, DAN CRISAN AND FRANÇOIS DELARUE picard[ k ]() is defined as follows. $’&’% ˜ X k,jt “ ξ ` ş tr k b ` ˜ X k,js , ˜ Y k,js , r ˜ X k,js , ˜ Y k,js s ˘ d s ` ş tr k σ ` ˜ X k,js , r ˜ X k,js s ˘ d W s , ˜ Y k,jt “ picard[ k ` ]( ˜ X k,j ´ r k ` ) ´ ş r k ` t ˜ Z k,js ¨ d W s ` ş r k ` t f ` ˜ X k,j ´ s , ˜ Y k,js , ˜ Z k,js , r ˜ X k,j ´ s , ˜ Y k,js s ˘ d s , (17)for j ě and with ´ ˜ X k, s “ ξ ` ş tr k b p X k, s , , r X k, s , sq d s ` ş tr k σ p X k, s , r X k, s sq d W s ¯ r k ď s ď r k ` ,and p ˜ Y k, s “ q r k ď s ď r k ` . We then define picard[ k ]( ξ ) : “ Y k,Jr k and (cid:15) k p ξ q : “ Y k,Jr k ´ U p r k , ξ, r ξ sq , where J ě is the number of Picard iterations.At level N ´ , which is the last level for our recursive algorithm, the Picard iterationscheme is given by $’’’&’’’% ˜ X N ´ ,jt “ ξ ` ş tr N ´ b ` ˜ X N ´ ,js , ˜ Y N ´ ,js , r ˜ X N ´ ,js , ˜ Y N ´ ,js s ˘ d s ` ş tr N ´ σ ` ˜ X N ´ ,js , r ˜ X N ´ ,js s ˘ d W s , ˜ Y N ´ ,jt “ g p ˜ X N ´ ,j ´ T , r ˜ X N ´ ,j ´ T sq ´ ş Tt ˜ Z N ´ ,js ¨ d W s ` ş Tt f ` ˜ X N ´ ,j ´ s , ˜ Y N ´ ,js , ˜ Z N ´ ,js , r ˜ X N ´ ,j ´ s , ˜ Y N ´ ,js s ˘ d s . (18)Here, the terminal condition g is known and the error comes from the fact that the Picarditeration is stopped. It is then natural to set, for ξ P L p F T q , picard[ N ]( ξ ) “ g p ξ, r ξ sq and (cid:15) N p ξ q “ . (19)Practical implemention. As already noticed in 3) above, it is not possible to solve thebackward and forward equations in (17) perfectly, even though the system is decoupled.Hence, we need to introduce an approximation that can be implemented in practice.Given a continuous adapted input process X “ p X s q r k ď s ď r k ` such that E r sup r k ď s ď r k ` | X s | s ă8 and η P L p Ω , F r k ` , P ; R q , we thus would like to solve ˜ X t “ X r k ` ş tr k b ` ˜ X s , ˜ Y s , r ˜ X s , ˜ Y s s ˘ d s ` ş tr k σ ` ˜ X s , r ˜ X s s ˘ d W s ˜ Y t “ η ` ş r k ` t f ` X s , ˜ Y s , ˜ Z s , r X s , ˜ Y s s ˘ d s ´ ş r k ` t ˜ Z s ¨ d W s , for t P r r k , r k ` s .Let π be a discrete time grid of r , T s such that (cid:60) Ă π , π : “ t t : “ ă ¨ ¨ ¨ ă t n : “ T u and | π | : “ max i ă n p t i ` ´ t i q . (20)For ď k ď N ´ , we note π k : “ t t P π | r k ď t ď r k ` u and for later use, we define theindices p j k q ď k ď N as follows π k “ t t j k : “ r k ă ¨ ¨ ¨ ă t i ă ¨ ¨ ¨ ă r k ` “ : t j k ` u , for all k ă N . So, instead of a perfect solver for an iteration of the Picard scheme (17),we assume that we are given a numerical solver, denoted by solver[ k ]( ¯ X , η , f ) , whichcomputes an approximation of the process p ˜ X s , ˜ Y s , ˜ Z s q r k ď s ď r k ` on π k for a discretiza-tion p ¯ X t q t P π k of the time continuous process p X s q r k ď s ď r k ` . The output is denoted by p ¯ X t , ¯ Y t , ¯ Z t q t P π k . In parallel, we call input the triplet formed by the random variable η ,the discrete-time process p ¯ X t q t P π k and the driver f of the backward equation. In short,the output is what the numerical solver returns after one iteration in the Picard schemewhen the discrete input is p η, ¯ X , f q . Pay attention that, in contrast with b and σ , we UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 9 shall allow f to vary; this is the rationale for regarding it as an input. However, whenthe value of f is clear, we shall just regard the input as the pair p η, p ¯ X t q t P π k q .The full convergence analysis, including the discretization error, will be discussedin the next section in the following two cases: first for a generic (or abstract) solver solver[](,,) and second for an explicit solver, as given in the example below. Example 2.
This example is the prototype of the solver solver[](,,) . We consideran approximation of the Brownian motion obtained by quantization of the Brownian in-crements. At every time t P π , we denote by ¯ W t the value at time t of the discretizedBrownian motion. It may expressed as ¯ W t i : “ i ´ ÿ j “ ∆ ¯ W j , where ∆ ¯ W j : “ h j (cid:36) j , (cid:36) j : “ Γ d ´ h ´ j ` W t j ` ´ W t j ˘¯ , Γ d mapping R d onto a finite grid of R d . Importantly, Γ d is assumed to be bounded by Λ and each (cid:36) j is assumed to be centered and to have the identity matrix as covariancematrix. Of course, this is true if Γ d is of the form Γ d ` w , ¨ ¨ ¨ , w d ˘ : “ ` Γ p w q , ¨ ¨ ¨ , Γ p w d q ˘ , p w , ¨ ¨ ¨ , w d q P R d , where Γ is a bounded odd function from R onto a finite subset of R with a normalizedsecond order moment under the standard Gaussian measure. In practice, Γ d is intendedto take a small number of values. Of course, the typical example is the so-called binomialapproximation, in which case Γ is the sign function.On each interval r r k , r k ` s , given a discrete-time input process ¯ X and a terminal condi-tion η , we thus implement the following scheme (below, E t i is the conditional expectationgiven F t i ):(1) For the backward component:(a) Set as terminal condition, p ¯ Y t jk ` , ¯ Z t jk ` q “ p η, q .(b) For j k ď i ă j k ` , compute recursively ¯ Y t i “ E t i “ ¯ Y t i ` ` p t i ` ´ t i q f ` ¯ X t i , ¯ Y t i , ¯ Z t i , r ¯ X t i , ¯ Y t i s ˘‰ , ¯ Z t i “ E t i „ ∆ ¯ W i t i ` ´ t i ¯ Y t i ` . (2) For the forward component:(a) Set as initial condition, ¯ X t jk “ ¯ X r k .(b) For j k ă i ď j k ` , compute recursively ¯ X t i ` “ ¯ X t i ` b ` ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i s ˘ p t i ` ´ t i q ` σ ` ¯ X t i , r ¯ X t i s ˘ ∆ ¯ W i . Full algorithm for solver[]() . Using solver[](,,) , for each level, we can now give acompletely implementable algorithm for solver[]() . Its description is as follows.The value solver[ k ]( ξ ) , i.e. the value of the solver at level k with initial condition ξ P L p F r k q , is obtained through:(1) Initialize the backward component at ¯ Y k, t “ for t P π k and regard p ¯ X k, t q t P π k as the forward component of solver[ k ]( ξ , , ) (2) for ď j ď J (a) compute ¯ Y k,jr k ` “ solver[ k ` ]( ¯ X k,j ´ r k ` ) . (b) compute p ¯ X k,j , ¯ Y k,j , ¯ Z k,j q “ solver[ k ]( ¯ X k,j ´ , ¯ Y k,jr k ` , f ) (3) return ¯ Y k,Jr k ` .Following (19), we let solver[ N ]( ξ ) “ g p ξ, r ξ sq . (21)We first explain the initialization step. The basic idea is to set the backward componentto and then to solve the forward component as an approximation of the autonomousMcKean-Vlasov diffusion process in which the backward entry is null. Of course, thismay be solved by means of any standard method, but to make the notation shorten, wefelt better to regard the underlying solver as a specific case of a forward-backward solverwith null coefficients in the backward equation. We specify in the analysis below theconditions that this initial solver solver[](, , ) must satisfy.It is also worth noting that each Picard iteration used to define the solver at level k calls the solver at level k ` . This is a typical feature of the way the continuation methodmanifests from the algorithmic point of view. In particular, the total complexity is oforder O p J N K q , where K is the complexity of the solver solver[](,,) . In this regard,it must be stressed that, for a given length T , N is fixed, regardless of the time step | π | .Also, J is intended to be rather small as the Picard iterations are expected to convergegeometrically fast, see the numerical examples in Section 4 in which we choose J “ .However, it must be noticed that the complexity increases exponentially fast when T tends to , which is obviously the main drawback of this method. Again, we refer toSection 4 for numerical illustrations.Useful notations. Throughout the paper, } ¨ } p denotes the L p norm on p Ω , F , P q . Also, p ˆΩ , ˆ F , ˆ P q stands for a copy of p Ω , F , P q . It is especially useful to represent the Lions’derivative of a function of a probability measure and to distinguish the (somewhat artifi-cial) space used for representing these derivatives from the (physical) space carrying theWiener process. For a random variable X defined on p Ω , F , P q , we shall denote by x X y its copy on p ˆΩ , ˆ F , ˆ P q .We shall use the notations C Λ , c Λ for constants only depending on Λ (and possibly onthe dimension as well). They are allowed to increase from line to line. We shall use thenotation C for constants not depending upon the discretization parameters. Again, theyare allowed to increase from line to line. In most of the proofs, we shall just write C for C Λ , even if we use the more precise notation C Λ in the corresponding statement.2.2. A first analysis with no discretization error.
To conclude this section, we wantto understand how the error propagates through the solvers used at different levels in theideal case where the Picard iteration in (17) can be perfectly computed or equivalentlywhen the solver is given by solver[ k ]() “ picard[ k ]() . For j ě , we then denote by p ˜ X k,j , ˜ Y k,j , ˜ Z k,j q , the solution on r r k , r k ` s of (17).The main result of the section, see Theorem 5, is an upper bound for the error whenwe use picard[ ¨ ]( ¨ ) to approximate U . The proof of this theorem requires the followingproposition, which gives a local error estimate for each level. Proposition 3.
Let us define, for j P t , ¨ ¨ ¨ , J u , k P t , ¨ ¨ ¨ , N ´ u , ∆ jk : “ ››› sup t Pr r k ,r k ` s ` ˜ Y k,jt ´ U p t, ˜ X k,jt , r ˜ X k,jt sq ˘››› UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 11 then, there exist constants C Λ , c Λ such that, for ¯ δ : “ C Λ δ ă c Λ , ∆ jk ď ¯ δ j ∆ k ` j ÿ (cid:96) “ ¯ δ (cid:96) ´ e ¯ δ ›› (cid:15) k ` p ˜ X k,j ´ (cid:96)r k ` q ›› . (22)We recall that (cid:15) k p ξ q stands for the error term: (cid:15) k p ξ q “ picard[ k ]( ξ ) ´ U p r k , ξ, r ξ sq , with (cid:15) N p ξ q “ . Remark 4.
A careful inspection of the proof shows that, whenever σ depends on Y or b depends on Z , the same result holds true but with a constant C Λ depending on N . As N is fixed in practice, this might still suffice to complete the analysis of the discretizationscheme in that more general setting. Proof.
We suppose that the full algorithm is initialized at some level k P t , ¨ ¨ ¨ , N ´ u ,with an initial condition ξ P L p F r k q . As the value of the index k is fixed throughout theproof, we will drop it in the notations p ˜ X k,j , ˜ Y k,j , ˜ Z k,j q and ∆ jk .Applying Ito’s formula for functions of a measure argument, see [11, 20], we have d U p t, ˜ X jt , r ˜ X jt sq “ ˆ b p ˜ X jt , ˜ Y jt , r ˜ X jt , ˜ Y jt sq ¨ B x U p t, ˜ X jt , r ˜ X jt sq`
12 Tr “ a ` ˜ X jt , r ˜ X jt s ˘ B xx U p t, ˜ X jt , r ˜ X jt sq ‰ ` ˆ E ” b px ˜ X jt y , x ˜ Y jt y , r ˜ X jt , ˜ Y jt sq ¨ B µ U p t, ˜ X jt , r ˜ X jt sq ı ` ˆ E ”
12 Tr “ a ` x ˜ X jt y , r ˜ X jt s ˘ B υ B µ U p t, ˜ X jt , r ˜ X jt sq ‰ı ` B t U p t, ˜ X jt , r ˜ X jt sq ˙ d t ` B x U p t, ˜ X jt , r ˜ X jt sq ¨ ` σ ` ˜ X jt , r ˜ X jt s ˘ d W t ˘ . Expressing the integral in (1) as expectations on p ˆΩ , ˆ F , ˆ P q and combining with (1) and(17), we obtain d r ˇ Y jt ´ ˜ Y jt s “ ´(cid:32) b ` ˜ X jt , ˜ Y jt , r ˜ X jt , ˜ Y jt s ˘ ´ b ` ˜ X jt , ˇ Y jt , r ˜ X jt , ˇ Y jt s ˘( ¨ B x U p t, ˜ X jt , r ˜ X jt sq` p E ”(cid:32) b ` x ˜ X jt y , x ˜ Y jt y , r ˜ X jt , ˜ Y jt s ˘ ´ b ` x ˜ X jt y , x ˇ Y jt y , r ˜ X jt , ˇ Y jt s ˘( ¨ B µ U p t, ˜ X jt , r ˜ X jt sq ı ` f ` ˜ X jt , ˜ Y jt , ˜ Z jt , r ˜ X jt , ˜ Y jt s ˘ ´ f ` ˜ X jt , ˇ Y jt , ˇ Z jt , r ˜ X jt , ˇ Y jt s ˘¯ d t ` r ˇ Z jt ´ ˜ Z jt s ¨ d W t , where ˇ Y jt : “ U p t, ˜ X jt , r ˜ X jt sq and ˇ Z jt : “ B x u p t, ˜ X jt , r ˜ X jt sq σ p ˜ X jt , r ˜ X jt sq . Observe that thisargument is reminiscent of the four-step scheme, see [32].Using standard arguments from BSDE theory and p H q – p H q , we then compute ∆ j ď e Cδ ›› U p r k ` , ˜ X jr k ` , r ˜ X jr k ` sq ´ ˜ Y jr k ` ›› ď e Cδ ´›› (cid:15) k ` p ˜ X j ´ r k ` q ›› ` ›› U p r k ` , ˜ X jr k ` , r ˜ X jr k ` sq ´ U p r k ` , ˜ X j ´ r k ` , r ˜ X j ´ r k ` sq ›› ¯ , recalling ˜ Y jr k ` “ picard[ k ` ]( ˜ X j ´ r k ` ) and (16). Since U is Lipschitz, we have ∆ j ď e Cδ ´›› (cid:15) k ` p ˜ X j ´ r k ` q ›› ` L ›› ˜ X jr k ` ´ ˜ X j ´ r k ` ›› ¯ . (23) We also have that ˜ X jt ´ ˜ X j ´ t “ ż tr k (cid:32) b ` ˜ X js , ˜ Y js , r ˜ X jt , ˜ Y jt s ˘ ´ b ` ˜ X j ´ s , ˜ Y j ´ s , r ˜ X j ´ t , ˜ Y j ´ t s ˘( d s ` ż tr k (cid:32) σ ` ˜ X js , r ˜ X js s ˘ ´ σ ` ˜ X j ´ s , r ˜ X j ´ s s ˘( d W s . Using usual arguments (squaring, taking the sup, using Bürkholder-Davis-Gundy in-equality), we get, since b and σ are Lipschitz continuous, ››› sup t Pr r k ,r k ` s | ˜ X jt ´ ˜ X j ´ t | ››› ď C ´ δ ››› sup t Pr r k ,r k ` s | ˜ Y jt ´ ˜ Y j ´ t | ››› ` δ ››› sup t Pr r k ,r k ` s | ˜ X jt ´ ˜ X j ´ t | ››› ¯ . Observing that | ˜ Y js ´ ˜ Y j ´ s | ď | ˜ Y js ´ U p s, ˜ X js , r ˜ X js sq| ` | ˜ Y j ´ s ´ U p s, ˜ X j ´ s , r ˜ X j ´ s sq|` Λ p| ˜ X j ´ s ´ ˜ X js | ` } ˜ X j ´ s ´ ˜ X js } q , we obtain, for δ small enough, ››› sup t Pr r k ,r k ` s | ˜ X jt ´ ˜ X j ´ t | ››› ď Cδ p ∆ j ` ∆ j ´ q . (24)Combining the previous inequality with (23), we obtain, for δ small enough, ∆ j ď e Cδ ›› (cid:15) k ` p ˜ X j ´ r k ` q ›› ` Cδ ∆ j ´ , which by induction leads to ∆ j ď p Cδ q j ∆ ` j ÿ (cid:96) “ p Cδ q (cid:96) ´ e Cδ ›› (cid:15) k ` p ˜ X j ´ (cid:96)r k ` q ›› , and concludes the proof. l We now state the main result of this section, which explains how the local error inducedby the fact that the Picard iteration is stopped at rank J propagates through the variouslevels k “ N ´ , ¨ ¨ ¨ , . Theorem 5.
We can find two constants C Λ , c Λ ą and a continuous non-decreasingfunction B : R ` Ñ R ` matching in , only depending on Λ , such that, for ¯ δ : “ C Λ δ ă min p c Λ , q and β ě B p ¯ δ q satisfying p J ´ q Λ¯ δ J e βC Λ T e β ¯ δ ´ ď (25) where J is the number of Picard iterations in a period, it holds, for any period k Pt , ¨ ¨ ¨ , N u and ξ P L p F r k q , } solver[ k ]( ξ ) ´ U p r k , ξ, r ξ sq} ď Λ e βC Λ T β ¯ δ J ´ ` ` ›› P ‹ r k ,T p ξ q ›› ˘ , (26) where P r k ,t p ξ q is the solution at time t of the stochastic differential equation dX s “ b ` X s , , r X s , s ˘ d s ` σ ` X s , r X s s ˘ d W s , with X r k “ ξ as initial condition, and P ‹ r k ,t p ξ q “ sup s Pr r k ,t s | P r k ,s p ξ q| . UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 13
Of course, it is absolutely straightforward to bound ››› P ‹ r k ,T p ξ q ››› by C p ` } ξ } q in (26).Theorem 5 may be restated accordingly, but the form used in the statement is morefaithful to the spirit of the proof. Proof.
We prove the claim by an induction argument. We show below that for all k P t , . . . , N u , ››› (cid:15) k p ξ q ››› “ } solver[ k ]( ξ ) ´ U p r k , ξ, r ξ sq} ď θ k ´ ` ›› P ‹ r k ,T p ξ q ›› ¯ , (27)where p θ k q k “ , ¨¨¨ ,N ´ is defined by the following backward induction: θ N : “ , recall (19),and for k P t , ¨ ¨ ¨ , N ´ u , θ k : “ Λ¯ δ J ` e β ¯ δ θ k ` , (28)where β is such that ˆ γ ` γ ¯ δe γ ¯ δ p γ ` Λ1 ´ ¯ δ q ˙ ď e β ¯ δ , with γ : “ e ¯ δ ´ ¯ δ . (29)With this definition, we have, for all k P t , ¨ ¨ ¨ , N u , θ k “ Λ¯ δ J N ´ k ´ ÿ j “ e jβ ¯ δ ď Λ¯ δ J e βC Λ T e β ¯ δ ´ , (30)which gives the expected result.We now prove (27). Observe that it is obviously true for the last step N . Assume nowthat it holds true at step k ` , for k ă N , and that (30) holds true for θ k ` . Then,using (25), we have θ k ` j ď , for all j ď J ´ . (31)From Proposition 3, we have ∆ jk ď ¯ δ j ∆ k ` j ÿ (cid:96) “ ¯ δ (cid:96) ´ e ¯ δ } (cid:15) k ` p ˜ X k,j ´ (cid:96)r k ` q} . (32)Using the induction hypothesis (27), we compute ∆ jk ď ¯ δ j ∆ k ` e ¯ δ ´ ¯ δ θ k ` ` e ¯ δ θ k ` j ´ ÿ (cid:96) “ ¯ δ j ´ ´ (cid:96) ›› P ‹ r k ` ,T ` ˜ X k,(cid:96)r k ` ˘›› . (33)We study the last sum. Observe that for (cid:96) P t , ¨ ¨ ¨ , j ´ u , ›› P ‹ r k ` ,T ` ˜ X k,(cid:96)r k ` ˘›› ď ›› P ‹ r k ` ,T ` ˜ X k, r k ` ˘›› ` (cid:96) ÿ i “ ››› P ‹ r k ` ,T ` ˜ X k,ir k ` ˘ ´ P ‹ r k ` ,T ` ˜ X k,i ´ r k ` ˘››› . We observe that P r k ` ,t p ˜ X k, r k ` q “ P r k ,t p ˜ X k, r k q “ P r k ,t p ξ q , for t P r r k ` , T s . Hence, P ‹ r k ` ,T p ˜ X k, r k ` q ď P ‹ r k ,T p ξ q . Also, it is well-checked that there exists a constant C Λ suchthat each P ‹ t,T is C Λ -Lipschitz continuous from L p F t q into L p F T q . Then, j ´ ÿ (cid:96) “ ¯ δ j ´ ´ (cid:96) ›› P ‹ r k ,T ` ˜ X k,(cid:96)r k ` ˘›› ď C Λ j ´ ÿ (cid:96) “ ¯ δ j ´ ´ (cid:96) (cid:96) ÿ i “ ››› ˜ X k,ir k ` ´ ˜ X k,i ´ r k ` ››› ` j ´ ÿ (cid:96) “ ¯ δ (cid:96) ›› P ‹ r k ,T p ξ q ›› . Using (24) in the proof of Proposition 3 and changing the definition of ¯ δ , we obtain j ´ ÿ (cid:96) “ ¯ δ j ´ ´ (cid:96) ›› P ‹ r k ,T ` ˜ X k,(cid:96)r k ` ˘›› ď ¯ δ j ´ ÿ i “ p ∆ ik ` ∆ i ´ k q j ´ ÿ (cid:96) “ i ¯ δ j ´ ´ (cid:96) ` j ´ ÿ (cid:96) “ ¯ δ (cid:96) ›› P ‹ r k ,T p ξ q ›› . (34)Observing that, for all i ď j ´ , ř j ´ (cid:96) “ i ¯ δ j ´ ´ (cid:96) ď ´ ¯ δ , we get j ´ ÿ (cid:96) “ ¯ δ j ´ (cid:96) ›› P ‹ r k ,T ` ˜ X k,(cid:96)r k ` ˘›› ď δ ´ ¯ δ S j ´ k ` ´ ¯ δ ›› P ‹ r k ,T p ξ q ›› , (35)where S nk : “ ř ni “ ∆ ik . Inserting the previous estimate into (33) and changing ¯ δ into δ ,we obtain ∆ jk ď ¯ δ j ∆ k ` e ¯ δ ´ ¯ δ θ k ` ` ` ›› P ‹ r k ,T p ξ q ›› ˘ ` θ k ` ¯ δe ¯ δ ´ ¯ δ S j ´ k . (36)We note that ∆ k ď Λ p ` } P ‹ r k ,T p ξ q} q . Recalling γ in (29), equation (36) leads to ∆ jk ď a j ` γθ k ` ¯ δ S j ´ k . (37)where we set a j : “ p Λ¯ δ j ` γθ k ` qp ` } P ‹ r k ,T p ξ q} q . We have S jk ´ S j ´ k “ ∆ jk ď a j ` γθ k ` ¯ δ S j ´ k , and then S jk ď e γθ k ` ¯ δj S k ` j ÿ (cid:96) “ e γθ k ` ¯ δ p j ´ (cid:96) q a (cid:96) . (38)We compute j ÿ (cid:96) “ a (cid:96) ď ´ jγθ k ` ` Λ¯ δ ´ ¯ δ ¯` ` ›› P ‹ r k ,T p ξ q ›› ˘ , which combined with the properties (31) and (38) leads to, for all j ď J ´ , S jk ď e γ ¯ δ ˆ γ ` Λ1 ´ ¯ δ ˙ ` ` ›› P ‹ r k ,T p ξ q ›› ˘ , where we recall that S k “ ∆ k ď Λ p ` } P ‹ r k ,T p ξ q} q . We insert the previous inequalityinto (37) for j “ J and get ∆ Jk ď ˆ Λ¯ δ J ` ˆ γ ` γ ¯ δe γ ¯ δ p γ ` Λ1 ´ ¯ δ q ˙ θ k ` ˙ ´ ` ›› P ‹ r k ,T p ξ q ›› ¯ . Using (29), this rewrites ∆ Jk ď ´ Λ¯ δ J ` e β ¯ δ θ k ` ¯ ´ ` ›› P ‹ r k ,T p ξ q ›› ¯ , and validates (28) and thus (30). We then obviously have that (27) holds true. l UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 15 Convergence Analysis
Error analysis in the generic case.
We now study the convergence of a genericimplementable solver solver[]() , based upon the local solver solver[](,,) as de-scribed above, as long as the output of the local solver solver[ k ](,,) satisfies someconditions, which are shown to be true for Example 2.In order to define the required assumption, we use the same letters Λ and α as in p H q and p H q , except that, without any loss of generality, we assume that α is greater than1. For the same coefficients as in the equation (6), and in particular for the same driver f , we then ask solver[ k ](,,) to satisfy the following three conditions. p A q sup t P π k ›› U p t, ¯ X t , r ¯ X t sq ´ ¯ Y t ›› α ď e Λ δ ›› U p r k ` , ¯ X r k ` , r ¯ X r k ` sq ´ ¯ Y r k ` ›› α ` Λ max j k ď i ă j k ` ›› ¯ X t i ´ ¯ X t i ›› α ` D p| π |q ` D p| π |q ` ` } ξ } α α ˘ , p A q sup t P π k ›› ¯ X t ´ ¯ X t ›› α ď Λ δ sup t P π k ›› ¯ Y t ´ ¯ Y t ›› α , p A q ›› U p r k ` , ¯ X r k ` , r ¯ X r k ` sq ´ ¯ Y r k ` ›› α α ď Λ ›› U p r k ` , ¯ X r k ` , r ¯ X r k ` sq ´ ¯ Y r k ` ›› α , where p ¯ X, ¯ Y , ¯ Z q : “ solver[ k ]( ¯ X , η , f ) , for f as before, and p ¯ X , ¯ Y , ¯ Z q : “ solver[ k ]( ¯ X , η , f ) ,for another f either equal to f or , are two output values of solver[](,,) associatedto two input processes ¯ X , ¯ X , with the same initial condition ¯ X r k “ ¯ X r k “ ξ , and to twodifferent terminal conditions η and η . For i P t , u , the function D i : r ,
8q Ñ r , isa discretization error associated to the use of the grid π , which satisfies lim h Ó D i p h q “ .Importantly, both D and D are independent of ¯ X , ¯ η , J and N .In full analogy with the discussion right below Theorem 5, we shall also need someconditions on the solver solver[ k ](, , ) used to initialize the algorithm at each step.Following the definition of p P r k ,t q ď t ď T introduced in the statement of Theorem 5, we letby induction, for a given k P t , ¨ ¨ ¨ , N ´ u : P r k ,t p ξ q “ ` solver[ k ]( ξ , , ) ˘ t , t P π k , ξ P L p F r k q , where we recall that ` solver[ k ]( ξ , , ) ˘ is the forward component of the algorithm’soutput, and, for k ď N ´ , P r k ,t p ξ q “ P r (cid:96) ,t ` P r k ,r (cid:96) p ξ q ˘ , t P π (cid:96) , k ă (cid:96) ď N ´ , and then P ‹ r k ,T p ξ q “ max s P π,s Pr r k ,T s | P r k ,s p ξ q| , for ξ P L p F r k q . It then makes sense toassume p A q ›› P ‹ r k ,T p ξ q ´ P ‹ r k ,T p ξ q ›› α ď Λ ›› ξ ´ ξ ›› α p A q ›› P ‹ r k ,T p ξ q ›› α ď Λ ` ` ›› ξ ›› α ˘ where ξ, ξ P L α p F r k q and k P t , ¨ ¨ ¨ , N ´ u . Remark 6.
The main challenging assumption (and maybe the most surprising one) is p A q . It is obviously satisfied when α “ as long as Λ is assumed to be greater than1. We refer to [13] and [17, Chap. 12] for sets of conditions under which this is indeedtrue. When α ą , Assumption p A q is checked provided we have an a priori boundon } U p r k ` , ¯ X r k ` , r ¯ X r k ` sq ´ ¯ Y r k ` } α , see Lemma 10. This permits to invoke the resultproven in our previous paper [20] , which holds true in a weaker setting than the solvabilityresults obtained in [13] and [17, Chap. 12] . Theorem 7.
We can find two constants C Λ , c Λ ą and a continuous non-decreasingfunction B : R ` Ñ R ` matching in , only depending on Λ , such that, for ¯ δ : “ C Λ δ ă min p c Λ , q and β ě B p ¯ δ q satisfying p J ´ q ´ Λ¯ δ J ` e β ¯ δ D p| π |q ¯ e βC Λ T e β ¯ δ ´ ď , (39) where J is the number of Picard iterations in a period, it holds, for any period k Pt , ¨ ¨ ¨ , N u and ξ P L p F r k q , ›› solver[ k ]( ξ ) ´ U p r k , ξ, r ξ sq ›› α ď C ` ¯ δ J ´ ` p N ´ k q D p| π |q ˘ ` ` } ξ } α α ˘ ` C p N ´ k q D p| π |q , for a constant C independent of the discretization parameters. Proof.
The proof will follow closely the proof of Theorem 5 but we now have to takeinto account the discretization error. We will first show that for all k “ t , ¨ ¨ ¨ , N u , ›› (cid:15) k p ξ q ›› α ď θ k ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` ϑ k D p| π |q , (40)where (cid:15) k p ξ q “ solver[ k ]( ξ ) ´ U p r k , ξ, r ξ sq , and p θ k , ϑ k q k “ , ¨¨¨ ,N is defined by the following backward induction: p θ N , ϑ N q : “ p , q ,recall (21), and for k P t , ¨ ¨ ¨ , N ´ u , θ k : “ Λ¯ δ J ` e β ¯ δ t θ k ` ` D p| π |qu and ϑ k : “ e β ¯ δ p ϑ k ` ` q , (41) β being defined as in equation (48).Assume for a while that thids holds true. Then, we have, for all k “ t , ¨ ¨ ¨ , N ´ u , θ k ď ` Λ¯ δ J ` e β ¯ δ D p| π |q ˘ e β ¯ δ p N ´ k q ´ e β ¯ δ ´ and ϑ k ď e β ¯ δ e β p N ´ k q ¯ δ ´ e β ¯ δ ´ . (42)Recalling that ¯ δN “ C Λ T , we get the announced inequality.We now prove (40). Obviously, it holds true for the last step N . Assume now that itis true at step k ` , for k ă N and that (42) holds for θ k ` and ϑ k ` .In particular, using (39), we observe that θ k ` j ď , for all j ď J ´ . (43) First Step.
For j P t , . . . , J u , let ¯∆ jk : “ sup t P π k ›› U p t, ¯ X k,jt , r ¯ X k,jt sq ´ ¯ Y k,jt ›› α . Under p A q´p A q , we will prove in this first step an upper bound for ¯∆ jk , for j “ , ¨ ¨ ¨ , J ,similar to the one obtained in Proposition 3.Using p A q and p H q and the fact that ¯ Y k,jr k ` “ U ` r k ` , ¯ X k,j ´ r k ` , r ¯ X k,j ´ r k ` s ˘ ` (cid:15) k ` p ¯ X k,j ´ r k ` q , we observe that ¯∆ jk ď e Λ δ ”›› U ` r k ` , ¯ X k,jr k ` , r ¯ X k,jr k ` s ˘ ´ U ` r k ` , ¯ X k,j ´ r k ` , r ¯ X k,j ´ r k ` s ˘›› α ` ›› (cid:15) k ` p ¯ X k,j ´ r k ` q ›› α ı ` Λ max j k ď i ă j k ` ›› ¯ X k,jt i ´ ¯ X k,j ´ t i ›› α ` D p| π |q ` D p| π |q ` ` ›› ξ ›› α α ˘ ď C Λ max t P π k ›› ¯ X k,jt ´ ¯ X k,j ´ t ›› α ` e Λ δ ›› (cid:15) k ` p ¯ X k,j ´ r k ` q ›› α (44) ` D p| π |q ` D p| π |q ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ . UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 17
Using p A q , we also have sup t P π k ›› ¯ X k,jt ´ ¯ X k,j ´ t ›› α ď Λ δ sup t P π k ”›› ¯ Y k,jt ´ U p t, ¯ X k,jt , r ¯ X k,jt sq ›› α ` Λ ›› ¯ X k,jt ´ ¯ X k,j ´ t ›› α ` ›› U p t, ¯ X k,j ´ t , r ¯ X k,j ´ t sq ´ ¯ Y k,j ´ t ›› α ı ď C Λ δ ´ ¯∆ jk ` ¯∆ j ´ k ¯ , for δ small enough. Inserting the previous inequality in (44), we get ¯∆ jk ď C Λ δ ¯∆ j ´ k ` e C Λ δ ›› (cid:15) k ` p ¯ X k,j ´ r k ` q ›› α ` D p| π |q ` D p| π |q ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ , ď ¯ δ j ¯∆ k ` e ¯ δ j ´ ÿ (cid:96) “ ¯ δ (cid:96) ›› (cid:15) k ` p ¯ X k,j ´ ´ (cid:96)r k ` q ›› α ` D p| π |q ´ ¯ δ ` D p| π |q ´ ¯ δ ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ , with ¯ δ : “ C Λ δ . We note that compared to (22), there is a new term, namely p D p| π |q ` D p| π |qp ` } P ‹ r k ,T p ξ q} α α q{p ´ ¯ δ q , which is due to the discretization. Second Step.
Using (40) at the previous step k ` and noting that ¯∆ k ď Λ p `} P ‹ r k ,T p ξ q} α q ď p ` } P ‹ r k ,T p ξ q} α α q , we claim that ¯∆ jk ď ` δ j ` γ D p| π |q ˘ ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` γ p ϑ k ` ` q D p| π |q` e ¯ δ θ k ` j ´ ÿ (cid:96) “ ¯ δ j ´ ´ (cid:96) ´ ` ›› P ‹ r k ` ,T p ¯ X k,(cid:96)r k ` q ›› α α ¯ , (45)where γ : “ e ¯ δ {p ´ ¯ δ q .This corresponds to equation (33) adapted to our context. By p A q , we have, for (cid:96) ď J ´ , ›› P ‹ r k ` ,T p ¯ X k,(cid:96)r k ` q ´ P ‹ r k ` ,T p ¯ X k, r k ` q ›› α ď C Λ sup t P π k ›› ¯ X k,(cid:96)t ´ ¯ X k, t ›› α . (46)Using p A q , we then compute, recalling that ¯ Y k, “ , sup t P π k ›› ¯ X k,(cid:96)t ´ ¯ X k, t ›› α ď Λ δ sup t P π k ›› ¯ Y k,(cid:96)t ›› α ď Λ δ ˆ ¯∆ (cid:96)k ` Λ sup t P π k ›› ¯ X k,(cid:96)t ´ ¯ X k, t ›› α ` Λ ` ` ›› ξ ›› α ˘˙ ď C Λ δ ¯∆ (cid:96)k ` C Λ δ ` ` ›› ξ ›› α ˘ , where for the last inequality we used the fact that δ is small enough. Observing that ›› ξ ›› α ď ›› P ‹ r k ,T p ξ q ›› α and combining the previous inequality with (46), we obtain ›› P ‹ r k ` ,T p ¯ X k,(cid:96)r k ` q ´ P ‹ r k ` ,T p ¯ X k, r k ` q ›› α ď C Λ δ ¯∆ (cid:96)k ` C Λ δ ` ` ›› P ‹ r k ,T p ξ q ›› α ˘ . So that, by using the fact that P ‹ r k ` ,T p ¯ X k, r k ` q ď P ‹ r k ,T p ξ q together with a convexityargument, ›› P ‹ r k ` ,T p ¯ X k,(cid:96)r k ` q ›› α α ď ´ C Λ δ ¯∆ (cid:96)k ` ` ` C Λ δ ˘` ` ›› P ‹ r k ,T p ξ q ›› α ˘¯ α , ď ` ` C Λ δ ˘ α ´ ´ C Λ δ ` ¯∆ (cid:96)k ˘ α ` ` ` C Λ δ ˘›› P ‹ r k ,T p ξ q ›› α α ¯ , Appealing to p A q and redefining ¯ δ , we get ›› P ‹ r k ` ,T p ¯ X k,(cid:96)r k ` q ›› α α ď ¯ δ ¯∆ (cid:96)k ` e ¯ δ ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ , which may be rewritten as j ´ ÿ (cid:96) “ ¯ δ j ´ ´ (cid:96) ›› P ‹ r k ` ,T p ¯ X k,(cid:96)r k ` q ›› α α ď ¯ δ j ´ ÿ (cid:96) “ ¯ δ j ´ ´ (cid:96) ¯∆ (cid:96)k ` e ¯ δ ´ ¯ δ ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ . Recalling the notation γ “ e ¯ δ {p ´ ¯ δ q and letting ¯ S nk : “ ř ni “ ¯ δ n ´ i ¯∆ ik , we obtain a newversion of (37), namely ¯∆ jk ď Λ¯ δ j ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` ¯ a ` θ k ` γ ¯ δ ¯ S j ´ k , (47)where we changed the constant in (45) into Λ as we changed the value of ¯ δ , andwhere we put ¯ a “ ` γ θ k ` ` γ D p| π |q ˘` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` γ p ϑ k ` ` q D p| π |q . We straightforwardly deduce that ¯ S jk “ ¯∆ jk ` ¯ δ ¯ S j ´ k ď Λ¯ δ j ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` ¯ a ` ` ` γθ k ` ˘ ¯ δ ¯ S j ´ k ď e γθ k ` j ¯ δ j ¯ S k ` j ´ ÿ (cid:96) “ e γθ k ` (cid:96) ¯ δ (cid:96) ´ Λ¯ δ j ´ (cid:96) ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` ¯ a ¯ , which yields ¯ S jk ď Λ p j ` q ¯ δ j e γθ k ` p j ´ q ` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` ¯ a ´ e γθ k ` ¯ δ , where we used ¯ S k ď p ` } P ‹ r k ,T p ξ q} α α q . Thanks to (47), we get ¯∆ Jk ď Λ¯ δ J ´ ` ¯ δγ p J ` q θ k ` e γθ k ` p J ´ q ¯` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` ¯ a ´ e γθ k ` ¯ δ . Recalling that p J ´ q θ k ` ď , we deduce that, for ¯ δ small enough, ¯∆ Jk ď ` Λ¯ δ J ` e β ¯ δ t θ k ` ` D p| π |qu ˘` ` ›› P ‹ r k ,T p ξ q ›› α α ˘ ` e β ¯ δ p ϑ k ` ` q D p| π |q , provided that β satisfies γ ´ e γθ k ` ¯ δ ď e β ¯ δ . (48)This validates (41) and concludes the proof. l Convergence error for the implemented scheme.
We now analyse the globalerror of our method when the numerical algorithm is given by our benchmark Example2, see Section 4.1.
Lemma 8. (Scheme stability) Condition p A q holds true for the scheme given in Example2. Proof.
For k ď N ´ , we consider p ¯ X, ¯ Y , ¯ Z q : “ solver[ k ]( ¯ X , η , f ) and p ¯ X , ¯ Y , ¯ Z q : “ solver[ k ]( ¯ X , η , f ) with ¯ X r k “ ¯ X r k “ ξ . Letting ∆ X i “ ¯ X t i ´ ¯ X t i and ∆ Y i “ ¯ Y t i ´ ¯ Y t i ,we observe | ∆ X i ` | ď ˇˇˇˇ i ÿ (cid:96) “ j k p t (cid:96) ` ´ t (cid:96) q ∆ b (cid:96) ˇˇˇˇ ` ˇˇˇˇ i ÿ (cid:96) “ j k ∆ σ (cid:96) ∆ ¯ W (cid:96) ˇˇˇˇ , for i P t j k , ¨ ¨ ¨ , j k ` u , where ∆ b (cid:96) : “ b p ¯ X t (cid:96) , ¯ Y t (cid:96) , r ¯ X t (cid:96) , ¯ Y t (cid:96) sq ´ b p ¯ X t (cid:96) , ¯ Y t (cid:96) , r ¯ X t (cid:96) , ¯ Y t (cid:96) sq and,similarly, ∆ σ (cid:96) : “ σ p ¯ X t (cid:96) , r ¯ X t (cid:96) sq ´ σ p ¯ X t (cid:96) , r ¯ X t (cid:96) sq .Invoking Cauchy-Schwartz inequality for the first term and the Bürkholder-Davis-Gundy UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 19 inequality for discrete martingales for the second term and appealing to the Lipschitzproperty of b and σ , we get ›› ∆ X i ` ›› α ď Cδ max (cid:96) “ j k , ¨¨¨ ,i `›› ∆ Y (cid:96) ›› α ` ›› ∆ X (cid:96) ›› α ˘ ` C ›››› i ÿ (cid:96) “ j k | ∆ σ (cid:96) | ¨ | ∆ ˆ W (cid:96) | ›››› α ď Cδ max (cid:96) “ j k , ¨¨¨ ,i `›› ∆ Y (cid:96) ›› α ` ›› ∆ X (cid:96) ›› α ˘ ` C ˆ i ÿ (cid:96) “ j k p t (cid:96) ` ´ t (cid:96) q ›› ∆ X (cid:96) ›› α ˙ ď Cδ max (cid:96) “ j k , ¨¨¨ ,i `›› ∆ Y (cid:96) ›› α ` ›› ∆ X (cid:96) ›› α ˘ ` Cδ { max (cid:96) “ j k , ¨¨¨ ,i `›› ∆ X (cid:96) ›› α ˘ , where we used the identity t (cid:96) ` ´ t (cid:96) “ δ {p j k ` ´ j k q . For δ small enough (taking the supin the sum), we then obtain max j k ď i ď j k ` ›› ∆ X i ›› α ď Cδ max j k ď i ď j k ` ›› ∆ Y i ›› α , (49)which concludes the proof. l We now turn to the study of the approximation error.
Lemma 9.
Assume that p H q - p H q are in force. Then, condition p A q holds true forthe scheme given in Example 2 with D p| π |q ď C a | π | and D p| π |q ď C a | π | . Proof.
First Step.
Given the scheme defined in Example 2, we introduce its piecewisecontinuous version, which we denote by p ¯ X s q ď s ď T . For i ă n , t i ă s ă t i ` , ¯ X s : “ ¯ X t i ` b i p s ´ t i q ` σ i ? s ´ t i (cid:36) i , (cid:36) i : “ ? t i ` ´ t i ∆ ¯ W i , with p b i , σ i q : “ p b p ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i sq , σ p ¯ X t i , r ¯ X t i sqq . In preparation for the proof, we alsointroduce a piecewise càd-làg version, denoted by p ¯ X p λ q s q ď s ď T , where λ is a parameterin r , q . For i ă n , t i ă s ă t i ` , ¯ X p λ q s : “ ¯ X t i ` b i p s ´ t i q ` λσ i ? s ´ t i (cid:36) i . For the reader’s convenience, we also set ¯ U s : “ U ` s, ¯ X s , r ¯ X s s ˘ , ¯ V xs : “ B x U ` s, ¯ X s , r ¯ X s s ˘ , ¯ V µs : “ B µ U ` s, ¯ X s , r ¯ X s s ˘ px ¯ X s yq , ¯ V x, s : “ B x U ` s, ¯ X p q s , r ¯ X s s ˘ . Applying the discrete Itô formula given in Proposition 14, and using the PDE solved by U , recall (1), we compute ¯ U t i ` “ ¯ U t i ` ż t i ` t i ¯ V xs ¨ (cid:32) b ` ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i s ˘ ´ b ` ¯ X t i , ¯ U t i , r ¯ X t i , ¯ U t i s ˘( d s ` ż t i ` t i ˆ E “ ¯ V µs ¨ (cid:32) x b ` ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i s ˘ ´ b ` ¯ X t i , ¯ U t i , r ¯ X t i , ¯ U t i s ˘( y ‰ d s ´ p t i ` ´ t i q f ´ ¯ X t i , ¯ U t i , σ : ` ¯ X t i , r ¯ X t i s ˘ ¯ V xt i , r ¯ X t i , ¯ U t i s ¯ ` ¯ V xt i ¨ ´a t i ` ´ t i σ ` ¯ X t i , r ¯ X t i s ˘ (cid:36) i ¯ ` R wi ` R fi ` R bxi ` R bµi ` R σxi ` R σµi ` δ M p t i , t i ` q ` δ T p t i , t i ` q , with R wi : “ ż t i ` t i p ¯ V x, s ´ ¯ V x, t i q ¨ σ p ¯ X p q t i , r ¯ X t i sq (cid:36) i ? s ´ t i d s , R fi : “ ż t i ` t i ! f ´ ¯ X s , ¯ U s , σ : ` ¯ X s , r ¯ X s s ˘ ¯ V xs , r ¯ X s , ¯ U s s ¯ ´ f ´ ¯ X t i , ¯ U t i , σ : ` ¯ X t i , r ¯ X t i s ˘ ¯ V xt i , r ¯ X t i , ¯ U t i s ¯) d s , R bxi : “ ż t i ` t i ¯ V xs ¨ (cid:32) b ` ¯ X t i , ¯ U t i , r ¯ X t i , ¯ U t i s ˘ ´ b ` ¯ X s , ¯ U s , r ¯ X s , ¯ U s s ˘( d s , R bµi : “ ż t i ` t i ˆ E “ ¯ V µs ¨ (cid:32) x b ` ¯ X t i , ¯ U t i , r ¯ X t i , ¯ U t i s ˘ ´ b ` ¯ X s , ¯ U s , r ¯ X s , ¯ U s s ˘ y (‰ d s , and R σxi “ ż t i ` t i ż ∆ x p s, λ q d λ d s , R σµi “ ż t i ` t i ż ∆ µ p s, λ q d λ d s , (50)where ∆ x p s, λ q : “ Tr ! B xx U ` s, ¯ X p λ q s , r ¯ X s s ˘ a p ¯ X t i , r ¯ X t i sq ´ B xx U ` s, ¯ X s , r ¯ X s s ˘ a p ¯ X s , r ¯ X s sq ) ∆ µ p s, λ q : “ ˆ E ” Tr ! B v B µ U ` s, ¯ X s , r ¯ X s s ˘ px ¯ X p λ q s yqx a p ¯ X t i , r ¯ X t i sqy´ B v B µ U ` s, ¯ X s , r ¯ X s s ˘ px ¯ X s yqx a p ¯ X s , r ¯ X s sqy )ı . Also, δ M p t i , t i ` q is a martingale increment satisfying E “ | δ M p t i , t i ` q| α | F t i ‰ {p α q ď Ch i and } δ T p t i , t i ` q} α ď C Λ h i , recall Proposition 14. Second Step.
Denoting (cid:36) i : “ (cid:36) i {? t i ` ´ t i and δb i : “ h i ż t i ` t i ¯ V xs ¨ (cid:32) b ` ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i s ˘ ´ b ` ¯ X t i , ¯ U t i , r ¯ X t i , ¯ U t i s ˘( d s ` h i ż t i ` t i ˆ E “ ¯ V µs ¨ (cid:32) x b ` ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i s ˘ ´ b ` ¯ X t i , ¯ U t i , r ¯ X t i , ¯ U t i s ˘ y (‰ d s , UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 21 the previous equation reads ¯ U t i ` “ ¯ U t i ` ζ i ` h i ” δb i ´ f ´ ¯ X t i , ¯ U t i , σ : ` ¯ X t i , r ¯ X t i s ˘ ¯ V xt i , r ¯ X t i , ¯ U t i s ¯ ` ¯ V xt i ¨ ` σ p ¯ X t i , r ¯ X t i sq (cid:36) i ˘ı , (51)where ζ i : “ R wi ` R fi ` R bxi ` R bµi ` R σxi ` R σµi ` δ M p t i , t i ` q ` δ T p t i , t i ` q . On the other hand, the scheme can be rewritten as ¯ Y t i “ ¯ Y t i ` ` h i f ` ¯ X t i , ¯ Y t i , ¯ Z t i , r ¯ X t i , ¯ Y t i s ˘ ´ h i ¯ Z t i ¨ (cid:36) i ´ ∆ M i , (52)where ∆ M i satisfies E t i r ∆ M i s “ , E t i r (cid:36) i ¨ ∆ M i s “ and E “ | ∆ M i | ‰ ă 8 . (53)Denoting ∆ ¯ Y i “ ¯ Y t i ´ ¯ U t i , ∆ ¯ Z i “ ¯ Z t i ´ σ : p ¯ X t i , r ¯ X t i sq ¯ V xt i , and adding (51) and (52), weget ∆ ¯ Y i “ ∆ ¯ Y i ` ` h i p δb i ` δf i q ` ζ i ´ h i ∆ ¯ Z i ¨ (cid:36) i ´ ∆ M i , (54)where δf i “ f ` ¯ X t i , ¯ Y t i , ¯ Z t i , r ¯ X t i , ¯ Y t i s ˘ ´ f ´ ¯ X t i , ¯ U t i , σ : ` ¯ X t i , r ¯ X t i s ˘ ¯ V xt i , r ¯ X t i , ¯ U t i s ¯ . For later use, we observe that | δb i | ` | δf i | ď C Λ ` | ∆ ¯ Y i | ` ›› ∆ ¯ Y i ›› ` | ∆ ¯ Z i | ˘ . (55)Summing the equation (54) from i to j k ` ´ , we obtain ∆ ¯ Y i ` j k ` ´ ÿ (cid:96) “ i t h (cid:96) ∆ ¯ Z (cid:96) ¨ (cid:36) (cid:96) ` ∆ M (cid:96) u “ ∆ ¯ Y j k ` ` j k ` ´ ÿ (cid:96) “ i h (cid:96) p δb (cid:96) ` δf (cid:96) q ´ j k ` ´ ÿ (cid:96) “ i ζ (cid:96) . Squaring both sides and taking expectation, we compute, using (53) for the left side andYoung’s and conditional Cauchy-Schwarz inequality for the right side, E t q “ | ∆ ¯ Y i | ‰ ` j k ` ´ ÿ (cid:96) “ i h (cid:96) E t q “ | ∆ ¯ Z (cid:96) | ‰ ď E t q « p ` Cδ q| ∆ ¯ Y j k ` | ` C j k ` ´ ÿ (cid:96) “ i h (cid:96) | δb (cid:96) ` δf (cid:96) | ` Cδ ˆ j k ` ´ ÿ (cid:96) “ i ζ (cid:96) ˙ ff , for i ě q ě j k . Combining (55) and Young’s inequality, this leads to E t q “ | ∆ ¯ Y i | ‰ ` j k ` ´ ÿ (cid:96) “ i h (cid:96) E t q “ | ¯ Z (cid:96) | ‰ ď E t q « e Cδ | ∆ ¯ Y j k ` | ` C j k ` ´ ÿ (cid:96) “ i h (cid:96) | ∆ ¯ Y (cid:96) | ` Cδ ˆ j k ` ´ ÿ (cid:96) “ i ζ (cid:96) ˙ ff . Using the discrete version of Gronwall’s lemma and recalling that ř j k ` ´ (cid:96) “ j k h (cid:96) “ δ , weobtain, for i “ q , | ∆ ¯ Y i | ď E t i « e Cδ | ∆ ¯ Y j k ` | ` Cδ max j k ď i ď j k ` ´ ˆ j k ` ´ ÿ (cid:96) “ i ζ (cid:96) ˙ ff , and then, ∆ Y : “ max j k ď i ď j k ` ›› ∆ ¯ Y i ›› α ď e Cδ ›› ∆ ¯ Y j k ` ›› α ` Cδ ›››› max j k ď i ď j k ` ´ ˆ j k ` ´ ÿ (cid:96) “ i ζ (cid:96) ˙›››› α . (56) Third Step.
To conclude, we need an upper bound for the error } max j k ď i ď j k ` ´ p ř j k ` ´ (cid:96) “ i ζ (cid:96) q} α where ζ (cid:96) is defined in (52). To do so, we study each term in (52) separately. We alsodefine ∆ X : “ max t P π k ›› ¯ X t ´ ¯ X t ›› α and we recall that ¯ X r k “ ξ . Third Step a.
We first study the contribution of R fi to the global error term and notethat ›››› max j k ď i ď j k ` ˆ j k ` ´ ÿ (cid:96) “ i R f(cid:96) ˙›››› α ď C δ | π | j k ` ´ ÿ (cid:96) “ j k ›› R f(cid:96) ›› α . (57)We will upper bound this last term.Let us first observe, that, for t i ď s ď t i ` , | ¯ V xs ´ ¯ V xt i | ď ˇˇ B x U ` s, ¯ X s , r ¯ X s s ˘ ´ B x U ` t i , ¯ X t i , r ¯ X t i s ˘ˇˇ ď C ˆ | ¯ X s ´ ¯ X t i | ` W pr ¯ X s s , r ¯ X t i sq ` h i ` ` | ¯ X t i | ` } ¯ X t i } ˘˙ , where we used the Lipschitz property of B x U given in p H q , together with (8) and (12).Hence, ›› ¯ V xs ´ ¯ V xt i ›› α ď C ´›› ¯ X s ´ ¯ X t i ›› α ` h i ` ` } ¯ X t i } α ˘¯ . (58)From the boundedness of σ and the Lipschitz property of b and U , we compute ›› ¯ X s ´ ¯ X t i ›› α ď C Λ ´ h i ` h i ›› ¯ U t i ´ Y t i ›› α ` h i ›› ¯ X t i ›› α ¯ . (59)Using Lemma 15 from the appendix below, we obtain ›› ¯ V xs ´ ¯ V xt i ›› α ď C ´ h i ` ` } ξ } α ˘ ` h i ∆ Y ¯ . From the boundedness of B x U , σ and the lipschitz property of σ , we obtain ›› σ : ` ¯ X s , r ¯ X s s ˘ ¯ V xs ´ σ : ` ¯ X t i , r ¯ X t i s ˘ ¯ V xt i ›› α ď C ` h i ` ` } ξ } α ˘ ` h i ∆ Y ˘ , where we used the same argument as above to handle the difference between the two σ terms. Combining the previous inequality with the Lipschitz property of f and replicatingthe analysis to handle the difference between the ¯ U terms, we deduce ›› R fi ›› α ď Ch i ´ ∆ X ` h i ` ` } ξ } α ˘ ` h i ∆ Y ¯ . (60) Third Step b.
Combining the Lipschitz property of b , the fact that | ¯ V xs | ` ˆ E r| ¯ V µs | s ď C and Cauchy-Schwarz inequality, we get ›› R bxi ›› α ` ›› R bµi ›› α ď Ch i ›› ¯ U s ´ ¯ U t i ›› α ` ›› ¯ X s ´ ¯ X t i ›› α . (61)Arguing as in the previous step, we easily get ›› R bi ›› α ď Ch i ´ h i ` ` } ξ } α ˘ ` h i ∆ Y ¯ . (62) Third Step c.
We now study the contribution of the terms R wi to the global error. Fromthe independance property of p (cid:36) i q i “ , ¨¨¨ ,n ´ , we may regard each R w(cid:96) as a martingale UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 23 increment. By Burkholder-Davies-Gundy inequalities for discrete martingales, we firstcompute, using the fact that each (cid:36) i is uniformly bounded, ›››› max j k ď i ď j k ´ ˆ j k ` ´ ÿ (cid:96) “ i R w(cid:96) ˙›››› α ď C ›››››› j k ` ´ ÿ (cid:96) “ j k ˇˇˇˇˇż t i ` t i σ : p ¯ X p q t i , r ¯ X t i sq ¯ V ,xs ´ ¯ V ,xt i ? s ´ t i d s ˇˇˇˇˇ ›››››› α ď C j k ` ´ ÿ (cid:96) “ j k h i ˆ h i ` ` } ξ } α ˘ ` ››› sup s Pr t i ,t i ` s | ¯ X p q s ´ ¯ X t i | ››› α ˙ . Since | ¯ X p q s ´ ¯ X t i | ď h i | b p ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i sq| , for s P r t i , t i ` s , so that } ¯ X p q s ´ ¯ X t i } α ď C Λ h i p ` } ¯ X t i } α ` } ¯ Y t i } α q ď C Λ h i p ` } ¯ X t i } α ` ∆ Y q , the previous inequality, togetherwith Lemma 15, leads to ›››› max j k ď i ď j k ´ ˆ j k ` ´ ÿ (cid:96) “ i R w(cid:96) ˙›››› α ď Cδ | π | ´ ` } ξ } α ` | π | ∆ Y ¯ . Similarly, ›››› max j k ď i ď j k ´ ˆ j k ` ´ ÿ (cid:96) “ i δ M p t (cid:96) , t (cid:96) ` q ˙›››› α ď C j k ` ´ ÿ (cid:96) “ j k ››ˇˇ δ M p t (cid:96) , t (cid:96) ` q ˇˇ ›› α ď C δ | π | . Hence, ›››› max j k ď i ď j k ´ ˆ j k ` ´ ÿ (cid:96) “ i R w(cid:96) ˙›››› α ` ›››› max j k ď i ď j k ´ ˆ j k ` ´ ÿ (cid:96) “ i δ M p t (cid:96) , t (cid:96) ` q ˙›››› α ď Cδ | π | ` ` } ξ } α ˘ ` Cδ | π | ∆ Y . (63) Third Step d. (i) We study the contribution of R σxi . We observe that | ∆ x p s, λ q| ď ˇˇ B xx U ` s, ¯ X p λ q s , r ¯ X s s ˘ ´ B xx U ` s, ¯ X s , r ¯ X s s ˘ˇˇ ¨ | a p ¯ X t i , r ¯ X t i sq|` ˇˇ B xx U p s, ¯ X s , r ¯ X s sq ˇˇ ¨ ˇˇ a p ¯ X t i , r ¯ X t i sq ´ a p ¯ X s , r ¯ X s sq ˇˇ , for s P r t i , t i ` s . Using the boundedness and Lipschitz continuity of B xx U and σ , we get,from the previous expression, ›› ∆ x p s, λ q ›› α ď C ´›› ¯ X p λ q s ´ ¯ X s ›› α ` ›› ¯ X s ´ ¯ X t i ›› α ¯ . (64)Observing that ›› ¯ X p λ q s ´ ¯ X s ›› α ď C ? h i , we obtain using (59), for t i ď s ď t i ` ›› ∆ x p s, λ q ›› α ď Ch i ` ` h i ›› ¯ U t i ´ Y t i ›› α ` h i ›› ¯ X t i ›› α ˘ . which leads, using Lemma 15 again, to ›› R σxi ›› α ď Ch i ´ h i ` h i ` ∆ Y ` } ξ } α ˘¯ . (65)(ii) To study R σµi , we first observe that | ∆ µ p s, λ q| ď C ˆ E ”ˇˇ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p λ q s yq ´ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq ˇˇı (66) ` ˆ E “ˇˇ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq ˇˇ ¨ ˇˇ x a p ¯ X t i , r ¯ X t i sq ´ a p ¯ X s , r ¯ X s sqy ˇˇ‰ . For the last term, we combine Cauchy-Schwarz inequality (10) and boundedness andLipschitz continuity of σ to get ˆ E “ˇˇ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq ˇˇ ¨ ˇˇ x a p ¯ X t i , r ¯ X t i sq ´ a p ¯ X s , r ¯ X s sqy ˇˇ‰ ď C ›› ¯ X t i ´ ¯ X s ›› ď C ›› ¯ X t i ´ ¯ X s ›› α . Recalling from (59) that ›› ¯ X s ´ ¯ X t i ›› α ď C Λ p h i ` h i p ∆ Y ` } ¯ X t i } α qq , we obtain, usingLemma 15, that ˆ E “ˇˇ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq ˇˇ ¨ ˇˇ x a p ¯ X t i , r ¯ X t i sq ´ a p ¯ X s , r ¯ X s sqy ˇˇ‰ ď C Λ h i ˆ ` h i t ∆ Y ` } ξ } α u ˙ . (67)For the first term in (66), we use p H q equation (13) to get |B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p λ q s yq ´ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq|ď C ! ` |x ¯ X p λ q s y| α ` |x ¯ X s y| α ` ›› ¯ X s ›› α ) |x ¯ X p λ q s y ´ x ¯ X s y| . By Cauchy Schwarz inequality, we obtain ˆ E ” |tB υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p λ q s yq ´ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yqu| ı ď C a h i ´ ` } ¯ X p λ q s } α α ` } ¯ X s } α α ¯ . (68)We then observe that ›› ¯ X p λ q s ›› α ` ›› ¯ X s ›› α ď C ´›› ¯ X t i ›› α ` h i ›› ¯ U t i ´ ¯ Y t i ›› α ` a h i ¯ ď C p ` } ξ } α ` δ ∆ Y q , where we used lemma 15 for the last inequality. Combining the last inequality with (68)and using also (67), we compute | R σµi | ď Ch i p ` } ξ } α ` δ ∆ Y q , and then ›››› j k ` ´ ÿ (cid:96) “ j k | R σµ(cid:96) | ›››› α ď C | π | δ ` ` δ ∆ Y ` } ξ } α ˘ . (69)4. Collecting the estimates (60), (62) and (65), we compute ¨˝ j k ` ´ ÿ (cid:96) “ j k ›› R f(cid:96) ` R b(cid:96) ` R σx(cid:96) ›› α ˛‚ ď Cδ ` ∆ X ` | π |t ` } ξ } α u ` | π | ∆ Y ˘ . Observing that ¨˝ j k ` ´ ÿ (cid:96) “ j k ›› δ T p t i , t i ` q ›› α ˛‚ ď Cδ | π | , and combining the previous inequality with (69), (63) and (56), we obtain ∆ Y ď e Cδ } ¯ U r k ` ´ ¯ Y r k ` } α ` C ´ δ ∆ X ` | π | ` ` } ξ } α ˘ ` | π | δ ∆ Y ¯ , which concludes the proof for δ small enough. l UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 25
Lemma 10.
Assume that g and f p¨ , , , r¨ , sq are bounded. Then p A q is satisfied what-ever the value of α . Proof.
It suffices to prove that U is bounded on the whole space and that ¯ Y is boundedindependently of the discretization parameters.We refer to [20] for the proof of the boundedness of U .The bound for ¯ Y may obtained by squaring (52) and then by taking the conditionalexpectation exactly as done in the second step of the proof of Lemma 9. l Assumptions p A q and p A q are easily checked. It suffices to observe that p P r k ,t p ξ qq t P π,t ě r k coincides with the solution of the discrete Euler scheme: ¯ X t i ` “ ¯ X t i ` p t i ` ´ t i q b ` X t i , , r X t i , s ˘ ` a t i ` ´ t i σ ` X t i , r X t i s ˘ (cid:36) i , with ¯ X r k “ ξ as initial condition.Combining Lemma 9, Lemma 8 and Lemma 10 with Theorem 7, we have the followingresult. Corollary 11.
Under p H q - p H q , assuming (39) , the following holds ›› solver[ k ]( ξ ) ´ U p r k , ξ, r ξ sq ›› α ď C ´ p Cδ q J ´ ` | π | δ ´ p ` } ξ } α q ¯ , for ¯ δ small enough. The first term in the right hand side is connected with the local Picard iterationson a step of length δ . As expected, it decreases geometrically fast with the number ofiterations. The second term is due to the propagation of the error along the mesh. Theleading term | π | is consistent with that observed for classical forward-backward systems,see for instance [21, 22]. The normalization by δ is due to the propagation of the errorthrough the successive local solvers.4. Numerical applications
In practice, we would like to approximate the value of U p , ¨q at some point p x, µ q P R d ˆ P p R d q . In the first section below, we explain how to retrieve such approximationusing the approximation of U p , ξ, r ξ sq given by the algorithm solver[ ]() , for some ξ „ µ . In a second part, we discuss the numerical results obtained by implementing solver[ ]() with two levels, i.e. N “ . In particular, we show that it is more efficientthan an algorithm based simply on Picard Iterations.4.1. Approximation of U p , x, µ q . The goal of this section is to show how to obtainan approximation of U p , x, r ξ sq with ξ „ µ and x P supp p µ q . We will assume that wethus have at hand a discrete valued random variable ξ | π | „ µ | π | “ ř M(cid:96) “ p (cid:96) δ x(cid:96) such that µ | π | is a good approximation of µ for the Wasserstein distance. For instance, such anapproximation can be constructd by using quantization techniques. Then, we can use solver[ ]( ξ | π | ) to obtain an approximation of U p , ξ | π | , r ξ | π | sq .Note that solver[ ]( ξ | π | ) is a discrete random variable as the algorithm is initialisedby a discrete random variable as well. In practice, this means that each point x (cid:96) will bethe root of a tree and will be associated to an output value y (cid:96) “ U p , x (cid:96) , r ξ | π | sq and then solver[ ]( ξ | π | ) „ ř M(cid:96) “ p (cid:96) δ y (cid:96) . It is important to remark that the computations on thetrees are connected via the McKean-Vlasov interaction. Using the Lipschitz continuity of U , one easily obtains | U p , x, µ q ´ U p , x ¯ (cid:96) , µ | π | q| ď C ˆ min y P supp pr ξ | π | sq | y ´ x | ` W p µ | π | , µ q ˙ “ : E p| π | , ξ q , (70)where x ¯ (cid:96) is a point in the support of µ | π | realising the minimum in the first line. Remark 12.
In many cases, it will be easy to have x P supp p µ | π | q and thus reduce theabove error to the term W p µ | π | , µ q . This is obviously the case if ξ is deterministic. As mentioned above, the approximation of U p , x ¯ (cid:96) , µ | π | q is obtained by running solver[ ]( ξ | π | ) and by taking its value on the tree initiated at x ¯ (cid:96) , precisely we have U p , x ¯ (cid:96) , µ | π | q “ y ¯ (cid:96) .The corresponding pointwise error is given by E p| π | , δ, ξ q : “ | y ¯ (cid:96) ´ U p , x ¯ (cid:96) , r ξ | π | sq| . (71)Of a course, this might be estimated by E p| π | , δ, ξ q ď p ¯ (cid:96) ››› U p , ξ | π | , r ξ | π | sq ´ solver[ ]( ξ | π | ) ››› , but this is very poor when the initial distribution µ is diffuse and accordingly when µ | π | has a large support, in which case p ¯ (cid:96) is expected to be small.To bypass this difficulty, we must regard E p| π | , δ, ξ q as a conditional error. Somehow,it is the error of the numerical scheme conditional on the initial root of the tree. It requiresa new analysis, but it should not be so challenging: Now that we have investigated theerror for the McKean-Vlasov component, we can easily revisit the proof of Theorem 7 inorder to derive a bound for this conditional error.Instead of revisiting the whole proof, we can argue by doubling the variables. For ξ and x as above, we can regard the four equations (2), (3), (4) and (5) as a single forward-backward system of the McKean-Vlasov type. The forward component of such a doubledsystem is X “ p X ,x,µ , X ,ξ q and the backward components are Y “ p Y ,x,µ , Y ,ξ q and Z “ p Z ,x,µ , Z ,ξ q . Except for the fact that the dimension of X is no longer equal to thedimension of the noise, which we assumed to be true for convenience only, and for thefact that Y takes values in R , the setting is exactly the same as before, namely p X , Y , Z q can be regarded as the solution of a McKean-Vlasov forward-backward SDE in whichthe mean field component reduces to the marginal law of p X ,ξ , Y ,ξ q . We observe inparticular that Y ,x,µt “ U p t, X ,x,µt , r X ,ξt sq , Y ,ξt “ U p t, X ,ξt , r X ,ξt sq , t P r , T s , with similar relationships for Z ,x,µ and Z ,ξ . Hence, Y t (and Z t ) can be represented asa function of X , which was the key assumption in our analysis. For sure, the fact that Y takes values in dimension 2 is not a limitation for duplicating the arguments used toprove Theorem 7.Numerically speaking, the tree initiated at root x ¯ (cid:96) under the initial distribution µ | π | provides an approximation of U p , x ¯ (cid:96) , r ξ | π | sq , which is equal to Y ,x ¯ (cid:96) , r ξ | π | s . So our numer-ical (implemented) scheme is in fact a numerical for the whole process p X , Y , Z q .This leads us to the following result. Theorem 13.
Let y ¯ (cid:96) be the approximation of U p , x, µ q obtained by calling solver[ ]( ξ | π | ) ,where ¯ (cid:96) is defined in (70) . Then, the following holds | U p , x, µ q ´ y ¯ (cid:96) | ď E p| π | , ξ q ` E p| π | , δ, ξ q , UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 27 where E p| π | , δ, ξ q can be estimated by Corollary 11, with p ` } ξ } α q replaced by p `| x ¯ (cid:96) | ` } ξ } α q . Numerical illustration.
In this section, we will prove empirically the convergenceof the approximation obtained by the solver solver[]() . In particular, we will comparethe output of our algorithm solver[]() , when implemented with two levels, i.e. N “ (we simply call it two-level algorithm ), with the output of a basic algorithm based onlyon Picard iterations, which can be seen as a solver solver[]() , but with only one level,i.e. N “ (we simply call it one-level algorithm ). In both cases, we use Example 2 asdiscretization scheme, with a standard Bernoulli quantization of the normal distribution, d being equal to . In the numerical studies below, we show that the two-level algorithm converges in case when the one-level algorithm fails.4.2.1. The example of a linear model.
In this part, we compare the output of both algo-rithms for the following linear model where a closed-form solution is available: d X t “ ´ ρ E r Y s t d t ` σ d W t , X “ x , d Y t “ ´ aY t d t ` Z t d W t , and Y T “ X T , for ρ, a ą , and the true solution for E r X s “ m is given by Y “ m e aT ` ρa p e aT ´ q . The errors for various time steps and for both algorithms are shown on the log-log errorplot of Figure . The parameters are fixed as follows: ρ “ . , a “ . , σ “ , T “ and x “ . Moreover, the two-level algorithm uses Picard Iterations per level, and the one-level algorithm computes Picard Iterations.4.2.2.
Efficiency of the solver[]() algorithm.
In this section, we compare the two-levelalgorithm and the one-level algorithm on two models, for which existence and uniquenessto the master equation (or the FBSDE system) hold true for any arbitrary terminal time T and Lipschitz constant L of the coefficients function. Nevertheless, as stated in thetheorems above, the convergence of the algorithms is guaranted only for a periods oftime which are controlled by L and T . Here, we fix the terminal date T and allow L to vary with the use of a coupling parameter ρ , see equations (72) (for a case withoutMcKean-Vlasov interaction) and (73) (for a case with McKean-Vlasov interaction). Wewill see below that, as expected, the two-level algorithm converges for a larger range ofcoupling parameter than the one-level algorithm .An example with no McKean-Vlasov interaction. .Here, the model is the following d X t “ ρ cos p Y t q d t ` σ d W t , X “ x ,Y t “ E t r sin p X T qs . (72)On Figure , we plot the output of the two-level and one-level algorithm along witha proxy of the true solution computed by usual BSDE approximation method (after aGirsanov transform) and with a very high-level of precision. On the graph, the value Y stands for the approximation of U p , x q : There is no dependence upon the initial measureas there is no MKV interaction in this example. The parameters are fixed as follows: σ “ , T “ and x “ . Moreover, the two-level algorithm uses Picard Iterations perlevel, and the one-level algorithm computes Picard Iterations.
Figure 1.
Convergence of the algorithms: log-log error plot for the samedata as in the text. We can observe that both algorithms return the samevalue which is close to the true value. This validates the convergence ofboth methods in this simple linear setting.An example from large population stochastic control. .For this part, the model is given by d X t “ ´ ρY t d t ` d W t , X “ x , d Y t “ atan p E r X t sq d t ` Z t d W t and Y T “ G p X T q : “ atan p X T q . (73)coming from Pontryagin principle applied to MFG inf α E „ G p X αt q ` ż T ˆ ρ α t ` X αt atan p E r X αt sq ˙ d t with d X αt “ α t d t ` d W t , see e.g. [14].We do not know the exact solution for this model and it is not possible to obtain easilyan approximation as in the previous example. We plot on Figure , the output valueof the one-level algorithm and two-level algorithm . On the graph, the value Y standsfor the approximation of U p , x, δ x q . The parameters are fixed as follows: σ “ , T “ and x “ . Moreover, the two-level algorithm uses Picard Iterations per level, and the one-level algorithm computes Picard Iterations.5.
Appendix
A discrete Itô formula.
We consider the following Euler scheme on the discretetime grid π of the interval r , T s , recall (20), ¯ X t i ` “ ¯ X t i ` b i p t i ` ´ t i q ` σ i a t i ` ´ t i (cid:36) i , (74) UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 29
Figure 2.
Comparison of algorithms’ output for different value of thecoupling parameter and for the same data as in Example (72): two-level(black star), one-level (blue cross), true value (red line). The two-levelalgorithm converges for larger coupling parameter than the one-level al-gorithm . It is close to the true solution up to parameter ρ “ , thediscrepancy for large coupling parameter coming most probably from thediscrete-time error. Interestingly, the one-level algorithm shows bifurca-tions.where p (cid:36) i q i ď n are i.i.d. centered R d -valued random variables such that the covariancematrix E r (cid:36) i (cid:36) : i s is the identity matrix and } (cid:36) i } α ď Λ h i , and p b i , σ i q P L p F t i q , for all i ď n .We also introduce a piecewise continuous version of the previous scheme, for i ă n , t i ď s ă t i ` and λ P r , s , the process p ¯ X p λ q t q ď t ď T , ¯ X p λ q s “ ¯ X t i ` b i p s ´ t i q ` σ i λ ? s ´ t i (cid:36) i (75)and ¯ X p λ q t n “ ¯ X t n . Following the notation used in the proof of Lemma 9, we just write p ¯ X s q ď s ď T for p ¯ X p q s q ď s ď T , which defines a continuous version of the Euler scheme givenin (74). Figure 3.
Algorithms’ output for the same data as in Example (73): one-level algorithm (blue line), two-level algorithm (black line). We observethe same phenomenon as in the previous model: The two-level algorithm converges to a unique value for a larger range of coupling parameter thanthe one-level algorithm , which exhibits a bifurcation. Observe that the two-level algorithm fails to converge at some points: One should add alevel of computation to shorten the time period δ . Proposition 14.
For any i P t , ¨ ¨ ¨ , n ´ u , the following holds true: U p t i ` , ¯ X t i ` , r ¯ X t i ` sq “ U p t i , ¯ X t i , r ¯ X t i sq ` ż t i ` t i B t U p s, ¯ X s , r ¯ X s sq d s ` ż t i ` t i ˆ B x U p s, ¯ X s , r ¯ X s sq ¨ b i ` ż Tr “ B xx U p s, ¯ X p λ q s , r ¯ X s sq a i ‰ d λ ˙ d s ` ż t i ` t i ˆ E “ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq ¨ x b i y ‰ ds ` ż ˆ E ” Tr “ B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p λ q s yqx a i y ‰ d λ ı d s ` ż t i ` t i B x U p s, ¯ X p q s , r ¯ X s sq σ i (cid:36) i ? s ´ t i d s ` δ M p t i , t i ` q ` δ T p t i , t i ` q , where a i is here equal to σ i σ : i , and δ M p t i , t i ` q is a martingale increment satisfying } δ M p t i , t i ` q} α ď C Λ h i and } δ T p t i , t i ` q} α ď C Λ h i . Proof.
By writing ¯ X t i ` “ ¯ X t i ` ż t i ` t i ` b i ` σ i (cid:36) i ? s ´ t i ˘ ds, UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 31 and by using the standard chain rule for continuously differentiable functions on a Hilbertspace, we get U p t i ` , ¯ X t i ` , r ¯ X t i ` sq “ U p t i , ¯ X t i , r ¯ X t i sq ` ż t i ` t i B t U p s, ¯ X s , r ¯ X s sq d s ` ż t i ` t i ˆ B x U p s, ¯ X s , r ¯ X s sq ¨ ` b i ` σ i (cid:36) i ? s ´ t i ˘ d s ` ˆ E „ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq ¨ x b i ` σ i (cid:36) i ? s ´ t i y ˙ d s . Now we observe that, B x U p s, ¯ X s , r ¯ X s sq “ B x U p s, ¯ X p q s , r ¯ X s sq ` ? s ´ t i ż B xx U p s, ¯ X p λ q s , r ¯ X s sq σ i (cid:36) i d λ “ B x U p s, ¯ X p q s , r ¯ X s sq ` ? s ´ t i B xx U p s, ¯ X p q s , r ¯ X s sq σ i (cid:36) i ` ? s ´ t i T p s q , where T p s q is a random variable defined on p Ω , F , P q such that } T p s q} α ď Ch i , and B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X s yq “ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p q s yq` ? s ´ t i ż B υ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p λ q s yqx σ i (cid:36) i y d λ “ B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p q s yq` ? s ´ t i B v B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p q s yqx σ i (cid:36) i y ` ? s ´ t i T p s q , where T p s q is a random variable on the enlarged space p Ω ˆ ˆΩ , F b ˆ F , P b ˆ P q such that ˆ E “ | T p s q| α ‰ {p α q ď Ch i .We insert these expansions back into the identity we obtained for the term U p t i ` , ¯ X t i ` , r ¯ X t i ` sq .We let δ M p t i , t i ` q “ ż t i ` t i ” B xx U p s, ¯ X p q s , r ¯ X s sq σ i (cid:36) i ¨ ` σ i (cid:36) i ˘ ´ E t i ” B xx U p s, ¯ X p q s , r ¯ X s sq σ i (cid:36) i ¨ ` σ i (cid:36) i ˘ıı d s ,δ T p t i , t i ` q “ ż t i ` t i ` T p s q ` T p s q ˘ ¨ σ i (cid:36) i d s . It defines a martingale increment satisfying E t i “ | δ M p t i , t i ` q| α ‰ {p α q ď Ch i . Observingthat for t i ď s ď t i ` , ˆ E ” B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p q s yq ¨ x σ i (cid:36) i y ı “ , E t i ” B xx U p s, ¯ X p q s , r ¯ X s sq σ i (cid:36) i ¨ ` σ i (cid:36) i ˘ı “ E t i ” Tr ` B xx U p s, ¯ X p q s , r ¯ X s sq a i ˘ı , E t i ” B xx U p s, ¯ X p q s , r ¯ X s sq σ i (cid:36) i ¨ ` σ i (cid:36) i ˘ı “ E t i ” Tr ` B xx U p s, ¯ X p q s , r ¯ X s sq a i ˘ı , ˆ E ” B v B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p q s yq? s ´ t i x σ i (cid:36) i y ¨ x σ i (cid:36) i y ı “ ˆ E ” Tr ` B v B µ U p s, ¯ X s , r ¯ X s sqpx ¯ X p q s yqx a i y ˘ı , we complete the proof. l Estimates for the scheme given in Example 2.Lemma 15.
Under p H q - p H q , the following holds for the forward component of thescheme given in Example 2 and its continuous version, max t P π k ›› ¯ X t ›› α ď C Λ ˆ ` ›› ¯ X r k ›› α ` δ max t P π k ›› U p t, ¯ X t , r ¯ X t sq ´ Y t ›› α ˙ , (76) Proof.
We introduce d i : “ | U p t i , ¯ X t i , r ¯ X t i sq ´ ¯ Y t i | and observe from the Lipschitz prop-erty of b and U that ˇˇ b ` ¯ X t i , ¯ Y t i , r ¯ X t i , ¯ Y t i s ˘ˇˇ ď C Λ ` ` | ¯ X t i | ` ›› ¯ X t i ›› α ` d i ` ›› d i ›› α ˘ . (77)Recall that the scheme for the forward component reads ¯ X t i ` “ ¯ X r k ` i ÿ (cid:96) “ j k b ` ¯ X t (cid:96) , ¯ Y t (cid:96) , r ¯ X t (cid:96) , ¯ Y t (cid:96) s ˘ p t (cid:96) ` ´ t (cid:96) q ` i ÿ (cid:96) “ j k σ ` ¯ X t (cid:96) , r ¯ X t (cid:96) s ˘ ∆ ¯ W (cid:96) . Squaring the previous inequality, using Cauchy-Schwarz inequality for the first sum andthe martingale property for the second sum, we obtain ›› ¯ X t i ` ›› α ď C ›› ¯ X r k ›› α ` C i ÿ (cid:96) “ j k h (cid:96) ´ δ ›› b p ¯ X t (cid:96) , ¯ Y t (cid:96) , r ¯ X t (cid:96) , ¯ Y t (cid:96) sq ›› α ` } σ p ¯ X t (cid:96) , r ¯ X t (cid:96) sq} α ¯ , where we used again Bürkholder-Davis-Gundy inequality for discrete martingales.Combining (77) with the boundedness of σ , we then have ›› ¯ X t i ` ›› α ď C ´›› ¯ X r k ›› α ` δ ` δ max j k ď i ă j k ` } d i } α ` Cδ i ÿ (cid:96) “ j k h (cid:96) ›› ¯ X t (cid:96) ›› α ¯ . Using the discrete version of Gronwall’s lemma, the result easily follows. l References [1]
Y. Achdou and I. Capuzzo-Dolcetta (2010) Mean field games: numerical methods
SIAM J.Numer. Anal , 48, pp. 1136-1162.[2]
Y. Achdou, F. Camilli and I. Capuzzo-Dolcetta (2013) Mean field games: convergence ofa finite difference method,
SIAM J. Numer. Anal. , 51, pp. 2585-2612.[3]
Y. Achdou and A. Porretta (2016) Convergence of a Finite Difference Scheme to WeakSolutions of the System of Partial Differential Equations Arising in Mean Field Games,
SIAM J.Numer. Anal. , 54, pp. 161-186.[4]
S. Alanko
Regression-based Monte Carlo methods for solving nonlinear PDEs. , PhD disserta-tion, New York University, 2015.[5]
Bayraktar, E. and Budhiraja, A. and Cohen, A. (2016)
Rate Control under Heavy Trafficwith Strategic Servers preprint , http://arxiv.org/abs/1605.09010 [6] J.D. Benamou and G. Carlier (2015) Augmented Lagrangian Methods for Transport Opti-mization, Mean Field Games and Degenerate Elliptic Equations,
Journal of Optimization Theoryand Applications , 167, pp. 1-26.[7]
C. Bender and J. Zhang (2008) Time discretization and Markovian iteration for coupledFBSDEs,
Ann. Appl. Probab. , 18 (1), pp. 143-177.[8]
A. Bensoussan, J. Frehse and P. Yam
Mean Field Games and Mean Field Type ControlTheory , Springer Verlag, 2013.[9]
A. Bensoussan, J. Frehse and P. Yam (2015) The Master equation in mean field theory,
Journal de Mathématiques Pures et Appliquées , 103, pp. 1441-1474.[10]
A. Bensoussan, J. Frehse and P. Yam (2016) On the interpretation of the master equation,
Stochastic Processes and their Applications , doi = "10.1016/j.spa.2016.10.004"[11]
R. Buckdahn, J. Li, S. Peng and C. Rainer (2014) Mean-field stochastic differential equa-tions and associated PDEs, forthcoming in
Annals of Probability . UMERICAL METHOD FOR FBSDES OF MCKEAN-VLASOV TYPE 33 [12]
P. Cardaliaguet (2012) Notes from P.L. Lions’ lectures at the Collège de France, notes , „ cardalia/MFG100629.pdf .[13] P. Cardaliaguet, F. Delarue, J.-M. Lasry and P.-L. Lions (2015) The master equationand the convergence problem in mean field games, preprint , http://arxiv.org/abs/1509.02505 .[14] R. Carmona and F. Delarue (2013) Probabilistic Analysis of Mean Field Games,
SIAMJournal on Control and Optimization , 51 (4), pp. 2705-2734.[15]
R. Carmona, and F. Delarue (2013) Mean Field Forward-Backward Stochastic DifferentialEquations,
Electronic Communications in Probability .[16]
R. Carmona and F. Delarue
The master equation for large population equilibriums, In:Stochastic Analysis and Applications, Springer Verlag, pp. 77 Ð 128, 2014.[17]
R. Carmona and F. Delarue
Probabilistic Theory of Mean Field Games: vol. I, Mean FieldFBSDEs, Control, and Games , Springer Verlag, 2017.[18]
R. Carmona and F. Delarue
Probabilistic Theory of Mean Field Games: vol. II, Mean FieldGames with Common Noise and Master Equations , Springer Verlag, 2017.[19]
R. Carmona, and F. Delarue, and A. Lachapelle (2013) Control of McKean-Vlasov versusMean Field Games,
Mathematics and Financial Economics , 7, pp. 131-166.[20]
J.-F. Chassagneux, D. Crisan and F. Delarue (2014) A Probabilistic approach to classicalsolutions of the master equation for large population equilibria, preprint .[21]
F. Delarue and S. Menozzi (2006) A forward-backward stochastic algorithm for quasi-linearPDEs,
Ann. Appl. Probab. , 16 (1), pp 140-184.[22]
F. Delarue and S. Menozzi (2008). An Interpolated Stochastic Algorithm for Quasi-LinearPDEs.
Mathematics of Computation , 77, pp. 125–158.[23]
D.A. Gomes, L. Nurbekyan and E. Pimentel , Economic Models and Mean-field GamesTheory , Publicaões Matemáticas, IMPA, Rio, Brazil, 2015.[24]
D.A. Gomes, E. Pimentel and V. Voskanyan
Regularity Theory for Mean-Field Game Sys-tems , Springer International Publishing Switzerland, 2016.[25]
O. Guéant (2012) New numerical methods for mean field games with quadratic costs,
Networksand Heterogeneous Media , 2, pp. 315-336.[26]
A. Lachapelle, J. Salomon and G. Turinici (2010) Computation of mean field equilibria ineconomics,
Mathematical Models and Methods in Applied Sciences , 20, pp. 567-588.[27]
J.M. Lasry and P.L. Lions (2006) Jeux à champ moyen I. Le cas stationnaire,
Comptes Rendusde l’Académie des Sciences de Paris, ser. A , 343 (9).[28]
J.M. Lasry and P.L. Lions (2006) Jeux à champ moyen II. Horizon fini et contrôle optimal,
Comptes Rendus de l’Académie des Sciences de Paris, ser. A , 343 (10).[29]
J.M. Lasry and P.L. Lions (2007) Mean Field Games,
Japanese Journal of Mathematics , 2(1), pp. 229-260.[30]
P.L. Lions (2014) Estimées nouvelles pour les équations quasil-inéaires
Seminar in Applied Mathematics at the Collège de France , [31] J. Ma, H. Yin and J. Zhang (2012) On non-Markovian forward-backward SDEs and backwardstochastic PDEs.
Stochastic Processes and their Applications , 122, pp. 3980-4004.[32]
J. Ma, P. Protter and J. Yong (1994) Solving forward-backward stochastic differential equa-tions explicitly – a four step scheme.
Probab. Theory Related Fields , 98, pp. 339-359.[33]
H. Pham, and X. Wei (2015) Bellman equation and viscosity solutions for mean field stochasticcontrol problem, preprint ,,