A numerical scheme for a mean field game in some queueing systems based on Markov chain approximation method
Erhan Bayraktar Amarjit Budhiraja Asaf Cohen


On the repeated inversion of a covariance matrix
    M. de Jong

  NWO-I, Nikhef, PO Box 41882, Amsterdam, 1098 DB Netherlands 

  Leiden University, Leiden Institute of Physics, PO Box 9504, Leiden, 2300 RA Netherlands

    
===========================================================================================================================================================================




We use the Markov chain approximation method to construct approximations for the solution of the mean field game (MFG) with reflecting barriers studied in <cit.>. The  MFG is formulated in terms of a controlled reflected diffusion with a cost function that depends on the reflection terms in addition to the standard variables: state, control, and the mean field term.  This MFG arises from the asymptotic analysis of an N-player 
game for single server queues with strategic servers. By showing that our scheme is an almost contraction, we establish the convergence of this numerical scheme over a small time interval. 




AMS Classification: 65M12, 60K25, 91A13, 60K35, 93E20,

65M12, 60F17.



Keywords: Numerical scheme, mean field games, Nash equilibrium, rate control, reflected diffusions, heavy traffic limits, queuing systems, Markov chain approximation method.






§ INTRODUCTION



The theory of mean field games (MFGs) was initiated a decade ago in the seminal work of Lasry and Lions <cit.>, and Huang, Malhamé, and Caines <cit.>. For theoretical study and applications of this theory see <cit.> and the references therein. MFGs are control problems that approximate many player games with weak interaction between the players that is given in terms of the empirical distribution of the players' states. 

In these control problems the empirical distribution that governs the interaction is replaced by a deterministic flow of measures. A solution of the MFG is a probability measure on the path space of the single player state 
that is the distribution of the state process under the optimal control for the control problem associated with the 
flow of measures given by the (time-)marginal distributions of this probability measure.



A standard (probablistic) method to prove the existence of a MFG solution is by solving a fixed point theorem on the space of probability measures on certain path spaces. 
A probability measure on the path space is fixed and a stochastic control problem is formulated in terms of the flow of time marginals of this probability measure.
Then a `best reply' to the probability measure  is found by solving this  control problem. The distribution of the state process under the best reply is another measure on the path space. 
Solution of the MFG is the fixed point of this map that takes a probability measure to its implied `best reply' distribution.
There are other ways to describe a MFG solution, for example the seminal papers of Lasry and Lions <cit.> represent a MFG solution through two coupled nonlinear partial differential equations; one is an equation of Hamilton-Jacobi-Bellman (HJB) type while the second takes the form of a Kolmogorov forward equation, and  recent works of Carmona, Delarue, and Lacker <cit.>, using probabilistic methods,  characterize the MFG solution as a solution to certain forward backward stochastic differential equations. 


In general closed form solutions for MFGs are not available and thus one needs numerical approximations. In our work we study one such procedure that uses the Markov chain approximation method (<cit.>) and establish convergence of the scheme over a small time interval.

In recent years there have been several works on numerical schemes for MFGs, most of which are based on the PDE system of <cit.>. 
Achdou and Capuzzo-Dolcetta <cit.> were the first to suggest a finite difference method for approximating the PDE system relying on monotone approximations of the Hamiltonians and a weak formulation of the forward equation. Together with Camilli, the same authors proved in <cit.> the convergence of the scheme. In <cit.>, Achdou and Porretta showed that the solutions of a certain discrete system converges to a weak solution of the PDE system. In <cit.>, Lachapelle, Salomon, and Turinici provided an iterative scheme using a discrete Markov decision problem. Taking advantage of the structure of the problem (in particular, the problem is linear-quadratic in the control), they used the monotonic algorithm method introduced in <cit.> and iteratively constructed a value function, control, and a measure by using finite differences based on the forward-backward system. Guéant studies numerical schemes when the Hamiltonians are quadratic, see <cit.>. Semi-Lagrangian schemes were studied by Carlini and Silva in <cit.>. In a recent paper, Chassagneux, Crisan, and Delarue <cit.> used the master equation and by making smoothness assumptions on this infinite dimensional PDE, they proposed an algorithm based on Picard iterations and the continuation method.  
The master equation is a parabolic partial differential equation with a terminal condition. Its variables are time, state, and measure and its solution approximates the value function of the MFG, see e.g., <cit.>. 

Our method in contrast to above methods is purely probabilistic. We do not make smoothness assumptions as in <cit.>. We use an iterative Markov chain approximation method (see <cit.>) to construct numerical solutions of the MFG. Specifically, we discretize time and space and for a fixed measure on the path space we define a Markov decision problem that is suggested by the MFG. In the first step of the iteration, the law of the solution of the MDP is computed. Then we take this law as the starting point to formulate the MDP for the second iteration and repeat the process. Unfortunately, it is not clear that the map defined by such iterations is in general a contraction. We instead show that the map is an almost contraction over a small time interval with length independent of the discretization parameter. By an almost contraction we roughly mean that the map is a contraction
up to  an additional term  that vanishes as the discretization parameter approaches 0. 
The proof of this almost contraction property relies on the construction of a coupling between certain controlled reflected Markov chains (see proof of Proposition <ref>) which we believe is of independent interest. 
Using the above almost contraction property, tightness of relevant processes, and weak convergence arguments, we show the convergence of the laws obtained from the iteration scheme to the solution of the MFG over a small time interval.
Proving the convergence of a Markov chain based approximation method of the form considered in this work over an arbitrary time interval is for now a challenging open problem. 









The paper is organized as follows. In Section <ref> we present the MFG and summarize the results from <cit.>. In Section <ref> we provide the numerical scheme and present our main convergence result (Theorem <ref>). Section <ref> provides proofs of some auxiliary results from Section <ref>. Finally Section <ref> gives a numerical example.












 §.§ Preliminaries
 
We use the following notation.
For every t∈(0,) and
 f:[0,)→^d, let f_t≐sup_[0,t]f. In case that d=1, we often use |f|_t. 
For any two metric spaces _1, _2 denote by (_1:_2) the space of continuous functions mapping _1 to _2. When _2 =, we use the notation (_1). For a Polish space , the space ([0,T]: ) will be equipped with the
 uniform topology. We will denote by ([0,T]: ) the space of functions mapping [0,T] to  that are right continuous and have left limits (RCLL) defined on [0,T]. This space is equipped with the usual Skorohod topology.






Denote by () the space of probability measures on . We endow () with the topology of weak convergence of measures. Convergence in distribution of S valued random variable X_n to X will be denoted as
X_n ⇒ X. For T, L ∈ (0,∞), the space (([0,T]:[0,L])) will be denoted as _T,L.
The Wasserstein distance of order 1 on (), where  is a compact metric space,
 is defined as

    W_1(η',η)=inf{[∫_d(x,y)dπ(x,y)] : π∈(×)   with marginals η' and η},

 where η,η'∈().



For ϕ∈^1,2([0,T]× [0,L]), D_tϕ, Dϕ, D^2ϕ will denote the time derivative and the first two space derivatives of ϕ, respectively.
For x ∈, δ_x ∈() denotes the
Dirac measure at x.




Throughout the paper we will make extensive use of the Skorohod map, which for the particular setting of interest here is recalled below.  Fix T,L>0.


Given ψ∈([0,T]:) such that ψ(0)∈ [0,L], we say the triplet of functions (,ζ_1,ζ_2)∈([0,T]: ^3) solve the Skorohod problem for ψ if the following properties are satisfied:

(i) For every t∈[0,T],  (t)=ψ(t)+ζ_1(t)-ζ_2(t)∈ [0,L].

(ii) ζ_i are nonnegative and nondecreasing, ζ_1(0)=ζ_2(0)=0, and

    ∫_[0,T]1_(0,L]((s))dζ_1(s)=∫_[0,T]1_[0,L)((s))dζ_2(s)=0.

We denote by Γ(ψ)=(Γ_1,Γ_2,Γ_3)(ψ)≐ (,ζ_1,ζ_2) and refer to Γ as the Skorohod map.

It is known that there is a unique solution to the Skorohod problem for every ψ∈([0,T]: ) and so the Skorohod map in Definition <ref> is well defined.
The Skorohod map has the following Lipschitz property (see <cit.>).

There exists c_S ∈ (0,∞) such that
for all ,∈([0, T]:) with (0), (0) ∈ [0,L],

    ∑_i=1^3_i() - _i()_T≤ c_S-_T.












§ THE MFG AND RELATED RESULTS
 
We now provide a precise description of the MFG that was studied in <cit.> and state some relevant results from there. 


 §.§ Description of the MFG

Fix L,T>0. Here T denotes the terminal time of our finite time horizon and [0,L] will be the state space of the controlled process X. Also, let U be a compact subset of  representing the control space. Let (Ω,,{_t},) be a filtered probability space that supports a one dimensional standard _t-Brownian motion B. We will refer to the collection (Ω,,{_t},, B)
as a system and denote it by Ξ.
Given (x,t,ν) ∈ [0,L]× [0,T]×_T,L, we denote by (Ξ, t,x,ν) the collection of
all pairs (α, Z) where α = {α(s)}_0≤ s ≤ T-t is a U-valued
_s-progressively measurable process, Z = {Z(s)}_0≤ s ≤ T-t is a [0,L]×_+×_+ valued _s-adapted continuous process such that, Z=(X,Y,R) and

    Z(s) = (X,Y,R)(s)=Γ(x+∫_0^·b̅(u)du+σ B(·))(s),   s∈[0,T-t],

where 

    b̅(u)≐ b(t+u,ν(t+u),X(u), (u)),  u ∈ [0, T-t],

b:[0,T]×([0,L])×[0,L]× U, ν(s) is the marginal of ν at time instant s and σ is a (strictly) positive constant.

Given ν∈_T,L, t ∈ [0,T], x ∈ [0,L], and a system
Ξ as above, let (α,Z) ∈(Ξ, t,x,ν). The cost functon is given by,
 
    J_ν(t,x,,Z)   ≐[∫_0^T-t f(s+t,ν(s+t),X(s),(s) )ds + g(ν(T),X(T-t))
         +∫_0^T-ty(s+t,ν(s+t))dY(s)+∫_0^T-tr(s+t,ν(s+t))dR(s)],

and the value function is

    V_ν(t,x)=inf_Ξinf_(,Z)∈(Ξ,t,x,ν)J_ν(t,x,,Z).

Conditions on f,g,y,r will be specified below. We now introduce the notion of a solution to the MFG associated with
(<ref>)–(<ref>).

	
A solution to the MFG, 	associated with
	(<ref>)–(<ref>), with initial condition x ∈ [0,L] is defined to be a ν∈_T,L such that
there exist a system Ξ and an (α,Z) ∈(Ξ,0,x,ν) such that Z=(X,Y,R)
satisfies ∘ X^-1 = ν and

    V_ν(0,x) = J_ν(0, x, α,Z).

If there exists a unique such ν, we refer to V_ν(0,x) as the value of the MFG with initial condition
x.




 §.§ Background results


  
The following conditions were used in <cit.> in order to characterize  the value function V_ν and the optimal control.



  (a) There exists c_L ∈ (0,∞) such that for every (t,η,x,),(t',η',x',')∈[0,T]×([0,L])×[0,L]× U,

    |f(t,η,x,)-f(t',η',x',')| + |g(η,x) - g(η',x')|+|b(t,η,x,)-b(t',η',x',')|
          +|y(t,η)-y(t',η')|+|r(t,η)-r(t',η')|
           ≤ c_L (|t-t'|+W_1(η,η')+|x-x'|+|-'|).


  (b)
For every (t,η, x,p) ∈ [0,T]×([0,L])× [0,L]×, there is a unique
(t,η,x,p) ∈ U such that

    (t,η,x,p)=u∈ Umin h(t,η,x,u,p) ,

where 

    h(t,η,x,u,p)=f(t,η,x,u)+b(t,η,x,u)p.




As argued in <cit.>, Berge's maximum theorem (see <cit.>) together with part (b) of the above assumption
implies that 
 is continuous.
Also note that (<ref>) implies that b,f,g,y,r are bounded functions, in particular,

    sup_(t,η,x,u)∈ [0,T]×([0,L])× [0,L]× U|b(t,η,x,u)| ≐ c_B < ∞.

For further discussion about the assumption, see <cit.>.







The Hamilton-Jacobi-Bellman equation for the value function V_ν(t,x) is given as follows.

    -D_t ϕ-H(t,ν(t),x,D ϕ)-1/2σ^2 D^2ϕ=0,     (t,x)∈[0,T]×[0,L],

with the boundary conditions

    ϕ(T,x)=g(ν(T),x),  Dϕ(t,0)=-y(t,ν(t)), and  Dϕ(t,L)=r(t,ν(t)),  t∈[0,T],

where H is the Hamiltonian given as

    H(t,η,x,p)=inf_u∈ U h (t,η,x,u,p).





The following class of Hölder continuous ν∈_T,L plays a key role in the analysis.

    _0:={ν∈_T,L:sup_0≤ s<t≤ TW_1(ν(t),ν(s))/(t-s)^1/2<}.


Fix  ν∈_0 and suppose that Assumption <ref> holds. Then V_ν is continuously differentiable w.r.t. t and twice continuously differentiable with respect to (w.r.t.) x. It is the unique solution of the Hamilton-Jacobi-Bellman 
equation (<ref>) with 



the boundary conditions (<ref>).





Furthermore, with α̂ as introduced in Assumption <ref>,



the map (s,x')↦(s,ν(s),x',D V_ν(s,x')) is continuous and the
feedback control γ̂(u, x')≐(u+t,ν(u+t),x', D V_ν(u+t,x')) is an optimal feedback control for (<ref>) for every t ∈ (0,T). Moreover, any optimal control  for (<ref>) satisfies (u,ω)=γ̂(u,X(u,ω)), _T^t⊗ almost surely (a.s.), where λ_T^t denotes the Lebesgue measure on [0,T-t].

Using Proposition <ref>, <cit.> proves the existence of a solution of MFG under Assumption <ref>.
In order to establish  uniqueness of the solution, we need an additional condition. Fix η_0 ∈([0,L]).

For every (t,η,x,u)∈[0,T]×([0,L])×[0,L]× U,

    b(t,η,x,u)   =b(t,η_0,x,u),  
    f(t,η,x,u)=f_0(t,η,x)+f_1(t,x,u),
    
    y(t,η)   =y(t,η_0),  
    r(t,η)=r(t,η_0).

Moreover, for every t∈[0,T] and η,η'∈([0,L]), f_0 and g satisfy the following monotonicity property
 
    ∫_0^L [f_0(t,η,x)-f_0(t,η',x))]d(η-η')(x)   ≥ 0,
    
     ∫_0^L (g(η,x)-g(η',x))d(η-η')(x)   ≥ 0.


Abusing notation, when Assumption <ref> holds, we will write b(t,x,u) =b(t,η_0,x,u), y(t)=y(t,η_0), and r(t)=r(t,η_0).
The following is one of the main results from <cit.>.




Under Assumption <ref>, there exists a solution of the MFG. If in addition Assumption <ref> holds then there is a unique MFG solution.



 §.§ Rate control in queues with strategic servers



The MFG described above arises from the heavy traffic analysis of a large queuing system that consists of many symmetric strategic servers that are weakly interacting. Consider a collection of n critically loaded single server queues. Given past information, each server controls the arrival and service rate associated with its own queue. In addition the rates depend on time, the individual queue length, and the empirical measures of all the queue states. The servers aim to minimize individual costs, that in particular account for the scaled idleness and rejection processes. The cost also depends on the individual queue state, the control action and the state of the overall system given through the empirical measure of states of all queues.  The main goal is to find asymptotic Nash equilibrium in this game as the system approaches criticality (i.e. heavy traffic limit) and the number of queues approach ∞, simultaneously.  It is shown in <cit.> that given a solution of the MFG and an optimal control associated with it of the form in Section <ref>, one can construct an asymptotic (in number of players and in heavy traffic limit) Nash equilibrium for the n-player game  such that the solution of the MFG and its associated value function approximate the empirical distribution of the states of the queues and the value function of each server. The current work provides a numerical approximation for the solution of the MFG that is needed in order for constructing the above n-player asymptotic Nash equilibrium.






































































































§ NUMERICAL SCHEME FOR THE MFG


In this section we will use the Markov chain approximation method (<cit.>) to construct numerical solutions of the MFG. The main result of the paper
Theorem <ref> is given here. The numerical scheme is composed of two main steps. First, in Section <ref>, given a probability measure in _T,L,
 we construct a finite state, discrete time, controlled Markov chain and provide a numerical scheme to construct a measure over ([0,T]:[0,L]). Then, in Section <ref> we show that, under assumptions that include the existence of a unique solution of the MFG, the measure  constructed from the chain converges to the solution of the MFG over a small time interval.
 
Throughout the section we assume that Assumption <ref> is satisfied and that (<ref>) holds. 
Note that we do not assume (<ref>) or the monotonicity condition in Assumption <ref>, however
we will introduce additional assumptions as needed.




























We now introduce the  controlled Markov chain constructed on some probability space
(, , ) that will be used to approximate the solution of the MFG.




 §.§ Approximating controlled Markov chains



Fix a discretization parameter h>0 such that L is an integer multiple of h. Denote the h-grid {-h,0,h,…,L+h} by ^h. This is a discretized version of the state space [0,L]. Since 0 and L
are reflecting barriers for the state process X, we will consider two
 types of transition steps for the approximating chain. The first, which occurs when the chain is away from the boundary,
will be referred to as the rate control step and the second occurs at the end points L+h and -h and is referred to as the reflection step.

Rate control step. For every t∈ [0,T], η∈([0,L]), u∈ U, and x∈^h_0≐^h∖{-h,L+h} let

    q^h(t,η,u; x,x± h)   ≐± hb(t,η,x,u)+σ^2/2σ^2,



Note that

    ∑_y ∈{x± h} q^h(t,η,u; x,y) = 1

and that for 0<h<σ^2/c_B, the transition probabilities are positive. Hereafter, these inequalities on h are in force.
 Also, define

    Δ^h≐h^2/σ^2.

This will be used to define the continuous time interpolation of the controlled Markov chain. Denote the Δ^h-grid {0,Δ^h,2Δ^h,…,T-Δ^h} by 𝕋^h.
Note that
Δ^h→ 0 as h→ 0.
One can verify that
the following local consistency conditions (cf. <cit.>) hold for every x∈^h_0,

    m^h_0(t,η,u,x)   ≐∫_^h(x̃-x)q^h(t,η,u; x,dx̃) = b(t,η,x,u)Δ^h,
    
    (σ^h_0)^2(t,η,u,x)   ≐∫_^h(x̃-x-m^h_0)^2q^h(t,η,u; x,dx̃) =σ^2Δ^h-(b(t,η,u,x)Δ^h)^2.


Reflection step. Such a step occurs only when x∈{-h,L+h}.
For every t∈[0,T], η∈([0,L]), and u∈ U let

    q^h(t,η,u; L+h,L)=q^h(t,η,u; -h,0)=1.








We will now define a controlled Markov chain {X_n^h,ν}_n ∈_0 associated with the parameter h, a measure ν∈_T,L and an initial condition x_0∈ [0,L]. 
We choose to work with a deterministic initial state for simplicity of presentation.  The results continue to hold  when the  initial state is random. In that case, in Construction <ref>,
one needs an additional initialization step at which one takes a random draw x_0 from the initial distribution and then sets x^h_0=⌊ x_0/h⌋ h, as in the construction below.


We will assume that
I(h) ≐ T/Δ^h and L/h are integers. 


	The numerical scheme that we develop will be based on controlled Markov chains associated with the probability kernel q^h. Such controlled Markov chain based schemes are  closely related to explicit finite difference schemes for parabolic PDE. Although not studied in the current work, one can also consider Markov chain approximation schemes that have behavior similar to that of implicit finite difference schemes. One of the important steps in convergence proofs
 of finite difference schemes is the identification of appropriate stability conditions for space-time discretizations. 
For the Markov chain approximation method, the analogue of such stability conditions are the local consistency condition of the form in (<ref>)–(<ref>) which form the heart of our convergence proof.  These local consistency requirements in particular imply for our scheme the space-time scaling of the form in (<ref>).  
	





	
	
	
		

	
  * Define X_0^h,ν = x^h_0=⌊ x_0/h⌋ h, set t^h,ν_0=0 and let α_-1^h,ν be a fixed element of U.
	
  * Having defined for i=0, 1, …, n time instants t_i^h,ν<T and random variables X_i^h,ν, α^h,ν_i-1 with values in ^h and U respectively, let _i^h,ν≐σ{X^h,ν_j, α^h,ν_j-1: j = 0, 1, … i}.
	
  * Choose the control α_n^h,ν for the n-th step that is a U valued _n^h,ν measurable random variable and
	let X^h,ν_n+1 be such that its conditional distribution 
	
	given ^h,ν_n is
	q^h(t_n^h,ν,ν(t_n^h,ν), α_n^h,ν,  X_n^h,ν, ·), where ν(t) denotes the marginal distribution of ν
	at time instant t. Also define
	
    t^h,ν_n+1≐ t_n^h,ν + Δ^h 1_{X_n^h,ν∉{-h, L+h}}

	where the indicator in the above definition will ensure that when we do a continuous time interpolation of the chain, reflection steps `occur instantaneously'. Note that the choice of α^h,ν_n is irrelevant if
	X^h,ν_n∈{-h, L+h}.


If α_n^h,ν = ϑ(t_n^h,ν, X_n^h,ν) for some ϑ: 𝕋^h ×𝕊^h → U 
then the function ϑ is referred to as a feedback control.

Some auxiliary processes. We will now introduce some processes that will be useful in the analysis of the h-th Markov chain. 
Consider the piecewise constant processes

    (X^h,ν(t), α^h,ν(t)) ≐ (X^h,ν_n^h,ν(t), α^h,ν_n^h,ν(t)),   t ∈ [0,T],

where

    n^h,ν(t) ≐
    	 max{n: t^h,ν_n = j Δ^h},   t ∈ [jΔ^h, (j+1)Δ^h),   j = 0, … I(h)-1.

Let F^h,ν(0)=B^h,ν(0)=Y^h,ν(0)=R^h,ν(0)=0 and for every t∈[0,T] let

    F^h,ν(t)   ≐∑_j=0^n^h,ν(t)-1[X^h,ν_j+1-X^h,ν_j|^h,ν_j]1_{X^h,ν_j∉{-h,L+h}},
    
    B^h,ν(t)   ≐1/σ∑_j=0^n^h,ν(t)-1(X^h,ν_j+1-X^h,ν_j-[X^h,ν_j+1-X^h,ν_j|^h,ν_j])1_{X^h,ν_j∉{-h,L+h}},
    
    Y^h,ν(t)   ≐∑_j=0^n^h,ν(t)-1(X^h,ν_j+1-X^h,ν_j)1_{X^h,ν_j=-h} = h∑_j=0^n^h,ν(t)-11_{X^h,ν_j=-h},
    
    R^h,ν(t)   ≐∑_j=0^n^h,ν(t)-1(X^h,ν_j-X^h,ν_j+1)1_{X^h,ν_j=L+h} = h∑_j=0^n^h,ν(t)-11_{X^h,ν_j=L+h}.


One can verify that the following representation holds

    (X^h,ν,Y^h,ν,R^h,ν)(t)=Γ(x^h_0+F^h,ν(·)+σ B^h,ν(·))(t),   t∈[0,T].

Also, from (<ref>) and (<ref>) it follows that, on the set {X^h,ν_n∉{-h,L+h}} 

    [X^h,ν_n+1-X^h,ν_n|^h,ν_n]=b(t^h,ν_n,ν(t^h,ν_n),X^h,ν(t^h,ν_n),^h,ν(t^h,ν_n))Δ^h
    ,
    
       [[X^h,ν_n+1-X^h,ν_n-[X^h,ν_n+1-X^h,ν_n|^h,ν_n] ]^2|^h,ν_n ]=
    σ^2Δ^h+o(Δ^h).


Define ^h,ν_t ≐^h,ν_n^h,ν(t). Then since n^h,ν(t) for each fixed t is a {^h,ν_j} stopping time, we have by optional sampling theorem that B^h,ν(·)
is a {^h,ν_t} martingale. Also, from the above,

    F^h,ν(t)   =∫_0^t b(l^h(s),ν(l^h(s)),X^h,ν(s),^h,ν(s))ds,

where

    l^h(s)≐⌊ s/Δ^h⌋Δ^h,     s∈[0,T].








Cost function for the MDP.
For 
every (t,x)∈𝕋^h×𝕊^h and any admissible control α^h,ν used to construct the h-th controlled Markov chain, define the associated cost

    J^h,ν(t,x,^h,ν)   ≐[∫_t^T f(l^h(s),ν(l^h(s)),X^h,ν(s),^h,ν(s))ds+g(ν(T),X^h,ν(T)).
    
           .+∫_t^Ty(s)dY^h,ν(s)+∫_t^Tr(s)dR^h,ν(s)| X^h,ν(t)=x].




The value function associated with the above cost is given by,

    V^h_ν(t,x)≐inf_J^h,ν(t,x,),

where the infimum is taken over all admissible controls.


We now provide properties of the value function V^h_ν and the optimal strategy in the h-th  MDP. 
For every (t,x,ν)∈𝕋^h×𝕊^h_0×, define the h-th finite difference of the value function w.r.t. x, as follows

    ^h_xV^h_ν(t,x)   ≐1/2h(V^h_ν(t+Δ^h,x+h)-V^h_ν(t+Δ^h,x-h)),

where 

    V^h_ν(t+Δ^h,L+h)   ≐ r(t+Δ^h)h+V^h_ν(t+Δ^h,L),
    
     V^h_ν(t+Δ^h,-h)   ≐ y(t+Δ^h)h+V^h_ν(t+Δ^h,0).



The optimal control in the h-th  MDP is given in state feedback form as

    ϑ^h,ν(t,x)= (t,ν(t),x,^h_xV^h_ν(t,x)),   (t,x)∈𝕋^h×𝕊^h.

Letting
 
    ^h,ν(t,x)≐ϑ^h,ν(l^h(t),x),   (t,x)∈ [0,T]×𝕊^h,

 there exists a constant c_d(T)∈(0,), such that for every (t,ν)∈𝕋^h× and for every h, one has,
 
    V^h_ν(t,X^h,ν(t))+ σ∑_s∈𝕋^h, s≥ t^h_xV^h_ν(s,X^h,ν(s))(B^h,ν(s+Δ^h)-B^h,ν(s))
         =g(ν(T),X^h,ν(T))+∫_t^Tf(l^h(s),ν(l^h(s)),X^h,ν(s),^h,ν(s,X^h,ν(s)))ds
           +∫_t^Ty(s)dY^h,ν(s)+∫_t^Tr(s)dR^h,ν(s),

where (X^h,ν,B^h,ν,Y^h,ν,R^h,ν) are as in (<ref>)-(<ref>) with {α^h,ν_n} replaced with the optimal feedback control
ϑ^h,ν and for all (t,x)∈𝕋^h×𝕊^h

    |^h_xV^h_ν(t,X^h,ν(t))|≤ c_d(T).



We note that Lemma <ref> gives, in an explicit form, the finite difference scheme associated with the dynamic programming equation for the cost function (<ref>). Indeed, recall that the function  is given in (<ref>), that the gradient ^h_xV^h_ν(t,x) from (<ref>) is calculated based on the values of the value function at time t+Δ^h, and that the integrals in (<ref>) can be written as finite sums over s∈𝕋^h∩[t,T]. Now, by taking (conditional) expected values in (<ref>), the finite difference scheme follows from a backwards induction.  
The proof of the lemma is given in Section <ref>.

The induced measure Φ^h(ν). Recall from (<ref>) that for j = 0, 1, … I(h), n^h,ν(j Δ^h) = max{i: t_i^h,ν = j Δ^h}. Let {X̂^h,ν(t)}_t∈ [0,T] be the continuous stochastic process which is linear on [jΔ^h, (j+1)Δ^h]
	and equals X^h,ν_n^h,ν(j Δ^h) at t = jΔ^h, for j=0, … I(h)-1, where {X^h,ν_n} is the controlled Markov chain constructed using the
	optimal feedback control ^h,ν.
	
	Let
	
    Φ^h(ν)≐∘(X̂^h,ν)^-1.

	Next, we show that,  under suitable conditions,  Φ^h is a contraction up to an (h^2) term, over a small time interval.  This `almost-contraction' property lies at the heart of our main result, Theorem <ref>.

An almost-contraction property.
Recall that we assume  that Assumption <ref> is satisfied.
In addition, we will make the following assumption on a Lipschitz property of the  function  from (<ref>).
 
There exists c_∈ (0,∞) such that for every t∈[0,T], η,η'∈([0,L]), x,x'∈[0,L], and p,p'∈,

    |(t,η,x,p)-(t,η',x',p')|≤ c_(W_1(η,η')+|x-x'|+|p-p'|).



The following lemma gives a sufficient condition for Assumption <ref> to hold.




Suppose that the drift and the cost functions satisfy the following properties.


  (a) For every (t,η,x,)∈[0,T]×([0,L])×[0,L]× U,

    b(t,η,x,)=b_1(t,η,x)+b_2(t).


  (b) There exists c_m∈ (0,∞) such that for every t∈[0,T], η,η'∈([0,L]), x,x'∈[0,L], and ,'∈ U,

    f(t,η,x,')-f(t,η,x,)-('-)f_(t,η,x,)≥ c_m('-)^2.


  (c)
The map α↦ f(t,η,x,,p) is continuously differentiable for every (t,η,x,p) ∈ [0,T]×([0,L])×[0,L]× and
there exists c_l ∈ (0,∞) such that for every t∈[0,T], η,η'∈([0,L]), x,x'∈[0,L], ∈ U, and p,p'∈,

    |f_(t,η,x,,p)-f_(t,η',x',,p')|≤ c_l(W(η,η')+|x-x'|+|p-p'|).


Then Assumption <ref> is satisfied.

The proof of the lemma is deferred to Section <ref>.


The conditions in the above lemma are not new to the literature of MFG. For example, Assumptions (A.1), (A.2), and (A.3) in <cit.> are stronger. Also, parts (b) and (c) above that concern the running cost are imposed by <cit.>, which studies a rate control problem (part (a) is irrelevant for that model).
A basic example that satisfies parts (a)–(c) in the Lemma, in addition to Assumptions <ref> and <ref>, is the following

    b(t,η,x,)   =b_1(t,x)+b_2(t),
    
    f(t,η,x,)   =a_1(t,x)+a_2(t,x)k()+a_3(t)(c_1+a_4(x))∫_0^La_4(y)dη(y),
    
    g(η,x)   =(c_2+a_5(x))∫_0^La_5(y)dη(y),
    
    y(t,η)   =a_6(t),     r(t,η)=a_7(t),

where b_1,a_1,a_2:[0,T]×[0,L]→,  b_2,a_3,a_6, a_7:[0,T]→, a_4,a_5:[0,L]→ are Lipschitz functions,  c_1,c_2∈,  k: U → is a C^2-strictly convex function (e.g.
k(α) = (α-α_0)^2 for some α_0 ∈),
and
a_2(t,x)≥ c_m>0 for all (t,x)∈ [0,T]× [0,L].

We note that, although Assumption <ref> is not explicitly imposed in the current work, we will assume later in the section that the MFG has a unique solution (see Assumption <ref>) which from Proposition <ref> holds under Assumptions
<ref> and <ref>. For this reason we presented an example that satisfies all three assumptions (i.e., Assumptions <ref>, <ref> and <ref>). From a modeling perspective, by choosing positive and nondecreasing a_3 and a_4 and a positive a_3, the system planner penalizes all servers collectively for congestion when the empirical measure has high
a_3 and a_4-moments and in addition it penalizes  individual servers for long queues. Also, when a_7>0, rejections of jobs by an individual server are disincentivized and when a_6<0, idleness is being rewarded. Finally a convex nondecreasing k assigns costs for increasing the rates.  








The next result plays an important role in the proof of Theorem <ref>. Recall Assumption <ref>
and (<ref>) are in force.

Suppose that Assumption  <ref> is satisfied. Then there exist T̂>0, ĥ>0 and q∈(0,1), such that for every T≤T̂, ν,ν'∈ and h∈(0,ĥ∧T̂),

    W^2_1(Φ^h(ν),Φ^h(ν')) ≤ q (h^2+ W^2_1(ν,ν')).


The proof of the proposition is given in Section <ref>.





 §.§  Approximating the solution of the MFG

We now provide the numerical scheme that approximates the solution of the MFG.
 
	
	Let T̂, ĥ be as in Proposition <ref>.
Fix T<T̂ and (x^h_0,ν^1,h)∈𝕊^h×× (0,ĥ∧T̂). Let {X^h,ν^1_n} be the h-th Markov chain from Construction <ref> associated with the optimal control ^h,ν^1. Having defined for m∈ the process {X^h,ν^m_n}
, set ν^m+1≐Φ^h(ν^m) and let {X^h,ν^m+1_n} be the h-th Markov chain from Construction <ref> associated with the optimal control ^h,ν^m+1.



With q∈ (0,1) as in Proposition <ref> we get that for every h as in Construction <ref> and every k∈,

    W^2_1(Φ^h(ν^k),ν^k)=W^2_1(Φ^h(ν^k),Φ^h(ν^k-1))≤ q(h^2+W^2_1(ν^k,ν^k-1)).

By iterating this bound we obtain

    W^2_1(Φ^h(ν^k),ν^k)≤q/1-qh^2+q^k-1W^2_1(ν^2,ν^1)).

Set

    k_h≐min{k∈ : W^2_1(Φ^h(ν^k),ν^k)≤2q/1-qh^2}

and

    ν_h≐ν^k_h.

We note that k_h depends also on ν^1, however, it plays no role in the sequel and is therefore omitted from the notation. Processes (X^h,ν_h,Y^h,ν_h,R^h,ν_h,B^h,ν_h) are defined as in Construction <ref>,
replacing ν with ν_h and α_n^h,ν with α̂_n^h,ν_h.


As an immediate consequence of the definition of ν_h, we get the following proposition, which is key to the proof of the approximation result in Theorem <ref> below.

Suppose that Assumption  <ref> is satisfied. Then with T̂ as in Proposition <ref>,
 for every T≤T̂

    lim_h→ 0 W^2_1(Φ^h(ν_h),ν_h)=0.


For the main result of this section (Theorem <ref>), in addition to Assumptions <ref>, <ref> and the property in (<ref>), we also need the following assumption.

	 
	There is a unique ν̅∈_T,L that solves the MFG with initial condition x.


In order to formalize the main result we introduce the notion of relaxed controls. The reason is that we need to argue the tightness of control sequences in an appropriate space.
For this, we borrow a relaxed control formulation from <cit.>. Consider the relaxation of the stochastic control problem in (<ref>)–(<ref>) where the control space U is replaced by
(U), the drift function b is replaced by the function b_:[0,T]× [0,L]×(U) →
defined as

    b_(t,x, r) ≐∫_U b(t, x, u) r(du),
    ,

and the running cost f is replaced by f_: [0,T]×([0,L])× [0,L]×(U) →, defined as

    f_(t,η, x, r) ≐∫_U f(t, η, x, u) r(du)
    .

Finally, we replace the class of admissible controls (Ξ, t,x, ν̅) by _(Ξ, t,x, ν̅) of pairs (α_, Z) that are similar to pairs 
(α, Z) introduced above (<ref>) except that α_ is (U) valued rather than U valued and in (<ref>) we replace
b̅(u) = b(u, X(u), α(u)) with b_(u, X(u), α_(u)). The corresponding cost function
J_ν̅, is defined by (<ref>) with f replaced by f_. The value function in this relaxed formulation, denoted as V_ν̅,,
is given by (<ref>) with  replaced by _. Define the function h_ by (<ref>), replacing (f,b) with (f_,b_).
Then, from Assumption <ref>(b),

    H(t,η,x,p) = inf_u ∈ U h(t,η, x, u,p) = inf_r ∈(U) h_(t,η, x, r,p).

To see the last equality note that

    inf_r ∈(U) h_(t,η, x, r,p)=inf_r∈(U)∫ h(t,η, x, u,p)dr(u)≥inf_u∈ Uh(t,η, x, u,p)
 and on the other hand 

    inf_u∈ Uh(t,η, x, u,p)= inf_u∈ U∫ h(t,η, x, a,p)dδ_u(a)≥inf_r ∈(U) h_(t,η, x, r,p).


Therefore, V_ν and V_ν, are both solutions of the partial differential equation (<ref>)-(<ref>). In view of the uniqueness result given in Proposition <ref>, V_ν = V_ν,.



	
Recall from Proposition <ref> that 
Assumption <ref> is satisfied if in addition to Assumption <ref>,
Assumption <ref> holds.
Also, from Assumptions <ref> and <ref> it follows from arguments as in the proof of Proposition 3.1 in <cit.> (see also the statement of Proposition <ref> here),
that there is a continuous map 
γ:[0,T] × [0,L] → U such that if there exist a system Ξ and an (α_,Z) ∈_(Ξ,0,x,ν) such that

    Z(t)≡ (X,Y,R)(t)=Γ(x+∫_0^·
    b_(t, ν(t), X(t), α_(t))dt +σ B(·) )(t),  t ∈ [0,T],

 ∘ X^-1 =  ν, and

    V_ν(0,x) = J_ν(0, x, α_,Z).

then α_(t, ω) = δ_γ(t, X(t, ω)), λ_T^0⊗-a.s., where recall that λ_T^0 denotes the Lebesgue measure on [0,T], and ν = ν̅.











We now present the main result of the paper. 
Let (U × [0,T]) be the space of finite measures on U × [0,T] equipped with the topology of weak convergence. 
Define (U× [0,T]) valued random variable m̂^h,ν_h as

    m̂^h,ν_h(du  ds) ≐δ_α̂^h,ν_h(s,X^h,ν_h(s))(du) ds.
  Also, recall Assumption <ref>
and (<ref>) are in force throughout this section.

	
Suppose that T≤T̂ where T̂ is as in Proposition <ref>.
Also suppose that  Assumptions
<ref> and <ref> are satisfied. Recall the processes (X^h,ν_h, Y^h,ν_h, R^h,ν_h, B^h,ν_h) introduced below Construction <ref> and
consider a sequence {h}→ 0.
Then
the sequence 
    (X^h,ν_h, Y^h,ν_h, R^h,ν_h, B^h,ν_h, m̂^h,ν_h),
  converges in distribution to
(X,Y, R, B, m) in ([0,T]:^4)×(U × [0,T]) and

    lim_h→0ν_h= ν

 in _T,L where the
limit processes defined on some probability space (Ω, , ) satisfy the following.

	
  (a) B is a _t ≐σ{ B(s), X(s), Y(s), R(s), m(A×[0,s]): s ≤ t, A ∈(U)} Brownian motion and so
	 Ξ≐ (Ω, , {_t}, , B) is a system.
	
  (b) Disintegrating m(du  ds) = m_s(du) ds, the following relationship holds a.s.
	
    Z(t)≐ (X(t), Y(t), R(t)) = Γ(x + ∫_0^·b_(s, ν(s), X(s), m_s) ds + σ B(·))(t),   t ∈ [0,T].

	
  (c) ∘ (X)^-1 = ν.
	
	
  (d) The pair (m,Z) ∈_(Ξ,0,x,ν) and V_ν(0,x) = J_ν, (0, x, m,Z). In particular, with γ as in Remark <ref>,
	m(du  ds) = δ_γ(s, X(s))(du) ds, and ν = ν̅, the unique solution of the MFG.






















 
Proof.
Since many steps in the proof are quite standard we will only provide details where appropriate. Using properties
(<ref>) and (<ref>) of the controlled transition probability kernel it can be argued (cf. Proof of <cit.>) that {(F^h,ν_h, B^h,ν_h)}_h>0 are tight in ([0,T]:^2). Using this tightness property along with the continuity of the Skorohod map (Lemma <ref>) it now follows that {(X^h,ν_h, Y^h,ν_h, R^h,ν_h, B^h,ν_h)}_h>0 is tight in
([0,T]:^4). In fact, this sequence is -tight.
Next note that, since m̂^h,ν_h(U× [0,T]) =T and U is compact, the sequence {m̂^h,ν_h}_h>0
is tight in (U× [0,T]). Also, recalling the definition of the interpolated processes (right before (<ref>)) it can be checked that

    |X^h,ν_h - X̂^h,ν_h|_T → 0,  h→ 0.

Combining this with the tightness of {X^h,ν_h} and the fact that Φ^h(ν_h) = ∘ (X̂^h,ν_h)^-1 gives the relative compactness of {Φ^h(ν_h)} in _T,L. Suppose now that along a subsequence (relabeled again as {h})

    (X^h,ν_h, Y^h,ν_h, R^h,ν_h, B^h,ν_h, m̂^h,ν_h) ⇒ (X,Y, R, B, m)

	and
	
    Φ^h(ν_h) →ν.

Then by (<ref>) we also have that

    ν_h→ν.

By standard martingale methods it follows that B is a {_t} Brownian motion (see e.g. Proof of <cit.>) proving part (a) of the theorem.

Using (<ref>),
F^h,ν_h converges, along with the above processes, in distribution to

    ∫_U× [0,·] b(s, ν(s), X(s), u) m(du  ds) = ∫_0^· b_(s, ν(s), X(s), m_s) ds,

where m(du ds) = m_s(du) ds. Using the continuity property of the Skorohod map we now get (b).

Also, part (c) is immediate on using (<ref>), recalling that Φ^h(ν_h) is the probability law of X̂^h,ν_h, and by (<ref>).

Clearly (m,Z) ∈_(Ξ, 0,x, ν). We will now argue the first statement in (d), namely

    V_ν(0,x) = J_ν, (0, x, m,Z).

This together with Remark  <ref> will prove the second statement in (d) for the subsequence.

Since the convergent subsequence was arbitrary, we will get the convergence asserted in the statement of the theorem and complete the proof.


Proof of (<ref>). For the rest of the proof we consider the subsequence along which the convergence in (<ref>) and (<ref>) holds.
Using arguments similar to those in Proposition 4.2 in <cit.> it can be checked that

    lim_h→ 0 J^h,ν_h(α̂^h,ν_h) = J_ν, (0, x,m,Z),

where we have suppressed (0,x_0^h) from the notation in J^h,ν_h.
Next, let (β, Z) ∈(Ξ̅, 0, x, ν) for some system Ξ̅= (Ω̅, , {_t}, , B̅).
	We will now show that for every _0 > 0 there is a sequence of controls {β^h,ν_h_i} for the controlled Markov chains defined as in Construction <ref> with the sequence {ν_h}

, such that
	
    lim sup_h→ 0 J^h,ν_h(β^h,ν_h) ≤ J_ν(0, x, β, Z) + _0.

	Note that from the optimality property of α̂^h,ν_h
	
    lim sup_h→ 0 J^h,ν_h(α̂^h,ν_h) ≤lim sup_h→ 0 J^h,ν_h(β^h,ν_h).

	Combining the above inequality with (<ref>) and (<ref>) we now get that
	
    V_ν(0,x) ≤ J_ν, (0, x, m,Z) = lim_h→ 0 J^h,ν_h(α̂^h,ν_h) ≤lim sup_h→ 0 J^h,ν_h(β^h,ν_h) ≤ J_ν(0, x, β, Z) + _0 .

	Since _0>0, the system Ξ̅, and (β, Z) ∈(Ξ̅, 0, x, ν) are arbitrary we have (<ref>), completing the proof of the theorem.
	
	We now prove (<ref>). Using arguments as in the proof of <cit.>
	it can be shown that there is a θ_1:[0,∞) → [0,∞) such that θ_1(κ)→ 0 as κ→ 0
	and for every >0
	there is a system Ξ^≐ (Ω^,^,{_t^},^, B^)
	and (β^, Z^) ∈(Ξ^, 0,x,ν)
	with the following properties
	
		
  * Z^ satisfies the following equation for t ∈ [0,T].
		
    Z^(t) ≡ (X^,Y^,R^)(t)=Γ(x+∫_0^·
    		b(t, ν(s), X^(s), β^(s))ds +σ B^(·) )(t).

		
  * For some δ > 0, β^ is piecewise constant on intervals of the form [lδ,(l+1)δ), l=0,1,…, T/δ. For some finite set U^⊂ U, β^(s) takes values in U^ for every
		s ∈ [0,T].
		
  * For some θ>0, for each u ∈ U^
		
    ^(β^(lδ)=u|B^(s),s≤ lδ, β^(jδ),j<l)   =^(β^(lδ)=u|B^(pθ),pθ≤ lδ, β^(jδ),j<l)
       = F^_u(B^(pθ),pθ≤ lδ, β^(jδ),j<l),

		where for suitable Î,L̂∈, F^_u: ^Î× (U^)^L̂→ [0,1] is a measurable function such that
		F_u(·, 𝐮) is continuous on ^Î for every 𝐮∈ (U^)^L̂.
		
  * Letting m^(du dt) ≐δ_β^(t)(du) dt, m(du dt) ≐δ_β(t)(du) dt, as → 0,
		
    (X^,Y^, R^, B^, m^) ⇒ (X,Y, R, B, m)

		in ([0,T]:^4)×(U× [0,T]).
		
  * J_ν(0,x,β^,Z^) ≤ J_ν(0,x,β,Z) + θ_1().
	
	







We will now use the piecewise constant control β^ to construct a collection of control sequences {β^h,ν_h_i} as stated above (<ref>) and for which

    lim_h→ 0 J^h,ν_h(β^h,ν_h, ν_h) = J_ν(0, x, β^, Z^).

Note that since θ_1()→ 0 as → 0, this will prove (<ref>) and complete the proof of the theorem. The construction is carried out as follows.

	
  * Define X^h,ν_h_0 = x^h,ν_h, t_0^1=0 and let β_-1^h,ν_h be a fixed element of U^.
	
  * Having defined for i=0, 1, …, n time instants t_i^h,ν_h<T and random variables X_i^h,ν_h, β^h,ν_h_i-1 with values in ^h and U^ respectively, let _i^h,ν_h≐σ{X^h,ν_h_j, β^h,ν_h_j-1: j = 0, 1, … i}.
	
  * Choose the control β_n^h,ν_h for the n-th step that is a U^ valued _n^h,ν_h measurable random variable as follows:
	
		
  * If t_n-1^h,ν_h and t_n^h,ν_h both lie in [jδ, (j+1)δ) for some j, set β^h,ν_h_n = β^h,ν_h_n-1.
		
  * If t_n-1^h,ν_h< lδ≤ t_n^h,ν_h < (l+1)δ for some l, choose β^h,ν_h_n according the the conditional distribution
		
    P(β^h,ν_h_n = u| X_i^h,ν_h, β^h,ν_h_i-1, 0 ≤ i ≤ n) = F_u^(B^h,ν_h(pθ),pθ≤ lδ, β^h,ν_h(jδ),j<l),

		where B^h,ν_h is defined as in (<ref>), β^h,ν_h(s) = β^h,ν_h_n^h,ν_h(s), and n^h,ν_h(·) is as in (<ref>).
	
	
  * Let X^h,ν_h_n+1 be such that the conditional distribution of X^h,ν_h_n+1 
	
	given ^h,ν_h_n equals 
	q^h(t_n^h,ν_h, ν_h(t_n^h,ν_h), β_n^h,ν_h, X_n^h,ν_h, ·), where ν_h is as introduced above Theorem <ref>.
	 Also define
	
    t^h,ν_h_n+1≐ t_n^h,ν_h + Δ^h 1_{X_n^h,ν_h∉{-h, L+h}}.


Now define processes X^h,ν_h, F^h,ν_h, Y^h,ν_h, R^h,ν_h as in (<ref>)–(<ref>). Exactly as in the first part of the proof we now have that {(X^h,ν_h, Y^h,ν_h, R^h,ν_h, B^h,ν_h)}_h is -tight in
([0,T]:^4) and the sequence {m^h,ν_h}_k≥ 1, where m^h,ν_h(du ds) = δ_β^h,ν_h(s)(du) ds
is tight in (U× [0,T]). Arguing as before, if along a further subsequence
the convergence (<ref>) holds (with m̂^h,ν_h replaced with m^h,ν_h) then parts (a) and (b) as in the statement of Theorem <ref> are satisfied with ν as in (<ref>).
Using the continuity property of F_u and the fact that the control is piecewise constant with values in a finite set it follows that
(cf. Proof of <cit.>) (B, m) has the same distribution as (B^, m^). By unique solvability of (<ref>), that follows from the
Lipschitz property of b (Assumption <ref>) and the Lipschitz property of the Skorohod map (Lemma <ref>) we now have that
(Z^, m^) has the same law as (Z, m) for every limit point of the chosen further subsequence. Since the chosen further subsequence was arbitrary, this proves the weak convergence of (Z^h,ν_h, m^h,ν_h) to (Z^,m^) along the subsequence fixed above (<ref>) and
 arguing once again as in the proof of Proposition
4.2 in <cit.> we have the convergence of costs
as in (<ref>), completing the proof of the theorem.






§ PROOFS OF RESULTS FROM SECTION <REF>
 
In this section we provide the  proofs of Lemmas <ref>, <ref> and Proposition <ref>.



Proof of Lemma <ref>. 
We start by analyzing the evolution of the value function V^h_ν. Clearly, V^h_ν(T,x) = g(ν(T),x)
for all x ∈𝕊^h_0. Using backwards induction, we get that for any 
(t,x,ν)∈𝕋^h×𝕊^h_0×, one has
from the definition of  (cf. (<ref>)), q^h (see (<ref>)) and (<ref>)-(<ref>),

    V^h_ν(t,x)   =min_u∈ U{Δ^hf(t,ν(t),x,u)
               +q^h(t,ν(t),u;x,x+h)V^h_ν(t+Δ^h,x+h)
    
    
               +q^h(t,ν(t),u;x,x-h)V^h_ν(t+Δ^h,x-h)}
       =Δ^hf(t,ν(t),x,^h,ν(t,x))
         +q^h(t,ν(t),^h,ν(t,x);x,x+h)V^h_ν(t+Δ^h,x+h)
    
         +q^h(t,ν(t),^h,ν(t,x);x,x-h)V^h_ν(t+Δ^h,x-h),

where ^h,ν is as in (<ref>).
The above identity in particular shows that ϑ^h,ν gives an optimal feedback control.

We now use (<ref>) to prove (<ref>).
Define the following h-th finite differences for (t,x)∈𝕋^h×𝕊^h_0,

    ^h_tV^h_ν(t,x)   ≐1/Δ^h(V^h_ν(t+Δ^h,x)-V^h_ν(t,x)),
    
    ^h_xxV^h_ν(t,x)   ≐1/h^2(V^h_ν(t+Δ^h,x+h)-2V^h_ν(t+Δ^h,x)+V^h_ν(t+Δ^h,x-h)).

Simplifying (<ref>) by using (<ref>), we get that
 
    ^h_tV^h_ν(t,x) =
            -H(t,ν(t),x,^h_x(t,x)) -1/2σ^2^h_xxV^h_ν(t,x).




Notice also that X^h,ν(t+Δ^h) X^h,ν(t) if and only if R^h,ν(t+Δ^h)-R^h,ν(t)=Y^h,ν(t+Δ^h)-Y^h,ν(t)=0, in which case, one has

    V^h_ν(t+Δ^h,X^h,ν(t+Δ^h))-V^h_ν(t,X^h,ν(t))
         =^h_tV^h_ν(t,X^h,ν(t))Δ^h+^h_xV^h_ν(t,X^h,ν(t))(X^h,ν(t+Δ^h)-X^h,ν(t))+1/2h^2^h_xxV^h_ν(t,X^h,ν(t)).

This can be easily verified by considering separately the cases X^h,ν(t+Δ^h)=X^h,ν(t)± h.
In case that R^h,ν(t+Δ^h)-R^h,ν(t)=h, that is X^h,ν(t+Δ^h)= X^h,ν(t)=L, (<ref>) implies that

    V^h_ν(t+Δ^h,X^h,ν(t+Δ^h))-V^h_ν(t,X^h,ν(t))
          =V^h,ν(t+Δ^h,L)-V^h,ν(t+Δ^h,L+h)+V^h,ν(t+Δ^h,L+h)-V^h,ν(t,L)
         =-r(t+Δ^h)h+^h_tV^h_ν(t,X^h,ν(t))Δ^h+^h_xV^h_ν(t,X^h,ν(t))h+1/2h^2^h_xxV^h_ν(t,X^h,ν(t))

and finally, in case that Y^h,ν(t+Δ^h)-Y^h,ν(t)=h,

one similarly has

    V^h_ν(t+Δ^h,X^h,ν(t+Δ^h))-V^h_ν(t,X^h,ν(t))
    
         =-y(t+Δ^h)h+^h_tV^h_ν(t,X^h,ν(t))Δ^h+^h_xV^h_ν(t,X^h,ν(t))(-h)+1/2h^2^h_xxV^h_ν(t,X^h,ν(t)).

From the last two equalities, we get that,

    V^h_ν(t+Δ^h,X^h,ν(t+Δ^h))-V^h_ν(t,X^h,ν(t))
         = -H(t,ν(t),X^h,ν(t),^h_xV^h_ν(t,X^h,ν(t)))Δ^h + ^h_xV^h_ν(t,X^h,ν(t))(X^h,ν(t+Δ^h)-X^h,ν(t))
           -(y(t+Δ^h)+
    ^h_xV^h_ν(t,X^h,ν(t))
    )(Y^h,ν(t+Δ^h)-Y^h,ν(t))
           -(r(t+Δ^h)-^h_xV^h_ν(t,X^h,ν(t)))(R^h,ν(t+Δ^h)-R^h,ν(t))
         = -f(t,ν(t),X^h,ν(t),^h,ν(t,X^h,ν(t)))Δ^h + σ^h_xV^h_ν(t,X^h,ν(t))(B^h,ν(t+Δ^h)-B^h,ν(t))
           -y(t+Δ^h)(Y^h,ν(t+Δ^h)-Y^h,ν(t))-r(t+Δ^h)(R^h,ν(t+Δ^h)-R^h,ν(t)).

Summing up the terms over t∈𝕋^h, one gets (<ref>).
We will postpone the proof of (<ref>) to the end of the paper. 


Proof of Lemma <ref>.
Fix t∈[0,T], η,η'∈([0,L]), x,x'∈[0,L], and p,p'∈. Denote =(t,η,x,p) and '=(t,η',x',p'). 
Recall the definition of h from (<ref>). By (<ref>), (<ref>), and the definition of ', we get that

    h(t,η',x',,p')≥ h(t,η',x',',p')≥ h(t,η',x',,p')+('-)h_(t,η',x',,p')+c_m|'-|^2.

From the minimizing property of α we see  that ('-)h_(t,η,x,,p)≥ 0. Subtracting this term from the right side of the above

    c_m|'-|^2   ≤ |'-|·|h_(t,η',x',,p')-h_(t,η,x,,p)|
       ≤ c |'-|(W_2(η,η')+|x-x'|+|p-p'|) ,

where c=c_l+sup_t∈[0,T]|b_2(t)| and the second inequality follows by (<ref>) and (<ref>) . Result follows on dividing both sides by |'-|/c_m.


Proof of Proposition <ref>.  We begin by introducing a coupling between two optimally controlled chains, one associated with ν and the other with ν'.

Coupling. Fix x∈[0,L] and ν,ν'∈. Let {X^ν_n} and {X^ν'_n} be the Markov chains from Construction <ref> associated with the parameter h and the optimal strategies given by (<ref>). Denote by Σ^ν=(n^ν,F^ν,B^ν,X^ν,Y^ν,R^ν) and Σ^ν'=(n^ν',F^ν',B^ν',X^ν',Y^ν',R^ν') the processes that were defined immediately after Construction <ref>, where we suppressed the index h since it is fixed in the rest of the proof.  Also, denote

    b^ν(t)≐ b(l^h(t),ν(l^h(t)),X^ν(t),^ν(t,X^ν(t))),

where recall that l^h(t)=⌊ t/Δ^h⌋Δ^h. Similarly define b^ν'. We now define a coupling of the chains through a time change of an underlying Markov chain {(Z^ν_n,Z^ν'_n)}. The main idea in the construction of the latter Markov chain is to keep track of the proper time. Whenever an `instantaneous jump' occurs for only one of the Z-processes, the other process has a degenerate step, that is, it remains at the same position. Therefore, we use two sequences of times. The first, which we refer as time instants, {t_n} has the same role as in (<ref>). The second is referred as time steps and denoted by {(N^ν_n,N^ν'_n)}. Each of the components counts how many non-degenerate steps the respective Z process has taken so far.


Set Z^ν_0=Z^ν'_0=x_0, t_0=0, and N^ν_0= N^ν'_0=0. Having defined for i=0,1,…,n time instants t^ν_i<T, time steps N^ν_i, N^ν'_i∈, and random variables Z^ν_i,Z^ν'_i with values in 𝕊^h, define them for the (n+1)-th step as follows.

	
  * If Z^ν_n,Z^ν'_n∉{-h,L+h}, then

    (Z^ν_n+1,Z^ν'_n+1)=(Z^ν_n,Z^ν'_n) +
      
        (h,h),    w.p.  (hmin{b^ν(t_n),b^ν'(t_n)}+σ^2)/(2σ^2), 
    
         (h,-h),    w.p.  h( b^ν(t_n)- b^ν'(t_n))^+/(2σ^2), 
    
            (-h,h),    w.p.  h(b^ν(t_n)- b^ν'(t_n))^-/(2σ^2),
    
            (-h,-h),    w.p.  (-hmax{b^ν(t_n),b^ν'(t_n)}+σ^2)/(2σ^2),

where w.p. stands for `with probability', and

    t_n+1=t_n+Δ^h,    
     N^ν_n+1= N^ν_n+1,       and      N^ν'_n+1= N^ν'_n+1,

where x^+=max{0,x} and x^-=max{0,-x}.


  * 
If Z^ν_n∉{-h,L+h} and Z^ν'_n∈{-h,L+h} then

    (Z^ν_n+1,Z^ν'_n+1)=(Z^ν_n,Z^ν'_n)+
     (0, h(-1)^_{Z^ν'_n=L+h})

and

    t_n+1=t_n,    
     N^ν_n+1= N^ν_n,       and      N^ν'_n+1= N^ν'_n+1.


The transition probabilities when Z^ν'_n∉{-h,L+h} and Z^ν_n∈{-h,L+h}  are defined similarly.


  * If Z^ν_n,Z^ν'_n∈{-h,L+h}, then with probability 1,

    (Z^ν_n+1,Z^ν'_n+1)=(Z^ν_n,Z^ν'_n)+(h(-1)^_{Z^ν_n=L+h}, h(-1)^_{Z^ν'_n=L+h})

and

    t_n+1=t_n,    
     N^ν_n+1= N^ν_n+1,       and      N^ν'_n+1= N^ν'_n+1.



  * For every n∈, set

    X^ν_n≐ Z^ν_M^ν_n    and     X^ν'_n≐ Z^ν'_M^ν'_n,

where M^ν_n≐max{m : N^ν_m≤ n} and  M^ν'_n≐max{m : N^ν'_m≤ n}.




With the above construction {X^ν_n} and {X^ν'_n} are controlled Markov chains constructed using the optimal feedback controls ^ν and
^ν' respectively, given on the same probability space. Also relationships (<ref>)–(<ref>) are satisfied by
Σ^ν and Σ^ν'. The above coupling of the two processes gives the joint evolution  (X^ν(t),X^ν'(t))_0≤ t≤ T as follows.  X^ν(0)=X^ν'(0)=x_0 and for every t∈𝕋^h,














    (X^ν(t+Δ^h),X^ν'(t+Δ^h))- (X^ν(t),X^ν'(t))
                 =
      (h_{X^ν(t) L},h_{X^ν'(t) L}),    w.p.  (hmin{b^ν(t),b^ν'(t)}+σ^2)/(2σ^2), 
    (h_{X^ν(t) L},-h_{X^ν'(t) 0}),    w.p.  h(Δ b(t))^+/(2σ^2), 
    (-h_{X^ν(t) 0},h_{X^ν'(t) L}),    w.p.  h(Δ b(t))^-/(2σ^2),
    (-h_{X^ν(t) 0},-h_{X^ν'(t) 0}),    w.p.  (-hmax{b^ν(t),b^ν'(t)}+σ^2)/(2σ^2),

where

    Δ b(t)≐ b^ν(t)-b^ν'(t),     t∈[0,T].




We also define the corresponding `unconstrained' increment as

    (Z^ν(t+Δ^h),Z^ν'(t+Δ^h))-(X^ν(t),X^ν'(t))
               ≐(h,h),    w.p.  (hmin{b^ν(t),b^ν'(t)}+σ^2)/(2σ^2), 
    (h,-h),    w.p.  h(Δ b(t))^+/(2σ^2), 
    (-h,h),    w.p.  h(Δ b(t))^-/(2σ^2),
    (-h,-h),    w.p.  (-hmax{b^ν(t),b^ν'(t)}+σ^2)/(2σ^2).

Bounding W^2_1(Φ^h(ν),Φ^h(ν')).








Denote Δ X(t)≐ X^ν(t)-X^ν'(t). The processes Δ B, Δ R and Δ Y are defined similarly.
Note that

    W^2_1(Φ^h(ν),Φ^h(ν'))≤[|Δ X|^2_T].

 
 
We now estimate [|Δ X|^2_T].
Recall that Δ X(0)=0. From (<ref>) and Lemma <ref>,

    |Δ X|_T≤ |Δ X|_T + |Δ Y|_T + |Δ R|_T
    ≤ c_S(∫_0^T|Δ b(s)|ds+σ|Δ B|_T).

Therefore,

    [|Δ X|^2_T]≤ 2c^2_S[(∫_0^T|Δ b(s)|ds)^2+σ^2 |Δ B|^2_T].

We now estimate the second term on the right side. By using the martingale property
of B^h,ν(t)- B^h,ν'(t) and
Doob's inequality,

    [|Δ B|_T^2]
    ≤ 4[∑_s∈𝕋^h(Δ B(s+Δ^h)-Δ B(s))^2].

From (<ref>) and (<ref>),






    σ |Δ B(s+Δ^h)- Δ B(s)|   ≤|(Δ X+Δ R-Δ Y)(s+Δ^h)- (Δ X+Δ R-Δ Y)(s)|+Δ^h |Δ b(s)|.

If (Z^ν(t+ Δ^h) - X^ν(t))(Z^ν'(t+ Δ^h) - X^ν'(t)) >0, i.e.,
 the unconstrained increments are of the same sign, then

    |(Δ X+Δ R-Δ Y)(s+Δ^h)- (Δ X+Δ R-Δ Y)(s)|
         =
    |(X^ν+R^ν-Y^ν)(s+Δ^h)- (X^ν+R^ν-Y^ν)(s) .
               .-[(X^ν'+R^ν'-Y^ν')(s+Δ^h)- (X^ν'+R^ν'-Y^ν')(s)]|
    =0.

If the signs are different, i.e., (Z^ν(t+ Δ^h) - X^ν(t))(Z^ν'(t+ Δ^h) - X^ν'(t)) <0, then

    |(Δ X+Δ R-Δ Y)(s+Δ^h)- (Δ X+Δ R-Δ Y)(s)|≤ 2h.

Hence,

    σ|(Δ B(s+Δ^h)-Δ B(s))|≤ 2h_Ê_s+|Δ b(s)|Δ^h,

where

    Ê_s≐{ω : (Z^ν(t+ Δ^h) - X^ν(t))(Z^ν'(t+ Δ^h) - X^ν'(t)) <0}.

From (<ref>) we now have that,

    (Ê_s|^h_s)≤ h|Δ b(s)|/(2σ^2),

where {^h_t} is the filtration generated by the process (X^ν(t),X^ν'(t))_0≤ t≤ T.

As a consequence,

    [(σ(Δ B(s+Δ^h)-Δ B(s)))^2|^h_s]
       ≤[(
    2h_Ê_s+|Δ b(s)|Δ^h
    )^2|^h_s]
       ≤2h |Δ b(s)|(Δ^h+ChΔ^h),

where in the above expression, and in the rest of the proof, C refers to a finite positive constant that is independent of h and s, ν, ν' and which  can change from one line to the next. Applying the above bound to (<ref>) and taking h sufficiently small such that Ch≤ 1/2, we get that for sufficiently small h, 

    σ^2[|Δ B|^2_T]   ≤ 12Δ^h∑_s∈𝕋^hh[|Δ b(s)|]=12[∫_0^Th|Δ b(s)|ds].

Combining this with (<ref>) and using the inequality,

    h∫_0^T|Δ b(s)|ds
    
    
       ≤1/2[h^2T^1/2+T^-1/2(∫_0^T|Δ b(s)|ds)^2]
    ≤1/2T^1/2[h^2 + ∫_0^T|Δ b(s)|^2ds]

we get that

    [|Δ X|^2_T]   ≤ C[(∫_0^T|Δ b(s)|ds)^2+h∫_0^T |Δ b(s)|ds]
       ≤ CT^1/2h^2+C(T^1/2+T)[∫_0^T(Δ W(s))^2ds+∫_0^T(Δ X(s))^2ds+∫_0^T(Δ^h_xV^h(s))^2ds
     ],

where for t∈[0,T]

    Δ^h_xV^h(t)≐^h_xV^h_ν(l^h(t),X^ν(t))-^h_xV^h_ν'(l^h(t),X^ν(t)),  Δ W(t) ≐ W_1(ν(l^h(t)),ν'(l^h(t)))

and
the above inequality also uses the Lipschitz property of b (Assumption <ref>), the Lipschitz property of  (Assumption <ref>) and (<ref>).

We now consider the last term on the right side of (<ref>). We will
 show that for some T̂ that does not depend on h, ν, ν' and all T ≤T̂,

    [∫_0^T(Δ^h_xV^h(s))^2ds]≤ C(h^2+ sup_0≤ s≤ T(Δ W(s))^2
    +[|Δ X|_T^2]).




Define,

    Δ V^h(s)   ≐ V^h_ν(s,X^ν(s))-V^h_ν'(s,X^ν'(s)), Δ g(T)≐Δ V^h(T) = g(ν(T),X^ν(T))-g(ν'(T),X^ν'(T)),
    
    
    Δ f(s)   ≐ f(l^h(s),ν(l^h(s)),X^ν(s),^h,ν(s,X^ν(s)))-f(l^h(s),ν'(s),X^ν'(l^h(s)),^h,ν'(s,X^ν'(s))).


From (<ref>) we get that

    Δ V^h(0)+ σ∑_s∈𝕋^hΔ^h_xV^h(s)(B^ν(s+Δ^h)-B^ν(s))
         =- σ∑_s∈𝕋^h^h_xV^h_ν'(s,X^ν'(s))(Δ B(s+Δ^h)-Δ B(s))+Δ g(T)+∫_0^TΔ f(s)ds
            +∫_0^Ty(s)d(Δ Y(s))+∫_0^Tr(s)d(Δ R(s)).

By squaring both sides and taking expectations
, we have  

    (Δ V^h(0))^2+ [(σ∑_s∈𝕋^hΔ^h_xV^h(s)(B^ν(s+Δ^h)-B^ν(s)))^2]
         ≤ 5{[(σ∑_s∈𝕋^h^h_xV^h_ν'(s,X^ν'(s))(Δ B(s+Δ^h)-Δ B(s)))^2]+[(Δ g(T))^2].
    
               +.[(∫_0^TΔ f(s)ds)^2]+  [(∫_0^Ty(s)d(Δ Y(s)))^2]+[(∫_0^Tr(s)d(Δ R(s)))^2]}.

Here we have used the fact that B^h,ν is a {^h,ν_t}-martingale  and therefore has mean 0. Also, we use the elementary inequality that (∑_i=1^5 a_i)^2 ≤ 5 ∑_i=1^5 a_i^2.
We now bound the terms in the inequality above. For any s∈𝕋^h, we get by (<ref>) and (<ref>) that

    [(Δ^h_xV^h(s)(B^ν(s+Δ^h)-B^ν(s))^2|^h_s]=(Δ^h_xV^h(s))^2(Δ^h+C(Δ^h)^2).

Using once more the martingale property of B^ν
 we get that, for  sufficiently small h,

    [(∑_s∈𝕋^hΔ^h_xV^h(s)(B^ν(s+Δ^h)-B^ν(s)))^2]
    
       =∑_s∈𝕋^h[(Δ^h_xV^h(s)(B^ν(s+Δ^h)-B^ν(s)))^2]
       ≥1/2[∫_0^T(Δ^h_xV^h(s))^2ds].


Recall that from (<ref>), |^h_xV^h_ν'(·,X^ν'(·))|_T≤ c_d(M) whenever T ≤ M. Henceforth we will only consider T ∈ [0,M].
Using (<ref>), we get that for sufficiently small h,

    [(σ∑_s∈𝕋^h^h_xV^h_ν'(s,X^ν'(s))(Δ B(s+Δ^h)-Δ B(s)))^2]
         =∑_s∈𝕋^h[(σ^h_xV^h_ν'(s,X^ν'(s))(Δ B(s+Δ^h)-Δ B(s)))^2]
    
         ≤ 3(c_d(M))^2[∫_0^Th|Δ b(s)|ds]
         ≤ CT^1/2[h^2+∫_0^T(Δ W(s))^2ds+∫_0^T(Δ X(s))^2ds+∫_0^T(Δ^h_xV^h(s)))^2ds
     ],

where the last inequality follows by a similar inequality as in (<ref>).

From (<ref>), it is easy to see that

    [(Δ g(T))^2]   ≤ C((Δ W(T))^2+[(Δ X (T))^2]),

and that

    [(∫_0^TΔ f(s)ds)^2]
       ≤
    C{(∫_0^TΔ W(s)ds)^2+ (∫_0^TΔ X(s)ds)^2 +(∫_0^TΔ^h_xV^h(s)ds)^2}
       ≤
    CT[∫_0^T(Δ W(s))^2ds+∫_0^T(Δ X(s))^2ds+∫_0^T(Δ^h_xV^h(s)))^2ds
     ].


Using integration by parts and the boundedness of y














    [   (∫_0^Ty(s)d(Δ Y(s)))^2]
       ≤ C[ |Δ Y|^2_T
    ]
    
    
       ≤ C[|F^ν+σ B^ν-F^ν'-σ B^ν'|^2_T]
    
       ≤ C[(∫_0^T|Δ b(s)|ds)^2+σ^2|Δ B|_T^2]
       ≤ CT^1/2h^2+C(T^1/2+T)[∫_0^T(Δ W(s))^2ds+∫_0^T(Δ X(s))^2ds+∫_0^T(Δ^h_xV^h(s)))^2ds
     ],

where the inequality on the third line is from Lemma <ref>, the fourth line is from (<ref>) and the last line uses (<ref>) and (<ref>).
A similar bound holds for [(∫_0^Tr(s)d(Δ R(s)))^2].



From (<ref>)–(<ref>), we have

    (Δ V^h(0))^2+ 1/2σ^2[∫_0^T(Δ^h_xV^h(s))^2ds]
         ≤ C( (Δ W(T))^2+[(Δ X(T))^2]) 
           +CT^1/2h^2+C(T^1/2+T)[∫_0^T(Δ W(s))^2ds+∫_0^T(Δ X(s))^2ds+∫_0^T(Δ^h_xV^h(s)))^2ds
     ].

Thus we can find a T̂_1 ∈ (0,M) and ĥ_1>0 such that for all T ≤T̂_1 and h ≤ĥ_1 ∧T̂_1,
 (<ref>) is satisfied. Together with (<ref>) we now get that there exist T̂∈ (0,M), ĥ>0 and q∈(0,1) such that for every T≤T̂ and   h∈(0,ĥ∧T̂),

    W^2_1(Φ^h(ν),Φ^h(ν'))
       ≤[|Δ X|^2_T]≤ q(h^2+sup_0≤ s≤ T(Δ W(s))^2)≤ q(h^2+W^2_1(ν,ν')).









We finally prove the last statement in Lemma <ref>, namely the inequality in (<ref>).


Proof of (<ref>). Fix ν∈, x x' in 𝕊^h_0, and t_0∈𝕋^h, which will be regarded as the initial time. As in the proof of Proposition <ref>, one can define a coupling of two processes on the same h-grid, both of which are driven by the same ν. The first one is denoted as (X(s),Y(s),R(s),B(s),(s))_t_0≤ s≤ T, where its components are defined in (<ref>)–(<ref>) with X(t_0)=x and  is the optimal policy for this process. The second process, denoted as (X'(s),Y'(s),R'(s),B'(s),(s))_t_0≤ s≤ T is also given by (<ref>)–(<ref>) using the same control process {α(s)} as for the first one, except that the second one  starts at x', i.e., X'(t_0)=x'.

For every s∈[t_0,T], let Δ X(s)≐ X(s)-X'(s) and

    Δ b(s)≐ b(l^h(s),ν(l^h(s)),X(s),(s))-b(l^h(s),ν(l^h(s)),X'(s),(s)).

Processes Δ Y(s) and Δ R(s) are defined in a similar manner.
By definition Δ X(0)=x-x'. The arguments that lead to (<ref>) can also be applied here, but in fact they are simpler here since Δ W=0 and the controls are the same for both processes.
Specifically,

    [ sup_t_0≤ s≤ T(|Δ X(s)| + |Δ R(s)| + |Δ Y(s)|)^2]
         ≤ C(x-x')^2+C_1(T-t_0)^1/2h^2+C((T-t_0)^1/2+(T-t_0))[∫_t_0^T(Δ X(s))^2ds
     ]
    
    
         ≤ C(x-x')^2(1+T^1/2)+C(T^1/2+T)[∫_t_0^T(Δ X(s))^2ds
     ].

 The second inequality is a consequence of the fact that since x x', h≤ |x-x'|.
Therefore,

    [ sup_t_0≤ s≤ T(Δ X(s))^2]
    ≤ C(x-x')^2(1+T^1/2)+C(T^1/2+T)[∫_t_0^T(Δ X(s))^2ds
     ].

By 
Grönwall's inequality,
[ sup_t_0≤ s≤ T(Δ X(s))^2] ≤ c_T|x-x'|^2,







where 
    c_T≐ C(1+T^1/2)exp{C_1(T^3/2+T^2)}.
 Using the above bound  we have that the left side in (<ref>) can be bounded above by

    C(x-x')^2(1+T^1/2)+C(T^1/2+T)T [sup_t_0≤ s≤ T(Δ X(s))^2
     ]
     ≤c̃_T |x-x'|^2,

where c̃_T:=(4CL^2(1+T^1/2)+C(T^1/2+T)Tc_T)^1/2, and so,

    [ sup_t_0≤ s≤ T( |Δ X(s)| + |Δ Y(s)| + |Δ R(s)| )]≤ (c̃_T)^1/2 |x-x'|,

Consequently, using integration by parts as in (<ref>) and Lipschitz property of f and g, we get that


    V^h_ν(t_0,x')-V^h_ν(t_0,x)
         ≤ J^h,ν(t_0,x',)-J^h,ν(t_0,x,)
         ≤[|g(ν(T),X'(T))-g(ν(T),X(T))|+ |∫_t_0^Ty(s)d(Δ Y(s))|+
    |∫_t_0^Tr(s)d(Δ Y(s))|
    
                +∫_t_0^T|f(l^h(s),ν(l^h(s)),X'(s),(s))-f(l^h(s),ν(l^h(s)),X(s),(s))|ds]
         ≤c̅_T|x-x'|,

where c̅_T depends on the parameters c_L,c_S, the bounds on y and r,  and the terminal time T.
By reversing the roles of the processes we get that for every t_0∈𝕋^h, 

    |V^h_ν(t_0,x)-V^h_ν(t_0,x')|≤c̅_T|x-x'|

and the result follows.










 

§ NUMERICAL STUDY
 
In this section we present a numerical example. 

We set the parameters L=1, T=0.4, σ=1, and U={-0.75,0.25}. Also, b(t,η,x,)=2x+7, f(t,η,x,)=(4x-5η̅)^2+^2, g(η,x)=(4x-5η̅)^2, y(t,x)=0, and r(t,x)=15, where η̅ is the mean of η. Assumption <ref> obviously holds and therefore, by Proposition <ref> the MFG admits a unique solution and Assumption  <ref> holds. 
 The initial state of the MFG is taken to be X(0)=0.5 and in the numerical scheme x^h=X^h(0)=⌊ x/h ⌋· h. We choose the initial function ν^h(t)=δ_⌊ x/h ⌋· h, t∈[0,T].

 We implemented the algorithm described in Construction <ref> by computing 15 iterations of the map Φ^h for each of the h's taken from the set {1/5,1/10,1/15,1/20,1/25}. For each h we calculated the value function of the MDP after each iteration V^h(x^h). 
 Since our example depends on ν^h only through its mean, we also calculated the mean of ν^h(·), which we denote by ϖ^h.
Figure <ref> illustrates the convergence of the value functions of the MDP's to the value function of the MFG. The convergence of the means ϖ^h is illustrated in Figure <ref>. Finally in Figure <ref> we present the distribution taken from the last iteration, ν^1/25. Here we provide a bird's eye view of this distribution, where the darker areas represent greater density.

 
















































 

 








 








 As can be seen by Figure <ref>, shortly after time zero, the means tend to increase w.r.t. the time. The reason for this is that there are two opposing “forces” in our example. The reflection cost r=15 “pushes” the process X^h downwards when it is close to the boundary L and then the optimal control   is close to -0.75. When X^h is relatively far away from the boundary, the control =-0.75 is too costly and then the optimal control  is close to 0.25. As we approach the terminal time, the rejection cost has less impact and therefore the distribution has higher expectation. 


Acknowledgement. We are thankful to the anonymous referee and AE for their suggestions, which helped us to improve the presentation of the paper.




abbrv 
 


Erhan Bayraktar

Department of Mathematics

University of Michigan

Ann Arbor, MI 48109, USA

email: erhan@umich.edu

web: www-personal.umich.edu/∼ erhan/


Amarjit Budhiraja 

Department of Statistics and Operations Research

University of North Carolina

Chapel Hill, NC 27599, USA

email: budhiraj@email.unc.edu
web: http://www.unc.edu/∼ budhiraj/


Asaf Cohen

Department of Statistics

University of Haifa

Haifa 31905, Israel

email: shloshim@gmail.com

web: https://sites.google.com/site/asafcohentau/



