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

Alexis Anagnostakis, Antoine Lejay, Denis Villemonais

► To cite this version:

Alexis Anagnostakis, Antoine Lejay, Denis Villemonais. General diffusion processes as the limit of time-space Markov chains. The Annals of Applied Probability, 2023, 33 (5), pp.3620-3651. 10.48550/arXiv.2206.03713 . hal-02897819v2

HAL Id: hal-02897819

https://inria.hal.science/hal-02897819v2

Submitted on 8 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# General diffusion processes as the limit of time-space Markov chains

Alexis Anagnostakis*¹, Antoine Lejay†¹ and Denis Villemonais‡¹

¹Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France

## Abstract

We prove the convergence of the law of grid-valued random walks, which can be seen as time-space Markov chains, to the law of a general diffusion process. This includes processes with sticky features, reflecting or absorbing boundaries and skew behavior. We prove that the convergence occurs at any rate strictly inferior to $(1/4) \wedge (1/p)$ in terms of the maximum cell size of the grid, for any $p$-Wasserstein distance. We also show that it is possible to achieve any rate strictly inferior to $(1/2) \wedge (2/p)$ if the grid is adapted to the speed measure of the diffusion, which is optimal for $p \le 4$. This result allows us to set up asymptotically optimal approximation schemes for general diffusion processes. Last, we experiment numerically on diffusions that exhibit various features.

**Subject classification:** 60J60, 60F17

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

## 1 Introduction

In the diffusion process literature, the most well-studied and straightforward way to approximate diffusion processes is the Euler scheme. While such approximations work well for non-degenerate stochastic differential equations, this is not the case for more general diffusion

<footnote>* alexis.anagnostakis@univ-lorraine.fr</footnote>

<footnote>† antoine.lejay@univ-lorraine.fr</footnote>

<footnote>‡ denis.villemonais@university.edu</footnote>
processes [22]. The Euler scheme is also not well-defined for processes that exhibit sticky features, skew behavior [20] or slowly reflecting boundaries.

First studied by Feller [17, 18], general diffusion processes are used to replicate a wide range of phenomena like semi-permeable layers [24], principal agent problem dynamics [28, 32] and interest rates behavior close to 0 [27]. For more theoretical results and applications of general diffusions, we refer the reader to [7, 13, 14, 21, 25]. 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 [3], 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 [27] and [26], 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 [19], where the authors use a Continuous Time Markov Chain to identify events in genomics evolution. In [5, 6], 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 [29, 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 [25] 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 [6] and [15]. 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) \wedge (2/p)$ in terms of the maximum cell size for all $p$-Wasserstein distances. This convergence rate is optimal for $p \le 4$ according the Donsker invariance principle, as this is the rate simple randoms walk converges to the standard Brownian motion [12].

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 [3, 12, 15] are all special cases of it.

We make heavy use of the speed measure/scale function characterization of diffusions: from Proposition 3.14 of [29, Chapter VII], the law of a one-dimensional diffusion process with state-space $\mathbb{I}$ an open interval of $\mathbb{R}$ is entirely determined by an increasing continuous function $s$ and a positive locally finite measure $m$ defined on $\mathbb{I}$. The function $s$ called the scale function, is defined as the unique<fnref n="1" /> function such that,

$$P_x(\tau_b < \tau_a) = \frac{s(x) - s(a)}{s(b) - s(a)},$$

for $a < x < b \in \mathbb{I}$ and with $\tau_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(\tau_{ab}) = \int_{(a,b)} G_{a,b}(x, y)m(dy),$$

for $a < x < b \in \mathbb{I}$ and where $G_{a,b}(x, y)$ is the Green function defined as

$$G_{a,b}(x, y) = \begin{cases} \frac{(s(x) - s(a))(s(b) - s(y))}{s(b) - s(a)}, & \text{for } x \le y, \\ \frac{(s(y) - s(a))(s(b) - s(x))}{s(b) - s(a)}, & \text{for } x > y. \end{cases} \quad (1)$$

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

$$L = D_m D_s, \\ \mathrm{dom}(L) = \{f \in C_b^0(\mathbb{I}, \mathbb{R}) : L f \in C_b^0(\mathbb{I}, \mathbb{R})\},$$

where $D_m$ and $D_s$ are defined by

$$D_m f(x) = \lim_{h \to 0; h > 0} \frac{f(x+h) - f(x)}{m(x, x+h)}, \\ D_s f(x) = \lim_{h \to 0; h > 0} \frac{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(\tau_{ab}^k 1_{\tau_b < \tau_a})$ for $k \in \mathbb{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

<footnote n="1">¹up to an affine transformation</footnote>
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 $\infty$) 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

## 2.1 The approximation scheme

In this section we define the approximation process for a one-dimensional diffusion process on natural scale<fnref n="2" /> with state space $\mathbb{I}$, an open interval of $\mathbb{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 $\mathbb{I}$. We introduce incrementally this notion as follows, which we illustrate by Figure 1. Let $\mathbb{I}$ be an interval of $\mathbb{R}$.

* A *grid* over $\mathbb{I}$ is a countable subset of $\mathbb{I}$ with no accumulation points within $\mathbb{I}$.

* A *cell* $\mathbf{c}$ of a grid $\mathbf{g}$ is an open interval with endpoints in $\mathbf{g}$ with a single element of $\mathbf{g}$ in its interior, i.e.,

$$ \mathrm{Card} (\mathbf{g} \cap \mathbf{c}) = 1. $$

We denote with $\mathcal{C}(\mathbf{g})$ the set of all cells of the grid $\mathbf{g}$.

* We call $x \in \mathbb{I}$ the *center* of the cell $\mathbf{c} \in \mathcal{C}(\mathbf{g})$ iff

$$ \mathbf{c} \cap \mathbf{g} = \{x\}. $$

* Finally, a *covering grid* of $\mathbb{I}$ is a grid $\mathbf{g}$ such that $\mathbb{I} = \bigcup_{\mathbf{c} \in \mathcal{C}(\mathbf{g})} \mathbf{c}$.

Examples of covering grids of $(0, \infty)$ and $(-\infty, \infty)$ are $\{1/n; n \in \mathbb{N}\}$ and $\mathbb{Z}$ respectively.

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

$$ |\mathbf{g}| = \sup_{\mathbf{c} \in \mathcal{C}(\mathbf{g})} |\mathbf{c}|, \quad |\mathbf{g}|_X = \sup_{\mathbf{c} \in \mathcal{C}(\mathbf{g})} \{s(\mathbf{c})m(\mathbf{c})\}, \qquad (2) $$

where $|\mathbf{c}| = (b - a)$ and $s(\mathbf{c}) = s(b) - s(a)$ with $a$ and $b$ being the endpoints of $\mathbf{c}$. The convergence results will be expressed in the latter metric.

<footnote n="2">²Which means that $s(x) = x$.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: The covering grid $\mathbf{g} = \mathbb{Z}$ of $(-\infty, \infty)$ along with some of its cells $C(\mathbf{g})$.</figcaption></figure>

Let $X$ be a diffusion process on natural scale with state-space $\mathbb{I}$ and speed measure $m$. Let $X$ be defined on the family of filtered probability spaces $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, (P_x)_{x \in \mathbb{I}})$, where $P_x$ is the law of $X$ such that, for any $x \in \mathbb{I}$, $P_x\{X_0 = x\} = 1$. For any covering grid $\mathbf{g}$ of $\mathbb{I}$, we define the process $\tilde{X}^{\mathbf{g}} = (\tilde{X}_t^{\mathbf{g}})_{t \ge 0}$ as the asymmetric random walk with:

* state-space $\mathbf{g}$,

* initial distribution equal to the distribution of $X$ the first time it touches the grid,

* the same transition probabilities as $X$ over $\mathbf{g}$,

* conditional transition times that match the conditional expected transition times of $X$ over $\mathbf{g}$.

Thus, under $P_x$, if $a$ and $b$ are respectively the closest lower and upper elements to $x$ of $\mathbf{g}$,

$$ \tilde{X}_0^{\mathbf{g}} = \begin{cases} a, & \text{with probability } P_x(\tau_b < \tau_a), \\ b, & \text{with probability } P_x(\tau_a < \tau_b) = 1 - P_x(\tau_b < \tau_a), \end{cases} \quad (3) $$

where $\tau_a := \inf\{t > 0 : X_t = a\}$. For the rest of the trajectory, we define $\tau_{ab} := \tau_a \wedge \tau_b$ and $(T^{\mathbf{g}}(n))_{n \ge 0}$ as the consecutive jumping times of $(\tilde{X}_t^{\mathbf{g}})_{t \ge 0}$. Then, for all $k \in \mathbb{N}_0$ and $a < x < b$ adjacent points of $\mathbf{g}$,

$$ P(\tilde{X}_{T^{\mathbf{g}}(k+1)}^{\mathbf{g}} = b | \tilde{X}_{T^{\mathbf{g}}(k)}^{\mathbf{g}} = x) = P_x(\tau_b < \tau_a), \quad (4) $$

and

$$ T^{\mathbf{g}}(k+1) - T^{\mathbf{g}}(k) = \begin{cases} E_x(\tau_{ab} | \tau_b < \tau_a), & \text{on } \{\tilde{X}_{T^{\mathbf{g}}(k+1)}^{\mathbf{g}} = b\} \cap \{\tilde{X}_{T^{\mathbf{g}}(k)}^{\mathbf{g}} = x\}, \\ E_x(\tau_{ab} | \tau_a < \tau_b), & \text{on } \{\tilde{X}_{T^{\mathbf{g}}(k+1)}^{\mathbf{g}} = a\} \cap \{\tilde{X}_{T^{\mathbf{g}}(k)}^{\mathbf{g}} = x\}. \end{cases} \quad (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 $\mathbf{c}_x$ be the cell of the grid $\mathbf{g}$ containing $x$, i.e., $\mathbf{c}_x \in C(\mathbf{g})$ and $x = \mathbf{c}_x \cap \mathbf{g}$. From (4), (5) and Bayes' rule, if $\mathbf{c}_x = (a, b)$, both $P(\tilde{X}_{T^{\mathbf{g}}(k+1)}^{\mathbf{g}} = b | \tilde{X}_{T^{\mathbf{g}}(k)}^{\mathbf{g}} = x)$ and $T^{\mathbf{g}}(k+1) - T^{\mathbf{g}}(k)$ only
depend on $x$. Thus, if we know the quantities

$$ 
\begin{align} 
p^+[x, (a, b)] &= P_x(\tau_b < \tau_a), & T^+[x, (a, b)] &= E_x[\tau_{ab} | \tau_b < \tau_a], \tag{6} \\ 
p^-[x, (a, b)] &= P_x(\tau_a < \tau_b), & T^-[x, (a, b)] &= E_x[\tau_{ab} | \tau_a < \tau_b], 
\end{align} 
$$

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 (6). Practical examples are given in Sections 6 and 7.

This algorithm has been first introduced in [15] in the situation of SDE solutions with measurable coefficients, where the speed measure of the process satisfies

$$ c \, dx \le m(dx) \le C \, dx. \tag{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 [15] 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 $\tilde{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 $\tilde{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) + \rho \delta_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[\tau_{ab} \mathbf{1}_{\tau_b < \tau_a}] = \int_{(a,b)} G_{a,b}(x, \zeta) v_0(\zeta) m_c(d\zeta) + \rho G_{a,b}(x, 0) v_0(0), \tag{8} $$

where $v_0(x) = P_x(\tau_b < \tau_a)$.

For the convergence, we make the further assumption that the diffusion process $(X_t)_{t \ge 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) \ge k_1 \frac{1}{1+x^2} dx, \tag{9} $$

for all $x \in \mathbb{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
**Algorithm 1** Space-Time Markov Chain Approximation (STMCA) Algorithm

**Input:** $x$ initial value, $T$ time horizon, $\mathbf{g} = \{x_j\}_{j \in J}$ grid on $\mathbb{I}$

**Output:** $(\hat{\mathbf{X}}[k])_k$ consecutive values taken by the approximation process, $(\hat{t}[k])_k$ consecutive transition times

**Initialization:**

$$ \hat{t}[0] = 0, n = 0 \\ j = \arg \min_{i \in J} \{|x_i - x|\} \\ U \sim \text{Bernoulli}(p^+[x, (x_j, x_{j+1})]) $$

if $U == 1$ then

$$ j = j + 1 $$

end if

$$ \hat{\mathbf{X}}[0] = x_j $$

**Main loop:**

while $\hat{t}[n] < T$ do

$$ U \sim \text{Bernoulli}(p^+[x_j, (x_{j-1}, x_{j+1})]) $$

if $U == 1$ then

$$ j = j + 1 \\ \hat{t}[n+1] = \hat{t}[n] + T^+[x_j, (x_{j-1}, x_{j+1})] $$

else

$$ j = j - 1 \\ \hat{t}[n+1] = \hat{t}[n] + T^-[x_j, (x_{j-1}, x_{j+1})] $$

end if

$$ \hat{\mathbf{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 $\mu$ and $\nu$ of processes with càdlàg<fnref n="3" /> paths is defined as

$$ W_p[\mu, \nu] = \inf_{(\zeta, \xi) \sim \Gamma(\mu, \nu)} \left\| \|\zeta - \xi\|_{\infty} \right\|_{L^p}, $$

where by $\Gamma(\mu, \nu)$ we denote the collection of all measures with marginals $\mu$ and $\nu$.

**Theorem 2.1.** Let $X$ be a diffusion process with state-space $\mathbb{I}$ an interval of $\mathbb{R}$, on natural scale, whose speed measure satisfies Condition (9) for some constant $k_1 > 0$. Let $\mathbf{g}$ be a covering grid of $\mathbb{I}$. Then, for all $p \ge 1$, $\delta \in (0, \frac{1}{4} \wedge \frac{1}{p})$, $T > 0$ and $x \in \mathbb{I}$ there exists a constant $C > 0$ such that

$$ W_p \left[ \text{Law}((\tilde{\mathbf{X}}_t^{\mathbf{g}})_{t \in [0,T]}), \text{Law}((X_t)_{t \in [0,T]}) \right] \le C |\mathbf{g}|_X^\delta, \quad (10) $$

<footnote n="3">³This stands for “continue à droite avec une limite à gauche”, that is right-continuous with left-limit.</footnote>
where $|g|_X = \sup_{c \in C(g)} \{|c|m(c)\}$.

*Remark 3.* In the case where $m(dx) \ge k_1 dx$, the constant$^4$ $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 (10) can be replaced by $C|g|^{2\delta}$.

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

**Corollary 2.2.** For all $T > 0$, the processes $(\tilde{X}_t^g)_{t \in [0,T]}$ converges in law to $(X_t)_{t \in [0,T]}$ in the Skorokhod space $D([0,T], \mathbb{I})$ as $|g|_X \to 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 [6, 15]). 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.

## 2.2 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 $\mathbb{I}$ an open interval of $\mathbb{R}$, scale function $s$ and speed measure $m$. We assume that

* $s$ belongs to the Sobolev space $W^{1,1}(\mathbb{I})$, so from Theorem 8.2 of [9], as $s$ is continuous,

$$s(x) - s(y) = \int_y^x s'(t) dt,$$

for all $y, x$ in $\mathbb{I}$.

* there exists a $k_1 > 0$ such that for all $x \in \mathbb{I}$,

$$m(dx) \ge k_1 \frac{s'(x)}{1 + (s(x))^2} dx, \quad (11)$$

* the inverse of $s$ is $\alpha$-Hölder continuous, i.e., there exists a constant $C > 0$ such that for all $x \neq y \in \mathbb{I}$,

$$\frac{|s^{-1}(\bar{x}) - s^{-1}(\bar{y})|}{|\bar{x} - \bar{y}|^\alpha} \le C. \quad (12)$$

<footnote>$^4$This results in the bound of Theorem 2.1 not depending on the starting point of the diffusion.</footnote>
Given a grid $\mathbf{g}$, we consider the random walk $\tilde{X}_t^g$ 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 $\mathbf{g}$ be a covering grid over the state-space $\mathbb{I}$ of $X$. Then, for all $p \ge 1$, $\delta \in (0, \frac{1}{4} \wedge \frac{1}{p})$, $T > 0$ and $x \in \mathbb{I}$ there exists positive constants $C_1$, and $C_2$ such that

$$ \mathcal{W}_p[\text{Law}((\tilde{X}_t^g)_{t \in [0,T)}), \text{Law}((X_t)_{t \in [0,T)})] \le C_1 |\mathbf{g}|_X^\delta, $$

where $|\mathbf{g}|_X$ is defined in (2).

*Proof.* We define the proxy process $Y = (s(X_t))_{t \ge 0}$ which has state-space $s(\mathbb{I})$, scale function $s_Y(x) = x$ and speed measure $m_Y(dx) = m \circ s^{-1}(dx)$. From condition (11) and a change of variables, we get that $(Y_t)_{t \ge 0}$ satisfies condition (9) for the same constant $k_1$. We also define $\tilde{Y}_t^{s(\mathbf{g})}$ evolving according to Algorithm 1, with covering grid $s(\mathbf{g}) = \{s(x); x \in \mathbf{g}\}$. It can be defined on the canonical space of $X$ so that $s(\tilde{X}_t^g) = \tilde{Y}_t^{s(\mathbf{g})}$ almost surely. Thus, Condition (12) implies that

$$ |X_t - \tilde{X}_t^g| \le C|Y_t - \tilde{Y}_t^{s(\mathbf{g})}|^\alpha. $$

Along with the fact that

$$ |\mathbf{g}|_Y = \sup_{\mathbf{c} \in \mathcal{C}(\mathbf{g})} \{|s(\mathbf{c})|m_Y(s(\mathbf{c}))\}, $$

Theorem 2.1 implies Corollary 2.3. $\square$

## 2.3 Grid tuning

We observe that for all $\epsilon > 0$, in the case of a Brownian motion, Theorem 2.1 yields a convergence rate of $O(|\mathbf{g}|^{\frac{1-\epsilon}{2}})$ as $|\mathbf{g}| \to 0$, which is optimal from Donsker's invariance principle [12]. 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 $\mathbb{R}$, defined through $s$ and $m$ with

$$ s(x) = x, \quad m(dx) = 2dx + \rho\delta_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 $\mathbf{g}$ of $\mathbb{R}$ and $\mathbf{c} \in \mathcal{C}(\mathbf{g})$ with $\mathbf{c} \neq \mathbf{c}_0$ and $\mathbf{c}_0 \in \mathcal{C}(\mathbf{g})$ being the cell containing 0,

$$2|\mathbf{c}|^2 = m(\mathbf{c})|\mathbf{c}| \leq \sup_{\mathbf{c} \in \mathcal{C}(\mathbf{g})} m(\mathbf{c})|\mathbf{c}|, \\ \rho|\mathbf{c}_0| + 2|\mathbf{c}_0|^2 = m(\mathbf{c}_0)|\mathbf{c}_0| \leq \sup_{\mathbf{c} \in \mathcal{C}(\mathbf{g})} m(\mathbf{c})|\mathbf{c}|. \qquad (13)$$

From (13), for a uniform grid of step size h,

$$2h\rho < |\mathbf{g}|_X.$$

Thus, for any $\epsilon > 0$, the convergence rate given by Theorem 2.1 is $O(h^{\frac{1}{4}-\epsilon})$ as $|\mathbf{g}| \to 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\rho, h^2/2\rho)$, i.e.

$$\mathbf{g} = \left\{ \bigcup_{k \in \mathbb{Z}_+} \left\{ -\frac{h^2}{2\rho} - k\frac{h}{2} \right\} \right\} \cup \{0\} \cup \left\{ \bigcup_{k \in \mathbb{Z}_+} \left\{ \frac{h^2}{2\rho} + k\frac{h}{2} \right\} \right\}. \qquad (14)$$

As the approximation process is a random walk, for every $k$ steps it makes, it spends $O(\sqrt{k})$ steps in the cell containing 0 (see [10]). Thus, running the algorithm on either grid yields the same algorithmic complexity, while the convergence rate is improved to $O(h^{\frac{1}{2}-\epsilon})$ for the adapted grid. Numerical examples are given in Section 7.1.2.

The general case is covered by the following Corollary:

**Corollary 2.5.** Let $X$ be a diffusion process and $\tilde{X}^g$ the approximation process defined by Algorithm 1. Then, if $\mathbf{g}$ is a grid such that

$$|\mathbf{g}|_X \leq C|\mathbf{g}|^2, \qquad (15)$$

we can bound the $p$-Wasserstein distance between the laws of $(\tilde{X}_t^g)_{t \in [0,T)}$ and $(X_t)_{t \in [0,T)}$ in Theorem 2.1 by $|\mathbf{g}|^{2\delta}$ instead of $|\mathbf{g}|_X^{\delta}$. Thus, for all $\epsilon > 0$, the law of the random walk converges in any $p$-Wasserstein distance at the rate $O(|\mathbf{g}|^{(\frac{1}{2} \wedge \frac{2}{p})-\epsilon})$ instead of $O(|\mathbf{g}|^{(\frac{1}{4} \wedge \frac{1}{p})-\epsilon})$ as $|\mathbf{g}| \to 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 (7) and (11) 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 $\rho > 0$, one needs to consider the points $\{-h^2/\rho, 0, h^2/\rho\}$ to have a tuned grid,

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

Let $(\mu, \sigma)$ be a pair of functions satisfying the following condition:

**Condition 1.** The functions are measurable $\mathbb{R} \mapsto \mathbb{R}$ mappings and the SDE

$$dX_t = \mu(X_t) dt + \sigma(X_t) dB_t, \qquad (16)$$

has a unique weak solution, where $B$ is a standard Brownian motion.

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

$$s'(x) = e^{-\int_a^y \frac{2\mu(\zeta)}{\sigma^2(\zeta)} d\zeta}, \quad \text{and} \quad m(dx) = \frac{1}{s'(x)} \frac{2}{\sigma^2(x)} dx, \qquad (17)$$

with $s'$ being the right-derivative of $s$. Then, if $g = \{x_j\}_{j \in J}$ is a covering grid such that for a constant $C > 0$ and every $j \in J$

$$\left| (s(x_{j+1}) - s(x_{j-1})) \int_{x_{j-1}}^{x_{j+1}} \frac{2}{s'(\zeta)\sigma^2(\zeta)} d\zeta \right| \wedge (x_{j+1} - x_{j-1})^2 \le Ch^2, \qquad (18)$$

it satisfies (15) and $|g| = Ch$. Thus, from Corollary 2.5, for any $\epsilon > 0$ using such grids give us a convergence rate of $O(|g|^{1-\epsilon})$ 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:

$$\left| (s(y) - s(x_{j-1})) \int_{x_{j-1}}^{y} \frac{2}{s'(\zeta)\sigma^2(\zeta)} d\zeta \right| \le h^2/2. \qquad (19)$$

Then the next element of the grid is defined as

$$x_j = \begin{cases} x_j^{(0)} & \text{if } x_j^{(0)} - x_{j-1} \le h, \\ x_{j-1} + h & \text{if } x_j^{(0)} - x_{j-1} > h. \end{cases} \qquad (20)$$

The problem (19) can be solved numerically using a fixed point algorithm. Examples of STMCA simulations using tuned grids computed solving (19)-(20) are given in Figures 4 and 5.
## 2.4 The case of diffusions with boundary conditions

When presenting the results and the structure of the scheme, we considered only processes where $\mathbb{I}$ is an open set, thus considering diffusion with unreachable boundaries. Our results also adapt to the situation where either $\ell$ 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 $\mathbb{I} = [0, +\infty)$ (the adaptation to $\mathbb{I} = (\ell, r]$ or $\mathbb{I} = [\ell, r]$ or $\mathbb{I} = [\ell, r)$ with $\ell \in \mathbb{R} \cup \{-\infty\}$ and $r \in \mathbb{R} \cup \{+\infty\}$ is straightforward).

It is well known (see for instance Section 5.11 of Itô's book [23]) that the finite boundary 0 can be of four types. Setting, for some fixed $c > 0$,

$$ \mathcal{I} = \iint_{0<y<x<c} m(dx) \, dy, \quad \mathcal{II} = \iint_{0<y<x<c} m(dy) \, dx, $$

then

* 0 is an exit boundary if $\mathcal{I} < \infty$ and $\mathcal{II} = \infty$,

* 0 is a regular boundary if $\mathcal{I} < \infty$ and $\mathcal{II} < \infty$,

* 0 is a natural boundary if $\mathcal{I} = \infty$ and $\mathcal{II} = \infty$,

* 0 is an entrance boundary if $\mathcal{I} = \infty$ and $\mathcal{II} < \infty$.

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 $\mathbb{I}$ containing 0 and by adapting the dynamics of $\tilde{X}^g$ as follows. The dynamic of $\tilde{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 $\tilde{X}^g$ when it reaches 0;

* if 0 is a reflecting regular boundary, then the process $\tilde{X}^g$ jumps from 0 to $b := \min g \setminus \{0\}$ with probability 1 and after a time $\int_{[0,b]} (b - \zeta) m(d\zeta)$. We emphasize that in this configuration, 0 may be a sticky boundary (i.e., with $m(0) \in (0, +\infty)$).

In both cases of reflection and absorption at 0, the boundary is attainable. Thus, 0 must be a point of any covering grid of $\mathbb{I} = [0, \infty)$. 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 $\mathbb{R}$ on $\mathbb{I}$ with endpoint in $g$. For example $g_0 = \mathbb{Z}_+$ is a covering grid of $[0, \infty)$ and $C(g_0) = \{[0, 1), (n-1, n+1)_{n \in \mathbb{N}}\}$.
<figure><img src="image_3.png" /><figcaption>Figure 2: The covering grid $\mathbf{g} = \mathbb{Z}^+$ of $[0, \infty)$ along with its first cells $\mathcal{C}(\mathbf{g})$.</figcaption></figure>

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,

$$E_0[\tau_b] = \int_{[0,b)} (b - \zeta) m(d\zeta).$$

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.

## 2.5 Markovian embedding

The consecutive values of the process $(\tilde{X}_t)_{t \ge 0}$ defined through (3)-(5) form a Markov chain with, by construction, the same transition probabilities as $(X_t)_{t \ge 0}$ on $\mathbf{g}$. We define the embedding times of $(X_t)_{t \ge 0}$ in $\mathbf{g}$ as,

$$\begin{cases} \tau_0^{\mathbf{g}} = 0, \\ \tau_k^{\mathbf{g}} = \inf \{t > \tau_{k-1}^{\mathbf{g}} : X_t \in \mathbf{g} \setminus \{X_{\tau_{k-1}^{\mathbf{g}}}\}\}, & \forall k \ge 1. \end{cases} \quad (21)$$

As both $\tilde{X}_{Tg(n)}$ and $X_{\tau_n^g}$ are both Markov chains with the same transition probabilities with $\tilde{X}_0$ forced to be equal in law to $X_{\tau_1^g}$ (see Section 2.1), the following equality in law holds,

$$\mathrm{Law}(\tilde{X}_{Tg(n)}; n \ge 0) = \mathrm{Law}(X_{\tau_n^g}; n \ge 1). \quad (22)$$

We define $K^g(t)$ as the inverse of $T^g(n)$, i.e.,

$$K^g(t) = \inf \left\{ n \in \mathbb{N} : \sum_{k=1}^{n} E \left[ \tau_k^g - \tau_{k-1}^g \mid X_{\tau_{k-1}^g}, X_{\tau_k^g} \right] > t \right\}. \quad (23)$$

Thus, we get the following Proposition.
**Proposition 2.6.** Let $(X_t)_{t \ge 0}$ be a diffusion process, $g$ a grid defined over its state space $\mathbb{I}$ and $\tilde{X}_t^g$ be the approximation process defined in (4) and (5). Then, if $(\tau_n^g)_{n \ge 0}$ are the embedding times of $(X_t)_{t \ge 0}$ in $g$, the following equality in law holds,

$$ \mathrm{Law}(\tilde{X}_t; t \ge 0) = \mathrm{Law}(X_{\tau_{Kg(t)}}; t \ge 0), $$

where $K^g(t)$ is the random index defined in (23).

## 3 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(\tau_b < \tau_a)$ and conditional transition times $E_x(\tau_{ab}|\tau_b < \tau_a)$. In this section we show that quantities of the form $v_k(x) = E_x(\tau_{ab}^k \mathbb{1}_{\tau_b < \tau_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 \ge 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(\tau_b < \tau_a), \quad v_1(x) = E_x(\tau_{ab} \mathbb{1}_{\tau_b < \tau_a}) \quad \text{and} \quad \bar{v}_1(x) = E_x(\tau_{ab} \mathbb{1}_{\tau_a < \tau_b}). $$

The quantities of (6) are then

$$ \begin{aligned} p^+[x, (a, b)] &= v_0(x), & T^+[x, (a, b)] &= \frac{v_1(x)}{v_0(x)}, \\ p^-[x, (a, b)] &= 1 - v_0(x), & T^-[x, (a, b)] &= \frac{\bar{v}_1(x)}{1 - v_0(x)}. \end{aligned} \quad (24) $$

**Proposition 3.1.** The function $v_0(x) = P_x(\tau_b < \tau_a)$ is solution to the problem with Dirichlet boundary conditions

$$ \left\{ \begin{array}{l} L u = 0, \quad x \in (a, b), \\ u(a) = 0, \\ u(b) = 1, \end{array} \right. \quad (25) $$

where $L$ is the infinitesimal generator of $(X_t)_{t \ge 0}$, which also implies that $v_0 \in \mathrm{dom}(L)$.

*Proof.* Let $x \in (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 \frac{s(\cdot) - s(a)}{s(b) - s(a)} = D_m \frac{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_b^0$, we deduce that $v_0 \in \text{dom}(L)$ and $L v_0 = 0$. Under $P_b$, the stopping time $\tau_b$ equal 0 a.s. and the process has a.s. continuous trajectories, hence $\tau_a > 0$ a.s., i.e.,

$$v_0(b) = P_b(\tau_b < \tau_a) = P_b(0 < \tau_a) = 1.$$

This, along with the symmetrical argument, allow us to retrieve the boundary conditions of (25). $\square$

**Proposition 3.2.** For every $k \in \mathbb{N}$, let $v_k$ be the function defined for every $x \in (a, b)$ by $v_k(x) = E_x(\tau_{ab}^k 1_{\tau_b < \tau_a})$. Then,

$$v_k(x) = k \int_{(a,b)} G_{a,b}(x, \zeta) v_{k-1}(\zeta) m(d\zeta). \qquad (26)$$

*Proof.* Since $\int_0^{\tau_{ab}} (\tau_{ab} - t)^{k-1} dt = \tau_{ab}^k/k$,

$$v_k(x) = k E_x \left[ 1_{\tau_b < \tau_a} \int_0^{\tau_{ab}} (\tau_{ab} - t)^{k-1} dt \right] = k E_x \left[ 1_{\tau_b < \tau_a} \int_0^{\infty} 1_{t \le \tau_{ab}} (\tau_{ab} - t)^{k-1} dt \right].$$

From the Markov property, by conditioning on $\mathcal{F}_t$ and as $1_{t \le \tau_{ab}}$ is $\mathcal{F}_t$-measurable,

$$v_k(x) = k E_x \left[ \int_0^\infty 1_{t \le \tau_{ab}} E \left[ 1_{\tau_b < \tau_a} (\tau_{ab} - t)^{k-1} \mid \mathcal{F}_t \right] dt \right] = k E_x \left[ \int_0^{\tau_{ab}} E_{X_t} \left[ 1_{\tau_b < \tau_a} \tau_{ab}^{k-1} \right] dt \right].$$

The equality (26) results by applying directly Green's formula. $\square$

**Proposition 3.3.** The function $v_k(x) = E_x(\tau_{ab}^k 1_{\tau_b < \tau_a})$ is solution to the problem with Dirichlet boundary conditions

$$\begin{cases} L u = -kv_{k-1}, & x \in (a, b), \\ u(a) = 0, \\ u(b) = 0. \end{cases}$$

**Lemma 3.4.** Let $g(x) = \int_{(a,b)} G_{a,b}(x, y) f(y)m(dy)$, where $f \in C_b^0(a, b)$ and $G_{a,b}(x, y)$ is the Green function defined in (1). Then $g \in \text{dom}(L)$ and

$$L g(x) = -f(x), \quad \forall x \in (a, b).$$

*Proof.* Let $x \in (a, b)$. Using the $D_m D_s$ factorization of L and the dominated convergence theorem we get

$$\begin{aligned} L g(x) &= D_m D_s \int_{(a,b)} \left[ 1_{y<x} \frac{(s(y)-s(a))(s(b)-s(x))}{s(b)-s(a)} + 1_{y \ge x} \frac{(s(x)-s(a))(s(b)-s(y))}{s(b)-s(a)} \right] f(y)m(dy) \\ &= -D_m \int_{(a,x)} v_0(y)f(y)m(dy) + D_m \int_{[x,b)} (1-v_0(y))f(y)m(dy) \\ &= -v_0(x)f(x) - (1-v_0(x))f(x) = -f(x). \end{aligned}$$
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 \in \text{dom}(L)$ and that on $(a, b)$ we have $Lg = -f$. $\square$

*Proof of Proposition 3.3.* From Proposition 3.2,

$$v_k(x) = \int_{(a,b)} G_{a,b}(x, y)kv_{k-1}(y)m(dy).$$

As $v_0 \in C_b^0(a, b)$, from Lemma 3.4 we deduce iteratively that $v_k \in C_b^0(a, b)$ for all $k \in \mathbb{N}$ and that $L v_k = -kv_{k-1}$ on $(a, b)$. For the boundary conditions, we observe that $\tau_{ab} = \tau_a \wedge \tau_b = 0$ a.s. under $P_a$, so for $k \ge 1$

$$v_k(a) = E_a(\tau_{ab}^k 1_{\tau_b < \tau_a}) \le E_a(\tau_{ab}^k) = 0.$$

With the same argument we show that $v_k(b) \le 0$ and as they are obviously positive quantities $v_k(a) = v_k(b) = 0$. $\square$

# 4 Convergence of the embedding times

In order to prove the convergence of the process $(\tilde{X}_t^g)_{t \ge 0}$ we need to control quantities of the form $E_x[\sup_{t \le T} |t - \tau_{Kg(T)}^g|^p]$, where $\tau_k^g$ 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 \in \mathcal{C}(g)} \{s(c)m(c)\} \qquad (27)$$

defined in (2). If there exists a constant $K > 0$ such that $m(c) \le K$ for all $c \in \mathcal{C}(g)$, then, $|c|m(c) \le K|c|$ and $|g|_X \le 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.

## 4.1 Bounds on the conditional moments of the exit times

**Lemma 4.1.** Let $v_k(x) = E_x(\tau_{ab}^k 1_{\tau_b < \tau_a})$, then for all $k \in \mathbb{N}$

$$\left\| \frac{v_k}{v_{k-1}} \right\|_{\infty} \le k|g|_X. \qquad (28)$$
Proof. We first observe that

$$G_{a,b}(x, y) \frac{v_0(y)}{v_0(x)} = \begin{cases} \frac{(y-a)(b-x)}{(b-a)} \frac{y-a}{x-a}, & x > y, \\ \frac{(x-a)(b-y)}{(b-a)} \frac{y-a}{x-a}, & x \le y. \end{cases}$$

As for $x > y$ the ratio $\frac{y-a}{x-a} < 1$,

$$G_{a,b}(x, y) \frac{v_0(y)}{v_0(x)} \le \begin{cases} \frac{(y-a)(b-x)}{(b-a)}, & x > y, \\ \frac{(y-a)(b-y)}{(b-a)}, & x \le y, \end{cases} \quad (29)$$

which is bounded by $(b-a)$ in both cases as $x, y \in (a, b)$. Thus,

$$\int_a^b G_{a,b}(x, y) v_0(y) m(dy) = v_0(x) \int_a^b G_{a,b}(x, y) \frac{v_0(y)}{v_0(x)} m(dy) \le v_0(x)(b-a)m((a, b)). \quad (30)$$

From Proposition 3.2, $v_k(x) = k \int_{(a,b)} G_{a,b}(x, y)v_{k-1}(y)m(dy)$ and

$$v_k(x) = k! \int_{(a,b)} G_{a,b}(x, x_k) \int_{(a,b)} G_{a,b}(x_k, x_{k-1}) \cdots \int_{(a,b)} G_{a,b}(x_2, x_1) v_0(x_1) m(dx_1) \cdots m(dx_k). \quad (31)$$

From (30) and (31),

$$v_k(x) \le k(b-a)m((a,b))v_{k-1}(x).$$

Since $(b-a)m((a,b)) \le |g|_X$, we get the desired result on $v_k/v_{k-1}$. $\square$

**Corollary 4.2.** Let $k \in \mathbb{N}$. Then, we have the following bound for $v_k(x) = E_x(\tau_{ab}^k 1_{\tau_b < \tau_a})$,

$$\|v_k\|_{\infty} \le k! |g|_X^k.$$

**Lemma 4.3.** Let $X$ be a diffusion process with state space $\mathbb{I}$ an interval of $\mathbb{R}$, $m$ the speed measure of $X$, $a, x, b \in \mathbb{I}$ such that $a < x < b$, $c_x = (a, b)$, $|c_x|_X = m(c_x)|b-a|$ and $\lambda > 0$ such that $\lambda|c_x|_X \in (0, 1)$. Then,

$$E_x(e^{\lambda \tau_{ab}}) \le \exp \left( \frac{\lambda}{1 - \lambda |c_x|_X} E_x(\tau_{ab}) \right). \quad (32)$$

Proof. Developing the exponential series,

$$E_x(e^{\lambda \tau_{ab}}) = 1 + \sum_{N=1}^{\infty} \frac{\lambda^N}{N!} E_x(\tau_{ab}^N) = 1 + E_x(\tau_{ab}) \sum_{N=1}^{\infty} \frac{\lambda^N}{N!} \frac{E_x(\tau_{ab}^N)}{E_x(\tau_{ab})}.$$

From Corollary 4.2, we can bound the ratio of expected values by $(N-1)!|c_x|_X^{N-1}$. As $\lambda|c_x|_X \in (0, 1)$,

$$E_x(e^{\lambda \tau_{ab}}) \le 1 + E_x(\tau_{ab}) \sum_{N=1}^{\infty} \frac{\lambda}{N} (\lambda |c_x|_X)^{N-1} \le 1 + E_x(\tau_{ab}) \frac{\lambda}{1 - \lambda |c_x|_X}.$$

Thus, we only need to apply the inequality $1 + x \le e^x$ to get (32). $\square$
**Lemma 4.4.** Let $t, M > 0$ and $\lambda > 0$ such that $\lambda|g|_X \in (0, 1)$. Then,

$$P_x(\tau_{Kg(t)}^g > M) \le e^{-\lambda M} e^{\lambda \frac{t+|g|_X}{1-\lambda|g|_X}}. \quad (33)$$

*Proof.* From Markov's inequality,

$$P_x(\tau_{Kg(t)}^g > M) \le e^{-\lambda M} E_x \left[ e^{\lambda \tau_{Kg(t)}} \right] = e^{-\lambda M} E_x \left[ e^{\lambda \sum_{k=1}^{Kg(t)} (\tau_k - \tau_{k-1})} \right]. \quad (34)$$

Conditioning on the $\sigma$-algebra $\mathcal{B}$ generated by the trajectory of $X_t$ on the grid $g$, i.e., $\mathcal{B} = \sigma\{X_{\tau_k}; k \in \mathbb{N}_0\}$, and as $K^g(t)$ is $\mathcal{B}$-measurable,

$$E_x \left[ e^{\lambda \sum_{k=1}^{Kg(t)} (\tau_k - \tau_{k-1})} \right] = E_x \left[ \prod_{k=1}^{Kg(t)} E \left[ e^{\lambda(\tau_k - \tau_{k-1})} \mid \mathcal{B} \right] \right]. \quad (35)$$

From the definition of $|g|_X$, $\lambda|c|_X \le \lambda|g|_X \in (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(\tau_{Kg(t)}^g > M) \le e^{-\lambda M} E_x \left[ \exp \left( \frac{\lambda}{1-\lambda|g|_X} \sum_{k=1}^{Kg(t)} E(\tau_k - \tau_{k-1}|\mathcal{B}) \right) \right] \quad (36)$$

From (23),

$$t < \sum_{k=1}^{Kg(t)} E(\tau_k - \tau_{k-1}|\mathcal{B}) \le t + E(\tau_{Kg(t)+1} - \tau_{Kg(t)}|\mathcal{B})$$

Thus, from (28),

$$\sum_{k=1}^{Kg(t)} E(\tau_k - \tau_{k-1}|\mathcal{B}) \le t + \|v_1/v_0\|_{\infty} \le t + |g|_X, \quad (37)$$

From (36) and (37), we get (33). $\square$

## 4.2 Convergence of the embedding times

For this section, let $X$ be a diffusion process with state-space $\mathbb{I}$ an interval of $\mathbb{R}$, $g$ a covering grid of $\mathbb{I}$, $\tau_k^g$ the embedding times of $X$ in $g$ as defined in (21) and $K^g(t)$ as defined in (23).

**Lemma 4.5.** For any $T > 0$,

$$\sum_{k=1}^{Kg(T)} \mathrm{Var} \left( \tau_k^g - \tau_{k-1}^g \mid X_{\tau_{k-1}^g}, X_{\tau_k^g} \right) \le 2|g|_X(T + |g|_X),$$

where $|\cdot|_X$ is defined in (2).
Proof. For all $x \in \mathbf{g}$, let $c_x$ be the cell of $\mathbf{g}$ containing $x$. Then,

$$
\begin{align*}
\sum_{k=1}^{K^{\mathbf{g}}(T)} \mathrm{Var}(\tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}} | X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}})
&= \sum_{k=1}^{K^{\mathbf{g}}(T)} \mathrm{E} \left[ (\tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}})^2 \middle| X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right] - \left( \mathrm{E} \left[ \tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}} \middle| X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right] \right)^2 \\
&\le \sup_{k \le K^{\mathbf{g}}(t)} \left\{ \frac{\mathrm{E} \left[ (\tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}})^2 \middle| X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right]}{\mathrm{E} \left[ \tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}} \middle| X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right]} \right\} \sum_{k=1}^{K^{\mathbf{g}}(T)} \mathrm{E} \left[ (\tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}}) \middle| X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right].
\end{align*}
$$

So from Lemma 4.1 and the definition of $K^{\mathbf{g}}(t)$,

$$
\begin{align*}
\sum_{k=1}^{K^{\mathbf{g}}(T)} \mathrm{Var}(\tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}} | X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}})
&\le \left\| \frac{v_2}{v_1} \right\|_{\infty} \sum_{k=1}^{K^{\mathbf{g}}(T)} \mathrm{E} \left[ (\tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}}) | X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right] \\
&\le \left\| \frac{v_2}{v_1} \right\|_{\infty} \left( T + \left\| \frac{v_1}{v_0} \right\|_{\infty} \right) \le 2 |\mathbf{g}|_X (T + |\mathbf{g}|_X),
\end{align*}
$$

which is the desired inequality.

<figure><img src="image_4.png" /></figure>

**Proposition 4.6.** Let $(\mathcal{F}_t)_{t \ge 0}$ be the canonical filtration of $X$. Let also $\mathcal{A}_n = \mathcal{F}_{\tau_n^{\mathbf{g}}}$, $\mathcal{B} = \sigma((X_{\tau_k^{\mathbf{g}}})_{k \in \mathbb{N}_0})$ and $\Delta \tau_k^{\mathbf{g}} = \tau_k^{\mathbf{g}} - \tau_{k-1}^{\mathbf{g}}$. If we define the augmented filtration $\mathcal{G}_n = \mathcal{A}_n \vee \mathcal{B}$, then the process

$$
M_n = \sum_{k=1}^{n} \Delta \tau_k^{\mathbf{g}} - \mathrm{E} \left[ \Delta \tau_k^{\mathbf{g}} | X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right],
$$

is a $\mathcal{G}_n$-martingale.

Proof. For $m \le n$,

$$
\begin{align*}
\mathrm{E} [M_n | \mathcal{G}_m] &= \mathrm{E} [M_n | \mathcal{A}_m, \mathcal{B}] \\
&= \mathrm{E} \left[ \sum_{k=1}^{m} \Delta \tau_k^{\mathbf{g}} - \mathrm{E} \left[ \Delta \tau_k^{\mathbf{g}} | X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right] + \sum_{k=m+1}^{n} \Delta \tau_k^{\mathbf{g}} - \mathrm{E} \left[ \Delta \tau_k^{\mathbf{g}} | X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right] \middle| \mathcal{A}_m, \mathcal{B} \right].
\end{align*}
$$

From Lemma A.1 and as $\tau_k^{\mathbf{g}}$ is $\mathcal{A}_k$ measurable,

$$
\mathrm{E} [M_n | \mathcal{G}_m] = M_m + \sum_{k=m+1}^{n} \mathrm{E}_{X_{\tau_m^{\mathbf{g}}}} \left[ \Delta \tau_k^{\mathbf{g}} - \mathrm{E} \left[ \Delta \tau_k^{\mathbf{g}} | X_{\tau_{k-1}^{\mathbf{g}}}, X_{\tau_k^{\mathbf{g}}} \right] | \mathcal{B} \right] = M_m.
$$

This proves the result.

<figure><img src="image_5.png" /></figure>
**Theorem 4.7.** For any $T > 0$,

$$ E \left[ \sup_{t \in [0, T]} |\tau_{Kg(t)}^{g} - t|^2 \right] \le 2|g|_X (4(T + |g|_X) + 1), \quad (38) $$

where $|\cdot|_X$ is defined in (2).

*Proof.* The convexity inequality $(a+b)^p \le 2^{p-1}(a^p + b^p)$ yields for $p = 2$

$$ \begin{aligned} E \left[ \sup_{t \in [0, T]} |\tau_{Kg(t)}^{g} - t|^2 \right] \le 2 E \left[ \sup_{t \in [0, T]} \left| \tau_{Kg(t)}^{g} - \sum_{k=1}^{Kg(t)} E \left[ \Delta \tau_k^{g} | X_{\tau_{k-1}^{g}}, X_{\tau_k^{g}} \right] \right|^2 \right] \\ + 2 E \left[ \sup_{t \in [0, T]} \left| \sum_{k=1}^{Kg(t)} E \left[ \Delta \tau_k^{g} | X_{\tau_{k-1}^{g}}, X_{\tau_k^{g}} \right] - t \right|^2 \right]. \end{aligned} \quad (39) $$

From Proposition 4.6, the term $\sum_{k=1}^n \Delta \tau_k^g - E[\Delta \tau_k^g | X_{\tau_{k-1}^g}, X_{\tau_k^g}]$ is a $G_n$-martingale, where $G_n = \mathcal{F}_{\tau_n^g} \vee \mathcal{B}$ and $\mathcal{B} = \sigma((X_{\tau_k^g})_{k \in \mathbb{N}_0})$. Thus, from Doob's $L^p$ inequality,

$$ E \left[ \sup_{t \in [0, T]} |M_{Kg(t)}|^2 \right] = E \left[ \sup_{k \le Kg(T)} |M_k|^2 \right] \le 2 E \left[ |M_{Kg(T)}|^2 \right]. $$

By conditioning on $\mathcal{B}$, from Lemma A.1,

$$ \begin{aligned} & E \left[ \sup_{t \in [0, T]} \left| \sum_{k=1}^{Kg(t)} \Delta \tau_k^g - E \left[ \Delta \tau_k^g | X_{\tau_{k-1}^g}, X_{\tau_k^g} \right] \right|^2 \right] \\ & \le 2 E \left[ \left( \sum_{k=1}^{Kg(T)} \Delta \tau_k^g - E \left[ \Delta \tau_k^g | X_{\tau_{k-1}^g}, X_{\tau_k^g} \right] \right)^2 \right] \\ & = 2 E \left[ \sum_{k=1}^{Kg(T)} \mathrm{Var} \left( \Delta \tau_k^g | X_{\tau_{k-1}^g}, X_{\tau_k^g} \right) \right], \end{aligned} $$

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 (23) and since $K^g(t) \ge 1$ for any $t > 0$,

$$ \begin{aligned} & \sum_{k=1}^{Kg(t)} E \left[ \tau_k^g - \tau_{k-1}^g | X_{\tau_{k-1}^g}, X_{\tau_k^g} \right] - t \le E \left[ \tau_{Kg(t)}^g - \tau_{Kg(t)-1}^g | X_{\tau_{Kg(t)-1}^g}, X_{\tau_{Kg(t)}^g} \right] \\ & \le \left\| \frac{v_1}{v_0} \right\|_{\infty} \le |g|_X. \end{aligned} $$

So having bounded both additive parts of the right hand side of (39), we get (38). $\square$
# 5 Convergence rate of the Markov chain

## 5.1 Moment bounds

In this section we prove that Theorem 3.1 of [5] 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 $\mathbb{R}$ and a speed measure $m_X$ that satisfies Condition (9). Then, for all $p \ge 2$, there exist two constants $C$, $C' > 0$ such that

$$E_x \left[ \sup_{t \in [0, T]} |X_t - x|^p \right] \le C' [1 + |x|^p] e^{CT}, \quad (40)$$

where $k_1$ is a constant such that Condition (9) is satisfied, $C \le 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) = \mathbb{1}_{|x|<1} \frac{k_1}{2} dx + \mathbb{1}_{|x|\ge 1} \frac{k_1}{2x^2} dx.$$

We note that $\frac{k_1}{2x^2} \le \frac{k_1}{1+x^2}$ for all $x \ge 1$. The dynamic of $Z$ can be shown to be

$$dZ_t = \begin{cases} \frac{2}{\sqrt{k_1}} dB_t, & \text{for } |Z_t| < 1, \\ \frac{2}{\sqrt{k_1}} Z_t dB_t, & \text{for } |Z_t| \ge 1, \end{cases}$$

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 \ge 0$, $X_t = B_{\gamma_X(t)}$ and $Z_t = W_{\gamma_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 $\gamma_X(t)$, $\gamma_Z(t)$ being the respective right-inverses<fnref n="5" /> of

$$A_X(t) = \frac{1}{2} \int_{\mathbb{I}} L_t^x(B) m_X(dx) \quad \text{and} \quad A_Z(t) = \frac{1}{2} \int_{\mathbb{I}} L_t^x(W) m_Z(dx).$$

Using the same underlying Brownian motion in these definitions, we have $A_Z(t) \le A_X(t)$, and hence $\gamma_X(t) \le \gamma_Z(t)$. Thus,

$$\begin{aligned} E_x \left[ \sup_{t \in [0, T]} |X_t - x|^p \right] &= E_x \left[ \sup_{t \in [0, \gamma_X(T)]} |B_t - x|^p \right] \\ &\le E_x \left[ \sup_{t \in [0, \gamma_Z(T)]} |B_t - x|^p \right] = E_x \left[ \sup_{t \in [0, T]} |Z_t - x|^p \right]. \quad (41) \end{aligned}$$

$$f^{-1}(x) = \inf\{\zeta \ge 0 : f(\zeta) > x\}.$$

<footnote n="5">⁵The right-inverse of a function $f$ is given by,</footnote>
Thus, from Doob's $L^p$ and convexity inequalities,

$$ E_x \left[ \sup_{t \in [0, T]} |X_t - x|^p \right] \le \frac{2^{p-1}p}{p-1} \left( E_x \left[ |Z_T|^p \right] + |x|^p \right). \quad (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$,

$$ \begin{aligned} & E_x \left[ (|Z_t| - 1)^{3q} 1_{|Z_t| \ge 1} + 1 \right] \\ & \le (|x| - 1)^{3q} 1_{|x| \ge 1} + 1 + \int_0^t \frac{3q}{2} (3q - 1) E_x \left[ (|Z_s| - 1)^{3q-2} 1_{|Z_s| \ge 1} \frac{4Z_s^2}{k_1} \right] ds. \end{aligned} $$

Using the inequality $|x - 1|^{3q-2}x^2 \le 4|x - 1|^{3q} + 4$ for all $x \ge 1$,

$$ \begin{aligned} & E_x \left[ (|Z_t| - 1)^{3q} 1_{|Z_t| \ge 1} + 1 \right] \\ & \le (|x| - 1)^{3q} 1_{|x| \ge 1} + 1 + \int_0^t \frac{24q}{k_1} (3q - 1) E_x \left[ (|Z_s| - 1)^{3q} 1_{|Z_s| \ge 1} + 1 \right] ds. \end{aligned} $$

Using Gronwall's Lemma, we deduce that, for all $t \ge 0$,

$$ E_x \left[ (|Z_t| - 1)^{3q} 1_{|Z_t| \ge 1} + 1 \right] \le \left[ (|x| - 1)^{3q} 1_{|x| \ge 1} + 1 \right] e^{24q(3q-1)t/k_1}. $$

Let $C'_q > 0$ be a constant such that $(|x| - 1)^{3q} 1_{|x| \ge 1} + 1 \ge C'_q |x|^{3q}$ for all $x \in \mathbb{R}$. Thus,

$$ E_x \left[ |Z_t|^{3q} \right] \le C'_q \left[ (|x| - 1)^{3q} 1_{|x| \ge 1} + 1 \right] e^{24q(3q-1)t/k_1}, \quad \forall t \ge 0. $$

Hence, for all $^6 p > 2$, we have

$$ E_x \left[ |Z_t|^p \right] \le C'_{p/3} \left[ 1 + |x|^p \right] e^{8p(p-1)t/k_1}, \quad \forall t \ge 0. \quad (43) $$

Replacing (43) in (42), we get the bound (40) for $C' = \frac{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 $\mathbb{I}$ an interval of $\mathbb{R}$, on natural scale and with a speed measure $m_X$ that satisfies Condition (9). Then, for each $T > 0$ and $\gamma \in (0, \frac{1}{2})$, there exists a constant $C > 0$ such that

$$ \left\| \sup_{s \neq t \leq T} \frac{|X_t - X_s|}{|t - s|^{\gamma}} \right\|_{L^p(\mathbb{P}_x)} \leq C(1 + |x|). \quad (44) $$

The result also holds in the presence of a reflecting boundary $\zeta \in \mathbb{I}$.

<footnote n="6">⁶The extension to $p = 2$ is straightforward.</footnote>
Proof. For the non-reflecting case, the proof works using the same arguments as in the proof of Theorem 3.1 in [5].

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 $\zeta \in \mathbb{I}$. We observe that

$$X = |X^{\circ} - \zeta| + \zeta \qquad (45)$$

in law, where $X^{\circ}$ is the non-reflecting diffusion on natural scale with speed measure

$$m_{X^{\circ}}(dx) = 1_{x \ge \zeta} m(dx) + 1_{x < \zeta} m(2\zeta - dx)$$

which also satisfies (11) for the same constant $k_1$. From (45) and the triangle inequality,

$$\frac{|X_t - X_s|}{|t - s|^{\gamma}} = \frac{|X_t^{\circ} - \zeta| - |X_s^{\circ} - \zeta|}{|t - s|^{\gamma}} \le \frac{|X_t^{\circ} - X_s^{\circ}|}{|t - s|^{\gamma}}. \qquad (46)$$

The diffusion process $X^{\circ}$ is non-reflecting and its speed measure $m_{X^{\circ}}$ satisfies (11). Thus, (44) holds for $X^{\circ}$ for a constant $C > 0$ and from (46),

$$\left\| \sup_{s \neq t \leq T} \frac{|X_t - X_s|}{|t - s|^{\gamma}} \right\|_{L^p(P_x)} \leq \left\| \sup_{s \neq t \leq T} \frac{|X_t^{\circ} - X_s^{\circ}|}{|t - s|^{\gamma}} \right\|_{L^p(P_x)} \leq C(1 + |x|). \quad \square$$

## 5.2 Proof of the convergence rate for a process on natural scale

*Proof of Theorem 2.1.* From Proposition 2.6,

$$\begin{aligned} \mathcal{W}_p \left[ \text{Law}((\tilde{X}_t^g)_{t \in [0,T]}), \text{Law}((X_t)_{t \in [0,T]}) \right] \\ = \inf \left\{ \left\| \zeta - \xi \right\|_{L^p}; \zeta \sim \text{Law}((\tilde{X}_t^g)_{t \in [0,T]}), \xi \sim \text{Law}((X_t)_{t \in [0,T]}) \right\} \\ \le \left\| \sup_{t \in [0,T]} |X_{\tau_{Kg(t)}}^g - X_t| \right\|_{L^p(P_x)}. \end{aligned}$$

Let $M > T$. From Minkowski inequality,

$$\left\| \sup_{t \in [0, T]} |X_{\tau_{Kg(t)}}^g - X_t| \right\|_{L^p(P_x)} \le \left\| \mathbb{1}_{\tau_{Kg(T)}^g \le M} \sup_{t \in [0, T]} |X_{\tau_{Kg(t)}}^g - X_t| \right\|_{L^p(P_x)} \\ + \left\| \mathbb{1}_{\tau_{Kg(T)}^g > M} \sup_{t \in [0, T]} |X_{\tau_{Kg(t)}}^g - X_t| \right\|_{L^p(P_x)} \quad (47)$$
For the first additive term of (47), for any $\gamma > 0$, by multiplying and dividing by $|\tau_{Kg(t)}^g - t|^\gamma$,

$$
\begin{aligned}
& \left\| 1_{\tau_{Kg(T)}^g} \le M \sup_{t \in [0,T]} |X_{\tau_{Kg(T)}^g} - X_t| \right\|_{L^p(P_x)} \\
& = \left\| 1_{\tau_{Kg(T)}^g} \le M \sup_{t \in [0,T]} \left\{ \frac{|X_{\tau_{Kg(t)}^g} - X_t|}{|\tau_{Kg(t)}^g - t|^\gamma} |\tau_{Kg(t)}^g - t|^\gamma \right\} \right\|_{L^p(P_x)} \\
& \le \left\| 1_{\tau_{Kg(T)}^g} \le M \sup_{t \in [0,T]} \left\{ \frac{|X_{\tau_{Kg(t)}^g} - X_t|}{|\tau_{Kg(t)}^g - t|^\gamma} \right\} \sup_{t \in [0,T]} |\tau_{Kg(t)}^g - t|^\gamma \right\|_{L^p(P_x)}.
\end{aligned}
$$

As $\tau_{Kg(t)}^g$ is increasing with respect to $t$,

$$
\begin{aligned}
& \left\| 1_{\tau_{Kg(T)}^g} \le M \sup_{t \in [0,T]} |X_{\tau_{Kg(T)}^g} - X_t| \right\|_{L^p(P_x)} \\
& \le \left\| \sup_{s \ne t \le M} \left\{ \frac{|X_s - X_t|}{|s - t|^\gamma} \right\} \sup_{t \in [0,T]} |\tau_{Kg(t)}^g - t|^\gamma \right\|_{L^p(P_x)} \\
& = \left( E_x \left[ \sup_{s \ne t \le M} \left\{ \frac{|X_s - X_t|}{|s - t|^\gamma} \right\} \sup_{t \in [0,T]} |\tau_{Kg(t)}^g - t|^\gamma \right] \right)^{1/p}.
\end{aligned}
$$

From Hölder's inequality for $q \ge 1$ and $q/(q-1)$ conjugates exponents and (44), for every $\gamma \in (0, \frac{1}{2})$, there exists a constant $C_1 = C_1(M, \gamma, p(q-1)/q, x) > 0$ such that

$$
\begin{aligned}
& \left\| 1_{\tau_{Kg(T)}^g} \le M \sup_{t \in [0,T]} |X_{\tau_{Kg(T)}^g} - X_t| \right\|_{L^p(P_x)} \\
& \le \left( E_x \left[ \sup_{s \ne t \le M} \frac{|X_s - X_t|}{|s - t|^\gamma} \right]^{p(q-1)/q} \right)^{q/p(q-1)} \\
& \qquad \left( E_x \left[ \sup_{t \in [0,T]} |\tau_{Kg(t)}^g - t|^\gamma \right]^{pq} \right)^{1/pq} \\
& \le C_1^{q/p(q-1)} \left( E_x \left[ \sup_{t \in [0,T]} |\tau_{Kg(t)}^g - t|^{\gamma pq} \right] \right)^{1/pq}.
\end{aligned}
$$

Choosing $q = 2/\gamma p$, from (38), for any $\gamma \in (0, \frac{1}{2} \wedge \frac{2}{p})$

$$
\left\| 1_{\tau_{Kg(T)}^g} \le M \sup_{t \in [0,T]} |X_{\tau_{Kg(T)}^g} - X_t| \right\|_{L^p(P_x)} \le C_1^{1/p(1-\gamma p)} (2|g|_X (4T + |g|_X))^{\gamma/2}. \quad (48)
$$

For the second additive term of (47),

$$
E_x \left[ 1_{\tau_{Kg(T)}^g} > M \sup_{t \in [0,T]} |X_{\tau_{Kg(T)}^g} - X_t|^p \right] = \sum_{m=M}^{\infty} E_x \left[ 1_{\tau_{Kg(T)}^g} \in [m, m+1) \sup_{t \in [0,T]} |X_{\tau_{Kg(T)}^g} - X_t|^p \right]. \quad (49)
$$
For each term of the sum in (49), from Hölder's inequality,

$$
\begin{aligned}
& E_x \left[ \mathbb{1}_{\tau_{K\mathbf{g}(T)} \in [m,m+1]} \sup_{t \in [0,T]} |X_{\tau_{K\mathbf{g}(t)}} - X_t|^p \right] \\
& \le \left[ P_x(\tau_{K\mathbf{g}(T)} > m) \right]^{1/q'} \left[ E_x \left[ \mathbb{1}_{\tau_{K\mathbf{g}(T)} \in [m,m+1]} \sup_{t \in [0,T]} |X_{\tau_{K\mathbf{g}(t)}} - X_t|^p \right] \right]^{q'/(q'-1)} \\
& \le \left[ P_x(\tau_{K\mathbf{g}(T)} > m) \right]^{1/q'} \left[ E_x \left[ \mathbb{1}_{\tau_{K\mathbf{g}(T)} < m+1} \sup_{t \in [0,T]} |X_{\tau_{K\mathbf{g}(t)}} - X_t|^p \right] \right]^{q'/(q'-1)}.
\end{aligned}
$$

As each $m$ in the sum in (49) satisfies $m + 1 > M > T$, from Minkowski's inequality,

$$
\begin{aligned}
& E_x \left[ \mathbb{1}_{\tau_{K\mathbf{g}(T)} < m+1} \sup_{t \in [0,T]} |X_{\tau_{K\mathbf{g}(t)}} - X_t|^p \right] \\
& \le E_x \left[ \mathbb{1}_{\tau_{K\mathbf{g}(T)} < m+1} 2^{p-1} \sup_{t \in [0,T]} \left\{ |X_{\tau_{K\mathbf{g}(t)}} - x|^p + |X_t - x|^p \right\} \right] \\
& \le 2^p E_x \left[ \sup_{t \in [0,m+1]} |X_t - x|^p \right],
\end{aligned}
$$

which from Lemma 5.1 is bounded by $C_2[1 + |x|^p]e^{C_3(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 $\lambda > 0$ chosen such that $\alpha = \lambda|\mathbf{g}|_X < 1$,

$$
E_x \left[ \mathbb{1}_{\tau_{K\mathbf{g}(T)} \in [m,m+1]} \sup_{t \in [0,T]} |X_{\tau_{K\mathbf{g}(t)}} - X_t|^p \right] \le \left[ C_2 [1 + |x|^{pq'}] e^{C_3(m+1)} \right]^{1/q'} \left[ e^{-\lambda m} e^{\lambda \frac{T+|\mathbf{g}|_X}{1-\lambda|\mathbf{g}|_X}} \right]^{(q'-1)/q'}
$$

where $C_2(pq') \le 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$,

$$
\begin{aligned}
& E_x \left[ \mathbb{1}_{\tau_{K\mathbf{g}(T)} \in [m,m+1]} \sup_{t \in [0,T]} |X_{\tau_{K\mathbf{g}(t)}} - X_t|^p \right] \\
& \le C_4 \exp \left\{ \frac{1}{q'} \left[ C_3 m + (q' - 1) \left( \lambda \frac{T + |\mathbf{g}|_X}{1 - \lambda |\mathbf{g}|_X} - \lambda m \right) \right] \right\} \\
& = C_4 \exp \left\{ \frac{1}{q'} \left[ C_3 m + (q' - 1) \frac{\alpha}{|\mathbf{g}|_X} \left( \frac{T + |\mathbf{g}|_X}{1 - \alpha} - m \right) \right] \right\} \\
& = C_5 \exp \left\{ \frac{1}{q'} \left[ C_3 m + (q' - 1) \frac{\alpha}{|\mathbf{g}|_X} \left( \frac{T}{1 - \alpha} - m \right) \right] \right\},
\end{aligned}
$$
where $C_5 := C_4 \exp(\frac{q'-1}{q'} \frac{\alpha}{1-\alpha})$. If we choose $q' > 1$ such that$^7 A := C_3 - (q' - 1)\alpha/|g|_X < 0$,

$$
\begin{align*}
\sum_{m=M}^{\infty} \mathrm{E}_x \left[ \mathbb{1}_{\tau_{Kg}(T) \in [m,m+1]} \sup_{t \in [0,T]} |X_{\tau_{Kg}(t)} - X_t|^p \right] & \\
&\le C_5 \exp \left\{ \frac{(q'-1)\alpha T}{(1-\alpha)|g|_X} \right\} \sum_{m=M}^{\infty} (e^A)^m = C_5 \exp \left\{ \frac{(q'-1)\alpha T}{(1-\alpha)|g|_X} \right\} \frac{e^{AM}}{1-e^A} \\
&= C_5 \exp \left\{ \frac{(q'-1)\alpha T}{(1-\alpha)|g|_X} \right\} \frac{\exp \left\{ \left( C_3 - \frac{(q'-1)\alpha}{|g|_X} \right) M \right\}}{1 - \exp \left\{ C_3 - \frac{(q'-1)\alpha}{|g|_X} \right\}}.
\end{align*}
$$

If $M$ is chosen such that $M > T/(1 - \alpha)$, then

$$
\sum_{m=M}^{\infty} \mathrm{E}_x \left[ \mathbb{1}_{\tau_{Kg}(T) \in [m,m+1]} \sup_{t \in [0,T]} |X_{\tau_{Kg}(t)} - X_t|^p \right] \le \frac{C_5 e^{C_3 M}}{1-e^A} \exp \left\{ \frac{(q'-1)\alpha}{|g|_X} \left( \frac{T}{1-\alpha} - M \right) \right\}, (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),

$$
\left\| \sup_{t \in [0,T]} |X_{\tau_{Kg}(t)} - X_t| \right\|_{L^p(P_x)} \le C_1^{1/p(1-\gamma p)} (2|g|_X (4T + |g|_X))^{\gamma/2} + C^{(1)} e^{-C^{(2)}/|g|_X}.
$$

As $|g|_X^\gamma$ and $e^{-1/|g|_X}$ are both $O(|g|_X^{\gamma/2})$, there exists a constant $C > 0$, such that,

$$
\left\| \sup_{t \in [0, T]} |X_{\tau_{Kg}(t)} - X_t| \right\|_{L^p(P_x)} \le C |g|_X^{\gamma/2},
$$

which is (10).

<figure><img src="image_6.png" /></figure>

6 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 $(\mu, \sigma)$ two real-valued functions that satisfy Condition 1 and $X$ be the diffusion that solves

$$
dX_t = \mu(X_t) dt + \sigma(X_t) dB_t,
$$

<footnote>$^7$We remark that if this is satisfied for a grid $g$, then it is satisfied for all grids $g'$ such that $|g'|_X \le |g|_X$.</footnote>
where $B$ is a standard Brownian motion. Let also $\mathbb{I}$ be the state-space of $X$. A straightforward computation using Itô's formula gives us the infinitesimal generator of $X$,

$$ (L, \mathrm{dom}(L)) = \begin{cases} L f(x) = \mu(x)f'(x) + \frac{1}{2}\sigma^2(x)f''(x), & \forall f \in \mathrm{dom}(L), \\ \mathrm{dom}(L) = (f \in C_b(\mathbb{I}) : L f \in C(\mathbb{I})). \end{cases} \quad (51) $$

In particular, if $\mu$ and $\sigma$ are continuous, then, $\mathrm{dom}(L) = C^2(\mathbb{I})$. From Proposition 3.1 and (51) we can infer that,

$$ s'(x) = e^{-\int_a^y \frac{2\mu(\zeta)}{\sigma^2(\zeta)} d\zeta} \quad \text{and} \quad m(dx) = \frac{1}{s'(x)} \frac{2}{\sigma^2(x)} dx. \quad (52) $$

Thus, from Proposition 3.1,

$$ v_0(x) = \frac{\int_a^x e^{-\int_a^y \frac{2\mu(\zeta)}{\sigma^2(\zeta)} d\zeta} d\zeta}{\int_a^b e^{-\int_a^y \frac{2\mu(\zeta)}{\sigma^2(\zeta)} d\zeta} d\zeta}. \quad (53) $$

From (26) and (52),

$$ \begin{aligned} v_1(x) &= \int_a^b G_{a,b}(x, \zeta) v_0(\zeta) \frac{1}{s'(x)} \frac{2}{\sigma^2(x)} d\zeta, \\ \bar{v}_1(x) &= \int_a^b G_{a,b}(x, \zeta) (1 - v_0(\zeta)) \frac{1}{s'(x)} \frac{2}{\sigma^2(x)} d\zeta. \end{aligned} \quad (54) $$

The scale functions and the speed measures are defined up to a multiplicative constant: for $\alpha \in \mathbb{R}$ and $\lambda > 0$, the pairs $(s, m)$ and $(\alpha + \lambda s, \lambda^{-1}m)$ are associated to the same diffusion. In particular, as we are only concerned with points $x \in [a, b]$, we could use $v_0$ for the scale function. The speed measure shall be adapted accordingly. For $x, \zeta \in [a, b]$, the Green function in (1) takes the simpler form,

$$ G_{a,b}(x, \zeta) = v_0(x \wedge \zeta) (1 - v_0(x \vee \zeta)). $$

Thus, expressions (54) become

$$ \begin{aligned} v_1(x) &= \int_a^b v_0(x \wedge \zeta) (1 - v_0(x \vee \zeta)) \frac{v_0(\zeta)}{v'_0(\zeta)} \frac{2}{\sigma^2(\zeta)} d\zeta, \\ \bar{v}_1(x) &= \int_a^b v_0(x \wedge \zeta) (1 - v_0(x \vee \zeta)) \frac{1 - v_0(\zeta)}{v'_0(\zeta)} \frac{2}{\sigma^2(\zeta)} d\zeta. \end{aligned} \quad (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.
# 7 Numerical experiments

## 7.1 Approximation in distribution

### 7.1.1 Standard Brownian motion

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

$$s(x) = x \quad \text{and} \quad m(dx) = 2 \, dx.$$

Let $v_k(x) := E_x[\tau_{ab}^k 1_{\tau_b < \tau_a}]$ and $\bar{v}_k(x) := E_x[\tau_{ab}^k 1_{\tau_a < \tau_b}]$ for all $k \in \mathbb{N}_0$, where $c_x = (a, b)$. Then, from the definition of the scale function,

$$v_0(x) = \frac{x-a}{b-a}.$$

From Proposition 3.2,

$$v_1^{BM}(x) = (x-a)(b-x) \left( \frac{2}{3} \frac{(x-a)^2}{(b-a)^2} + \frac{b-x}{b-a} - \frac{2}{3} \frac{(b-x)^2}{(b-a)^2} \right),$$

$$\bar{v}_1^{BM}(x) = (x-a)(b-x) \left( \frac{2}{3} \frac{(b-x)^2}{(b-a)^2} + \frac{x-a}{b-a} - \frac{2}{3} \frac{(x-a)^2}{(b-a)^2} \right),$$

where by $BM$ we mean these quantities are associated with the Brownian motion. Thus, from (6), we have all the necessary quantities we need to implement the algorithm.

### 7.1.2 Sticky Brownian motion

The Brownian motion with a sticky point at 0 of stickiness $\rho > 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 \quad \text{and} \quad m(dx) = 2 \, dx + \rho \delta_0(dx).$$

From the definition of the scale function,

$$v_0(x) = \frac{x-a}{b-a}.$$

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

$$v_1(x) = v_1^{BM}(x) + \rho 1_{0 \in (a,b)} G_{(a,b)}(x, 0) v_0(0),$$

$$\bar{v}_1(x) = \bar{v}_1^{BM}(x) + \rho 1_{0 \in (a,b)} G_{(a,b)}(x, 0) (1 - v_0(0)),$$

where $v_1^{BM}(x)$ and $\bar{v}_1^{BM}(x)$ 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]).
<figure><img src="image_7.png" /><figcaption>Figure 3: (a): Histogram of simulated values at $T = 1$ of a sticky Brownian motion of parameter $\rho = 0.7$ with initial value $x_0 = 0$ using Algorithm 1 with the tuned grid (14) of size-criteria $h = 0.01$.</figcaption></figure>

<caption>(b): Histogram of simulated values at $T = 1$ of a skew Brownian motion of parameter $\beta = 0.9$ with initial value $x_0 = 0$ using Algorithm 1 with a uniform grid of step-size $h = 0.01$.</caption>

### 7.1.3 Skew Brownian motion

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

$$s(x) = \begin{cases} \frac{x}{\beta}, & x \ge 0, \\ \frac{x}{1-\beta}, & x \le 0 \end{cases} \quad \text{and} \quad m(dx) = \begin{cases} 2\beta \, dx, & x > 0, \\ 2(1-\beta) \, dx, & x < 0. \end{cases}$$

From the definition of the scale function,

$$v_0(x) = \frac{s(x) - s(a)}{s(b) - s(a)}.$$

As the $\beta$ and $(1-\beta)$ 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_1^{BM}(x) \text{ and } \bar{v}_1(x) = \bar{v}_1^{BM}(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].
### 7.1.4 Ornstein-Uhlenbeck process

Let $X$ be the Ornstein-Uhlenbeck with mean reversion force $\theta > 0$, long-term mean $\mu$ and diffusion parameter $\sigma > 0$ define through its path-wise description

$$dX_t = \theta(\mu - X_t) dt + \sigma dB_t, \quad X_0 \in \mathbb{R},$$

where $B$ is a standard Brownian motion. We will see that $v_0$ can be expressed in terms of the Gaussian imaginary error function as $s(x) = \text{erfi}(\sqrt{\frac{\theta}{\sigma^2}}(\mu - x))$,

$$v_0(x) = \frac{\int_a^x e^{-\int_a^y \frac{2\mu(\zeta)}{\sigma^2(\zeta)} d\zeta} dy}{\int_a^b e^{-\int_a^y \frac{2\mu(\zeta)}{\sigma^2(\zeta)} d\zeta} dy} = \frac{\text{erfi}(\sqrt{\frac{\theta}{\sigma^2}}(\mu - x)) - \text{erfi}(\sqrt{\frac{\theta}{\sigma^2}}(\mu - a))}{\text{erfi}(\sqrt{\frac{\theta}{\sigma^2}}(\mu - b)) - \text{erfi}(\sqrt{\frac{\theta}{\sigma^2}}(\mu - a))},$$

where $\text{erfi}(x) = \frac{2}{\pi} \int_0^x e^{t^2} dt = \frac{2}{\pi} e^{x^2} D_+(x)$, with $D_+(x)$ being the Dawson function$^8$. Thus, for an Ornstein-Uhlenbeck process (54) becomes

$$v_1(x) = \int_a^b v_0(x \wedge \zeta) (1 - v_0(x \vee \zeta)) v_0(\zeta) c \exp\left(\frac{\theta(\zeta - \mu)^2}{\sigma^2}\right) \frac{2}{\sigma^2} dx,$$

$$\bar{v}_1(x) = \int_a^b v_0(x \wedge \zeta) (1 - v_0(x \vee \zeta)) (1 - v_0(\zeta)) c \exp\left(\frac{\theta(\zeta - \mu)^2}{\sigma^2}\right) \frac{2}{\sigma^2} dx,$$

$$\text{where } c = \text{erfi}\left(\sqrt{\frac{\theta}{\sigma^2}}(\mu - b)\right) - \text{erfi}\left(\sqrt{\frac{\theta}{\sigma^2}}(\mu - a)\right).$$

### 7.1.5 Cox-Ingersoll-Ross process

The Cox-Ingersoll-Ross process or CIR process [11], introduced first by W. Feller [16], is the diffusion that solves the SDE

$$dX_t = \theta(\mu - X_t) dt + \sigma \sqrt{X_t} dB_t, \quad X_0 > 0,$$

where $B_t$ is a standard Brownian motion. The parameter $\theta > 0$ expresses its mean reversion speed, $\mu$ its long term speed and $\sigma > 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\theta\mu > \sigma^2$. A large body of work have been devoted to the simulation of the CIR and related process, see e.g. [1].

From (53) and (55),

$$v_0(x) = \frac{\int_a^x y^{-\frac{2\theta\mu}{\sigma^2}} e^{\frac{2\theta}{\sigma^2}y} dy}{\int_a^b y^{-\frac{2\theta\mu}{\sigma^2}} e^{\frac{2\theta}{\sigma^2}y} dy}$$

<footnote>$^8$We use the latter representation in our numerical results.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 4: Histogram of simulated values at $T = 1$ of a CIR process of $(\theta, \mu, \sigma) = (5, 5, 1)$ with initial value $x_0 = 1$ using Algorithm 1 with:</figcaption></figure>

(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).

and

$$v_1(x) = \int_a^b v_0(x \wedge \zeta) (1 - v_0(x \vee \zeta)) \frac{v_0(\zeta)}{v'_0(\zeta)} \frac{2}{\sigma^2 \zeta} d\zeta,$$

$$\bar{v}_1(x) = \int_a^b v_0(x \wedge \zeta) (1 - v_0(x \vee \zeta)) \frac{1 - v_0(\zeta)}{v'_0(\zeta)} \frac{2}{\sigma^2 \zeta} d\zeta.$$

As the scale function yields no satisfactory closed formula, we perform a numerical approximation of both $v_0$ and the couple $(v_1, \bar{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 and 5 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, \bar{v}_1)$. We also superpose on the same figures the corresponding transition kernels.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Figure 5: Histogram of simulated values at $T = 1$ of a CIR process of $(\theta, \mu, \sigma) = (1, 2, 1)$ with initial value $x_0 = 5$ using Algorithm 1 with:</figcaption></figure>

(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 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).
<figure><img src="image_13.png" /><figcaption>Figure 6: Histogram of simulated values at $T = 1$ of the Skew Bessel process of parameters $(\delta, \beta) = (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.</figcaption></figure>

### 7.1.6 Skew or reflected Bessel process

The Bessel process is defined as the solution of the SDE

$$dX_t = \frac{\delta - 1}{2X_t} dt + dB_t,$$

where $B$ is a standard Brownian motion. If $\delta \in (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 $\delta \in (0, 2)$ and skew $\beta \in (0, 1)$ is the diffusion process with state-space $\mathbb{R}$ defined through $s$ and $m$, where (see [2])

$$s(x) = \begin{cases} \frac{1}{\beta} \frac{x^{2-\delta}}{2-\delta}, & x > 0, \\ -\frac{1}{1-\beta} \frac{|x|^{2-\delta}}{2-\delta}, & x \le 0, \end{cases} \qquad m(dx) = \begin{cases} 2\beta x^{\delta-1} dx, & x > 0, \\ 2(1-\beta)|x|^{\delta-1} dx, & x \le 0. \end{cases}$$

This yields the following expressions for the quantities we need to compute in order to
<figure><img src="image_14.png" /><figcaption>Figure 7: (a): histogram of simulated values at $T = 1$ of a Reflected Bessel process of parameter $\delta = 1.1$ with initial value $x_0 = 5$ using Algorithm 1 with a uniform grid of step-size $h = 0.01$.</figcaption></figure>

(b): histogram of simulated values at $T = 1$ of a Bessel process of parameter $\delta = 1.8$ with initial value $x_0 = 5$ using Algorithm 1 with a uniform grid of step-size $h = 0.01$. The quantities (56) were approximated using 100-step Riemann approximations.

implement the algorithm,

$$ 
\begin{aligned}
v_0(x) &= \frac{s(x) - s(a)}{s(b) - s(a)}, \\
v_1(x) &= \int_a^b G_{a,b}(x, \zeta)v_0(\zeta)2|\zeta|^{\delta-1} d\zeta, \\
\bar{v}_1(x) &= \int_a^b G_{a,b}(x, \zeta)(1 - v_0(\zeta))2|\zeta|^{\delta-1} d\zeta,
\end{aligned} 
\quad (56) $$

where $G_{a,b}(x, \zeta)$ is defined in (1).

In Figures 6 and 7 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 [2].

## 7.2 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
orders of convergence of localized path-sensitive functionals. An example is the stickiness parameter estimation in [4]. It is shown there that, if $X$ is a sticky Brownian motion of stickiness $\rho > 0$, $\alpha \in (0, 1/2)$ and $g$ is an integrable function vanishing in the vicinity of 0, the statistics

$$ \hat{L}_n[g, \alpha] := \frac{n^\alpha}{n} \sum_{i=1}^{[nt]} g(n^\alpha X_{\frac{i-1}{n}}) \quad (57) $$

and

$$ \hat{\rho}_n[g, \alpha] := 2 \frac{\lambda(g)}{\hat{L}_n[g, \alpha]} \frac{1}{n} \sum_{i=1}^{[nt]} \mathbb{1}\{X_{\frac{i-1}{n}}=0\} \quad (58) $$

converge in probability to $L_t^0(X)$ and $\rho$ respectively. This qualifies $\hat{\rho}_n[g, \alpha]$ 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 $\alpha \in (0, 1)$.

* The convergence speed of (58) seems to increase in terms of $\alpha \in (0, 1)$.

* The higher the $\alpha$, 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 $\alpha$ 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\} \cup \{\pm(h^2/\rho + kh); k \in \mathbb{N}\}, \quad g_1(h) = \{\pm x_k(h); k \ge 0\}, \quad (59) $$

where $\{x_j(h)\}_{j \ge 0}$ is defined recursively by

$$ x_0 = 0, \quad x_j = x_{j-1} + \left( \frac{h^2}{\rho} \frac{1}{x_{j-1} + 1} + h \left( 1 - \frac{1}{x_{j-1} + 1} \right) \right) \mathbb{1}_{\{x_{j-1} < 1\}} + h \mathbb{1}_{\{x_{j-1} \ge 1\}}. \quad (60) $$

We observe that

$$ |g_1| = |g_0| \quad \text{and} \quad |g_1|_X = |g_0|_X. \quad (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 \le N_{MC}$ we simulate the path of an approximation process $X_t^j$. To assess the quality of each Monte Carlo estimation we use the following metrics:
<figure><img src="image_15.png" /><figcaption>Figure 8: Stickiness parameter estimations histograms using (58) with $n = 100 000$ and $\alpha = 0.3$ (true value $\rho = 1$).</figcaption></figure>

<table><caption>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.</caption><thead><tr><th>α</th><th>n</th><th>ρ̂<sub>MC</sub></th><th>Ŝ<sup>2</sup><sub>MC</sub></th><th>σ̂<sub>MC</sub></th><th>acĉ</th><th>rej/N<sub>MC</sub></th></tr></thead><tbody><tr><td>0.3</td><td>100000</td><td>1.334</td><td>0.740</td><td>0.860</td><td>0.125</td><td>0/2000</td></tr><tr><td>0.4</td><td>100000</td><td>1.061</td><td>0.142</td><td>0.378</td><td>0.046</td><td>0/2000</td></tr><tr><td>0.5</td><td>100000</td><td>1.314</td><td>0.162</td><td>0.403</td><td>0.011</td><td>0/2000</td></tr><tr><td>0.55</td><td>100000</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2000/2000</td></tr><tr><td>0.6</td><td>100000</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2000/2000</td></tr><tr><td>0.65</td><td>100000</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2000/2000</td></tr></tbody></table>

* $(\hat{\rho}_{MC}, \hat{S}_{MC}^2, \hat{\sigma}_{MC})$: Monte Carlo estimation, variance and standard deviation of the stickiness parameter estimations,

* $\hat{acc}$: average number of path-values observed by $g$, i.e

$$ \hat{acc} = \frac{1}{N_{MC}} \sum_{j=1}^{N_{MC}} \frac{1}{n} \sum_{i=1}^{n} \mathbb{1}_{g(X_{(i-1)/n}^{j}) \neq 0} \quad (62) $$
<figure><img src="image_16.png" /><figcaption>Figure 9: Stickiness parameter estimations histograms using (58) with $n = 100 000$ and $\alpha = 0.55$ (true value $\rho = 1$).</figcaption></figure>

<table><caption>Table 2: Stickiness parameter estimations using the grid $g_1$ for $h = 0.01$. Computation time (single-core): 67 seconds.</caption><thead><tr><th>α</th><th>n</th><th>ρ̂<sub>MC</sub></th><th>Ŝ<sup>2</sup><sub>MC</sub></th><th>σ̂<sub>MC</sub></th><th>acĉ</th><th>rej/N<sub>MC</sub></th></tr></thead><tbody><tr><td>0.3</td><td>100000</td><td>1.252</td><td>0.629</td><td>0.793</td><td>0.1272</td><td>0/2000</td></tr><tr><td>0.4</td><td>100000</td><td>1.068</td><td>0.187</td><td>0.432</td><td>0.0437</td><td>0/2000</td></tr><tr><td>0.5</td><td>100000</td><td>1.023</td><td>0.072</td><td>0.268</td><td>0.0140</td><td>0/2000</td></tr><tr><td>0.55</td><td>100000</td><td>1.018</td><td>0.067</td><td>0.259</td><td>0.0079</td><td>0/2000</td></tr><tr><td>0.6</td><td>100000</td><td>1.012</td><td>0.031</td><td>0.177</td><td>0.0044</td><td>0/2000</td></tr><tr><td>0.65</td><td>100000</td><td>1.015</td><td>0.026</td><td>0.163</td><td>0.0025</td><td>2/2000</td></tr></tbody></table>

* rej: percentage of trajectories where the local time estimation equals 0, i.e.,

$$ \text{rej} = \#\left\{j \le N_{\text{MC}} : \frac{n^\alpha}{n} \sum_{i=1}^{[nt]} g(n^\alpha X_{\frac{i-1}{n}}^j) = 0\right\}. \quad (63) $$

We use the test function $g(x) = 1_{1<x<5}/8$ which satisfies the conditions of Corollary 1.3 of [4]. Within each table (Tables 1-2) we use the same simulated STMCA trajectories of the sticky Brownian motion of parameter $\rho = 1$.

We observe that the usage of grid $g_1$ yields far superior results than $g_0$ (see Tables 1-2).
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 $\alpha$ and we can achieve higher orders of convergence (see Figures 8-9, and Tables 1-2). 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 $\mathbb{I}$ an interval of $\mathbb{R}$ defined on a family of probability spaces $\mathcal{P} = (\Omega, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P}_x)_{x \in \mathbb{I}}$ such that for every $x \in \mathbb{I}$, $\mathbb{P}_x(X_0 = x) = 1$. Let also $g$ be a covering grid of $\mathbb{I}$, $\{\tau_j\}_{j \in \mathbb{Z}_+}$ the sequence of embedding times of $X$ in $g$ defined in (21) and $\mathcal{B}$ the sigma-algebra defined by $\mathcal{B} = \sigma\{X_{\tau_j}; j \in \mathbb{N}_0\}$. Then, for any measurable path-functional $F: C^0(\mathbb{R}_+, \mathbb{I}) \to \mathbb{R}$ and $j \ge 1$ and $x \in \mathbb{I}$,

$$ \mathrm{E}_x(F((X_t^{\tau_{j-1}, \tau_j})_{t \ge 0}) | \mathcal{B}) = \mathrm{E}_x(F((X_t^{\tau_{j-1}, \tau_j})_{t \ge 0}) | X_{\tau_{j-1}}, X_{\tau_j}), \quad (64) $$

where $X^{\tau_{j-1}, \tau_j}$ is the process defined for every $t \ge 0$ by $X_t^{\tau_{j-1}, \tau_j} = X_{(\tau_{j-1}+t) \wedge \tau_j}$.

*Proof.* Let us fix $x_1, x_2, \dots$ be a sequence of points in the grid. Let us define

$$ Q(x; x_1, x_2, \dots) := \mathrm{P}_x \left( X_{\tau_1} = x_1, X_{\tau_2} = x_2, \dots \right) $$

By the strong Markov property,

$$ \begin{aligned} \mathrm{P}_x \left( X_{\tau_{i+1}} = x_{i+1}, X_{\tau_{i+2}} = x_{i+2}, \dots \middle| \mathcal{F}_{\tau_i} \right) &= \mathrm{P}_{X_{\tau_i}} \left( X_{\tau_1} = x_{i+1}, X_{\tau_{i+2}} = x_{i+2} \right) \\ &= Q(X_{\tau_i}; x_{i+1}, x_{i+2}, \dots). \end{aligned} $$

Using the strong Markov property twice, first by conditioning first with respect to $\mathcal{F}_{\tau_{j+1}}$ and then with respect to $\mathcal{F}_{\tau_j}$,

$$ \begin{aligned} \mathrm{E}_x \left( \mathbb{1}_{X_{\tau_1} = x_1, \dots, X_{\tau_j} = x_j} F((X_t^{\tau_{j-1}, \tau_j})_{t \ge 0}) \mathbb{1}_{X_{\tau_j} = x_j, X_{\tau_{j+1}} = x_{j+1}, \dots} \right) \\ &= \mathrm{E}_x \left( \mathbb{1}_{X_{\tau_1} = x_1, \dots, X_{\tau_j} = x_j} F((X_t^{\tau_{j-1}, \tau_j})_{t \ge 0}) \mathbb{1}_{X_{\tau_j} = x_j} Q(X_{\tau_j}) \right) \\ &= \mathrm{E}_x \left( \mathbb{1}_{X_{\tau_1} = x_1, \dots, X_{\tau_{j-1}} = x_{j-1}} (X_t^{\tau_{j-1}, \tau_j})_{t \ge 0} \right) Q(x_j) = \mathrm{E}_x \left( \mathbb{1}_{X_{\tau_1} = x_1, \dots, X_{\tau_j} = x_{j-1}} \right) R(x_{j-1}, x_j) Q(x_j) \end{aligned} $$

with $R(x, y) := \mathrm{P}_x(F((X_t^{0, \tau_1})_{t \ge 0}) \mathbb{1}_{X_{\tau_1} = y})$. Therefore, using the definition of the conditional expectation,

$$ \begin{aligned} \mathrm{E}_x \left( F((X_t^{\tau_{j-1}, \tau_j})_{t \ge 0}) \middle| X_{\tau_1} = x_1, X_{\tau_2} = x_2, \dots \right) \\ &= \frac{R(x_{j-1}, x_j)}{\mathrm{P}_{x_{j-1}}(X_{\tau_1} = x_j)} = \mathrm{E}_{x_{j-1}} \left( F((X_t^{0, \tau_1})_{t \ge 0}) \middle| X_{\tau_1} = x_j \right). \end{aligned} $$

This is sufficient to prove the result. $\square$
**References**

[1] Aurélien Alfonsi. On the discretization schemes for the CIR (and Bessel squared) processes. *Monte Carlo Methods Appl.*, 11(4):355-384, 2005.

[2] Larbi Alili and Andrew Aylwin. On the semi-group of a scaled skew Bessel process. *Statist. Probab. Lett.*, 145:96-102, 2019.

[3] Madjid Amir. Sticky Brownian motion as the strong limit of a sequence of random walks. *Stochastic Process. Appl.*, 39(2):221-237, 1991.

[4] Alexis Anagnostakis. Functional convergence to the local time of a sticky diffusion, 2022.

[5] Stefan Ankirchner, Thomas Kruse, and Mikhail Urusov. Wasserstein convergence rates for coin tossing approximations of continuous markov processes, 03 2019. Preprint arXiv:1903.07880.

[6] Stefan Ankirchner, Thomas Kruse, and Mikhail Urusov. A functional limit theorem for coin tossing Markov chains. *Ann. Inst. Henri Poincaré Probab. Stat.*, 56(4):2996-3019, 2020.

[7] Richard F. Bass. A stochastic differential equation with a sticky point. *Electron. J. Probab.*, 19:no. 32, 22, 2014.

[8] Andrei N. Borodin and Paavo Salminen. *Handbook of Brownian motion—facts and formulae*. Probability and its Applications. Birkhäuser Verlag, Basel, 1996.

[9] Haim Brezis. *Functional analysis, Sobolev spaces and partial differential equations*. Universitext. Springer, New York, 2011.

[10] K. L. Chung and G. A. Hunt. On the zero $\sum_{1}^{n} \pm 1$. *Ann. Math. (2)*, 50:385-400, 1949.

[11] John C. Cox, Jonathan E. Ingersoll, Jr., and Stephen A. Ross. A theory of the term structure of interest rates. *Econometrica*, 53(2):385-407, 1985.

[12] Monroe D. Donsker. An invariance principle for certain probability limit theorems. *Mem. Amer. Math. Soc.*, 6:12, 1951.

[13] Andreas Eberle and Raphael Zimmer. Sticky couplings of multidimensional diffusions with different drifts. *Ann. Inst. Henri Poincaré Probab. Stat.*, 55(4):2370-2394, 2019.

[14] Hans-Jürgen Engelbert and Goran Peskir. Stochastic differential equations for sticky Brownian motion. *Stochastics*, 86(6):993-1021, 2014.
[15] Pierre Étoré and Antoine Lejay. A Donsker theorem to simulate one-dimensional processes with measurable coefficients. *ESAIM Probab. Stat.*, 11:301–326, 2007.

[16] William Feller. Two singular diffusion problems. *Ann. of Math. (2)*, 54:173–182, 1951.

[17] William Feller. The parabolic differential equations and the associated semi-groups of transformations. *Ann. of Math. (2)*, 55:468–519, 1952.

[18] William Feller. Generalized second order differential operators and their lateral conditions. *Ill. J. Math.*, 1:459–504, 1957.

[19] Anna Ferrer-Admetlla, Christoph Leuenberger, Jeffrey Jensen, and Daniel Wegmann. An approximate markov model for the wright-fisher diffusion and its application to time series data. *Genetics*, 203, 04 2016.

[20] Noufel Frikha. On the weak approximation of a skew diffusion by an Euler-type scheme. *Bernoulli*, 24(3):1653–1691, 2018.

[21] Hatem Hajri, Mine Caglar, and Marc Arnaudon. Application of stochastic flows to the sticky Brownian motion equation. *Electron. Commun. Probab.*, 22:Paper No. 3, 10, 2017.

[22] Martin Hutzenthaler and Arnulf Jentzen. *Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients*, volume 1112 of *Mem. Am. Math. Soc*. Providence, RI: American Mathematical Society (AMS), 2015.

[23] Kiyosi Itō. *Essentials of stochastic processes*, volume 231. American Mathematical Soc., 2006.

[24] Antoine Lejay. Monte Carlo methods for fissured porous media: a gridless approach. *Monte Carlo Methods Appl.*, 10(3-4):385–392, 2004.

[25] Antoine Lejay. On the constructions of the skew Brownian motion. *Probab. Surv.*, 3: 413–466, 2006.

[26] Christian Meier, Lingfei Li, and Zhang Gongqiu. Markov chain approximation of one-dimensional sticky diffusions. *10 2019*. Preprint arXiv:1910.14282.

[27] Yutian Nie and Vadim Linetsky. Sticky reflecting ornstein-uhlenbeck diffusions and the vasicek interest rate model with the sticky zero lower bound. *Stochastic Models*, 0(0): 1–19, 2019.

[28] Tomasz Piskorski and Mark M. Westerfield. Optimal dynamic contracts with moral hazard and costly monitoring. *J. Econom. Theory*, 166:242–281, 2016.

[29] Daniel Revuz and Marc Yor. *Continuous martingales and Brownian motion*, volume 293 of *Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]*. Springer-Verlag, Berlin, third edition, 1999.
[30] L. C. G. Rogers and David Williams. *Diffusions, Markov processes, and martingales. Vol. 2*. Cambridge Mathematical Library. Cambridge University Press, Cambridge, 2000. Itô calculus, Reprint of the second (1994) edition.

[31] Cédric Villani. *Optimal transport. Old and new*, volume 338 of *Grundlehren Math. Wiss*. Berlin: Springer, 2009.

[32] John Y. Zhu. Optimal contracts with shirking. *Rev. Econ. Stud.*, 80(2):812–839, 2013.