Parabolic Approximation of Brownian Motion

• Piecewise parabolic approximation of Brownian motion is a method that uses quadratic polynomials to match both endpoints and the time-integral, enhancing convergence over linear schemes.
• The method employs Gaussian increments and area variables to construct continuous quadratic splines, enabling efficient simulation and tight coupling with finer discretizations.
• Applications include variance reduction in SDE simulations and control variate construction, yielding improved error constants and computational efficiency.

Piecewise parabolic approximation of Brownian motion refers to the representation or simulation of Brownian paths using quadratic (parabolic) polynomials on subintervals of a time discretization. This technique arises out of both optimal polynomial expansions for Brownian motion and practical implementations for efficient stochastic simulation, strong/weak approximation, control variate construction, and ε-strong simulation. The essential feature is that the parabolic interpolant matches not only endpoints but also the time-integral (area) of the Brownian path over each mesh interval, yielding improved convergence and error properties compared to piecewise linear schemes.

1. Mathematical Construction and Characterization

Given a standard Brownian motion $W_t$ on $[0,T]$ and a discretization $0=t_0 < t_1 < \dots < t_N=T$ with mesh size $h = T/N$, the piecewise parabolic approximation replaces each segment of the path by a unique quadratic interpolant $\widetilde{W}_t$ on $[t_k, t_{k+1}]$ that satisfies:

• $\widetilde{W}_{t_k} = W_{t_k}$,
• $\widetilde{W}_{t_{k+1}} = W_{t_{k+1}}$,
• $$\displaystyle \int_{t_k}^{t_{k+1}} \widetilde{W}_u\,du = \int_{t_k}^{t_{k+1}} W_u\,du$$.

Defining $\Delta W_k := W_{t_{k+1}} - W_{t_k}$ and the area variable $$H_k := -\frac{1}{2}\Delta W_k + \frac{1}{h}\int_{t_k}^{t_{k+1}}(W_u-W_{t_k}) du$$, one obtains $H_k\sim N(0,h/12)$ independent of $\Delta W_k$. For $t\in[t_k, t_{k+1}]$, setting $\theta=(t-t_k)/h$ yields

$$\widetilde{W}_t = W_{t_k} + \theta\,\Delta W_k + 6\theta(1-\theta)H_k.$$

This construction is conditional-Gaussian: the pair $(\Delta W_k, H_k)$ encodes both the martingale increment and area statistics and supports tight coupling to finer discretizations via the Gaussian conditioning theorem (Garnier et al., 11 Nov 2025, Foster et al., 2019). On adaptive or uniform grids, the same structure applies, and for finer-scale representations, parabolic coefficients can be determined directly from collections of normal increments.

2. Polynomial Expansions and Optimal Approximations

Orthogonal polynomial expansions of Brownian motion use eigenfunctions of the covariance operator for the Brownian bridge, leading to global representations in terms of Jacobi polynomials $P_n^{(-1,-1)}$ and weighted $L^2$ error minimization by Mercer’s theorem (Foster et al., 2019):

• The Karhunen–Loève expansion expresses $W(t) = W_1 t + \sum_{k=1}^\infty I_k e_k(t)$, where $e_k$ is a polynomial eigenfunction and $I_k$ are independent Gaussians with variance $\lambda_k = 1/[k(k+1)]$.
• Truncating at $N$ terms yields an optimal global polynomial approximant, with weighted mean-square error decaying as $O(N^{-1/2})$.

Locally, parabolic approximation on $[t_i, t_{i+1}]$ can be built from Brownian increments and rescaled area variables, yielding a quadratic spline path $\widehat{W}$ that matches endpoint values and time-integrals (Foster et al., 2019). This paradigm admits both global and local error analyses, with strong (pointwise) error $O(h)$ and global supremum-norm error $O(\sqrt{h\log(1/h)})$.

3. Algorithmic Implementation and Coupling

The quadratic spline construction is fully explicit and requires only i.i.d. Gaussian samples:

• For $M$ intervals on $[0,T]$:
  • Sample $\Delta W_i \sim N(0,h)$ and $H_i \sim N(0,h/12)$.
  • Compute $W(t_{i+1}) = W(t_i) + \Delta W_i$.
  • Define $\wideparen{W}^{(i)}(u) = W(t_i) + \frac{u-t_i}{h}\Delta W_i + \frac{6(u-t_i)(t_{i+1}-u)}{h^2} H_i$ for $u\in[t_i, t_{i+1}]$.
• The total number of independent Gaussian variables used is $2M$, and the output is a continuous quadratic spline.

In control variate or multilevel Monte Carlo contexts (Garnier et al., 11 Nov 2025), primary fine-scale discretizations (using Brownian increments $g_i'$) can be efficiently coupled to secondary coarse parabolic approximations. Conditioned on fine increments, the pair $(\Delta W_k, H_k)$ used in the parabolic interpolant is Gaussian with the correct joint law, enabling strong coupling. Computational complexity per coarse step is $O(q)$ (sum over fine increments), reducible to $O(1)$ with streaming updates. Storage requirements are $O(N)$ or single-interval in constant-memory approaches.

4. Error Estimates and Convergence Properties

Quantitative error estimates demonstrate the advantages over piecewise linear (Wong–Zakai) approximation:

Approximation Type	Strong Order	Weak Order	Pointwise MSE	Area Matching
Piecewise Linear	$1/2$	$1$	$O(h)$	Endpoints only
Piecewise Parabolic	$1$	$2$	$O(h)$	Endpoints + Area
Global Polynomial (KL)	$-$	$-$	$O(N^{-1/2})$	Weighted $L^2$

• On each interval, sup-norm error is $O(h)$; global supremum scales as $O(\Delta t)$ (Garnier et al., 11 Nov 2025), with $O((\Delta t)^2)$ weak error for smooth functionals.
• In SDE integration (e.g., inhomogeneous geometric Brownian motion), the parabola–ODE method exhibits strong error rate $1.0$ with error constant $\sim7\times$ smaller than linear; log-ODE methods reach strong order $1.5$ and weak order $2.0$ (Foster et al., 2019).

A plausible implication is that parabolic schemes, while not improving the order for pathwise error over linear, achieve substantially tighter error constants and enable higher-order convergence for integrated SDE functionals.

5. ε-Strong Approximation via Heat-Ball Skeletons

An alternative approach to parabolic path approximation leverages ε-strong simulation via Brownian skeletons constructed from heat-ball domains (Deaconu et al., 2020). Given desired accuracy $\epsilon>0$:

• The time-space domain $D_k$ is defined so that the boundary $$\varphi_\epsilon(t) = \sqrt{t\ln(\epsilon^2 e / t)}$$ is compactly supported and $|\varphi_\epsilon(t)| \leq \epsilon$.
• Sequential exit times and positions are determined via sampling $A_k \sim \text{Gamma}(3/2,2)$, $Z_k \in \{\pm1\}$:
  • $U_k^\epsilon = \epsilon^2 e^{1-A_k}$,
  • $T_k = T_{k-1} + U_k^\epsilon$,
  • $X_k = X_{k-1} + Z_k \varphi_\epsilon(U_k^\epsilon)$.
• The piecewise parabolic interpolant is then $$\widehat{B}_t = X_{k-1} + Z_k \varphi_\epsilon(t-T_{k-1})$$ for $t \in [T_{k-1}, T_k]$.

This algorithm guarantees $$\sup_{0 \leq t \leq T} |B_t - \widehat{B}_t| \leq \epsilon$$ almost surely, with average number of steps $N_T^\epsilon \leq (3\sqrt{3}\delta T)/\epsilon^2$ for arbitrary $\delta>1$, and overall computational complexity $O(T/\epsilon^2)$ (Deaconu et al., 2020). These strong error bounds are rigorous and independent of path complexity.

6. Practical Application in SDE Simulation and Monte Carlo

Piecewise parabolic approximation schemes are widely applied in the numerical simulation of SDEs, especially in the context of control variate construction and variance reduction (Garnier et al., 11 Nov 2025), as well as strong/weak convergence analysis (Foster et al., 2019):

• Control variate construction couples fine and coarse discretizations via common Brownian increments, and the parabolic coarse path acts as a high-correlation control variate, yielding minimal bias.
• In SDE integrators, combining the parabolic path with high-order ODE solvers (parabola–ODE, log–ODE) enhances accuracy and efficiency; parabola–ODE can be 5–10 times faster than piecewise linear methods.
• For applications requiring strict uniform error control (e.g., rare-event simulation, numerical analysis with guaranteed pathwise bounds), $\epsilon$-strong skeleton methods provide almost sure approximation with explicit sampling and sharp step complexity.

A plausible implication is that practitioners should select parabolic or global polynomial schemes when area (integral) matching or weak order accuracy is paramount, and employ skeleton-based approaches when uniform error thresholds are required.

7. Advantages, Limitations, and Recommendations

Piecewise parabolic approximation confers several advantages:

• Accuracy: Strong order one for sup-norm errors, weak order two for expectations of smooth functionals, and tighter constants than linear schemes (Garnier et al., 11 Nov 2025, Foster et al., 2019).
• Explicitness: All constructions use elementary Gaussian sampling; polynomial expansions are optimally weighted in Mercer’s sense.
• Coupling: Conditioning on fine increments yields nearly perfect correlation for control variate strategies.
• Computational Cost: Minor overhead relative to SDE drift/diffusion evaluation when used for variance reduction; storage can be $O(N)$ or constant-memory if past intervals are not revisited.

Drawbacks include the need to code area-matching and parabola construction, and the fact that strong order improvement over linear interpolation is only in constants, not in order, unless further combined with advanced solvers (e.g., log–ODE). Linear schemes suffice when only endpoint matching is required, but suffer from lower accuracy and variance in coupled estimators.

For SDE applications, use the parabola–ODE or log–ODE method for optimal convergence. For path simulation with guaranteed uniform error, apply the heat-ball skeleton algorithm. For high-efficiency Monte Carlo, favor parabolic control variates and couple to fine-scale discretizations via Gaussian conditioning.

---

Piecewise parabolic approximation unifies optimal polynomial expansion, spline-based simulation, control variate strategies, and strong simulation in a systematic framework, providing explicit, high-accuracy representations for Brownian motion and associated diffusions. The technique is foundational for advanced Monte Carlo, multilevel, and high-order SDE integration methodologies.