Hybrid Scheme for Rough Process Simulation

• The paper demonstrates that the hybrid scheme accurately approximates the kernel’s singularity to reduce bias in roughness estimation.
• It partitions simulation into a singular region governed by a power law and a regular region using step-function approximations, balancing precision and computational speed.
• Monte Carlo analyses confirm the method’s optimal convergence rates with improved performance in rough volatility and spatial VMMA models.

A hybrid scheme, in the sense introduced by Bennedsen, Lunde, and Pakkanen, refers to a simulation methodology for stochastic integral representations of rough processes—specifically Brownian semistationary (BSS) processes and spatial volatility modulated moving average (VMMA) fields—where the kernel’s singular local behavior is approximated separately from its regular behavior away from zero. This approach yields a discretization that sharply reduces bias in path roughness estimation and other statistical features, while retaining computational efficiency on par with traditional Riemann-sum methods (Bennedsen et al., 2015, Heinrich et al., 2017).

1. Mathematical Framework and Motivation

The hybrid scheme’s central object is the BSS process

$$X(t)=\int_{-\infty}^t g(t-s)\,\sigma(s)\,dW(s),\quad t\in\mathbb{R},$$

where $W$ is (two-sided) standard Brownian motion and $\sigma$ is a predictable volatility process with $\mathbb{E}[\sigma(t)^2]<\infty$. The kernel $g:(0,\infty)\to[0,\infty)$ is required to be square-integrable and regularly varying at zero: $g(x)=x^\alpha L_g(x)$ for $x\in(0,1]$, with $L_g$ slowly varying and $\alpha\in(-\tfrac12,\tfrac12)\setminus\{0\}$.

The local singularity of $g$ at zero governs the process’s roughness; specifically, the variogram behaves as $V_X(h)\sim C(\alpha) h^{2\alpha+1}L_g(h)^2$ as $h\downarrow0$, with $C(\alpha)$ described explicitly in terms of $\sigma$’s variance and an integral over the kernel power law. Standard Riemann-sum approaches to discretizing the stochastic integral systematically underestimate process roughness, as the step-function approximation fails to capture the local singularity.

2. Hybrid Discretization Principles

The hybrid scheme partitions the simulation into two regions:

• Singular region ($k=1,\dots,\kappa$): The kernel is approximated by its dominant power law form, preserving the singularity:

$$g(t-s)\approx (t-s)^{\alpha}L_g\left(\frac{k}{n}\right),\quad s\in[t-kh,\,t-(k-1)h].$$

• Regular (tail) region ($k\ge\kappa+1$): Away from zero, a step-function approximation suffices:

$$g(t-s)\approx g\left(\frac{b_k}{n}\right),\quad s\in[t-kh,\,t-(k-1)h],$$

with $b_k\in[k-1,k]$ chosen, for instance, to minimize mean square error.

This leads to an overall discretization

$X_n(t) = \check X_n(t) + \hat X_n(t),$

where the singular part $\check X_n(t)$ involves direct simulation of Wiener integrals of the power function at small lags, and the tail part $\hat X_n(t)$ can be evaluated using FFT-accelerated convolutions (Bennedsen et al., 2015).

3. Error Analysis and Optimality Properties

The hybrid scheme achieves a sharp asymptotic mean square error,

$$\mathbb{E}\left[|X(t)-X_n(t)|^2\right] \sim J(\alpha,\kappa,\{b_k\})\,\mathbb{E}[\sigma(0)^2]\,n^{-(2\alpha+1)} L_g(1/n)^2,$$

for $n\to\infty$, under minimal regularity conditions on the long-range decay of $g$, the truncation (with $N_n \sim n^{1+\gamma}$), and the smoothness of $\sigma$. The constant

$$J(\alpha,\kappa,\{b_k\}) = \sum_{k=\kappa+1}^\infty \int_{k-1}^k [y^\alpha - b_k^\alpha]^2 dy$$

can be minimized by an explicit choice of $b_k$ in each cell. For $\kappa=0$, $b_k=k$, the scheme reduces to the standard Riemann-sum error (Bennedsen et al., 2015).

In the VMMA setting in $d$ dimensions, similar structure obtains with

$$X^n(t) = \sum_{z\in K_\kappa}L(\|z\|/n)\,\sigma(t+z/n)\,W^n_{t+z} +\sum_{z\in\bar K_\kappa} g(z/n)\,\sigma(t+z/n)\,W^n_{t+z},$$

where the hybrid partition is geometrically defined and $W^n_{t+z}$ are appropriately correlated Gaussian integrals (Heinrich et al., 2017). The asymptotic $L^2$-error is

$$n^{2(\alpha+1)} L(1/n)^{-2}\, \mathbb{E}[\bigl(X(t)-X^n(t)\bigr)^2] \to \mathbb{E}[\sigma(0)^2]\,J(\alpha,\kappa)$$

with

$$J(\alpha,\kappa) = \sum_{z\notin K_\kappa} \int_{[0,1]^d} (\|u\|^\alpha - \|z+u\|^\alpha)^2 du.$$

4. Algorithmic Implementation and Complexity

The hybrid scheme consists of:

• Calculating the local (singular) terms ($\kappa$ terms per location), involving direct Wiener integrals, with precomputed covariances (small Cholesky decompositions).
• Simulating the regular part via fast convolution. In one dimension, this is an $O(n\log n)$ operation; for $d$-dimensional VMMAs, $O(n^d \log n)$.
• The overall work per realization is dominated by the FFT convolution, matching the Riemann-sum complexity.

For practical settings, $\kappa=1$ or $2$ suffices to recover roughness and suppress bias, and major speed gains accrue when multiple paths are simulated by vectorizing the FFT and sharing local covariance computations (Heinrich et al., 2017, Bennedsen et al., 2015). Table 1 illustrates computational times, with the hybrid scheme competitive with circulant embedding methods for many replicates:

Method	Single Sim (s)	100 Sims (s)
Hybrid $\kappa=0$	12.6	51.0
Hybrid $\kappa=1$	13.2	61.3
Hybrid $\kappa=2$	14.3	72.6
Hybrid $\kappa=3$	15.3	77.7
Circulant Embedding	0.8	75.6
Riemann-sum	1.2	32.5

The fixed overhead for local covariance precomputation is amortized in batch settings.

5. Statistical and Applied Performance

Monte Carlo analyses demonstrate the scheme’s efficacy:

• For estimation of the roughness index $\alpha$, bias and standard deviation in the COF-estimator collapse rapidly as $h\to0$, with negligible bias for hybrid paths at $\kappa\ge1$ versus substantial negative bias for Riemann-sum schemes, especially at $\alpha$ near $-0.5$ (Bennedsen et al., 2015).
• In rough volatility modeling (e.g., rough Bergomi model), hybrid-scheme-driven simulations of the volatility process reproduce the implied-volatility surface accurately, matching exact Cholesky-sampled smiles even for extreme strikes. Standard Riemann sums mis-level the smile (Bennedsen et al., 2015).
• For spatial VMMAs, the estimated Hausdorff dimension matches the theoretical $d-\alpha$ almost perfectly for hybrid $\kappa\ge1$. Riemann-sum approaches systematically underestimate this roughness parameter (Heinrich et al., 2017).

6. Practical Guidelines, Variants, and Related Schemes

Best practices from empirical and theoretical investigation include:

• $\kappa=1$ nearly always sufficient, with minimal benefit from larger $\kappa$.
• Choose the lowest $\gamma$ meeting technical restrictions to limit convolution cost.
• The choice of block evaluation points $b_k$ can be optimized for the MSE constant but in practice, midpoints are nearly optimal.
• For higher-dimensional applications ($d>2$), direct implementation proceeds analogously, though local covariance integrals may lack closed forms.
• The hybrid strategy is robust to kernel shapes—anisotropic kernels can be accommodated by local blockwise constant approximations, with an increase in local covariance costs (Heinrich et al., 2017).

Hybrid schemes are entirely time or space domain, in contrast to some earlier Fourier-based approaches, but exhibit similar computational scalability. The Riemann-sum scheme is recovered as the degenerate case with $\kappa=0$.

7. Impact and Extensions

The hybrid simulation methodology provides a crucial tool for generating high-fidelity rough paths and fields in a statistically principled, computationally feasible manner. Its principal advantage is the restoration of the kernel singularity—essential for accurate inference on roughness and for correct pricing in rough volatility models—without requiring full exact simulation. Extensions to spatial fields (VMMAs) have been developed by adapting the hybrid principle to higher-dimensional grids, with similar accuracy guarantees and robust empirical performance (Heinrich et al., 2017).

The hybrid scheme is now integrated into numerical studies of rough volatility and spatially inhomogeneous stochastic models, including joint calibration algorithms for VIX and SPX derivatives via the rough Bergomi model, where it underpins FFT-accelerated simulation in practical calibration pipelines (Jacquier et al., 2017). Its convergence rates, proven and observed, are at or near the theoretical optimum for discretization methods respecting local kernel behavior.