Alexis Anagnostakis
email: alexis.anagnostakis@univ-lorraine.fr

Antoine Lejay
email: antoine.lejay@univ-lorraine.fr

Denis Villemonais
email: denis.villemonais@university.edu

General diffusion processes as the limit of time-space Markov chains

Keywords: classification: 60J60, 60F17 Markov chain approximation, Markov chain approximation, random walk, singular diffusion, sticky, skew, slow reflection, Donsker's invariance principle, Wasserstein distance

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the diffusion process literature, the most well-studied and straightforward way to approximate diffusion processes is the Euler scheme. While such approximations works well for non-degenerate stochastic differential equations, this is not the case for more general diffusion processes [START_REF] Hutzenthaler | Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF]. The Euler scheme is also not well-defined for processes that exhibit sticky features, skew behavior [START_REF] Frikha | On the weak approximation of a skew diffusion by an Euler-type scheme[END_REF] or slowly reflecting boundaries.

First studied by Feller [START_REF] Feller | The parabolic differential equations and the associated semi-groups of transformations[END_REF][START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF], general diffusion processes are used to replicate a wide range of phenomena like semi-permeable layers [START_REF] Lejay | Monte Carlo methods for fissured porous media: a gridless approach[END_REF], principal agent problem dynamics [START_REF] Piskorski | Optimal dynamic contracts with moral hazard and costly monitoring[END_REF][START_REF] Zhu | Optimal contracts with shirking[END_REF] and interest rates behavior close to 0 [START_REF] Nie | Sticky reflecting ornstein-uhlenbeck diffusions and the vasicek interest rate model with the sticky zero lower bound[END_REF]. For more theoretical results and applications of general diffusions, we refer the reader to [START_REF] Bass | A stochastic differential equation with a sticky point[END_REF][START_REF] Eberle | Sticky couplings of multidimensional diffusions with different drifts[END_REF][START_REF] Engelbert | Stochastic differential equations for sticky Brownian motion[END_REF][START_REF] Hajri | Application of stochastic flows to the sticky Brownian motion equation[END_REF][START_REF] Lejay | On the constructions of the skew Brownian motion[END_REF]. The resurgence of interest in these processes and the need to have a universal simulation method for diffusions have motivated the search for new kinds of approximation processes.

Several works aim at overcoming the shortcomings of the Euler scheme and allow us to approximate the law of more general diffusion processes. In [START_REF] Amir | Sticky Brownian motion as the strong limit of a sequence of random walks[END_REF], the author proposes to approximate the sticky Brownian motion with a simple random walk that stops for a fixed amount of time when it hits 0. In [START_REF] Nie | Sticky reflecting ornstein-uhlenbeck diffusions and the vasicek interest rate model with the sticky zero lower bound[END_REF] and [START_REF] Meier | Markov chain approximation of one-dimensional sticky diffusions[END_REF], Continuous Time Markov Chains are used to approximate slowly reflected SDE solutions, where the jumping intensities are computed using approximated discretizations of the infinitesimal generator of the diffusion. Another work where such processes are defined is [START_REF] Ferrer-Admetlla | An approximate markov model for the wright-fisher diffusion and its application to time series data[END_REF], where the authors use a Continuous Time Markov Chain to identify events in genomics evolution. In [START_REF] Ankirchner | Wasserstein convergence rates for coin tossing approximations of continuous markov processes[END_REF][START_REF] Ankirchner | A functional limit theorem for coin tossing Markov chains[END_REF], the authors develop a numerical scheme to approximate diffusions on natural scale as long as a mild non-explosion condition is satisfied. They use symmetric random walks with fixed-time step whose magnitude depends on the average local behavior of the target process speed measure. Choosing the step-size this way allows the approximation process to replicate non-boundary sticky features.

In this paper, we prove the convergence in law of grid-valued random walks to any one-dimensional general diffusion process at an asymptotically optimal rate. General diffusion processes are regular one-dimensional strong Markov processes with continuous trajectories, see for instance [START_REF] Revuz | Continuous martingales and Brownian motion, volume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF]Chapter 7.3] where they are defined as linear continuous Markov processes. This convergence result allows us to set up approximation schemes that, while make it straightforward to take into account for sticky points, can also be applied to any diffusion process that satisfy a mild non-explosion condition. This includes processes with boundary conditions like absorption, reflection or slow-reflection as well as the skew diffusions such as the Skew Brownian motion [START_REF] Lejay | On the constructions of the skew Brownian motion[END_REF] and its generalizations. The values taken by the random walk correspond to values taken by the target process at random times, allowing us to classify it as an embeddable scheme along with [START_REF] Ankirchner | A functional limit theorem for coin tossing Markov chains[END_REF] and [START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF]. We prove that for a grid adapted to the speed measure of the diffusion process, the laws of the random walks converges at any rate strictly inferior to (1/2) ∧ (2/p) in terms of the maximum cell size for all p-Wasserstein distances. This convergence rate is optimal for p ≤ 4 according the Donsker invariance principle, as this is the rate simple randoms walk converges to the standard Brownian motion [START_REF] Donsker | An invariance principle for certain probability limit theorems[END_REF].

Besides the asymptotic optimal convergence rate, the usage of such an approximation process yield several advantages. Firstly, the static character of the grid makes involved quantities good candidates for numerical approximation (see Sections 6 and 7). Moreover, this scheme makes it straightforward to take into account potential sticky points of the diffusion. Finally, its universality is further validated by the fact that the Donsker invariance principle and [START_REF] Amir | Sticky Brownian motion as the strong limit of a sequence of random walks[END_REF][START_REF] Donsker | An invariance principle for certain probability limit theorems[END_REF][START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF] are all special cases of it.

We make heavy use of the speed measure/scale function characterization of diffusions: from Proposition 3.14 of [START_REF] Revuz | Continuous martingales and Brownian motion, volume 293 of Grundlehren der Mathematischen Wissenschaften[END_REF]Chapter VII], the law of a one-dimensional diffusion process with state-space I an open interval of R is entirely determined by an increasing continuous function s and a positive locally finite measure m defined on I. The function s called the scale function, is defined as the unique1 function such that, (

We can also express the infinitesimal generator of the diffusion L in terms of m and s:

L = D m D s , dom(L) = {f ∈ C 0 b (I, R) : L f ∈ C 0 b (I, R)},
where D m an D s are defined by

D m f (x) = lim h→0;h>0 f (x + h) -f (x) m(x, x + h] , D s f (x) = lim h→0;h>0 f (x + h) -f (x) s(x + h) -s(x) .
These results give us analytic formulations of quantities of the form v k (x) = E x (τ k ab 1 τ b <τa ) for k ∈ N 0 and allow us to bound these quantities in terms of the size of the interval (a, b). The latter will be particularly useful for proving the convergence results of the algorithm.

Outline. In Section 2, we present the approximation scheme along with its properties. In Section 3, we give analytical characterizations of the quantities that determine the law of the random walk defined by the algorithm, allowing us to implement it. Section 4 is dedicated to proving the convergence of embedding times. In Section 5, we prove the main convergence result in terms of the maximum cell size of the grid. The case of solution of an SDE is studied in Section 6. Section 7 is dedicated to numerical experiments. In Section 7.1, we apply the scheme for several types of diffusions (reflection, stickiness, skew-behavior, singularity at 0 or ∞) and illustrate its convergence. In Section 7.2, we exhibit the flexibility of STMCAs on a local time approximation problem.

2 The Space-Time Markov Chain Approximation and its properties

The approximation scheme

In this section we define the approximation process for a one-dimensional diffusion process on natural scale2 with state space I, an open interval of R. The general case is obtained by a change of scale, as detailed in Section 2.2.

The possible values taken by the approximation process are given as input of the scheme and must form a covering grid of I. We introduce incrementally this notion as follows, which we illustrate by Figure 1. Let I be an interval of R.

• A grid over I is a countable subset of I with no accumulation points within I.

• A cell c of a grid g is an open interval with endpoints in g with a single element of g in its interior, i.e., Card g ∩ c = 1.

We denote with C(g) the set of all cells of the grid g.

• We call x ∈ I the center of the cell c ∈ C(g) iff c ∩ g = {x}.
• Finally, a covering grid of I is a grid g such that I = c∈C(g) c.

Examples of covering grids of (0, ∞) and (-∞, ∞) are {1/n; n ∈ N} and Z respectively.

For any covering grid g of I and diffusion process X with state-space I, defined through (s, m), let |g|, |g| X be the grid metrics:

|g| = sup c∈C(g) |c|, |g| X = sup c∈C(g) {s(c)m(c)}, (2)
where |c| = (ba) and s(c) = s(b)s(a) with a and b being the endpoints of c. The convergence results will be expressed in the latter metric.

-5 -4 -3 -2 -1 0 1 2 3 4 5

Figure 1: The covering grid g = Z of (-∞, ∞) along with some of its cells C(g).

Let X be a diffusion process on natural scale with state-space I and speed measure m. Let X be defined on the family of filtered probability spaces (Ω, F, (F t ) t≥0 , (P x ) x∈I ), where P x is the law of X such that, for any x ∈ I, P x {X 0 = x} = 1. For any covering grid g of I, we define the process X g = ( X g t ) t≥0 as the asymmetric random walk with: • state-space g, • initial distribution equal to the distribution of X the first time it touches the grid, • the same transition probabilities as X over g, • conditional transition times that match the conditional expected transition times of X over g. Thus, under P x , if a and b are respectively the closest lower and upper elements to x of g,

X g 0 = a, with probability P x (τ b < τ a ), b, with probability P x (τ a < τ b ) = 1 -P x (τ b < τ a ), (3)
where τ a := inf{t > 0 : X t = a}. For the rest of the trajectory, we define τ ab := τ a ∧ τ b and (T g (n)) n≥0 as the consecutive jumping times of ( X g t ) t≥0 . Then, for all k ∈ N 0 and a < x < b adjacent points of g,

P X g T g (k+1) = b X g T g (k) = x = P x (τ b < τ a ), (4)
and

T g (k + 1) -T g (k) = E x (τ ab |τ b < τ a ), on { X g T g (k+1) = b} ∩ { X g T g (k) = x}, E x (τ ab |τ a < τ b ), on { X g T g (k+1) = a} ∩ { X g T g (k) = x}. (5)
As proved in Section 3, the quantities that appear on the right hand side of (4) and ( 5) are explicit functionals of the speed measure m.

Let c x be the cell of the grid g containing x, i.e., c x ∈ C(g) and x = c x ∩ g. From (4), [START_REF] Ankirchner | Wasserstein convergence rates for coin tossing approximations of continuous markov processes[END_REF] and Bayes' rule, if c x = (a, b), both P X g T g (k+1) = b X g T g (k) = x and T g (k + 1) -T g (k) only depend on x. Thus, if we know the quantities

p + [x, (a, b)] = P x (τ b < τ a ), T + [x, (a, b)] = E x τ ab τ b < τ a , p -[x, (a, b)] = P x (τ a < τ b ), T -[x, (a, b)] = E x τ ab τ a < τ b , (6)
for any adjacent a < x < b in g, we can simulate the random walk using Algorithm 1. We discuss in Section 3 on how to compute the quantities in [START_REF] Ankirchner | A functional limit theorem for coin tossing Markov chains[END_REF]. Practical examples are given in Sections 6 and 7.

This algorithm has been first introduced in [START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF] in the situation of SDE solutions with measurable coefficients, where the speed measure of the process satisfies

c dx ≤ m(dx) ≤ C dx. (7)
Our main contribution is that we allow non-elliptic speed measures with vanishing and unbounded density (as in e.g. the Bessel process case, see Section 7.1), speed measures with singular part (as in e.g. the sticky Brownian motion setting, see Sections 7.1 and 7.2), scale functions not in C 1 (as in e.g. the skew Bessel process, see Section 7.1) and non-trivial boundary behaviors (see Section 2.4). The probabilistic arguments we use to prove our results allow for greater flexibility, while the proofs of [START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF] are based on elliptic PDE theory. This allows us to handle degenerate diffusions and to perform grid tuning and achieve higher orders of convergence (see Section 2.3).

Remark 1. We observe that the process X g is not a one-dimensional Markov chain. It is though a Markov chain in space and time since the joint law of the next position of X g on the grid and the next transition time are both determined by the current position on the grid. Hence the terms: space-time Markov chains and Space Time Markov Chain Approximation (STMCA).

Remark 2. In the case of sticky diffusions, where the speed measure m has the form m(dx) = m c (dx) + ρδ 0 (dx), the transition probabilities and transition times (6) can be directly inferred from the ones of the diffusion without the sticky term. Indeed, Proposition 3.2 yields

E x τ ab 1 τ b <τa = (a,b) G a,b (x, ζ)v 0 (ζ)m c (dζ) + ρ G a,b (x, 0)v 0 (0), (8)
where v 0 (x) = P x (τ b < τ a ).

For the convergence, we make the further assumption that the diffusion process (X t ) t≥0 satisfies the following non-explosion condition: there exists a k 1 > 0 such that the speed measure of the diffusion process satisfies

m(dx) ≥ k 1 1 1 + x 2 dx, (9)
for all x ∈ I. Practically, this means that the process does not move faster than a log-normal process for large values of X t . We may now express the convergence result in terms of

Initialization: t[0] = 0, n = 0 j = arg min i∈J {|x i -x|} U ∼ Bernoulli p + [x, (x j , x j+1 )] if U == 1 then j = j + 1 end if X[0] = x j Main loop: while t[n] < T do U ∼ Bernoulli p + [x j , (x j-1 , x j+1 )] if U == 1 then j = j + 1 t[n + 1] = t[n] + T + [x j , (x j-1 , x j+1 )] else j = j -1 t[n + 1] = t[n] + T -[x j , (x j-1 , x j+1 )] end if X[n + 1] = x j n = n + 1 end while
the step-size of the grid in terms of p-Wasserstein distances. In the following result, the p-Wasserstein distance W p between two laws µ and ν of processes with càdlàg3 paths is defined as

W p µ, ν = inf (ζ,ξ)∼Γ(µ,ν) ∥ζ -ξ∥ ∞ L p ,
where by Γ(µ, ν) we denote the collection of all measures with marginals µ and ν.

Theorem 2.1. Let X be a diffusion process with state-space I an interval of R, on natural scale, whose speed measure satisfies Condition (9) for some constant k 1 > 0. Let g be a covering grid of I. Then, for all p ≥ 1, δ ∈ (0, 1 4 ∧ 1 p ), T > 0 and x ∈ I there exists a constant C > 0 such that

W p Law ( X g t ) t∈[0,T ] , Law (X t ) t∈[0,T ] ≤ C|g| δ X , (10)
where |g| X = sup c∈C(g) {|c|m(c)}.

Remark 3. In the case where m(dx) ≥ k 1 dx, the constant4 C > 0 in Theorem 2.1 does not depend on the starting point of the diffusion.

Remark 4. If X is a diffusion process on natural scale such that (7) holds, the bound in [START_REF] Chung | On the zero n 1 ±1[END_REF] can be replaced by C|g| 2δ .

The convergence of the Wasserstein distances implies the convergence in law [31, p. 109].

Corollary 2.2. For all T > 0, the processes ( X g t ) t∈[0,T ] converges in law to (X t ) t∈[0,T ] in the Skorokhod space D([0, T ], I) as |g| X -→ 0.

In Section 2.5, we observe that, by suitable choice of the probability space, the values taken by the approximation process correspond to values taken by X. We call this class of approximation schemes embeddable schemes (other embeddable schemes are the ones developed in [START_REF] Ankirchner | A functional limit theorem for coin tossing Markov chains[END_REF][START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF]). Proving the convergence of an embeddable scheme usually involves: embedding the approximation process in the trajectory of the target diffusion process and controlling the embedding times, controlling the speed at which the process moves.

Convergence rate for the general case

The convergence results established in the previous section are proven in the case of a diffusion process on natural scale. In this section, we show how more general results can be inferred. Let X be a diffusion process with state-space I an open interval of R, scale function s and speed measure m. We assume that • s belongs to the Sobolev space W 1,1 (I), so from Theorem 8.2 of [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF], as s is continuous,

s(x) -s(y) = x y s ′ (t) dt,
for all y, x in I.

• there exists a k 1 > 0 such that for all x ∈ I,

m(dx) ≥ k 1 s ′ (x) 1 + s(x) 2 dx, (11)
• the inverse of s is α-Hölder continuous, i.e., there exists a constant C > 0 such that for all x ̸ = y ∈ I,

|s -1 (x) -s -1 (ȳ)| |x -ȳ| α ≤ C. (12)
Given a grid g, we consider the random walk X g t defined by Algorithm 1, where the transition probabilities and transition times in ( 6) can be computed using the formulas derived in Section 3. We obtain the following corollary of Theorem 2.1.

Corollary 2.3. Let X be a diffusion process with scale function and speed measure satisfying the above conditions. Let also g be a covering grid over the state-space I of X. Then, for all p ≥ 1, δ ∈ (0, 1 4 ∧ 1 p ), T > 0 and x ∈ I there exists positive constants C 1 , and C 2 such that

W p Law ( X g t ) t∈[0,T ) , Law (X t ) t∈[0,T ) ≤ C 1 |g| δ X ,
where |g| X is defined in (2).

Proof. We define the proxy process Y = (s(X t )) t≥0 which has state-space s(I), scale function

s Y (x) = x and speed measure m Y (dx) = m • s -1 (dx).
From condition [START_REF] Cox | A theory of the term structure of interest rates[END_REF] and a change of variables, we get that (Y t ) t≥0 satisfies condition (9) for the same constant k 1 . We also define Y s(g) t evolving according to Algorithm 1, with covering grid s(g) = {s(x); x ∈ g}. It can be defined on the canonical space of X so that s(

X g t ) = Y s(g) t
almost surely. Thus, Condition [START_REF] Donsker | An invariance principle for certain probability limit theorems[END_REF] implies that

|X t -X g t | ≤ C|Y t -Y s(g) t | α .
Along with the fact that

|g| Y = sup c∈C(g) {|s(c)|m Y (s(c))},
Theorem 2.1 implies Corollary 2.3.

Grid tuning

We observe that for all ϵ > 0, in the case of a Brownian motion, Theorem 2.1 yields a convergence rate of O(|g| 1 2 -ϵ ) as |g| -→ 0, which is optimal from Donsker's invariance principle [START_REF] Donsker | An invariance principle for certain probability limit theorems[END_REF]. We would like to have this result for all diffusion processes, but the following example illustrates that this is not the case. We then show how we can remediate to this by using a custom grid and extrapolate this method to the general case via Corollary 2.5.

Example 2.4. Let X be the diffusion process with state-space R, defined through s and m with s

(x) = x, m(dx) = 2 dx + ρδ 0 (dx).
This process is called the sticky Brownian motion and is the "most elementary" sticky diffusion process. As such, it spends a positive amount of time at 0 and the Euler scheme is known to not be well defined for these processes. We observe that, for any covering grid g of R and c ∈ C(g) with c ̸ = c 0 and c 0 ∈ C(g) being the cell containing 0,

2|c| 2 = m(c)|c| ≤ sup c∈C(g) m(c)|c|, ρ|c 0 | + 2|c 0 | 2 = m(c 0 )|c 0 | ≤ sup c∈C(g) m(c)|c|. (13)
From ( 13), for a uniform grid of step size h,

2hρ < |g| X .
Thus, for any ϵ > 0, the convergence rate given by Theorem 2.1 is O(h 1 4 -ϵ ) as |g| -→ 0. This means that there are functionals of the trajectory for which the convergence rate is much slower for this process in comparison with a standard Brownian motion. In order to remediate to this, we propose a preliminary step to the approximation scheme that involves finding a grid that is "adapted" to the speed measure of the process. In the case of the Brownian motion with a sticky point at 0, such a grid can be defined as one that has uniform non-adjacent cells to 0 of size h and with the cell of center 0 being (-h 2 /2ρ, h 2 /2ρ), i.e.

g = k∈Z + - h 2 2ρ -k h 2 ∪ 0 ∪ k∈Z + h 2 2ρ + k h 2 . ( 14
)
As the approximation process is a random walk, for every k steps it makes, it spends O( √ k) steps in the cell containing 0 (see [START_REF] Chung | On the zero n 1 ±1[END_REF]). Thus, running the algorithm on either grid yields the same algorithmic complexity, whilte the convergence rate is improved to O(h The general case is covered by the following Corollary: Corollary 2.5. Let X be a diffusion process and X g the approximation process defined by Algorithm 1. Then, if g is a grid such that

|g| X ≤ C|g| 2 , ( 15
)
we can bound the p-Wasserstein distance between the laws of (

X g t ) t∈[0,T ) and (X t ) t∈[0,T ) in Theorem 2.1 by |g| 2δ instead of |g| δ X .
Thus, for all ϵ > 0, the law of the random walk converges in any p-Wasserstein distance at the rate O(|g|

( 1 2 ∧ 2 p )-ϵ ) instead of O(|g| ( 1 4 ∧ 1 p )-ϵ ) as |g| -→ 0.
We now show how one can create grids such that (15) holds in the case of homogeneous SDEs. Considering sticky and/or skew points is straightforward.

• For a process whose speed measure satisfies [START_REF] Bass | A stochastic differential equation with a sticky point[END_REF] and [START_REF] Cox | A theory of the term structure of interest rates[END_REF] with one or more skew points, no grid modification is required,

• For a process whose speed measure satisfies ( 7) and ( 11) and has a sticky point at 0 of stickiness ρ > 0, one needs to consider the points {-h 2 /ρ, 0, h 2 /ρ} to have a tuned grid,

• The case of a reflection at a boundary is treated in Section 2.4.

Let (µ, σ) be a pair of functions satisfying the following condition:

Condition 1. The functions are measurable R → R mappings and the SDE

dX t = µ(X t ) dt + σ(X t ) dB t , (16)
has a unique weak solution, where B is a standard Brownian motion.

Let X be the diffusion that solves ( 16), I its state-space, s its scale function and m its speed measure given by [8, p. 17]

s ′ (x) = e -y a 2µ(ζ) σ 2 (ζ)
dζ , and m(dx

) = 1 s ′ (x) 2 σ 2 (x) dx, (17)
with s ′ being the right-derivative of s. Then, if g = {x j } j∈J is a covering grid such that for a constant C > 0 and every j ∈ J s(x j+1 )s(x j-1 )

x j+1

x j-1

2 s ′ (ζ)σ 2 (ζ) dζ ∧ (x j+1 -x j-1 ) 2 ≤ Ch 2 , (18)
it satisfies [START_REF] Étoré | A Donsker theorem to simulate one-dimensional processes with measurable coefficients[END_REF] and |g| = Ch. Thus, from Corollary 2.5, for any ϵ > 0 using such grids give us a convergence rate of O(|g|

1 2 -ϵ ) in Theorem 2.
1. Generating such grids numerically can be done choosing a starting point x 0 and adding points x j to the grid iteratively as follows: given x j-1 , let x j be chosen such that:

s(y) -s(x j-1 ) y x j-1 2 s ′ (ζ)σ 2 (ζ) dζ ≤ h 2 /2. ( 19
)
Then the next element of the grid is defined as

x j = x (0) j if x (0) j -x j-1 ≤ h, x j-1 + h if x (0) j -x j-1 > h. (20)
The problem [START_REF] Ferrer-Admetlla | An approximate markov model for the wright-fisher diffusion and its application to time series data[END_REF] can be solved numerically using a fixed point algorithm. Examples of STMCA simulations using tuned grids computed solving ( 19)- [START_REF] Frikha | On the weak approximation of a skew diffusion by an Euler-type scheme[END_REF] are given in Figures 4 and5.

The case of diffusions with boundary conditions

When presenting the results and the structure of the scheme, we considered only processes where I is an open set, thus considering diffusion with unreachable boundaries. Our results also adapt to the situation where either ℓ and/or r are reachable, and in this case some adjustments are needed, depending on the nature of finite boundaries and on the condition at regular boundaries. In order to keep the presentation simple, we assume that the process is on natural scale and that I = [0, +∞) (the adaptation to

I = (ℓ, r] or I = [ℓ, r] or I = [ℓ, r) with ℓ ∈ R ∪ {-∞} and r ∈ R ∪ {+∞} is straightforward).
It is well known (see for instance Section 5.11 of Itô's book [START_REF] Itō | Essentials of stochastic processes[END_REF]) that the finite boundary 0 can be of four types. Setting, for some fixed c > 0,

I = 0<y<x<c m(dx) dy, II = 0<y<x<c m(dy) dx, then • 0 is an exit boundary if I < ∞ and II = ∞, • 0 is a regular boundary if I < ∞ and II < ∞, • 0 is a natural boundary if I = ∞ and II = ∞, • 0 is an entrance boundary if I = ∞ and II < ∞.
The entrance type been excluded for a finite boundary of a diffusion process on natural scale, and the natural type been considered in the settings of Theorem 2.1, this leaves us with two possible types for the boundary 0: exit or regular. If 0 is an exit boundary, then the diffusion process X is absorbed at the boundary 0. If 0 is a regular boundary, then the diffusion process can either be absorbed or reflected at 0. In these cases, the convergence result of Theorem 2.1 can be extended by considering a grid g on I containing 0 and by adapting the dynamics of X g as follows. The dynamic of X g is the same as in Algorithm 1, up to the time when it reaches 0, then:

• if 0 is an absorbing boundary (exit or regular), then the result can be immediately extended by stopping X g when it reaches 0;

• if 0 is a reflecting regular boundary, then the process X g jumps from 0 to b := min g \{0} with probability 1 and after a time

[0,b) (b -ζ) m(dζ).
We emphasize that in this configuration, 0 may be a sticky boundary (i.e., with m(0) ∈ (0, +∞)).

In both cases of reflection and absorption at 0, the boundary is attainable. Thus, 0 must be a point of any covering grid of I = [0, ∞). This leads to the following adaptation of the notion of grid cells. The cells C(g) of such grid g are the open intervals for the induced topology of R on I with endpoint in g. For example The proof of the convergence in these situations is omitted here, since it is a straightforward adaptation of the proof of Theorem 2.1, using in particular the fact that, in the case of a reflecting boundary,

g 0 = Z + is a covering grid of [0, ∞) and C(g 0 ) = {[0, 1), (n -1, n + 1) n∈N }.
E 0 τ b = [0,b) (b -ζ) m(dζ).
The case of killing boundaries, and in general of a process with non-zero killing measure, leads to additional non-trivial difficulties. Devising an algorithm and a similar result as Theorem 2.1 for such processes remains an active area of research.

Markovian embedding

The consecutive values of the process ( X t ) t≥0 defined through (3)-( 5) form a Markov chain with, by construction, the same transition probabilities as (X t ) t≥0 on g. We define the embedding times of (X t ) t≥0 in g as,

τ g 0 = 0, τ g k = inf t > τ g k-1 : X t ∈ g \ {X τ g k-1 } , ∀k ≥ 1. (21)
As both X T g (n) and X τ g n are both Markov chains with the same transition probabilities with X 0 forced to be equal in law to X τ g 1 (see Section 2.1), the following equality in law holds,

Law( X T g (n) ; n ≥ 0) = Law(X τ g n ; n ≥ 1). ( 22
)
We define K g (t) as the inverse of T g (n), i.e.,

K g (t) = inf n ∈ N : n k=1 E τ g k -τ g k-1 X τ g k-1 , X τ g k > t . (23)
Thus, we get the following Proposition.

Proposition 2.6. Let (X t ) t≥0 be a diffusion process, g a grid defined over its state space I and X g t be the approximation process defined in (4) and (5). Then, if (τ g n ) n≥0 are the embedding times of (X t ) t≥0 in g, the following equality in law holds, Law X t ; t ≥ 0 = Law X τ g K g (t) ; t ≥ 0 , where K g (t) is the random index defined in [START_REF] Itō | Essentials of stochastic processes[END_REF].

Moment characterization of conditional exit times

The law of the approximation process defined in the previous section was shown to be determined by the transition probabilities P x (τ b < τ a ) and conditional transition times E x (τ ab |τ b < τ a ). In this section we show that quantities of the form

v k (x) = E x (τ k ab 1 τ b <τa )
yield an integral formulation with respect to the speed measure of the diffusion and involving the scale function (we do not assume that the diffusion is on natural scale in the present section). We also show that this results in them being solutions to Dirichlet problems where the differential operator is the infinitesimal generator L of the diffusion. This allows us to simulate such processes via Algorithm 1 and thus to approximate the law of the target diffusion process (X t ) t≥0 .

In terms of Algorithm 1, we need to compute for three adjacent points a, x, b of the grid the quantities

v 0 (x) = P x (τ b < τ a ), v 1 (x) = E x (τ ab 1 τ b <τa ) and v 1 (x) = E x (τ ab 1 τa<τ b ).
The quantities of (6) are then

p + [x, (a, b)] = v 0 (x), T + [x, (a, b)] = v 1 (x) v 0 (x) , p -[x, (a, b)] = 1 -v 0 (x), T -[x, (a, b)] = v 1 (x) 1 -v 0 (x) . (24)
Proposition 3.1. The function v 0 (x) = P x (τ b < τ a ) is solution to the problem with Dirichlet boundary conditions

     L u = 0, x ∈ (a, b), u(a) = 0, u(b) = 1, ( 25
)
where L is the infinitesimal generator of (X t ) t≥0 , which also implies that v 0 ∈ dom(L).

Proof. Let x ∈ (a, b), from the definition of the scale function and the factorization of the infinitesimal generator

L = D m D s L v 0 = D m D s s(•) -s(a) s(b) -s(a) = D m 1 s(b) -s(a)
.

which equals 0 as m(dx) is a positive measure. As v 0 and L v 0 = 0 are both functions in C 0 b , we deduce that v 0 ∈ dom(L) and L v 0 = 0. Under P b , the stopping time τ b equal 0 a.s. and the process has a.s. continuous trajectories, hence τ a > 0 a.s., i.e.,

v 0 (b) = P b (τ b < τ a ) = P b (0 < τ a ) = 1.
This, along with the symmetrical argument, allow us to retrieve the boundary conditions of [START_REF] Lejay | On the constructions of the skew Brownian motion[END_REF]. Proposition 3.2. For every k ∈ N, let v k be the function defined for every x ∈ (a, b) by

v k (x) = E x (τ k ab 1 τ b <τa ). Then, v k (x) = k (a,b) G a,b (x, ζ)v k-1 (ζ)m(dζ). (26)
Proof. Since

τ ab 0 (τ ab -t) k-1 dt = τ k ab /k, v k (x) = k E x 1 τ b <τa τ ab 0 (τ ab -t) k-1 dt = k E x 1 τ b <τa ∞ 0 1 t≤τ ab (τ ab -t) k-1 dt .
From the Markov property, by conditioning on F t and as

1 t≤τ ab is F t -measurable, v k (x) = k E x ∞ 0 1 t≤τ ab E 1 τ b <τa (τ ab -t) k-1 F t dt = k E x τ ab 0 E Xt 1 τ b <τa τ k-1 ab dt .
The equality [START_REF] Meier | Markov chain approximation of one-dimensional sticky diffusions[END_REF] results by applying directly Green's formula.

1 y<x (s(y) -s(a))(s(b) -s(x)) s(b) -s(a) + 1 y≥x (s(x) -s(a))(s(b) -s(y)) s(b) -s(a) f (y)m(dy) = -D m (a,x) v 0 (y)f (y)m(dy) + D m [x,b) (1 -v 0 (y))f (y)m(dy) = -v 0 (x)f (x) -(1 -v 0 (x))f (x) = -f (x).
The continuity of g is a consequence of Lebesgue's convergence theorem for integrals. Moreover as G a,b (x, y) is bounded by s(b)s(a), m(dx) is locally finite and f is bounded, g is also bounded. So, we deduce that f ∈ dom(L) and that on (a, b) we have L g = -f .

Proof of Proposition 3.3. From Proposition 3.2,

v k (x) = (a,b)
G a,b (x, y)kv k-1 (y)m(dy).

As v 0 ∈ C 0 b (a, b), from Lemma 3.4 we deduce iteratively that v k ∈ C 0 b (a, b) for all k ∈ N and that L v k = -kv k-1 on (a, b).
For the boundary conditions, we observe that τ ab = τ a ∧ τ b = 0 a.s. under P a , so for k ≥ 1

v k (a) = E a (τ k ab 1 τ b <τa ) ≤ E a (τ k ab ) = 0.
With the same argument we show that v k (b) ≤ 0 and as they are obviously positive quantities

v k (a) = v k (b) = 0.

Convergence of the embedding times

In order to prove the convergence of the process ( X g t ) t≥0 we need to control quantities of the form E x sup t≤T |tτ g K g (T ) | p , where τ g k are the embedding times of the process X in the grid g. In this section, we show the existence of such bounds in terms of the metric

|g| X = sup c∈C(g) {s(c)m(c)} (27)
defined in [START_REF] Alili | On the semi-group of a scaled skew Bessel process[END_REF]. If there exists a constant K > 0 such that m(c) ≤ K for all c ∈ C(g), then, |c|m(c) ≤ K|c| and |g| X ≤ K|g|. Thus, we can bound the quantities of interest in terms of |g| instead of |g| X . But, in doing so, we do not track correctly the convergence rates of the approximation process. For example, for the standard Brownian motion, |g| B = |g| 2 . Moreover, as shown in Section 2.3, such bounds give us a direct way to adapt the grid to the speed measure in order to accelerate the convergence of the scheme.

Bounds on the conditional moments of the exit times

Lemma 4.1. Let v k (x) = E x (τ k ab 1 τ b <τa ), then for all k ∈ N v k v k-1 ∞ ≤ k|g| X . (28)
Proof. We first observe that

G a,b (x, y) v 0 (y) v 0 (x) = (y-a)(b-x) (b-a) y-a x-a , x > y, (x-a)(b-y) (b-a) y-a
x-a , x ≤ y.

As for x > y the ratio y-a x-a < 1,

G a,b (x, y) v 0 (y) v 0 (x) ≤ (y-a)(b-x) (b-a) , x > y, (y-a)(b-y) (b-a) , x ≤ y, (29)
which is bounded by (ba) in both cases as x, y ∈ (a, b). Thus,

b a G a,b (x, y)v 0 (y)m(dy) = v 0 (x) b a G a,b (x, y) v 0 (y) v 0 (x) m(dy) ≤ v 0 (x)(b -a)m((a, b)). (30) From Proposition 3.2, v k (x) = k (a,b) G a,b (x, y)v k-1 (y)m(dy) and v k (x) = k! (a,b) G a,b (x, x k ) (a,b) G a,b (x k , x k-1 )• • • (a,b) G a,b (x 2 , x 1 )v 0 (x 1 )m(dx 1 ) . . . m(dx k ).
(31) From ( 30) and [START_REF] Villani | Optimal transport. Old and new[END_REF],

v k (x) ≤ k(b -a)m((a, b))v k-1 (x). Since (b -a)m((a, b)) ≤ |g| X , we get the desired result on v k /v k-1 .
Corollary 4.2. Let k ∈ N. Then, we have the following bound for

v k (x) = E x (τ k ab 1 τ b <τa ), ∥v k ∥ ∞ ≤ k!|g| k X . Lemma 4.3.
Let X be a diffusion process with state space I an interval of R, m the speed measure of X, a, x, b ∈ I such that a < x < b, c x = (a, b), |c x | X = m(c x )|b -a| and λ > 0 such that λ|c x | X ∈ (0, 1). Then,

E x (e λτ ab ) ≤ exp λ 1 -λ|c x | X E x (τ ab ) . (32)
Proof. Developing the exponential series,

E x (e λτ ab ) = 1 + ∞ N =1 λ N N ! E x (τ N ab ) = 1 + E x (τ ab ) ∞ N =1 λ N N ! E x (τ N ab ) E x (τ ab )
.

From Corollary 4.2, we can bound the ratio of expected values by

(N -1)!|c x | N -1 X . As λ|c x | X ∈ (0, 1), E x (e λτ ab ) ≤ 1 + E x (τ ab ) ∞ N =1 λ N (λ|c x | X ) N -1 ≤ 1 + E x (τ ab ) λ 1 -λ|c x | X .
Thus, we only need to apply the inequality 1 + x ≤ e x to get [START_REF] Zhu | Optimal contracts with shirking[END_REF].

Lemma 4.4. Let t, M > 0 and λ > 0 such that λ|g| X ∈ (0, 1). Then,

P x (τ g K g (t) > M ) ≤ e -λM e λ t+|g| X 1-λ|g| X . (33)
Proof. From Markov's inequality,

P x (τ g K g (t) > M ) ≤ e -λM E x e λτ g K g (t) = e -λM E x e λ K g (t) k=1 (τ k -τ k-1 ) . (34)
Conditioning on the σ-algebra B generated by the trajectory of X t on the grid g, i.e., B = σ{X τ k ; k ∈ N 0 }, and as K g (t) is B-measurable,

E x e λ K g (t) k=1 (τ k -τ k-1 ) = E x K g (t) k=1 E e λ(τ k -τ k-1 ) B . ( 35
)
From the definition of |g| X , λ|c| X ≤ λ|g| X ∈ (0, 1) for each cell c of the grid g. Thus, applying Lemma 4.3 on each term of the product in (35),

P x (τ g K g (t) > M ) ≤ e -λM E x exp λ 1 -λ|g| X K g (t) k=1 E(τ k -τ k-1 |B) (36)
From [START_REF] Itō | Essentials of stochastic processes[END_REF],

t < K g (t) k=1 E(τ k -τ k-1 |B) ≤ t + E(τ K g (t)+1 -τ K g (t) |B)
Thus, from (28),

K g (t) k=1 E(τ k -τ k-1 |B) ≤ t + ∥v 1 /v 0 ∥ ∞ ≤ t + |g| X , (37)
From ( 36) and (37), we get (33).

Convergence of the embedding times

For this section, let X be a diffusion process with state-space I an interval of R, g a covering grid of I, τ g k the embedding times of X in g as defined in ( 21) and K g (t) as defined in [START_REF] Itō | Essentials of stochastic processes[END_REF].

Lemma 4.5. For any T > 0,

K g (T ) k=1 Var τ g k -τ g k-1 X τ g k-1 , X τ g k ≤ 2|g| X T + |g| X , where |.| X is defined in (2).
Proof. For all x ∈ g, let c x be the cell of g containing x. Then,

K g (T ) k=1 Var τ g k -τ g k-1 X τ g k-1 , X τ g k = K g (T ) k=1 E (τ g k -τ g k-1 ) 2 X τ g k-1 , X τ g k -E τ g k -τ g k-1 X τ g k-1 , X τ g k 2 ≤ sup k≤K g (t) E (τ g k -τ g k-1 ) 2 X τ g k-1 , X τ g k E τ g k -τ g k-1 X τ g k-1 , X τ g k K g (T ) k=1 E (τ g k -τ g k-1 ) X τ g k-1 , X τ g k .
So from Lemma 4.1 and the definition of K g (t),

K g (T ) k=1 Var τ g k -τ g k-1 X τ g k-1 , X τ g k ≤ v 2 v 1 ∞ K g (T ) k=1 E (τ g k -τ g k-1 ) X τ g k-1 , X τ g k ≤ v 2 v 1 ∞ T + v 1 v 0 ∞ ≤ 2|g| X (T + |g| X ),
which is the desired inequality.

Proposition 4.6. Let (F t ) t≥0 be the canonical filtration of X. Let also

A n = F τ g n , B = σ (X τ g k ) k∈N 0 and ∆τ g k = τ g k -τ g k-1
. If we define the augmented filtration G n = A n ∨ B, then the process

M n = n k=1 ∆τ g k -E ∆τ g k X τ g k-1 , X τ g k , is a G n -martingale. Proof. For m ≤ n, E M n G m = E M n A m , B = E m k=1 ∆τ g k -E ∆τ g k X τ g k-1 , X τ g k + n k=m+1 ∆τ g k -E ∆τ g k X τ g k-1 , X τ g k A m , B .
From Lemma A.1 and as

τ g k is A k measurable, E M n G m = M m + n k=m+1 E X τ g m ∆τ g k -E ∆τ g k X τ g k-1 , X τ g k B = M m .
This proves the result.

Theorem 4.7. For any T > 0,

E sup t∈[0,T ] |τ g K g (t) -t| 2 ≤ 2|g| X 4 T + |g| X + 1 , (38)
where |.| X is defined in (2).

Proof. The convexity inequality (a + b) p ≤ 2 p-1 (a p + b p ) yields for p = 2

E sup t∈[0,T ] |τ g K g (t) -t| 2 ≤ 2 E sup t∈[0,T ] τ g K g (t) - K g (t) k=1 E ∆τ g k X τ g k-1 , X τ g k 2 + 2 E sup t∈[0,T ] K g (t) k=1 E ∆τ g k X τ g k-1 , X τ g k -t 2 . (39)
From Proposition 4.6, the term

n k=1 ∆τ g k -E ∆τ g k X τ g k-1 , X τ g k is a G n -martingale, where G n = F τ g n ∨ B and B = σ (X τ g k ) k∈N 0 . Thus, from Doob's L p inequality, E sup t∈[0,T ] |M K g (t) | 2 = E sup k≤K g (T ) |M k | 2 ≤ 2 E |M K g (T ) | 2 .
By conditioning on B, from Lemma A.1,

E sup t∈[0,T ] K g (t) k=1 ∆τ g k -E ∆τ g k X τ g k-1 , X τ g k 2 ≤ 2 E K g (T ) k=1 ∆τ g k -E ∆τ g k X τ g k-1 , X τ g k 2 = 2 E K g (T ) k=1 Var ∆τ g k X τ g k-1 , X τ g k ,
which from Lemma 4.5 is bounded by 4|g| X T + |g| X . For the second term on the right hand side of (39), from [START_REF] Itō | Essentials of stochastic processes[END_REF] and since K g (t) ≥ 1 for any t > 0,

K g (t) k=1 E τ g k -τ g k-1 X τ g k-1 , X τ g k -t ≤ E τ g K g (t) -τ g K g (t)-1 X τ g K g (t)-1 , X τ g K g (t) ≤ v 1 v 0 ∞ ≤ |g| X .
So having bounded both additive parts of the right hand side of (39), we get (38).

5 Convergence rate of the Markov chain

Moment bounds

In this section we prove that Theorem 3.1 of [START_REF] Ankirchner | Wasserstein convergence rates for coin tossing approximations of continuous markov processes[END_REF] holds also for reflected processes and for a sharper constant. This result, combined with the bound (38) is used to prove the convergence of the approximation process in Section 5.2.

Lemma 5.1. Let X be a diffusion process on natural scale with state-space R and a speed measure m X that satisfies Condition [START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]. Then, for all p ≥ 2, there exist two constants C,

C ′ > 0 such that E x sup t∈[0,T ] |X t -x| p ≤ C ′ [1 + |x| p ]e CT , ( 40
)
where k 1 is a constant such that Condition (9) is satisfied, C ≤ 8p(p -1)/k 1 and C ′ > 0 is a constant that depends only on p.

Proof. Let Z be the diffusion process on natural scale with speed measure,

m Z (dx) = 1 |x|<1 k 1 2 dx + 1 |x|≥1 k 1 2x 2 dx. We note that k 1 2x 2 ≤ k 1 1+x 2 for all x ≥ 1.
The dynamic of Z can be shown to be

dZ t =      2 √ k 1 dB t , for |Z t | < 1, 2 √ k 1 Z t dB t , for |Z t | ≥ 1,
where B is a standard Brownian motion. We also assume that P x (Z 0 = x) = 1. As X and Z are on natural scale, they can be expressed as time-changed Brownian motion [30, Theorem 47.1, p. 277], i.e., for every t ≥ 0, X t = B γ X (t) and Z t = W γ Z (t) , where B and W are two standard Brownian motions with P x (B 0 = x) = P x (W 0 = x) = 1, respective local times L x (B), L x (W ) and with γ X (t), γ Z (t) being the respective right-inverses5 of

A X (t) = 1 2 I L x t (B)m X (dx) and A Z (t) = 1 2 I L x t (W )m Z (dx).
Using the same underlying Brownian motion in these definitions, we have A Z (t) ≤ A X (t), and hence γ X (t) ≤ γ Z (t). Thus,

E x sup t∈[0,T ] |X t -x| p = E x sup t∈[0,γ X (T )] |B t -x| p ≤ E x sup t∈[0,γ Z (T )] |B t -x| p = E x sup t∈[0,T ] |Z t -x| p . (41)
Thus, from Doob's L p and convexity inequalities,

E x sup t∈[0,T ] |X t -x| p ≤ 2 p-1 p p -1 E x |Z T | p + |x| p . ( 42
)
Using Itô's formula, followed by a classical localization argument, Fatou's Lemma along with the standard dominated convergence theorem, one obtains that for all q > 2/3,

E x (|Z t | -1) 3q 1 |Zt|≥1 + 1 ≤ (|x| -1) 3q 1 |x|≥1 + 1 + t 0 3q 2 (3q -1) E x (|Z s | -1) 3q-2 1 |Zs|≥1 4Z 2 s k 1 ds. Using the inequality |x -1| 3q-2 x 2 ≤ 4|x -1| 3q + 4 for all x ≥ 1, E x (|Z t | -1) 3q 1 |Zt|≥1 + 1 ≤ (|x| -1) 3q 1 |x|≥1 + 1 + t 0 24q k 1 (3q -1) E x (|Z s | -1) 3q 1 |Zs|≥1 + 1 ds.
Using Gronwall's Lemma, we deduce that, for all t ≥ 0,

E x (|Z t | -1) 3q 1 |Zt|≥1 + 1 ≤ [(|x| -1) 3q 1 |x|≥1 + 1]e 24q(3q-1)t/k 1 .
Let C ′ q > 0 be a constant such that (|x| -1) 3q 1 |x|≥1 + 1 ≥ C ′ q |x| 3q for all x ∈ R. Thus,

E x |Z t | 3q ≤ C ′ q [(|x| -1) 3q 1 |x|≥1 + 1]e 24q(3q-1)t/k 1 , ∀t ≥ 0.
Hence, for all 6 p > 2, we have

E x |Z t | p ≤ C ′ p/3 [1 + |x| p ]e 8p(p-1)t/k 1 , ∀t ≥ 0. (43)
Replacing (43) in (42), we get the bound (40) for

C ′ = p p-1 2 p-1 (C ′ p/3 +1) and C = 8p(p-1)/k 1 .
Proposition 5.2. Let X be a diffusion process with state-space I an interval of R, on natural scale and with a speed measure m X that satisfies Condition (9). Then, for each T > 0 and γ ∈ (0, 1 2 ), there exists a constant C > 0 such that

sup s̸ =t≤T |X t -X s | |t -s| γ L p (Px) ≤ C(1 + |x|). (44)
The result also holds in the presence of a reflecting boundary ζ ∈ I. 6 The extension to p = 2 is straightforward.

Proof. For the non-reflecting case, the proof works using the same arguments as in the proof of Theorem 3.1 in [START_REF] Ankirchner | Wasserstein convergence rates for coin tossing approximations of continuous markov processes[END_REF].

For the reflecting case: Let X be a diffusion process on natural scale, with speed measure m satisfying (11) for a constant k 1 > 0 and a reflecting boundary at ζ ∈ I. We observe that

X = |X • -ζ| + ζ (45)
in law, where X • is the non-reflecting diffusion on natural scale with speed measure

m X • (dx) = 1 x≥ζ m(dx) + 1 x<ζ m(2ζ -dx)
which also satisfies [START_REF] Cox | A theory of the term structure of interest rates[END_REF] for the same constant k 1 . From (45) and the triangle inequality,

|X t -X s | |t -s| γ = |X • t -ζ| -|X • s -ζ| |t -s| γ ≤ |X • t -X • s | |t -s| γ . ( 46
)
The diffusion process X • is non-reflecting and its speed measure m X • satisfies [START_REF] Cox | A theory of the term structure of interest rates[END_REF]. Thus, (44) holds for X • for a constant C > 0 and from (46),

sup s̸ =t≤T |X t -X s | |t -s| γ L p (Px) ≤ sup s̸ =t≤T |X • t -X • s | |t -s| γ L p (Px) ≤ C(1 + |x|).

Proof of the convergence rate for a process on natural scale

Proof of Theorem 2.1. From Proposition 2.6,

W p Law ( X g t ) t∈[0,T ] , Law (X t ) t∈[0,T ] = inf ∥ζ -ξ∥ L p ; ζ ∼ Law ( X g t ) t∈[0,T ] , ξ ∼ Law (X t ) t∈[0,T ] ≤ sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px)
.

Let M > T . From Minkowski inequality,

sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) ≤ 1 τ g K g (T ) ≤M sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) + 1 τ g K g (T ) >M sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) (47)
For the first additive term of (47), for any γ > 0, by multiplying and dividing by |τ g K g (t) -t| γ ,

1 τ g K g (T ) ≤M sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) = 1 τ g K g (T ) ≤M sup t∈[0,T ] |X τ g K g (t) -X t | |τ g K g (t) -t| γ |τ g K g (t) -t| γ L p (Px) ≤ 1 τ g K g (T ) ≤M sup t∈[0,T ] |X τ g K g (t) -X t | |τ g K g (t) -t| γ sup t∈[0,T ] |τ g K g (t) -t| γ L p (Px)
.

As τ g K g (t) is increasing with respect to t,

1 τ g K g (T ) ≤M sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) ≤ sup s̸ =t≤M |X s -X t | |s -t| γ sup t∈[0,T ] |τ g K g (t) -t| γ L p (Px) = E x sup s̸ =t≤M |X s -X t | |s -t| γ sup t∈[0,T ] |τ g K g (t) -t| γ p 1/p
. From Hölder's inequality for q ≥ 1 and q/(q -1) conjugates exponents and (44), for every γ ∈ (0, 1 2 ), there exists a constant C 1 = C 1 (M, γ, p(q -1)/q, x) > 0 such that

1 τ g K g (T ) ≤M sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) ≤ E x sup s̸ =t≤M |X s -X t |
|s -t| γ p(q-1)/q q/p(q-1)

E x sup t∈[0,T ] |τ g K g (t) -t| γ pq 1/pq ≤ C q/p(q-1) 1 E x sup t∈[0,T ] |τ g K g (t) -t| γpq 1/pq .
Choosing q = 2/γp, from (38), for any γ ∈ (0, 1 2 ∧ 2 p )

1 τ g K g (T ) ≤M sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) ≤ C 1/p(1-γp) 1 2|g| X (4T + |g| X ) γ/2 . ( 48
)
For the second additive term of (47),

E x 1 τ g K g (T ) >M sup t∈[0,T ] |X τ g K g (t) -X t | p = ∞ m=M E x 1 τ g K g (T ) ∈[m,m+1) sup t∈[0,T ] |X τ g K g (t) -X t | p . ( 49
)
For each term of the sum in (49), from Hölder's inequality,

E x 1 τ g K g (T ) ∈[m,m+1) sup t∈[0,T ] |X τ g K g (t) -X t | p ≤ P x (τ g K g (T ) > m) 1/q ′ E x 1 τ g K g (T ) ∈[m,m+1) sup t∈[0,T ] |X τ g K g (t) -X t | p q ′ /(q ′ -1) ≤ P x (τ g K g (T ) > m) 1/q ′ E x 1 τ g K g (T ) <m+1 sup t∈[0,T ] |X τ g K g (t) -X t | p q ′ /(q ′ -1)
.

As each m in the sum in (49) satisfies m + 1 > M > T , from Minkowski's inequality, m+1) , where C 2 and C 3 are positive constants depending on p. This, along with Lemma 4.4 and Hölder's inequality gives us for λ > 0 chosen such that α = λ|g| X < 1,

E x 1 τ g K g (T ) <m+1 sup t∈[0,T ] |X τ g K g (t) -X t | p ≤ E x 1 τ g K g (T ) <m+1 2 p-1 sup t∈[0,T ] |X τ g K g (t) -x| p + |X t -x| p ≤ 2 p E x sup t∈[0,m+1] |X t -x| p , which from Lemma 5.1 is bounded by C 2 [1 + |x| p ]e C 3 (
E x 1 τ g K g (T ) ∈[m,m+1) sup t∈[0,T ] |X τ g K g (t) -X t | p ≤ C 2 [1 + |x| pq ′ ]e C 3 (m+1) 1/q ′ e -λm e λ T +|g| X 1-λ|g| X (q ′ -1)/q ′
where C 2 (pq ′ ) ≤ 2pq ′ (pq ′ -1)/c and C 2 (pq ′ ) > 0 a positive constant depending only on pq ′ .

Thus, setting C 4 (pq

′ ) := C 2 (pq ′ ) [1 + |x| pq ′ ] 1/q ′ e C 3 /q ′ > 0, E x 1 τ g K g (T ) ∈[m,m+1) sup t∈[0,T ] |X τ g K g (t) -X t | p ≤ C 4 exp 1 q ′ C 3 m + (q ′ -1) λ T + |g| X 1 -λ|g| X -λm = C 4 exp 1 q ′ C 3 m + (q ′ -1) α |g| X T + |g| X 1 -α -m = C 5 exp 1 q ′ C 3 m + (q ′ -1) α |g| X T 1 -α -m ,
where C 5 := C 4 exp( q ′ -1

q ′ α 1-α ). If we choose q ′ > 1 such that 7 A := C 3 -(q ′ -1)α/|g| X < 0, ∞ m=M E x 1 τ g K g (T ) ∈[m,m+1) sup t∈[0,T ] |X τ g K g (t) -X t | p ≤ C 5 exp (q ′ -1)αT (1 -α)|g| X ∞ m=M (e A ) m = C 5 exp (q ′ -1)αT (1 -α)|g| X e AM 1 -e A = C 5 exp (q ′ -1)αT (1 -α)|g| X exp C 3 -(q ′ -1)α |g| X M 1 -exp C 3 -(q ′ -1)α |g| X . If M is chosen such that M > T /(1 -α), then ∞ m=M E x 1 τ g K g (T ) ∈[m,m+1) sup t∈[0,T ] |X τ g K g (t) -X t | p ≤ C 5 e C 3 M 1 -e A exp (q ′ -1)α |g| X T 1 -α -M , ( 50
)
where the bound is O(e -1/|g| X ), and can be rewritten as C (1) e -C (2) /|g| X . From (47), ( 48) and (50),

sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) ≤ C 1/p(1-γp) 1 2|g| X 4T + |g| X γ/2 + C (1) e -C (2) /|g| X .
As |g| γ X and e -1/|g| X are both O(|g| γ/2 X ), there exists a constant C > 0, such that,

sup t∈[0,T ] |X τ g K g (t) -X t | L p (Px) ≤ C|g| γ/2 X ,
which is [START_REF] Chung | On the zero n 1 ±1[END_REF].

Computations for the classical SDE case and beyond

In order to implement Algorithm 1 for simulating paths of a diffusion X, one needs two things. First, a grid g adapted to the scale function and speed measure of X. Second, good approximations of the transition probabilities and conditional transition times of X over g. The first point was covered in Section 2.3.

In this section, we show how to compute the quantities (6) in the pure SDE case. The extension to SDE solutions with point-wise singularities is straightforward. This allows us via Algorithm 1 to simulate all such processes. Let (µ, σ) two real-valued functions that satisfy Condition 1 and X be the diffusion that solves

dX t = µ(X t ) dt + σ(X t ) dB t ,
where B is a standard Brownian motion. Let also I be the state-space of X. A straightforward computation using Itô's formula gives us the infinitesimal generator of X,

(L, dom(L)) = L f (x) = µ(x)f ′ (x) + 1 2 σ 2 (x)f ′′ (x), ∀f ∈ dom(L), dom(L) = f ∈ C b (I) : L f ∈ C(I) . (51)
In particular, if µ and σ are continuous, then, dom(L) = C 2 (I). From Proposition 3.1 and (51) we can infer that,

s ′ (x) = e -y a 2µ(ζ) σ 2 (ζ) dζ and m(dx) = 1 s ′ (x) 2 σ 2 (x) dx. (52)
Thus, from Proposition 3.1,

v 0 (x) =
From ( 26) and ( 52),

v 1 (x) = b a G a,b (x, ζ)v 0 (ζ) 1 s ′ (x) 2 σ 2 (x) dζ, v 1 (x) = b a G a,b (x, ζ) 1 -v 0 (ζ) 1 s ′ (x) 2 σ 2 (x) dζ. (54)
The scale functions and the speed measures are defined up to a multiplicative constant: for α ∈ R and λ > 0, the pairs (s, m) and (α + λs, λ -1 m) are associated to the same diffusion.

In particular, as we are only concerned with points x ∈ [a, b], we could use v 0 for the scale function. The speed measure shall be adapted accordingly. For x, ζ ∈ [a, b], the Green function in (1) takes the simpler form,

G a,b (x, ζ) = v 0 (x ∧ ζ) 1 -v 0 (x ∨ ζ) .
Thus, expressions (54) become

v 1 (x) = b a v 0 (x ∧ ζ) 1 -v 0 (x ∨ ζ) v 0 (ζ) v ′ 0 (ζ) 2 σ 2 (ζ) dζ, v 1 (x) = b a v 0 (x ∧ ζ) 1 -v 0 (x ∨ ζ) 1 -v 0 (ζ) v ′ 0 (ζ) 2 σ 2 (ζ) dζ. (55)
The quantities (53) and ( 54) or (55) can be computed analytically as we do for the Ornstein-Uhlenbeck process in Section 7.1.4 or approximated numerically as for the Cox-Ingersoll-Ross process in Section 7.1.4.

Skew Brownian motion

The skew Brownian motion at 0 with parameter β ∈ (0, 1) can be defined (see [8, p. 126]) as the diffusion process with scale function and speed measure

s(x) =    x β , x ≥ 0, x 1 -β , x ≤ 0 and m(dx) = 2β dx, x > 0, 2(1 -β) dx, x < 0.
From the definition of the scale function,

v 0 (x) = s(x) -s(a) s(b) -s(a) .
As the β and (1β) terms between the speed measure and the scale function compensate themselves in the expressions giving v k in Proposition 3.2,

v 1 (x) = v BM 1 (x) and v 1 (x) = v BM 1 (x).
In Figure 3-(b) are simulation results of the skew Brownian motion using Algorithm 1 and a plot of its probability transition kernel [8, p. 126].

Cox-Ingersoll-Ross process

The Cox-Ingersoll-Ross process or CIR process [START_REF] Cox | A theory of the term structure of interest rates[END_REF], introduced first by W. Feller [START_REF] Feller | Two singular diffusion problems[END_REF], is the diffusion that solves the SDE

dX t = θ(µ -X t ) dt + σ X t dB t , X 0 > 0,
where B t is a standard Brownian motion. The parameter θ > 0 expresses its mean reversion speed, µ its long term speed and σ > 0 is its diffusivity parameter. This equation has a diffusion coefficients that degenerates at 0. It however remains non-negative given X 0 > 0 and almost surely never hit 0 when 2θµ > σ 2 . A large body of work have been devoted to the simulation of the CIR and related process, see e.g. [START_REF] Alfonsi | On the discretization schemes for the CIR (and Bessel squared) processes[END_REF]. 19)-( 20) with Newton's method and (250, 200)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 46.5 sec). and

v 1 (x) = b a v 0 (x ∧ ζ) 1 -v 0 (x ∨ ζ) v 0 (ζ) v ′ 0 (ζ) 2 σ 2 ζ dζ, v 1 (x) = b a v 0 (x ∧ ζ) 1 -v 0 (x ∨ ζ) 1 -v 0 (ζ) v ′ 0 (ζ) 2 σ 2 ζ dζ.
As the scale function yields no satisfactory closed formula, we perform a numerical approximation of both v 0 and the couple (v 1 , v 1 ).

These functions may be computed numerically. One may choose a suitable grid when the process is close to 0, leading to improved convergence results. It is also noteworthy that the approximating process does not cross 0, a problem which arise when using Euler type schemes.

In Figures 4 and5 we plot the transition kernels and histograms of simulated values of CIR processes at T = 1 using Algorithm 1 with uniform and tuned grids and approximation precisions of the quantities (v 0 , v 1 , v1 ). We also superpose on the same figures the corresponding transition kernels. 19)-( 20) with Newton's method and (100, 100)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 12.8 seconds). (c): Same as Figure 5a but with a (250, 200)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 11.5 seconds). (d): Same as Figure 5b but with a (250, 200)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 12.2 seconds). 1.2, 0.8) with initial value x 0 = 0 using Algorithm 1 with a uniform grid of step-size h = 0.01. The quantities (56) were approximated using 100-step Riemann approximations.

Skew or reflected Bessel process

The Bessel process is defined as the solution of the SDE

dX t = δ -1 2X t dt + dB t ,
where B is a standard Brownian motion. If δ ∈ (0, 2), then 0 is an attainable boundary (see Section 2.4). In this case, we can then define the following behaviors at 0: reflection or skew/partial reflection. In particular, the Skew-Bessel process of dimension δ ∈ (0, 2) and skew β ∈ (0, 1) is the diffusion process with state-space R defined through s and m, where (see [START_REF] Alili | On the semi-group of a scaled skew Bessel process[END_REF])

s(x) =        1 β x 2-δ 2 -δ , x > 0, - 1 1 -β |x| 2-δ 2 -δ , x ≤ 0, m(dx) = 2βx δ-1 dx, x > 0, 2(1 -β)|x| δ-1 dx, x ≤ 0.
This yield the following expressions for the quantities we need to compute in order to orders of convergence of localized path-sensitive functionals. An example is the stickiness parameter estimation in [START_REF] Anagnostakis | Functional convergence to the local time of a sticky diffusion[END_REF]. It is shown there that, if X is a sticky Brownian motion of stickiness ρ > 0, α ∈ (0, 1/2) and g is an integrable function vanishing in the vicinity of 0, the statistics

L n [g, α] := n α n [nt] i=1 g(n α X i-1 n ) (57)
and ρ n [g, α] := 2 λ(g)

L n [g, α]

1 n

[nt] i=1

1 X i-1 n =0 (58)
converge in probability to L 0 t (X) and ρ respectively. This qualifies ρ n [g, α] as a consistent estimator of the stickiness parameter.

It is also observed that, as long as n is high enough:

• The convergence of (58) seems to hold for any α ∈ (0, 1).

• The convergence speed of (58) seems to increase in terms of α ∈ (0, 1).

• The higher the α, the less observations of the trajectories are observed via g.

Algorithm 1 gives us the flexibility to remediate to the latter by using grids of higher precision around the point of stickiness. Thus, the statistics (57) and (58) are relevant for larger values of α and we achieve higher orders of convergence without significant increase in the numerical complexity.

Let g 0 and g 1 be two grids defined by: g 0 (h) = {0} ∪ ± (h 2 /ρ + kh); k ∈ N , g 1 (h) = ± x k (h); k ≥ 0 ,

where {x j (h)} j≥0 is defined recursively by

x 0 = 0, x j = x j-1 + h 2 ρ 1 x j-1 + 1 + h 1 - 1 x j-1 + 1 1 x j-1 <1 + h1 x j-1 ≥1 . (60)
We observe that

|g 1 | = |g 0 | and |g 1 | X = |g 0 | X . (61)
Thus, from Theorem 2.1, the rate of convergence of the STMCA is O(h 1/2 ) for both grids.

Simulation results

In this section we present Monte Carlo simulation results. The integer N MC will be the Monte Carlo simulation size. For every j ≤ N MC we simulate the path of an approximation process X j t . To assess the quality of each Monte Carlo estimation we use the following metrics:

Using g 1 we have an abundance of simulated path-wise observations close to the point of stickiness. The statistic (58) remains thus relevant for large values of α and we can achieve higher orders of convergence (see Figures 89, and Tables 12). It is worth noting that simulation with h = 0.001

A Conditioning on the embedded path

Lemma A.1. Let X be a diffusion process with state-space I an interval of R defined on a family of probability spaces P = (Ω, (F t ) t≥0 , P x ) x∈I such that for every x ∈ I, P x (X 0 = x) = 1.

Let also g be a covering grid of I, {τ j } j∈Z + the sequence of embedding times of X in g defined in [START_REF] Hajri | Application of stochastic flows to the sticky Brownian motion equation[END_REF] and B the sigma-algebra defined by B = σ{X τ j ; j ∈ N 0 }. Then, for any measurable path-functional F : C 0 (R + , I) → R and j ≥ 1 and x ∈ I, E x F ((X τ j-1 ,τ j t

) t≥0 ) B = E x F ((X

τ j-1 ,τ j t ) t≥0 ) X τ j-1 , X τ j , (64)
where X τ j-1 ,τ j is the process defined for every t ≥ 0 by X τ j-1 ,τ j t = X (τ j-1 +t)∧τ j .

Proof. Let us fix x 1 , x 2 , . . . be a sequence of points in the grid. Let us define Q(x; x 1 , x 2 , . . . ) := P x X τ 1 = x 1 , X τ 2 = x 2 , . . .

By the strong Markov property, P x X τ i+1 = x i+1 , X τ i+2 = x i+2 , . . . F τ i = P Xτ i X τ 1 = x i+1 , , X τ i+2 = x i+2 = Q(X τ i ; x i+1 , x i+2 , . . . ).

Using the strong Markov property twice, first conditioning first with respect to F τ j+1 and then with respect to F τ j , E x 1 Xτ 1 =x 1 ,...,Xτ j =x j F ((X τ j-1 ,τ j t ) t≥0 )1 Xτ j =x j ,Xτ j+1 =x j+1 ,... = E x 1 Xτ 1 =x 1 ,...,Xτ j =x j F ((X τ j-1 ,τ j t ) t≥0 )1 Xτ j =x j Q(X τ j ) = E x 1 Xτ 1 =x 1 ,...,Xτ j-1 =x j-1 (X τ j-1 ,τ j t ) t≥0 Q(x j ) = E x 1 Xτ 1 =x 1 ,...,Xτ j =x j-1 R(x j-1 , x j )Q(x j )

with R(x, y) := P x F ((X 0,τ 1 t ) t≥0 )1 Xτ 1 =y . Therefore, using the definition of the conditional expectation, E x F ((X τ j-1 ,τ j t ) t≥0 ) X τ 1 = x 1 , X τ 2 = x 2 , . . .

=

R(x j-1 , x j ) P x j-1 X τ 1 = x j = E x j-1 F ((X 0,τ 1 t ) t≥0 X τ 1 = x j .

This is sufficient to prove the result.

P

  x (τ b < τ a ) = s(x)s(a) s(b)s(a) , for a < x < b ∈ I and with τ a being the hitting time of a. The measure m, called the speed measure, is the unique positive locally finite measure such that, E x (τ ab ) = (a,b) G a,b (x, y)m(dy), for a < x < b ∈ I and where G a,b (x, y) is the Green function defined as G a,b (x, y) )s(a) s(b)s(y) s(b)s(a) , for x ≤ y, s(y)s(a) s(b)s(x) s(b)s(a) , for x > y.

1 2

 1 -ϵ ) for the adapted grid. Numerical examples are given in Section 7.1.2.

Figure 2 :

 2 Figure 2: The covering grid g = Z + of [0, ∞) along with its first cells C(g).

Proposition 3 . 3 .Lemma 3 . 4 .

 3334 The function v k (x) = E x (τ k ab 1 τ b <τa ) is solution to the problem with Dirichlet boundary conditions      L u = -kv k-1 , x ∈ (a, b), u(a) = 0, u(b) = 0. Let g(x) = (a,b) G a,b (x, y)f (y)m(dy), where f ∈ C 0 b (a, b) and G a,b (x, y) is the Green function defined in (1). Then g ∈ dom(L) and L g(x) = -f (x), ∀x ∈ (a, b). Proof. Let x ∈ (a, b). Using the D m D s factorization of L and the dominated convergence theorem we get L g(x) = D m D s (a,b)

Figure 3 :

 3 Figure 3: (a): Histogram of simulated values at T = 1 of a sticky Brownian motion of parameter ρ = 0.7 with initial value x 0 = 0 using Algorithm 1 with the tuned grid (14) of size-criteria h = 0.01. (b): Histogram of simulated values at T = 1 of a skew Brownian motion of parameter β = 0.9 with initial value x 0 = 0 using Algorithm 1 with a uniform grid of step-size h = 0.01.

From

Figure 4 :

 4 Figure 4: Histogram of simulated values at T = 1 of a CIR process of (θ, µ, σ) = (5, 5, 1) with initial value x 0 = 1 using Algorithm 1 with: (a): a uniform grid of step-size h = 0.01 and (250, 200)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 44.5 sec). (b): a tuned grid of size-criteria h = 0.01 computed solving numerically (19)-(20) with Newton's method and (250, 200)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 46.5 sec).

Figure 5 :

 5 Figure 5: Histogram of simulated values at T = 1 of a CIR process of (θ, µ, σ) = (1, 2, 1) with initial value x 0 = using Algorithm 1 with: (a): a uniform grid of step-size h = 0.015 and (100, 100)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 10.8 seconds). (b): a tuned grid of size-criteria h = 0.015 computed solving numerically (19)-(20) with Newton's method and (100, 100)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 12.8 seconds). (c): Same as Figure5abut with a (250, 200)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 11.5 seconds). (d): Same as Figure5bbut with a (250, 200)-step Riemann approximation of (v 0 , v 1 ) (simulation time: 12.2 seconds).

Figure 6 :

 6 Figure 6: Histogram of simulated values at T = 1 of the Skew Bessel process of parameters (δ, β) = (1.2, 0.8) with initial value x 0 = 0 using Algorithm 1 with a uniform grid of step-size h = 0.01. The quantities (56) were approximated using 100-step Riemann approximations.

  Algorithm 1 Space-Time Markov Chain Approximation (STMCA) Algorithm Input: x initial value, T time horizon, g = {x j } j∈J grid on I Output: ( X[k]) k consecutive values taken by the approximation process, ( t[k]) k consecutive transition times

up to an affine transformation

Which means that s(x) = x.

This stands for "continue à droite avec une limite à gauche", that is right-continuous with left-limit.

This results in the bound of Theorem 2.1 not depending on the starting point of the diffusion.

The right-inverse of a function f is given by,f -1 (x) = inf{ζ ≥ 0 : f (ζ) > x}.

We remark that if this is satisfied for a grid g, then it is satisfied for all grids g ′ such that |g ′ | X ≤ |g| X .

Acknowledgments The PhD thesis of A. Anagnostakis is supported by a scholarship from the Grand-Est Region (France).

The standard Brownian motion can be defined as (see [8, p. 119]) the diffusion process with scale function and speed measure From Proposition 3.2,

, where by BM we mean these quantities are associated with the Brownian motion. Thus, from [START_REF] Ankirchner | A functional limit theorem for coin tossing Markov chains[END_REF], we have all the necessary quantities we need to implement the algorithm.

Sticky Brownian motion

The Brownian motion with a sticky point at 0 of stickiness ρ > 0, already presented in Example 2.4, can be defined [8, p. 123] as the diffusion process with scale function and speed measure s(x) = x and m(dx) = 2 dx + ρδ 0 (dx).

From the definition of the scale function,

From Proposition 3.2 we may deduce the following expressions for the conditional exit times of c x = (a, b),

are the analogous quantities for the standard Brownian motion.

In Figure 3-(a) are simulation results of the sticky Brownian motion using algorithm 1 and its probability transition kernel (see [8, p. 124]).

where G a,b (x, ζ) is defined in [START_REF] Alfonsi | On the discretization schemes for the CIR (and Bessel squared) processes[END_REF].

In Figures 6 and7 we plot the transition kernels and histograms of simulated values of: a skew and two reflected Bessel processes. The plotted corresponding transition kernels are computed in [START_REF] Alili | On the semi-group of a scaled skew Bessel process[END_REF].

Local time approximation

The following example illustrates the flexibility of Space-Time Markov Chain Approximation (STMCA) generated via Algorithm 1. One feature of such approximation processes is that they are defined on a given grid. With a suitable choice of grid it is possible to achieve higher Table 1: Stickiness parameter estimations using the grid g 0 for h = 0.01. The missing values in the table corresponds to cases where the statistic (57) is observed to be 0. Computation time (single-core): 4 seconds.

• ρ MC , S 2 MC , σ MC : Monte Carlo estimation, variance and standard deviation of the stickiness parameter estimations,

• acc: average number of path-values observed by g, i.e

We use the test function g(x) = 1 1<|x|<5 /8 which satisfies the conditions of Corollary 1.3 of [START_REF] Anagnostakis | Functional convergence to the local time of a sticky diffusion[END_REF]. Within each table (Tables 12) we use the same simulated STMCA trajectories of the sticky Brownian motion of parameter ρ = 1.

We observe that the usage of grid g 1 yields far superior results than g 0 (see Tables 12).


Belal Hmedan

Dorilys Kilgus

Humbert Fiorino

Aurélie Landry

Damien Pellier

Adapting Cobot Behavior to Human Task Ordering Variability for Assembly Tasks

In this paper, we address flexible assembly plans generation to accommodate human task variability. Unlike existing approaches, the proposed approach is based on a free-style human-robot interaction (HRI) that does not impose any task order on the participants and furthermore can accommodate their errors and help them to correct the errors without stopping the whole assembly process. Our approach is implemented in a real robotic architecture that combines sensory-motor modules with Hierarchical Task Networks (HTN) to endow the cobot with the necessary adaptability to correspond to human actions dynamically. We show experimentally with 56 participants on a simulated industrial assembly task that the cobot increases the task performance (reduction in the number of errors and of gestures) without increasing the participants' cognitive load.

Introduction

Manufacturing assembly automation is a permanent challenge for the industry as innovation is a real lever in order to keep production performance in a worldwide economic market. However, this automation has obliged the skilled operators, who cannot be completely substituted by robots, to follow specific operation paradigms to cope with the robots' operation procedures. This type of constrained operation increases the cognitive load of the operators [START_REF] Kalakoski | Effects of a cognitive ergonomics workplace intervention (cogerg) on cognitive strain and well-being: a cluster-randomized controlled trial. a study protocol[END_REF][START_REF] Karasek | Job demands, job decision latitude, and mental strain: Implications for job redesign[END_REF], reduces the room for maneuver in their task performance and causes musculoskeletal disorders [START_REF] Caroly | Innovation in the occupational health physician profession requires the development of a work collective to improve the efficiency of msd prevention. 8th World congress on Ergonomics -Designing a sustainable future[END_REF][START_REF] Srinivasan | Motor variability in occupational health and performance[END_REF]. This is due to a lack of human variability consideration in the design of the automation process. It is therefore important to take into account human variability, e.g., level of experience or fatigue, size, etc. very early in the robotic architecture design, and to endow the robot with the ability to adapt to human beings in order to prevent these risks.

Among the human variabilities that must be taken into account, human task variability is one of the most important. It is a question of being able for a robot to adapt the order in which it performs the tasks to the order chosen by a human and not the other way around. The aim is to give the human a great deal of freedom in carrying out a collaborative task with a robot and not to constrain him/her more than necessary and to alleviate his/her cognitive load by not imposing specific task orders, leaving the human free to delegate more or fewer tasks to the robot.

Accommodating robot's behaviors to human task variability is a challenging issue due to the uncertainty in human behavior [START_REF] Faroni | A layered control approach to human-aware task and motion planning for human-robot collaboration[END_REF] and the probability of conflicts between human and cobot actions [START_REF] Cheng | Humanaware robot task planning based on a hierarchical task model[END_REF]. To deal with this issue, we propose in this paper a robotic architecture bridging sensory and action elements with task planning reasoning capabilities [START_REF] Pellier | PDDL4J: a planning domain description library for java[END_REF] implemented on the cobot YuMi to ensure a smooth workflow in collaborative tasks where humans and cobots work together in close vicinity and share the same environment. Unlike other approaches e.g., [START_REF] Faroni | A layered control approach to human-aware task and motion planning for human-robot collaboration[END_REF][START_REF] Gabler | A game-theoretic approach for adaptive action selection in close proximity human-robot-collaboration[END_REF][START_REF] Paxton | Evaluating methods for end-user creation of robot task plans[END_REF][START_REF] Roncone | Transparent role assignment and task allocation in human robot collaboration[END_REF], we do not impose any order on the operator. Moreover, the robot and the operator can simultaneously act in the same area to accomplish a shared assembly task.

Our contribution in this paper is twofold: (1) a free-style human-robot interaction (HRI) that does not impose the task-order on the participants and furthermore can accommodate their errors and help them to correct the errors without stopping the whole assembly process, and (2) an extensive evaluation methodology that highlights the benefits of integrating a cobot into the assembly process using not only a classical robotics metric to evaluate task performance but also metrics from ergonomics. This methodology handles the problem from different aspects, including the subjective cognitive load measured through two phases with and without the cobot, of experiments performed with 56 participants.

The rest of the paper is organized as follows: Section II presents the related work; section III presents the architecture proposed; Finally, section IV introduces the experimental setup and the results showing that our approach reduces the cognitive load of the operators and increases the task performance.

RELATED WORK

Robotic architecture bridging sensory and action elements with task planning reasoning capabilities in the context of industrial human-robot collaboration was discussed mainly from the performance point of view, where the time of planning and scheduling is the main criteria for evaluation [START_REF] Faroni | A layered control approach to human-aware task and motion planning for human-robot collaboration[END_REF][START_REF] Hawkins | Anticipating human actions for collaboration in the presence of task and sensor uncertainty[END_REF][START_REF] Johannsmeier | A hierarchical human-robot interaction-planning framework for task allocation in collaborative industrial assembly processes[END_REF][START_REF] Zhu | Hierarchical planning for long-horizon manipulation with geometric and symbolic scene graphs[END_REF]. In contrast, human aspects such as acceptability, usability and cognitive load were less discussed [START_REF] Cheng | Humanaware robot task planning based on a hierarchical task model[END_REF][START_REF] Gabler | A game-theoretic approach for adaptive action selection in close proximity human-robot-collaboration[END_REF][START_REF] Roncone | Transparent role assignment and task allocation in human robot collaboration[END_REF]. Furthermore, most works limit the freedom of ac-Figure 1: The proposed modular architecture tion of humans in the execution of their task, for example by assigning specific tasks to humans [START_REF] Gabler | A game-theoretic approach for adaptive action selection in close proximity human-robot-collaboration[END_REF][START_REF] Hawkins | Anticipating human actions for collaboration in the presence of task and sensor uncertainty[END_REF][START_REF] Roncone | Transparent role assignment and task allocation in human robot collaboration[END_REF], or limiting their choices in terms of assembly zones and assembled block colors [START_REF] Faroni | A layered control approach to human-aware task and motion planning for human-robot collaboration[END_REF].

To relax these constraints it is necessary to grant the human a higher degree of liberty to make the interaction with the cobot smooth, natural, and intuitive. Hence, it is better not to impose task order on the human, but instead, adapt the task plans with respect to the human actions. A number of works have tried to deal with the problem of adapting a robot's behavior to the order of execution of a task by a human and many techniques have been proposed to generate behaviors more or less dynamically. These are the techniques (in non-exhaustive list): Planning and Scheduling [START_REF] Faroni | A layered control approach to human-aware task and motion planning for human-robot collaboration[END_REF][START_REF] Le | Hierarchical human-motion prediction and logic-geometric programming for minimal interference human-robot tasks[END_REF][START_REF] Lee | Adaptive task planner for performing home service tasks in cooperation with a human[END_REF][START_REF] Pan | A general task and motion planning framework for multiple manipulators[END_REF], AND-OR Graphs [START_REF] Hawkins | Anticipating human actions for collaboration in the presence of task and sensor uncertainty[END_REF], Behavior-Tree [START_REF] Paxton | Evaluating methods for end-user creation of robot task plans[END_REF][START_REF] Rovida | Extended behavior trees for quick definition of flexible robotic tasks[END_REF], Hierarchical Task Networks (HTN) [START_REF] Cheng | Humanaware robot task planning based on a hierarchical task model[END_REF][START_REF] Grigore | Preference-based assistance prediction for human-robot collaboration tasks[END_REF][START_REF] Kuhner | Closed-loop robot task planning based on referring expressions[END_REF], Game Theory [START_REF] Gabler | A game-theoretic approach for adaptive action selection in close proximity human-robot-collaboration[END_REF], Markov Decision Process (MDP), or a combination of them [START_REF] Roncone | Transparent role assignment and task allocation in human robot collaboration[END_REF]. Among these techniques, task planning and more precisely hierarchical planning is widely used and very effective to generate symbolic plans in a similar way of the human decomposition of complex tasks into simpler actions [START_REF] Erol | Htn planning: Complexity and expressivity[END_REF].

Subsequently, the challenge of the tasks order adaption is tackled in this paper with a real industrial cobot YuMi to best imitate an industrial assembly task and not only on simulation [START_REF] Faroni | A layered control approach to human-aware task and motion planning for human-robot collaboration[END_REF][START_REF] Le | Hierarchical human-motion prediction and logic-geometric programming for minimal interference human-robot tasks[END_REF]. The robot's behavior is made explicit to the human, through an interface so as to consolidate his/her understanding and his/her acceptability of its behavior [START_REF] Milliez | Using human knowledge awareness to adapt collaborative plan generation, explanation and monitoring[END_REF]. An explanation of the cobot next action, and human errors are provided for error correction while the assembly keeps uninterrupted.

Adaptive Robotic Architecture

Our approach in tackling human task variability stands on a modular architecture consisting of four modules (Fig 1 ): 1. Perception module to perceive the workplace (including the human task variability) and encode the perceived information into a discrete representation.

2.

Planning module to generate the adaptive behavior that corresponds to the human task variability.

3.

Execution module to execute the intelligent behavior on the cobot.

4.

Controller module to manage the interaction between the operator and the cobot, and explain the intelligent behav-

PLANNING MODULE

To adapt the cobot behavior to the human task variability we take advantage of the flexibility and the intelligence that Automated Planning (AP) can generate. Specifically, Hierarchical Task Networks (HTN) technique allows to compute action sequences i.e. plans with respect to given world states, and to provide explanations on the causal links of these actions at different levels of abstractions [START_REF] Milliez | Using human knowledge awareness to adapt collaborative plan generation, explanation and monitoring[END_REF].

For this purpose we have developed a conceptual model composed of a domain, and a problem in the Hierarchical Domain Definition Language (HDDL) [START_REF] Höller | HDDL: an extension to PDDL for expressing hierarchical planning problems[END_REF]. The planner used in our architecture is the Totally Ordered Fast Downward (TFD) Planner [START_REF] Pellier | Totally and partially ordered hierarchical planners in PDDL4J library[END_REF]. The domain in our model is static and encompasses the objects in our scenario and the relations (primitive actions) between them, as well as the higher-level tasks, and different methods to decompose the tasks into primitive actions according to triggering preconditions. On the other hand, the problem in our case is dynamic in the sense that the initial state is obtained from the discrete representation generated by the vision system, which varies along the cobot actions, and the operator variable actions (Fig 3). Re-planning is done after each perception phase sequentially for each block until a plan is found for a given block, otherwise if there is no plan found for all blocks, hence planning stops and the user is informed to do the task alone.

The predicates in HDDL define the status of each object, for instance the gripper position, orientation, fingers status (open/close), and whether it is loaded or not. Furthermore, the points are described in terms of occupancy and neighboring points. Finally, blocks orientation, gripping position, and whether the block is attached to the workspace or not. The tasks in our domain are limited to Pick, and Place, and the parameters are the gripper which will be used, the block to be picked/placed, and the direction of the placement:

(:task pick :parameters (?g -Gripper ?x -Block ?d -Direction)) (:task place :parameters (?g -Gripper ?x -Block ?pos -Point_Workspace))

These tasks are combined with the necessary primitive actions, and the hierarchical decomposition methods to provide a flexible way to execute the same task in various available approaches. An example of the actions done here is holding a block with the gripper, this action changes the status of the cobot fingers from open to close, and the gripper becomes loaded with the block:

(:action Hold_gripper :parameters (?x -Block ?g -Gripper) :precondition (and (is_open ?g) (unloaded ?g)) :effect (and (loaded ?g ?x)

(not (is_open ?g)) (not(unloaded ?g))))

The method to decompose the pick task in the case of 2x2 block will have the selected gripper, the block to be picked, the direction of the placement, and the waypoint where the gripper will move to after picking the block as parameters. This method decomposes the pick task into a series of subtasks/actions in the case that all the pre-conditions are satisfied:

(:method pick_2x2 :parameters(?g -Gripper ?x -Block_2x2 ?d -Direction ?pa -Point_Air) :task (pick ?g ?x ?d) :precondition (and (not (left_direction ?d))

(not (right_direction ?d)) :ordered-subtasks (and (Open_gripper ?g) (Rotate_empty_gripper_V ?g) (Move_open_gripper_V_To_Stock ?g ?pa)

EXECUTION MODULE

The execution of the adapted behavior represented in the generated plan is done in our scenario on the cobot YuMi (for the sake of simplicity we used one arm only). Pick/Place locations are passed from the plan to the low-level language (RAPID) to execute the corresponding motion. We have developed a special Application Programming Interface (API) to use the Robot Web Services (RWS) platform provided by the robot manufacturer in order to execute the plan on the cobot, and get the execution status (Fig 4).

USER INTERFACE MODULE

It is important to express the cobot behavior to the operator to ensure workflow smoothness in the collaborative joint task, hence a user interface is used to explain: i) Assembly progress. ii) Cobot's next movement, to ensure human awareness of the cobot behavior in order to avoid conflicts. iii) Operator mistakes, to make the cobot adaptation clear and understandable for the operator and to increase the operator acceptance. Therefore, we have developed a user interface to demonstrate the cobot behavior, where the elements inside this interface emulate the situation of the workplace. It contains three distinctive areas (Fig 5 ): 1. Assembly pattern selection and control area. 2. Assembly and swap zone area. 3. Messages area.

The first area is used to control the experiments, while the rest are interactive areas, area (2) visualises the location and the color signals, while area (3) is used to display the interactive messages.

The color signals used in the graphical interface are as follows: I) Green-colored squares to represent the green base-plate to indicate unoccupied positions. II) Red-blinking squares to represent the positions that the cobot will use to place a block. III) Yellow/Blue colored squares to represent the occupied positions. IV) Red cross on top of the Yellow/Blue squares to represent any incorrect placement.

Moreover, the interactive messages are: i) "I'm taking a photo to understand the situation." ii) "I'm placing a block in the blinking area." iii) "I can't help you anymore, I'll let you finish." iv) "The task is solved, thank you very much"

Experiments

The objective of this experimentation is to evaluate the impacts of integrating intelligent cobot equipped with an adaptive action-selection mechanism and free-style user interaction on runtime and operators' well-being. The former (runtime) represents productivity, while the latter (well-being) is interpreted as the physical and cognitive demands of the task. Thus, we consider the performance with an extended vision that takes into account not only runtime but also human ergonomics.

The experiments were performed in two phases, the first with 33 participants, and the second with 23 participants different from the first phase. They all were informed about the experiments and signed a consent form beforehand.

Experimental Scenario

The experiments are meant to simulate a collaborative assembly task in an industrial environment, hence LEGO blocks were used to construct four patterns ( In this workspace divided into four spaces (Fig 2 ), the participant is supposed to use the stock from his/her dedicated area, and the swap area, while the cobot is programmed to use the stock from its dedicated area. Initially, the swap and the assembly areas are empty, while the stock of both the participant and the cobot is filled with blocks of two different sizes 2 × 2 and 2 × 4, and two different colors: yellow, and blue. The user interface is used to start the experiment, where the cobot starts by perceiving the environment, and the user is informed through a message displayed on the user interface "I'm taking a photo to understand the situation.", the perceived environment is interpreted into a discrete representation exploited by our conceptual model to generate the adaptive response. After that the cobot executes the sub-task (single block assembly) and the participant is informed of the cobot action and place of execution by the user interface: the message "I'm placing a block in the blinking area." appears on the user interface, and the red-blinking squares, which reflects the positions that the cobot needs to do the placement, appear on the screen (Fig 5). The user at this stage is supposed to do the assembly at the same time as the cobot, with respect to the red-blinking positions that the cobot will use. This loop of perception, planning, collaborative action is repeated until the cobot cannot place any block because of the size of the robot's fingers. Then the user is asked to complete the task alone in the user interface: "I can't help you anymore, I'll let you finish.", and this message is presented until the perception module verifies the task completion, thus the last message appears: "The task is solved, thank you very much". The assembly task complexity stems from the fact that the task is changing according to the selected pattern, thus the size and the color of the blocks vary from one pattern to another making the task cognitively demanding.

Moreover, this scenario reflects actual collaborative assembly tasks in manufacturing where the operator and the cobot execute teamwork. The code used is released online1 for experiment reproducibility.

Evaluation Criteria

Four factors were used for evaluation in our case: time of the task completion, to evaluate the impact on the productivity, errors committed during the task (any misplacement is considered an error, even if it was corrected later), to evaluate the impact on the quality, count of gestures per task (each movement of the participant's hand to do a sub-task is counted as a gesture, even if he did not place a block), which in turn reflects the physical demands of the task, and cognitive load measured using the NASA-TLX [START_REF] Hart | Development of nasa-tlx (task load index): Results of empirical and theoretical research[END_REF] form.

Methodology

The experiments were conducted in two separate phases with different groups of participants in each phase: in the first phase, one person fulfills the block assembly task, this constitutes our baseline. In the second phase, the assembly task is jointly realised by the cobot and a human operator. It is noteworthy that the experiments focus on the Human-Robot adaptation in joint task completion, not manipulation skills. From this perspective, comparable setups can be met in the industry without significant changes because the assembly process is similar, though the size and the shape of the assembled parts vary with the considered application.

First Phase: Human Alone

The first phase of the experiments was performed with humans only (without the cobot), in order to compare the results with the second phase. This phase of the experiments consists of four steps:

1. Introduction of the tasks to be carried out: a short introduction of the tasks and the framework is given to familiarize the participant with the experiment. The setup is as in Fig. 2, except that the cobot, the camera, and the user interface are not present in the first phase. 2. A demonstration of assembly pattern: a demonstration of "Model D" assembly is carried out to eliminate any ambiguity, then the participant is asked if he has any further questions about the task assembly to be clarified. 3. The assembly tasks: The user is asked to assemble the other three models (A, B and C) according to a given document in a specific order 2 one by one, and the experiments are video recorded to be analyzed later. 4. The survey: Finally, the participant is asked to answer the NASA-TLX test to measure his/her cognitive load.

Second Phase: Human and Cobot

The second phase involves the cobot with the human to achieve a shared task. Furthermore, this phase has three additional components: the cobot, the camera to capture the world, and the user interface projected on the external screen beside the cobot. The second phase of the experiment is executed in four steps as follows:

1. Introduction of the tasks to be carried out: a short introduction of the cobot, the tasks, and the framework is given to familiarize the participant with the experiment. The setup that the participant is introduced to consists of:

• The user interface (Fig 5 ) including the signals and the messages that the participant will encounter during the experiment. • The camera, and the captured workspace. And the workplace: (a) The cobot abilities and the high safety constraints to make the participant feels confident that there is no harm in working in collaboration with the cobot. (b) The pattern which has to be assembled in collaborative manner with the cobot. (c) The stock to be used by the participant. (d) The stock of the cobot, which the participant should not use. (e) The assembly zone where the tasks are carried out. (f) The swap zone, where the cobot handles a block to the participant to assist him in doing the task in case the user is out of stock of a specific size/color. 2 The sequence introducing the patterns is randomized to avoid any bias in the results.

2. A demonstration of assembly pattern: a demonstration of "Model D" assembly is carried out with assistance of the cobot and the user interface to eliminate any ambiguity or panic, then the participant is asked if he/she has any further questions about the task assembly, the cobot behavior, or the user interface to be clarified.

The assembly tasks:

The user is asked to assemble the other three models (A, B and C) in a specific sequence 2 one by one with assistance of the cobot, and the experiments are video recorded to be analyzed later. Note that at this phase, the start/finish signals of the tasks were replaced by the messages emitted from the user interface, and not expressed by the participant as in the first phase.

4. The survey: Finally, the participant is asked to answer the NASA-TLX test to measure his/her cognitive load.

EXPERIMENTAL RESULTS

The preliminary results (Table 1) indicate a significant decrease in the average number of errors (60%), the average number of gestures (31%), and the average subjective cognitive load (2%). On the contrary, the time for task completion has increased (150%) due to the simplification 3

Discussion

In this work, we proposed a modular architecture to adapt cobot behavior to human task variability. This architecture makes use of HTN to generate flexible plans that account for human task variability, allows free-style human-robot interaction, and assists the human partner in correcting their errors without hindering the progress of the assembly process. These preliminary results show that, from an ergonomic perspective, this architecture and the automated planning are adequate to address human-cobot interactions effectively in industrial workplaces. The intelligent behavior acceptability is enhanced by a user interface that makes explicit the cobot behavior to the user. The errors decrease significantly indicating higher quality in performing the task, the number of gestures also decreases, which reflects less physical load, and the cognitive load has significantly decreased.

Task completion time increase highlights the importance of utilizing an enhanced module for the motion and the coordination between the robot and the human partner. We do believe there is room to enhance this architecture's runtime performance and converge towards human-alone effectiveness.

Figure 2 :

 2 Figure 2: Workplace zones

Figure 3 :

 3 Figure 3: Planning module

Figure 4 :

 4 Figure 4: Execution module

Figure 5 :

 5 Figure 5: Graphical interface module

  Fig 6) each of them is 4 × 4 blocks of two different colors yellow/blue (Fig 6), pattern (D) is used for demonstration purpose, while the rest are used for the experiments with the participants. All the experiments have been designed by ergonomists (coauthors) expert in industrial workplaces to specifically investigate free-style human-robot interaction in realistic situations.

Figure 6 :

 6 Figure 6: Assembly patterns

Table 1 :

 1 in the framework, hence finer optimization has to be considered in the further studies to remediate this drawback. Experimental Results. SD -Standard Deviation, Avg -Average value. Phase 1: 33 participants. Phase 2: 23 participants different from Phase 1.

			Phase 1
		Min	Max	SD	Avg
	Time (minutes)	1.1	4.4	0.92	2.73
	Errors	0	16	4	5
	Gestures	41	124	20	74
	Cognitive load	15.2% 59.07% 12.91% 37.49%
			Phase 2
		Min	Max	SD	Avg
	Time (minutes)	4	10.68	1.55	6.82
	Errors	0	6	1.83	2
	Gestures	34	72	10	51
	Cognitive load 11.33% 62.53% 15.32% 35.52%

https://github.com/bhomaidan1990/Legos

Only left arm of the cobot was used, and the speed was reduced for the sake of comfort human-robot interaction.