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

Technical Summary: A Modified Euler-Maruyama Method to Simulate a One-Dimensional Sticky Diffusion

Problem Statement
Sticky diffusions are stochastic processes capable of spending a non-zero amount of time on a lower-dimensional boundary (e.g., a point in one dimension). While mathematically well-studied since Feller's work in the 1950s, simulating these processes presents a unique challenge: standard numerical methods for Stochastic Differential Equations (SDEs), such as the Euler-Maruyama (EM) scheme, take random steps in continuous space and thus almost never land exactly on the boundary. Consequently, they fail to capture the singular measure (the "stickiness") inherent to the process. Existing alternatives, such as continuous-time Markov chain approximations on pre-defined grids, face scalability issues in higher dimensions and difficulties handling curved boundaries. Furthermore, recent time-change methods based on reflected diffusions are limited to cases where the drift and diffusion coefficients are identical in the interior and on the boundary.

Methodology
The authors propose a modified Euler-Maruyama algorithm designed to simulate one-dimensional sticky diffusions that can jump exactly to and from a sticky point (taken as $x=0$). The method is a minimal modification of the standard reflected EM scheme, introducing a probabilistic mechanism to handle the boundary layer.

The algorithm divides the domain into two regions based on the time step $h$:

1. Exterior ($\Omega_e$): When the process is far from the boundary, it evolves using a standard reflected Euler-Maruyama update.
2. Boundary Layer ($\Omega_b$): When the process is close to the boundary (specifically within a distance scaling with $\sqrt{h}$), the algorithm performs a mixed update. With probability $\lambda(x)$, it performs a reflected EM step. With probability $1 - \lambda(x)$, it jumps exactly to the sticky point ($x=0$).

The core of the methodology lies in the derivation of the jump probability function $\lambda(x)$. The authors derive $\lambda(x)$ by requiring that the transition operator of the discrete scheme approximates the transition operator of the continuous sticky diffusion over a single time step $h$. This is achieved by:

• Expanding the generator of the sticky diffusion and the discrete scheme using Taylor series.
• Utilizing the sticky boundary condition (a relation between the first and second spatial derivatives at the boundary) to express the first derivative in terms of the second.
• Solving for $\lambda(x)$ such that the coefficients of the drift and diffusion terms in the local error expansion match the continuous dynamics up to a specified tolerance.

For a Sticky Brownian Motion (SBM) with drift $b=0$ and diffusion $\sigma=1$, the derived probability is:
$$\lambda(x) = \frac{x^2 + 2\kappa x + h(\frac{\kappa}{\sqrt{3h}} + 1)}{x^2 + h + \kappa\sqrt{3h}}$$
where $\kappa$ is the stickiness parameter. The method is subsequently generalized to arbitrary one-dimensional sticky diffusions with drift $b(x)$ and diffusion $\sigma(x)$, resulting in a more complex expression for $\lambda(x)$ that accounts for the state-dependent coefficients.

Key Contributions

1. Algorithm Design: The paper introduces a simple, computationally efficient algorithm that allows a discrete process to reside exactly on the boundary, overcoming the limitation of standard EM methods.
2. Probability Derivation: The authors provide a rigorous derivation for the jump probability $\lambda(x)$, ensuring the discrete scheme is consistent with the generator of the desired sticky diffusion dynamics.
3. Convergence Analysis: The paper proves that the proposed scheme converges weakly with order 1 (global error $O(h)$).
   • The proof decomposes the global error into local errors and the number of steps spent in each region.
   • In the exterior, the local error is $O(h^2)$.
   • In the boundary layer, the local error is $O(h^{3/2})$.
   • Crucially, the authors prove that the number of steps spent in the boundary layer scales as $O(h^{-1/2})$. The product of the local error and the step count yields the global $O(h)$ convergence.
4. Regularity Conditions: The paper establishes sufficient conditions on the terminal function $\phi$ (compatibility conditions at the boundaries) to ensure the regularity of the backward equation solution required for the convergence proofs.

Results

• Theoretical: The authors prove weak convergence of order 1 for both Sticky Brownian Motion and general sticky diffusions with state-dependent drift and diffusion. The proofs rely on bounding the local error and demonstrating that the occupation time of the boundary layer is controlled.
• Numerical: Numerical experiments validate the theoretical findings. Simulations of SBM and general sticky diffusions with various drift and diffusion coefficients demonstrate empirical convergence rates close to 1. The results show that the method effectively captures the time spent near the sticky boundary, with trajectories exhibiting longer dwell times as the stickiness parameter $\kappa$ increases.

Significance and Claims
The authors position this work as a significant step toward a general method for simulating arbitrary sticky diffusions. They claim that their method is advantageous because:

• It is a small modification of standard simulation methodologies, making it easier to incorporate into existing frameworks (e.g., those handling interacting particle systems or chemical reactions).
• It does not require the process to be a time-change of a reflected diffusion, thereby allowing for different drift and diffusion coefficients on the boundary compared to the interior.
• It avoids the computational bottlenecks of grid-based Markov chain approximations in higher dimensions.

The paper concludes by noting that while the current work focuses on one dimension, the heuristics and convergence tools described provide a pathway for extending these ideas to multi-dimensional sticky diffusions, though ensuring the validity of the switching probabilities in higher dimensions remains a technical challenge. The authors also suggest future work could explore higher-order schemes and applications to processes on manifolds.