Randomly shifted lattice rules for SDEs with constant diffusion: error bounds uniform in time discretization and dimension

Introduction

Stochastic differential equations (SDEs) model systems subject to random perturbations. They are widely used in biology, chemistry, mechanics, and finance [9]. Many applications require computing expectations of functions of SDE solutions. For example, such expectations are used to price financial derivatives [6] and to solve parabolic partial differential equations [19].

In this paper, we study the numerical approximation of such expectations, focusing on SDEs with constant diffusion. Specifically, let $I = \mathbb{E}[g(X_T)]$, where $X$ solves

Here $T > 0$ is fixed, $\boldsymbol{X}_t \in\mathbb{R}^{d_x}$, and $\boldsymbol{W}$ is a standard $d_w$-dimensional Brownian motion. The drift is $b : [0,T] \times\mathbb{R}^{d_x} \to\mathbb{R}^{d_x}$, the diffusion matrix $\Sigma\in\mathbb{R}^{d_x \times d_w}$ is constant, and $g : \mathbb{R}^{d_x} \to\mathbb{R}$ is the payoff function. The initial state $x_0 \in\mathbb{R}^{d_x}$ is deterministic. SDEs of this form arise in various applications, including overdamped Langevin dynamics [26].

Explicit solutions of SDEs are generally unavailable, so numerical schemes are used to approximate their sample paths. We use the Euler–Maruyama (EM) scheme [10] on the uniform grid $\{t_i = ih : i = 0, \ldots, m\}$, with $h = T/m$, and set $s = d_wm$. For $z = (z_1, \ldots, z_m) \in(\mathbb{R}^{d_w})^m \cong\mathbb{R}^s$, set $\widetilde{X}_{t_0}(z) = x_0$ and define

Let $Z \sim N(0,I_{d_s})$ be an $s$-dimensional standard Gaussian vector. Then $\widetilde{X}_{t_m}(Z)$ is the EM approximation at time $T$. The corresponding approximation to $\mathcal{I}$ can be written as

where $\Phi$ is the standard Gaussian distribution function and $\Phi^{-1}$ acts componentwise. Under suitable regularity and growth conditions, $\mathcal{I}_m$ approximates $\mathcal{I}$ with the classical weak error bound $|\mathcal{I}_m-\mathcal{I}| \leq Cm^{-1}$, where $C$ is independent of $m$; see [10 Theorem 14.5.1]. Thus, increasing the number of time steps $m$ reduces the discretization error.

To compute $\mathcal{I}_m$, however, we typically approximate the integral in (3) numerically, since it generally has no closed-form expression. In this paper, we focus on the error introduced by numerical integration, referred to as the sampling error, as distinct from the time discretization error $|\mathcal{I}_m-\mathcal{I}|$. The integration dimension $s=d_wm$ grows as the time grid is refined and becomes arbitrarily large as $m \to\infty$. This growth poses a central challenge in obtaining sampling error bounds that remain uniform in $m$. Monte Carlo (MC) provides a straightforward approach. For independent uniform points $u_0,\ldots,u_{N-1} \in(0,1)^s$, the MC estimator is

The estimator achieves a root mean square error of $O(N^{-1/2})$, with an implied constant independent of $m$ provided that the variance of the integrand is uniformly bounded. This convergence rate is nevertheless slow. Quasi-Monte Carlo (QMC) replaces independent random points with low-discrepancy point sets and achieves nearly first-order convergence under bounded-variation assumptions [2, 17, 4]. However, $\Phi^{-1}$ in (3) can make the integrand unbounded near the cube boundary, placing it outside the scope of the classical theory. To address this limitation, approaches based on boundary growth conditions [23], projection [22], and weighted function spaces [16] have been developed. Nevertheless, the boundary-growth and projection-based approaches do not directly resolve the dimension dependence.

The weighted-space framework of Nichols and Kuo [16] allows dimension-independent bounds, but applying it here requires mixed derivative estimates and summability conditions on the function-space weights. The nested composition structure of the EM integrand makes such estimates difficult to obtain. It remains unclear whether the dependence of this integrand on the Gaussian inputs permits weights satisfying these conditions uniformly in $m$. Recent studies offer insights into possible low-dimensional structure that may help address this question. For EM approximations, Ouyang, Pan, and He [20] derive randomized QMC (RQMC) error bounds of the form $C_{\varepsilon,s}N^{-1+\varepsilon}$, but the constants grow with $s$. Ouyang and Wang [21] later establish effective-dimension bounds uniform in $m$. Numerical experiments in both studies show advantages of RQMC over MC. These findings suggest that EM integrands may possess low-dimensional structure favorable to QMC, raising the possibility that suitable weights satisfy the required summability conditions.

It is also worth noting that QMC performance depends on how Brownian paths are constructed from Gaussian inputs. Different constructions can be represented by replacing $\Phi^{-1}(u)$ in (1.3) with $U\Phi^{-1}(u)$, where $U \in\mathbb{R}^{s \times s}$ is orthogonal. Such transformations leave the standard MC error unchanged but alter the structure of the integrand, and hence can affect the QMC error [12, 27]. In particular, we focus on the Brownian bridge (BB) construction, which first generates $W_T$ and then samples conditional midpoint values on successively finer intervals [15 Section 1.1.2]. Such a construction assigns the largest-scale variation of the Brownian path to the Gaussian inputs generating $W_T$, while each subsequent refinement introduces Gaussian inputs that affect smaller-scale fluctuations with decreasing amplitude. This hierarchical structure may lead to low-dimensional structure in the integrand that is favorable to QMC.

In this paper, we use randomly shifted rank-1 lattice rules, an RQMC method whose points are generated by an integer generating vector and a uniform random shift. Combined with the Brownian bridge construction, these rules are used to approximate $I_m$ for SDEs with constant diffusion. Within the weighted function-space framework of Nichols and Kuo [16], we establish mixed derivative estimates and verify the required summability conditions on the weights. To handle the nested composition structure of the EM integrand, we use the multivariate Faà di Bruno formula [14 Theorem 1] (see also [3, 8]) to expand its mixed derivatives. We organize the resulting terms and their propagation through the EM recursion using a rooted-tree representation, following the standard terminology in [13 Section 4.2]. The decay and localized supports of the BB basis functions then allow us to obtain summability estimates uniform in $m$.

Since the resulting weights are not available in an explicitly computable form, a component-by-component (CBC) construction [18, 16] is impractical for these weights. To avoid this difficulty, we choose generating vectors at random and use the median of independent lattice estimates [7, 28]. Specifically, we independently generate $R = 2k - 1$ lattice estimates, each using a uniformly random generating vector and an independent uniform shift, and take their median. This approach requires no explicit knowledge of the weights. Taking the median also improves stability by reducing the influence of unusually large errors and increasing the probability of obtaining an accurate estimate.

Our analysis of the randomly shifted lattice rule and its median estimator yields the following main results on sampling error.

• For an $N$-point randomly shifted lattice rule with a uniformly random generating vector, we prove a mean absolute sampling error bound of $O(N^{-1+\varepsilon})$ for every $\varepsilon\in(0,1/2)$ and every prime $N$. The implied constant is independent of the number of time steps $m$. See Theorem 3.1.

• The median of $R = 2k - 1$ independent lattice estimates achieves an $L^q$ sampling error of $O(N^{-1+\varepsilon})$ for every $1 \le q < k/(1-\varepsilon)$, with an implied constant independent of $m$. In particular, the median of three estimates achieves the same rate in root mean square error. See Theorem 3.1.

• Under an additional summability condition on the columns of the diffusion matrix and suitable uniform bounds on the model data, the implied constants in the error bounds for both estimators are further independent of the state and Brownian path dimensions. See Corollary 3.3.

The remaining time discretization error is covered by the classical strong and weak convergence theory for Euler–Maruyama approximations; see, for example, [10 Chapters 10 and 14] and [9 Section 5]. We retain a general assumption on the Euler–Maruyama time discretization error and verify it under two sets of time regularity conditions. Theorems 3.6 and 3.7 establish bias rates of order $m^{-1/2}$ and $m^{-1}$, respectively. The sampling error estimates and the Euler–Maruyama time discretization error bound together give the total error bounds for both estimators in Corollary 3.8.

The remainder of the paper is organized as follows. § 2 states the assumptions on the SDE (1.1) and payoff $g$ and describes the Brownian bridge construction. § 3 states the sampling and total error results. § 4 establishes the mixed derivative summability result, and § 5 derives the lattice sampling error bounds. § 6 presents the numerical experiments. § 7 concludes the paper. Technical proofs are deferred to the appendices.

Notation and preliminaries

This section states the assumptions used in our analysis and reviews the Brownian bridge construction, the Euler–Maruyama scheme, the weighted function spaces of Nichols and Kuo [16], and randomly shifted lattice rules. We also introduce the median lattice estimator following [7, 28]. Before proceeding, we fix some notation.

We write $\mathbb{Z}$ for the integers, $\mathbb{N} = \{1,2,\ldots\}$ for the positive integers, and $[s] = \{1,\ldots,s\}$ for $s \in\mathbb{N}$. For a finite set $u$, we denote its cardinality by $|u|$ and the set of its partitions into nonempty blocks by $\mathcal{P}(u)$, with $\mathcal{P}(\emptyset) = \{\emptyset\}$. For $\pi\in\mathcal{P}(u)$, define the set of labels occurring as singleton blocks of $\pi$ by

For an interval $I \subseteq\mathbb{R}$, $|I|$ denotes its length. Throughout the paper, $d_x$ and $d_w$ denote the state and Brownian path dimensions, respectively. An unnamed $C$ denotes a positive constant that may change from line to line. Named constants retain the definitions and parameter dependence specified when they are introduced. We write $\phi(z) = (2\pi)^{-1/2}e^{-z^2/2}$ and $\Phi$ for the standard Gaussian density and distribution function, respectively, and $\phi_s(\boldsymbol{z}) = \prod_{j=1}^{s}\phi(z_j)$ for the standard Gaussian density on $\mathbb{R}^s$. The inverse $\Phi^{-1}$ is applied componentwise to vector arguments.

Vectors are denoted by bold symbols and their scalar components by unbolded symbols. The $r$th standard basis vector is denoted by $\boldsymbol{e}_r$, with its dimension determined by the context. For $d \in\mathbb{N}$, $\mathrm{Id}_d$ denotes the $d \times d$ identity matrix. For $\boldsymbol{z} = (z_1,\ldots,z_s)$ and $u \subseteq[s]$, we write $\boldsymbol{z}_u = (z_j)_{j\in u}$ and $\boldsymbol{z}_{-u} = (z_j)_{j\notin u}$. For a sufficiently smooth scalar-valued function $f : \mathbb{R}^s \to\mathbb{R}$, we define

Thus $\partial_u$ differentiates once with respect to each coordinate indexed by $u$.

For $\boldsymbol{z} \in\mathbb{R}^s$ and a matrix $A = (A_{ij}) \in\mathbb{R}^{r\times w}$, we use the Euclidean vector norm, the induced matrix 2-norm, and the Hilbert–Schmidt matrix norm: