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

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.

• The paper introduces a modified reflected Euler–Maruyama scheme that jumps exactly to the sticky boundary with a generator-consistent probability, avoiding the boundary-missing problem of standard EM methods.
• The analysis proves an O(h) weak error bound for sticky Brownian motion and general one-dimensional sticky diffusions by controlling boundary-layer occupation and local truncation errors.
• Numerical experiments confirm near-first-order convergence, with fitted rates from 0.908 to 1.027, while the method remains easier to integrate into particle simulations than many CTMC alternatives.

Motivation and contribution

Sticky diffusions—processes that spend positive Lebesgue time on a lower-dimensional boundary—have re-emerged as useful models in finance (interest rates pinned near zero), colloidal physics with short-ranged adhesion, cell biology, epidemiology, ecology, and queueing/network models. Simulating them is difficult because standard Euler–Maruyama (EM) schemes take random spatial steps and essentially never land exactly on the boundary, while the alternative continuous-time Markov chain (CTMC) approximations of Bou-Rabee–Holmes-Cerfon and Meier et al., though effective in low dimensions, require rate constants that scale poorly with dimension, need a valid spatial discretization, and are awkward to integrate into existing particle-simulation frameworks.

This paper by Jiang and Holmes-Cerfon proposes a simple modification of the reflected EM scheme for a one-dimensional sticky diffusion on $[0,\infty)$ that can jump exactly to the sticky point at the origin. The main result is a proof of weak convergence of order 1: for bounded test functions $\phi$ whose backward solution is sufficiently regular,

$\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$

where $\{X_k\}$ is the numerical chain and $Y_t$ the true sticky process.

Background: characterizations of one-dimensional sticky diffusions

The paper first collects equivalent formal descriptions of a sticky diffusion with drift $b$, diffusion coefficient $\sigma$ ($a=\sigma^2/2$), and stickiness parameter $\kappa>0$. The process solves the SDE system

$$dX_t = b(X_t)1_D(X_t)\,dt + \sigma(X_t)1_D(X_t)\,dW_t + dL_t,\qquad \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)\,dt = \kappa\, dL_t,$$

with local time $\phi$0. Its generator has domain restricted by the Wentzell-type sticky boundary condition

$\phi$1

The authors also derive, via probability conservation, a forward equation for a measure of the form $\phi$2: the density $\phi$3 evolves under the Fokker–Planck operator while the atom exchanges mass with the interior through the flux $\phi$4 evaluated at zero, closed by the condition $\phi$5. They verify formally that this forward formulation is adjoint to the backward equation, and note that the stationary distribution carries an atom proportional to $\phi$6 at the boundary. These characterizations motivate the design principle of the scheme: match the discrete transition operator to the generator including its boundary condition.

The modified Euler–Maruyama scheme

For the sticky Brownian motion (SBM), the domain is split into a boundary layer $\phi$7 and exterior $\phi$8. In the exterior the chain takes standard EM increments with uniform noise $\phi$9 (chosen over Gaussians to simplify analysis). In the boundary layer, the chain performs a reflected EM update $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$0 with probability $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$1, and jumps exactly to $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$2 with probability $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$3:

$\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$4

The function $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$5 is derived by requiring the one-step transition operator to approximate the generator: Taylor-expanding $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$6 and using the sticky boundary condition $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$7 to replace $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$8 by $\left|E[\phi(X_N)] - E[\phi(Y_T)]\right| \le Ch,$9 up to $\{X_k\}$0 error, one obtains a moment-matching condition solved exactly by the formula above. A notable subtlety the authors report: dropping the intermediate term $\{X_k\}$1 in the boundary-condition expansion yields a $\{X_k\}$2 that fails to be a valid probability; retaining it restores validity. Indeed, $\{X_k\}$3 is shown to satisfy $\{X_k\}$4 on the layer, increasing in $\{X_k\}$5 and decreasing in $\{X_k\}$6, with $\{X_k\}$7.

Convergence analysis for the SBM

The global error is decomposed into local errors weighted by occupation counts in the two regions. Under uniform bounds on the first four spatial derivatives and second time derivative of the backward solution $\{X_k\}$8, the local error satisfies

$\{X_k\}$9

A naive worst-case bound would give only $Y_t$0 globally, since the boundary-layer local error is $Y_t$1 and there are $Y_t$2 steps. The key technical lemma shows instead that the expected occupation measure of the layer is uniformly bounded:

$Y_t$3

independent of $Y_t$4 and initial condition. Since points in the layer are $Y_t$5-close to the origin, this implies only $Y_t$6 effective visits, yielding global order 1. The proof constructs an exponential supersolution $Y_t$7 with $Y_t$8 piecewise linear (flat near the boundary, decaying away from it), verifies a one-step drift inequality case-by-case against the transition operator, and closes via backward induction—a Lyapunov/supersolution argument in the spirit of Leimkuhler–Sharma–Tretyakov's analysis of Robin-reflected diffusions.

Because the convergence theorem is stated in terms of regularity of $Y_t$9 rather than of $b$0, the paper supplies a sufficient condition: if $b$1 satisfies the compatibility conditions $b$2, $b$3, $b$4, $b$5 on a compact interval with reflecting endpoint, then $b$6 for the Wentzell generator $b$7, and semigroup theory gives the required derivative bounds. The argument extends verbatim to general Wentzell conditions at $b$8.

Extension to general drift and diffusion

For a general sticky diffusion, the boundary layer becomes state-dependent, $b$9, which may be disconnected; the authors show it is contained in $\sigma$0 for small $\sigma$1 under local $\sigma$2 regularity of $\sigma$3 near 0, $\sigma$4, and linear growth $\sigma$5. The jump probability generalizes to

$\sigma$6

derived from the boundary condition $\sigma$7 at 0, assuming the nondegeneracy $\sigma$8 so that $\sigma$9 stays bounded. Validity of $a=\sigma^2/2$0 as a probability is proven for small $a=\sigma^2/2$1 via careful sign analysis across three regimes of $a=\sigma^2/2$2 versus $a=\sigma^2/2$3.

Two additional complications are handled explicitly. First, negative drift creates a thin reflection strip $a=\sigma^2/2$4 where the scheme performs a pure reflected EM step with no jump to the boundary; the first conditional moment there still approximates $a=\sigma^2/2$5 to order $a=\sigma^2/2$6, preserving the $a=\sigma^2/2$7 local error. Second, when $a=\sigma^2/2$8 the function $a=\sigma^2/2$9 degenerates; a separate argument exploiting the fact that the sticky boundary condition then forces $\kappa>0$0 recovers the same local error rate. The occupation bound for the enlarged layer $\kappa>0$1 is obtained by rescaling so $\kappa>0$2, showing the diffusion jump probability differs from the SBM one by $\kappa>0$3, coupling the transition operators, and adapting the exponential supersolution with a cutoff that avoids imposing global boundedness on the coefficients. The exterior local error is $\kappa>0$4, controlled via a fourth-moment Gronwall bound, giving the global estimate $\kappa>0$5.

Numerical validation

Monte Carlo experiments confirm first-order weak convergence. For the SBM with $\kappa>0$6, $\kappa>0$7, comparing $\kappa>0$8 against a finely resolved finite-difference solution of the backward equation, least-squares fits give slopes of approximately 0.957 ($\kappa>0$9) and 1.011 ($$dX_t = b(X_t)1_D(X_t)\,dt + \sigma(X_t)1_D(X_t)\,dW_t + dL_t,\qquad \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)\,dt = \kappa\, dL_t,$$0). For general diffusions with four combinations of drift and diffusion—including strongly negative drift $$dX_t = b(X_t)1_D(X_t)\,dt + \sigma(X_t)1_D(X_t)\,dW_t + dL_t,\qquad \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)\,dt = \kappa\, dL_t,$$1 with weak diffusion—the fitted slopes range from 0.908 to 1.027, with the strongly drifted case converging more slowly but exhibiting absolute errors several orders of magnitude smaller. Notably, the test functional $$dX_t = b(X_t)1_D(X_t)\,dt + \sigma(X_t)1_D(X_t)\,dW_t + dL_t,\qquad \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)\,dt = \kappa\, dL_t,$$2 does not satisfy the compatibility conditions required by the sufficient regularity proposition, yet empirical order 1 persists—an observation the authors flag but do not resolve theoretically.

Limitations and open questions

Several restrictions are acknowledged. The method is developed and proven only in one dimension; extension to multidimensional sticky diffusions—with genuinely different dynamics along the boundary—is anticipated but unproven, and the authors note that ensuring switching probabilities remain valid probabilities was already nontrivial even in 1D. The derivation assumes $$dX_t = b(X_t)1_D(X_t)\,dt + \sigma(X_t)1_D(X_t)\,dW_t + dL_t,\qquad \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)\,dt = \kappa\, dL_t,$$3 for the generic case, with the degenerate case treated separately but only under positivity of the drift. The convergence constants may depend on the initial position for general diffusions, and random initial conditions require integrability of the resulting bound. The regularity assumptions on the backward solution are verified only on compact intervals; a comparable regularity theory for the sticky/Wentzell problem on the half-line is not readily available in the literature. Open questions raised include higher-order schemes, characterization of the full class of admissible jump probabilities (including deliberately simpler ones achieving only $$dX_t = b(X_t)1_D(X_t)\,dt + \sigma(X_t)1_D(X_t)\,dW_t + dL_t,\qquad \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)\,dt = \kappa\, dL_t,$$4), long-time behavior and accuracy of stationary distributions, application to manifolds and stratified boundaries, and probabilistic solution of parabolic PDEs with sticky boundary conditions.

Conclusion

The paper provides a minimal modification of reflected Euler–Maruyama that faithfully simulates one-dimensional sticky diffusions by allowing exact jumps to the boundary with a generator-consistent probability. With rigorous proofs of $$dX_t = b(X_t)1_D(X_t)\,dt + \sigma(X_t)1_D(X_t)\,dW_t + dL_t,\qquad \frac{\sigma^2(X_t)}{2}1_{\partial D}(X_t)\,dt = \kappa\, dL_t,$$5 weak convergence for both sticky Brownian motion and general sticky diffusions, supported by extensive Monte Carlo evidence, it offers a practical alternative to CTMC discretizations that integrates naturally with existing constrained-diffusion simulation methodology.