Amarjit Budhiraja
email: budhiraj@email.unc.edu

Jiang Chen
email: jiangc@live.unc.edu

Sylvain Rubenthaler
email: rubentha@unice.fr

A Numerical Scheme for Invariant Distributions of Constrained Diffusions

Keywords: Reflected Diffusions, Heavy Traffic Theory, Stochastic Networks, Skorohod Problem, Invariant Measures, Stochastic Algorithms. MSC2000 Subject Classification: Primary: 60J60, Secondary: 60J70, 60K25. OR/MS subject classification: Primary: Probability: diffusion, applications, Secondary: Probability: Markov processes, stochastic model applications

Reflected diffusions in polyhedral domains are commonly used as approximate models for stochastic processing networks in heavy traffic. Stationary distributions of such models give useful information on the steady state performance of the corresponding stochastic networks and thus it is important to develop reliable and efficient algorithms for numerical computation of such distributions. In this work we propose and analyze a Monte-Carlo scheme based on an Euler type discretization of the reflected stochastic differential equation using a single sequence of time discretization steps which decrease to zero as time approaches infinity. Appropriately weighted empirical measures constructed from the simulated discretized reflected diffusion are proposed as approximations for the invariant probability measure of the true diffusion model. Almost sure consistency results are established that in particular show that weighted averages of polynomially growing continuous functionals evaluated on the discretized simulated system converge a.s. to the corresponding integrals with respect to the invariant measure. Proofs rely on constructing suitable Lyapunov functions for tightness and uniform integrability and characterizing almost sure limit points through an extension of Echeverria's criteria for reflected diffusions. Regularity properties of the underlying Skorohod problems play a key role in the proofs. Rates of convergence for suitable families of test functions are also obtained. A key advantage of Monte-Carlo methods is the ease of implementation, particularly for high dimensional problems. A numerical example of a eight dimensional Skorohod problem is presented to illustrate the applicability of the approach.

1. Introduction Reflected diffusion processes in polyhedral domains have been proposed as approximate models for critically loaded stochastic processing networks. Starting with the influential paper of Reiman [START_REF] Reiman | Open queueing networks in heavy traffic[END_REF], there have been many works [START_REF] Peterson | A heavy traffic limit theorem for networks of queues with multiple customer types[END_REF][START_REF] Dai | A multiclass station with Markovian feedback in heavy traffic[END_REF][START_REF] Mandelbaum | State-dependent stochastic networks. Part I: Approximations and applications with continuous diffusion Limits[END_REF][START_REF] Yamada | Diffusion approximation for open state-dependent queueing networks in the heavy traffic situation[END_REF][START_REF] Kushner | Heavy traffic analysis of controlled queueing and communication networks[END_REF][START_REF] Williams | Diffusion approximations for open multiclass queueing networks: Sufficient conditions involving state space collapse[END_REF] that justify approximations via reflected diffusions rigorously by establishing a limit theorem under appropriate heavy traffic assumptions. Many performance measures for stochastic networks are formulated to capture the long term behavior of the system and a key object involved in the computation of such measures is the corresponding steady state distribution. Although classical heavy traffic limit theorems only justify approximations of the network behavior through the associated diffusion limit over any fixed finite time horizon, there are now several results [START_REF] Gamarnik | Validity of heavy traffic steady-state approximations in open queueing networks[END_REF][START_REF] Budhiraja | Stationary distribution convergence for generalized Jackson networks in heavy traffic[END_REF][START_REF] Budhiraja | Stability of constrained Markov modulated diffusions[END_REF] that prove, for certain generalized Jackson network models, the convergence of steady state distributions of stochastic networks to those of the associated limit diffusions. Such limit theorems then lead to the important question: How does one compute the stationary distributions of reflected diffusions? Indeed, one of the main motivations for introducing diffusion approximations in the study of stochastic processing systems is the expectation that diffusion models are easier to analyze than their stochastic network counterparts. Classical results of Harrison and Williams [START_REF] Harrison | Multidimensional reflected Brownian motions having exponential stationary distributions[END_REF] show that under certain geometric conditions on the underlying problem data, stationary densities of reflected Brownian motions have explicit product form expressions. However, once one moves away from this special family of models there are no explicit formulas and thus one needs to use numerical procedures.

The objective of the current work is to propose and study the performance of one such numerical proce-dure for computing stationary distributions of reflected diffusions in polyhedral domains. For diffusions in R m there are two basic approaches for computation of invariant distributions: PDE methods and Monte-Carlo methods. PDE approaches are based on the well known basic property that invariant densities of diffusions can be characterized as solutions of certain stationary Fokker-Planck equations. For reflected Brownian motions in polyhedral domains the papers [START_REF] Dai | Characterization of the stationary distribution for a semi-martingale reflecting Brownian motion in a convex polyhedron[END_REF][START_REF] Kurtz | A control formulation for constrained Markov processes[END_REF][START_REF] Dai | Reflected Brownian motion in an orthant: Numerical methods for steady-state analysis[END_REF] develop similar characterization results. The characterization in this case is formulated for the invariant density together with certain boundary densities and is given in terms of the second order differential operator describing the underlying unconstrained dynamics and a collection of first order operators corresponding to the boundary reflections. Using this characterization as a starting point Dai and Harrison [START_REF] Dai | Reflected Brownian motion in an orthant: Numerical methods for steady-state analysis[END_REF] develop an approximation scheme for the stationary density by constructing projections on to certain finite dimensional Hilbert spaces that are described in terms of the above collection of differential operators. Although PDE methods such as above are quite efficient for settings where the state dimension m is small, one finds that Monte-Carlo methods, based on the use of the ergodic theorem, have advantages in higher dimensions. With this in mind, we will propose and study here a Monte-Carlo method for the computation of stationary distributions. Approximations of invariant distributions of diffusions in R m using simulation of paths have been studied in several works [START_REF] Basak | Weak convergence of recursions[END_REF][START_REF] Pelletier | Weak convergence rates for stochastic approximation with application to multiple targets and simulated annealing[END_REF][START_REF] Talay | Stochastic Hamiltonian systems: Exponential convergence to the invariant measure, and discretization by the implicit Euler scheme, Markov Process[END_REF][START_REF] Talay | Second order discretization schemes of stochastic differential systems for the computation of the invariant law[END_REF][START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion[END_REF]. One of the key difficulties in using simulation methods is that paths of diffusions cannot be simulated exactly and so one has to contend with two sources of errors: Discretization of the SDE and finite time empirical average approximation for the steady state behavior. In particular, the long term behavior of the discretized SDE could, in general, be quite different from that of the original system and thus a performance analysis of such Monte-Carlo schemes requires a careful understanding of the stability properties of the underlying systems.

The Monte-Carlo approach studied in the current work is inspired by the papers [START_REF] Basak | Weak convergence of recursions[END_REF], [START_REF] Pelletier | Weak convergence rates for stochastic approximation with application to multiple targets and simulated annealing[END_REF], [START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion[END_REF] which have analyzed the properties of weighted empirical measures constructed from a Euler scheme, based on a single sequence of time discretization steps decreasing to zero, for diffusions in R m . For multi-dimensional diffusions with reflection one first needs to describe a suitable analog of an 'Euler discretization step'. In order to do so, we begin with a precise description of the stochastic dynamical system of interest.

Let G ⊂ R m be the convex polyhedral cone in R m with the vertex at origin given as the intersection of half spaces G i , i = 1, . . . , N . Let n i be the unit vector associated with G i via the relation

G i = {x ∈ R m : x, n i ≥ 0}.
Denote the boundary of a set S ⊂ R m by ∂S. We will denote the set {x ∈ ∂G : x, n i = 0} by F i . For x ∈ ∂G, define the set, n(x), of unit inward normals to G at x by n(x) . = {r : |r| = 1, r, xy ≤ 0, ∀y ∈ G}.

With each face F i we associate a unit vector d i such that d i , n i > 0. This vector defines the direction of constraint associated with the face

F i . For x ∈ ∂G define d(x) . =    d ∈ R m : d = i∈In(x) α i d i ; α i ≥ 0; |d| = 1    , where In(x) . = {i ∈ {1, 2, • • • , N } : x, n i = 0}
. Roughly speaking, the set d(x) represents the set of permissible directions of constraint available at a point x ∈ ∂G. In a typical stochastic network setting this set valued function is determined from the routing structure of the network and governs the precise constraining mechanism that is used. This mechanism specifies how a RCLL trajectory ψ with values in R m is constrained to form a new trajectory with values in G, through the associated Skorohod problem, which is defined as follows.

Let D([0, ∞) : R m ) denote the set of functions mapping [0, ∞) to R m that are right continuous and have left limits. We endow D([0, ∞) : R m ) with the usual Skorokhod topology. Let

D G ([0, ∞) : R m ) . = {ψ ∈ D([0, ∞) : R m ) : ψ(0) ∈ G}.
(i) φ(t) = ψ(t) + η(t); (ii) φ(t) ∈ G; (iii) |η|(t) < ∞; (iv) |η|(t) = [0,t] I {φ(s)∈∂G} d|η|(

s);

(v) There exists Borel measurable γ : [0, ∞) → R m such that γ(t) ∈ d(φ(t)), d|η|-almost everywhere and

η(t) = [0,t]
γ(s)d|η|(s).

In the above definition φ represents the constrained version of ψ and η describes the correction applied to ψ in order to produce φ. On the domain D ⊂ D G ([0, ∞) : R m ) on which there is a unique solutions to the Skorokhod problem we define the Skorokhod map (SM) Γ as Γ(ψ) . = φ, if (φ, ψφ) is the unique solution of the Skorokhod problem posed by ψ. We will make the following assumption on the regularity of the Skorokhod map defined by the data

{(d i , n i ); i = 1, 2, • • • , N }. Condition 1.1 The Skorokhod map is well defined on all of D G ([0, ∞) : R m ), that is, D = D G ([0, ∞) : R m
) and the SM is Lipschitz continuous in the following sense. There exists a K < ∞ such that for all

φ 1 , φ 2 ∈ D G ([0, ∞) : R m ), sup 0≤t<∞ |Γ(φ 1 )(t) -Γ(φ 2 )(t)| < K sup 0≤t<∞ |φ 1 (t) -φ 2 (t)|.
We will also make the following assumption on the problem data.

Condition 1.2 For every x ∈ ∂G, there is a n ∈ n(x) such that d, n > 0 for all d ∈ d(x).
The above condition is equivalent to the assumption that the N × N matrix with (i, j) th entry d i , n j is complete-S (see [START_REF] Dupuis | On Lipschitz continuity of the solution mapping to the Skorohod problem, with applications[END_REF][START_REF] Reiman | A boundary property of semimartingale reflecting Brownian motions[END_REF]). When G = R m + and N = m, it is known that Condition 1.1 implies Condition 1.2 (see [START_REF] Taylor | Existence and uniqueness of semimartingale reflecting Brownian motions in an orthant[END_REF]). An important consequence of Condition 1.2 that will be used in our work is the following result from [START_REF] Budhiraja | An ergodic control problem for constrained diffusion processes: Existence of optimal markov control[END_REF] (see also [START_REF] Dai | Existence and uniqueness of semmartingale reflecting Brownian motions in convex polyhedrons[END_REF]).

Lemma 1.1 Suppose that Condition 1.2 holds. Then there exists a

g ∈ C 2 b (R m ) such that ∇g(x), d i ≥ 1 ∀x ∈ F i , i ∈ {1, ..., N }. (1)
We remark here that the function constructed in [START_REF] Budhiraja | An ergodic control problem for constrained diffusion processes: Existence of optimal markov control[END_REF] is defined only on G, however a minor modification of the construction there gives a C 2 extension to all of R m .

We refer the reader to [START_REF] Harrison | Reflected Brownian motion on an orthant[END_REF][START_REF] Dupuis | On Lipschitz continuity of the solution mapping to the Skorohod problem, with applications[END_REF][START_REF] Dupuis | Convex duality and the Skorokhod Problem. I[END_REF] for sufficient conditions under which Condition 1.1 and Condition 1.2 hold. For example, the paper [START_REF] Dupuis | Convex duality and the Skorokhod Problem. I[END_REF] shows that if G = R m + , N = m and the square matrix

D = [d 1 , ..., d m ] is of the form D = M (I -V ),
where M is a diagonal matrix with positive diagonal entries, V is off diagonal and the spectral radius of |V | is less than 1, then both Conditions 1.1 and 1.2 hold. Here |V | represents the matrix with entries (|V ij |), where V ij is the (i, j)-th entry of V .

We now describe the constrained diffusion process that will be studied in this paper. Let (Ω, F , P) be a complete probability space on which is given a filtration {F t } t≥0 satisfying the usual hypotheses. Let (W (t), F t ) be a m-dimensional standard Wiener process on the above probability space. For x ∈ G, denote by X x the unique solution to the following stochastic integral equation,

X x (t) = Γ x + • 0 σ(X x (s))dW (s) + • 0 b(X x (s))ds (t), (2)
where σ : G → R m×m and b : G → R m are maps satisfying the following condition.

Condition 1.3 There exists a 1 ∈ (0, ∞) such that |σ(x) -σ(y)| + |b(x) -b(y)| ≤ a 1 |x -y| ∀x, y ∈ G (3) and |σ(x)| ≤ a 1 , |b(x)| ≤ a 1 , ∀x ∈ G. (4)
Unique solvability of (2) can be shown using the above condition and the regularity assumption on the Skorokhod map. In fact, the classical method of Picard iteration gives the following:

Theorem 1.1 For each x ∈ G there exists a unique pair of continuous {F t } adapted process (X x (t), k(t)) t≥0 and a progressively measurable process (γ(t)) t≥0 such that the following hold:

(i) X x (t) ∈ G, for all t ≥ 0, a.s.

(ii) For all t ≥ 0,

X x (t) = x + t 0 σ(X x (s))dW (s) + t 0 b(X x (s))ds + k(t), (5)
a.s.

(iii) For all T ∈ [0, ∞), |k|(T ) < ∞ a.s.

(iv) Almost surely, for every t ≥ 0,

|k|(t) = t 0 I {X x (s)∈∂G} d|k|(s), k(t) = t 0 γ(s)d|k|(s), and γ(s) ∈ d(X x (s)) a.e, [d|k|].
In this work we are interested in the invariant distributions of the strong Markov process {X x }. One of the basic results due to Harrison and Williams [START_REF] Harrison | Brownian models of open queueing networks with homogeneous customer populations[END_REF] (see also [START_REF] Budhiraja | Simple necessary and sufficient conditions for the stability of constrained processes[END_REF]) on invariant distributions of such Markov processes says that if b and σ are constants and σ is invertible, then X x has a unique invariant probability measure if b ∈ C o (the interior of C), where

C . = - N i=1 α i d i : α i ≥ 0; i ∈ {1, • • • , N } .
This result was extended to a setting with state dependent coefficients in [START_REF] Atar | On positive recurrence of constrained diffusion process[END_REF] as follows. We introduce the following two additional assumptions. For δ ∈ (0, ∞), define C(δ) . = {v ∈ C : dist(v, ∂C) ≥ δ}.

Condition 1.4 There exists a δ ∈ (0, ∞) such that for all x ∈ G, b(x) ∈ C(δ).

Condition 1.5 There exists σ ∈ (0, ∞) such that for all x ∈ G and

α ∈ R m , α ′ (σ(x)σ ′ (x))α ≥ σα ′ α.
The following is the main result of [START_REF] Atar | On positive recurrence of constrained diffusion process[END_REF].

Theorem 1.2 Assume that Conditions 1.1-1.5 hold. Then the strong Markov process {X x (•); x ∈ G} is positive recurrent and has a unique invariant probability measure.

We remark that in [START_REF] Atar | On positive recurrence of constrained diffusion process[END_REF] a somewhat weaker assumption than Condition 1.4 is used, which says that b(x) ∈ C(δ) for all x outside a bounded set. In the current work, for simplicity we will use the stronger form as in Condition 1.4. Conditions 1.1-1.5 will be assumed to hold for the rest of this work and will not be explicitly noted in the statements of various results.

We now summarize some of the notation that will be used in this work. For a Polish space S, P(S) will denote the space of probability measures, and M F (S) the space of finite measures on S endowed with the usual topology of weak convergence. For a closed set

G ⊂ R m , we say f ∈ C 2 b (G), [respectively f ∈ C 2 c (G)] if f is defined on some open set O ⊃ G
and f is a twice continuously differentiable on O with bounded first two derivatives [respectively compact support]. For ν ∈ P(S) and a ν-integrable f : S → R, we write S f dν as f, ν or ν(f ) interchangeably. We will use the symbol "⇒" or " --→ denote convergence in probability and L p respectively. Denote by || • || ∞ the supremum norm. A vector v ∈ R m is said to be nonnegative (and we write v ≥ 0) if it is componentwise nonnegative.

1.1 Numerical Scheme and Main Results Throughout this work, the unique invariant measure for the Markov process {X x } will be denoted by ν. The goal of this work is to develop a convergent numerical procedure for approximating ν. We now describe this procedure.

Let {λ k } k≥1 be a sequence of positive real numbers such that

λ k → 0, as k → ∞ and letting Λ n := n k=1 λ k , Λ n → ∞ as n → ∞. (6)
Note the condition is satisfied if

λ n = 1 n θ with θ ∈ (0, 1]. Define the map S : G × R m → G by the relation S(x, v) = Γ(x + vi)(1), (7)
where i : [0, ∞) → [0, ∞) is the identity map. The map S will be used to construct an Euler discretization of the stochastic dynamical system described by [START_REF] Budhiraja | Ergodic control for constrained diffusions: Characterization using HJB equations[END_REF]. We now introduce the noise sequence that will be used in the Euler discretization of (5).

Let {U k,j ; k ∈ N, j = 1, ..., m} be an array of mutually independent R valued random variables, given on some probability space (Ω, F , P), such that EU k,j = 0 and EU 2 k,j = 1, for all k ∈ N, j = 1, ..., m. We denote the R m valued random variable (U k,1 , ..., U k,m ) ′ by U k . We will make the following assumption on the array {U k,j }.

Condition 1.6 For some α ∈ (0, ∞), Ee λU k,j ≤ e αλ 2 for all k ∈ N, j = 1, ..., m, λ ∈ R.

The above condition is clearly satisfied when U k,j ∼ N (0, 1). Also, using well known concentration inequalities it can be checked that the condition also holds if supp(U k,j ) is uniformly bounded (see Appendix for a proof of the latter statement). Condition 1.6 will be assumed to hold throughout this work.

The Euler scheme is given as follows. Define iteratively, sequences

{X k } k∈N0 , {Y k } k∈N0 of G and R m valued random variables, respectively, as follows. Fix x 0 ∈ G.      X 0 = x 0 , Y k+1 = X k + b(X k )λ k+1 + σ(X k ) λ k+1 U k+1 , X k+1 = S(X k , Y k+1 -X k ) . (8)
Note that {X k } is a sequence of G valued random variables. The last equation of the above display describes a projection for the Euler step that is consistent with the Skorohod problem associated with the problem data.

Define a sequence of P(G) valued random variables as

ν n = 1 Λ n n k=1 λ k δ X k-1 , n ∈ N.
The above random measures define our basic sequence of approximations for the invariant measure ν. In particular, they yield an approximation for any integral of the form G f (x)dν(x) through the corresponding weighted averages:

1 Λ n n k=1 λ k f (X k-1 ). ( 9
)
The following is the first main result of this work.

Theorem 1.3 As n → ∞, ν n converges weakly to ν, almost surely.

The above result ensures that ( 9) gives an almost surely consistent approximation for ν(f ) for any bounded and continuous f . In fact we have a substantially stronger statement as follows:

Theorem 1.4 There exists a ζ ∈ (0, ∞) such that for all continuous f : G → R satisfying lim sup x→∞ e -ζ|x| |f (x)| = 0, we have ν n (f ) → ν(f ), a.s.
The key ingredient in the proof of the above almost sure limit theorems is a certain Lyapunov function that was introduced in [START_REF] Budhiraja | Long time asymptotics for constrained diffusions in polyhedral domains[END_REF] to study geometric ergodicity properties of reflected diffusions. Using this Lyapunov function we establish a.s. bounds on exponential moments of ν n that are uniform in n. These bounds in particular guarantee tightness of {ν n (ω), n ≥ 1}, for a.e. ω. Then the remaining work, for proving the above theorems, lies in the characterization of the limit points of ν n (ω). For this we use an extension of the well known Echeverria criterion for invariant distributions of Markov processes that has been developed in [START_REF] Dai | Characterization of the stationary distribution for a semi-martingale reflecting Brownian motion in a convex polyhedron[END_REF][START_REF] Kurtz | A control formulation for constrained Markov processes[END_REF] (see also [START_REF] Budhiraja | An ergodic control problem for constrained diffusion processes: Existence of optimal markov control[END_REF]). Verification of this criteria (stated as Theorem 2.1 in the current work) for a typical limit point ν 0 of {ν n } requires showing that, ν 0 along with a certain collection {µ i 0 , i = 1, • • • N } of finite measures supported on various parts of the boundary of G satisfy a relation of the form in [START_REF] Pelletier | Weak convergence rates for stochastic approximation with application to multiple targets and simulated annealing[END_REF]. The measures {µ i 0 } are obtained by taking weak limits of certain finite measures constructed from the Euler scheme. Although these pre-limit measures may place positive mass away from the boundary of the domain, we argue using the regularity properties of the Skorohod map (a key ingredient here is Lemma 1.1), that in the limit these finite measures are supported on the correct parts of the boundary.

Under additional assumptions, one can obtain rates of convergence as follows. For α > 0, set

Λ (α) n = λ α 1 + ... + λ α n .
Denote the normal distribution with mean a and variance b 2 by N (a, b 2 ). For φ ∈ C 3 (G) (space of three times continuously differentiable functions on G) and v ∈ R m , let

D 3 φ(x)(v) ⊗3 = i,j,k D 3 i,j,k φ(x)v i v j v k . For f ∈ C 2 c (G), define Af : G → R and D i f : G → R; i = 1, ..., N as Af (x) = b(x) • ∇f (x) + 1 2 σ ′ (x)D 2 f (x)σ(x), x ∈ G, D i f (x) = d i • ∇f (x), x ∈ G,
where ∇ is the gradient operator and D 2 is the m × m Hessian matrix.

∇φ(x), d i = 0, ∀x ∈ F i , ∀i; D 2 φ(x)d i = 0, ∀x ∈ F i , ∀i. ( 10
)
Then the following CLT holds:

Λ n ν n (Aφ) L -→ N 0, G |σ T ∇φ| 2 dν . (b) Slowly decreasing step. Suppose that lim n→∞ (1/ √ Λ n )Λ (3/2) n = λ ∈ (0, +∞], φ ∈ C 3 (G) and D 3 φ is bounded and Lipschitz. Further suppose that      ∇φ(x), d i = 0, ∀x ∈ F i , ∀i; D 2 φ(x)d i = 0, ∀x ∈ F i , ∀i; D 3 •jk φ(x) • d i = 0, ∀x ∈ F i , ∀i, j, k. ( 11
)
Then we have

Λ n ν n (Aφ) L -→ N λ m, G |σ T ∇φ| 2 dν if λ < ∞, (12)
Λ n Λ (3/2) n ν n (Aφ) P -→ m if λ = +∞, ( 13
)
where

m = - 1 6 G R m D 3 φ(x)(σ(x)u) ⊗3 µ(du)ν(dx).
Note that when

λ k = 1 k α , Λ (3/2) n / √ Λ n converges to 0 [resp. ∞, λ ∈ (0, +∞)], if α > 1/2 [resp. α < 1/2, α = 1/2]
. Also note that if φ is a smooth function supported in the interior of G then it automatically satisfies [START_REF] Cottle | The linear complementarity problem[END_REF] and [START_REF] Dai | Reflected Brownian motion in an orthant: Numerical methods for steady-state analysis[END_REF].

Proof of Theorem 1.5 is quite similar to that of Theorem 9 in [START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion[END_REF], the main difference is in the treatment of the reflection terms for which once more we appeal to regularity properties of the Skorohod map and an estimate based on Lemma 1.1 (see proof of [START_REF] Taylor | Existence and uniqueness of semimartingale reflecting Brownian motions in an orthant[END_REF] which is crucially used in proofs of Section 3).

A key step in the implementation of the Euler scheme in [START_REF] Budhiraja | Stationary distribution convergence for generalized Jackson networks in heavy traffic[END_REF] is the evaluation of the one time step Skorohod map S(x, v). In Section 4.1 we describe one possible approach to this evaluation that uses relationships between Skorohod problems and Linear Complementarity problems(LCPs). There are many well developed numerical codes for solving LCPs (for example in MATLAB) and we will describe in Section 4.2 some results from numerical experiments that use a quadratic programming algorithm for LCPs (cf. [START_REF] Cottle | The linear complementarity problem[END_REF]) in implementing the scheme in [START_REF] Budhiraja | Stationary distribution convergence for generalized Jackson networks in heavy traffic[END_REF]. As remarked earlier, one of the advantages of Monte-Carlo methods is the ease of implementation, particularly for high dimensional problems. To illustrate this, in Section 4.2 we present numerical results for a eight dimensional Skorohod problem.

The paper is organized as follows. In Section 2 we prove Theorem 1.3 and 1.4. Theorem 1.3 is proved in two steps. Section 2.1 shows the tightness of the random measures {ν n }, and Section 2.2 characterizes the limit of the measures {ν n } as the invariant measure of the constrained diffusion in (2). Section 2.3 gives the proof of Theorem 1.4. Rate of convergence theorem (Theorem 1.5) is proved in Section 3. Finally we conclude by describing some numerical results in Section 4.

Proofs of Theorems 1.3 and 1.4

The proof of Theorem 1.3 proceeds by showing that for a.e. ω, the sequence of random probability measures {ν n (ω)} n≥1 is tight and then characterizing the limit points of the sequence using a generalization of Echeverria's criteria. Tightness is argued in Section 2.1 while the limit points are characterized in Section 2.2. Finally in Section 2.3, we give the proof of Theorem 1.4.

Tightness

We begin by presenting a Lyapunov function introduced in [START_REF] Atar | On positive recurrence of constrained diffusion process[END_REF] that plays a key role in the stability analysis of constrained diffusion processes of the form studied here (see [START_REF] Atar | On positive recurrence of constrained diffusion process[END_REF][START_REF] Budhiraja | Long time asymptotics for constrained diffusions in polyhedral domains[END_REF][START_REF] Budhiraja | Stationary distribution convergence for generalized Jackson networks in heavy traffic[END_REF][START_REF] Budhiraja | Stability of constrained Markov modulated diffusions[END_REF][START_REF] Budhiraja | An ergodic control problem for constrained diffusion processes: Existence of optimal markov control[END_REF][START_REF] Budhiraja | Ergodic control for constrained diffusions: Characterization using HJB equations[END_REF][START_REF] Budhiraja | Exit time and invariant measure asymptotics for small noise constrained diffusions[END_REF]).

Throughout this work we will fix a δ > 0 as in Condition 1.4.

For x ∈ G, let A(x) be the collection of all absolutely continuous functions z : [0, ∞) → R m defined via

z(t) . = Γ x + • 0 v(s)ds (t), t ∈ [0, ∞), (14)
for some v : [0, ∞) → C(δ) which satisfies

t 0 |v(s)|ds < ∞, for all t ∈ [0, ∞). (15)
Define T : G → [0, ∞) by the relation

T (x) . = sup z∈A(x) inf{t ∈ [0, ∞) : z(t) = 0}, x ∈ G. ( 16
)
The function T has the following properties (see [START_REF] Atar | On positive recurrence of constrained diffusion process[END_REF]).

Lemma 2.1 There exist constants c, C ∈ (0, ∞) such that the following hold:

(i) For all x, y ∈ G, |T (x) -T (y)| ≤ C|x -y|.
(ii) For all x ∈ G, T (x) ≥ c|x|. Thus, in particular, for all M ∈ (0, ∞) the set {x ∈ G :

T (x) ≤ M } is compact. (iii) Fix x ∈ G and let z ∈ A(x).
Then for all t > 0,

T (z(t)) ≤ (T (x) -t) + .
We next present an elementary lemma that will be used in obtaining moment estimates. For k ∈ N, let

F k = σ(U 1 , ..., U k ). Set F 0 = {∅, Ω}.
Lemma 2.2 There exist c 1 , c 2 ∈ (1, ∞) for which the following holds. Let {v i } i∈N be a sequence of R m valued random variables such that v i is F i-1 measurable for all i ≥ 1 and ess sup

ω |v i (ω)| ≡ |v i | ∞ < ∞. Let S n = n i=1 v i • U i , n ∈ N. Then for every r ≥ 0 and n ≥ 1, E max 1≤i≤n e r|Si| ≤ c 1 e c2r 2 n i=1 |vi| 2 ∞ .
Proof. We will only give the proof for the case m = 1. The general case is treated similarly. From Doob's maximal inequalities for submartingales, we have

E max 1≤i≤n e r|Si| ≤ 4Ee r|Sn| ≤ 4 Ee rSn + Ee -rSn
From Condition 1.6, it follows that for every r ∈ R,

E e rSn |F n-1 ≤ e rSn-1 e αr 2 v 2 n ≤ e rSn-1 e αr 2 |vn| 2 ∞
The result now follows by a successive conditioning argument.

2

Define λ : [0, ∞) → [0, ∞) and j : [0, ∞) → N 0 as λ(s) = Λ k ; j(s) = k, if Λ k ≤ s < Λ k+1 , k ∈ N 0 ;
where we define Λ 0 = 0. Define piecewise linear R m valued stochastic process as follows,

Ŵ (t) = i≤j(t) λ i U i + t -λ(t) λ j(t)+1 U j(t)+1 , t ≥ 0.
Let X(t) be the solution of the following integral equation

X(t) = Γ x 0 + • 0 b( X(λ(s)))ds + • 0 σ( X(λ(s)))d Ŵ (s) (t), t ≥ 0.
Clearly, X(λ(t)) = X j(t) for all t ≥ 0.

Fix ρ ∈ (0, 1]. Define ̟ = 1 2(1 + ρ)L , ∆ = 4λ 0 + 16L ln(c 1 ), (17)
where

L = c 2 a 2 1 C 2 K 2 and λ 0 = sup i≥1 λ i . Let V : G → R + be defined as V (x) = e ̟T (x) , x ∈ G. ( 18
)
Lemma 2.3 There exist β ∈ (0, 1) and φ ∈ [0, ∞) such that for each ζ ∈ [0, ρ] and for all t ≥ 0,

E(V 1+ζ (X j(t+∆) )|F j(t) ) ≤ (1 -β)V 1+ζ (X j(t) ) + φ (19) Proof. Fix t ≥ 0 and ζ ∈ [0, ρ]. Define ξ : [λ(t), ∞) → G as ξ(s) = Γ X j(t) + λ(t)+• λ(t) b(X j(u) )du (s -λ(t)), s ≥ λ(t).
Using the Lipschitz property of the Skorokhod map (Condition 1.1), we have sup

λ(t)≤s≤λ(t)+∆+λ0 | X(s) -ξ(s)| ≤K sup λ(t)≤s≤λ(t)+∆+λ0 s λ(t) σ( X(λ(u)))d Ŵ (u) =:K ν(t, ∆). Note that ∆ -λ 0 ≤ λ(t + ∆) -λ(t) ≤ ∆ + λ 0 .
Using this observation along with Lemma 2.1 (i) and (iii),

T ( X(λ(t + ∆))) ≤T (ξ(λ(t + ∆))) + CK ν(t, ∆) ≤(T ( X(λ(t))) -(λ(t + ∆) -λ(t))) + + CK ν(t, ∆) ≤(T ( X(λ(t))) -(∆ -λ 0 )) + + CK ν(t, ∆).
From the above estimate and the definition of V (x), we now have

E(V ( X(λ(t + ∆))) 1+ζ |F j(t) ) V ( X(λ(t))) 1+ζ ≤ E exp(̟(1 + ζ) (T ( X(λ(t))) -(∆ -λ 0 )) + + CK ν(t, ∆)) F j(t) × exp(-̟(1 + ζ))T ( X(λ(t))). ( 20
)
Letting, for q ∈ N 0 , σ q = σ(X q ), we have, for any s ∈ [λ(t), λ(t)

+ ∆ + λ 0 ], s λ(t) σ( X(λ(u)))d Ŵ (u) ≤      j(s) q=j(t) σ q λ q+1 U q+1 , if σ j(s) U j(s) ≥ 0 j(s)-1 q=j(t) σ q λ q+1 U q+1 , if σ j(s) U j(s) < 0 ,
which can be bounded by max

j(t)≤j≤j(s) j q=j(t) σ q λ q+1 U q+1 . Similarly, - s λ(t) σ( X(λ(u)))d Ŵ (u) ≤ max j(t)≤j≤j(s) - j q=j(t) σ q λ q+1 U q+1 .
And therefore ν(t, ∆) = sup

λ(t)≤s≤λ(t)+∆+λ0 s λ(t) σ( X(λ(u)))d Ŵ (u) ≤ max j(t)≤j≤j * t j q=j(t) σ q λ q+1 U q+1 ,
where j * t = j(λ(t) + ∆ + λ 0 ). Using Lemma 2.2, we now have that, with

m 0 = ̟(1 + ζ)CK, E e m0 ν(t,∆) F j(t) ≤ c 1 e c2m 2 0 a 2 1 j * t q=j(t) λq+1 ≤ c 1 e c2m 2 0 a 2 1 (∆+2λ0) . (21)
In the case T ( X(λ(t))) ≥ ∆λ 0 , we have from [START_REF] Harrison | Reflected Brownian motion on an orthant[END_REF] and ( 21) that

E(V ( X(λ(t + ∆))) 1+ζ |F j(t) ) ≤ V ( X(λ(t))) 1+ζ e -̟(1+ζ)(∆-λ0) × c 1 e c2m 2 0 a 2 1 (∆+2λ0) .
Recalling the choice of ̟ and ∆, we now see that

E(V ( X(λ(t + ∆))) 1+ζ |F j(t) ) ≤ (1 -β)V ( X(λ(t))) 1+ζ ,
where β = 1e -3 ln c1 .

In the case T ( X(λ(t))) < ∆λ 0 , we have

E(V ( X(λ(t + ∆))) 1+ζ |F j(t) ) ≤ E e m0 ν(t,∆) |F j(t) ≤ c 1 e c2m 2 0 a 2 1 (∆+2λ0) ≤ c 1 e 1 4L (∆+2λ0) ≡ φ.
Combining the two cases, we have [START_REF] Hall | Martingale limit theory and its application[END_REF]. 2

The following lemma follows from Lemma 2.3 through a recursive argument.

Lemma 2.4 There exists a 2 ∈ (0, ∞) such that

sup t E(V ( X(λ(t))) 1+ρ ) ≤ a 2 . ( 22
)
Proof.

For any t ∈ (∆, ∞), we can find t ′ ∈ (0, ∆] and j ∈ N such that t = t ′ + j∆. By a recursive argument using [START_REF] Hall | Martingale limit theory and its application[END_REF], we then have

E(V ( X(λ(t))) 1+ρ ) ≤ E(V ( X(λ(t ′ ))) 1+ρ ) + φ β . Thus sup t E(V ( X(λ(t))) 1+ρ ) ≤ sup 0≤t<∆ E(V ( X(λ(t))) 1+ρ ) + φ β .
The supremum on the right side is bounded by max j≤j(∆+λ0) E(V (X j ) 1+ρ ), which is finite using Condition 1.6, boundedness of b, σ and the Lipschitz property of Γ. 2

Now we can prove the following lemma.

Lemma 2.5 For a.e. ω, sup n V, ν n (ω) < ∞. Consequently, the sequence {ν n (ω)} n≥1 is tight for a.e. ω.

Proof. Let n 0 be such that Λ n0 > ∆. Then it suffices to consider the supremum in the above display over all n ≥ n 0 . For i ∈ N 0 , define s(i) = inf{j ∈ N 0 : Λ j ≥ i∆}. Then s(⌊Λ n /∆⌋) ≤ n and therefore, for n ≥ n 0 ,

ν n (V ) = 1 Λ n n k=1 λ k V (X k-1 ) = 1 Λ n Λn 0 V ( X(λ(t)))dt ≤ 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 V ( X(λ(t)))dt. Using Lemma 2.3 with ζ = 0, we have 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 V ( X(λ(t)))dt ≤ 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt + 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t)))(1 -β) + φ]dt.
Thus, rearranging terms,

β Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 V ( X(λ(t)))dt ≤ 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt + φ (⌊Λ n /∆⌋ + 1)∆ Λ s(⌊Λn/∆⌋) .
Next note that Λ s(⌊Λn/∆⌋) ≥ ⌊Λ n /∆⌋∆, and, for n ≥ n 0 ,

Λ s(⌊Λn/∆⌋) ≥ Λ s(1) ≥ λ 1 . Thus sup n≥n0 φ (⌊Λ n /∆⌋ + 1)∆ Λ s(⌊Λn/∆⌋) ≤ φ(1 + ∆ λ 1 ) < ∞.
To prove the lemma, it is now enough to show

sup n≥n0 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt < ∞, a.e. ω. (23)
The above expression can be split into two terms:

1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt = 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t))) -V ( X(λ(t + ∆)))]dt + 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t + ∆))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt ≡ T 1 + T 2
Consider the first term:

T 1 = 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 V ( X(λ(t)))dt - 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 V ( X(λ(t + ∆)))dt = 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 V ( X(λ(t)))dt - 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆+∆ ∆ V ( X(λ(t)))dt = 1 Λ s(⌊Λn/∆⌋) ∆ 0 V ( X(λ(t)))dt - 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆+∆ (⌊Λn/∆⌋+1)∆ V ( X(λ(t)))dt ≤ 1 Λ s(⌊Λn/∆⌋) k:Λ k-1 ≤∆ λ k V (X k-1 ). Let Z = k:Λ k-1 ≤∆ λ k V (X k-1
). Then from ( 22), we have EZ ≤ a 2 (∆ + λ 0 ). Combining this with the fact that for n ≥ n 0 , Λ s(⌊Λn/∆⌋) ≥ λ 1 , we have that

sup n≥n0 T 1 (ω) < ∞, a.e. ω. (24)
Next, consider T 2 :

T 2 = 1 Λ s(⌊Λn/∆⌋) (⌊Λn/∆⌋+1)∆ 0 [V ( X(λ(t + ∆))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt = 1 Λ s(⌊Λn/∆⌋) ⌊Λn/∆⌋ i=0 (i+1)∆ i∆ [V ( X(λ(t + ∆))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt.
From Kronecker's Lemma (see page 63 of [START_REF] Durrett | Probability: theory and examples[END_REF]), the last sum is bounded in n a.s. (in fact converges to 0) if the following series is summable a.s.

∞ i=1 1 Λ s(i) (i+1)∆ i∆ [V ( X(λ(t + ∆))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt.
Consider the sum over even and odd terms separately. For even terms, the sum can be written as

∞ k=1 1 Λ s(2k) (2k+1)∆ 2k∆ [V ( X(λ(t + ∆))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt. ( 25
)
Let

ξ k+1 = 1 Λ s(2k) (2k+1)∆ 2k∆ [V ( X(λ(t + ∆))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt
and G k = F j(2k∆) , then we have E(ξ i+1 |G i ) = 0. Also note that ξ i+1 is G i+1 measurable. Thus S n = n i=1 ξ i is a martingale with respect to the filtration {G n }. Consequently, by Chow's Theorem (see Theorem 2.17 of [START_REF] Hall | Martingale limit theory and its application[END_REF]), the series in [START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion[END_REF]

is a.s. summable if ∞ k=1 E(|ξ k | 1+ρ ) < ∞. Now note that E|ξ k | 1+ρ = E   1 Λ s(2k) (2k+1)∆ 2k∆ [V ( X(λ(t + ∆))) -E(V ( X(λ(t + ∆)))|F j(t) )]dt 1+ρ   ≤ 2 1+ρ ∆ 1+ρ Λ 1+ρ s(2k) sup t E(V ( X(t)) 1+ρ ) ≤ 2 1+ρ ∆ 1+ρ a 2 Λ 1+ρ s(2k)
,

where the last inequality follows form Lemma 2.4. Since Λ s(k) ≥ k∆, we have that

∞ k=1 1 Λ 1+ρ s(k) ≤ 1 ∆ 1+ρ ∞ k=1 1 k 1+ρ < ∞.
This proves that the series in ( 25) is summable. The odd terms are treated in a similar manner. Thus we have proved sup

n≥n0 T 2 (ω) < ∞, a.e. ω. (26)
Now ( 23) is an immediate consequence of ( 24) and ( 26), which proves the lemma. 2

2.2 Identification of the limit In this section we will complete the proof of Theorem 1.3 by arguing that for a.e. ω, every weak limit point of ν n (ω) equals ν. For this we will use the following extension of the Echeverria Criteria (see [START_REF] Kurtz | A control formulation for constrained Markov processes[END_REF][START_REF] Weiss | Invariant measures of diffusions in bounded domains[END_REF], see also Theorem 5.7 of [START_REF] Budhiraja | An ergodic control problem for constrained diffusion processes: Existence of optimal markov control[END_REF]).

Theorem 2.1 Let ν 0 ∈ P(G) and µ i 0 ∈ M F (F i ), i = 1, ..., N be such that for all f ∈ C 2 c (G), ν 0 (Af ) + N i=1 µ i 0 (D i f ) = 0. ( 27
)
Then ν 0 = ν.

In order to apply the above theorem to show convergence of ν n to ν, we will consider a sequence of finite measure {µ i n } n∈N ; i = 1, ..., N , which, roughly speaking, correspond to the prelimit versions of the measures {µ i 0 } that appear in the theorem above. We now describe this sequence.

For u ∈ R m , v ∈ G, r ∈ (0, ∞), define, for t ∈ [0, 1], z(u, v, r|t) ≡ z(t) = v + (b(v)r + σ(v) √ ru)t, x(u, v, r|t) ≡ x(t) = Γ(z)(t), y(u, v, r|t) ≡ y(t) = x(t) -z(t).
Then, one can represent the trajectory y as

y(t) = N i=1 d i t 0 α i (s)d|y|(s); t ∈ [0, 1], (28)
where α i (s) ≡ α i (u, v, r|s) ∈ [0, 1] and α i (s) > 0 only if x(s) ∈ F i . Also, let, for t ∈ [0, 1]

Π t (u, v, r) = z(1) + t(x(1) -z(1)), L i (u, v, r) = 1 0 α i (t)d|y|(t), i = 1, ..., N.
Finally for k ∈ N 0 , let

Π t k = Π t (U k+1 , X k , λ k+1 ), L i k = L i (U k+1 , X k , λ k+1 ).
For k ∈ N 0 and i = 1, ..., N , define a M F (R m ) valued random variable m i k by the relation

ψ, m i k = 1 0 E X k [ψ(Π t k )L i k ]dt, ψ ∈ BM + (R m ), (29)
where E X [Z] denotes E[Z|X], and BM + (R m ) is the space of nonnegative bounded measurable functions on R m .

For n ∈ N and i = 1, ..., N , let µ i n be a M F (R m ) valued random variable defined as

µ i n (A) = 1 Λ n n-1 k=0 m i k (A); A ∈ B(R m ).
The following lemma relates the above family of random measures with our approximation scheme. Recall the definition of the filtration {F k } in Section 2.1.

Lemma 2.6 For every f ∈ C 2 b (R m ), there exists a sequence of real random variables {ξ f n } n∈N such that 1 Λ n n-1 k=0 E[f (X k+1 ) -f (X k )|F k ] = N i=1 µ i n (D i f ) + ν n (Af ) + ξ f n , (30)
and sup n ξ f n (ω) < ∞ a.s. Furthermore if f has compact support then ξ f n → 0 a.s. as n → ∞.

Proof. Fix (u, v, r) ∈ R m × G × (0, ∞). Using the notation introduced above, we have from Taylor's theorem,

f (z(1)) -f (v) = ∇f (v), η + 1 2 η ′ D 2 f (v)η + R 2 (v, z(1)) where R 2 (x, y) = f (y) -f (x) -∇f (x), y -x - 1 2 (y -x) T D 2 f (x)(y -x) and η ≡ η(u, v, r) = b(v)r + σ(v) √ ru. Define r 2 (x, y) = 1 2 sup t∈(0,1) ||D 2 f (x + t(y -x)) -D 2 f (x)||, then we have |R 2 (x, y)| ≤ r 2 (x, y)|x -y| 2 . Also f (x(1)) -f (z(1)) = 1 0 df (z(1) + t(x(1) -z(1))) dt dt = 1 0 ∇f (z(1) + t(x(1) -z(1))) • (x(1) -z(1))dt = N i=1 1 0 ∇f (z(1) + t(x(1) -z(1)))dt • d i 1 0 α i (t)d|y|(t).
Fix a k ∈ N and let v = X k , u = U k+1 and r = λ k+1 . Then

E[f (X k+1 ) -f (X k )|F k ] = E[f (x(1)) -f (z(1)) + f (z(1)) -f (v)|F k ].
From the definition of m i k in [START_REF] Reiman | Open queueing networks in heavy traffic[END_REF] and observing that {X k } is a Markov chain (with respect to the filtration {F k }) and U k+1 is independent of F k , it follows that

E[f (x(1)) -f (z(1))|F k ] = N i=1 m i k (D i f ), Using independence of U k+1 from F k once more, E[f (z(1)) -f (v)|F k ] =λ k+1 ∇f (X k ), b(X k ) + 1 2 λ k+1 σ(X k ) ′ D 2 f (X k )σ(X k ) + 1 2 λ 2 k+1 b(X k ) ′ D 2 f (X k )b(X k ) + E[R 2 (X k , X k + η k )|F k ] =λ k+1 Af (X k ) + ξ f (k), where ξ f (k) = 1 2 λ 2 k+1 b(X k ) ′ D 2 f (X k )b(X k ) + E[R 2 (X k , X k + η k )|F k ] and η k = η(U k+1 , X k , λ k+1 ). Thus we have 1 Λ n n-1 k=0 E[f (X k+1 ) -f (X k )|F k ] = 1 Λ n n-1 k=0 N i=1 m i k (D i f ) + λ k+1 Af (X k ) + ξ f (k) = N i=1 µ i n (D i f ) + ν n (Af ) + 1 Λ n n-1 k=0 ξ f (k).
Equality in [START_REF] Reiman | A boundary property of semimartingale reflecting Brownian motions[END_REF] follows on taking

ξ f n = 1 Λn n-1 k=0 ξ f (k). We now show that sup n ξ f n (ω) < ∞ a.s. Write ξ f (k) = 1 2 λ 2 k+1 b(X k ) ′ D 2 f (X k )b(X k ) + E[R 2 (X k , X k + η k )|F k ] ≡ ξ f 1 (k) + ξ f 2 (k).
The term 1

Λn n-1 k=0 ξ f 1 (k) converges to zero because of the boundedness of b and D 2 f . Consider now the contribution from ξ f 2 (k). Let for p ∈ R + ,

h(p) = 1 2 sup x1,x2∈R m |x1-x2|≤p ||D 2 f (x 2 ) -D 2 f (x 1 )||. Then |R 2 (X k , X k + η k )| ≤ h(η k )|η k | 2
and so for some

κ 1 ∈ (0, ∞), |ξ f 2 (k)| ≤ E[h(η k )|η k | 2 |F k ] ≤ ||h|| ∞ κ 1 λ k+1 .
Thus sup n ξ f n (ω) < ∞ a.s. This completes the first part of the lemma. Finally if f in addition has compact support, we have h(p) → 0 as p → 0. Fix ǫ > 0. Since b, σ are bounded, we can find for each θ ∈ (0, ∞),

k θ ∈ N such that for every k ≥ k θ , |h(b(x k )λ k+1 + σ(x k ) λ k+1 U k+1 )|1 |U k+1 |≤θ ≤ ǫ. Also, for some l η ∈ (0, ∞), for all k ∈ N, E[|η k | 2 1 |U k+1 |≥θ |F k ] ≤ l η (λ 3/2 k + λ k E[|U 1 | 2 1 |U1|≥θ ]) a.s., E[|η k | 2 |F k ] ≤ l η λ k a.s. Choose θ 0 ∈ (0, ∞) such that E[|U 1 | 2 1 |U1|≥θ0 ] ≤ ǫ. Then 1 Λ n n-1 k=k θ 0 |ξ f 2 (k)| ≤ ǫl η 1 Λ n n-1 k=k θ 0 λ k + ||h|| ∞ l η ( 1 Λ n n-1 k=k θ 0 λ 3/2 k + ǫ Λ n n-1 k=k θ 0 λ k ).
Thus,

1 Λ n n-1 k=0 |ξ f 2 (k)| ≤ 1 Λ n k θ 0 -1 k=0 ξ f 2 (k) + ǫl η (1 + ||h|| ∞ ) + ||h|| ∞ l η 1 Λ n n-1 k=k θ 0 λ 3/2 k .
Sending n → ∞ and then ǫ → 0, we now see that

1 Λn n-1 k=0 |ξ f 2 (k)| → 0 as n → ∞. The result follows. 2
The following lemma shows that the left side of the expression in [START_REF] Reiman | A boundary property of semimartingale reflecting Brownian motions[END_REF] converges to 0 as n → ∞.

Lemma 2.7 For every f ∈ C 2 b (G), 1 Λ n n-1 k=0 E[f (X k+1 ) -f (X k )|F k ] → 0 a.s., as n → ∞.
Proof. We can split the sum into two terms:

1 Λ n n-1 k=0 E[f (X k+1 ) -f (X k )|F k ] = 1 Λ n n-1 k=0 (E[f (X k+1 )|F k ] -f (X k+1 )) + 1 Λ n n-1 k=0 (f (X k+1 ) -f (X k )) =T 1 + T 2 .
Note that,

|T 2 | = 1 Λ n |f (X n ) -f (X 0 )| → 0,
as n → ∞, since f is bounded and Λ n → ∞. Also, using Kronecker's Lemma,

T 1 = 1 Λ n n-1 k=0 (E[f (X k+1 )|F k ] -f (X k+1 ))
will converge to 0 once the martingale

M f n := n-1 k=1 1 Λ k (E[f (X k+1 )|F k ] -f (X k+1 )) converges a.s. Finally observing that E(f (X k+1 )|F k ) minimizes the L 2 distance from f (X k+1 ) among F k measurable square integrable random variables, E M f ∞ = k≥1 ( 1 Λ k ) 2 E (f (X k+1 ) -E(f (X k+1 )|F k )) 2 ≤ k≥1 ( 1 Λ k ) 2 E (f (X k+1 ) -f (X k )) 2 ≤ ||Df || ∞ k≥1 ( 1 Λ k ) 2 E (X k+1 -X k ) 2 ≤ κ 1 k≥1 λ k+1 Λ 2 k < ∞
for some constant κ 1 , where the last inequality follows from the observation that for a positive sequence λ k , k≥1 λ k+1 /Λ 2 k < ∞. The lemma follows. 2

Next we consider the limit of the first term on the right side of [START_REF] Reiman | A boundary property of semimartingale reflecting Brownian motions[END_REF]. We can regard µ i n to be a finite measure on the one point compactificaion of R m , denoted as Rm . In order to show that {µ i n } is a.s. a precompact sequence in M F ( Rm ), it suffices to show that µ i n (R m ) is an a.s. bounded sequence of R + valued random variables. This is shown in the following lemma.

Lemma 2.8 For

i = 1, ..., N , sup n µ i n (R m ) < ∞, a.s.
Proof. Let g ∈ C 2 b (R m ) be as in Lemma 1.1. Then for fixed (u, v, r) ∈ R m × G × (0, ∞) and with notation as introduced above Lemma 2.6,

g(x(1)) =g(v) + 1 0 [∇g(x(s)) • (b(v)r + σ(v) √ ru)]ds + N i=1 1 0 d i • ∇g(x(s))α i (s)d|y|(s) (31)
Since α i (s) is nonzero only when x(s) ∈ F i , and ∇g(x),

d i ≥ 1, for all x ∈ F i , i ∈ {1, ..., N }, we have N i=1 L i (v, u, r) = N i=1 1 0 α i (s)d|y|(s) ≤ N i=1 1 0 d i • ∇g(x(s))α i (s)d|y|(s) ≤ |g(x(1)) -g(v)| + ||∇g|| ∞ |b(v)r + σ(v) √ ru| ≤ ||∇g|| ∞ |x(1) -v| + ||∇g|| ∞ |b(v)r + σ(v) √ ru| ≤ ||∇g|| ∞ (K + 1)|b(v)r + σ(v) √ ru|, (32)
where the second inequality uses [START_REF] Talay | Second order discretization schemes of stochastic differential systems for the computation of the invariant law[END_REF], and the last inequality uses the Lipschitz property of the Skorokhod map.

Let κ 1 = ||∇g|| ∞ (K + 1)a 1 , then from [START_REF] Talay | Stochastic Hamiltonian systems: Exponential convergence to the invariant measure, and discretization by the implicit Euler scheme, Markov Process[END_REF] we have for i ∈ {1, ..., N },

L i k ≤ κ 1 λ k+1 |U k+1 | + λ k+1 . (33)
Also note that,

sup t∈[0,1] |x k (t) -X k | ≤ K|b(X k )λ k+1 + σ(X k ) λ k+1 U k+1 | ≤ Ka 1 λ k+1 |U k+1 | + Ka 1 λ k+1 , (34)
and for t ∈ [0, 1],

|Π t k -X k | ≤ t|x k (1) -v| + (1 -t)|z k (1) -v| ≤ (K + 1)a 1 λ k+1 + (K + 1)a 1 λ k+1 |U k+1 |, (35)
where 33)- [START_REF] Williams | Diffusion approximations for open multiclass queueing networks: Sufficient conditions involving state space collapse[END_REF] we have that

x k (t) = x(U k+1 , X k , λ k+1 |t), z k (t) = z(U k+1 , X k , λ k+1 |t). Combining (
E X k (|Π t k -x k (s i k )|L i k ) ≤ (2K + 1)a 1 κ 1 mλ k+1 + ϕ(λ k+1 )λ k+1 , (36)
where ϕ : (0, ∞) → (0, ∞) is a bounded function satisfying ϕ(α) → 0 as α → 0.

Next note that L i k is not equal to 0 only if there exists s

∈ [0, 1] such that α k i (s) > 0, i.e., x k (s) ∈ F i , where α k i (t) ≡ α i (U k+1 , X k , λ k+1 |t). And in that case, D i g(Π t k ) ≥ D i g(x k (s)) -||D 2 g|| ∞ |Π t k -x k (s)| ≥ 1 -||D 2 g|| ∞ |Π t k -x k (s)|. Let A i k = {ω : there exists s ∈ [0, 1] such that α k i (s) > 0} and
s i k (ω) = inf{s ∈ [0, 1] : α k i (s) > 0} if ω ∈ A i k , 1 if ω / ∈ A i k . Then, from (36), E X k [D i g(Π t k )L i k 1 A i k ] ≥ E X k [L i k 1 A i k ] -||D 2 g|| ∞ E X k [|Π t k -x k (s i k )|L i k 1 A i k ] ≥ E X k [L i k ] -||D 2 g|| ∞ ((2K + 1)a 1 κ 1 mλ k+1 + ϕ(λ k+1 )λ k+1 ) Thus we have D i g, m i k = 1 0 E X k [D i g(Π t k )L i k ]dt = 1 0 E X k [D i g(Π t k )L i k 1 A i k ]dt ≥ 1, m i k -||D 2 g|| ∞ ((2K + 1)a 1 mκ 1 λ k+1 + ϕ(λ k+1 )λ k+1 ) . Rearranging the terms, we have 1, m i k ≤ D i g, m i k + ||D 2 g|| ∞ ((2K + 1)a 1 mκ 1 λ k+1 + ϕ(λ k+1 )λ k+1
) . Summing over k from 0 to n -1 and i from 1 to N , we obtain

N i=1 1, µ i n ≤ N i=1 D i g, µ i n + N ||D 2 g|| ∞ ((2K + 1)a 1 κ 1 m + |ϕ| ∞ ) . (38)
Using Lemma 2.6

N i=1 µ i n (D i g) = 1 Λ n n-1 k=0 E[g(X k+1 ) -g(X k )|F k ] -ν n (Ag) -ξ g n . Since g ∈ C 2 b (R m
), the second term on the right side is bounded. Also from Lemma 2.7, the first term converges to 0 as n → ∞. Finally from Lemma 2.6, the third term is bounded, a.s.

From this it follows that

sup n N i=1 µ i n (D i g) < ∞ a.s.
Result follows on using this observation in (38). 2

The following lemma will be used to show that for a.e. ω, any limit point of

µ i n (ω) is supported on F i , i = 1, ..., N . Lemma 2.9 Fix i ∈ {1, ..., N }. Let ψ ∈ C 2 c (R m ) be such that ψ(x) ≥ 0 for all x ∈ R m . Suppose that there is a ǫ > 0, such that ψ(x) = 0 if dist(x, F i ) ≤ ǫ. Then ψ(x)µ i n (dx) → 0, a.s. as n → ∞.
Proof. We have

ψ, µ i n = 1 Λ n n-1 k=0 ψ, m i k = 1 Λ n n-1 k=0 1 0 E X k [ψ(Π t k )L i k ]dt = 1 Λ n n-1 k=0 1 0 E X k 1 0 ψ(Π t k )α k i (s)d|y k | s dt ≤ |ψ| ∞ Λ n n-1 k=0 E X k 1 0 1 0 1 {|Π t k -x k (s)|>ǫ} α k i (s)d|y k | s dt, (39)
where, recall that

α k i (s) ≡ α i (U k+1 , X k , λ k+1 |s) and x k (s) = x(U k+1 , X k , λ k+1 |s), y k (s) = y(U k+1 , X k , λ k+1 |s).
The last inequality in the above display follows from noting that α k i (s) > 0 only when x k (s) ∈ F i and if for such a s,

|Π t k -x k (t)| ≤ ǫ,we have by our choice of ψ that ψ(Π t k ) = 0. Next note that {(t, s, ω) : |Π t k -x k (s)| > ǫ} ⊂ {(t, s, ω) : |x k (1) -x k (s)| > ǫ} ∪ {(t, s, ω) : |z k (1) -x k (s)| > ǫ}, where recall that z k (t) = z(U k+1 , X k , λ k+1 |t).
Also, from the Lipschitz property of the Skorokhod map,

|x k (1) -x k (s)| ≤ Ka 1 λ k+1 + Ka 1 λ k+1 |U k+1 |, and
|z k (1) -x k (s)| ≤ |z k (1) -X k | + |X k -x k (s)| ≤ (K + 1)a 1 λ k+1 + (K + 1)a 1 λ k+1 |U k+1 |. Thus {ω : |Π t k -x k (s)| > ǫ for some t, s ∈ [0, 1]} ⊂ {ω : |U k+1 (ω)| ≥ p k }, where p k = ǫ/((K+1)a1)-λ k+1 √ λ k+1
. Using this observation in (39), we have

ψ, µ i n ≤ |ψ| ∞ Λ n n-1 k=0 E X k 1 0 1 0 1 {|U|≥p k } α k i (s)d|y k | s dt ≤ |ψ| ∞ Λ n n-1 k=0 E X k 1 {|U|≥p k } 1 0 α k i (s)d|y k | s ≤ |ψ| ∞ Λ n n-1 k=0 E X k 1 0 α k i (s)d|y k | s 2 P(|U | ≥ p k ).
From [START_REF] Talay | Stochastic Hamiltonian systems: Exponential convergence to the invariant measure, and discretization by the implicit Euler scheme, Markov Process[END_REF] it follows that for some

κ 1 ∈ (0, ∞), sup k E X k 1 0 α k i (s)d|y k | s 2 ≤ κ 1 . Also using Condition 1.6, E|U | j < ∞ for all j ≥ 1. Choose k 0 large enough so that λ k+1 ≤ ǫ 2(K+1)a1 for all k ≥ k 0 . Fix j > 4, then ψ, µ i n ≤ |ψ| ∞ Λ n √ κ 1 k 0 + |ψ| ∞ Λ n n-1 k=k0 √ κ 1 (E|U | j ) 1/2 p -j/2 k
.

The result now follows on observing that for some κ 2 ∈ (0, ∞), p

-j/2 k ≤ κ 2 λ j/4 k+1 for all k ≥ k 0 . 2
We are now ready to complete the proof of Theorem 1.3.

Proof.

[Proof of Theorem 1.3] Fix f ∈ C 2 c (G).
Then such a function can be extended to a function in C 2 c (R m ). We denote this function once more by f . Then from Lemma 2.6,

1 Λ n n-1 k=0 E[f (X k+1 ) -f (X k )|F k ] = N i=1 µ i n (D i f ) + ν n (Af ) + ξ f n . (40)
From Lemmas 2.5, 2.6, 2.7 and 2.8, there exists Ω 0 ∈ F such that P(Ω 0 ) = 1 and for every ω ∈ Ω 0 ,

• {ν n (ω)} n is precompact in P(G),

• {µ i n (ω)} n is precompact in M F ( Rm ), for every i = 1, ..., N , • Left hand side of (40) converges to 0, • ξ f n (ω) converges to 0. Fix a ω ∈ Ω 0 and let ν ∞ (ω), µ i ∞ (ω), i = 1, ..., N , be a subsequential limit of ν n (ω) and µ i n (ω), respectively. Then from (40) and the above observations, we have ( suppressing ω )

ν ∞ (Af ) + N i=1 µ i ∞ (D i f ) = 0.
To complete the proof, in view of Theorem 2.1, it suffices to argue that

R m 1 F c i (x)µ i ∞ (ω)(dx) = 0. ( 41
)
By convergence of µ i n to µ i ∞ , we have for every ψ as in Lemma 2.9, [START_REF] Durrett | Probability: theory and examples[END_REF]. Fix ζ ∈ (0, ̟c). We will prove the theorem with such a choice of ζ. Consider an f as in the statement of the theorem. Then there exists constant κ 1 such that |f (x)| ≤ κ 1 e ζ|x| Without loss of generality, we assume f ≥ 0.

R m ψ(x)µ i ∞ (ω)(dx) = 0. Therefore R m 1 F ǫ,r i (x)µ i ∞ (dx) = 0 ∀ǫ,
From Theorem 1.3, for any L > 0, we have

(f ∧ L)dν n → (f ∧ L)dν a.s.
In order to prove the theorem, it suffices to show that

(f ∧ L)dν n → f dν n , and (f ∧ L)dν → f dν, as L → ∞. First, consider sup n f dν n -(f ∧ L)dν n ≤ sup n 1 f >L f dν n ≤ sup n ν 1/p n (f > L)[ν n (f q )] 1/q ,
where p, q ∈ (1, ∞) are such that p -1 + q -1 = 1 and the last inequality follows from Hölder's inequality. Choose q > 1 such that ζq < ̟c, then from Lemma 2.5 we have

sup n [ν n (f q )] 1/q ≤ κ 1 sup n [ e ζq|x| ν n (dx)] 1/q ≤ κ 1 sup n ν 1/q n (V ) < ∞, a.s. ( 42
)
Using Markov's Inequality, we have

ν 1/p n (f > L) ≤ ν 1/p n (f ) L 1/p ,
which using (42) converges to 0 as L goes to infinity. Combining the above three displays, we have

sup n f dν n -(f ∧ L)dν n → 0, a.s. as L → ∞. (43)
Also, from Fatou's lemma we have, for a.e. ω,

f dν -(f ∧ L)dν = (f -f ∧ L)dν ≤ lim inf n (f -f ∧ L)dν n ≤ sup n (f -f ∧ L)dν n .
Using (43) the last expression converges to 0 as L → ∞. The result follows.

3. Proof of Theorem 1. [START_REF] Budhiraja | Ergodic control for constrained diffusions: Characterization using HJB equations[END_REF] We begin with a few preliminary lemmas.

Lemma 3.1 If φ ∈ C 2 (G), then Λ n ν n (Aφ) = n k=1 λ k Aφ(X k-1 ) = Z (0) n -(N n + 4 i=1 Z (i) n + 4 i=1 Y (i) n )
with

Z (0) n = φ(X n ) -φ(X 0 ), N n = n k=1 λ k ∇φ(X k-1 ), σ(X k-1 )U k , Z (1)
n = 1 2 n k=1 λ 2 k b(X k-1 ) T D 2 φ(X k-1 )b(X k-1 ), Z (2)
n = n k=1 λ 3/2 k b(X k-1 ) T D 2 φ(X k-1 )σ(X k-1 )U k , Z (3) n
= 1 2 n k=1 λ k [(σ(X k-1 )U k ) T D 2 φ(X k-1 )(σ(X k-1 )U k ) -E((σ(X k-1 )U k ) T D 2 φ(X k-1 )(σ(X k-1 )U k )|F k-1 )], Z (4)
n = n k=1 R 2 (X k-1 , X k ), and
Y (1) n = n k=1 ∇φ(X k-1 ), y k-1 , Y (2)
n = 1 2 n k=1 y T k-1 D 2 φ(X k-1 )y k-1 , Y (3)
n = n k=1 λ k b(X k-1 ) T D 2 φ(X k-1 )y k-1 , Y (4)
n = n k=1 λ k y T k-1 D 2 φ(X k-1 )σ(X k-1 )U k , where R 2 (x, y) = φ(y) -φ(x) -∇φ(x), y -x -1 2 (y -x) T D 2 φ(x)(y -x), and
y k = y(U k+1 , X k , λ k+1 |1). Proof. Denote δφ(X k ) = φ(X k ) -φ(X k-1 ) and δX k = X k -X k-1 . We deduce from (8) that δφ(X k ) = ∇φ(X k-1 ), δX k + 1 2 δX T k D 2 φ(X k-1 )δX k + R 2 (X k-1 , X k ) = ∇φ(X k-1 ), y k-1 + λ k Aφ(X k-1 ) + λ k ∇φ(X k-1 ), σ(X k-1 )U k + 1 2 y T k-1 D 2 φ(X k-1 )y k-1 + 1 2 λ 2 k b(X k-1 ) T D 2 φ(X k-1 )b(X k-1 ) + 1 2 λ k [(σ(X k-1 )U k ) T D 2 φ(X k-1 )(σ(X k-1 )U k ) -E((σ(X k-1 )U k ) T D 2 φ(X k-1 )(σ(X k-1 )U k )|F k-1 )] + λ k b(X k-1 ) T D 2 φ(X k-1 )y k-1 + λ 3/2 k b(X k-1 ) T D 2 φ(X k-1 )σ(X k-1 )U k + λ k y T k-1 D 2 φ(X k-1 )σ(X k-1 )U k + R 2 (X k-1 , X k ).
The lemma follows by summing the above equality over k = 1, ..., n and rearranging the terms.

2 Lemma 3.2 Let W : G → R be a continuous function such that sup n∈N ν n (W ) < ∞, a.s. Let φ ∈ C 1 (G), be such that lim |x|→∞ |∇φ(x)| 2 /W (x) = 0. Then 1 √ Λ n n k=1 λ k ∇φ(X k-1 ), σ(X k-1 )U k L -→ N 0, G |σ T ∇φ| 2 dν .
Proof. This lemma follows from Theorem 1.3 using the martingale central limit theorem, along the lines of Proposition 2 of [START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion[END_REF]. Details are left to the reader.

(3/2) n P -→ 1 6 G R m D 3 φ(x)(σ(x)u) ⊗3 µ(du)ν(dx), (4) n Λ
as n → ∞.

Proof. The proof is similar to that of Lemma 10 of [START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion[END_REF] except for the treatment of reflection terms. Using the notation above Theorem 1.5 and in Lemma 3.1, we have

R 2 (x, y) = 1 6 D 3 φ(x)(y -x) ⊗3 + R 4 (x, y), (44)
with

|R 4 (x, y)| ≤ L 6 |y -x| 4 ,
where L is the Lipschitz constant for

D 3 φ. Hence R 2 (X k-1 , X k ) = 1 6 D 3 φ(X k-1 )(δX k ) ⊗3 + r k , (45)
with

|r k | ≤ L 6 |δX k | 4 ≤ κ 1 λ 2 k (1 + |U k | 4 ), k ∈ N, for some κ 1 ∈ (0, ∞). Since E|U k | 4 := µ 4 < ∞ from Condition 1.6, we have E n k=1 |r k | ≤ κ 1 (1 + µ 4 ) n k=1 λ 2 k .
From the assumption lim n→∞

(1/ √ Λ n ) n k=1 λ 3/2 k = λ ∈ (0, +∞], we deduce that lim n→∞ n k=1 λ 3/2 k = +∞ and lim n→∞ n k=1 λ 2 k /Λ (3/2) n = 0. (46)
Therefore, 1

Λ (3/2) n n k=1 r k L 1 --→ 0. ( 47
)
Now consider the first term on the right side of (45).

D 3 φ(X k-1 )(δX k ) ⊗3 =D 3 φ(X k-1 )(λ k b(X k-1 ) + λ k σ(X k-1 )U k + y k-1 ) ⊗3 =λ 3/2 k D 3 φ(X k-1 )( λ k b(X k-1 ) + σ(X k-1 )U k ) ⊗3 + f (1) k (X k-1 , U k ) =λ 3/2 k D 3 φ(X k-1 )(σ(X k-1 )U k ) ⊗3 + f (2) k (X k-1 , U k ) + f (1) k (X k-1 , U k ), where f (1) k (X k-1 , U k ) and f (2)
k (X k-1 , U k ) are defined through the second and third equalities, respectively. Next observe that • From the assumptions, we have

|λ k b(X k-1 ) + λ k σ(X k-1 )U k | ≤ a 1 λ k (|U k | + λ 0 ).
• From ( 28) and ( 33), we have

y k-1 = N d i L i
k-1 and for some κ 2 ∈ (0, ∞),

L i k-1 ≤ κ 2 λ k (|U k | + 1), for all k ∈ N • The term L i k-1 is non zero only if there exists s ∈ [0, 1] such that x k-1 (s) ∈ F i , where x k-1 (s) = x(U k , X k-1 , λ k |s).
And in that case, we have from [START_REF] Weiss | Invariant measures of diffusions in bounded domains[END_REF], the Lipschitz property of D 3 φ and ( 11) that, for some κ 3 ∈ (0, ∞),

|D 3 •jk φ(X k-1 ) • d i | ≤ κ 3 λ k (|U k | + 1), ∀j, k.
Combining these estimates, we see that

E n k=1 |f (1) k (X k-1 , U k )| ≤ κ 4 n k=1 λ 2 k . Using (46) we now have 1 Λ (3/2) n n k=1 f (1) k (X k-1 , U k ) L 1 --→ 0. (48)
For the term f

k (X k-1 , U k ), using the boundedness of D 3 φ, b, and σ, it can be easily checked that E|f

(2) k (X k-1 , U k )| ≤ κ 5 λ 2 k . Thus E n k=1 |f b (X k-1 , U k )| ≤ κ 5 n k=1 λ 2 k ,
and so using (46) once again, we have 1

Λ (3/2) n n k=1 f b (X k-1 , U k ) P -→ 0. (49) Let Θ(X k-1 , U k ) = D 3 φ(X k-1 )(σ(X k-1 )U k ) ⊗3 . Since sup k E|Θ(X k-1 , U k )| 2 < ∞ and lim n→∞ Λ (3) n /(Λ (3/2) n ) 2 = 0, we have 1 Λ (3/2) n n k=1 λ 3/2 k [Θ(X k-1 , U k ) -E(Θ(X k-1 , U k )|F k-1 )] L 2 --→ 0. ( 50
) Observe that E(Θ(X k-1 , U k )|F k-1 ) = J(X k-1
), where J is given by

J(x) := R m D 3 φ(x)(σ(x)u) ⊗3 µ(du). Since Λ (3/2) n → ∞ as n → ∞, we can apply Theorem 1.3 to the measure νn = 1 Λ (3/2) n n k=1 λ 3/2 k δ X k-1 .
Since J is continuous and bounded, we have lim n→∞ νn (J) = Jdν a.s., and the lemma follows on combining this fact with (44)-(50).

2

We are now ready to prove Theorem 1.5.

Proof of Theorem 1.5 The proof is similar as the proof of Theorem 9 of [START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion[END_REF], once again the main difference is in the treatment of reflection terms. Using the notation of Lemma 3.1, we first observe that, for any sequence of positive numbers {a n } n∈N such that lim n→∞ a n = ∞, we have Z (0) n /a n → 0 Suppose that |J| = p and J = {i 1 , ..., i p }. Define a p × p matrix R J be the relation R

J (k, l) = (π J d i l ) i k , for k, l = 1, ...p. Let u J , v J ∈ R p be the solution of LCP for (R J , π J θ), i.e., (u J , v J ) = (L 1 p (R J , π J θ), L 2 p (R J , π J θ)).
Once again unique solvability of LCP for (R J , π J θ) is assured from Condition 1.1. Denote u J = (η 1 , ..., η p ) and define x 1

(t) = x + θt + t p j=1 η j d ij . Let τ 1 = inf{t ≥ 0|In(x 1 (t)) = In(x)}.
We define τ 1 = ∞ if the above set is empty. Then Γ(x + θi)(t) = x 1 (t) for all t < τ 1 . If τ 1 < ∞ set the initial point to be x 1 = x 1 (τ 1 ) and define the trajectory {x 2 (t)} t≥0 in a similar way as {x 1 (t)} by replacing x with x 1 . Set τ 2 = inf{t ≥ 0|In(x 2 (t)) = In(x 1 )}. Then Γ(x + θi)(τ 1 + t) = Γ(x 1 + θi)(t) = x 2 (t), for all t < τ 2 .

Define now recursively trajectory {x j (t)} with time points τ j and end points x j (τ j ), j = 3, 4, .... Let j 0 be such that Thus the evaluation of S(x, θ) can be carried out by recursively solving a sequence of LCP problems.

One difficulty in implementing the above scheme is the possibility that ∞ i=1 τ i ≤ 1. However using regularity property of the Skorokhod map, we see that this occurs only when S(x, θ) is zero. Thus in the practical implementation of the algorithm we fix a finite threshold L and carry out the above recursive procedure at most L times and set S(x, θ) = 0 if L i=1 τ i < 1. It can be checked that the spectral radius of Q is less than 1, and so Conditions 1.1 and 1.2 hold. Take the drift function b(x) = [-1/2, -1/2, -1/2] T and σ(x) = I, x ∈ R 3 + . The stationary distribution ν for this example is of product form (see [START_REF] Harrison | Brownian models of open queueing networks with homogeneous customer populations[END_REF]): exp(1.1667) ⊗ exp(1.0938) ⊗ exp(0.8537), where exp(µ) is the exponential distribution with parameter µ. In implementing the above numerical scheme, we set our initial point to X 0 = [1, 1, 1] T , and simulate {X k } n k=1 defined by equation ( 8), taking U k ∼ N (0, I), λ k = 1/ √ k and n = 10 7 . Figure 1 shows the comparison between the exact distribution with the first-coordinate marginal of the measure ν n . This example was considered in [START_REF] Dai | Reflected Brownian motion in an orthant: Numerical methods for steady-state analysis[END_REF]. We consider the first moment of the x 1 -coordinate. The exact value for this moment is known to be 0.5. We consider λ k = k -α and examine the influence of the choice of α on the numerical performance. The results are given in Figure 2. We find that α = 0.5 gives the best numerical convergence.

An 8-d symmetric SRBM.

A SRBM is said to be symmetric if its covariance matrix Γ, drift vector µ and reflection matrix R are symmetric in the following sense: Γ ij = Γ ji = ρ for 1 ≤ i < j ≤ d, µ i = -1 for 1 ≤ i ≤ d and R ij = R ji = -r for 1 ≤ i < j ≤ d, where r ≥ 0 . The positiveness of Γ implies -1/(d -1) < ρ < 1 and the completely-S condition of R implies r(d -1) < 1. In this case, It is known (see [START_REF] Dai | Reflected Brownian motion in an orthant: Numerical methods for steady-state analysis[END_REF]) that, the first moment of each of the component is the same, and is given by the following formula We consider time step sequence λ n = n -α with different choice of α and study the influence of α on numerical convergence. The thin solid line, the dotted line, the thick solid line, the dash-dot line, and the dashed line correspond to α =0.1, 0.3, 0.5, 0.7, and 0.9 respectively. The x-axis shows the value of n while the y-axis corresponds to x 1 ν n (dx).

Here we take d = 8. Then the conditions on the data yield -1/7 < ρ < 1 and 0 ≤ r < 1/7. Letting ρ range through {-0.1, -0.05, 0, 0.2, 0.9}, and r take value 0.1, we obtain estimates of m 1 using algorithm in this work. We take λ k = k -α , α = 0.5 and n = 10 7 . The results are shown in Table 1. The results show that as the correlation coefficient ρ approaches 1, the performance of the algorithm deteriorates.

For

  η ∈ D([0, ∞) : R m ) let |η|(T ) denote the total variation of η on [0, T ] with respect to the Euclidean norm on R m . Definition 1.1 Let ψ ∈ D G ([0, ∞) : R m ) be given. Then the pair (φ, η) ∈ D([0, ∞) : R m ) × D([0, ∞) : R m ) solves the Skorokhod problem (SP) for ψ with respect to G and d if and only if φ(0) = ψ(0), and for all t ∈ [0, ∞)

L

  -→" to denote convergence in distribution. Let R m denote the set of m-dimensional real vectors. Euclidean norm will be denoted by | • | and the corresponding inner product by •, • . The symbols, P -→, L p

Theorem 1 . 5 √

 15 Assume that U i 's are i.i.d with common distribution µ. There exists a ζ ∈ (0, ∞) such that whenever φ ∈ C 2 (G) satisfies lim |x|→∞ e -ζ|x| |∇φ(x)| 2 = 0, we have the following: (a) Fast-decreasing step. Suppose lim n→∞ Λ (3/2) n Λn = 0, D 2 φ is bounded and Lipschitz, and

2 2. 3

 23 r > 0, where F ǫ,r i = {x ∈ R m | dist(x, F i ) ≥ ǫ and |x| ≤ r}. The equality in (41) now follows on sending ǫ → 0 and r → ∞. Proof of Theorem 1.4 Recall c from Lemma 2.1 and ̟ from

2 Lemma 3 . 3

 233 Under the assumptions of Theorem 1.5(b), we have, Z

j0 i=1 τ i < 1 ≤

 1 j0+1 i=1 τ i . Then S(x, θ) = Γ(x + θi)(1) = Γ(x j0 + θi)(1 -j0 i=1 τ i ).

4. 2

 2 Results. 4.2.1 A 3-d Example with Product Form Stationary Distribution. Let m = 3 and suppose that the reflection matrix is of the form R = I + Q, where I is the identity matrix, and Q is given as

4. 2 . 2

 22 Effect of Choice of {λ k }. Consider a two-dimensional SRBM with covariance matrix σ(x) = I, drift vector b(x) = [-1, 0] T and reflection matrix

m 1 = 1 -

 11 (d -2)r + (d -1)rρ 2(1 + r).

Figure 1 :

 1 Figure 1: Comparison between the nnumerically computed distribution with exact distribution. The left figure shows the comparison between the empirical cumulative distribution function (cdf) and the exact cdf. The middle figure makes a comparison between the estimated density function and the exact density function. And the right figure is the qq-plot of the empirical quantiles versus the exact quantiles.

Figure 2 :

 2 Figure 2:We consider time step sequence λ n = n -α with different choice of α and study the influence of α on numerical convergence. The thin solid line, the dotted line, the thick solid line, the dash-dot line, and the dashed line correspond to α =0.1, 0.3, 0.5, 0.7, and 0.9 respectively. The x-axis shows the value of n while the y-axis corresponds to x 1 ν n (dx).

A. Budhiraja, J. Chen and S. Rubenthaler: Numerical Scheme for Invariant Distributions of Constrained Diffusions Mathematics of Operations Research xx(x), pp. xxx-xxx, c 200x INFORMS

in probability. This is because, from Lemma 2.4, the sequence {X n } n∈N is tight, and consequently so is {φ(X n )} n∈N as well.

We also derive from the definitions of Z [START_REF] Atar | On positive recurrence of constrained diffusion process[END_REF] n , Z

n and Z [START_REF] Budhiraja | An ergodic control problem for constrained diffusion processes: Existence of optimal markov control[END_REF] n the inequalities E|Z (1) n

and

for some κ 1 ∈ (0, ∞), for all n ≥ 1.

(a) Now assume that lim n→∞ (1/ √ Λ n )Λ

(3/2) n = 0. We then have lim n→∞ n k=1 λ 2 k / √ Λ n = 0, and it follows from (51) that Z

(1)

--→ 0. We also deduce from (52), that Z

n . Denoting the Lipschitz norm of D 2 φ by L, we have

where the second inequality follows from the Lipschitz property of the Skorokhod map (Condition 1.1). Thus, there exists κ 2 ∈ (0, ∞) such that, for all n ≥ 1,

and therefore Z

We now, consider Y

n , for j = 1, 2,

Y

From (33), we have

. Also, for any fixed i, L i k-1 is not equal to 0 only if there exists

and in that case, using Taylor's theorem and the Lipschitz property of D 2 φ, there exists κ 4 ∈ (0, ∞), such that,

Combining this with (10), we have

Thus we have

E|Y (1) n | ≤ κ 5

for some constant κ 5 . Using similar arguments as above, we obtain:

Thus we have that Y

From Lemma 2.5, and recalling the definition of V (see ( 18)) we have that, for every ζ ∈ (0, c̟),

For such a ζ, under the assumption that lim |x|→∞ e -ζ|x| |∇φ(x)| 2 = 0, applying Lemma 3.2, we now have

This completes the proof of part (a). --→ 0, and from (52) that

--→ 0, for j = 2, 3.

Under the assumptions of part (b) (i.e. that D 3 φ is bounded, Lipschitz and (11) holds), we have, using similar arguments as in part (a), for some κ 6 ∈ (0, ∞),

Thus we have that Y

--→ 0, for j = 1, 2, 3, 4.

Applying Lemma 3.2 once again, we have, for φ satisfying lim |x|→∞ e -ζ|x| |∇φ(x)| 2 = 0,

Also from Lemma 3.3

Now, if λ < +∞, we have from the above observations that Z

n / √ Λ n P -→ 0, for j = 1, 2, 3, 4 and

The statement in [START_REF] Dai | A multiclass station with Markovian feedback in heavy traffic[END_REF] now follows on combining this with (57).

Finally, if λ = +∞, we have Z

Numerical Results

Evaluation of the Euler Time

Step. A key step in simulating the sequence {X k } in ( 8) is the evaluation of S(X k , Y k+1 -X k ), where S : G × R m → G is the time-1 Skorokhod map defined in [START_REF] Budhiraja | Long time asymptotics for constrained diffusions in polyhedral domains[END_REF]. In this section we describe a procedure for computing S(x, v) that uses well known relationships between Skorokhod problems and linear complementary problems (LCP). We restrict ourselves to a setting where N = m and G = R m + . We begin by recalling the basic formulation of the LCP (see [START_REF] Cottle | The linear complementarity problem[END_REF]). For j ∈ N, a j × j matrix R and a j-dimensional vector θ, the LCP for (R, θ) is to find vectors u, v ∈ R j such that

It is well known (see [START_REF] Dupuis | Convex duality and the Skorokhod Problem. I[END_REF] and [START_REF] Budhiraja | Simple necessary and sufficient conditions for the stability of constrained processes[END_REF]) that with R = [d 1 , ..., d m ], under Condition 1.1, for every θ ∈ R m , the LCP for (R, θ) admits a unique solution (u, v) ≡ (L 1 m (R, θ), L 2 m (R, θ)), and furthermore L 2 m (R, θ) = S(0, θ). Thus the evaluation of S(0, θ) reduces to solving the above LCP for which numerous algorithms are available. In the examples considered in the current work we used a quadratic programming algorithm. Evaluation of S(x, θ) for x = 0 can be carried out using a localization procedure as follows.

Fix x ∈ G and let J = In(x) = {j ∈ {1, ..., m}| x, e j = 0}. Let P J = {z ∈ R m | z, e j = 0, ∀j ∈ J c }. Let π J : R m → P J be the orthogonal projection:

z, e j e j . Using the convexity of the function e λx , we have

Taking expectations in the above inequality and using Taylor's expansion, we have

2 .

The lemma then follows on taking α = 1 2 . 2


Mayer-Rokitansky-Küster-Hauser Syndrome Karine Morcel

Bruno Dallapiccola

Laurent Pasquier

Tanguy Watrin

Laura Bernardini

Dr Daniel Guerrier
email: daniel.guerrier@univ-rennes1.fr

Rokitansky-Küster-Hauser Syndrome

Karine Morcel

CUGC for Mayer

Disease characteristics 1.1 Name of the Disease (Synonyms):

Mayer-Rokitansky-Küster-Hauser (MRKH) syndrome, MRKH type I or isolated MRKH or Rokitansky sequence, MRKH type II or MURCS association (Müllerian duct aplasia, Renal dysplasia and Cervical Somite anomalies), Congenital Absence of the Uterus and Vagina (CAUV), Genital Renal Ear Syndrome (GRES), Müllerian Aplasia (MA).

OMIM# of the Gene(s):

Putative candidate genes : 189907 (TCF2) [1,2], 601999 (LHX1) [2][START_REF] Cheroki | Genomic imbalances associated with mullerian aplasia[END_REF][4][START_REF] Ledig | Recurrent aberrations identified by array-CGH in patients with Mayer-Rokitansky-Kuster-Hauser syndrome[END_REF], 312865 (SHOX) [6], 602427 (TBX6) [2], 609783 (ITIH5) [START_REF] Morcel | Utero-vaginal aplasia (Mayer-Rokitansky-Kuster-Hauser syndrome) associated with deletions in known DiGeorge or DiGeorge-like loci[END_REF]. Putative candidate locus : 274000 (TAR) [START_REF] Cheroki | Genomic imbalances associated with mullerian aplasia[END_REF][START_REF] Ledig | Recurrent aberrations identified by array-CGH in patients with Mayer-Rokitansky-Kuster-Hauser syndrome[END_REF].

Mutational Spectrum:

Maximal deletion in 1q21.1: 398.5 kb [START_REF] Cheroki | Genomic imbalances associated with mullerian aplasia[END_REF][START_REF] Ledig | Recurrent aberrations identified by array-CGH in patients with Mayer-Rokitansky-Kuster-Hauser syndrome[END_REF]; in 4q34-qter: 8 Mb [START_REF] Bendavid | Phenotypic variability of a 4q34-->qter inherited deletion: MRKH syndrome in the daughter, cardiac defect and Fallopian tube cancer in the mother[END_REF]; in 8p23.1: 1.2 Mb [START_REF] Morcel | Utero-vaginal aplasia (Mayer-Rokitansky-Kuster-Hauser syndrome) associated with deletions in known DiGeorge or DiGeorge-like loci[END_REF]; in 10p14-15: 0.2 Mb [START_REF] Morcel | Utero-vaginal aplasia (Mayer-Rokitansky-Kuster-Hauser syndrome) associated with deletions in known DiGeorge or DiGeorge-like loci[END_REF]; in 16p11.2: 600kb [2]; in 17q12: 1.5 Mb [2,[START_REF] Cheroki | Genomic imbalances associated with mullerian aplasia[END_REF][START_REF] Ledig | Recurrent aberrations identified by array-CGH in patients with Mayer-Rokitansky-Kuster-Hauser syndrome[END_REF]

Analytical Methods:

Search for microrearrangements by means of Multiplex Ligation-dependent Probe Amplification (MLPA) using the SALSA P023 DiGeorge MLPA kit, Comparative Genomic Hybridization array (CGH array), Duplex PCR/Liquid Chromatography (DP/LC), and/or FISH

Analytical Validation

Validation of MLPA results or CGH array results by DP/LC or FISH 1.8 Estimated Frequency of the Disease (Incidence at birth or "birth prevalence" or population prevalence):

1 in 4500 female births [START_REF] Folch | Mullerian agenesis: etiology, diagnosis, and management[END_REF] 1.9 If applicable, prevalence in the ethnic group of investigated person: Comment: WNT4 syndrome is close but different from MRKH. It differs from this latter by an hyperandrogenism and gonadal affection [START_REF] Biason-Lauber | WNT4 and sex development[END_REF]. It thus needs to be considered in regards to differential diagnosis [START_REF] Morcel | MRKH) syndrome[END_REF].

$ $ $ $

Test characteristics 2.1 Analytical Sensitivity

(proportion of positive tests if the genotype is present)

Nearly 100% using analytical methods described above.

Analytical Specificity

(proportion of negative tests if the genotype is not present)

Nearly 100% using analytical methods described above.

Clinical Sensitivity (proportion of positive tests if the disease is present)

The clinical sensitivity can be dependent on variable factors such as age or family history. In such cases a general statement should be given, even if a quantification can only be made case by case.

The synthesis of the literature [2-5, 9, 10] on the genome-wide analysis on MRKH patients, was used to estimate that only ~1% show a 1q21.1 rearrangement, ~1% in 16p11.2, ~6% in 17q12 and ~4% in 22q11.21, from a total of 153 cases.

Clinical Specificity

(proportion of negative tests if the disease is not present)

The clinical specificity can be dependent on variable factors such as age or family history. In such cases a general statement should be given, even if a quantification can only be made case by case.

100% of control subjects tested by MLPA (DiGeorge kit) showed negative results.

Positive clinical predictive value

(life time risk to develop the disease if the test is positive). not applicable

Negative clinical predictive value

(Probability not to develop the disease if the test is negative). Assume an increased risk based on family history for a non-affected person. Allelic and locus heterogeneity may need to be considered.

Index case in that family had been tested: not applicable Index case in that family had not been tested: not applicable

Clinical Utility

(Differential) diagnosis: The tested person is clinically affected

(To be answered if in 1.10 "A" was marked)

Can a diagnosis be made other than through a genetic test?

No.

(continue with 3.

How is the cost effectiveness of alternative diagnostic methods to be judged?

Alternative diagnosis methods described above are, at present, the only way to accurately diagnose MRKH syndrome since the phenotype-genotype correlations still cannot be established. Clinical examination, imaging, and biological measurements, will in any case remain necessary to evaluate the level of care for patients.

Will disease management be influenced by the result of a genetic test?

No.

Yes. Therapy (please describe) Prognosis (please describe) Management (please describe) The results of genetic tests will influence genetic counselling in case of surrogate pregnancy demand.

Predictive Setting: The tested person is clinically unaffected but carries an increased risk based on family history

(To be answered if in 1.10 "B" was marked)

Will the result of a genetic test influence lifestyle and prevention?

If the test result is positive (please describe)

If the test result is negative (please describe)

Which options in view of lifestyle and prevention does a person at-risk have if no genetic test has been done (please describe)?

$ $ $ $ $

Genetic risk assessment in family members of a diseased person

(To be answered if in 1.10 "C" was marked)

Does the result of a genetic test resolve the genetic situation in that family?

Can improve genetic counselling In rare familial forms

Prenatal diagnosis

(To be answered if in 1.10 "D" was marked)

If applicable, further consequences of testing

Please assume that the result of a genetic test has no immediate medical consequences. Is there any evidence that a genetic test is nevertheless useful for the patient or his/her relatives? (Please describe)

1. 2 1 . 3

 213 OMIM# of the Disease: 277000 (MRKH, CAUV), 601076 (MURCS) Name of the Analysed Genes or DNA/Chromosome Segments: 1q21.1, 4q34-qter, 8p23.1, 10p14-15, 16p11.2, 17q12, 22q11.21, Xpter-p22.32

  and in 22q11.21: 3 Mb [2, 3, 5, 7, 9, 10]. Maximal duplication in 1q21.1: 200 kb [3]. Partial duplication of the Xpter pseudoautosomal region 1 [6].

3. 3 . 2 No 3 . 3 . 3

 32333 Can a genetic test in the index patient save genetic or other tests in family members? Does a positive genetic test result in the index patient enable a predictive test in a family member?

3. 4 . 1

 41 Does a positive genetic test result in the index patient enable a prenatal diagnosis?

2 Describe the burden of alternative diagnostic methods to the patient

  Initial clinical examination and transabdominal ultrasonography must be the first investigations in evaluating patients with suspected Müllerian aplasia. If necessary, magnetic resonance imaging (MRI) affords to clearly precise the malformation. A full check-up (transabdominal ultrasonography, spine radiography, and/or heart echography, audiogram) should be undertaken to search for associated malformations in order to determine the MRKH type (I or II). Testosterone dosage should help orientating the diagnosis towards MRKH or WNT4 syndrome.

	1.4)
	Yes,
	clinically.
	imaging.
	endoscopy.
	biochemistry.
	electrophysiology.
	other (please describe): Testosterone dosage
	3.1.

Acknowledgement

This work was supported by EuroGentest, an EU-FP6 supported NoE, contract number 512148 (EuroGentest Unit 3: "Clinical genetics, community genetics and public health", Workpackage 3.2)

Conflict of interest

The authors declare no conflict of interest

Clinical utility gene card for:

Karine Morcel 1,2 , Bruno Dallapiccola 3 , Laurent Pasquier 4 , Tanguy Watrin 1,5 , Laura Bernardini 6 and Daniel Guerrier 1,5