A control variate method based on polynomial approximation of Brownian path

TL;DR

This work tackles efficiently estimating expectations of SDE solutions by introducing a control variate built from a coarse, parabolic Brownian motion discretization that is conditioned on the fine Brownian increments used in the primary simulation. The main innovation is strong coupling between the fine-path estimator and the coarse-path control variate, allowing a budget-driven optimization that yields a quadratic error decaying as ${\cal E} \lesssim {\cal C}^{-\frac{4\alpha\gamma+2\alpha}{4\alpha\gamma+2\alpha+1}}$, outperforming standard MC rates (e.g., ${\cal E} \lesssim {\cal C}^{-2/3}$ for Euler). The method is demonstrated across Euler and RK discretizations, including multiplicative-noise SDEs and a double-well potential, with explicit asymptotic rates such as ${\cal E} \lesssim {\cal C}^{-6/7}$ for Euler+parabolic BM and ${\cal E} \lesssim {\cal C}^{-12/13}$ for SRA1, confirming substantial variance reduction. A discussion of a multilevel Monte Carlo integration shows how the CV approach can be embedded in MLMC to further reduce computational cost in many-query settings, highlighting practical impact for finance and physics applications where SDE expectations arise frequently.

Abstract

We present a novel control variate technique for enhancing the efficiency of Monte Carlo (MC) estimation of expectations involving solutions to stochastic differential equations (SDEs). Our method integrates a primary fine-time-step discretization of the SDE with a control variate derived from a secondary coarse-time-step discretization driven by a piecewise parabolic approximation of Brownian motion. This approximation is conditioned on the same fine-scale Brownian increments, enabling strong coupling between the estimators. The expectation of the control variate is computed via an independent MC simulation using the coarse approximation. We characterize the minimized quadratic error decay as a function of the computational budget and the weak and strong orders of the primary and secondary discretization schemes. We demonstrate the method's effectiveness through numerical experiments on representative SDEs.

Figures (5)

• Figure 1: Empirical quadratic error $\mathcal{E}$\ref{['eq:empirical_quadratic_error']} ($M" = 1000$) versus cost $\mathcal{C}$ for the MC estimation of a QoI associated with the SDE \ref{['eq:sde1']}. The fine discretisation is based on the Euler-Maruyama method ($\alpha=1$). The parabola ODE method ($\gamma=1$) is used for the CV method. Left: Comparison between the optimized standard MC expressed in \ref{['eq:st_MC']} - \ref{['eq:st_optimized_parameters']} - \ref{['eq:st_optimized_error']} and the optimized CV MC expressed in \ref{['eq:defIMMp']} - \ref{['eq:cv_optimized_parameters']} - \ref{['eq:cv_optimized_error']}. The steepest slope $-6/7$ corresponds to the CV strategy and the other $-2/3$ to the standard optimized MC estimation. Right: Comparison between different time step sizes and numbers of MC samples in the CV strategy where $h \sim \mathcal{C}^{-x}$, $h' \sim \mathcal{C}^{-y}$, $M \sim \mathcal{C}^{1-x}$, and $M' \sim \mathcal{C}^{1-y}$. We examine six configurations. S1) $x = \frac{1}{7}$, $y = \frac{2}{7}$ (blue squares); S2) $x = \frac{1}{7}$, $y = \frac{3}{7}$ (red squares); S3) $x = \frac{1}{7}$, $y = \frac{4}{7}$ (black squares); C1) $x = \frac{2}{7}$, $y = \frac{3}{7}$ (red circles); C2) $x = \frac{2}{7}$, $y = \frac{4}{7}$ (black circles); and T1) $x = \frac{3}{7}$, $y = \frac{4}{7}$ (black triangles). The corresponding slopes are shown in the legend. These are numerical values obtained by log regression of the empirical quadratic errors. The steepest slope is observed for $x = \frac{1}{7}$, $y = \frac{3}{7}$, which is consistent with the theory.
• Figure 2: Empirical quadratic error $\mathcal{E}$ versus cost $\mathcal{C}$ for the MC estimation of a QoI related to the SDE \ref{['eq:sde2']}. Same configuration as Figure \ref{['fig:BSIGMAGENERAL1']}. The steepest slope is again observed for $x=1/7$, $y=3/7$, which is consistent with the theory.
• Figure 3: Empirical quadratic error $\mathcal{E}$\ref{['eq:empirical_quadratic_error']} versus cost $\mathcal{C}$ for the MC estimation of a QoI associated with the SDE \ref{['eq:sde3']}. The fine discretization is based on the stochastic RK method (SRA1) ($\alpha=2$). The parabola ODE method ($\gamma=1$) is used for the CV method. Left: Comparison between the optimized standard MC expressed in \ref{['eq:st_MC']} - \ref{['eq:st_optimized_parameters']} - \ref{['eq:st_optimized_error']} and the optimized CV MC expressed in \ref{['eq:defIMMp']} - \ref{['eq:cv_optimized_parameters']} - \ref{['eq:cv_optimized_error']}. The steepest slope $-12/13$ corresponds to the CV strategy and the other $-4/5$ to the standard optimized MC estimation. Right: Comparison between different time step sizes and number of MC samples in the CV strategy where $h \sim \mathcal{C}^{-x}$, $h' \sim \mathcal{C}^{-y}$, $M \sim \mathcal{C}^{1-x}$, and $M' \sim \mathcal{C}^{1-y}$. We examine six configurations. S1) $x = \frac{1}{13}$, $y = \frac{2}{13}$ (blue squares); S2) $x = \frac{1}{13}$, $y = \frac{3}{13}$ (red squares); S3) $x = \frac{1}{13}$, $y = \frac{4}{13}$ (black squares); C1) $x = \frac{2}{13}$, $y = \frac{3}{13}$ (red circles); C2) $x = \frac{2}{13}$, $y = \frac{4}{13}$ (black circles); and T1) $x = \frac{3}{13}$, $y = \frac{4}{13}$ (black triangles). The corresponding slopes are shown in the legend. The steepest slope is observed for $x = \frac{1}{13}$, $y = \frac{3}{13}$, which is consistent with the theory.
• Figure 4: Strong convergence error $E(h) = \mathbb{E} [ | \check{X}_1^{h} - {X}_1|^2 ]^{1/2}$ for the parabola method and the IGBM case. We numerically find that the strong error is $\epsilon h^\gamma$ with $\gamma=1$ (as predicted by the general theory) and $\epsilon=\exp(-4)$.
• Figure 5: This figure illustrates the comparison between the minimized theoretical errors and the empirical errors observed in simulation. Left:$x^\star$ below $0.2$ and $y^\star$ above $0.4$ versus $\log(\mathcal{C})$. These parameters minimize the right-hand-side of \ref{['epsE']}, where $h,h',M$ and $M'$ take the form $h \sim \mathcal{C}^{-x^\star}, h' \sim \mathcal{C}^{-y^\star}, M \sim \mathcal{C}^{1-x^\star}$ and $M' \sim \mathcal{C}^{1-y^\star}$. Right:$-\log(\mathcal{E}^{\rm th})/\log(\mathcal{C})$ where $\mathcal{E}^{\rm th}$ is the minimized theoretical quadratic error involving the parameter $\epsilon$. The black dots linked by straight lines correspond to $-\log(\mathcal{E}^{\rm emp})/\log(\mathcal{C})$ versus $\log(\mathcal{C})$ where the empirical error $\mathcal{E}^{\rm emp}$ similar to \ref{['eq:empirical_quadratic_error']} is estimated with $h \sim \mathcal{C}^{-x^\star}, h' \sim \mathcal{C}^{-y^\star}, M \sim \mathcal{C}^{1-x^\star}$ and $M' \sim \mathcal{C}^{1-y^\star}$. Due to the computational time limitation, we only show $\log (\mathcal{C}) \in [15,20]$. In both subfigures, the range of values for $\log(\mathcal{C})$ is $[10, 100]$. The curves marked with circles corresponds to $\epsilon = 10^{-2}$, while the curves marked with squares corresponds to $\epsilon = 1$. The solid line is the asymptotic limit for $\mathcal{C}$ large.

Theorems & Definitions (7)

• Remark 1
• Remark 2
• Proposition 3
• Remark 4
• Remark 5
• Remark 6
• Lemma 7