Abstract

A sticky diffusion is a process that can stick to and detach from a lower-dimensional boundary. A challenge in simulating such a process is in capturing the change in dimension in a dynamically consistent way. We introduce a numerical algorithm to simulate a one-dimensional sticky diffusion, which sticks to and detaches from a point. Our method is a simple modification of the standard Euler-Maruyama scheme, which chooses with some probability between a reflected Euler-Maruyama update and a jump to the sticky point. We show how to choose this probability to be consistent with the generator of the desired dynamics, and we prove that our scheme converges weakly to a sticky diffusion with order 1.

1 Introduction↩︎

Sticky diffusions are diffusion processes that can spend a nonzero amount of time on a lower-dimensional boundary. The simplest example is a sticky Brownian motion, which behaves as a Brownian motion away from the origin, but is slowed down at the origin in such a way that it spends a positive measure of time there. A more general sticky diffusion solves a stochastic differential equation (SDE) in an ambient space, but it can stick to a boundary, where it can move according to an SDE that can be different from that in the ambient space.

Sticky diffusions have been studied mathematically since the pioneering work of Feller in the 1950s [1], [2], who aimed to classify all possible boundary behaviours of a one-dimension diffusion process. His work spurred subsequent developments in both probability and PDE theory [3]–[7]; see also the mathematical history by Peskir [8].

Somewhat neglected for a few decades, sticky diffusions have regained attention recently because they have been found to be useful models in a wide variety of applications. In finance, they are used to model financial variables whose values may linger near certain threshholds [9], such as Japanese interest rates during 1995-2004 [10], US interest rates following the 2008 economic crisis [11], and US Treasury bill yields in the 1930s [12], which all spent long times hovering near zero. In physics and chemistry, sticky diffusion processes can model the dynamics of small particles in a fluid with short-ranged adhesive interactions, where particles may form transient bonds [13]–[18]. Biology has found a range of applications of sticky diffusions, including in cell biology, where they have been used to study the capture and release of molecules diffusing on a cell surface [19], [20]; epidemiology, where the evolution of pathogen load can have a sticky boundary at zero to capture the nonzero probability of an individual being not infected [21], [22]; and ecology, where they have recently been evoked as a mechanism to sustain the distribution of biomass and in particular the abundance of rare species [23]. Operations researchers have used sticky diffusions in queuing models and storage processes [24], and have recently generalized such ideas to networks, modelling e.g. traffic flow or data transmission, where sticky boundary conditions at the vertices of a network can capture congestion effects [25]. Sticky diffusions have also found a range of applications in mathematical coupling techniques, used to prove convergence rates of sampling algorithms and solutions to stochastic differential equations [26]–[28].

Given their burgeoning areas of application, it is natural to ask how to simulate a sticky diffusion. Yet, doing so is a challenge because of the need to capture a singular measure on the boundary. Standard methods for simulating SDEs, such as the Euler-Maruyama method [29], are not straightforward to adapt, as they take random steps in space so they will never lie exactly on the boundary. Inspired by interacting particle systems, one can approximate a sticky boundary condition by applying a strong, short-ranged force to keep a process near the boundary [30]–[32], but such a method is not efficient as resolving the force requires time steps much smaller than the timescales of interest.

An alternative idea, and the one which has seen the most development recently, is to take random steps in time, and move on a set of pre-defined points in space, leading to a continuous-time Markov chain approximation [33]. This idea was first proposed to simulate a sticky Brownian motion in [32], and it was later extended to more general one-dimensional [34] and multidimensional [35] sticky diffusions with flat boundaries. While this method works effectively for low-dimensional problems, there remain a few drawbacks, particularly for higher-dimensional problems. One is that to choose a direction to move in, the method requires evaluating a number of rate constants which scales at least linearly but often quadratically or higher with dimension. Two, it can be hard to find a discretization of space which leads to a valid continuous-time Markov chain. The authors in [35] propose to compute the discrete moves using the eigenvectors of the diffusion covariance matrix, but in high-dimensional problems such an eigendecomposition can become infeasible. Three, it is not clear how feasible such a method will be when applied to problems with curved boundaries, such as arise when considering distance constraints between point particles. Finally, such a method is significantly different from how diffusion processes and in particular interacting particle systems are normally simulated, making it hard to incorporate into existing simulation methodologies [36]–[38].

A recent study made steps toward developing a perhaps more natural discretization based on the Euler-Maruyama method [39]. This method takes a (discrete) step based on a discrete version of an Euler-Maruyama proposal step, but if a boundary is crossed then it reflects the step as for a reflecting diffusion [40] while simultaneously adjusting the time step to account for the slowdown near the sticky boundary. This method can simulate a sticky diffusion where the drift and diffusion are the same in the interior and on the boundary, because such a process can be constructed from a time-change of a reflected diffusion process. Unfortunately, it does not seem easy to generalize to allow for different drift and diffusion on the boundary, because such a process cannot be constructed as an aforementioned time-changed process. It also seems a disadvantage that the method never moves exactly to the boundary, as it fails to leverage the many methods developed for simulating processes with constraints [41]–[45].

Our contribution. In this paper we make an important step toward developing a general method to simulate arbitrary sticky diffusions, by proposing a modified Euler-Maruyama method to simulate a one-dimensional sticky reflecting diffusion process, which can jump exactly to and from a boundary. The method is only a small modification of the standard reflected Euler-Maruyama scheme, but it allows the process to lie exactly on the boundary, without requiring it to be a time-change of a reflecting diffusion. The method works as follows: when the process is far enough away from its sticky point (which we take to be at \(x=0\)), the algorithm proposes updates as for the regular Euler-Maruyama (EM) algorithm with reflection. When the process is close enough to 0, it performs a reflected EM update with some probability \(\lambda(x)\) depending on its current position \(x\), and it moves exactly to 0 with some probability \(1-\lambda(x)\). The function \(\lambda\) is chosen to make the transition operator over one step of the discrete scheme, a good approximation for the transition operator over a timestep \(h\) of the continuous dynamics.

We prove our method converges weakly with order \(O(h)\). The global error arises from a combination of the local error, i.e. the error over one timestep, and the number of points in each state: near the boundary, or away from the boundary. Specifically, away from the boundary, the local error of the EM update is \(O(h^2)\); the number of points in such a state is \(O(h^{-1})\) leading to \(O(h)\) global convergence. Near the boundary, the local error of our scheme is \(O(h^{3/2})\), however the number of points near the boundary is \(O(h^{-1/2})\) contributing to the same global error. The most technical part of the proof is bounding the number of points near the boundary.

While we develop the method in detail here for a one-dimensional process, we anticipate the ideas can be extended to a multi-dimensional process by following the heuristics and rigorous convergence tools that we describe in some detail here. An advantage of such an extension compared with other methods, is that we propose standard moves, with the only modification being that we choose from the discrete set of moves with probabilities consistent with a sticky diffusion process. Thus, we expect such a method to be easy to incorporate into existing simulation methodologies, which allow for such discrete choices of dynamics, for example when incorporating chemical reactions.

Overview of the paper. In Section 2 we briefly recall the properties of one-dimensional sticky diffusion processes. We also outline some formal tools for deriving their associated forward and backward equations and for understanding their properties. In Section 3 we present our method for a sticky Brownian motion (SBM). This section gives an overview of the method, a formal overview of how to derive the function \(\lambda\), and some numerical experiments, before presenting our proofs of convergence. We focus on an SBM first, because the ideas are easier to follow, and because the more general case reuses some of the results developed for an SBM. A general one-dimensional sticky diffusion is considered in Section 4, following the same structure as for SBM. The main convergence results and the lemmas used to prove them are presented at the end of the section, but the proofs are saved for the Supplement as they are technically involved but use standard techniques.

2 A brief introduction to one-dimensional sticky diffusion processes↩︎

In this section we give an overview of the key equations characterizing one-dimensional sticky diffusions. Our overview will be mainly formal, with the goal of providing accessible ways of writing down the equations characterizing with such diffusions.

2.1 Stochastic differential equation↩︎

A sticky diffusion process is the weak solution to a particular stochastic differential equation (SDE). Though not the main tool we use to characterize such processes, we present the SDE here for completeness.

We consider a process which, on domain \(D=(0,\infty)\), behaves as a diffusion process with drift \(b(x)\) and diffusion \(\sigma(x)\). At the boundary, \(\partial D = \{0\}\), the process is delayed before it is reflected, in such a way that it cumulatively spends a finite, nonzero amount of time there. The process can be constructed as the solution to the following system of stochastic differential equations [46]: \[\tag{1} \begin{align} dX_t &= b(X_t)1_D(X_t)dt + \sigma(X_t)1_D(X_t) dW_t + dL_t,\tag{2}\\ \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)dt &= \kappa dL_t. \tag{3} \end{align}\] Here \(1_A(\cdot)\) is the indicator function for set \(A\), \(W_t\) is a one-dimensional Brownian motion, and \(L_t\) is the local time process, a nondecreasing process which increases only when \(X_t\) is on the boundary, as \[L_t = \int_0^t 1_{\partial D}(X_s)dL_s.\] The constant \(\kappa>0\) characterizes how “sticky” the origin is. Larger \(\kappa\) implies the origin is stickier, while \(\kappa=0\) is not at all sticky and corresponds to a purely reflected process.

The existence and uniqueness of weak solutions to 1 is shown in [46], where a sufficient condition is shown to be that \(\sigma,b\) are bounded and Lipschitz continuous on \((0,\infty)\) with \(\sigma(0)^2>0\). (This reference shows existence and uniqueness for a general \(d\)-dimensional sticky diffusion on a half space but we specialize the equation and remarks here to a one-dimensional process.)

2.2 Backward equation↩︎

The key object we will work with to show weak convergence of our numerical scheme is the backward equation associated with 1 . We start by recalling the generator for 1 is the differential operator \(\mathcal{L}\) given by \[\label{eq:L} \mathcal{L}f(x) = b(x)\partial_x f + a(x) \partial_{xx}f.\tag{4}\] As usual \(a(x) = \frac{1}{2}\sigma^2(x)\). The sticky boundary implies a restriction on the domain of the generator, formulated as the sticky boundary condition [46] \[\label{eq:BC} a\partial_x f = \kappa \mathcal{L}f \qquad \text{at } x=0.\tag{5}\] The backward equation is thus characterized by the following.

Proposition 1. Let \(\phi:[0,\infty)\to \mathbb{R}\) be a bounded function, let \(T>0\), and suppose the backward equation \[\label{eq:back} \begin{align} \partial_s u + \mathcal{L} u &= 0, & (x\geq 0) \label{back95sde}\\ a\partial_x u &= \kappa \mathcal{L}u & (x=0) \label{back95bc}\\ u(T,x) &= \phi(x). \label{back95terminal} \end{align}\] {#eq: sublabel=eq:eq:back,eq:back95sde,eq:back95bc,eq:back95terminal} has a solution \(u:[0,T]\times[0,\infty)\to \mathbb{R}\) which is such that \(u,\partial_tu,\partial_xu,\partial_{xx}u\) are continuous and bounded on \([0,T]\times[0,\infty)\), and \(\sigma\partial_xu\) is bounded on that same domain. Then \(u\) has the representation \(u(t,x) = \mathbb{E}[\phi(X_T)|X_t=x]\), where \((X_t)_{t\geq 0}\) is a weak solution to 1 .

Proof. Let \(u\) solve ?? . Using Itô’s formula, rewriting \(1_D\mathcal{L} u = \mathcal{L} u-1_{\partial_D}\mathcal{L} u\), and using 3 , we evaluate: \[\begin{gather} u(T,X_T) \\= u(t,X_t) + \int_t^T (\partial_s u + 1_D\mathcal{L} u)(s,X_s)ds + \int_t^T (\sigma 1_D\partial_x u)(s,X_s)dW_s + \int_t^T \partial_x u(s,X_s)dL_s\\ = u(t,X_t) + \int_t^T (\partial_s u + \mathcal{L} u)(s,X_s)ds + \int_t^T (\sigma 1_D\partial_x u)(s,X_s)dW_s + \int_t^T 1_{\partial D}(-\mathcal{L}u + \frac{a}{\kappa}\partial_x u)(s,X_s))ds. \end{gather}\] The deterministic integrals exist because of the regularity assumptions on \(u\) and its derivatives, and the stochastic integral is a martingale because of the regularity assumption on \(\sigma \partial u\).

The first integrand is identically zero, by ?? . The last integrand is identically zero, by the sticky boundary condition ?? . Taking \(\mathbb{E}[\cdot | X_t = x]\) shows that \(\mathbb{E}[u(X_T,T) | X_t = x] = u(x,t)\). Then we apply the terminal condition ?? to obtain the result. ◻

2.3 Stationary distribution↩︎

One way to see the role of parameter \(\kappa\) in determining the stickiness of the process 1 is to compute the stationary distribution. We will show (formally) that a locally stationary distribution is \[\label{eq:pi} \pi(x) = \tilde{\pi}(x)(1 + \kappa \delta(x)),\tag{6}\] where \(\tilde{\pi}\) is a locally stationary distribution for a Markov process with generator 4 assuming reflecting boundary conditions. That is, \(\tilde{\pi}\) satisfies \(\mathcal{L}^*\tilde{\pi} = 0\), where \(\mathcal{L}^*\) is the formal adjoint of \(\mathcal{L}\): \[\label{eq:Lstar} \mathcal{L}^* p = \partial_x(-b(x)p + \partial_x(a(x)p)).\tag{7}\] Furthermore we have \(-b\pi +\partial_x(a\pi)|_{x=0,\infty}=0\) by the reflecting boundary condition on \(\tilde{\pi}\). Hence, the \(\kappa\) is directly proportional to the atom of probability on the boundary.

We use a weak formulation of the usual condition \(\mathcal{L}^*\pi=0\) for a stationary distribution, which asks to show that \(\langle \pi,\mathcal{L}f\rangle=0\) for all \(f\) in the domain of the generator, where \(\langle f ,g\rangle=\int_0^\infty fgdx\) is the \(L_2\) inner product on \(\bar D\). We make use of an easy calculation: \[\label{eq:LLstar} p\mathcal{L}f = f\mathcal{L}^* p + \partial_x\Big(f\big(pb - \partial_x(pa)\big) + pa\partial_x f\Big).\tag{8}\] Now we may calculate: \[\begin{align} \langle \pi,\mathcal{L}f\rangle &= \tilde{\pi}\kappa \mathcal{L}f\big|_{x=0} + \Big(\tilde{\pi} a \partial_x f + f(\tilde{\pi} b - \partial_x(\tilde{\pi} a) \Big)\Big|_{x=0}^\infty + \int_0^\infty f\mathcal{L}^*\tilde{\pi} dx\\ &= \tilde{\pi}(\kappa \mathcal{L}f-a \partial_x f )\\ &= 0. \end{align}\] The last two terms in the first line were zero because of the condition that \(\tilde{\pi}\) was a locally stationary distribution for \(\mathcal{L}^*\). The final terms vanished by the boundary condition 5 .

2.4 Derivation via conservation of probability↩︎

We’ll now consider an alternative formulation of the evolution equations for a sticky diffusion, in terms of the evolution of its probability measure. This formulation can be derived from simple considerations of conservation of probability. Consider a process that has probability density \(p_1(x,t)\) for \(x>0\), and which has an atom of probability \(p_0(t)\) at 0, so the total density of the process can be written as \[\rho(x,t) = p_1(x,t)1(x>0) + p_0(t)\delta(x).\] Suppose \(p_1\) evolves with the diffusion dynamics described by the operator \(\mathcal{L}^*\) in 7 : \[\partial_t p_1 = \mathcal{L}^* p_1.\] Let’s write \(j = bp_1 - \partial_x(ap_1)\) for the flux associated with \(\mathcal{L}^*\), so that \(\mathcal{L}^*p_1 = -\partial_x j\).

Probability is exchanged with the atom via flux from the exerior: \[\partial_t p_0 = -j|_{x=0}.\] One can verify this condition conserves the total probability \(\int_0^\infty \rho dx\).

We need one additional condition to close this system and obtain a well-posed system of equations for \(p_0,p_1\). One such condition is \[p_0(t) = \kappa p_1(0,t)\] for some constant \(\kappa>0\). This condition was shown in [17], [32] to be consistent with the stickiness near zero arising from the limit of a strong, short-ranged force pushing the diffusion process towards 0.

Putting this together and using that \(\partial_t p_0=\kappa\partial_t p_1|_{x=0} = \kappa \mathcal{L}^* p_1|_{x=0}\), shows the probability density is \[\label{eq:rho} \rho(x,t) = p_1(x,t)(1+\kappa \delta(x))\tag{9}\] where \(p_1\) solves \[\label{eq:FP} \partial_t p_1 = \mathcal{L}^* p_1 \quad (x\geq 0), \qquad \kappa \mathcal{L}^* p_1 = -j \quad (x=0),\tag{10}\] with \(j = bp_1 - \partial_x(ap_1)\).

2.5 Equivalence of the forward and backward equations↩︎

Finally, let’s show that the operator \(\mathcal{L}^*\) in 10 with the given boundary condition, is the (formal) adjoint of \(\mathcal{L}\) with boundary condition 5 . We wish to show that \[\label{eq:adjoint} \langle f,\partial_t \rho\rangle = \langle \mathcal{L}f,\rho\rangle\tag{11}\] for all test functions \(f\) satisfying 5 , and for all \(\rho\) of the form 9 where \(p_1\) is a test function (infinitely differentiable with compact support). This is a weak form of the usual condition \(\langle f,\mathcal{L}^*\rho\rangle = \langle \mathcal{L}f,\rho\rangle\), and is more suited to handling \(\delta\)-functions. We calculate \[\langle f,\partial_t \rho\rangle = \int_0^\infty f\partial_t p_1 dx \;\;+ f\kappa\partial_t p_1\Big|_{x=0},\] and, using 8 , \[\begin{align} \langle \mathcal{L}f,\rho\rangle &= \int_0^\infty f\mathcal{L}^* p_1dx \;\; - \Big(fj + p_1a\partial_x f \Big)\Big|_{x=0} \; + \kappa p_1\mathcal{L}f\Big|_{x=0}\\ &= \int_0^\infty f\mathcal{L}^* p_1dx - fj\Big|_{x=0}\;\; , \end{align}\] by the sticky boundary condition 5 . Thus, we can see that 11 holds, exactly when \(p_1\) satisfies 10 .

3 Sticky Brownian motion↩︎

In this section we present our algorithm for simulating a sticky Brownian motion (SBM), which has drift \(b(x)=0\) and diffusion \(\sigma(x)=1\). We first give a high-level description of the algorithm and our numerical experiments (Section 3.1), then we show how to formally derive the algorithm, using simple Taylor expansions (Section 3.2). In Section 3.3 we present our main convergence results, with proofs in Section 3.4. Section 3.5 proves a sufficient condition for regularity of the solution to the associated backward equation.

3.1 Overview of the simulation algorithm↩︎

Consider a SBM \(\{Y_t\}_{t\in [0,T]}\), where \(T\) is the total length of time. We construct a numerical approximation \(\{X_k\}_{k=1,\ldots,N}\), where \(X_k\) is the value of the approximation at time \(t_k=kh\), and \(h\) is our time step, so that \(N=T/h\) is the total number of steps. We always assume that \(h\) is chosen such that \(N\) is an integer.

The numerical scheme works as follows. We divide the domain into two regions: \[\begin{align} \text{boundary layer:}\quad & \Omega_b=[0,\sqrt{3h}],\\ \text{exterior:}\quad & \Omega_e = (\sqrt{3h},\infty). \end{align}\] The algorithm proposes different updates, depending on whether \(X_k\) is in the boundary layer or the exterior.

Given a point in the exterior, \(X_k \in \Omega_e\), it is updated with a standard Euler-Maruyama (EM) increment. Specifically, \[\label{eqn:cont95scheme95outlayer} X_{k+1} = X_k + Z_{k},\tag{12}\] where \(\{Z_k\}_{k=0,\ldots,N-1} \sim \mathcal{U}([-\sqrt{3h}, \sqrt{3h}])\) is a sequence of iid uniform random variables, which have mean 0 and variance \(\text{Var}(Z_k) = h\). We choose uniform random variables rather than normal random variables to simplify our later analysis. Notice that, because of the choice of uniform steps, if \(X_k\in \Omega_e\), then \(X_{k+1} > 0\), so we do not need to consider reflection at this stage.

Figure 1: Plot of \lambda(x) for x\in[0,\sqrt{3h}] with various \kappa (shown in the legends) and time step h=10^{-4} (left) and h=4\times 10^{-4} (right). We always have \lambda(\sqrt{3h}) = 1, while \lambda(0) is a positive constant for all \kappa which decreases to 0 as \kappa increases. Increasing h has the effect of increasing \lambda(\cdot). == — Figure 1: Plot of \(\lambda(x)\) for \(x\in[0,\sqrt{3h}]\) with various \(\kappa\) (shown in the legends) and time step \(h=10^{-4}\) (left) and \(h=4\times 10^{-4}\) (right). We always have \(\lambda(\sqrt{3h}) = 1\), while \(\lambda(0)\) is a positive constant for all \(\kappa\) which decreases to 0 as \(\kappa\) increases. Increasing \(h\) has the effect of increasing \(\lambda(\cdot)\). ==

Set \(X_0=x\). Generate \(Z_k \sim \mathcal{U}([-\sqrt{3h},\sqrt{3h}])\). Update \(X_{k+1}\) using 12 . (\(X_k \in [0,\sqrt{3h}]\)) Calculate \(\lambda(X_k)\) using 14 . Generate \(u \sim \mathcal{U}([0,1])\). Update \(X_{k+1} = |X_k + Z_k|\). Set \(X_{k+1} = 0\).

Given a point \(X_k\) in the boundary layer, \(X_k\in \Omega_b\), it does one of two things: with probability \(\lambda(X_k)\) it is updated with a reflected EM increment [40], and with probability \(1-\lambda(X_k)\) it jumps exactly to \(0\): \[\label{eqn:cont95scheme95inlayer} X_{k+1} = \begin{cases} \left|X_k + Z_{k}\right|, & \text{with probability } \lambda(X_k) \\ 0, & \text{with probability } 1-\lambda(X_k). \\ \end{cases}\tag{13}\] We remark that this same update is performed when \(X_k=0\); if the second case in 13 is chosen, the process simply stays at 0.

We will show momentarily that we achieve good convergence properties when \(\lambda\) is the following function: \[\label{eqn:lambda95cont95scheme} \lambda(x) = \frac{x^2 + 2\kappa x + h}{\left(\frac{\kappa}{\sqrt{3h}}+1\right)x^2 + h + \kappa\sqrt{3h}}.\tag{14}\] Figure 1 shows \(\lambda(x)\) for different values of \(\kappa\) and \(h\). It has many sensible properties: one, \(\lambda\) is an increasing function of \(x\), with \(\lambda(\sqrt{3h}) =1\); this is expected because when \(X_k\) is further away from 0, it should have a smaller probability to jump there. Two, as \(\kappa\) increases, \(\lambda\) decreases, which also makes sense because when the boundary is stickier, there should be a greater chance of jumping there.

Our algorithm is summarized in algorithm [alg:cont95scheme]. Some realizations from our algorithm are illustrated in Figure 2. We can see that as sticky parameter \(\kappa\) increases, realizations spend more time within the boundary layer.

Figure 2: Simulations of an SBM with X_0=0 using Algorithm [alg:cont95scheme] with sticky parameter \kappa=0.6 (left) and \kappa=3 (right), and with time steps h \in \{10^{-2}, 10^{-3}, 10^{-5}\} as shown in the legend. The same realizations for the random variables are used for both plots for the same h. Readers can observe how stickiness causes trajectories to spend more time near the sticky boundary, before beginning excursions. — Figure 2: Simulations of an SBM with \(X_0=0\) using Algorithm [alg:cont95scheme] with sticky parameter \(\kappa=0.6\) (left) and \(\kappa=3\) (right), and with time steps \(h \in \{10^{-2}, 10^{-3}, 10^{-5}\}\) as shown in the legend. The same realizations for the random variables are used for both plots for the same \(h\). Readers can observe how stickiness causes trajectories to spend more time near the sticky boundary, before beginning excursions.

Note that it is essential that \(\lambda\) be a valid probability, i.e. that \(0\leq \lambda(x)\leq 1\). This is not guaranteed in the derivation to follow, but is something that must be checked a posteriori.

Lemma 1. The jump probability \(\lambda\) defined in 14 is a valid probability: \(\lambda(x) \in [0,1]\) for all \(x \in [0, \sqrt{3h}]\).

Proof. By factoring part of the denominator, we can write \(\lambda\) as \[\lambda(x) = \frac{x^2+h+2\kappa x}{x^2+h + 2\kappa x+ \frac{\kappa}{\sqrt{3h}}(x-\sqrt{3h})^2 }.\] From this is it clear that \(\lambda(x)> 0\), since every term in the expression is positive. It is also clear that \(\lambda(x)\leq 1\), with equality only when \(x=\sqrt{3h}\) or when \(\kappa=0\). ◻

3.2 How to choose \(\lambda(x)\)↩︎

The key step of this algorithm is in choosing the function \(\lambda\) which governs the probability of jumping directly to 0, instead of performing an EM update. This is done by solving for a \(\lambda\) that makes the transition operator of the discrete scheme a good approximation to the transition operator of SBM over one timestep. In this section we derive a formula for \(\lambda\) formally, using Taylor expansions of the generator for SBM. Our intent is to show the key steps one must follow to apply these ideas more generally; we prove rigorous convergence statements in subsequent sections.

Let the transition operator of our numerical scheme with timestep \(h\) be denoted by \(P^h\); it characterizes the evolution of statistics via \[(P^h f)(x) = \mathbb{E}_{x} f(X_1).\] We write \(\mathbb{E}_x(\cdot):=\mathbb{E}(\cdot|X_{0}=x)\). We wish to find \(\lambda\) such that \(P^h\) is a good approximation for the transition operator of the true dynamics; i.e. such that \[\label{eq:genapprox0} (P^h f)(x) -f(x) \approx \mathbb{E}_x \int_0^h \mathcal{L}f(Y_s)ds\,,\tag{15}\] where \(\mathcal{L}f = \frac{1}{2}\partial_{xx}f\) for a SBM. We have that \[\mathbb{E}_x\int_0^h \mathcal{L}f(Y_s)ds\approx h(\mathcal{L}f)(x),\] and \[\begin{align} (P^h f)(x) - f(x) &= \int_0^{+\infty} p_h(y|x)\left(f(t,y)-f(t,x)\right)\;dy\\ &\approx \partial_xf(t,x)\cdot \mathbb{E}_x(\Delta X) + \frac{1}{2}\partial_{xx}f(t,x)\cdot \mathbb{E}_x(\Delta X^2). \end{align}\] Here \(p_h(y|x)\) is the transition probability to move from \(x\) to \(y\) in one step of the scheme, and \(\Delta X := X_{1} - X_0\). Thus, asking that the approximations for both sides of 15 are within some tolerance say \(\epsilon\), our task is to find a \(\lambda\) such that \[\label{eq:genapprox} |\mathbb{E}_x(\Delta X)(\partial_xf) + \frac{1}{2} \mathbb{E}_x(\Delta X^2)(\partial_{xx}f) - \frac{h}{2}(\partial_{xx}f)| \leq \epsilon.\tag{16}\] Later, we will show that choosing \(\epsilon =O(h^2+hx)\) gives good convergence properties.

We must choose \(\lambda\) such that the coefficients of \(\partial_x f, \partial_{xx}f\) on the left-hand side of 16 are sufficiently small. This is not yet possible since there is only one term involving \(\partial_xf\). Thus, our next step is to approximate \(\partial_xf\) in terms of \(\partial_{xx}f\), using the sticky boundary condition \(\partial_xf(t,0) = \kappa \partial_{xx}f(t,0)\). Taylor-expanding \(\partial_xf, \partial_{xx}f\) around \(x\) shows that \[\begin{align} \partial_xf(t,0) &= \partial_xf(t,x) - x\partial_{xx}f(t,x) + \frac{x^2}{2}f^{(3)}(t,z_x),\\ \partial_{xx}f(t,0) &= \partial_{xx}f(t,x) - xf^{(3)}(t,y_x) , \end{align}\] for some \(z_x, y_x \in [0,x]\). From the sticky boundary condition, we have \[\label{eq:taylorsticky} \partial_xf(t,x) = \kappa\partial_{xx}f(t,x) +x\partial_{xx}f(t,x)- \kappa xf^{(3)}(t,y_x) - \frac{x^2}{2}f^{(3)}(t,z_x).\tag{17}\] Therefore, an approximation of \(\partial_xf(t,x)\) using only \(\partial_{xx}f(t,x)\) is \[\label{eq:stickyapprox} \partial_xf(t,x) = (\kappa+x)\partial_{xx}f(t,x) + O(h+x).\tag{18}\] We remark that to obtain an error of order \(O(h+x)\), we do not need the additional term \(x\partial_{xx}f(t,x)\); in fact, this term has the same magnitude as some higher-order terms that we are dropping. However, we found that without including this term, we obtained a function \(\lambda\) that was not always a valid probability, i.e. \(\lambda\notin[0,1]\).

Substituting for \(\partial_xf\) into 16 , shows that we must find a function \(\lambda\) such that \[\label{eq:genapprox2} |(\kappa+x)\mathbb{E}_x(\Delta X) + \frac{1}{2}\mathbb{E}_x(\Delta X^2) - \frac{h}{2} |\approx 0.\tag{19}\]

At this point we may evaluate \(\mathbb{E}_x(\Delta X), \mathbb{E}_x(\Delta X^2)\), and then solve 19 for \(\lambda\), treating 19 as an equality. We have that \[p_h(y|x) = (1-\lambda)1_{\{0\}}(y) + \frac{\lambda}{\sqrt{3h}}1_{[0,x+\sqrt{3h}]}(y) + \frac{\lambda}{\sqrt{3h}}1_{[0,|x-\sqrt{3h}|]}(y).\] From this, we can calculate \[\begin{align} \mathbb{E}_x(\Delta X_k) &=-x + \frac{\lambda\sqrt{3h}}{2} + \frac{\lambda x^2}{2\sqrt{3h}}\,, \nonumber\\ \mathbb{E}_x(\Delta X_k^2) &= -\frac{\lambda x^3}{\sqrt{3h}} + (1+\lambda)x^2 - \sqrt{3h}\lambda x + \lambda h\,. \label{eq:DelX} \end{align}\tag{20}\] Substituting into 19 and solving for \(\lambda\) gives gives Formula 14 .

Our calculations above have demonstrated two Lemmas, which we will use shortly in our convergence proofs.

Lemma 2. Consider a function \(u(t,x)\) satisfying the sticky boundary condition at \(x=0\), i.e. such that \(\frac{1}{2}\partial_xu(t,0) = \kappa \partial_{xx}u(t,0)\). Given \(x \in \Omega_b\) and \(t\in [0,+\infty)\), there exists \(z_x,y_x \in [0,x]\) such that \[\label{eqn:sticky95bc64x} \partial_xu(t,x) = (\kappa + x)\cdot \partial_{xx}u(t,x) - \kappa x u^{(3)}(t,y_x) - \frac{x^2}{2}u^{(3)}(t,z_x).\tag{21}\]

Proof. This follows from our calculations deriving 17 . ◻

Lemma 3. For \(\lambda\) chosen as 14 , and given \(x\in \Omega_b\), \[(\kappa+x)\mathbb{E}_x(\Delta X) + \frac{1}{2}\mathbb{E}_x(\Delta X^2) - \frac{h}{2} =0.\]

Proof. This follows from our calculations following 19 , which showed that we chose \(\lambda\) such that the above equation holds. ◻

3.3 Primary convergence results↩︎

Now we give an overview of our main convergence results. We wish to show weak convergence of the numerical scheme with order \(h\). That is, we wish to show that, for \(h\) sufficiently small, there is a constant \(C>0\) such that \[|\mathbb{E}\phi(X_N) - \mathbb{E}\phi(Y_T)| \leq Ch,\] where \(\phi\) is from a suitable class of functions. This is shown by decomposing the error into a local, one-step error, which is different for the boundary and the interior, and a contribution related to the number of steps accumulating each type of error. Our proof strategy closely follows that used to show convergence of diffusions satisfying Robin boundary conditions [40].

At this point we make a small but important remark about the assumed initial condition. Throughout the global-error analysis in this subsection, we fix a deterministic initial condition, and assume \[X_0=Y_0=x_0.\] All expectations and conditional expectations are taken under this initial condition. We make this assumption mainly in order to simplify notation and the readability of proofs, as it allows us to suppress the dependence on the initial condition in our notation. Thus, we suppress the dependency of \(\mathbb{E}\) on \(x_0\), except when the dependence on the starting point needs to be emphasized.

The only place that the initial condition could affect our estimates is in the global error analysis, where we could have a bounding constant depending on \(x_0\). For an SBM, the constants obtained below are uniform in \(x_0\). For the sticky diffusion scheme in Section 4, the present argument may produce constants that depend on \(x_0\). If one instead allows a random initial condition \(X_0\), the same argument may be applied conditionally on \(X_0\), after which one integrates with respect to the law of \(X_0\); in that case one must assume that the resulting \(x_0\)-dependent bound is integrable.

As is typical in proofs of weak convergence [40], [47], our main tool will be the backward equation ?? , repeated here for an SBM for ease of reference: \[\begin{gather} \label{eq:backSBM} \partial_t u + \frac{1}{2} \partial^2_{xx}u = 0 \qquad (x\geq 0), \\ \partial_x u(t,0) = \kappa \cdot \partial^2_{xx}u(t,0), \quad u(T,x) = \phi(x).\nonumber \end{gather}\tag{22}\] We may observe that for \(X_0=x_0\), \[\mathbb{E}(\phi(Y_T)) = u(0,x_0), \qquad \mathbb{E}(\phi(X_N)) = \mathbb{E}(u(T,X_N)),\] where the latter comes from observing that \(u(T,X_N) = \mathbb{E}[\phi(Y_T) | Y_T=X_N] = \phi(X_N)\). Therefore, writing \(t_k=kh\), \(u_k=u(t_k,X_k)\), we have \[\label{eq:telescope} \mathbb{E}\left(\phi(X_N)\right) - \mathbb{E}\left(\phi(Y_T)\right) = \mathbb{E}\left(\sum_{k=0}^{N-1}u_{k+1}-u_k\right).\tag{23}\] This sum gives us a way to decompose the global error into a sum of local errors, via: \[\begin{align} |\mathbb{E}\phi(X_N) - \mathbb{E}\phi(Y_T)| &= \left|\mathbb{E}\left(\sum_{k=0}^{N-1} \mathbb{E}\left(u_{k+1}-u_k|X_k\right)\right)\right|\nonumber\\ &\leq \mathbb{E}\left(\sum_{k=0}^{N-1} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k\right)\right|\cdot \mathbf{1}_{X_k \in [0,\sqrt{3h}]}\right) \nonumber\\&\quad + \mathbb{E}\left(\sum_{k=0}^{N-1} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k\right)\right|\cdot \mathbf{1}_{X_k \notin [0,\sqrt{3h}]}\right). \label{eq:err95decompose} \end{align}\tag{24}\] The global error thus arises from two things: the local errors \(\left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k\right)\right|\) in the boundary and exterior, and the number of steps in each of these states.

To present our main results we need some assumptions.

Assumption 1. The solution \(u(t,x)\) of 22 has uniformly bounded spatial derivatives up to order \(4\). That is, \[\|\partial_x^j u\|_\infty := \sup_{(t,x)\in [0,T]\times[0,\infty)} \left|\partial_x^j u(t,x)\right| <\infty, \qquad j=1,2,3,4.\]

Assumption 2. The solution \(u(t,x)\) of 22 has time derivative up to order 2 where \(\left\|\partial_{tt} u\right\|_{\infty} < \infty\).

We further suppose that \(h\) is always chosen such that \(N := T/h \in \mathbb{Z}\).

Theorem 2. (Local Error of \(\{X_k\}_{k\in\mathbb{N}}\) in Algorithm [alg:cont95scheme])
Under Assumptions 1, 2, there exists a constant \(M\) depending on \(\left\| u^{(3)} \right\|_{\infty}\) \(\left\| u^{(4)} \right\|_{\infty}\) and \(\left\|\partial_{tt} u\right\|_{\infty}\), such that \[\label{eqn:local95err95cont95scheme} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| \leq \begin{cases} M( h^2 + hx), & \forall x\in [0,\sqrt{3h}]\\ Mh^2 & \forall x > \sqrt{3h}. \end{cases}\tag{25}\]

Thus, the local error of the scheme in the boundary layer is \(O(h^{3/2})\) and in the exterior is \(O(h^2)\).

A naive estimate for the global error would be \(O(h^{1/2})\), which comes from taking the worst-case local error of \(O(h^{3/2})\), and multiplying by the total number of steps which is \(O(h^{-1})\). This however is too pessimistic, as will turn out that the number of steps in the boundary layer is \(O(h^{-1/2})\), leading to a global error of \(O(h^{3/2}\cdot h^{-1/2}) + O(h^2\cdot h^{-1}) \sim O(h)\).

The specific type of estimate we must construct on the number of steps in the boundary layer arises from this corollary.

Corollary 1. Under the assumptions of Theorem 2, \[\mathbb{E}\left(\sum_{k=0}^{N-1} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k\right)\right|\cdot \mathbf{1}_{X_k \in [0,\sqrt{3h}]}\right) \leq M_1\left(h + h \cdot \mathbb{E}\left(\sum_{k=0}^{N-1} X_k \cdot \mathbf{1}_{X_k \in [0,\sqrt{3h}]} \right)\right).\]

Proof. This is simply the expectation applied to the first case of 25 , and summed over all steps. ◻

Our main lemma to estimate the number of points in the boundary layer is the following.

Lemma 4. For sufficiently small \(h\), there exists a constant \(C>0\), independent of \(h\) and initial position \(x_0\ge0\), such that \[\mathbb{E}\left(\sum_{k=0}^{N-1} X_k \cdot \mathbf{1}_{X_k \in [0,\sqrt{3h}]} \right) < C.\]

This lemma implies the number of points in the boundary layer is \(O(h^{-1/2})\), because \(|X_k| =O(h^{1/2})\) in the boundary layer. It is the most technical lemma to prove. Applying this Lemma to 24 gives the global error bound, our main convergence theorem.

Theorem 3. (Global Error of \(\{X_k\}_{k\in\mathbb{N}}\) in Algorithm [alg:cont95scheme])
Let \(\phi:[0,\infty)\to \mathbb{R}\) be such that the associated backward solution of 22 , \[u(t,x):=\mathbb{E}\!\left[\phi(Y_T)\mid Y_t=x\right],\] satisfies Assumptions 1 and 2. Then there exist constants \(\delta>0\) and \(C>0\), independent of \(h\) and initial position \(x_0\ge0\), such that for all \(h<\delta\) and \(x_0 \geq 0\), \[\label{eqn:global95err95cont95scheme} \left|\mathbb{E}\bigl(\phi(X_N)-\phi(Y_T)\bigr)\right|\le Ch.\tag{26}\]

Remark 1. The main convergence theorems are stated in terms of regularity assumptions on the backward solution \(u\), rather than directly on the terminal datum \(\phi\). A sufficient condition on \(\phi\) ensuring Assumptions 1 and 2 is given later in Section 3.5; see in particular Proposition 4. For simplicity, that sufficient condition is stated on a compact interval \([0,L]\) with a reflecting boundary condition at \(x=L\) (general Wentzell conditions). In fact, at \(0\) and \(L\), we can generalize the proof of Proposition 4 to any Wentzell condition.

Proof. This follows from 24 , applying Theorem 2, Corollary [corollary:local95err], and Lemma 4. ◻

Thus, our scheme is weakly convergent with order 1.

Numerical experiments. The convergence results are illustrated numerically in Figure 3. We run our algorithm with \(X_0=0, T=1\) and estimate \(\mathbb{E}[X_T]\). In our algorithm we additionally impose a reflecting boundary condition at \(x=L\), which is implemented in our algorithm using a reflected EM update (similar to the update in our algorithm at \(x=0\)). This is known to converge weakly with order 1 [40]. We compare our estimated mean, to the solution \(u(0,0)\) to the backward equation 22 with \(\phi(x) = x\). This solution is computed using a finite-difference scheme with with extremely small space and time steps; specifically, \(\Delta x = 10^{-3}\) and \(\Delta t = \frac{\Delta x^2}{8*5} = 2.5\times10^{-8}\).

Figure 3 demonstrates \(O(h)\) convergence for two different values of the sticky parameter. It is notable that our choice of \(\phi\) does not satisfy the compatibility conditions given in Proposition 4, which are used to show that Assumptions 1 and 2 hold for the solution to the backward equation. Nevertheless, this does not appear to affect our empirical order of convergence.

Figure 3: Log-log plot for the weak error e_h = |\mathbb{E}(X_T|X_0=0)-u(0,0)| of our simulation scheme [alg:cont95scheme] using \kappa = 1 (left) and \kappa=10 (right). Both comparisons used L=2,T=1. Error bars are \pm 1 empirical standard deviation. The left plot takes time step-size h \in \{1/64, 1/128, 1/256, 1/512, 1/1024\}. To keep the relative error bar similar, we quadruple the number of realizations every time we halve the time step size starting from 2\times 10^8 realization. The least square fit gives the slope of 0.956615. The right plot takes time step-size h \in \{1/32, 1/64, 1/128, 1/256, 1/512\}. Similarly, to keep the relative error bar similar, we start with 8 \times 10^8 realizations and we quadruple the number of realizations every time we halve the time step size. The least square fit for the last three data points has the slope of 1.011378. — Figure 3: Log-log plot for the weak error \(e_h = |\mathbb{E}(X_T|X_0=0)-u(0,0)|\) of our simulation scheme [alg:cont95scheme] using \(\kappa = 1\) (left) and \(\kappa=10\) (right). Both comparisons used \(L=2,T=1\). Error bars are \(\pm 1\) empirical standard deviation. The left plot takes time step-size \(h \in \{1/64, 1/128, 1/256, 1/512, 1/1024\}\). To keep the relative error bar similar, we quadruple the number of realizations every time we halve the time step size starting from \(2\times 10^8\) realization. The least square fit gives the slope of 0.956615. The right plot takes time step-size \(h \in \{1/32, 1/64, 1/128, 1/256, 1/512\}\). Similarly, to keep the relative error bar similar, we start with \(8 \times 10^8\) realizations and we quadruple the number of realizations every time we halve the time step size. The least square fit for the last three data points has the slope of 1.011378.

3.4 Convergence proofs↩︎

It remains to prove Theorem 2 and Lemma 4. We start with a couple of lemmas regarding the magnitudes of the moments of \(\Delta X_k\). In this section, we use notation \(\mathbb{E}_xf(\Delta X_k) \equiv \mathbb{E}[f(\Delta X_k)|X_k = x]\) for a function \(f\) of the increments.

Lemma 5. Given the sequence \(\{X_k\}\) generated from algorithm [alg:cont95scheme] with \(\kappa>0\), there is a constant \(C\) such that \[\sup_{x\in [0,\infty)} |\mathbb{E}_x(\Delta X_k^2)| \leq Ch\,, \qquad \sup_{x\in [0, \sqrt{3h}]} |\mathbb{E}_x(\Delta X_k)| \leq Ch.\]

The second part of the lemma is not trivial, as we might expect \(|\mathbb{E}_x(\Delta X_k)|\sim O(h^{1/2})\), since this is the size of a typical jump to 0.

Proof. The bound on \(\mathbb{E}_x(\Delta X_k^2)\) follows trivially from the fact that \(|\Delta X_k^2|\leq 3h\). The bound on \(\mathbb{E}_x(\Delta X_k)\) then follows from Lemma 3, which shows that \(\mathbb{E}_x(\Delta X_k)\) is a sum of terms which are all \(O(h)\) or higher provided \(\kappa >0\). ◻

Lemma 6. There exists a constant \(C\) such that \[\sup_{x\in [0, \sqrt{3h}]}\mathbb{E}_x(|\Delta X_k|^3)\leq C(h^2 + h x)\,.\]

Proof. Fix \(x\in[0,\sqrt{3h}]\), an increment \(\Delta X_k\) is supported on \([-x,\sqrt{3h}]\). We decompose \(\Delta X_k\) into its positive and negative parts, \[\Delta X_k^+=\max\{\Delta X_k,0\},\qquad \Delta X_k^-=\max\{-\Delta X_k,0\},\] so that \(|\Delta X_k|^3=(\Delta X_k^+)^3+(\Delta X_k^-)^3\) and therefore \[\mathbb{E}_x(|\Delta X_k|^3) =\mathbb{E}_x[(\Delta X_k^+)^3]+\mathbb{E}_x[(\Delta X_k^-)^3].\]

Positive part. Whenever \(\Delta X_k>0\) the increment comes from the continuous branch \(Z_k\) in Algorithm [alg:cont95scheme]. This means the largest positive increment is \(\sqrt{3h}\), hence \[\label{eq:3nd95mom95pospart} \mathbb{E}_x[(\Delta X_k^+)^3]\le \lambda(x)\,(3h)^{3/2}.\tag{27}\] To bound \(\lambda(x)\) in 14 , recall that for \(x\in[0,\sqrt{3h}]\) we have \(x^2\le 3h\). The denominator is \[\Bigl(\frac{\kappa}{\sqrt{3h}}+1\Bigr)x^2+h+\kappa\sqrt{3h} \;\ge\;\kappa\sqrt{3h},\] The numerator obeys \(x^2+2\kappa x+h\le 4h+2\kappa x\). Therefore there exist constants \(A,B>0\), independent of \(h\) and \(x\), such that \[\label{eq:3nd95mom95lambdabound} \lambda(x)\leq \frac{4\sqrt{h}}{\sqrt{3}\,\kappa} + \frac{2x}{\sqrt{3h}}\le A \sqrt{h} + B\frac{x}{\sqrt{h}}.\tag{28}\] Combining 28 with 27 shows there exist constants \(C_1,C_2>0\) for which \[\label{eq:3nd95mom95pos-final} \mathbb{E}_x[(\Delta X_k^+)^3]\le C_1 h^2 + C_2 h x.\tag{29}\]

Negative part. Because \(\Delta X_k\ge -x\), we have \(0\le \Delta X_k^-\le x\), and therefore \[(\Delta X_k^-)^3 \le x^3.\] Since \(x\in [0,\sqrt{3h}]\), \[\label{eq:3nd95mom95neg-final} \mathbb{E}_x[(\Delta X_k^-)^3] \le x^3 \le 3hx.\tag{30}\]

Combining 29 and 30 gives \[\mathbb{E}_x(|\Delta X_k|^3) \;\le\; C_1 h^2 + (C_2+3)hx.\] The bound holds for \(|\mathbb{E}_x(\Delta X_k^3)|\) with \(C=\max\{C_1,C_2+3\}\). This proves the lemma. ◻

We are now ready to prove the local error.

Proof of Theorem 2. Let \(u\) be the solution to the backward equation 22 . We perform a Taylor expansion of \(u\) to 3rd order centered at \((t_k, X_k)\). We write \(u_k=u(t_k,X_k), \partial_xu_k = \partial_xu(t_k,X_k)\) (and similarly for \(\partial_{xx}u_k, \partial_{xxx}u_k\)), and we write \(Y_k\) for a point such that \(Y_k\in[X_k,X_{k+1}]\) or \(Y_k\in[X_{k+1},X_{k}]\), and \(s_k\) will be a point in interval \(s_k\in[kh, (k+1)h]\)). Expanding gives \[\label{eq:taylor4} \begin{align} u_{k+1} - u_k &= \partial_t u_k \cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot \Delta X_k^2\\ &\quad + \frac{1}{2}\partial_{tt} u(s_k, X_k) \cdot h^2 + \frac{1}{6}\partial_{xxx}u(t_k, Y_{k})\cdot \Delta X_k^3. \end{align}\tag{31}\]

We are interested in \(\mathbb{E}_x(u_{k+1}-u_k)\). We start by assuming \(x\in \Omega_b\). We now bound the expectation of each of the terms in the Taylor-expansion 31 . We have \[\begin{align} \mathbb{E}_x|\partial_{tt} u(s_k, X_k) \cdot h^2| & \leq h^2\|\partial_{tt}u\|_{\infty}\\ \mathbb{E}_x|\partial_{xxx}u(t_k, Y_{k})\cdot \Delta X_k^3| &\leq C(h^2+hx)\|u^{(3)}\|_\infty, \end{align}\] where the second inequality follows from Lemma 6.

We also have, using the fact that \(u\) solves the backward equation 22 and using Lemma 2, \[\begin{align} & \mathbb{E}_x\left[\partial_t u_k\cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot (\Delta X_k^2) \right] \\ &\quad = -\frac{1}{2}\partial_{xx} u(t_k,x)\cdot h + (\kappa + x)\cdot\partial_{xx}u(t_k,x) \cdot \mathbb{E}_x(\Delta X_k) + \frac{1}{2}\partial_{xx}u(t_k,x)\cdot \mathbb{E}_x(\Delta X_k^2)\\ &\qquad \qquad - \mathbb{E}_x(\Delta X_k)\left(\kappa x\cdot u^{(3)}(t_k,y_x) + \frac{x^2}{2}\cdot u^{(3)}(t_k,z_x)\right) \\ &\quad = - \mathbb{E}_x(\Delta X_k)\left(\kappa x\cdot u^{(3)}(t_k,y_x) + \frac{x^2}{2}\cdot u^{(3)}(t_k,z_x)\right). \end{align}\] Here \(y_x,z_x\in [0,x]\), and the last step used Lemma 3.

Therefore, applying the triangle inequality to 31 , and then using Lemma 5, shows that there exist constants \(M_1,M_2\), independent of \(x,t,h\), such that \[\begin{align} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| & \leq \left|\mathbb{E}_x(\Delta X_k)\right|\cdot \left(\kappa x \|u^{(3)}\|_{\infty} + \frac{h}{2}\|u^{(3)}\|_{\infty}\right) \\ &\qquad\qquad+ h^2\|\partial_{tt}u\|_{\infty} + C(h^2+hx)\|u^{(3)}\|_\infty\\ &\leq M_1 h^2 + M_2 hx. \end{align}\]

We now check \(x \in \Omega_e\). Recall that, outside the boundary layer, the simulation scheme makes a standard Euler Maruyama jump, with moments \[\mathbb{E}_x(\Delta X_k) = \mathbb{E}_x(\Delta X_k^3) = 0\,, \quad \mathbb{E}_x(\Delta X_k^2) = h, \quad \mathbb{E}_x(\Delta X_k^4) = \frac{9}{5}h^2.\] We exploit the Taylor expansion up to fourth order, \[\label{eq:taylor495ext} \begin{align} u_{k+1} - u_k &= \partial_t u_k \cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot \Delta X_k^2\\ &\quad + \frac{1}{2}\partial_{tt} u(s_k, X_k) \cdot h^2 + \frac{1}{6}\partial_{xxx}u_k\cdot \Delta X_k^3 + \frac{1}{24}\partial_{xxxx}u(t_k, Y_{k})\cdot \Delta X_k^4. \end{align}\tag{32}\] Using the fact that \(u\) solves the backward equation 22 , we know, \[\begin{align} & \mathbb{E}_x\left[\partial_t u_k\cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot (\Delta X_k^2) \right] \\ &\quad = -\frac{1}{2}\partial_{xx} u(t_k,x)\cdot h + \partial_xu_k \cdot 0 + \frac{1}{2}\partial_{xx}u(t_k,x)\cdot h\\ &\quad = 0. \end{align}\] Apply the Taylor expansion 31 and the fact that \(\mathbb{E}_x(\Delta X_k^3) = 0\), \[\begin{align} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| &= \left|\frac{1}{2}\partial_{tt} u(s_k, x)\cdot h^2 + 0 + \frac{3}{40}\partial_{xxxx}u(t_k, Y_{k})\cdot h^2\right|\\ &\leq h^2 \left(\|\partial_{tt}u\|_{\infty} + \|\partial_{xxxx}u\|_{\infty}\right). \end{align}\] ◻

It remains to prove Lemma 4. For this, we first record a standard discrete backward equation for additive functionals of a Markov chain, proven here for completeness [40], [47].

Lemma 7. Let \(\{X_i\}_{i=0}^N\) be the Markov chain generated by our numerical scheme, with time grid \(t_i=ih\), \(i=0,\dots,N\), and \(t_N=T\). Let \(P\) denote the transition operator of the scheme, i.e. whenever the conditional expectation is well defined, \[PU(t_i,x):=\mathbb{E}\bigl[U(t_{i+1},X_{i+1})\,\big|\,X_i=x\bigr], \qquad i=0,\dots,N-1.\] Let \[g:\{t_0,\dots,t_{N-1}\}\times [0,\infty)\to \mathbb{R}\] be measurable, and assume that for every \(i=0,\dots,N\) and every \(x\in [0,\infty)\), \[\mathbb{E}\left(\sum_{k=i}^{N-1} |g(t_k,X_k)|\,\middle|\, X_i=x\right)<\infty.\] Define \[U(t_i,x):=\mathbb{E}\left(\sum_{k=i}^{N-1} g(t_k,X_k)\,\middle|\, X_i=x\right), \qquad i=0,\dots,N,\] with the convention that the sum is empty when \(i=N\). Then \(U\) satisfies \[\label{eq:V} \begin{align} PU(t_i,x)-U(t_i,x) &= -g(t_i,x), \qquad (t_i,x)\in \{t_0,\dots,t_{N-1}\}\times [0,\infty),\\ U(t_N,x) &= 0, \qquad x\in [0,\infty). \end{align}\tag{33}\]

Proof. By assumption, the random variable \(\sum_{k=i}^{N-1} g(t_k,X_k)\) is integrable under \(\mathbb{P}(\cdot\mid X_i=x)\), so \(U(t_i,x)\) is well defined. When \(i=N\), the sum is empty, hence \[U(t_N,x)=0, \qquad x\in [0,\infty).\] Now fix \(i\in\{0,\dots,N-1\}\). Since \(X_i=x\) on the conditioning event, \[\begin{align} U(t_i,x) &= g(t_i,x)+\mathbb{E}\left(\sum_{k=i+1}^{N-1} g(t_k,X_k)\,\middle|\, X_i=x\right). \end{align}\] Applying the tower property, \[\begin{align} \mathbb{E}\left(\sum_{k=i+1}^{N-1} g(t_k,X_k)\,\middle|\, X_i=x\right) &= \mathbb{E}\left( \mathbb{E}\left(\sum_{k=i+1}^{N-1} g(t_k,X_k)\,\middle|\, X_{i+1}\right) \,\middle|\, X_i=x\right)\\ &= \mathbb{E}\bigl[U(t_{i+1},X_{i+1})\,\big|\,X_i=x\bigr]\\ &= PU(t_i,x). \end{align}\] Therefore, \[U(t_i,x)=g(t_i,x)+PU(t_i,x).\] ◻

Now, to prove Lemma 4, it suffices to show that

\[U(t_0,x)=\mathbb{E}\left(\sum_{k=0}^{N-1} X_k \mathbf{1}_{\{X_k\in[0,\sqrt{3h}]\}} \,\middle|\, X_0=x\right)\] is uniformly bounded. From Lemma 7, we know that \(U\) solves 33 with \(g(t,x) = x \cdot \mathbf{1}_{x \in [0,\sqrt{3h}]}\). Therefore, it is enough to construct a uniformly bounded function \(V\) such that \[\begin{align} PV(t_i,x)-V(t_i,x) &\le -x\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}, \qquad i=0,\dots,N-1,\\ V(T,x) &\ge 0. \end{align}\] Indeed, these two inequalities imply \(U(t_i,x)\le V(t_i,x)\) for all grid times by backward induction (shown explicitly in the proof), and in particular \[U(t_0,x)\le V(t_0,x).\] Hence the uniform boundedness of \(V\) yields the desired bound on \(U\).

Proof of Lemma 4. We seek a bounded test function whose one-step drift is negative by an amount comparable to \(x\) whenever the chain lies in the boundary layer. Once such a function is found, a backward comparison argument will imply \(U\le V\), and hence the desired uniform bound. The function \(e^{w(x)}\) constructed below is maximal near the sticky boundary and decreases away from it, so each visit to \([0,\sqrt{3h}]\) forces a definite loss in one step. The exponential factor \(e^{K(T-t_i)}\) is then used later to absorb the order-\(h\) remainder terms in this drift estimate. Readers will see in 36 later that the function \(V(t,x)\) we construct has a structure similar to (3.27) in [40].

We introduce the following function, where \(M,s>16\): \[\label{eqn:chosen95wx95continuous} w(x) = \begin{cases} M , & x\in [0, \sqrt{3h}/2] \\ M - s\left(x-\frac{\sqrt{3h}}{2}\right), & x\in (\sqrt{3h}/2,+\infty) . \end{cases}\tag{34}\]

A direct case-by-case computation, provided in the supplement 6.1, shows that there exist constants \(C_1(s,M)>0\) and \(\delta>0\), independent of \(x\) and \(h\), such that for all \(h<\delta\), \[\label{eqn:wx95goal} Pe^{w(x)} \leq e^{w(x)}\left(1-f_h(x)+C_1h\right), \qquad \forall x\ge 0,\;\forall h<\delta,\tag{35}\] where \(f_h\) is defined piecewise by 66 , 69 , and 70 in the supplement, and satisfies \[f_h(x)\geq x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}.\] Now we introduce a new function \[\label{eqn:V95define} V(t_i,x) := e^{K(T-t_i)}e^{w(x)}\,,\quad (t,x)\in \{t_0,\dots,t_N\}\times[0,\infty).\tag{36}\]

Throughout the remainder of the proof, constants are allowed to depend on the fixed parameters \(s,M\) (and on \(K\) after \(K\) is chosen), but are independent of \(x\) and \(h\). After possibly decreasing \(\delta\), we may also assume that \[e^{-Kh}=1-Kh+r_K(h), \qquad |r_K(h)|\le C_2(K)h^2, \qquad \forall h<\delta.\] Hence, for every \(h<\delta\), \[\label{eq:PV95step1} \begin{align} PV(t_i,x)-V(t_i,x) &= e^{K(T-t_i-h)}Pe^{w(x)} - e^{K(T-t_i)}e^{w(x)}\\ &= e^{K(T-t_i)}\Bigl(e^{-Kh}Pe^{w(x)}-e^{w(x)}\Bigr)\\ &\le e^{K(T-t_i)}e^{w(x)} \Bigl[\bigl(1-Kh+C_2(K)h^2\bigr)\bigl(1-f_h(x)+C_1h\bigr)-1\Bigr]. \end{align}\tag{37}\] Expanding the bracket gives \[\label{eq:PV95step2} \begin{align} &\bigl(1-Kh+C_2(K)h^2\bigr)\bigl(1-f_h(x)+C_1h\bigr)-1\\ &\qquad = -f_h(x) + (C_1-K)h + K h f_h(x) + C_2(K)h^2\bigl(1-f_h(x)+C_1h\bigr)-K C_1 h^2. \end{align}\tag{38}\]

Next, note that \(f_h(x)=0\) for \(x>\sqrt{3h}\). On the other hand, when \(x\in[0,\sqrt{3h}]\), the explicit formulas 66 –70 in the Supplement together with \(0\le \lambda(x)\le 1\) imply that \[0\le f_h(x)\le C_3\sqrt{h}, \qquad \forall x\ge 0,\;\forall h<\delta,\] for some constant \(C_3>0\) depending only on \(s\). Therefore \[h f_h(x)\le C_3 h^{3/2}, \qquad \forall x\ge 0,\;\forall h<\delta.\] After decreasing \(\delta\) once more so that \(\delta\le 1\), the \(h^2\)-terms in 38 can be bounded by \(C_4(K)h^{3/2}\) for some constant \(C_4(K)>0\). Consequently, there exists a constant \(C_5(K)>0\) such that \(\forall x\ge 0,\;\forall h<\delta\), \[\label{eq:PV95step3} \bigl(1-Kh+C_2(K)h^2\bigr)\bigl(1-f_h(x)+C_1h\bigr)-1 \le -f_h(x)+(C_1-K)h + C_5(K)h^{3/2}.\tag{39}\] Substituting 39 into 37 , we obtain, \(\forall x\ge 0,\;\forall h<\delta\) \[\label{eq:PV95final} PV(t_i,x)-V(t_i,x) \le -\,e^{K(T-t_i)}e^{w(x)} \Bigl( f_h(x) + (K-C_1)h - C_5(K)h^{3/2} \Bigr).\tag{40}\] Now choose \(K:=C_1+1\). After decreasing \(\delta\) again if necessary, we may assume that \[(K-C_1)h-C_5(K)h^{3/2}\ge 0, \qquad \forall h<\delta.\] Hence 40 yields \[PV(t_i,x)-V(t_i,x)\le -\,e^{K(T-t_i)}e^{w(x)}f_h(x), \qquad \forall x\ge 0,\;\forall h<\delta.\] Finally, after decreasing \(\delta\) one last time so that \[M-\frac{s}{2}\sqrt{3h}\ge 0, \qquad \forall h<\delta,\] we have \(w(x)\ge 0\) for every \(x\in[0,\sqrt{3h}]\), and therefore \[e^{K(T-t_i)}e^{w(x)}\ge 1 \qquad \text{on } [0,\sqrt{3h}].\] Since also \(f_h(x)\ge x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}\), we conclude that \[\label{eq:PV95final95final} PV(t_i,x)-V(t_i,x)\le -\,x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}, \qquad \forall x\ge 0,\;\forall h<\delta.\tag{41}\] It remains to compare \(V\) with the function \[U(t_i,x):=\mathbb{E}\left(\sum_{k=i}^{N-1} X_k\mathbf{1}_{\{X_k\in[0,\sqrt{3h}]\}}\,\middle|\,X_i=x\right), \qquad i=0,\dots,N.\] By Lemma 7, \(U\) satisfies \[\label{eq:PU} PU(t_i,x)-U(t_i,x)=-x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}, \qquad i=0,\dots,N-1,\tag{42}\] with terminal condition \[U(T,x)=0.\] On the other hand, we have shown that for all sufficiently small \(h\), \[PV(t_i,x)-V(t_i,x)\le -\,x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}, \qquad i=0,\dots,N-1.\] Moreover, \[V(T,x)=e^{w(x)}\ge 0 = U(T,x).\]

We now prove by backward induction that \[U(t_i,x)\le V(t_i,x), \qquad \forall i=0,\dots,N,\;\forall x\ge 0.\] The claim is true for \(i=N\) by the terminal inequality above. Assume it holds at time \(t_{i+1}\), namely \(U(t_{i+1},x)\le V(t_{i+1},x)\), for all \(x\ge 0\). This implies \[PU(t_i,x) = \mathbb{E}\left(U(t_{i+1}, X_{i+1})\mid X_i=x\right)\le \mathbb{E}\left(V(t_{i+1}, X_{i+1})\mid X_i=x\right) = PV(t_i,x), \qquad \forall x\ge 0.\] Using 42 , the supersolution inequality 41 , and \(f_h(x) \geq x \cdot \mathbf{1}_{x \in [0,\sqrt{3h}]}\) we obtain \[\begin{align} U(t_i,x) &= PU(t_i,x)+x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}\\ &\le PV(t_i,x)+x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}\\ &\le V(t_i,x). \end{align}\] This closes the induction. Finally, since \(w(x)\le M\) for all \(x\ge 0\), we have \[V(t_i,x)=e^{K(T-t_i)}e^{w(x)}\le e^{KT+M}, \qquad \forall i=0,\dots,N,\;\forall x\ge 0.\] Hence \[U(t_0,x)\le V(t_0,x)\le e^{KT+M}, \qquad \forall x\ge 0.\] Therefore there exists a constant \(C = e^{KT+M}>0\), independent of \(h, x_0\), such that \[\mathbb{E}\left(\sum_{k=0}^{N-1} X_k\mathbf{1}_{\{X_k\in[0,\sqrt{3h}]\}}\right)\le C.\] ◻

As a side remark, we in fact use a special case of 3.19-3.21 in Leimkuhler’s paper [40] where \(q(\cdot) = 1\).

3.5 Sufficient conditions for backward equation regularity↩︎

The main convergence theorems are stated in terms of regularity assumptions on the backward solution \(u\), rather than directly on the terminal datum \(\phi\). It is therefore natural to ask for a concrete condition on \(\phi\) which guarantees that these regularity assumptions hold. On the half-line this requires a separate regularity theory for the sticky/Wentzell problem, which is not esaily available in the literature, so here we record a convenient sufficient condition on a compact interval \([0,L]\). For simplicity, we impose a reflecting boundary condition at \(x=L\); the same argument extends to a general Wentzell boundary condition, as noted in the remark after proving Proposition 4.

The key point is that the compatibility conditions are not ad hoc. The boundary conditions imposed on \(\phi\) ensure that \(\phi\) belongs to the domain of the spatial generator \(A\) of the SBM, while the additional conditions obtained by applying the same boundary conditions to \(\frac{1}{2}\phi''\) ensure that \(A\phi\) also belongs to the domain of \(A\), that is, \(\phi\in D(A^2)\). Once this is established, standard semigroup regularity for the associated boundary-value problem yields uniform bounds on the second time derivative and on the spatial derivatives up to order four. The next proposition makes this precise.

Proposition 4. Let \(\phi\in C^4([0,L])\), and let \(u(t,x)\) solve \[\label{eq:back95compact95reflecting} \begin{align} \partial_t u + \frac{1}{2} \partial_{xx}u &= 0, && (t,x) \in [0,T)\times[0,L], \label{eq:back95compact95reflecting95pde}\\ \partial_x u(t,0) &= \kappa\,\partial_{xx}u(t,0), && t\in[0,T), \label{eq:back95compact95reflecting95sticky}\\ \partial_x u(t,L) &= 0, && t\in[0,T), \label{eq:back95compact95reflecting95reflecting}\\ u(T,x) &= \phi(x), && x\in[0,L]. \label{eq:back95compact95reflecting95terminal} \end{align}\] {#eq: sublabel=eq:eq:back95compact95reflecting,eq:eq:back95compact95reflecting95pde,eq:eq:back95compact95reflecting95sticky,eq:eq:back95compact95reflecting95reflecting,eq:eq:back95compact95reflecting95terminal} Assume that \(\phi\) satisfies compatibility conditions \[\phi'(0)=\kappa \phi''(0), \qquad \phi'(L)=0,\] and \[\phi^{(3)}(0)=\kappa \phi^{(4)}(0), \qquad \phi^{(3)}(L)=0.\] Then \[\sup_{(t,x)\in[0,T]\times[0,L]} |\partial_x^j u(t,x)|<\infty, \qquad j=1,2,3,4,\] and \[\sup_{(t,x)\in[0,T]\times[0,L]} |\partial_{tt}u(t,x)|<\infty.\] In particular, Assumptions 1 and 2 hold on the compact strip \([0,T]\times[0,L]\).

Proof. Set \(v(s,x):=u(T-s,x)\) for \((s,x)\in[0,T]\times[0,L]\). Then \(v\) solves \[\begin{align} \partial_s v &= \frac{1}{2} \partial_{xx}v, && (s,x) \in (0,T]\times[0,L],\\ \partial_x v(s,0) &= \kappa\,\partial_{xx}v(s,0), && s\in (0,T],\\ \partial_x v(s,L) &= 0, && s\in (0,T],\\ v(0,x) &= \phi(x), && x\in[0,L]. \end{align}\]

Let \(X:=C([0,L])\), and define \(Af:=\frac{1}{2} f''\) on \[D(A):=\left\{f\in C^2([0,L]) : f'(0)=\kappa f''(0),\;f'(L)=0\right\}.\] By Theorem 1.1 in [48] for generalized Wentzell boundary conditions on \(C([0,L])\), \(A\) generates an analytic semigroup \(\{T(s)\}_{s\ge0}\) on \(X\). Hence \(v(s)=T(s)\phi\).

The first pair of conditions on \(\phi\) says exactly that \(\phi\in D(A)\). Since \(A\phi=\frac{1}{2} \phi''\), the second pair of conditions says exactly that \(A\phi\in D(A)\). Therefore \(\phi\in D(A^2)\). By the standard differentiability theorem for \(C_0\)-semigroups ([49] Section 2.4), since \(\phi\in D(A^2)\), we have \(v(s)\in D(A^2)\) for every \(s\in[0,T]\), and \[v_s(s)=T(s)A\phi, \qquad v_{ss}(s)=T(s)A^2\phi.\] Since \(\{T(s)\}_{s\ge0}\) is a strongly continuous semigroup on \(X=C([0,L])\), it is uniformly bounded on the compact time interval \([0,T]\). Thus, there exists \(M_T>0\) such that \[\|T(s)f\|_\infty \le M_T\|f\|_\infty, \qquad f\in X,\quad 0\le s\le T.\] Hence \[\|v_{ss}(s,\cdot)\|_\infty = \|T(s)A^2\phi\|_\infty \le M_T\|A^2\phi\|_\infty, \qquad 0\le s\le T.\] Hence, \(v_{ss}\) is uniformly bounded, \[\|v_{ss}(s,\cdot)\|_\infty \le M_T \|A^2\phi\|_\infty, \qquad 0\le s\le T.\] Next, since \(Af=\frac{1}{2} f''\) and \(A^2f=\frac{1}{4} f^{(4)}\) in the interior, we obtain \(v_{xx}\) and \(v_{xxxx}\) are uniformly bounded on \([0,T]\times[0,L]\) by, \[\|v_{xx}(s,\cdot)\|_\infty =2\|Av(s,\cdot)\|_\infty =2\|T(s)A\phi\|_\infty \le 2M_T\|A\phi\|_\infty,\] \[\|v_{xxxx}(s,\cdot)\|_\infty =4\|A^2v(s,\cdot)\|_\infty =4\|T(s)A^2\phi\|_\infty \le 4M_T\|A^2\phi\|_\infty.\]

To control the odd derivatives, we use the boundary conditions and integration. Since \(v(s,\cdot)\in D(A)\), we have \(\partial_x v(s,0)=\kappa\,\partial_{xx}v(s,0)\) and \[|\partial_x v(s,0)|\le \kappa \|\partial_{xx}v(s,\cdot)\|_\infty.\] As a result, \(\partial_x v\) is uniformly bounded, \[\begin{align} |\partial_x v(s,x)|&=|\partial_x v(s,0)|+\left|\int_0^x \partial_{xx}v(s,y)\,dy\right|\\ &\leq (\kappa+L)\|\partial_{xx}v(s,\cdot)\|_\infty. \end{align}\]

Similarly, since \(v(s,\cdot)\in D(A^2)\), we know \(Av(s,\cdot)\in D(A)\). Applying the sticky boundary condition to \(Av=\frac{1}{2} \partial_{xx}v\) and the uniform bound on fourth order derivative, we obtain \(|\partial_{xxx}v(s,0)|\le \kappa \|\partial_{xxxx}v(s,\cdot)\|_\infty\). Hence \(\partial_{xxx}v\) is uniformly bounded by similar reasoning above \[\begin{align} |\partial_{xxx}v(s,x)| &= |\partial_{xxx}v(s,0)|+\left|\int_0^x \partial_{xxxx}v(s,y)\,dy\right|\\ &\leq (\kappa+L)\|\partial_{xxxx}v(s,\cdot)\|_\infty. \end{align}\] We have now shown that \[\sup_{(s,x)\in[0,T]\times[0,L]} |\partial_x^j v(s,x)|<\infty, \qquad j=1,2,3,4,\] and \[\sup_{(s,x)\in[0,T]\times[0,L]} |\partial_{ss}v(s,x)|<\infty.\] Since \(u(t,x)=v(T-t,x)\), the same spatial bounds hold for \(u\), and \(\partial_{tt}u(t,x)=\partial_{ss}v(T-t,x)\). Therefore \[\sup_{(t,x)\in[0,T]\times[0,L]} |\partial_x^j u(t,x)|<\infty, \qquad j=1,2,3,4,\] and \[\sup_{(t,x)\in[0,T]\times[0,L]} |\partial_{tt}u(t,x)|<\infty.\] This completes the proof. ◻

Remark 2. The same argument works if the reflecting boundary condition at \(x=L\) in ?? is replaced by the general Wentzell boundary condition \[\frac{1}{2} \partial_{xx}u(t,L)+\beta_L\partial_x u(t,L)+\gamma_L u(t,L)=0.\] In that case one replaces the domain of \(A\) by \[D(A):=\left\{f\in C^2([0,L]) : f'(0)=\kappa f''(0),\;\frac{1}{2} f''(L)+\beta_L f'(L)+\gamma_L f(L)=0\right\},\] and the compatibility conditions at \(x=L\) become \[\frac{1}{2} \phi''(L)+\beta_L\phi'(L)+\gamma_L\phi(L)=0\] and \[\frac{1}{4} \phi^{(4)}(L)+\frac{\beta_L}{2}\phi^{(3)}(L)+\frac{\gamma_L}{2}\phi''(L)=0.\] The proof is unchanged.

4 General Sticky Diffusion↩︎

Now we consider a general one-dimensional sticky diffusion with drift \(b(x)\) and diffusion \(\sigma(x)\). For ease of reference we recall the generator of the process is \[\label{eq:Lgen} \mathcal{L}u = b \cdot \partial_x u + a\partial_{xx}u,\tag{43}\] with \(a = \frac{\sigma^2}{2}\), and with a sticky boundary condition at \(x=0\): \[\label{eqn:L95diffusion95BC} \left.a\partial_xu\right\vert_{x=0} = \left.\kappa\left(b\partial_xu + a\partial_{xx}u\right)\right\vert_{x=0}\, .\tag{44}\] The backward equation describing the process is \[\label{eqn:diffusion95bckwd95eqn} \partial_t u + \mathcal{L}u = 0, \quad x \in (0,+\infty), \qquad u(T,x) = \phi(x), \quad x \in (0,+\infty),\tag{45}\] with the sticky boundary condition 44 .

4.1 Overview of the simulation algorithm↩︎

The scheme proceeds in the same way as for a SBM: in a boundary layer, we propose a reflected EM jump with probability \(\lambda(\cdot)\), and otherwise we jump directly to the sticky boundary. A difference from a SBM is that the boundary layer is defined as a function of the variable diffusion \(\sigma(x)\), and hence may not be a continuous interval; this introduces some technicalities into the proofs. Further, the procedure for reflecting the process must be modified, to ensure the jump probability \(\lambda(\cdot)\) is indeed a valid probability. As reader will see later, subtleties arise in the scheme when the drift term is negative.

We generate a process \(\{X_k\}_{k=1,\dots,N}\) with time step \(h\). Given \(X_k=x\), the boundary layer and exterior are defined as \[\begin{align} \text{boundary layer:}\quad &\Omega_{b}=\{x: 0\leq \frac{x}{\sigma(x)} \leq \sqrt{3h}\},\label{eq:Omegab}\\ \text{exterior:}\quad & \Omega_e = \{x: \frac{x}{\sigma(x)} > \sqrt{3h}\}.\nonumber \end{align}\tag{46}\]

Remark 3. The reader will observe that the boundary layer as defined above can be disconnected. We’ll handle this by showing that, for small enough \(h\), it is contained in an interval \([0,\sigma_2\sqrt{3h}]\) for some \(\sigma_2>0\). We’ll then do a local error analysis that is valid in this larger interval, hence it will apply to our boundary layer, which is a subset of this interval. Furthermore, the local error in the part of the exterior which is in the interval \([0,\sigma_2\sqrt{3h}]\) will be controlled by standard error analyses for reflected EM jumps.

Given a point \(X_k\) in the exterior, \(X_k \in \Omega_e\), we update \(X_{k+1}\) using a standard Euler-Maruyama increment with reflection:\[\label{eqn:diffusion95cont95scheme95outlayer} X_{k+1} = \left|X_k + \sigma(X_k)Z_k + b(X_k)h\right|.\tag{47}\] Here \(\{Z_k\}\) is a sequence of iid uniform random variables on \([-\sqrt{3h}, \sqrt{3h}]\).

When the particle is within the boundary layer, \(X_k \in \Omega_b\), the scheme performs a type of reflected EM jump with probability \(\lambda(X_k)\), and otherwise it jumps to \(0\) with probability \(1-\lambda(X_k)\): \[\label{eqn:diffusion95cont95scheme95inlayer} X_{k+1} = \begin{cases} \Bigl| \left|X_k + \sigma(X_k)Z_k\right| + b(X_k)h \Bigr|, & \text{with probability } \lambda(X_k), \\ 0, & \text{with probability } 1-\lambda(X_k). \end{cases}\tag{48}\]

Remark 4. The reflected EM jump in 48 involves a nested absolute value. A more standard reflection would take a single absolute value of the entire expression. We attempted this more standard reflection first, but it led to an expression for \(\lambda(x)\) that we were unable to prove is a valid probability.

The jump probability function is \[\label{eqn:diffusion95lambda95cont95scheme} \lambda(x) = \frac{(\sigma^2(x)-2\kappa b(x))x^2 + 2\kappa \sigma^2(x) x + \sigma^4(x)h}{\left(\frac{\kappa\sigma(x)}{\sqrt{3h}} + \sigma^2(x) - 2\kappa b(x)\right)x^2 + \sigma^3(x) \kappa \sqrt{3h} + \sigma^4(x)h}.\tag{49}\] This expression is derived shortly via the method of undetermined coefficients. We prove in Section 4.3 that it is a valid probability for small enough \(h\).

An summary of the simulation scheme for a general sticky diffusion is given in Algorithm [alg:diffusion95cont95scheme].

Set \(X_0=x\). Generate \(Z_{k} \sim \mathcal{U}([-\sqrt{3h}, \sqrt{3h}])\). Calculate boundary threshold \(L_k = \sigma(X_k)\sqrt{3h}\). Update \(X_{k+1}\) using 47 . (\(X_k \in [0, L_k]\)) Calculate \(\lambda(X_k)\) using 49 . Generate \(u \sim \mathcal{U}([0,1])\). Update \(X_{k+1} = \bigl| |X_k + \sigma(X_k)Z_k| + b(X_k)h \bigr|\). Set \(X_{k+1} = 0\).

4.2 How to choose \(\lambda(x)\)↩︎

As for a SBM, we want to choose the jump probability \(\lambda\) so that the transition operator of our discrete scheme is a good approximation to the generator of sticky diffusion in the boundary layer. By Taylor expansion of \(P^h f - f\) and the generator \(\mathcal{L}\) defined in 43 , we want to find \(\lambda\) such that \[\label{eq:genapprox95diffusion} \left|(\partial_xf)\mathbb{E}_x(\Delta X) + \frac{1}{2}(\partial_{xx}f) \mathbb{E}_x(\Delta X^2) -hb(x)(\partial_xf) - ha(x)(\partial_{xx}f)\right| \leq \epsilon,\tag{50}\] for some sufficiently small \(\epsilon\); as before we’ll have \(\epsilon \leq C(h^2+hx)\). To replace \(\partial_xf\) by \(\partial_{xx}f\), we Taylor-expand both sides of the sticky boundary condition 44 : \[\begin{align} \sigma^2(0) \partial_xf(t,0) &= \sigma^2(x)\partial_x f(t,x) - x\sigma^2(x)\partial_{xx}f(t,x) + \mathcal{O}(h+x),\\ 2\kappa b(0)\partial_x f(t,0) + \kappa \sigma^2(0)\partial_{xx}f(t,0) &= 2\kappa b(x)\partial_x f(t,x) +\\& \qquad \left(- 2x\kappa b(x) + \kappa \sigma^2(x)\right)\partial_{xx}f(t,x) + \mathcal{O}(h+x). \end{align}\] Substituting these expansions back into the boundary condition 44 , and dropping the error terms, gives \[\label{eq:taylorsticky95diffusion} \begin{align} \partial_x f(t,x) &= \left(M(x) + x \right)\partial_{xx}f(t,x) + \mathcal{O}(h+x), \end{align}\tag{51}\] where \[M(x) = \frac{\kappa \sigma^2(x)}{\sigma^2(x)-2\kappa b(x)}.\] Approximation 51 assumes \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\) for \(h\) sufficiently small (if we assume continuity of \(b(x)\) and \(\sigma(x)\)). We will revisit the special case where this doesn’t hold when we do the local error analysis. Substituting 51 into 50 , dropping the remainder terms gives \[\label{eq:genapprox295diffusion} (M(x)+x)\left(\mathbb{E}_x(\Delta X)-hb(x)\right) + \frac{1}{2}\mathbb{E}_x(\Delta X^2) - ha(x) = \epsilon.\tag{52}\] We may now calculate \(\mathbb{E}_x(\Delta X)\), \(\mathbb{E}_x(\Delta X^2)\), as functions of \(\lambda\), and then solve 52 for \(\lambda\) (after setting \(\epsilon=0\)). For \(x-\sigma(x)\sqrt{3h}<0\), readers can verify that \[\begin{align} \mathbb{E}_x(\Delta X_{k}) &= \frac{\lambda x^2}{2\sqrt{3h} \cdot \sigma(x)} + \frac{\sigma(x) \lambda \sqrt{3h}}{2} + \lambda b(x)h - x + \mathcal{O}(hx + h^2)\,, \nonumber\\ \mathbb{E}_x(\Delta X_k^2) &= -\frac{\lambda}{\sigma(x)\sqrt{3h}}\cdot x^3 + (1+\lambda)x^2 - \lambda\sigma(x)\sqrt{3h}\cdot x + \lambda\sigma^2(x)\cdot h + \mathcal{O}(hx + h^2)\,. \label{eq:DelX95diffusion} \end{align}\tag{53}\] Substituting into 52 , and given \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), we have \[\lambda\cdot\left( \left( \frac{M}{\sigma\sqrt{3h}} + 1 \right)x^2 + M(x)\sigma\sqrt{3h} + \sigma^2 h + 2bMh \right) = \sigma^2 h + x^2 + 2Mx + 2bMh.\] Solving for \(\lambda\) gives 49 .

4.3 \(\lambda\) is a probability↩︎

To show that \(\lambda\) defined in 49 is a valid probability for small enough \(h\), we must impose assumptions on the drift and diffusion coefficients. The first is a local regularity assumption, which is used for Taylor expansions and coefficient bounds once the boundary layer has been localized near \(0\):

Assumption 3 (Local coefficient regularity near the sticky boundary). There exists \(r_0>0\) such that \(b,\sigma\in C^2([0,r_0])\), \(\sigma(0)>0\), and \[C_{\mathrm{loc}} := \max_{0\le j\le 2} \left( \sup_{x\in[0,r_0]} |\partial_x^j b(x)| + \sup_{x\in[0,r_0]} |\partial_x^j \sigma(x)| \right) <\infty.\] After decreasing \(r_0\) if necessary, we also assume that \[\sigma(x)\ge \sigma_* >0,\qquad x\in[0,r_0].\]

The second is a linear-growth assumption, which has two roles: it first ensures that the state-dependent boundary layer cannot extend far from the sticky boundary when \(h\) is small, and it later provides the standard moment control needed for the Euler–Maruyama part outside the boundary layer.

Assumption 4 (Linear growth of the coefficients). There exists \(L>0\) such that \[|b(x)|+|\sigma(x)|\le L(1+x),\qquad x\ge 0.\]

The following lemma shows that, once \(h\) is small, every point in the boundary layer lies in a compact interval on which \(b(x)\in [b_1,b_2]\) and \(\sigma(x)\in [\sigma_1,\sigma_2]\). In particular, \(\sigma_1\) and \(\sigma_2\) are simply uniform lower and upper bounds for the diffusion coefficient on the region when \(x \in \Omega_{br}\).

Lemma 8. Under Assumptions 3 and 4, there exists \(\delta>0\) and constants \(b_1,b_2,\sigma_1,\sigma_2\), independent of \(x\) and \(h\), with \(\sigma_1>0\), such that, for every \(h<\delta\), \[b(x)\in [b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2],\] whenever \[x-\sigma(x)\sqrt{3h}<0 \qquad\text{or}\qquad x-\sigma(x)\sqrt{3h}+b(x)h<0.\]

Proof. See Supplement. ◻

As a consequence, whenever \(h<\delta\), for any \(x\in \Omega_b\), we may apply the fixed boundary-layer bounds \[b(x)\in[b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2].\] Moreover, \[x<\sigma(x)\sqrt{3h}\le \sigma_2\sqrt{3h},\] so every point in the boundary layer is in fact \(O(\sqrt h)\) close to the sticky boundary, exactly as in the SBM case where the boundary layer is explicit.

Lemma 9. Under Assumptions 3 and 4, there exists \(\delta>0\) such that, for every \(h<\delta\) and every \(x\in \Omega_b\) the jump probability 49 satisfies \[\lambda(x)\in[0,1].\]

Proof. See Supplement. ◻

4.4 Primary convergence results↩︎

Now we present our main convergence results. We wish to show weak convergence of the numerical scheme with order \(h\) and \(X_0=x_0\), i.e. that \[|\mathbb{E}\phi(X_N) - \mathbb{E}\phi(Y_T)| \leq Ch\] for some constant \(C\), where \(\phi\) is from a suitable class of functions. Similar to our previous discussion in SBM, we let \(u(t,x)= \mathbb{E}[\phi(Y_T) | Y_t=x]\) be a function solving the backward equation for the sticky diffusion, 45 . Writing \(t_k=kh\), \(u_k=u(t_k,X_k)\), we have a similar decomposition of the error as in 24 : \[\begin{align} |\mathbb{E}\phi(X_N) - \mathbb{E}\phi(Y_T)| &\leq \mathbb{E}\left(\sum_{k=0}^{N-1} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k\right)\right|\cdot \mathbf{1}_{X_k \in \Omega_{br}}\right) \nonumber\\&\quad + \mathbb{E}\left(\sum_{k=0}^{N-1} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k\right)\right|\cdot \mathbf{1}_{X_k \notin \Omega_{br}}\right). \label{eq:err95decompose95diffusion} \end{align}\tag{54}\] Here \(\Omega_{br} = \Omega_{b} \cup \Omega_{r}\) with \[\begin{align} \Omega_{r} &:= \{x > \sigma(x)\sqrt{3h};\; x-\sigma(x)\sqrt{3h}+b(x)h\leq 0\}.\label{eq:Omegabi} \end{align}\tag{55}\] We include \(\Omega_{r}\) along with \(\Omega_b\) because both regions have the same form of local error, and we refer to \(\Omega_{br}\) henceforth as the boundary layer.

To present main results in the error analysis, we need the same assumptions we have for SBM and Assumptions 3, 4. Suppose that \(h\) is always chosen such that \(N := T/h \in \mathbb{Z}\).

Theorem 5. (Local Error of \(\{X_k\}_{k\in\mathbb{N}}\) in Algorithm [alg:diffusion95cont95scheme])
Under Assumptions 1, 2, 3, and 4, choose \(\epsilon=\mathcal{O}(hx+h^2)\) in 52 . Then there exists a constant \(M>0\), independent of \(x\) and \(h\), depending on the relevant derivative bounds of \(u\), the local boundary constants for \(b,\sigma\), the linear-growth constant of \(b,\sigma\), and \(\kappa\), such that for all sufficiently small \(h\), \[\label{ogkzhdpc} \left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k=x\right)\right| \leq \begin{cases} M(h^2+hx), & x\in \Omega_{br},\\ Mh^2(1+x^4), & x\in \Omega_e. \end{cases}\tag{56}\]

Hence, as for a SBM, the local error in the boundary layer is \(O(h^{3/2})\), while that in the exterior is \(O(h^2)\). Also as for a SBM, we will show the number of steps in the boundary layer is \(O(h^{-1/2})\), leading to a global error of \(O(h)\). We have \(Mh^2(1+x^4)\) because we analyze the local error outside the boundary layer by Taylor expansion up to fourth order, which contains products of drift and diffusion of order 4 and we apply the linear growth condition. Analogous to Corollary 1 and Lemma 4, we have

Corollary 2. Under the assumptions of Theorem 5, \[\mathbb{E}\left(\sum_{k=0}^{N-1} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k\right)\right|\cdot \mathbf{1}_{X_k \in \Omega_{br}}\right) \leq M\left(h + h \cdot \mathbb{E}\left(\sum_{k=0}^{N-1} X_k \cdot \mathbf{1}_{X_k \in \Omega_{br}} \right)\right).\]

Lemma 10. For sufficiently small \(h\), there exists a positive constant \(C\), independent of \(h\) and initial condition \(x_0\), such that \[\mathbb{E}\left(\sum_{k=0}^{N-1} X_k \cdot \mathbf{1}_{X_k \in \Omega_{br}} \right) < C.\]

The proof will recycle most of the construction we used to prove Lemma 4. If \(\sigma(0) = 1\), we can show the jump probability \(\lambda\) and the transition operator \(P^h\) for SBM and the sticky diffusion are sufficiently close.

The following lemma helps us analyze the error coming from the exterior where the scheme does a standard EM jump where we will transform the drift and diffusion into \((1+X_k)^4\) from the linear growth assumption.

Lemma 11. Under assumptions 3, 4, and suppose that \(X_0=x_0\) deterministically. Then, for \(T=Nh\), there exists a constant \(C_T>0\), independent of \(h\) and \(x_0\), such that \[\sup_{0\le k\le N}\mathbb{E}(1+X_k^4)\le C_T(1+x_0^4).\]

Applying these Lemmas to 54 gives the global error bound for the sticky diffusion, our main convergence theorem.

Theorem 6. (Global Error of \(\{X_k\}_{k\in\mathbb{N}}\) in Algorithm [alg:diffusion95cont95scheme])
Let \(\phi:[0,\infty)\to\mathbb{R}\) be such that the associated backward solution \[u(t,x):=\mathbb{E}\!\left[\phi(Y_T)\mid Y_t=x\right]\] satisfies Assumptions 1, 2, 3, and 4. Assume moreover that \(X_0=Y_0=x_0\) deterministically. Then there exist constants \(\delta>0\) and \(C>0\), independent of \(h\) and \(x_0\), such that for all \(h<\delta\), \[\label{rhoatlzi} \left|\mathbb{E}\bigl(\phi(X_N)-\phi(Y_T)\bigr)\right| \le C(1+x_0^4)h.\tag{57}\]

Proof. This follows from the error decomposition 54 . By Theorem 5, there exists a constant \(M>0\), independent of \(x\) and \(h\), such that \[\left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k=x\right)\right| \le \begin{cases} M(h^2+hx), & x\in\Omega_{br},\\ Mh^2(1+x^4), & x\in\Omega_e. \end{cases}\] Applying this bound to 54 , we obtain \[\begin{align} \left|\mathbb{E}\phi(X_N)-\mathbb{E}\phi(Y_T)\right| &\le M\mathbb{E}\left(\sum_{k=0}^{N-1}(h^2+hX_k)\mathbf{1}_{\{X_k\in\Omega_{br}\}}\right)\\ &\quad+ Mh^2\sum_{k=0}^{N-1}\mathbb{E}\left((1+X_k^4)\mathbf{1}_{\{X_k\in\Omega_e\}}\right). \end{align}\] For the boundary-layer contribution, Lemma 10 gives \[\mathbb{E}\left(\sum_{k=0}^{N-1}X_k\mathbf{1}_{\{X_k\in\Omega_{br}\}}\right)\le C_b,\] where \(C_b\) is independent of \(h\) and \(x_0\). Since \(N=T/h\), we also have \[\mathbb{E}\left(\sum_{k=0}^{N-1}h^2\mathbf{1}_{\{X_k\in\Omega_{br}\}}\right) \le Nh^2 =Th.\] Therefore, \[M\mathbb{E}\left(\sum_{k=0}^{N-1}(h^2+hX_k)\mathbf{1}_{\{X_k\in\Omega_{br}\}}\right) \le (MT+MC_b)h.\] For the exterior contribution, Lemma 11 yields \[\sup_{0\le k\le N}\mathbb{E}(1+X_k^4)\le C_T(1+x_0^4),\] where \(C_T\) is independent of \(h\) and \(x_0\). Hence \[\begin{align} Mh^2\sum_{k=0}^{N-1}\mathbb{E}\left((1+X_k^4)\mathbf{1}_{\{X_k\in\Omega_e\}}\right) &\le Mh^2\sum_{k=0}^{N-1}\mathbb{E}(1+X_k^4)\\ &\le Mh^2N\, C_T(1+x_0^4)\\ &= MTC_T(1+x_0^4)h. \end{align}\]

Combining the two estimates gives \[\left|\mathbb{E}\phi(X_N)-\mathbb{E}\phi(Y_T)\right| \le (MT+MC_b)h + MTC_T(1+x_0^4)h.\] Since \(1\le 1+x_0^4\), we may absorb the first term into the second and obtain \[\left|\mathbb{E}\phi(X_N)-\mathbb{E}\phi(Y_T)\right| \le C(1+x_0^4)h\] for some constant \(C>0\) independent of \(h\) and \(x_0\). ◻

Figure 4: Weak error e_h = |\mathbb{E}[X_T\mid x_0=0]-u(T,0)| of Algorithm [alg:diffusion95cont95scheme] for different choices of drift, diffusion, and sticky parameter: (a) b(x)=2+\sin(x), \sigma(x)=\frac{1}{2}+x, \kappa=4; (b) b(x)=2+\sin(x), \sigma(x)=2+x, \kappa=1; (c) b(x)=-1+\sin(x), \sigma(x)=2+x, \kappa=1; (d) b(x)=-2+\sin(x), \sigma(x)=\frac{1}{2}+x, \kappa=4. All experiments use L=4, T=1, and x_0=0. Error bars (commensurate with marker sizes) are \pm 1 empirical standard deviation. The respective least-squares fitted slopes are (a) 0.954162 (b) 1.026675, (c) 1.010317, (d) 0.907653 (fitted only through the last 2 data points in this case). — Figure 4: Weak error \(e_h = |\mathbb{E}[X_T\mid x_0=0]-u(T,0)|\) of Algorithm [alg:diffusion95cont95scheme] for different choices of drift, diffusion, and sticky parameter: (a) \(b(x)=2+\sin(x)\), \(\sigma(x)=\frac{1}{2}+x\), \(\kappa=4\); (b) \(b(x)=2+\sin(x)\), \(\sigma(x)=2+x\), \(\kappa=1\); (c) \(b(x)=-1+\sin(x)\), \(\sigma(x)=2+x\), \(\kappa=1\); (d) \(b(x)=-2+\sin(x)\), \(\sigma(x)=\frac{1}{2}+x\), \(\kappa=4\). All experiments use \(L=4\), \(T=1\), and \(x_0=0\). Error bars (commensurate with marker sizes) are \(\pm 1\) empirical standard deviation. The respective least-squares fitted slopes are (a) \(0.954162\) (b) \(1.026675\), (c) \(1.010317\), (d) \(0.907653\) (fitted only through the last 2 data points in this case).

4.5 Numerical experiments↩︎

Figure 4 empirically validates our convergence results using different choices for \(b,\sigma\) as given in the caption. As for a SBM, we impose a reflecting boundary condition at \(x=L\), in order to compare our stochastic estimates with a “ground truth” constructed by solving the backward equation numerically using finite difference with a very small time step. Also as before we choose \(\phi(x)=x\), which does not satisfy the boundary condition compatibility conditions that are typically required in PDE proofs of smoothness. Figure 4 (a,b,c) used timesteps \(h\in\{1/128,1/256,1/512,1/1024,1/2048\}\), and a number of realizations which started at \(5\times 10^7\) and was quadrupled each time \(h\) was halved. Figure 4 (d) used \(h\in\{1/1024,1/2048,1/4096,1/8192\}\), with initial number of realizations equal to \(10^9\).

All cases approach \(O(h)\) convergence, though it is notable that case (d) with strongly negative drift and weaker diffusion approaches it more slowly, requiring must smaller timesteps to approach a slope of 1. The error however in this case is several orders of magnitude smaller than for the others.

4.6 Key steps in the error analysis↩︎

We now show the key steps in performing a local error analysis, Theorem 5, and in proving Lemma 10. Compared to an SBM, an additional challenge arises when we have negative drift \(b(X_k)\), where our scheme does a reflected EM jump when \(x - \sigma(x)\sqrt{3h} + b(x)h< 0\). It is non-trivial to show this situation gives desired order of local accuracy. Compared to the SBM case, the challenge is that the near-boundary region is no longer exactly \([0,\sqrt{3h}]\). The variable diffusion changes the boundary threshold, and negative drift can create an additional thin reflection region. Thus, for the number-of-steps lemma, we first need to show that this enlarged region is still contained in an \(O(\sqrt h)\) neighborhood of the boundary, and then compare the resulting transition operator with that for an SBM.

This section collects the intermediate statements behind the convergence results in Section 4.4. All the proofs use standard techniques and are found in the Supplement.

4.6.1 Local error: Theorem 5 ↩︎

The role of this section is the same as that of Section 3.4 for sticky Brownian motion: we first isolate the local error estimate in the boundary layer, then record the occupation bound needed in the global error argument, and finally return to the weak error decomposition 54 . Unlike in Section 3.4, the error in \(\Omega_e\) is less straightforward because we do not want to impose uniform boundedness assumptions on the drift and diffusion. Instead, we assume linear growth and local regularity in the boundary layer.¹ The Euler-Maruyama jump and its local weak error is of order \(h^2\).

We start by concentrating on \(\Omega_{br}\). Recall this was divided into two separate pieces \(\Omega_{b},\Omega_{r}\) defined in 46 ,55 . There are three situations to treat separately:

\(x\in \Omega_{b}\) and \(b(0)\neq \frac{\sigma^2(0)}{2\kappa}\);
\(x\in \Omega_{r}\), which is the additional reflection strip created by negative drift;
\(x\in \Omega_{b}\) and \(b(0)= \frac{\sigma^2(0)}{2\kappa}\) (notice \(\Omega_{br} = \Omega_{b}\) and \(\Omega_{r} = \emptyset\) because\(b(x)>0\) for small \(h\)).

All lemmas and theorems below impose Assumptions 3-4 and the proofs can be found in the Supplements.

Similar to the proofs for SBM, we proceed to show that \(\lambda(x)\) is upper bounded by a multiple of \(\frac{x}{\sqrt{h}}+\sqrt{h}\). This not only tells us the jump probability to the sticky boundary is higher when the particle position \(x\) is close to the sticky boundary, but it also enables us to conclude that \(\lambda(x)\cdot h^{3/2} = \mathcal{O}(hx+h^2)\) in the local error analysis.

Lemma 12. There exists constant \(C, \delta > 0\), depending on \(\sigma\) and \(b\), such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), and for \(\lambda\) defined in 49 , \[0 \leq \lambda(x) \leq S_1\cdot \left(\frac{x}{\sqrt{h}} + \sqrt{h} \right).\]

Now we would like to prove a couple of lemmas to prepare for showing that our diffusion scheme has local error of order \(3/2\). We begin with the generic case \(x\in \Omega_{b}\) and \(b(0)\neq \frac{\sigma^2(0)}{2\kappa}\). The first lemma is the analogue of Lemma 2 from the SBM proof. It rewrites the sticky boundary condition in a form that lets us replace \(\partial_xu\) by \(\partial_{xx}u\) up to a remainder whose size is exactly what the local error proof can absorb.

Lemma 13. Let \(u\) be the solution of 45 satisfying corresponding boundary condition 44 and assume that \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\). Then, there exists \(\delta>0\) such that for all \(h<\delta\) and \(x\in \Omega_{b}\), \[\label{eqn:diffusion95BC95consequence} \begin{align} \partial_x u(t,x) &= \left(x+M(x) \right)\partial_{xx}u(t,x) + Re(x) \end{align}\tag{58}\] where \[|Re(x)|\leq C(h+x)=\mathcal{O}(h+x),\] and \(C\) depends on bounds of \(\partial_xu,\partial_{xx}u,u^{(3)}\) and on the local bounds of \(b,\sigma\) and their derivatives up to order \(2\) on the boundary neighborhood \([0,r_0]\), but is independent of \(x\) and \(h\).

The general idea of the proof is to Taylor expand both sides of the boundary condition 44 and bound remainder terms by regularity assumptions.

The following two lemmas record the first and second conditional moments of the increment in \(\Omega_{b}\). They are the diffusion analogues of the moment identities for SBM. The first and second moments both have an extra remainder term when the drift is negative because of the nested absolute value; we apply Lemma 12 to show that this term is of order \(hx+h^2\).

Lemma 14. There exists \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), we have, \[\label{eqn:first95mom95negative95drifted95diffusion95approx} \begin{align} \mathbb{E}_x(\Delta X_{k}) &= \frac{\lambda(x)x^2}{2\sqrt{3h} \cdot \sigma(x)} + \frac{\sigma(x) \lambda(x)\sqrt{3h}}{2} + \lambda(x)b(x)h - x + \mathcal{O}(hx + h^2).\\ \end{align}\tag{59}\]

Lemma 15. There exists \(\delta>0\) and \(C > 0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), we have, \[\mathbb{E}_x(\Delta X_k^2) = -\frac{\lambda(x)}{\sigma(x)\sqrt{3h}}\cdot x^3 + (1+\lambda(x))x^2 - \lambda(x)\sigma(x)\sqrt{3h}\cdot x + \lambda(x)\sigma^2(x)\cdot h + R_2(x),\] where \[|R_2(x)| \leq C\cdot (hx + h^2) = \mathcal{O}(hx + h^2).\]

The next lemma is the exact analogue of the moment identity used in Lemma 3. Its purpose is to isolate the leading-order cancellation built into the choice of \(\lambda(x)\). We apply the computation done in Lemma 14-15. They produce the remainder term in the following lemma.

Lemma 16. Assume that \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), there exist \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\) \[\label{eqn:1st952nd95mom95goal95cont95scheme95drifted} 2\left(\mathbb{E}_x(\Delta X_k) - hb(x)\right) \cdot (x + M(x)) + \mathbb{E}_x(\Delta X_k^2) = \sigma^2(x)h + R_3(x)\tag{60}\] where \(|R_3(x)| \leq K_3 \cdot (hx + h^2) = \mathcal{O}(hx + h^2)\) for some \(K_3\) independent of \(h\) but depending on the boundary-layer constants and the nondegeneracy constant from \(b(0)\neq\sigma^2(0)/(2\kappa)\).

To complete the local error proof we still need bounds on the first and third conditional moments themselves. The first-moment estimate is not immediate, because jumps to \(0\) are of size \(\mathcal{O}(\sqrt{h})\). Moreover, we will estimate first moment depending on \(x\in \Omega_{b}\) and \(x\in \Omega_{r}\) separately. Exactly as in the SBM argument, the previous cancellation lemma is what lowers the first moment to order \(h\).

Lemma 17. There exists \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), \[\label{eqn:first95mom95negative95drifted95bound95diffusion} \left|\mathbb{E}_x(\Delta X_k)\right| \leq C(h + hx+h^2),\tag{61}\] where \(C>0\) is independent of \(x\) and \(h\) and depends only on \(\kappa\), the local boundary bounds for \(b,\sigma\), and the nondegeneracy constant coming from \(b(0)\neq \sigma^2(0)/(2\kappa)\).

Lemma 18. There exists \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{br}\), \[\label{eqn:first95mom95drifted95bound} \left|\mathbb{E}_x(\Delta X_k^3)\right| \leq K_4(hx + h^2),\tag{62}\] where \(K_4\) is a fixed constant dependent on drift, diffusion, and \(\kappa\).

We now prove the local error estimate in \(\Omega_{b}\) for the generic case \(b(0)\neq \frac{\sigma^2(0)}{2\kappa}\). The proof follows the same Taylor-expansion template as in Section 3.5, but the cancellation now uses Lemma 13 and Lemma 16.

Theorem 7. (Local Error of \(\{X_k\}_{k\in\mathbb{N}}\) in algorithm [alg:diffusion95cont95scheme], situation 1)
Given \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), there exists \(\delta > 0\) and \(C\geq 0\) depending on the relevant bounds of derivatives of \(u\) and on the local boundary constants for \(b,\sigma\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), \[\label{eqn:local95err95cont95scheme95drifted} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| \leq C(h^2 + hx) = \mathcal{O}(hx + h^2)\tag{63}\]

Now, as mentioned before, if \(b(x) < 0\), then for some \(x\in \Omega_{r} \subseteq [0,\sigma_2\sqrt{3h}]\), we may have \(x - \sigma(x)\sqrt{3h} + b(x)h < 0\) and \(x - \sigma(x)\sqrt{3h} > 0\). According to our scheme (Algorithm [alg:diffusion95cont95scheme]), we do a reflected EM jump, and there is zero probability to jump to the sticky boundary. It is non-trivial to show this will give the desired order of local error because the first conditional moment of \(\Delta X_k\) might not be a good approximation of \(b(x)h\). The next lemma aims to resolve this doubt.

Lemma 19. There exists \(\delta > 0\) and a constant \(C\), which depends only on local boundary constants for \(b,\sigma\) coming from Lemma 8 such that \(\forall h < \delta\) and \(\forall x\in \Omega_{r}\), \[|E_x(\Delta X_k) - b(x)h| \leq C( hx + h^2) = \mathcal{O}(hx + h^2).\]

It is not hard to further compute that the second conditional moment of \(\Delta X_k\) is approximately \(\sigma^2(x)h\) with a remainder term of order \(hx + h^2\). Now, we are ready to show that the local error for \(x\) such that \(x - \sigma(x)\sqrt{3h} + b(x)h < 0\) and \(x - \sigma(x)\sqrt{3h}>0\),

Theorem 8. (Local Error of \(\{X_k\}_{k\in\mathbb{N}}\) in algorithm [alg:diffusion95cont95scheme], situation 2)
Given \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), there exists \(\delta > 0\) and \(C \geq 0\) depending on the relevant derivative bounds of \(u\) and on the local boundary constants for \(b\) and \(\sigma\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{r}\), \[\label{eqn:local95err95cont95scheme95drifted95ring} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| \leq C(h^2 + hx) = \mathcal{O}(hx + h^2).\tag{64}\]

So far, we have shown the order of local error for our scheme when \(x - \sigma(x)\sqrt{3h} < 0\) (where we need to consider absorption probability \(\lambda(x)\)) and \(\sigma(x)\sqrt{3h} \leq x < \sigma(x)\sqrt{3h} - b(x)h\) (where the negative drift can cause a reflected EM jump). Notice for the first situation, we assume that \(\sigma(0) \neq \frac{\sigma^2(0)}{2\kappa}\) so that \(M(x)\) can be a bounded function within the boundary layer. We need to compute the local error of our scheme when \(b(0) = \frac{\sigma^2(0)}{2\kappa}\) separately. We will show, similar to situation 2, that the first conditional moment of \(\Delta X_k\) is a good approximation of \(b(x)h\) within the boundary layer.

Lemma 20. When \(b(0) = \frac{\sigma^2(0)}{2\kappa}\), there exist \(\delta>0\) and a constant \(C\) which depends only on the boundary-layer bounds, such that \(\forall h< \delta\) and \(\forall x\in \Omega_{b}\), \[\left| \mathbb{E}_x\left(\Delta X_k\right) - hb(x) \right| \leq C(hx + h^2) = \mathcal{O}(hx + h^2).\]

We now proceed to show the order of local error for case 3, where \(x-\sigma(x)\sqrt{3h} < 0\) and \(b(0) = \frac{\sigma^2(0)}{2\kappa}\). This case must have a positive drift in \(\Omega_{br}\) for \(h\) small enough by its regularity assumption so there is no extra subtlety caused by negative drift.

Theorem 9. (Local Error of \(\{X_k\}_{k\in\mathbb{N}}\) in algorithm [alg:diffusion95cont95scheme], situation 3)
If \(b(0)= \frac{\sigma^2(0)}{2\kappa}\), there exist \(\delta > 0\) and \(C\geq 0\), independent of \(x\) and \(h\), depending only on the relevant bounds of derivatives of \(u\) and on the local boundary constants for \(b,\sigma\) such that \(\forall h<\delta\) and for all \(x\) such that \(x \in \Omega_{br}\), \[\label{eqn:local95err95cont95scheme95drifted95speical} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| \leq C(h^2 + hx) = \mathcal{O}(hx + h^2).\tag{65}\]

4.6.2 Bounding the number of steps: Lemma 10 ↩︎

To enable comparison with the SBM construction, we first normalize the boundary diffusion coefficient. Let \(\sigma_0:=\sigma(0)>0\) and define the rescaled process \[\hat{X}_t:=\frac{X_t}{\sigma_0}.\] Then \(\hat{X}\) is again a sticky diffusion of the same form, with coefficients \[\hat{b}(y):=\frac{b(\sigma_0 y)}{\sigma_0}, \qquad \hat{\sigma}(y):=\frac{\sigma(\sigma_0 y)}{\sigma_0},\] and sticky parameter \(\hat{\kappa}=\kappa/\sigma_0\). In particular, \(\hat{\sigma}(0)=1\). Therefore it is enough to prove the desired estimate in the normalized case \(\sigma(0)=1\), and we adopt this normalization below.

Our hope is to recycle our calculation for the global error analysis done in the proof of Theorem 3 to show similar theorems for the scheme with drift.

The first thing we notice is that when \(\sigma(x) = 1\), then the jump probability \(\lambda\) for the general diffusion, 49 , is \[\lambda(x) = \frac{(1-2\kappa b(x))x^2 + 2\kappa x + h}{\left(\frac{\kappa}{\sqrt{3h}}+1-2\kappa b(x)\right)x^2 + h + \kappa\sqrt{3h}},\] which differs from the jump probability for the SBM, 14 , only in terms involving \(b(x)\). Since these terms multiply \(x^2\), which is small in the boundary layer, and since if \(\sigma(0) = 1\) then \(\sigma(x)\approx 1\) in the boundary layer, we expect that these two jump probabilities will be not too different. The next lemma shows that the difference between the two is \(\mathcal{O}(\sqrt{h})\).

Lemma 21. Denote \(\tilde{\lambda}(\cdot)\) as the jump probability for SBM, 14 , and \(\lambda(\cdot)\) as the jump probability for the general diffusion, [alg:diffusion95cont95scheme]. If \(\sigma(0) =1\), then \[\left|\lambda(x) - \tilde{\lambda}(x)\right| \leq K_{b,\sigma} \sqrt{h},\quad \forall x\in [0,\sqrt{3h}],\] where \(K_{b,\sigma} > 0\) is a constant depending on the drift term \(b(\cdot)\), diffusion \(\sigma(\cdot)\), and sticky parameter \(\kappa\).

These results may be used to show Lemma 10, which is analogous to Lemma 4, and is the final step in proving the order of global error. The proof of Lemma 10 is in the Supplement.

5 Conclusion↩︎

This paper presented and analyzed a numerical scheme for simulating one-dimensional sticky diffusions. The scheme combines a reflected Euler-Maruyama increment, with a carefully designed jump probability to the sticky point, to handle the challenge of moving both in continuous space, and to a lower-dimensional boundary. Our scheme achieves a weak local error of order 3/2 within the boundary layer and a weak global error of order 1.

We provide rigorous theoretical analysis of our schemes, including detailed proofs of local and global error bounds. Numerical experiments corroborate our theoretical findings, showing empirical evidence of first-order global convergence for both sticky Brownian motion and sticky diffusion processes under various conditions.

A natural extension of our work would be to develop schemes for multidimensional sticky diffusions. A novel features of such diffusions, compared to the one-dimensional case, is that they can have nontrivial dynamics along their sticky boundaries, which can be different from the dynamics in the interior. We are optimistic that the method we have outlined can be adapted to such diffusions, because it provides a clear guideline for choosing how to switch between different types of dynamics (i.e. on sets of different dimensions), which are otherwise standard to simulate. A challenge however will be ensuring that the switching probabilities are indeed probabilities, as we have found this nontrivial to ensure even for a one-dimensional diffusion. Further extensions could include adapting these methods to processes on manifolds [43], and considering multiple intersecting sticky boundaries which form a stratification, as outlined in [17], [50].

Even for a one-dimensional diffusion there is further theoretical development possible. One could investigate higher-order schemes, which might achieve even better convergence rates, or could ask how to solve general parabolic-type PDEs with sticky boundary conditions [39], [51]. One could also ask for worse schemes, which achieve only \(O(h^{1/2})\) convergence – perhaps there is a wider scope of jump probabilities \(\lambda\) which allow for worse convergence, but are easier to either derive, or to implement. Such analysis could perhaps be applied to Monte-Carlo schemes which perform local moves to sample from a stationary distribution supported on manifolds of different dimensions, and show their convergence to a particular sticky diffusion [50]. Finally, it would be useful to perform a theoretical analysis of the long-time behavior of the numerical schemes and their ability to capture stationary distributions accurately.

In conclusion, this work contributes to the growing toolkit for numerical simulation of stochastic processes with complex boundary condition, offering both practical algorithms and theoretical guarantees for simulating certain sticky processes. As the study of such processes continues to grow in importance across various disciplines, we believe these methods will prove valuable tools for researchers and practitioners alike.

M.H.C. acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC), RGPIN-2023-04449 / Cette recherche a été financée par le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG).

6 Supplement↩︎

6.1 Additional proofs from Section 3↩︎

Verification of 35 in the proof of Lemma 4

As a reminder, we wish to show that there exist constants \(C_1(s,M)>0\) and \(\delta>0\), independent of \(x\) and \(h\), such that for all \(h<\delta\), \[Pe^{w(x)} \leq e^{w(x)}\left(1-f_h(x)+C_1h\right),\] where \(f_h\) satisfies \[f_h(x)\geq x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}}.\] We would like to calculate \(Pe^{w(x)}-e^{w(x)}\) for \(x\) in three cases: \([0,\frac{\sqrt{3h}}{2}]\), \([\frac{\sqrt{3h}}{2}, \sqrt{3h}]\) and \([\sqrt{3h},+\infty)\).

For \(x\in [0, \frac{\sqrt{3h}}{2}]\), we compute \[\begin{align} Pe^{w(x)} &= \frac{\lambda(x)}{2}e^M + \frac{\lambda(x)}{\sqrt{3h}}\cdot\int_{\frac{\sqrt{3h}}{2}}^{\sqrt{3h}-x} e^{M - s\left(z-\frac{\sqrt{3h}}{2}\right)}\; dz + \frac{\lambda(x)}{2\sqrt{3h}} \cdot \int_{\sqrt{3h}-x}^{\sqrt{3h}+x} e^{M- s\left(z-\frac{\sqrt{3h}}{2}\right)}\; dz \\&\quad + e^M(1-\lambda(x))\\ &= \lambda(x) \left(\frac{1}{2}e^M + e^M \cdot \frac{1}{2\sqrt{3h}s} \cdot \left( 2 - e^{-s\left(\frac{\sqrt{3h}}{2}-x\right)} - e^{-s\left(\frac{\sqrt{3h}}{2}+x\right)} \right)\right) + e^M(1-\lambda(x)). \end{align}\] In the first line, the first term corresponds to reflected jumps that land in the flat part of \(w\), the two integral terms correspond to reflected jumps that land in the sloping part of \(w\) (split according to whether the reflection map has one or two preimages), and the last term corresponds to the jump directly to \(0\).

Now we simplify this expression by Taylor expansion of exponentials with a remainder of order \(\mathcal{O}(h)\): \[\begin{align} \quad&Pe^{w(x)} - e^{w(x)}\\ &= e^M \cdot \frac{\lambda(x)}{2\sqrt{3h}s} \cdot \left( 2 - e^{-s\left(\frac{\sqrt{3h}}{2}-x\right)} - e^{-s\left(\frac{\sqrt{3h}}{2}+x\right)} \right) - e^M \cdot \frac{\lambda(x)}{2}\\ &= e^M\cdot \frac{\lambda(x)}{2\sqrt{3h}s} \cdot \left(\sqrt{3h}s - \frac{1}{2}s^2\left(\frac{3h}{2}+2x^2\right) + \mathcal{O}(h^{3/2})\right) - e^M \cdot \frac{\lambda(x)}{2}\\ &= e^M\cdot \frac{\lambda(x)}{2\sqrt{3h}s} \cdot \left(- \frac{1}{2}s^2\left(\frac{3h}{2}+2x^2\right) + \mathcal{O}(h^{3/2})\right)\\ &= e^M\cdot \left( -\frac{s}{4}\cdot \frac{\lambda(x)}{2}\sqrt{3h} - s \cdot \frac{\lambda(x)x^2}{2\sqrt{3h}} + \mathcal{O}(h)\right). \end{align}\]

In other words, we have shown that, for \(x \in [0,\sqrt{3h}/2]\), \[Pe^{w(x)} \leq e^{w(x)}\left(1-f_h(x) + \mathcal{O}(h)\right)\] where we define \[\label{eqn:fh95in95095to95sqrt403h41472} f_h(x):= \frac{s}{4}\cdot\frac{\lambda(x)}{2}\sqrt{3h} + s \cdot \frac{\lambda(x)\cdot x^2}{2\sqrt{3h}} \qquad \text{for } x\in [0, \frac{\sqrt{3h}}{2}].\tag{66}\]

Recall the expression for \(\mathbb{E}_x(\Delta X)\), 20 and substitute the expression for \(\lambda(x)\), 14 , \[\begin{align} \mathbb{E}_x(\Delta X_k) &=-x + \frac{\lambda(x)\sqrt{3h}}{2} + \frac{\lambda(x) x^2}{2\sqrt{3h}}\\ &= \frac{1}{2\sqrt{3h}}\cdot\frac{\left(x-\sqrt{3h}\right)^2(x^2+h)}{\left(\frac{\kappa}{\sqrt{3h}}+1\right)x^2+h+\kappa\sqrt{3h}}\\ &\geq 0. \end{align}\] Hence, \[f_h(x) - x > \mathbb{E}_x(\Delta X) > 0, \quad \forall x\in [0, \frac{\sqrt{3h}}{2}].\]
For \(x\in [\frac{\sqrt{3h}}{2}, \sqrt{3h}]\), the scheme does the same jump, but the expression for \(e^{w(x)}\) is different, which makes the calculation much more complex. It can be checked that \[\begin{align} Pe^{w(x)} &= \lambda(x) \left(e^M \cdot \frac{\sqrt{3h}-x}{\sqrt{3h}} + e^M \cdot \frac{\frac{\sqrt{3h}}{2}-(\sqrt{3h}-x)}{2\sqrt{3h}}\right)\\ &\quad + \frac{\lambda(x)}{2\sqrt{3h}}\cdot\int_{\frac{\sqrt{3h}}{2}}^{\sqrt{3h}+x}e^{M-s\left(z-\frac{\sqrt{3h}}{2}\right)}\; dz + (1-\lambda(x))e^M\\ &= \lambda(x) e^M \cdot \frac{3\sqrt{3h}-2sx + 2 - 2e^{-s\left(\frac{\sqrt{3h}}{2}+x\right)}}{4\sqrt{3h}s} + (1-\lambda(x))e^M. \end{align}\] Since \(e^{w(x)} = e^{M-s\left(x-\frac{\sqrt{3h}}{2}\right)}\), \[\begin{align} Pe^{w(x)} - e^{w(x)} &= A + B \end{align}\] where \[\begin{align} A&= \lambda(x) e^M \cdot \frac{3\sqrt{3h}-2sx + 2 - 2e^{-s\left(\frac{\sqrt{3h}}{2}+x\right)}}{4\sqrt{3h}s} - \lambda(x) \cdot e^{M-s\left(x-\frac{\sqrt{3h}}{2}\right)}\\ B &= (1-\lambda(x))\cdot e^M - (1-\lambda(x))\cdot e^{M-s\left(x-\frac{\sqrt{3h}}{2}\right)}. \end{align}\] Simplifying \(A,B\) by expanding the exponential, we have \[\label{eqn:A95simplify} \begin{align} A &= \frac{\lambda(x)}{4\sqrt{3h}s} \cdot e^{M-s\left(x-\frac{\sqrt{3h}}{2}\right)} \cdot \left( (3\sqrt{3h}s-2sx+2)e^{s\left(x-\frac{\sqrt{3h}}{2} \right)} - 2e^{-s\sqrt{3h}} - 4\sqrt{3h}s\right)\\ &= e^{w(x)} \cdot \frac{\lambda(x)}{4\sqrt{3h}s} \cdot \left( 3\sqrt{3h}s^2 x - \frac{9}{4}s^2\cdot 3h - s^2x^2 + \mathcal{O}(h^{3/2})\right)\\ &= \frac{\lambda(x)}{4}e^{w(x)} \left( 3sx - \frac{9}{4}\sqrt{3h}s - \frac{sx^2}{\sqrt{3h}} + \mathcal{O}(h) \right). \end{align}\tag{67}\] Meanwhile, \[\label{eqn:B95simplify} \begin{align} B &= e^{M-s\left(x-\frac{\sqrt{3h}}{2}\right)} \cdot (1-\lambda(x)) \left( e^{s\left(x-\frac{\sqrt{3h}}{2}\right)} - 1 \right)\\ &= e^{M-s\left(x-\frac{\sqrt{3h}}{2}\right)} \cdot (1-\lambda(x)) \left( s\left(x-\frac{\sqrt{3h}}{2}\right) + \mathcal{O}(h) \right). \end{align}\tag{68}\] Combining 67 , 68 together, we have \[Pe^{w(x)} - e^{w(x)} = A + B = -e^{w(x)} \left( f_h(x) + \mathcal{O}(h)\right)\] where we define \[\label{eqn:fh95in95sqrt403h4147295to95sqrt403h41} f_h(x) := \frac{s\lambda(x)}{16}\sqrt{3h} + \frac{s\lambda(x)x^2}{4\sqrt{3h}} + \frac{s\lambda(x)x}{4} + \frac{s\sqrt{3h}}{2} - sx\quad \text{for x\in [\frac{\sqrt{3h}}{2}, \sqrt{3h}]}.\tag{69}\] It can be checked by derivative that \(f_h(x)\) is decreasing for \(x\in [\frac{\sqrt{3h}}{2}, \sqrt{3h}]\), and \(f_h(\sqrt{3h}) = \frac{s}{16}\sqrt{3h} > \sqrt{3h}\), so we must have that \[f_h(x) - x > 0, \quad \forall x\in [\frac{\sqrt{3h}}{2}, \sqrt{3h}].\]
For \(x\in (\sqrt{3h},+\infty)\): in this region, the scheme does a standard EM jump. Using the fact that \(\mathbb{E}_x(\Delta X) = 0\), and then using the fact that \(w\) is a piecewise linear function with a constant piece near the boundary, so it smaller than its linear part: \(w(x+\Delta X) \leq M - s\left((x+\Delta X)-\frac{\sqrt{3h}}{2}\right)\), we have \[\begin{align} Pe^{w(x)} &= \mathbb{E}_x\left(e^{w(x+\Delta X)}\right)\\ &\leq \mathbb{E}_x\left(e^{w(x)- s\cdot \Delta X}\right)\\ &= e^{w(x)}\cdot (1 - s\cdot\mathbb{E}_x(\Delta X) + \mathcal{O}(h))\\ &= e^{w(x)}\cdot (1 + 0\cdot (-s) + \mathcal{O}(h))\\ &= e^{w(x)}\left(1 - f_h(x) + \mathcal{O}(h) \right)\\ \end{align}\] where \[\label{eqn:fh95in95sqrt403h4195to95inf} f_h(x) := 0 ,\qquad \text{for x\in (\sqrt{3h},+\infty).}\tag{70}\]

So far, we have shown that for \(x\in [0,+\infty)\), \[\label{eqn:wx95goal95supp} Pe^{w(x)} \leq e^{w(x)} \left(1 - f_h(x) + \mathcal{O}(h)\right)\tag{71}\] where \(f_h(x) \geq x \cdot \mathbf{1}_{x \in [0,\sqrt{3h}]}\) is defined according to 66 ,69 ,70 , and the bounding constant of the big-\(\mathcal{O}\) term depends on \(s\). In each of the three cases above, the remainder term \(\mathcal{O}(h)\) is uniform in \(x\) on the corresponding interval. Therefore, after taking the maximum of the three remainder constants and the minimum of the corresponding smallness thresholds, there exist constants \(C_1(s,M)>0\) and \(\delta>0\), independent of \(x\) and \(h\), such that for all \(h<\delta\), \[Pe^{w(x)} \le e^{w(x)}\bigl(1-f_h(x)+C_1h\bigr), \qquad x\ge 0.\] Together with the bounds established above, \[f_h(x)\ge x\,\mathbf{1}_{\{x\in[0,\sqrt{3h}]\}},\] which proves 35 .

6.2 Additional proofs from Section 4↩︎

Lemma 8. There exists \(\delta>0\) and constants \(b_1,b_2,\sigma_1,\sigma_2\) such that, for every \(h<\delta\), \[b(x)\in [b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2],\] whenever \[x-\sigma(x)\sqrt{3h}<0 \qquad\text{or}\qquad x-\sigma(x)\sqrt{3h}+b(x)h<0.\]

Proof. We first show that the boundary layer is contained in a fixed compact neighborhood of \(0\) for all sufficiently small \(h\). By Assumption 4, there exists \(L>0\) such that \[|b(x)|+|\sigma(x)|\leq L(1+x),\qquad x\geq 0.\] If \(x-\sigma(x)\sqrt{3h}<0\), then \[x<\sigma(x)\sqrt{3h}\leq |\sigma(x)|\sqrt{3h}\leq L(1+x)\sqrt{3h}.\] If \(x-\sigma(x)\sqrt{3h}+b(x)h<0\), then \[x<\sigma(x)\sqrt{3h}-b(x)h \leq |\sigma(x)|\sqrt{3h}+|b(x)|h \leq L(1+x)(\sqrt{3h}+h).\] Hence, after choosing \(\delta_0>0\) small enough so that \[L(\sqrt{3\delta_0}+\delta_0)\leq \frac{1}{2},\] both cases imply \[x\leq 2L(\sqrt{3h}+h),\qquad h<\delta_0.\] In particular, by decreasing \(\delta_0\) if necessary, every such \(x\) lies in the interval \([0,r_0]\) from Assumption 3. Moreover, \[|\partial_x b(y)|+|\partial_x\sigma(y)|\leq C_{\mathrm{loc}}, \qquad y\in[0,r_0].\] Now we sharpen this localization using the local boundedness of \(\partial_x b\) and \(\partial_x\sigma\) on \([0,r_0]\). Denote \[\sigma_0:=\sigma(0),\qquad b_0:=b(0),\] If \(x-\sigma(x)\sqrt{3h}+b(x)h<0\), then by the mean value theorem, \[\label{eqn:diffusion95boundary95layer95expansion} \begin{align} x &< \sigma(x)\sqrt{3h}-b(x)h\\ &= \left(\sigma_0+x\partial_x\sigma(z_x)\right)\sqrt{3h} -\left(b_0+x\partial_x b(y_x)\right)h \end{align}\tag{72}\] for some \(y_x,z_x\in[0,x]\). Therefore \[x\left(1-\partial_x\sigma(z_x)\sqrt{3h}+\partial_x b(y_x)h\right) < \sigma_0\sqrt{3h}-b_0h.\] After decreasing \(\delta_0\) if necessary, we have \[1-\partial_x\sigma(z_x)\sqrt{3h}+\partial_x b(y_x)h\geq \frac{1}{2}, \qquad h<\delta_0,\] and hence \[x<\sigma_0\sqrt{3h}+C_1h,\] where \(C_1\) depends only on \(b_0,\sigma_0\) and the local bound \(C_{\mathrm{loc}}\).

Similarly, if \(x-\sigma(x)\sqrt{3h}<0\), then \[\label{eqn:diffusion95boundary95layer95expansion952} x < \sigma(x)\sqrt{3h} = \left(\sigma_0+x\partial_x\sigma(z_x)\right)\sqrt{3h},\tag{73}\] so \[x\left(1-\partial_x\sigma(z_x)\sqrt{3h}\right)<\sigma_0\sqrt{3h}.\] After decreasing \(\delta_0\) if necessary, \[1-\partial_x\sigma(z_x)\sqrt{3h}\geq \frac{1}{2},\] and therefore \[x<\sigma_0\sqrt{3h}+C_2h,\] where \(C_2\) depends only on \(\sigma_0\) and \(C_{\mathrm{loc}}\). Let \(C:=C_1\vee C_2\). We have shown that for every \(h<\delta_0\), \[\left\{x\geq 0:\;x-\sigma(x)\sqrt{3h}<0 \;\text{or}\; x-\sigma(x)\sqrt{3h}+b(x)h<0\right\} \subseteq [0,\sigma_0\sqrt{3h}+Ch].\] Now choose \(\delta\in(0,\delta_0]\) small enough so that \[B_\delta:=[0,\sigma_0\sqrt{3\delta}+C\delta]\subseteq[0,r_0].\] Since \(b\) and \(\sigma\) are continuous on \(B_\delta\), define \[b_1:=\min_{x\in B_\delta} b(x),\qquad b_2:=\max_{x\in B_\delta} b(x),\] and \[\sigma_1:=\min_{x\in B_\delta}\sigma(x),\qquad \sigma_2:=\max_{x\in B_\delta}\sigma(x).\] By Assumption 3, after decreasing \(\delta\) if necessary, \(\sigma_1>0\). Therefore, for every \(h<\delta\) and every \(x\) satisfying either boundary-layer inequality, \[b(x)\in[b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2],\] as claimed. ◻

As a side remark, \(\sigma_i\) and \(b_i\) are different from \(\sigma_0 = \sigma(0)\) and \(b_0 = b(0)\) up to a term of order \(\sqrt{h}\) because we have shown the boundary layer of our scheme is included in set \(B_\delta\). Lemma 8 is useful in the sense that it holds when we want to do error analysis for our scheme whenever we do not make a standard EM jump.

Lemma 9. There exists \(\delta>0\) such that \(\forall h < \delta\), \(\forall x\) such that \(x - \sigma(x)\sqrt{3h} < 0,\) \[\lambda(x) \in [0,1]. \;\]

Proof. By Lemma 8, there exist \(\delta_1>0\) and constants \(b_1,b_2,\sigma_1,\sigma_2\), depending only on the local boundary constants and the linear-growth constant of \(b,\sigma\), such that for every \(h<\delta_1\) and every \(x\) with \(x-\sigma(x)\sqrt{3h}<0\), we have \[b(x)\in [b_1,b_2], \qquad \sigma(x)\in [\sigma_1,\sigma_2].\] Hence \[x<\sigma(x)\sqrt{3h}\le \sigma_2\sqrt{3h}.\] Therefore it suffices to prove the following algebraic statement: for all \((x,b,\sigma)\in[0,\sigma_2\sqrt{3h}]\times[b_1,b_2]\times[\sigma_1,\sigma_2]\), the quantity \(\Lambda(x,b,\sigma)\) defined below 75 satisfies \[0\le \Lambda(x,b,\sigma)\le 1.\] If \(b_2 \leq 0\), choose \(\delta_2 = 1\); if \(b_2 > 0\), \(\delta_2\) is chosen in such a way that the following inequalities hold for all \(h < \delta_2\), \[\begin{align} \frac{\kappa}{\sqrt{3h}} &> \frac{2\kappa b_2 \sigma_2}{\sigma_1^2},\\ \sqrt{3h} &< \frac{\sigma_2}{2b_2}, \end{align}\] which implies that for all \(b\in [b_1, b_2]\) and \(\sigma \in [\sigma_1, \sigma_2]\) \[\label{eqn:delta95lambda} \begin{align} \frac{\kappa}{\sqrt{3h}} &> \frac{(2\kappa b - \sigma^2)\sigma_2}{\sigma^2},\\ \frac{\kappa \sigma}{\sqrt{3h}} &> -(\sigma^2 - 2\kappa b).\\ \end{align}\tag{74}\] Consider, \(\forall h < \min\{\delta_1, \delta_2\}\), \[\label{eqn:Lambda} \Lambda(x, b, \sigma) = \frac{(\sigma^2-2\kappa b)x^2 + 2\kappa \sigma^2 x + \sigma^4h}{\left(\frac{\kappa\sigma}{\sqrt{3h}} + \sigma^2 - 2\kappa b\right)x^2 + \sigma^3 \kappa \sqrt{3h} + \sigma^4 h},\tag{75}\] and \[\label{eqn:one95minus95lambda} 1-\Lambda(x, b, \sigma) = \frac{{\frac{\kappa \sigma}{\sqrt{3h}}\left(x-\sigma\sqrt{3h}\right)^2}}{\left(\frac{\kappa\sigma}{\sqrt{3h}} + \sigma^2 - 2\kappa b\right)x^2 + \sigma^3 \kappa \sqrt{3h} + \sigma^4 h}.\tag{76}\] It suffices to conclude the lemma if we can show \(\forall (x,b, \sigma) \in [0,\sigma_2\sqrt{3h}]\times [b_1, b_2] \times [\sigma_1, \sigma_2]\), \[\label{eqn:sign95of95lambda} \Lambda(x, b, \sigma) \geq 0, \quad \text{and} \quad 1-\Lambda(x, b, \sigma) \geq 0.\tag{77}\] We consider three cases: (i) \(b < \frac{\sigma^2}{2\kappa}\) (ii) \(b > \frac{\sigma^2}{2\kappa}\) ; (iii) \(b = \frac{\sigma^2}{2\kappa}\).

For \(b < \frac{\sigma^2}{2\kappa}\), we notice that both numerator (denote as \(N_1(x,b)\)) and denominator (denote as \(D(x,b)\)) of \(\Lambda(x,b)\) are strictly positive. Besides, the numerator of \(1-\Lambda(x,b)\) is a square and its denominator is identical to \(D(x,b)\), so it is strictly positive. 77 is verified.
For \(b > \frac{\sigma^2}{2\kappa}\), we have \(2\kappa b-\sigma^2>0\), and in particular \(b_2>0\). From the first inequality in 74 , for every \(x\in [0,\sigma_2\sqrt{3h}]\), \[(2\kappa b-\sigma^2)x \le (2\kappa b-\sigma^2)\sigma_2\sqrt{3h} < \kappa\sigma^2.\] Hence \[2\kappa\sigma^2-(2\kappa b-\sigma^2)x>\kappa\sigma^2>0.\] Using this, the numerator of \(\Lambda\) satisfies \[\begin{align} N_1(x,b,\sigma) &= (\sigma^2-2\kappa b)x^2 + 2\kappa \sigma^2 x + \sigma^4h \\ &= \sigma^4h + x\bigl(2\kappa \sigma^2-(2\kappa b-\sigma^2)x\bigr) \\ &> \sigma^4h + \kappa\sigma^2 x \\ &> 0. \end{align}\] For the denominator, by the second inequality in 74 , \[\frac{\kappa\sigma}{\sqrt{3h}}+\sigma^2-2\kappa b>0.\] Therefore \[\begin{align} D(x,b,\sigma) &= \left(\frac{\kappa\sigma}{\sqrt{3h}}+\sigma^2-2\kappa b\right)x^2 +\sigma^3\kappa\sqrt{3h} +\sigma^4h \\ &\ge \sigma^3\kappa\sqrt{3h}+\sigma^4h \\ &>0. \end{align}\] Thus \(\Lambda(x,b,\sigma)\ge 0\) on \([0,\sigma_2\sqrt{3h}]\times [b_1,b_2]\times [\sigma_1,\sigma_2]\). We also have \(1-\Lambda(x,b,\sigma)\ge 0\). Hence 77 holds in this case.
For \(b = \frac{\sigma^2}{2\kappa}\), by simple algebra, \[\Lambda(x,b,\sigma) = \frac{2\sigma^2x + 2b\sigma^2h}{\frac{\sigma x^2}{\sqrt{3h}}+ \sigma^3\sqrt{3h} + 2b\sigma^2 h} \ge 0,\]

\[1-\Lambda(x,b,\sigma) = \frac{\frac{\sigma}{\sqrt{3h}}\left(x-\sigma\sqrt{3h}\right)^2}{\frac{\sigma x^2}{\sqrt{3h}}+ \sigma^3\sqrt{3h} + 2b\sigma^2 h} \ge 0.\] Hence 77 is verified in this case.

◻

Lemma 12. There exists constant \(C, \delta > 0\) dependent on \(\sigma(x)\) and \(b(x)\) such that \(\forall h < \delta\) and \(\forall x\in [0, \sigma_2\sqrt{3h}]\) such that \(x-\sigma(x)\sqrt{3h} < 0\), \[0 \leq \lambda(x) \leq C\cdot \left(\frac{x}{\sqrt{h}} + \sqrt{h} \right).\]

Proof. By Lemma 8, after decreasing \(\delta\) if necessary, there exist constants \(b_1,b_2,\sigma_1,\sigma_2\), independent of \(x\) and \(h\), such that whenever \(h<\delta\) and \(x-\sigma(x)\sqrt{3h}<0\), \[b(x)\in[b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2],\] and hence \[0\leq x<\sigma(x)\sqrt{3h}\leq \sigma_2\sqrt{3h}.\] Let \(B_M:=|b_1|\vee |b_2|\), by Lemma 9, after decreasing \(\delta\) if necessary, we already know that \(\lambda(x)\in [0,1]\). It remains to prove the sharper upper bound. Write \[\lambda(x)=\frac{N_h(x)}{D_h(x)},\] where \[N_h(x):=(\sigma^2(x)-2\kappa b(x))x^2+2\kappa\sigma^2(x)x+\sigma^4(x)h,\] and \[D_h(x):= \left(\frac{\kappa\sigma(x)}{\sqrt{3h}}+\sigma^2(x)-2\kappa b(x)\right)x^2 +\sigma^3(x)\kappa\sqrt{3h}+\sigma^4(x)h.\] Decrease \(\delta\) once more so that, for all \(h<\delta\), \[\frac{\kappa\sigma_1}{\sqrt{3h}}\geq 2\kappa B_M.\] Then, for all \(x\) under consideration, \[\frac{\kappa\sigma(x)}{\sqrt{3h}}+\sigma^2(x)-2\kappa b(x) \geq \frac{\kappa\sigma_1}{\sqrt{3h}}-2\kappa B_M \geq 0.\] Therefore \[D_h(x)\geq \sigma^3(x)\kappa\sqrt{3h}\geq \sigma_1^3\kappa\sqrt{3h}.\] On the other hand, \[|N_h(x)| \leq (\sigma_2^2+2\kappa B_M)x^2 +2\kappa\sigma_2^2x +\sigma_2^4h.\] Since \(x\leq\sigma_2\sqrt{3h}\), we have \(x^2\leq 3\sigma_2^2h\). Hence \[\label{eqn:diffusion95lambda95UB} \begin{align} \lambda(x) &\leq \frac{|N_h(x)|}{D_h(x)}\\ &\leq \frac{(\sigma_2^2+2\kappa B_M)\sigma_2^2\cdot 3h +2\kappa\sigma_2^2x+\sigma_2^4h}{\sigma_1^3\kappa\sqrt{3h}}\\ &= \frac{2\sigma_2^2}{\sqrt{3}\sigma_1^3}\cdot \frac{x}{\sqrt h} + \frac{4\sigma_2^4+6\kappa B_M\sigma_2^2}{\sqrt{3}\sigma_1^3\kappa}\cdot \sqrt h. \end{align}\tag{78}\] Choosing \[C:= \max\left\{ \frac{2\sigma_2^2}{\sqrt{3}\sigma_1^3}, \frac{4\sigma_2^4+6\kappa B_M\sigma_2^2}{\sqrt{3}\sigma_1^3\kappa} \right\}\] completes the proof. ◻

Lemma 13. Let \(u\) be the solution of 45 satisfying corresponding boundary condition 44 and assume that \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\). Then, there exists \(\delta\) small such that for all \(h<\delta\) and \(x\in \Omega_{b}\), \[\label{uphadkfv} \begin{align} \partial_x u(t,x) &= \left(x+M(x) \right)\partial_{xx}u(t,x) + Re(x) \end{align}\tag{79}\] where \[|Re(x)|\leq C(h+x)=\mathcal{O}(h+x),\] and \(C\) depends on bounds of \(\partial_xu,\partial_{xx}u,u^{(3)}\) and on the local bounds of \(b,\sigma\) and their derivatives up to order \(2\) on the boundary neighborhood \([0,r_0]\), but is independent of \(x\) and \(h\).

Proof. Similar to proof of Lemma 2, apply Taylor expansion with remainder to 44 centered at \(x\). By Lemma 8, after decreasing \(\delta\) if necessary, every \(x\in\Omega_{b}\) lies in \(B_\delta\subset[0,r_0]\). Hence all intermediate points in the Taylor expansions below also lie in \([0,r_0]\), where the local coefficient bounds from Assumption 3 apply. We have \[\begin{align} \sigma^2(0) \partial_xu(t,0) &= \sigma^2(x)\partial_x u(t,x) - x\sigma^2(x)\partial_{xx}u(t,x) -2x\sigma(x)\sigma'(x)\partial_xu(t,x) \\&\qquad+ \frac{x^2}{2}\partial_{xx}\left(\sigma^2(y_{1,x})\partial_x u(t,y_{1,x})\right)\\ &= \sigma^2(x)\partial_x u(t,x) - x\sigma^2(x)\partial_{xx}u(t,x) + \mathcal{O}(h+x) \end{align}\] and \[\begin{align} &\quad \quad 2\kappa b(0)\partial_x u(t,0) + \kappa \sigma^2(0)\partial_{xx}u(t,0)\\ &= 2\kappa \left(b\cdot\partial_xu - x\partial_x\left(b\cdot \partial_x u\right) + \frac{x^2}{2}\partial_{xx}\left(b\cdot \partial_x u\right)\right)\\ &\quad + \kappa \left(\sigma^2\cdot \partial_{xx}u - x\cdot\partial_x\left(\sigma^2\cdot \partial_{xx}u\right)\right)\\ &= 2\kappa b(x)\partial_x u(t,x) + \left(- 2x\kappa b(x) + \kappa \sigma^2(x)\right)\partial_{xx}u(t,x) + \mathcal{O}(h+x) \end{align}\] where we omit function arguments in expansion above. We choose \(\delta\) small such that \(\sigma^2(x) - 2\kappa b(x) \neq 0\) for all \(x\in [0,\sigma_2\sqrt{3h}]\) because we assume drift and diffusion terms to be continuous and \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\). This allows us to conclude that \((\sigma^2(x) - 2\kappa b(x))^{-1}\) to be bounded for \(x\in [0,\sigma_2\sqrt{3h}]\). This allows us to finish proving the lemma by equating two expressions above and rearranging terms. ◻

Lemma 14. There exists \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), we have, \[\label{vhtfjdmc} \begin{align} \mathbb{E}_x(\Delta X_{k}) &= \frac{\lambda(x)x^2}{2\sqrt{3h} \cdot \sigma(x)} + \frac{\sigma(x) \lambda(x)\sqrt{3h}}{2} + \lambda(x)b(x)h - x + \mathcal{O}(hx + h^2).\\ \end{align}\tag{80}\] where \[R_1(x) = \frac{\lambda(x)b(x)^2h^2}{\sigma(x)\sqrt{3h}} = \mathcal{O}(hx + h^2)\]

Proof. We first calculate \(\mathbb{E}_x(X_{k+1})\) where \(b(x)<0\). According to 48 , given \(X_k=x\), \(X_{k+1}\) can take four different forms

If \(x + \sigma(x)Z_{k+1} \geq 0\) and \(x + \sigma(x)Z_{k+1} + b(x)h \geq 0\), then \[X_{k+1} = x + \sigma(x)Z_{k+1} + b(x)h \;\;\; \text{where}\;\;\; Z_{k+1} \geq -\frac{1}{\sigma(x)}\left(b(x)h+x\right)\]
If \(x + \sigma(x)Z_{k+1} \geq 0\) and \(x + \sigma(x)Z_{k+1} + b(x)h \leq 0\), then \[X_{k+1} = -x - \sigma(x)Z_{k+1} - b(x)h \;\;\; \text{where}\;\;\; -\frac{x}{\sigma(x)}\leq Z_{k+1} \leq -\frac{1}{\sigma(x)}\left(b(x)h+x\right)\]
If \(x + \sigma(x)Z_{k+1} \leq 0\) and \(x + \sigma(x)Z_{k+1} + b(x)h \leq 0\), then \[X_{k+1} = x + \sigma(x)Z_{k+1} - b(x)h \;\;\; \text{where}\;\;\; \frac{1}{\sigma(x)}\left(b(x)h-x\right)\leq Z_{k+1} \leq -\frac{x}{\sigma(x)}\]
If \(x + \sigma(x)Z_{k+1} \leq 0\) and \(x + \sigma(x)Z_{k+1} + b(x)h \geq 0\), then \[X_{k+1} = -x - \sigma(x)Z_{k+1} + b(x)h \;\;\; \text{where}\;\;\; Z_{k+1} \leq \frac{1}{\sigma(x)}\left(b(x)h-x\right)\]

Therefore, \[\label{eqn:1st95mom95Xk43195negative95drift95diffusion} \begin{align} \mathbb{E}_x(X_{k+1}) &= \lambda(x)\mathbb{E}_x \left( \left| \left|x + \sigma(x)Z_{k+1}\right| + b(x)h \right| \right) + (1-\lambda(x))\cdot 0\\ &= \int_{-\frac{1}{\sigma(x)}\left(b(x)h+x\right)}^{\sqrt{3h}} (x+\sigma(x)y+b(x)h)\cdot \frac{\lambda(x)}{2\sqrt{3h}}\;dy\\ &\quad + \int_{-\frac{x}{\sigma(x)}}^{-\frac{1}{\sigma(x)}\left(b(x)h+x\right)} (-x-\sigma(x)y-b(x)h)\cdot \frac{\lambda(x)}{2\sqrt{3h}}\;dy\\ &\quad + \int_{\frac{1}{\sigma(x)}\left(b(x)h-x\right)}^{-\frac{x}{\sigma(x)}} (x+\sigma(x)y-b(x)h)\cdot \frac{\lambda(x)}{2\sqrt{3h}}\;dy\\ &\quad + \int_{-\sqrt{3h}}^{\frac{1}{\sigma(x)}\left(b(x)h-x\right)} (-x-\sigma(x)y+b(x)h)\cdot \frac{\lambda(x)}{2\sqrt{3h}}\;dy \\ &= \frac{\lambda(x)x^2}{2\sqrt{3h} \cdot \sigma(x)} + \frac{\sigma(x) \lambda(x)\sqrt{3h}}{2} + \lambda(x)b(x)h + \frac{\lambda(x)b(x)^2h^2}{\sigma(x)\sqrt{3h}}. \end{align}\tag{81}\] This gives that \[\label{eqn:first95mom95negative95drifted95diffusion} \begin{align} \mathbb{E}_x(\Delta X_{k}) &= \frac{\lambda(x)x^2}{2\sqrt{3h} \cdot \sigma(x)} + \frac{\sigma(x) \lambda(x)\sqrt{3h}}{2} + \lambda(x)b(x)h - x + \frac{\lambda(x)b(x)^2h^2}{\sigma(x)\sqrt{3h}}.\\ \end{align}\tag{82}\] When the drift term \(b(x)\) is positive calculation is simpler (as there is only absolute value on \(x+\sigma(x)Z_{k+1}\)), we can verify that: \[\label{eqn:first95mom95positive95drifted95diffusion} \mathbb{E}_x(\Delta X_{k}) = \frac{\lambda(x)x^2}{2\sqrt{3h} \cdot \sigma(x)} + \frac{\sigma(x) \lambda(x)\sqrt{3h}}{2} + \lambda(x)b(x)h - x.\tag{83}\] In other word, the proof is complete if the \(b(x) \geq 0\) and \(R_1 = 0\). Back to the case where \(b(x) < 0\), the first conditional moment of \(\Delta X_k\) is \[\mathbb{E}_x(\Delta X_k) = \eqref{eqn:first95mom95positive95drifted95diffusion} + \frac{\lambda(x)b(x)^2h^2}{\sigma(x)\sqrt{3h}},\] and the remainder \[\label{eqn:first95mom95diffusion95remainder} R_1(x) = \frac{\lambda(x)b(x)^2h^2}{\sigma(x)\sqrt{3h}} = \mathcal{O}(hx + h^2)\tag{84}\] because of Lemma 12 and the bounds \(b(x)\in[b_1,b_2]\), \(\sigma(x)\geq\sigma_1>0\) from Lemma 8. ◻

Lemma 15. There exists \(\delta>0\) and \(C > 0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), we have, \[\mathbb{E}_x(\Delta X_k^2) = -\frac{\lambda(x)}{\sigma(x)\sqrt{3h}}\cdot x^3 + (1+\lambda(x))x^2 - \lambda(x)\sigma(x)\sqrt{3h}\cdot x + \lambda(x)\sigma^2(x)\cdot h + R_2(x)\] where \[|R_2(x)| \leq C\cdot (hx + h^2) = \mathcal{O}(hx + h^2).\]

Proof. Again, by Lemma 8, after decreasing \(\delta\) if necessary, every \(x\in\Omega_{b}\) lies in \(B_\delta\subset[0,r_0]\). When drift is negative, recycle part of calculation in 81 and separate out remainder terms of order greater that \(3/2\), we have: (when drift is positive, \(\frac{b(x)^2h^2}{\sigma(x)\sqrt{3h}}\) will be removed from \(R_2(x)\), but we still reach same conclusion) \[\begin{align} & \mathbb{E}_x(\Delta X_k^2) \\ &= \lambda(x)\mathbb{E}_x\left(\left(\left| \left|x + \sigma(x) Z_{k+1}\right| + b(x)h \right| - x\right)^2\right) + (1-\lambda(x))x^2\\ &= \lambda(x)\mathbb{E}_x\left( \left( \left|x + \sigma(x) Z_{k+1}\right| + b(x)h \right)^2 - 2x\left| \left|x + \sigma(x)Z_{k+1}\right| + b(x)h \right| + x^2 \right) + (1-\lambda(x))x^2\\ &= \lambda(x)\mathbb{E}_x\left( \left|x + \sigma(x) Z_{k+1}\right|^2 - 2x \left| \left|x + \sigma(x) Z_{k+1}\right| + b(x)h \right| + x^2 \right) + (1-\lambda(x))x^2 + R_2'(x)\\ &= \lambda(x) \cdot \left( x^2 + \sigma^2(x) h - 2x \left( \frac{x^2}{2 \sigma(x) \sqrt{3h}} + \frac{\sigma(x)\sqrt{3h}}{2} + b(x)h + \frac{b(x)^2h^2}{\sigma(x)\sqrt{3h}}\right) + x^2 \right)\\ & \quad + (1-\lambda(x))x^2 + R_2'(x)\\ &= \lambda(x) \cdot \left( 2x^2 + \sigma^2(x)h - 2x \left( \frac{x^2}{2\sigma(x)\sqrt{3h}} + \frac{\sigma(x)\sqrt{3h}}{2} \right) \right)\\ &\quad + (1-\lambda(x))x^2 - 2x\lambda(x)\left(b(x)h + \frac{b(x)^2h^2}{\sigma(x)\sqrt{3h}}\right) + R_2'(x)\\ &= -\frac{\lambda(x)}{\sigma(x)\sqrt{3h}}\cdot x^3 + (1+\lambda(x))x^2 - \lambda(x)\sigma(x)\sqrt{3h}\cdot x + \lambda(x)\sigma^2(x)\cdot h + R_2(x) \end{align}\] where \[R_2'(x) = \lambda(x)\mathbb{E}_x\left(b(x)^2h^2 + 2b(x)h\left|x+\sigma(x)Z_{k+1}\right|\right)\] and \[\begin{align} R_2(x) &= -2x\lambda(x)\left(b(x)h + \frac{b(x)^2h^2}{\sigma(x)\sqrt{3h}}\right) + R_2'(x)\\ &= 2x\lambda(x)\left(b(x)h + \frac{b(x)^2h^2}{\sigma(x)\sqrt{3h}}\right) + \lambda(x)\mathbb{E}_x\left(b(x)^2h^2 + 2b(x)h\left|x+\sigma(x)Z_{k+1}\right|\right). \end{align}\] We can bound \(R_2(\cdot)\) using the local boundary-layer bounds from Lemma 8. Recall that \(B_M:=|b_1|\vee |b_2|\) and \(\sigma(x)\geq\sigma_1>0\) on the boundary layer: \[|R_2(x)|\leq 2B_Mhx + \frac{2B_M^2h^2x}{\sigma_1\sqrt{3h}} + B_M^2h^2 + 2B_Mh\mathbb{E}_x\left(\left|x+\sigma(x)Z_{k+1}\right|\right)\lambda(x).\] Hence, it suffices to complete the proof by showing that \(h\mathbb{E}_x\left(\left|x+\sigma(x)Z_{k+1}\right|\right)\lambda(x) = \mathcal{O}(hx+h^2)\). By Lemma 12, and using \(x\leq \sigma_2\sqrt{3h}\) and \(|Z_{k+1}|\leq \sqrt{3h}\), we have \[\begin{align} h\mathbb{E}_x\left(\left|x+\sigma(x)Z_{k+1}\right|\right)|\lambda(x)| &\leq h\left(x+\sigma_2\sqrt{3h}\right)|\lambda(x)|\\ &\leq C h^{3/2}\left(\frac{x}{\sqrt h}+\sqrt h\right)\\ &\leq C(hx+h^2). \end{align}\] ◻

Lemma 16. Assume that \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), there exist \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\) \[\label{ednfpvzr} 2\left(\mathbb{E}_x(\Delta X_k) - hb(x)\right) \cdot (x + M(x)) + \mathbb{E}_x(\Delta X_k^2) = \sigma^2(x)h + R_3(x)\tag{85}\] where \(|R_3(x)| \leq C(hx+h^2)=\mathcal{O}(hx+h^2)\) for some constant \(C\) independent of \(x\) and \(h\). The constant depends only on the boundary-layer bounds \(b_1,b_2,\sigma_1,\sigma_2\), \(\kappa\), and the lower bound of \(|\sigma^2(x)-2\kappa b(x)|\) in a sufficiently small neighborhood of \(0\).

Proof. Let \[\mathbb{E}_x(\Delta X_k)=A_h(x)+R_1(x),\] where \(A_h(x)\) denotes the principal term from 83 . By Lemma 12, \[\label{eq:R195bound95rewrite} |R_1(x)| \le C_1(hx+h^2)\tag{86}\] for some constant \(C_1\) independent of \(h\).

Likewise, by Lemma 15, \[\mathbb{E}_x(\Delta X_k^2)=B_h(x)+R_2(x),\] where \[B_h(x) = -\frac{\lambda(x)}{\sigma(x)\sqrt{3h}}x^3 +(1+\lambda(x))x^2 -\lambda(x)\sigma(x)\sqrt{3h}\,x +\lambda(x)\sigma^2(x)h,\] and \[\label{eq:R295bound95rewrite} |R_2(x)| \le C_2(hx+h^2)\tag{87}\] for some constant \(C_2\) independent of \(h\). We now substitute these expansions into the quantity of interest: \[2\bigl(\mathbb{E}_x(\Delta X_k)-hb(x)\bigr)(x+M(x)) +\mathbb{E}_x(\Delta X_k^2),\] where \[M(x)=\frac{\kappa \sigma^2(x)}{\sigma^2(x)-2\kappa b(x)}.\] This gives \[\label{eq:split95main95remainder} \begin{align} &2\bigl(\mathbb{E}_x(\Delta X_k)-hb(x)\bigr)(x+M(x)) +\mathbb{E}_x(\Delta X_k^2) \\ &\qquad= \underbrace{2\bigl(A_h(x)-hb(x)\bigr)(x+M(x))+B_h(x)}_{=: \mathcal{A}_h(x)} +2R_1(x)(x+M(x))+R_2(x). \end{align}\tag{88}\]

A direct rearrangement of \(\mathcal{A}_h(x)\) yields \[\label{eq:Ah95rearranged} \begin{align} \mathcal{A}_h(x) &= \lambda(x)\left( \left(\frac{M(x)}{\sigma(x)\sqrt{3h}}+1\right)x^2 +M(x)\sigma(x)\sqrt{3h} +\sigma^2(x)h +2b(x)M(x)h \right) \\ &\quad -x^2-2M(x)x-2b(x)M(x)h +2x(\lambda(x)-1)b(x)h. \end{align}\tag{89}\] By the definition 49 of \(\lambda(x)\), the bracketed terms are chosen precisely so that \[\begin{gather} \lambda(x)\left( \left(\frac{M(x)}{\sigma(x)\sqrt{3h}}+1\right)x^2 +M(x)\sigma(x)\sqrt{3h} +\sigma^2(x)h +2b(x)M(x)h \right)\\ -x^2-2M(x)x-2b(x)M(x)h = \sigma^2(x)h. \end{gather}\] Hence \[\label{eq:Ah95identity} \mathcal{A}_h(x)=\sigma^2(x)h+2x(\lambda(x)-1)b(x)h.\tag{90}\]

Substituting 90 into 88 , we obtain \[2\bigl(\mathbb{E}_x(\Delta X_k)-hb(x)\bigr)(x+M(x)) +\mathbb{E}_x(\Delta X_k^2) = \sigma^2(x)h+R_3(x),\] where \[\label{eq:R395definition95rewrite} R_3(x):= 2x(\lambda(x)-1)b(x)h +2R_1(x)(x+M(x)) +R_2(x).\tag{91}\]

It remains to bound \(R_3(x)\). Since \(b(0)\neq \sigma^2(0)/(2\kappa)\), continuity implies that, for sufficiently small \(\delta>0\), the quantity \[\sigma^2(x)-2\kappa b(x)\] is bounded away from \(0\) on \(x\in[0,\sigma_2\sqrt{3h}]\) whenever \(h<\delta\). Therefore \(M(x)\) is uniformly bounded on this region. Using also the boundary-layer bound \(b(x)\in[b_1,b_2]\), the fact that \(|\lambda(x)-1|\le 1\), the remainder bounds 86 –87 , and \(x\in[0,\sigma_2\sqrt{3h}]\), we conclude that \[|R_3(x)| \le C(hx+h^2)\] for some constant \(C\) independent of \(h\). This proves the claim. ◻

Lemma 17. Assume that \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), there exists \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), \[\label{pmojqelv} \left|\mathbb{E}_x(\Delta X_k)\right| \leq C(h + hx+h^2),\tag{92}\] where \(C>0\) is independent of \(x\) and \(h\) and depends only on \(\kappa\), the local boundary bounds for \(b,\sigma\), and the nondegeneracy constant coming from \(b(0)\neq \sigma^2(0)/(2\kappa)\).

Proof. By Lemma 8, after decreasing \(\delta\) if necessary, all \(x\in[0,\sigma_2\sqrt{3h}]\) lie in the fixed boundary neighborhood where the local coefficient bounds apply. Since \[b(0)\neq \frac{\sigma^2(0)}{2\kappa},\] continuity implies that, after decreasing \(\delta\) if necessary, there exists \(c_0>0\) such that \[|\sigma^2(x)-2\kappa b(x)|\ge c_0, \qquad x\in[0,\sigma_2\sqrt{3h}].\] Hence \(M(x)\) is uniformly bounded on this region. Moreover, since \[M(0)=\frac{\kappa\sigma^2(0)}{\sigma^2(0)-2\kappa b(0)}\neq 0,\] continuity also implies that there exists \(m_0>0\) such that \[|M(x)|\ge m_0,\qquad x\in[0,\sigma_2\sqrt{3h}]\] for all sufficiently small \(h\). Decreasing \(\delta\) once more so that \[\sigma_2\sqrt{3h}\le \frac{m_0}{2}, \qquad h<\delta,\] we obtain \[|x+M(x)|\ge |M(x)|-x\ge \frac{m_0}{2}, \qquad x\in[0,\sigma_2\sqrt{3h}].\] We derive from remainder estimation in Lemma 16 and the crude bound \(\mathbb{E}_x(\Delta X_k^2)\le Ch\) that \[\begin{align} |\mathbb{E}_x(\Delta X_k)| &= \left|\frac{\sigma^2(x)h + R_3(x) - \mathbb{E}_x(\Delta X_k^2)}{2\left(x+M(x)\right)} + hb(x)\right|\\ &\leq \frac{\sigma_2^2 h + |R_3(x)| + (\sigma_2 \sqrt{3h} + B_Mh)^2}{m_0} + B_M h\\ &\leq C(h+hx+h^2). \end{align}\] ◻

Lemma 18. There exists \(\delta>0\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{br}\), \[\label{fpsdxhbe} \left|\mathbb{E}_x(\Delta X_k^3)\right| \leq K_4(hx + h^2),\tag{93}\] where \(K_4\) is a fixed constant dependent on drift, diffusion, and \(\kappa\).

Similar to proof of Lemma 6, fix \(x\in[0,\sigma_2\sqrt{3h}]\), an increment \(\Delta X_k\) is supported on \([-x,\sigma_2\sqrt{3h} + B_Mh]\). We decompose \(\Delta X_k\) into its positive and negative parts, so that \(|\Delta X_k|^3=(\Delta X_k^+)^3+(\Delta X_k^-)^3\) and therefore \[\mathbb{E}_x(|\Delta X_k|^3) =\mathbb{E}_x[(\Delta X_k^+)^3]+\mathbb{E}_x[(\Delta X_k^-)^3].\]

Positive part. Whenever \(\Delta X_k>0\) the increment comes from the continuous branch \(Z_k\) in Algorithm [alg:diffusion95cont95scheme]. This means the largest positive increment is \(\sigma(x)\sqrt{3h} + b(x)h\), hence \[\label{eq:3nd95mom95pospart95diffusion} \mathbb{E}_x[(\Delta X_k^+)^3]\le \lambda(x)\,(\sigma_2\sqrt{3h} + B_Mh)^{3}.\tag{94}\] Recall by Lemma 12, for \(x\in[0,\sqrt{3h}]\), we have the following bound on \(\lambda(x)\), \[|\lambda(x)| \leq \frac{2 \sigma^2_2}{\sqrt{3}\sigma_1^3}\cdot \frac{x}{\sqrt{h}} + \frac{(4\sigma_2^4+6\kappa B_M\sigma_2^2)}{\sqrt{3}\sigma_1^3\kappa}\cdot \sqrt{h} \leq S_1\left(\frac{x}{\sqrt{h}} + \sqrt{h} \right).\] Combining with 94 gives \[\begin{align} \mathbb{E}_x[(\Delta X_k^+)^3] &\leq S_1\left(\frac{x}{\sqrt{h}} + \sqrt{h} \right)(\sigma_2\sqrt{3h} + B_Mh)^{3}\\ &\leq S_1\left( \sqrt{3}\sigma_2x + B_M\sqrt{h}x + \sqrt{3}\sigma_2 h + B_M h^{3/2} \right)(\sigma_2\sqrt{3h} + B_Mh)^{2}. \end{align}\] In other word, exist constant \(K_4^+ > 0\) that \[\label{eq:3nd95mom95pos-final95diffusion} \mathbb{E}_x[(\Delta X_k^+)^3]\le K^+_4 (hx + h^2).\tag{95}\]

Negative part. Because \(\Delta X_k\ge -x\), we have \(0\le \Delta X_k^-\le x\), and therefore \[(\Delta X_k^-)^3 \le x^3.\] Since \(x\in [0,\sigma_2\sqrt{3h}]\), \[\label{eq:3nd95mom95neg-final95diffusion} \mathbb{E}_x[(\Delta X_k^-)^3] \le x^3 \le 3\sigma_2^2hx.\tag{96}\]

Combining 95 and 96 gives \[\mathbb{E}_x(|\Delta X_k|^3) \;\le\; K_4^+ h^2 + (K_4^+ + 3\sigma_2^2)hx.\] Since \(|\mathbb{E}_x(\Delta X_k^3)|\le \mathbb{E}_x(|\Delta X_k|^3)\), the same bound holds for \(|\mathbb{E}_x(\Delta X_k^3)|\) with \(K_4=K_4^+ + 3\sigma_2^2\). We have shown the lemma for \(x\in\Omega_{b}\). It remains to consider \(x\in\Omega_{r}\). In this region the scheme performs the reflected Euler–Maruyama update \[X_{k+1}=\left|x+\sigma(x)Z_{k+1}+b(x)h\right|.\] Let \(\Delta X_k^{\mathrm{EM}}:=\sigma(x)Z_{k+1}+b(x)h\) denote the unreflected Euler–Maruyama increment. By symmetry of \(Z_{k+1}\) and the boundary-layer bounds on \(b\) and \(\sigma\), \[\left|\mathbb{E}_x\left[\left(\Delta X_k^{\mathrm{EM}}\right)^3\right]\right| = \left|3b(x)\sigma^2(x)h^2+b^3(x)h^3\right| \leq Ch^2.\] The reflected and unreflected increments differ only on \[A:=\left\{x+\sigma(x)Z_{k+1}+b(x)h<0\right\}.\] Since \(x\in\Omega_{r}\), Lemma 8 gives \(x=\mathcal{O}(\sqrt h)\), and the interval of \(Z_{k+1}\) values producing reflection has length \(\mathcal{O}(h)\). Because \(Z_{k+1}\) is uniform on an interval of length \(2\sqrt{3h}\), we have \[\mathbb{P}_x(A)=\mathcal{O}(\sqrt h).\] On \(A\), both \(|\Delta X_k|\) and \(|\Delta X_k^{\mathrm{EM}}|\) are \(\mathcal{O}(\sqrt h)\), so \[\left|\mathbb{E}_x\left[\Delta X_k^3-\left(\Delta X_k^{\mathrm{EM}}\right)^3\right]\right| \leq C h^{3/2}\mathbb{P}_x(A) =\mathcal{O}(h^2).\] Hence \[|\mathbb{E}_x(\Delta X_k^3)|\leq K_4'h^2\leq K_4'(hx+h^2),\] and the proof is complete. ◻

Theorem 7. Given \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), there exists \(\delta > 0\) and \(C\geq 0\) depending on the relevant bounds of derivatives of \(u\) and on the local boundary constants for \(b,\sigma\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{b}\), \[\label{grsuwacx} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| \leq C(h^2 + hx) = \mathcal{O}(hx + h^2)\tag{97}\]

Proof. We choose \(\delta\) small enough so that all lemmas we have proved in this section hold for \(h < \delta\). In particular, by Lemma 8, for all \(h<\delta\) and all \(x\in\Omega_{b}\), \[b(x)\in[b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2].\] These constants are local boundary-layer constants and are independent of \(x\) and \(h\). Let \(u\) be the solution to the backward equation 45 . We perform a Taylor expansion of \(u\) to 3rd order centered at \((t_k, X_k)\). We write \(u_k=u(t_k,X_k), \partial_xu_k = \partial_xu(t_k,X_k)\) (and similarly for \(\partial_{xx}u_k, \partial_{xxx}u_k\)), and we write \(Y_k\) for a point such that \(Y_k\in[X_k,X_{k+1}]\) or \(Y_k\in[X_{k+1},X_{k}]\), and \(s_k\) will be a point in interval \(s_k\in[kh, (k+1)h]\)). Expanding gives \[\label{eq:taylor495diffusion} \begin{align} u_{k+1} - u_k &= \partial_t u_k \cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot \Delta X_k^2\\ &\quad + \frac{1}{2}\partial_{tt} u(s_k, X_k) \cdot h^2 + \frac{1}{6}\partial_{xxx}u(t_k, Y_{k})\cdot \Delta X_k^3. \end{align}\tag{98}\] We are interested in \(\mathbb{E}_x(u_{k+1}-u_k)\) and know \(x\in \Omega_{b}\). We now bound the expectation of each of the terms in the Taylor-expansion 98 . We have \[\begin{align} \mathbb{E}_x|\partial_{tt} u(s_k, X_k) \cdot h^2| & \leq h^2\|\partial_{tt}u\|_{\infty}\\ \mathbb{E}_x|\partial_{xxx}u(t_k, Y_{k})\cdot \Delta X_k^3| &\leq C(h^2+hx)\|u^{(3)}\|_\infty, \end{align}\] where the second inequality follows from Lemma 18. We also have, using the backward equation 45 and Lemma 13, \[\begin{align} & \mathbb{E}_x\left[\partial_t u_k\cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot \Delta X_k^2 \right] \\ &\quad = \partial_xu(t_k,x)\left(\mathbb{E}_x(\Delta X_k)-b(x)h\right) +\frac{1}{2}\partial_{xx}u(t_k,x) \left(\mathbb{E}_x(\Delta X_k^2)-\sigma^2(x)h\right) \\ &\quad = \left[ \left(x+M(x)\right)\left(\mathbb{E}_x(\Delta X_k)-b(x)h\right) +\frac{1}{2}\left(\mathbb{E}_x(\Delta X_k^2)-\sigma^2(x)h\right) \right]\partial_{xx}u(t_k,x) \\ &\qquad +\left(\mathbb{E}_x(\Delta X_k)-b(x)h\right)Re(x) \\ &\quad = \frac{1}{2}\partial_{xx}u(t_k,x)R_3(x) +\left(\mathbb{E}_x(\Delta X_k)-b(x)h\right)Re(x). \end{align}\] The last step uses Lemma 16. Therefore, applying the triangle inequality to 98 , and then using Lemma 17 and remainder estimates in Lemma 13, 16 shows that there exist constants \(M_1,M_2\), independent of \(x,t,h\), such that \[\begin{align} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| & \leq \frac{1}{2}\|\partial_{xx}u\|_{\infty} \cdot |R_3(x)| + \left(\left|\mathbb{E}_x(\Delta X_k)\right| + B_Mh\right)\cdot |Re(x)| \\ &\qquad+ h^2\|\partial_{tt}u\|_{\infty} + C(h^2+hx)\|u^{(3)}\|_\infty\\ &\leq M_1 h^2 + M_2 hx \end{align}\] where \(B_M:=|b_1|\vee |b_2|.\) ◻

Lemma 19. There exists \(\delta > 0\) and constant \(C\) which depends only on the boundary-layer bounds such that \(\forall h < \delta\) and \(\forall x\in \Omega_{r}\), \[|E_x(\Delta X_k) - b(x)h| \leq C( hx + h^2) = \mathcal{O}(hx + h^2)\]

Proof. Recall that \[\Omega_{r}:=\left\{x\ge \sigma(x)\sqrt{3h}:\;x-\sigma(x)\sqrt{3h}+b(x)h\le 0\right\}.\] For \(x\in \Omega_{r}\) the scheme performs a reflected Euler–Maruyama step, \[X_{k+1}=\left|x+\sigma(x)Z_{k+1}+b(x)h\right|, \qquad Z_{k+1}\sim \mathrm{Unif}[-\sqrt{3h},\sqrt{3h}].\] Hence \[\begin{gather} \mathbb{E}_x(X_{k+1}) =\\ \frac{1}{2\sqrt{3h}} \int_{-\frac{x+b(x)h}{\sigma(x)}}^{\sqrt{3h}} \left(x+\sigma(x)z+b(x)h\right)\,dz + \frac{1}{2\sqrt{3h}} \int_{-\sqrt{3h}}^{-\frac{x+b(x)h}{\sigma(x)}} \left(-x-\sigma(x)z-b(x)h\right)\,dz. \end{gather}\] A direct calculation gives \[\mathbb{E}_x(X_{k+1}) = \frac{\sigma(x)\sqrt{3h}}{2} +\frac{(x+b(x)h)^2}{2\sigma(x)\sqrt{3h}}.\] Therefore \[\begin{align} \mathbb{E}_x(\Delta X_k)-b(x)h &= \frac{\sigma(x)\sqrt{3h}}{2} +\frac{(x+b(x)h)^2}{2\sigma(x)\sqrt{3h}} -x-b(x)h \\ &= \frac{\left(x+b(x)h-\sigma(x)\sqrt{3h}\right)^2}{2\sigma(x)\sqrt{3h}}. \end{align}\] In particular, \[\left|\mathbb{E}_x(\Delta X_k)-b(x)h\right| = \frac{\left(x+b(x)h-\sigma(x)\sqrt{3h}\right)^2}{2\sigma(x)\sqrt{3h}}.\]

Let \(B_M:=|b_1|\vee |b_2|\). By Lemma 8, after decreasing \(\delta\) if needed, we have \[\sigma(x)\ge \sigma_1>0, \qquad |b(x)|\le B_M\] for all \(x\in \Omega_{r}\) and all \(h<\delta\). Moreover, since \(x\in \Omega_{r}\), \[x\ge \sigma(x)\sqrt{3h}, \qquad x+b(x)h\le \sigma(x)\sqrt{3h}.\] Hence \[0\le \sigma(x)\sqrt{3h}-x-b(x)h \le -b(x)h \le B_M h,\] and therefore \[\left|\mathbb{E}_x(\Delta X_k)-b(x)h\right| \le \frac{B_M^2h^2}{2\sigma(x)\sqrt{3h}} \le \frac{B_M^2}{2\sigma_1\sqrt{3}}\,h^{3/2}.\] Finally, \(x\in \Omega_{r}\) also implies \[x\ge \sigma(x)\sqrt{3h}\ge \sigma_1\sqrt{3h},\] so \[h^{3/2}\le \frac{1}{\sigma_1\sqrt{3}}\,hx.\] Combining the last two bounds yields \[\left|\mathbb{E}_x(\Delta X_k)-b(x)h\right| \le \frac{B_M^2}{6\sigma_1^2}\,hx \le C(hx+h^2)\] for some constant \(C\) independent of \(h\) and \(x\). This completes the proof. ◻

Theorem 8. (Local Error of \(\{X_k\}_{k\in\mathbb{N}}\) in algorithm [alg:diffusion95cont95scheme], situation 2)
Given \(b(0) \neq \frac{\sigma^2(0)}{2\kappa}\), there exists \(\delta > 0\) and \(C \geq 0\) depending on the relevant derivative bounds of \(u\) and on the local boundary constants for \(b\) and \(\sigma\) such that \(\forall h < \delta\) and \(\forall x\in \Omega_{r}\), \[\label{yaiwxvse} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| \leq C(h^2 + hx) = \mathcal{O}(hx + h^2)\tag{99}\]

Proof. The proof is a simplified version of proof of Theorem 7. By Lemma 8, after decreasing \(\delta\) if necessary, for all \(h<\delta\) and all \(x\in\Omega_{r}\), \[b(x)\in[b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2].\] Let \(B_M:=|b_1|\vee |b_2|\). Recall from regularity assumption and Lemma 18, \[\begin{align} \mathbb{E}_x|\partial_{tt} u(s_k, X_k) \cdot h^2| & \leq h^2\|\partial_{tt}u\|_{\infty}\\ \mathbb{E}_x|\partial_{xxx}u(t_k, Y_{k})\cdot \Delta X_k^3| &\leq C(h^2+hx)\|u^{(3)}\|_\infty, \end{align}\] We also have, using the fact that \(u\) solves the backward equation 45 , \[\begin{align} & \mathbb{E}_x\left[\partial_t u_k\cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot (\Delta X_k^2) \right] \\ &\quad = \partial_xu(t_k,x)\cdot \left(\mathbb{E}_x(\Delta X_k) - b(x)h\right) + \frac{1}{2}\partial_{xx}u(t_k,x) \cdot \left(-\sigma^2(x)\cdot h + \mathbb{E}_x(\Delta X_k^2)\right). \end{align}\] Apply Lemma 19 and apply the lower bound of \(\mathbb{E}_x(X_{k+1})\) to the second conditional moment calculation, we have, \[|E_x(\Delta X_k) - b(x)h| \leq C_2( hx + h^2)\] \[|E_x(\Delta X_k^2) - \sigma^2(x)h| \leq C_3( hx + h^2)\] Therefore, applying the triangle inequality to the Taylor expansion 98 , and then using the first- and second-moment estimates above and Lemma 18, there exist constants \(M_1,M_2\), independent of \(x,t,h\), such that \[\begin{align} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| & \leq \|\partial_{x}u\|_{\infty} \cdot C_2( hx + h^2) + \frac{1}{2}\|\partial_{xx}u\|_{\infty} \cdot C_3( hx + h^2)\\ &\qquad+ h^2\|\partial_{tt}u\|_{\infty} + C_1(h^2+hx)\|u^{(3)}\|_\infty\\ &\leq M_1 h^2 + M_2 hx. \end{align}\] ◻

Lemma 20. When \(b(0) = \frac{\sigma^2(0)}{2\kappa}\), there exist \(\delta>0\) and and constant \(C\) which depends only on the boundary-layer bounds such that \(\forall h< \delta\) and \(\forall x\in \Omega_{b}\), \[\left| \mathbb{E}_x\left(\Delta X_k\right) - hb(x) \right| \leq C(hx + h^2) = \mathcal{O}(hx + h^2).\]

Proof. As before,we know there exists \(\delta > 0\) such that for all \(h < \delta\), \[b(x) \in [b_1, b_2]\,,\qquad \sigma(x) \in [\sigma_1, \sigma_2]\,,\qquad \lambda(x)\in[0,1], \qquad \forall x \in \Omega_{b}.\] It suffices to complete the proof if we show the inequality holds for all \(x\in [0, \sigma_2\sqrt{3h}]\). Since \(\sigma(0)>0\), the assumption \(b(0)=\sigma^2(0)/(2\kappa)\) implies \(b(0)>0\). By continuity, after further decreasing \(\delta\) if necessary, \(b(x)\ge 0\) for any \(x\in \Omega_{b}\) and \(\Omega_{r} = \emptyset\). We can check that when \(b(x) = \frac{\sigma^2(x)}{2\kappa}\), \[\mathbb{E}_x\left(\Delta X_k\right) = hb(x),\] we shall expect that \(=\) is replaced by \(\approx\) within the boundary layer. We use Taylor expansion to make the argument more accurate, \(\mathbb{E}_x\left(\Delta X_k\right)\) is a good approximation to \(hb(x)\). We know that, denote \(\sigma_0 := \sigma(0)\) and \(b_0 := b(0)\), \[\sigma^2(x) = \left(\sigma_0 + x\cdot\partial_x\sigma(z_x)\right)^2\,,\] \[b(x) = b_0 + x\cdot\partial_x b(y_x)\,.\] Therefore, \[\label{eqn:sigma95square95miunus952kappa95b95remainder} \begin{align} \sigma^2(x) - 2\kappa b(x) &= \sigma_0^2 - 2\kappa b_0 + \hat{R}(x) = \hat{R}(x) \end{align}\tag{100}\] where \[|\hat{R}(x)| = \left|\left(2\sigma_0\cdot\partial_x\sigma(z_x)-2\kappa\cdot\partial_xb(y_x)\right)x + \left(\partial_x\sigma(z_x)\right)^2 \cdot x^2 \right| \leq \hat{K}(x+x^2).\] where \(\hat{K}\) depends only on the local first-derivative bounds of \(b\) and \(\sigma\) on the boundary neighborhood, and is independent of \(x\) and \(h\). We substitute 100 into 49 , \[\begin{align} \lambda(x) &= \frac{\hat{R}(x)x^2 + 2\kappa \sigma^2(x)x + \sigma^2(x)h\cdot\left(2\kappa b(x) + \hat{R}(x)\right)}{\frac{\kappa \sigma(x)}{\sqrt{3h}}\cdot x^2 + \hat{R}(x)\cdot x^2 + \sigma^3(x)\kappa\sqrt{3h} + \sigma^2(x)h\cdot\left(2\kappa b(x) + \hat{R}(x)\right)}\\ \\ &= \frac{\kappa \sigma^2(x)\left(2x + 2b(x)h + R(x)\right)}{\kappa \sigma^2(x)\left(\frac{x^2}{\sigma(x)\sqrt{3h}} + \sigma(x)\sqrt{3h} + 2b(x)h + R(x)\right)}\\ &= \frac{2x + 2b(x)h + R(x)}{\frac{x^2}{\sigma(x)\sqrt{3h}} + \sigma(x)\sqrt{3h} + 2b(x)h + R(x)} \end{align}\] where \[R(x) = \frac{\sigma^2(x)h +x^2}{\kappa \sigma^2(x)}\cdot \hat{R}(x)\,.\] Recall \(x \le \sigma_2\sqrt{3h}\), \[\frac{\sigma^2(x)h +x^2}{\kappa \sigma^2(x)} \leq \frac{\sigma_2^2h + 3\sigma_2^2h}{k\sigma_1^2} = Ch.\] Hence, \[|R(x)| \leq Ch |\hat{R}(x)| \leq Ch(x + x^2) \leq C(hx+h^2).\] Since the drift is positive, substitute \(\lambda(x)\) back into 83 , we have \[\begin{align} \mathbb{E}_x(\Delta X_k) &= \lambda(x) \cdot \left( \frac{x^2}{2\sigma(x)\sqrt{3h}} + \frac{1}{2}\sigma(x)\sqrt{3h} + b(x)h + \frac{1}{2}R(x) - \frac{1}{2}R(x) \right) - x\\ &= x + b(x)h + \frac{1}{2}\left(1-\lambda(x)\right)R(x) - x\\ &= b(x)h + \frac{1}{2}\left(1-\lambda(x)\right)R(x) \end{align}\] Hence, \[\left| \mathbb{E}_x(\Delta X_k) - b(x)h\right| = \left|\frac{1}{2}\left(1-\lambda(x)\right)R(x)\right| \leq C(hx + h^2)\] because we know \(\lambda(x) \in [0,1]\), \(x\in [0,\sigma_2\sqrt{3h}]\), \(|\hat{R}(x)| \leq \hat{K}(x + x^2)\), and the coefficient \(\frac{\sigma^2(x)h +x^2}{\kappa \sigma^2(x)}\) we scale \(\hat{R}(\cdot)\) to \(R(\cdot)\) is of order \(h\) using the boundary-layer bounds \(\sigma_1\le \sigma(x)\le \sigma_2\) from Lemma 8. ◻

Theorem 9. (Local Error of \(\{X_k\}_{k\in\mathbb{N}}\) in algorithm [alg:diffusion95cont95scheme], situation 3)
If \(b(0)= \frac{\sigma^2(0)}{2\kappa}\), there exist \(\delta > 0\) and \(C\geq 0\), independent of \(x\) and \(h\), depending only on the relevant bounds of derivatives of \(u\) and on the local boundary constants for \(b,\sigma\) such that \(\forall h<\delta\) and for all \(x \in \Omega_{br}\), \[\label{sikgmvnj} \left|\mathbb{E}\left(u_{k+1}-u_k|X_k=x\right)\right| \leq C(h^2 + hx) = \mathcal{O}(hx + h^2).\tag{101}\]

Proof. The proof follows the same Taylor-expansion argument as in Theorem 7. Fix \(x\) such that \(x-\sigma(x)\sqrt{3h}<0\). By Lemma 8, after decreasing \(\delta\) if necessary, for all \(h<\delta\) and all \(x\) satisfying \(x-\sigma(x)\sqrt{3h}<0\), \[b(x)\in[b_1,b_2],\qquad \sigma(x)\in[\sigma_1,\sigma_2].\] These constants are independent of \(x\) and \(h\). By Taylor expansion of \(u\) centered at \((t_k,x)\), with \(Y_k\) between \(x\) and \(X_{k+1}\) and \(s_k\in[t_k,t_{k+1}]\), we have \[\begin{align} u_{k+1}-u_k &= \partial_tu_k\,h + \partial_xu_k\,\Delta X_k + \frac{1}{2} \partial_{xx}u_k\,\Delta X_k^2 \\ &\quad + \frac{1}{2} \partial_{tt}u(s_k,x)\,h^2 + \frac{1}{6} \partial_{xxx}u(t_k,Y_k)\,\Delta X_k^3 . \end{align}\] Taking conditional expectation given \(X_k=x\) and using the backward equation 45 , we obtain \[\label{eq:taylor495diffusion95expectation95s3} \begin{align} &\left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k=x\right)\right| \\ &\le \left|\partial_xu(t_k,x)\right|\cdot \left|\mathbb{E}_x(\Delta X_k)-b(x)h\right| + \frac{1}{2} \left|\partial_{xx}u(t_k,x)\right|\cdot \left|\mathbb{E}_x(\Delta X_k^2)-\sigma^2(x)h\right| \\ &\quad + \frac{h^2}{2}\|\partial_{tt}u\|_\infty + \frac{1}{6} \|u^{(3)}\|_\infty\,\left|\mathbb{E}_x(\Delta X_k^3)\right|. \end{align}\tag{102}\]

Because \(b(0)=\frac{\sigma^2(0)}{2\kappa}\), the sticky boundary condition implies \[\partial_{xx}u(t,0)=0.\] Hence, by Taylor expansion in space, \[\partial_{xx}u(t_k,x)=x\,u^{(3)}(t_k,y_x)\,, \qquad \text{for some y_x\in[0,x].}\] Therefore, \[|\partial_{xx}u(t_k,x)|\le \|u^{(3)}\|_\infty\,x.\] Next, by Lemma 20, \[\left|\mathbb{E}_x(\Delta X_k)-b(x)h\right|\le C_1(hx+h^2).\] Moreover, since \(x-\sigma(x)\sqrt{3h}<0\), the increment is supported in an interval of size \(O(\sqrt h)\), so \[\mathbb{E}_x(\Delta X_k^2)\le C_2 h\] for some constant \(C_2>0\) independent of \(x\) and \(h\) because on the boundary layer the update either jumps to \(0\) or takes a reflected Euler–Maruyama step. Finally, by the third-moment bound already established for the boundary layer, \[\left|\mathbb{E}_x(\Delta X_k^3)\right|\le C_3(hx+h^2).\] Substituting these estimates into the inequality 102 yields \[\begin{align} \left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k=x\right)\right| &\le \|\partial_xu\|_\infty \cdot C_1(hx+h^2) + \frac{1}{2} \|u^{(3)}\|_\infty x \cdot \bigl(\sigma^2(x)h+C_2h\bigr) \\ &\quad + \frac{h^2}{2}\|\partial_{tt}u\|_\infty + \frac{1}{6} \|u^{(3)}\|_\infty C_3(hx+h^2) \\ &\le M_1 h^2 + M_2 hx, \end{align}\] for suitable constants \(M_1,M_2\ge 0\) independent of \(x\) and \(h\), after possibly decreasing \(\delta\) so that all preceding lemmas hold simultaneously. ◻

Proof of Theorem 5. The estimate on \(\Omega_{br}\) follows from the three boundary-layer results already established. More precisely, when \(b(0)\neq \frac{\sigma^2(0)}{2\kappa}\), Theorems 7 and 8 cover \(\Omega_{b}\) and \(\Omega_{r}\), respectively. When \(b(0)=\frac{\sigma^2(0)}{2\kappa}\), we have \(b(0)>0\), so by continuity, after possibly decreasing \(\delta\), one has \(b(x)>0\) throughout the boundary layer; hence \(\Omega_{r}=\varnothing\), and Theorem 9 yields the estimate on all of \(\Omega_{br}=\Omega_{b}\). Therefore, after decreasing \(\delta>0\) finitely many times and enlarging the constant if needed, we may assume that \[\left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k=x\right)\right| \le M_b(h^2+hx), \qquad \forall x\in \Omega_{br},\;\forall h<\delta.\]

We now check \(x \in \Omega_e\). Recall that, outside the boundary layer, the simulation scheme makes a standard Euler–Maruyama jump. Since \[x-\sigma(x)\sqrt{3h}+b(x)h>0,\] the absolute value in 47 is inactive, and therefore \[\Delta X_k=\sigma(x)Z_{k+1}+b(x)h.\] Hence \[\mathbb{E}_x(\Delta X_k)=b(x)h,\qquad \mathbb{E}_x(\Delta X_k^2)=\sigma^2(x)h+b^2(x)h^2,\] \[\mathbb{E}_x(\Delta X_k^3)=3b(x)\sigma^2(x)h^2+b^3(x)h^3,\] and \[\label{eq:EM95fourth95moment95diffusion} \begin{align} \mathbb{E}_x(\Delta X_k^4) &= \mathbb{E}_x\left(\sigma(x)Z_{k+1}+b(x)h\right)^4\\ &= \frac{9}{5}\sigma^4(x)h^2 +6b^2(x)\sigma^2(x)h^3 +b^4(x)h^4. \end{align}\tag{103}\] The Taylor expansion upto fourth order is \[\label{eq:taylor495diffusion95ext} \begin{align} u_{k+1} - u_k &= \partial_t u_k \cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot \Delta X_k^2\\ &\quad + \frac{1}{2}\partial_{tt} u(s_k, X_k) \cdot h^2 + \frac{1}{6}\partial_{xxx}u_k\cdot \Delta X_k^3 + \frac{1}{24}\partial_{xxxx}u(t_k, Y_{k})\cdot \Delta X_k^4. \end{align}\tag{104}\] Using the fact that \(u\) solves the backward equation 45 , we know \[\begin{align} & \mathbb{E}_x\left[\partial_t u_k\cdot h + \partial_xu_k\cdot \Delta X_k + \frac{1}{2}\partial_{xx}u_k\cdot (\Delta X_k^2) \right] \\ &\quad = -b(x)\partial_xu(t_k,x)\cdot h -\frac{\sigma^2(x)}{2}\partial_{xx}u(t_k,x)\cdot h + \partial_xu_k \cdot b(x)h + \frac{1}{2}\partial_{xx}u(t_k,x)\cdot \bigl(\sigma^2(x)h+b^2(x)h^2\bigr)\\ &\quad = \frac{1}{2}b^2(x)\partial_{xx}u(t_k,x)\,h^2. \end{align}\] Taking conditional expectation in 104 , and using the previous calculation, gives \[\begin{align} \left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k=x\right)\right| &\le \frac{1}{2}|b(x)|^2\|\partial_{xx}u\|_\infty h^2 + \frac{h^2}{2}\|\partial_{tt}u\|_{\infty}\\ &\quad + \frac{1}{6}\|u^{(3)}\|_\infty \left|3b(x)\sigma^2(x)h^2+b^3(x)h^3\right|\\ &\quad +\frac{1}{24}\|u^{(4)}\|_\infty\,\left|\frac{9}{5}\sigma^4(x)h^2 +6b^2(x)\sigma^2(x)h^3 +b^4(x)h^4\right|. \end{align}\] By Assumption 4, there exists \(L>0\) such that \[|b(x)|+|\sigma(x)|\le L(1+x),\qquad x\ge 0.\] Hence, using \(h<1\) and \((1+x)^4\le C(1+x^4)\) for \(x\ge0\), we have \[|b(x)|^2+\left|b(x)\right|\sigma^2(x)+|b(x)|^3h \le C(1+x^3)\le C(1+x^4),\] and \[\sigma^4(x)+b^2(x)\sigma^2(x)h+b^4(x)h^2 \le C(1+x^4).\] Substituting these bounds and 103 into the previous estimate yields \[\left|\mathbb{E}\left(u_{k+1}-u_k\mid X_k=x\right)\right| \le M_eh^2(1+x^4),\qquad \forall x\in\Omega_e,\;\forall h<\delta,\] for some constant \(M_e\) independent of \(x\) and \(h\). ◻

Lemma 21. Denote \(\tilde{\lambda}(\cdot)\) as the jump probability we propose in algorithm [alg:cont95scheme] and \(\lambda(\cdot)\) as the one we propose in algorithm [alg:diffusion95cont95scheme], if for the diffusion \(\sigma(0)= 1\) for a sticky diffusion, then \[\left|\lambda(x) - \tilde{\lambda}(x)\right| \leq K_{b,\sigma} \sqrt{h},\quad \forall x\in [0,\sqrt{3h}]\] where \(K_{b,\sigma} > 0\) is a constant depending on the drift term \(b(\cdot)\), diffusion \(\sigma(\cdot)\), and sticky parameter \(\kappa\).

Proof. Recall that \[\lambda(x)=\frac{N_\sigma(x)}{D_\sigma(x)},\] where \[N_\sigma(x):=(\sigma^2(x)-2\kappa b(x))x^2+2\kappa\sigma^2(x)x+\sigma^4(x)h\] and \[D_\sigma(x):= \left(\frac{\kappa\sigma(x)}{\sqrt{3h}}+\sigma^2(x)-2\kappa b(x)\right)x^2 +\sigma^3(x)\kappa\sqrt{3h} +\sigma^4(x)h.\] Similarly, \[\tilde{\lambda}(x)=\frac{N_0(x)}{D_0(x)},\] where \[N_0(x):=x^2+2\kappa x+h, \qquad D_0(x):=\left(\frac{\kappa}{\sqrt{3h}}+1\right)x^2+h+\kappa\sqrt{3h}.\]

Using the same compactness argument as in Lemma 8, and using \(\sigma(0)=1\), there exist \(\delta_1>0\) and constants \(b_1,b_2,\sigma_1,\sigma_2\), with \(\sigma_1>0\), such that for every \(h<\delta_1\) and every \(x\in[0,\sqrt{3h}]\), \[b(x)\in[b_1,b_2], \qquad \sigma(x)\in[\sigma_1,\sigma_2].\] Set \[B_M:=|b_1|\vee |b_2|.\] Moreover, since \(\sigma(0)=1\) and \(\sigma\) is continuously differentiable near \(0\), after decreasing \(\delta_1\) if necessary, there exists a constant \(C_\sigma>0\), independent of \(x\) and \(h\), such that for \(n=1,2,3,4\), \[\label{eq:lambda95approx95sigma95power} |\sigma^n(x)-1|\le C_\sigma x\le C_\sigma\sqrt{3h}, \qquad \forall x\in[0,\sqrt{3h}],\quad h<\delta_1.\tag{105}\]

We also decrease \(\delta_1\) so that \(\delta_1\le 1\) and, for every \(h<\delta_1\), \[\frac{\kappa\sigma_1}{\sqrt{3h}}>\sigma_2^2+2\kappa B_M.\] Then, for every \(x\in[0,\sqrt{3h}]\), \[\label{eq:lambda95approx95D95lower} D_\sigma(x) \ge \sigma^3(x)\kappa\sqrt{3h} \ge \kappa\sigma_1^3\sqrt{3h}.\tag{106}\]

Now define \[\hat{\lambda}(x):=\frac{N_0(x)}{D_\sigma(x)}.\] We first compare \(\lambda(x)\) and \(\hat{\lambda}(x)\). Since \[N_\sigma(x)-N_0(x) = (\sigma^2(x)-1)x^2 -2\kappa b(x)x^2 +2\kappa(\sigma^2(x)-1)x +(\sigma^4(x)-1)h,\] using 105 and \(x\le \sqrt{3h}\) gives \[\label{eq:lambda95approx95N95bound} |N_\sigma(x)-N_0(x)|\le C_Nh\tag{107}\] for some constant \(C_N>0\) independent of \(x\) and \(h\). Hence, by 106 , \[|\lambda(x)-\hat{\lambda}(x)| = \frac{|N_\sigma(x)-N_0(x)|}{D_\sigma(x)} \le \frac{C_Nh}{\kappa\sigma_1^3\sqrt{3h}} \le C'_N\sqrt{h},\] where \(C'_N>0\) is independent of \(x\) and \(h\).

Next, we compare \(\hat{\lambda}(x)\) and \(\tilde{\lambda}(x)\). We have \[D_\sigma(x)-D_0(x) = \left(\frac{\kappa(\sigma(x)-1)}{\sqrt{3h}}+\sigma^2(x)-1-2\kappa b(x)\right)x^2 +\kappa(\sigma^3(x)-1)\sqrt{3h} +(\sigma^4(x)-1)h.\] Using 105 again, together with \(x\le \sqrt{3h}\) and \(|b(x)|\le B_M\), we obtain \[\label{eq:lambda95approx95D95bound} |D_\sigma(x)-D_0(x)|\le C_Dh\tag{108}\] for some constant \(C_D>0\) independent of \(x\) and \(h\). Moreover, \[\hat{\lambda}(x)-\tilde{\lambda}(x) = \frac{D_0(x)-D_\sigma(x)}{D_\sigma(x)}\cdot\frac{N_0(x)}{D_0(x)}.\] Using Lemma 1, which gives \(0\le \tilde{\lambda}(x)\le 1\) on \([0,\sqrt{3h}]\), together with 106 and 108 , we obtain \[|\hat{\lambda}(x)-\tilde{\lambda}(x)| \le \frac{|D_\sigma(x)-D_0(x)|}{D_\sigma(x)}|\tilde{\lambda}(x)| \le \frac{C_Dh}{\kappa\sigma_1^3\sqrt{3h}} \le C'_D\sqrt{h},\] where \(C'_D>0\) is independent of \(x\) and \(h\).

Combining the two estimates, for every \(h<\delta_1\) and every \(x\in[0,\sqrt{3h}]\), \[|\lambda(x)-\tilde{\lambda}(x)| \le |\lambda(x)-\hat{\lambda}(x)|+|\hat{\lambda}(x)-\tilde{\lambda}(x)| \le K_{b,\sigma}\sqrt{h}.\] Here \(K_{b,\sigma}:=C'_N+C'_D\) depends only on \(\kappa\), the local bounds \(b_1,b_2,\sigma_1,\sigma_2\), and the local derivative bound for \(\sigma\), but is independent of \(x\) and \(h\). This completes the proof with \(\delta=\delta_1\). ◻

Lemma 10. There exists \(\delta > 0\) and constant \(C\) such that for all \(h < \delta\), \[\mathbb{E}\left(\sum_{k=0}^{N-1} X_k \cdot \mathbf{1}_{\Omega_{br}} (X_k) \right) < C\]

Proof. Throughout this proof we use the normalization \(\sigma(0)=1\) made before the statement of the lemma. By the proof of Lemma 8, there exist \(\delta_B>0\) and a constant \(K_B>0\), depending only on local bounds for \(b\) and \(\sigma\) near \(0\), such that for every \(h<\delta_B\), \[\label{eq:keydiffusion95Bh} \Omega_{br}\subset B_h:= [0,\sqrt{3h}+K_Bh].\tag{109}\] Hence it is enough to prove a uniform bound with \(\mathbf{1}_{B_h}\) in place of \(\mathbf{1}_{\Omega_{br}}\). The rest of the proof is a discrete supersolution argument. We will construct a bounded barrier function \(\varphi_h=e^{w_h}\) whose one-step drift is negative by an amount comparable to \(x\) on the enlarged boundary layer \(B_h\). More precisely, the key estimate is \[P\varphi_h(x)\leq \varphi_h(x)(1-f_h(x)+C_1h), \qquad f_h(x)\geq x\mathbf{1}_{B_h}(x).\] Once this is proved, the time factor \(e^{K_T(T-t_i)}\) will absorb the harmless \(C_1h\) error, giving a supersolution for the discrete occupation equation. A backward comparison then bounds the expected occupation sum.

We now modify the test function used in Lemma 4 so that it is unchanged near the boundary but becomes constant away from the boundary. This avoids requiring \(b\) and \(\sigma\) to be globally bounded. Fix \(s>16\) and \(M>0\). Choose a fixed number \(R>1\), independent of \(h\). Later we decrease \(\delta_B\) if necessary so that \[B_h\subset [0,R/4].\]

Let \(\chi\in C^2(\mathbb{R})\) satisfy \[0\leq \chi\leq 1,\qquad \chi(r)=1\text{ for }r\leq 0,\qquad \chi(r)=0\text{ for }r\geq 1.\] Define \[\eta_R(x):=\int_0^x \chi\left(\frac{y-R}{R}\right)\,dy.\] Then \(\eta_R(x)=x\) for \(0\leq x\leq R\), while \(\eta_R(x)=A_R\) for \(x\geq 2R\), where \[A_R:=\int_0^{2R}\chi\left(\frac{y-R}{R}\right)\,dy.\] Moreover, \(0\leq \eta_R'(x)\leq 1\) and \(|\eta_R''(x)|\leq C_R\) for a constant \(C_R\) depending only on \(R\) and the cutoff \(\chi\).

Define \[\label{eq:keydiffusion95wh} w_h(x):= \begin{cases} M, & 0\leq x\leq \sqrt{3h}/2,\\ M-s\left(\eta_R(x)-\frac{\sqrt{3h}}{2}\right), & x>\sqrt{3h}/2. \end{cases}\tag{110}\] Then \(w_h\) agrees with the function used in Lemma 4 on \([0,R]\): \[w_h(x)= \begin{cases} M, & 0\leq x\leq \sqrt{3h}/2,\\ M-s\left(x-\frac{\sqrt{3h}}{2}\right), & \sqrt{3h}/2<x\leq R. \end{cases}\] Also, \(w_h\) is constant on \([2R,\infty)\). Let \(\varphi_h(x):=e^{w_h(x)}.\) Since \(w_h(x)\leq M\) for all \(x\geq 0\), \(\varphi_h(x)\leq e^M\) uniformly in \(x\) and \(h\). We first prove that there exist constants \(C_1>0\) and \(\delta_1>0\), independent of \(x\) and \(h\), such that for all \(h<\delta_1\), \[\label{eq:keydiffusion95Pphi95goal} P\varphi_h(x)\leq \varphi_h(x)\left(1-f_h(x)+C_1h\right), \qquad \forall x\geq 0,\tag{111}\] where \(f_h\) satisfies \[\label{eq:keydiffusion95fh95goal} f_h(x)\geq x\mathbf{1}_{B_h}(x).\tag{112}\] We introduce an auxiliary transition operator \(P_0\). The role of \(P_0\) is to separate the genuinely new diffusion effects from the SBM calculation. On \(B_h\), the normalization \(\sigma(0)=1\) implies \(\sigma(x)-1=O(\sqrt h)\) and the scheme makes uniformly distributed jump of size \(O(\sqrt h)\), so the continuous update in the diffusion scheme differs from the SBM reflected update by only \(O(h)\). The remaining complication is that \(P\) and \(P_0\) correspond to different boundary layer . Given \(Z\sim \mathcal{U}[-\sqrt{3h},\sqrt{3h}]\), \[P_0\varphi_h(x)= \begin{cases} \lambda(x)\mathbb{E}\left[\varphi_h(|x+Z|)\right]+(1-\lambda(x))\varphi_h(0), & 0\leq x\leq \sqrt{3h},\\ \mathbb{E}\left[\varphi_h(|x+Z|)\right], & x>\sqrt{3h}. \end{cases}\] Thus \(P_0\) is the same as the SBM transition operator, except that it uses the diffusion jump probability \(\lambda(x)\) instead of \(\tilde{\lambda}(x)\). We claim that for \(x\in B_h\), \[\label{eq:keydiffusion95P95P0} P\varphi_h(x)\leq P_0\varphi_h(x)+C_2h\varphi_h(x),\tag{113}\] where \(C_2\) is independent of \(x\) and \(h\). We prove this by coupling the true diffusion update and the auxiliary update using the same random variable \(Z\sim \mathcal{U}[-\sqrt{3h},\sqrt{3h}]\).

First, after decreasing \(\delta_1\) if necessary, \(B_h\subset[0,R/4]\). Since \(\sigma(0)=1\) and \(\sigma\) is locally Lipschitz near \(0\), there exists a constant \(C_\sigma>0\), independent of \(x\) and \(h\), such that \[|\sigma(x)-1|\leq C_\sigma x\leq C_\sigma(\sqrt{3h}+K_Bh)\leq C\sqrt h, \qquad x\in B_h.\] Similarly, \(b\) is locally bounded on a fixed compact neighborhood of \(0\), so \[|b(x)|\leq C,\qquad x\in B_h.\] Thus, for all possible values of \(Z\) and all \(x\in B_h\), \[|x+Z|\leq x+\sqrt{3h}\leq C\sqrt h,\] and \[\left||x+\sigma(x)Z|+b(x)h\right| \leq x+|\sigma(x)Z|+|b(x)|h \leq C\sqrt h.\] After decreasing \(\delta_1\) once more so that \(C\sqrt h<R\), all one-step positions appearing in the comparison below lie in \([0,R]\). We also record a bound on the Lipschitz constant of \(\varphi_h\). Since \(0\leq \eta_R'\leq 1\), the function \(\eta_R\) is \(1\)-Lipschitz. For \(h\) small enough that \(\sqrt{3h}/2<R\), we can write \[w_h(x)=M-s\left(\eta_R(x)-\frac{\sqrt{3h}}{2}\right)_+,\] Hence \(w_h\) is \(s\)-Lipschitz, uniformly in \(h\). Since \(w_h\leq M\), \[\label{eq:phi95lipschitz} |\varphi_h(y)-\varphi_h(z)| = |e^{w_h(y)}-e^{w_h(z)}| \leq e^M |w_h(y)-w_h(z)| \leq se^M|y-z|.\tag{114}\] Therefore \(\varphi_h\) is Lipschitz with constant \(L_\varphi:=se^M\), independent of \(x\) and \(h\). Finally, on \(B_h\) we also have a uniform lower bound on \(\varphi_h(x)\). Indeed, after decreasing \(\delta_1\) if necessary, \[\label{eq:lower95bound95phi} \varphi_h(x)=\exp(w_h(x)) \geq \exp\left(M-s\left(\frac{\sqrt{3h}}{2}+K_Bh\right)\right) \geq e^{M/2}, \qquad x\in B_h.\tag{115}\] Consequently, any estimate of the form \(Ch\) on \(B_h\) may be rewritten as \(Ch\varphi_h(x)\) after increasing the constant \(C\). We now compare the continuous branches. Let \[Y_0:=|x+Z|\] denote the continuous update in the auxiliary operator \(P_0\). For the true diffusion scheme, define the continuous update \[Y_1:= \left||x+\sigma(x)Z|+b(x)h\right|\] Using \(||a|-|b||\leq |a-b|\) gives \[\label{eq:keydiffusion95cont95branch95inlayer} \begin{align} |Y_1-Y_0| &= \left|\left||x+\sigma(x)Z|+b(x)h\right|-|x+Z|\right|\\ &\leq \left||x+\sigma(x)Z|+b(x)h-|x+Z|\right|\\ &\leq \left||x+\sigma(x)Z|-|x+Z|\right|+|b(x)|h\\ &\leq |(\sigma(x)-1)Z|+|b(x)|h\\ &\leq Ch. \end{align}\tag{116}\] Therefore, \[\label{eq:continuous95branch95comparison} |\varphi_h(Y_1)-\varphi_h(Y_0)| \leq L_\varphi |Y_1-Y_0| \leq Ch \leq Ch\varphi_h(x).\tag{117}\]

It remains to control the possible mismatch between the two boundary thresholds \(\sigma(x)\sqrt{3h}\) and \(\sqrt{3h}\). A mismatch occurs only when one of the inequalities \[x\leq \sigma(x)\sqrt{3h},\qquad x\leq \sqrt{3h}\] holds and the other fails. Equivalently, \(x\) lies between \(\sqrt{3h}\) and \(\sigma(x)\sqrt{3h}\). Hence \[\label{eq:keydiffusion95threshold95mismatch} |x-\sigma(x)\sqrt{3h}| \leq |\sigma(x)\sqrt{3h}-\sqrt{3h}| = |\sigma(x)-1|\sqrt{3h} \leq Ch.\tag{118}\] In this mismatch region, we also need that \(1-\lambda(x)\) is of order \(h\). Using the identity \[\label{eq:keydiffusion95one95minus95lambda} 1-\lambda(x) = \frac{\frac{\kappa\sigma(x)}{\sqrt{3h}}\left(x-\sigma(x)\sqrt{3h}\right)^2}{\left(\frac{\kappa\sigma(x)}{\sqrt{3h}}+\sigma^2(x)-2\kappa b(x)\right)x^2+\sigma^3(x)\kappa\sqrt{3h}+\sigma^4(x)h},\tag{119}\] and using the same denominator lower-bound argument as in Lemma 21, after decreasing \(\delta_1\) if necessary, there exists \(c>0\) such that \[\left(\frac{\kappa\sigma(x)}{\sqrt{3h}}+\sigma^2(x)-2\kappa b(x)\right)x^2+\sigma^3(x)\kappa\sqrt{3h}+\sigma^4(x)h \geq c\sqrt h, \qquad x\in B_h.\] Combining this lower bound with 118 and 119 , we obtain \[\label{eq:lambda95complement95estimate} 1-\lambda(x) \leq \frac{C\frac{1}{\sqrt h}\cdot h^2}{c\sqrt h} \leq Ch.\tag{120}\] We now consider the possible threshold configurations. The following four cases only account for which operator regards \(x\) as being in the sticky boundary layer. If both operators make the same decision, the comparison follows from the \(O(h)\) continuous-branch estimate. If their decisions differ, then \(x\) lies in the thin mismatch region, where \(1-\lambda(x)=O(h)\), so the extra sticky-jump contribution is also only \(O(h)\).

Suppose \(x\in\Omega_{b}\) and \(x\leq \sqrt{3h}.\) Then \[P\varphi_h(x) = \lambda(x)\mathbb{E}\left[\varphi_h(Y_1)\right]+(1-\lambda(x))\varphi_h(0),\] and \[\label{eq:P095boundary95layer} P_0\varphi_h(x) = \lambda(x)\mathbb{E}\left[\varphi_h(Y_0)\right]+(1-\lambda(x))\varphi_h(0).\tag{121}\] Thus, by 117 \[P\varphi_h(x)-P_0\varphi_h(x) = \lambda(x)\mathbb{E}\left[\varphi_h(Y_1)-\varphi_h(Y_0)\right] \leq Ch\varphi_h(x).\]
Suppose that \(x\in\Omega_{b}\) and \(x>\sqrt{3h}\), then \(P_0\) uses only the continuous branch, \[\begin{align} P\varphi_h(x)-P_0\varphi_h(x) &= \lambda(x)\mathbb{E}\left[\varphi_h(Y_1)-\varphi_h(Y_0)\right] +(1-\lambda(x))\left(\varphi_h(0)-\mathbb{E}\left[\varphi_h(Y_0)\right]\right)\\ &\leq Ch\varphi_h(x) +(1-\lambda(x))\cdot C\varphi_h(x)\\ &\leq Ch\varphi_h(x). \end{align}\] where we use 117 , 120 and the fact that \[|\varphi_h(0)-\varphi_h(Y_0)|\leq 2e^M \leq 2e^{M/2}\varphi_h(x)\leq C\varphi_h(x), \qquad x\in B_h,\] which follows from the upper bound \(\varphi_h\leq e^M\) and the lower bound 115 .
Suppose that \(x\notin\Omega_{b}\) and \(x\leq\sqrt{3h}\), then the true diffusion scheme uses only the continuous branch, while \(P_0\) is in its sticky layer. Hence \[\begin{align} P\varphi_h(x)-P_0\varphi_h(x) &= \mathbb{E}\left[\varphi_h(Y_1)-\varphi_h(Y_0)\right] +(1-\lambda(x))\left(\mathbb{E}\left[\varphi_h(Y_0)\right]-\varphi_h(0)\right)\\ &\leq Ch\varphi_h(x) +(1-\lambda(x))\cdot C\varphi_h(x)\\ &\leq Ch\varphi_h(x). \end{align}\]
Suppose that \(x\notin\Omega_{b}\) and \(x\geq\sqrt{3h}\), then both operators use only their continuous branch, and therefore \[P\varphi_h(x)-P_0\varphi_h(x) = \mathbb{E}\left[\varphi_h(Y_1)-\varphi_h(Y_0)\right] \leq Ch\varphi_h(x).\] Taking \(C_2\) large enough proves 113 .

We now estimate \(P_0\varphi_h\) on \([0,\sqrt{3h}]\). On this interval \(w_h\) agrees with the test function in Lemma 4. The only difference between \(P_0\) and the SBM transition operator is the replacement of \(\tilde{\lambda}(x)\) by \(\lambda(x)\). Let \(\tilde P_0\) denote the same auxiliary operator as \(P_0\), but with \(\lambda\) replaced by \(\tilde{\lambda}\). We wish to show \(P_0\varphi_h(x) \approx \tilde P_0\varphi_h(x)\) so that we can reuse calculation done in Lemma 4. For \(x\in[0,\sqrt{3h}]\), together with 121 , \[P_0\varphi_h(x)-\tilde P_0\varphi_h(x) = \left(\lambda(x)-\tilde{\lambda}(x)\right) \left( \mathbb{E}\left[\varphi_h(Y_0)\right]-\varphi_h(0) \right),\] where \(Y_0=|x+Z| \leq 2\sqrt{3h}\). Using the Lipschitz bound for \(\varphi_h\) 114 and the lower bound 115 , \[\left| \mathbb{E}\left[\varphi_h(Y_0)\right]-\varphi_h(0) \right| \leq \mathbb{E}\left|\varphi_h(Y_0)-\varphi_h(0)\right| \leq C\sqrt h \leq C\sqrt h\,\varphi_h(x).\] Apply Lemma 21, \[|P_0\varphi_h(x)-\tilde P_0\varphi_h(x)| \leq Ch\varphi_h(x), \qquad x\in[0,\sqrt{3h}].\] We now estimate \(P_0\varphi_h(x)-\varphi_h(x)\) by the case-by-case calculation in Lemma 4.

For \(x\in[0,\sqrt{3h}/2]\), \[P_0\varphi_h(x)\leq \varphi_h(x)\left(1-f_h(x)+C_3h\right),\] where \[f_h(x):= \frac{s}{4}\cdot\frac{\tilde{\lambda}(x)}{2}\sqrt{3h} +s\cdot\frac{\tilde{\lambda}(x)x^2}{2\sqrt{3h}}.\] The proof of Lemma 4 shows that, for \(s>16\), \[f_h(x)\geq x,\qquad x\in[0,\sqrt{3h}/2].\]
Likewise, for \(x\in[\sqrt{3h}/2,\sqrt{3h}]\), \[P_0\varphi_h(x)\leq \varphi_h(x)\left(1-f_h(x)+C_3h\right),\] where \[f_h(x):= \frac{s\tilde{\lambda}(x)}{16}\sqrt{3h} +\frac{s\tilde{\lambda}(x)x^2}{4\sqrt{3h}} +\frac{s\tilde{\lambda}(x)x}{4} +\frac{s\sqrt{3h}}{2} -sx.\] Again, the proof of Lemma 4 shows that, for \(s>16\), \[f_h(x)\geq x,\qquad x\in[\sqrt{3h}/2,\sqrt{3h}].\]
The preceding SBM comparison covers \([0,\sqrt{3h}]\). The enlarged set \(B_h\) contains one additional strip of width \(O(h)\), namely \([\sqrt{3h},\sqrt{3h}+K_Bh]\). This strip appears only because the diffusion boundary layer is not exactly the SBM boundary layer, so it must be checked directly. It remains to treat the thin strip \([\sqrt{3h},\sqrt{3h}+K_Bh]\). On this interval \(P_0\) performs standard EM jump without reflection because \(x\geq \sqrt{3h}\). Also \(x+Z\in[0,R]\) for \(h\) sufficiently small, so \(w_h\) agrees with the original function used in Lemma 4. A direct calculation gives, uniformly for \(x\in[\sqrt{3h},\sqrt{3h}+K_Bh]\), \[\label{eq:keydiffusion95thin95strip95calc} \begin{align} P_0\varphi_h(x)-\varphi_h(x) &= \left(\frac{3}{4}-\frac{x}{2\sqrt{3h}}\right)e^M +\frac{e^M}{2\sqrt{3h}}\int_{\sqrt{3h}/2}^{\sqrt{3h}+x}e^{-s(z-\sqrt{3h}/2)}\,dz -e^{w_h(x)}\\ &= -e^{w_h(x)} \left[ \frac{s}{4} \left( \frac{x^2}{\sqrt{3h}}-3x+\frac{9}{4}\sqrt{3h} \right) +\mathcal{O}(h) \right]. \end{align}\tag{122}\] Thus on this thin strip define \[f_h(x):=\frac{s}{4} \left( \frac{x^2}{\sqrt{3h}}-3x+\frac{9}{4}\sqrt{3h} \right).\] For \(x\in[\sqrt{3h},\sqrt{3h}+K_Bh]\), \[f_h'(x)=\frac{s}{4}\left(\frac{2x}{\sqrt{3h}}-3\right)\leq -\frac{s}{8}\] after decreasing \(\delta_1\) so that \(2K_Bh/\sqrt{3h}\leq 1/2\). Hence \(f_h\) is decreasing on the thin strip. Therefore \[\label{eq:keydiffusion95thin95strip95lower} \begin{align} f_h(x) &\geq f_h(\sqrt{3h}+K_Bh)\\ &=\frac{s}{16}\sqrt{3h}-\frac{sK_B}{4}h+\frac{sK_B^2}{4}\frac{h^2}{\sqrt{3h}}. \end{align}\tag{123}\] Since \(s>16\), after decreasing \(\delta_1\) if necessary, \[f_h(\sqrt{3h}+K_Bh)>\sqrt{3h}+K_Bh.\] Hence \[f_h(x)\geq x,\qquad x\in[\sqrt{3h},\sqrt{3h}+K_Bh].\]

Combining the estimates for \(P_0\varphi_h(x)-\varphi_h(x)\) on \(B_h\) with 113 , we obtain \[P\varphi_h(x)\leq \varphi_h(x)(1-f_h(x)+Ch), \qquad x\in B_h.\]

We now prove the corresponding estimate outside \(B_h\). For \(x\notin B_h\), set \[f_h(x):=0.\] Since \(\Omega_{br}\subset B_h\), if \(x\notin B_h\), then the true scheme is outside the boundary layer and the standard Euler–Maruyama update is nonnegative for every possible \(Z\). Thus \[X_{k+1}=x+\sigma(x)Z+b(x)h.\] By Assumption 4, there exists \(L>0\) such that \[|b(x)|+|\sigma(x)|\leq L(1+x),\qquad x\geq 0.\] If we can assume global boundedness on drift and diffusion, this part is straightforward because the usual Taylor expansion gives an order-\(h\) one-step bound for the test function under bounded coefficients. Without global boundedness assumption, we analyze \(P\varphi_h(x)\) in \[B_h^c\cap[0,R/2],\qquad [R/2,4R],\qquad [4R,\infty).\] The first region exploits the piecewise linear structure of \(w(x)\) and boundedness of \(\sigma(x),\, b(x)\) in the compact domain. The second region exploits the smoothness of \(\varphi\) and the boundedness of drift and diffusion. The third region, \(\varphi\) is constant by construction and we can make sure the next step is within the constant region by linear growth assumption.

Consider \(x\in B_h^c\cap[0,R/2]\). On \([0,R]\), \(b\) and \(\sigma\) are bounded. For \(h\) sufficiently small, every possible value of \(x+\sigma(x)Z+b(x)h\) lies in \([0,R]\). Since \(x\notin B_h\), we also have \(x>\sqrt{3h}/2\). On \([0,R]\), the function \(w_h\) is flat-linear, and for \(y\in[0,R]\) and \(x>\sqrt{3h}/2\), \[w_h(y)-w_h(x)\leq -s(y-x).\] Therefore, by Taylor expansion around \(0\), using first and second moment of uniform random variable and the boundedness of \(\sigma\), \(b\) in \([0,R]\), there exists \(C_4\) is independent of \(x\) and \(h\) such that \[\frac{P\varphi_h(x)}{\varphi_h(x)} \leq \mathbb{E}\left[e^{-s(\sigma(x)Z+b(x)h)}\right] \leq 1+C_4h.\]
Consider \(x\in[R/2,4R]\). For \(h\) sufficiently small, every possible value of \[Y:=x+\sigma(x)Z+b(x)h\] lies in \([R/4,5R]\). Also \(\sqrt{3h}/2<R/4\), so \(\varphi_h\) is \(C^2\) on \([R/4,5R]\). There exist constants \(c_{\varphi},C_{\varphi}>0\), depending on \(s,M,R\) and \(\chi\) but not on \(x\) or \(h\), such that \[0<c_{\varphi}\leq \varphi_h(y), \qquad |\varphi_h'(y)|+|\varphi_h''(y)|\leq C_{\varphi}, \qquad y\in[R/4,5R].\] Taylor expansion gives \[\varphi_h(Y) = \varphi_h(x)+\varphi_h'(x)(Y-x)+\frac{1}{2}\varphi_h''(\xi)(Y-x)^2\] for some \(\xi\) between \(x\) and \(Y\). Taking expectations and using \(\mathbb{E}Z=0\) and \(\mathbb{E}Z^2=h\), \[|P\varphi_h(x)-\varphi_h(x)| \leq C h.\] Since \(\varphi_h(x)\geq c_{\varphi}\) on this compact interval, this implies \[P\varphi_h(x)\leq \varphi_h(x)(1+C_5h), \qquad x\in[R/2,4R],\] where \(C_5\) is independent of \(x\) and \(h\).
Consider \(x\geq 4R\). By the linear-growth assumption, for all possible \(Z\), \[x+\sigma(x)Z+b(x)h \geq x-L(1+x)(\sqrt{3h}+h).\] Since \(R\geq 1\), \(1+x\leq \frac{5}{4}x\) for \(x\geq 4R\). After decreasing \(\delta_1\) so that \[\frac{5}{4}L(\sqrt{3h}+h)\leq \frac{1}{2},\] we get \[x+\sigma(x)Z+b(x)h\geq \frac{x}{2}\geq 2R.\] Thus both \(x\) and \(X_{k+1}\) lie in the region where \(w_h\) is constant. Therefore \[P\varphi_h(x)=\varphi_h(x),\qquad x\geq 4R.\]

Taking \(C_1\) to be larger than the constants appearing in the estimates above, and decreasing \(\delta_1\) if necessary, we obtain 111 for all \(x\geq 0\) and all \(h<\delta_1\). Moreover, by construction, \[f_h(x)\geq x\mathbf{1}_{B_h}(x).\] Since \(f_h\) is supported on \(B_h\) and the formulas above are all of order \(\sqrt h\), there exists \(C_6>0\), independent of \(x\) and \(h\), such that \[\label{eq:keydiffusion95fh95upper} 0\leq f_h(x)\leq C_6\sqrt h,\qquad \forall x\geq 0,\quad h<\delta_1.\tag{124}\] We now introduce the time-dependent supersolution \[\label{eq:keydiffusion95V95define} V(t_i,x):=e^{K_T(T-t_i)}\varphi_h(x), \qquad i=0,\dots,N.\tag{125}\] The constant \(K_T>0\) will be chosen below. Using 111 , \[\label{eq:keydiffusion95PV95calc} \begin{align} PV(t_i,x)-V(t_i,x) &=e^{K_T(T-t_i-h)}P\varphi_h(x)-e^{K_T(T-t_i)}\varphi_h(x)\\ &\leq e^{K_T(T-t_i)}\varphi_h(x) \left[ e^{-K_Th}(1-f_h(x)+C_1h)-1 \right]. \end{align}\tag{126}\] After decreasing \(\delta_1\) if necessary, \[e^{-K_Th}=1-K_Th+r_{K_T}(h),\qquad |r_{K_T}(h)|\leq C_7(K_T)h^2.\] Using this expansion in 126 , together with 124 , gives \[PV(t_i,x)-V(t_i,x) \leq -e^{K_T(T-t_i)}\varphi_h(x) \left( f_h(x)+(K_T-C_1)h-C_8(K_T)h^{3/2} \right).\] Choose \[K_T:=C_1+1.\] After decreasing \(\delta_1\) again, we may assume \[(K_T-C_1)h-C_8(K_T)h^{3/2}\geq 0.\] Hence \[\label{eq:keydiffusion95PV95super} PV(t_i,x)-V(t_i,x) \leq -e^{K_T(T-t_i)}\varphi_h(x)f_h(x).\tag{127}\] Finally, decrease \(\delta_1\) once more so that \[M-s\left(\frac{\sqrt{3h}}{2}+K_Bh\right)\geq 0.\] Then \(w_h(x)\geq 0\) on \(B_h\), so \(e^{K_T(T-t_i)}\varphi_h(x)\geq 1\) on \(B_h\). Combining this with 112 and 127 , we conclude that \[\label{eq:keydiffusion95PV95final} PV(t_i,x)-V(t_i,x) \leq -x\mathbf{1}_{B_h}(x), \qquad i=0,\dots,N-1,\quad x\geq 0.\tag{128}\] Define \[U_B(t_i,x):= \mathbb{E}\left( \sum_{k=i}^{N-1}X_k\mathbf{1}_{B_h}(X_k) \,\middle|\, X_i=x \right).\] By Lemma 7, applied with \(g(t_i,x)=x\mathbf{1}_{B_h}(x)\), we have \[PU_B(t_i,x)-U_B(t_i,x)=-x\mathbf{1}_{B_h}(x), \qquad i=0,\dots,N-1,\] with terminal condition \(U_B(T,x)=0\). On the other hand, \(V(T,x)=\varphi_h(x)\geq 0=U_B(T,x)\), and 128 shows that \(V\) is a supersolution. A backward induction, exactly as in the proof of Lemma 4, gives \[U_B(t_i,x)\leq V(t_i,x), \qquad i=0,\dots,N,\quad x\geq 0.\] Since \(w_h(x)\leq M\) for all \(x\geq 0\), \[V(t_i,x)\leq e^{K_TT+M}, \qquad i=0,\dots,N,\quad x\geq 0.\] Therefore \[U_B(t_0,x)\leq e^{K_TT+M}, \qquad x\geq 0.\] Using \(\Omega_{br}\subset B_h\), we get \[\mathbb{E}\left( \sum_{k=0}^{N-1}X_k\mathbf{1}_{\Omega_{br}}(X_k) \,\middle|\, X_0=x \right) \leq U_B(t_0,x) \leq e^{K_TT+M}.\] If \(X_0\) is random, taking expectation and using the tower property gives \[\mathbb{E}\left( \sum_{k=0}^{N-1}X_k\mathbf{1}_{\Omega_{br}}(X_k) \right) \leq e^{K_TT+M}.\] Thus the desired bound holds with \(C=e^{K_TT+M}\), which is independent of \(h\). ◻

Proof of Lemma 11. It suffices to prove the one-step estimate \[\mathbb{E}_x(X_{k+1}^4)\le (1+Ch)x^4+Ch,\] for a constant \(C\) independent of \(x\) and \(h\).

First consider \(x\in\Omega_e\). In this region the scheme makes a standard Euler–Maruyama jump, and the absolute value is inactive, so \[X_{k+1}=x+b(x)h+\sigma(x)Z_{k+1}, \qquad Z_{k+1}\sim \mathrm{Unif}[-\sqrt{3h},\sqrt{3h}].\] Using \(\mathbb{E}Z_{k+1}=0\), \(\mathbb{E}Z_{k+1}^3=0\), \(\mathbb{E}Z_{k+1}^2=h\), and \(\mathbb{E}Z_{k+1}^4=\frac{9}{5}h^2\), we have \[\begin{align} \mathbb{E}_x(X_{k+1}^4) &= (x+b(x)h)^4 +6(x+b(x)h)^2\sigma^2(x)h +\frac{9}{5}\sigma^4(x)h^2\\ &= x^4 +4x^3b(x)h +6x^2b^2(x)h^2 +4xb^3(x)h^3 +b^4(x)h^4\\ &\quad +6x^2\sigma^2(x)h +12xb(x)\sigma^2(x)h^2 +6b^2(x)\sigma^2(x)h^3 +\frac{9}{5}\sigma^4(x)h^2. \end{align}\] By the linear-growth assumption, \(|b(x)|+|\sigma(x)|\le L(1+x)\). Hence, using \(h<1\) and \((1+x)^4\le C(1+x^4)\) for \(x\ge0\), each term after \(x^4\) is bounded by \(Ch(1+x^4)\). For example, \[|x^3b(x)h|\le Chx^3(1+x) \le 2Ch(1+x^4),\] \[x^2b^2(x)h^2+x^2\sigma^2(x)h\le Chx^2(1+x)^2\le 3Ch(1+x^4),\] and the remaining terms are bounded similarly by \(Ch(1+x)^4\le Ch(1+x^4)\). Therefore, \[\mathbb{E}_x(X_{k+1}^4)\le x^4+Ch(1+x^4)\le (1+Ch)x^4+Ch,\qquad x\in\Omega_e.\]

Next consider \(x\in\Omega_{br}\). By Lemma 8, after decreasing \(\delta\) if necessary, \(x\le C\sqrt h\) and \(b,\sigma\) are bounded on the boundary layer. The scheme either jumps to \(0\) or makes a reflected Euler–Maruyama-type update, so in either case \[X_{k+1}\le x+\sigma_2\sqrt{3h}+B_Mh\le C\sqrt h.\] Hence \[\mathbb{E}_x(X_{k+1}^4)\le Ch^2\le Ch \le (1+Ch)x^4+Ch, \qquad x\in\Omega_{br}.\] Combining the two cases gives \[\mathbb{E}_x(X_{k+1}^4)\le (1+Ch)x^4+Ch, \qquad x\ge 0.\] Taking expectation with respect to \(X_k\) yields \[\mathbb{E}(X_{k+1}^4)\le (1+Ch)\mathbb{E}(X_k^4)+Ch.\] By the discrete Gronwall inequality, for \(0\le k\le N\) and \(Nh=T\), \[\mathbb{E}(X_k^4)\le e^{CT}\left(\mathbb{E}(X_0^4)+CT\right).\] Therefore, \[\sup_{0\le k\le N}\mathbb{E}(1+X_k^4)\le C_T,\] where \(C_T\) is independent of \(h\). ◻

References↩︎

[1]

W. Feller, “The parabolic differential equations and the associated semi-groups of transformations,” Annals of Mathematics, 1952, [Online]. Available: http://www.jstor.org/stable/1969644.

[2]

W. Feller, Generalized second order differential operators and their lateral conditions, vol. 1. Illinois Journal of Mathematics, 1957.

[3]

A. D. Venttsel, “On boundary conditions for multidimensional diffusion processes,” Theory of Probability & Its Applications, vol. 4, no. 2, pp. 164–177, 1959, doi: 10.1137/1104014.

[4]

N. Ikeda, “On the construction of two-dimensional diffusion processes satisfying Wentzell’s boundary conditions and its application to boundary value problems,” Memoirs of the College of Science, University of Kyoto. Series A: Mathematics, vol. 33, no. 3, pp. 367–427, 1961, doi: 10.1215/kjm/1250711995.

[5]

K. Itô and H. P. McKean, “Brownian motions on a half line,” Illinois Journal of Mathematics, vol. 7, no. 2, pp. 181–231, 1963, doi: 10.1215/ijm/1255644633.

[6]

D. W. Stroock and S. R. S. Varadhan, Gives sticky + oblique reflection condition for multidimensional diffusions. Proves uniqueness under certain conditions.“Diffusion processes with boundary conditions,” Communications on Pure and Applied Mathematics, vol. 24, no. 2, pp. 147–225, 1971, doi: 10.1002/cpa.3160240206.

[7]

S. Watanabe, “On stochastic differential equations for multi-dimensional diffusion processes with boundary conditions,” Journal of Mathematics of Kyoto University, vol. 11, no. 1, pp. 169–180, 1971, doi: 10.1215/kjm/1250523692.

[8]

G. Peskir, On Boundary Behaviour of One-Dimensional Diffusions: From Brown to Feller and Beyond. Selected Papers II, 2015.

[9]

A. Anagnostakis, “Pricing and hedging for a sticky diffusion,” Journal of Applied Probability, vol. 62, no. 4, pp. 1235–1259, 2025, doi: 10.1017/jpr.2025.2.

[10]

Y. Kabanov, M. Kijima, and S. Rinaz, “A positive interest rate model with sticky barrier,” Quantitative Finance, vol. 7, no. 3, pp. 269–284, 2007, doi: 10.1080/14697680600999351.

[11]

Y. Nie and V. Linetsky, “Sticky reflecting Ornstein-Uhlenbeck diffusions and the Vasicek interest rate model with the sticky zero lower bound,” Stochastic Models, vol. 36, no. 1, pp. 1–19, 2020, doi: 10.1080/15326349.2019.1630287.

[12]

F. A. Longstaff, “Multiple equilibria and term structure models,” Journal of Financial Economics, vol. 32, no. 3, pp. 333–344, 1992, doi: 10.1016/0304-405x(92)90031-r.

[13]

R. J. Baxter, “Percus–Yevick Equation for Hard Spheres with Surface Adhesion,” The Journal of Chemical Physics, vol. 49, no. 6, p. 2770, 1968, doi: 10.1063/1.1670482.

[14]

G. Stell, “Sticky spheres and related systems,” Journal of Statistical Physics, vol. 63, no. 5–6, pp. 1203–1221, 1991, doi: 10.1007/bf01030007.

[15]

M. A. Miller and D. Frenkel, “Phase diagram of the adhesive hard sphere fluid,” The Journal of Chemical Physics, vol. 121, no. 1, p. 535, 2004, doi: 10.1063/1.1758693.

[16]

G. Meng, N. Arkus, M. P. Brenner, and V. N. Manoharan, “The Free-Energy Landscape of Clusters of Attractive Hard Spheres,” Science, vol. 327, no. 5965, pp. 560–563, Jan. 2010, doi: 10.1126/science.1181263.

[17]

M. Holmes-Cerfon, S. J. Gortler, and M. P. Brenner, “A geometrical approach to computing free-energy landscapes from short-ranged potentials,” Proceedings of the National Academy of Sciences, vol. 110, no. 1, pp. E5–E14, 2013, doi: 10.1073/pnas.1211720110.

[18]

R. W. Perry, M. C. Holmes-Cerfon, M. P. Brenner, and V. N. Manoharan, “Two-Dimensional Clusters of Colloidal Spheres: Ground States, Excited States, and Structural Rearrangements,” Physical Review Letters, vol. 114, no. 22, pp. 228301–5, Jun. 2015, doi: 10.1103/physrevlett.114.228301.

[19]

A. Gandolfi, A. Gerardi, and F. Marchetti, “Association rates of diffusion-controlled reactions in two dimensions,” Acta Applicandae Mathematicae, vol. 4, no. 2–3, pp. 139–155, 1985, doi: 10.1007/bf00052459.

[20]

P. C. Bressloff, “Close encounters of the sticky kind: Brownian motion at absorbing boundaries,” Physical Review E, vol. 107, no. 6, p. 064121, 2023, doi: 10.1103/physreve.107.064121.

[21]

A. Calsina and J. Z. Farkas, “Steady states in a structured epidemic model with Wentzell boundary condition,” Journal of Evolution Equations, vol. 12, no. 3, pp. 495–512, Mar. 2012, doi: 10.1007/s00028-012-0142-6.

[22]

A. Pugliese and F. Milner, “A structured population model with diffusion in structure space,” Journal of Mathematical Biology, vol. 77, no. 6–7, pp. 2079–2102, 2018, doi: 10.1007/s00285-018-1246-6.

[23]

E. H. van Nes, D. G. F. Pujoni, S. A. Shetty, G. Straatsma, W. M. de Vos, and M. Scheffer, “A tiny fraction of all species forms most of nature: Rarity as a sticky state,” Proceedings of the National Academy of Sciences, vol. 121, no. 2, p. e2221791120, 2024, doi: 10.1073/pnas.2221791120.

[24]

J. M. Harrison and A. J. Lemoine, “Sticky Brownian motion as the limit of storage processes,” Journal of Applied Probability, vol. 18, no. 1, pp. 216–226, Jul. 2016, doi: 10.2307/3213181.

[25]

J. Berry and F. Camilli, “Stationary Mean Field Games on networks with sticky transition conditions,” ESAIM: Control, Optimisation and Calculus of Variations, vol. 31, p. 18, 2025, doi: 10.1051/cocv/2025007.

[26]

C. J. Howitt, “Stochastic Flows and Sticky Brownian motion,” PhD thesis, University of Warwick, 2007.

[27]

A. Eberle and R. Zimmer, “Sticky couplings of multidimensional diffusions with different drifts,” Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, vol. 55, no. 4, pp. 2370–2394, 2019, doi: 10.1214/18-aihp951.

[28]

A. Durmus, A. Eberle, A. Guillin, and K. Schuh, “Sticky nonlinear SDEs and convergence of McKean–Vlasov equations without confinement,” Stochastics and Partial Differential Equations: Analysis and Computations, vol. 12, no. 3, pp. 1855–1906, 2024, doi: 10.1007/s40072-023-00315-8.

[29]

P. E. Kloeden and E. Platen, Numerical solution of stochastic differential equations, vol. 23. Springer-Verlag, Berlin, 1992.

[30]

E. Peters and T. Barenbrug, “Efficient Brownian dynamics simulation of particles near walls. II. Sticky walls,” Physical Review E, vol. 66, no. 5, p. 056702, Nov. 2002, doi: 10.1103/physreve.66.056702.

[31]

R. Fantoni and P. Sollich, “Multicomponent adhesive hard sphere models and short-ranged attractive interactions in colloidal or micellar solutions,” Physical Review E, vol. 74, no. 5, p. 051407, Nov. 2006, doi: 10.1103/physreve.74.051407.

[32]

N. Bou-Rabee and M. C. Holmes-Cerfon, “Sticky Brownian Motion and Its Numerical Solution,” SIAM Review, vol. 62, no. 1, pp. 164–195, Feb. 2020, doi: 10.1137/19m1268446.

[33]

N. Bou-Rabee and E. Vanden-Eijnden, Continuous-time Random Walks for the Numerical Solution of Stochastic Differential Equations, vol. 256. 2018.

[34]

C. Meier, L. Li, and G. Zhang, “Markov Chain Approximation of One-Dimensional Sticky Diffusions,” Advances in Applied probability, vol. 53, no. 2, pp. 335–369, 2021.

[35]

C. Meier, L. Li, and G. Zhang, “Simulation of multidimensional diffusions with sticky boundaries via Markov chain approximation,” European Journal of Operational Research, vol. 305, no. 3, pp. 1292–1308, 2023, doi: 10.1016/j.ejor.2022.07.038.

[36]

J. A. Anderson, J. Glaser, and S. C. Glotzer, “HOOMD-blue: A Python package for high-performance molecular dynamics and hard particle Monte Carlo simulations,” Computational Materials Science, vol. 173, p. 109363, 2020, doi: 10.1016/j.commatsci.2019.109363.

[37]

A. P. Thompson et al., “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Computer Physics Communications, vol. 271, p. 108171, 2022, doi: 10.1016/j.cpc.2021.108171.

[38]

D. Frenkel and B. Smit, Understanding Molecular Simulation. Elsevier, 2023.

[39]

A. Sharma, “Weak approximation of stochastic differential equations with sticky boundary conditions,” arXiv, 2025, doi: 10.48550/arxiv.2508.06487.

[40]

B. Leimkuhler, A. Sharma, and M. V. Tretyakov, “Simplest random walk for approximating Robin boundary value problems and ergodic limits of reflected diffusions,” The Annals of Applied Probability, vol. 33, pp. 1904–1960, 2023, doi: 10.1214/22-aap1856.

[41]

E. Barth, K. Kuczera, B. Leimkuhler, and R. D. Skeel, “Algorithms for constrained molecular dynamics,” Journal of Computational Chemistry, vol. 16, no. 10, pp. 1192–1209, 1995, doi: 10.1002/jcc.540161003.

[42]

E. Vanden-Eijnden and G. Ciccotti, “Second-order integrators for Langevin equations with holonomic constraints,” Chemical Physics Letters, vol. 429, no. 1–3, pp. 310–316, Sep. 2006, doi: 10.1016/j.cplett.2006.07.086.

[43]

G. Ciccotti, T. Lelièvre, and E. Vanden-Eijnden, “Projection of diffusions on submanifolds: Application to mean force computation,” Communications on Pure and Applied Mathematics, vol. 61, no. 3, pp. 371–408, 2007, doi: 10.1002/cpa.20210.

[44]

N. Bou-Rabee, A. Donev, and E. Vanden-Eijnden, “Metropolis Integration Schemes for Self-Adjoint Diffusions,” Multiscale Modeling & Simulation, vol. 12, no. 2, pp. 781–831, Apr. 2014, doi: 10.1137/130937470.

[45]

B. Leimkuhler and C. Matthews, “Efficient molecular dynamics using geodesic integration and solvent–solute splitting,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 472, no. 2189, pp. 20160138–22, May 2016, doi: 10.1098/rspa.2016.0138.

[46]

N. Ikeda and S. Watanabe, Stochastic Differential Equations and Diffusion Processes. 1981.

[47]

G. N. Milstein and M. V. Tretyakov, “Stochastic Numerics for Mathematical Physics,” Scientific Computation, 2021, doi: 10.1007/978-3-030-82040-4.

[48]

K.-J. Engel, “The laplacian on \(C(\overline{\Omega})\) with generalized wentzell boundary conditions,” Archive for Rational Mechanics and Analysis, vol. 168, no. 3, pp. 247–261, 2003.

[49]

A. Pazy, Semigroups of linear operators and applications to partial differential equations. Springer Science & Business Media, 2012.

[50]

M. Holmes-Cerfon, “Simulating sticky particles: A Monte Carlo method to sample a stratification,” Journal of Chemical Physics, vol. 153, no. 16, p. 164112, Oct. 2020, doi: 10.1063/5.0019550.

[51]

B. Leimkuhler, A. Sharma, and M. V. Tretyakov, “Simplest random walk for approximating robin boundary value problems and ergodic limits of reflected diffusions,” The Annals of Applied Probability, vol. 33, no. 3, pp. 1904–1960, 2023.

If we can make assumptions like compact domain, making drift and diffusion uniformly bounded, the \(\Omega_e\) local error becomes straightforward to show.↩︎

A modified Euler-Maruyama method to simulate a one-dimensional sticky diffusion

Abstract

1 Introduction↩︎

2 A brief introduction to one-dimensional sticky diffusion processes↩︎

2.1 Stochastic differential equation↩︎

2.2 Backward equation↩︎

2.3 Stationary distribution↩︎

2.4 Derivation via conservation of probability↩︎

2.5 Equivalence of the forward and backward equations↩︎

3 Sticky Brownian motion↩︎

3.1 Overview of the simulation algorithm↩︎

3.2 How to choose \(\lambda(x)\)↩︎

3.3 Primary convergence results↩︎

3.4 Convergence proofs↩︎

3.5 Sufficient conditions for backward equation regularity↩︎

4 General Sticky Diffusion↩︎

4.1 Overview of the simulation algorithm↩︎

4.2 How to choose \(\lambda(x)\)↩︎

4.3 \(\lambda\) is a probability↩︎

4.4 Primary convergence results↩︎

4.5 Numerical experiments↩︎

4.6 Key steps in the error analysis↩︎

4.6.1 Local error: Theorem 5 ↩︎

4.6.2 Bounding the number of steps: Lemma 10 ↩︎

5 Conclusion↩︎

6 Supplement↩︎

6.1 Additional proofs from Section 3↩︎

6.2 Additional proofs from Section 4↩︎

References↩︎

Subjects

Updated on Academus

A modified Euler-Maruyama method to simulate a one-dimensional sticky diffusion

Abstract

1 Introduction↩︎

2 A brief introduction to one-dimensional sticky diffusion processes↩︎

2.1 Stochastic differential equation↩︎

2.2 Backward equation↩︎

2.3 Stationary distribution↩︎

2.4 Derivation via conservation of probability↩︎

2.5 Equivalence of the forward and backward equations↩︎

3 Sticky Brownian motion↩︎

3.1 Overview of the simulation algorithm↩︎

3.2 How to choose \(\lambda(x)\)↩︎

3.3 Primary convergence results↩︎

3.4 Convergence proofs↩︎

3.5 Sufficient conditions for backward equation regularity↩︎

4 General Sticky Diffusion↩︎

4.1 Overview of the simulation algorithm↩︎

4.2 How to choose \(\lambda(x)\)↩︎

4.3 \(\lambda\) is a probability↩︎

4.4 Primary convergence results↩︎

4.5 Numerical experiments↩︎

4.6 Key steps in the error analysis↩︎

4.6.1 Local error: Theorem 5↩︎

4.6.2 Bounding the number of steps: Lemma 10↩︎

5 Conclusion↩︎

6 Supplement↩︎

6.1 Additional proofs from Section 3↩︎

6.2 Additional proofs from Section 4↩︎

References↩︎

Subjects

Updated on Academus

4.6.1 Local error: Theorem 5 ↩︎

4.6.2 Bounding the number of steps: Lemma 10 ↩︎