Sampling Error Versus Time-Discretization Bias

Two Different Sources of Error

A Monte Carlo estimate based on simulated SDE paths can miss its target for two separate reasons. Even if the simulated quantity had the correct mean, a finite number of independent paths would produce sampling error. Separately, a time-stepping scheme replaces the continuous-time process by a discretized one, whose expected output may be biased. Increasing the number of paths reduces the first error but does not, by itself, remove the second.

Work on a filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_t),\mathbb{P})$ carrying a Brownian motion $B_t$. For an SDE whose solution is represented using an Ito integral, the integrand is taken to be adapted and progressively measurable, with $\mathbb{E}\int_0^T H_s^2\,ds<\infty$ on each finite time interval; this is the class $\mathcal{L}^2$. Suppose the target is $\theta=\mathbb{E}[F(X_T)]$, where $F(X_T)$ has a finite second moment. A numerical scheme with time step $h$ produces an output $Y_h$, also with a finite second moment. For example, $Y_h$ may be $F(X_T^h)$, where $X_T^h$ is the numerical approximation at terminal time.

Definition: Write $m_h=\mathbb{E}[Y_h]$, $v_h=\operatorname{Var}(Y_h)$, and $b_h=m_h-\theta$. The time-discretization bias is $b_h$. For $N$ independent simulated paths, let $\overline{Y}_{h,N}=N^{-1}\sum_{j=1}^N Y_h^{(j)}$, where the $Y_h^{(j)}$ are independent copies of $Y_h$. The mean-square error relative to the continuous-time target is $\mathbb{E}[(\overline{Y}_{h,N}-\theta)^2]$.

The Exact Mean-Square Error Decomposition

Theorem: If $Y_h$ has a finite second moment and the $N$ simulated outputs are independent and identically distributed, then $$ \mathbb{E}\big[(\overline{Y}_{h,N}-\theta)^2\big] =\frac{v_h}{N}+b_h^2. $$

Proof. The estimator has mean $\mathbb{E}[\overline{Y}_{h,N}]=m_h$. Independence gives $\operatorname{Var}(\overline{Y}_{h,N})=v_h/N$: in the variance of the sum, each of the $N$ terms contributes $v_h$, and every cross-covariance is zero. For any square-integrable random variable $U$, expanding $U-\theta=(U-\mathbb{E}U)+(\mathbb{E}U-\theta)$ and taking expectations shows that $\mathbb{E}[(U-\theta)^2]=\operatorname{Var}(U)+(\mathbb{E}U-\theta)^2$; the cross term vanishes because $\mathbb{E}[U-\mathbb{E}U]=0$. Apply this identity to $U=\overline{Y}_{h,N}$, whose mean is $m_h$, to obtain the formula. $\square$

The two terms have different behavior. At fixed $h$, the sampling variance $v_h/N$ decreases in proportion to $1/N$, so the standard error is $\sqrt{v_h/N}$. The bias $b_h$ is determined by the discretization and remains unchanged when more paths are simulated. Thus a standard error or confidence interval for $m_h$ describes uncertainty from sampling around the discretized mean; it does not automatically account for the distance from $m_h$ to $\theta$.

Worked Example: Euler Sampling and Bias for an Ornstein–Uhlenbeck Mean

Consider $dX_t=-X_t\,dt+dB_t$, $X_0=1$, on $[0,1]$, and estimate $\theta=\mathbb{E}[X_1]$ using Euler–Maruyama with four equal steps. The drift coefficient $-x$ is globally Lipschitz and has linear growth, and the diffusion coefficient $1$ is bounded; the usual strong solution exists. Its integral form is $X_t=1-\int_0^t X_s\,ds+\int_0^t 1\,dB_s$, with the adapted integrand $1\in\mathcal{L}^2$. The exact mean is $\theta=e^{-1}$.

Here $h=1/4$, and the Euler recursion is $X_{k+1}=(3/4)X_k+\Delta B_k$, where each $\Delta B_k$ has mean zero and variance $1/4$. Repeatedly taking expectations gives $m_h=(3/4)^4=81/256=0.3164$ to four significant digits. To compute the output variance, each independent Brownian increment contributes its variance multiplied by the square of its later recursion coefficient:

v_h=\frac14\left(1+\frac{9}{16}+\frac{81}{256}+\frac{729}{4096}\right) =\frac14\left(\frac{4096+2304+1296+729}{4096}\right) =\frac{8425}{16384}.

Thus $v_h=0.5142$ to four significant digits. The bias is $b_h=81/256-e^{-1}$, which is approximately $0.31640625-0.36787944=-0.05147$; its square is approximately $0.002649$. If $N=81$, the sampling variance is $(8425/16384)/81$, approximately $0.006348$. The exact decomposition therefore gives mean-square error approximately $0.006348+0.002649=0.008998$.

Check: the decomposition predicts a root mean-square error of approximately $\sqrt{0.008998}=0.09486$, larger than the sampling standard error $\sqrt{0.006348}=0.07968$, because it also includes bias. If the requested sampling standard error is at most $0.08=2/25$, then $N\ge v_h/(2/25)^2=(8425/16384)/(4/625)$, approximately $80.35$, so $N=81$ is the smallest integer meeting that sampling-only requirement. It does not eliminate the bias.

Using Weak-Order Information

For an output whose expected value converges with weak order $p>0$, a common asymptotic assumption is $m_h=\theta+c h^p+o(h^p)$ as $h$ tends to zero, for a constant $c$. This assumption concerns expectations, not pathwise error. It is useful for interpreting bias and comparing grid resolutions, but it need not describe a coarse grid accurately; it also allows $c=0$, in which case the leading bias term vanishes.

Proposition: Under the displayed weak-error expansion, the leading bias on the finer grid $h/2$ is estimated by $(m_h-m_{h/2})/(2^p-1)$, up to an $o(h^p)$ remainder.

Proof. The expansion at the two step sizes gives $m_h=\theta+c h^p+o(h^p)$ and $m_{h/2}=\theta+c(h/2)^p+o(h^p)$. Subtracting in the stated order yields $$ m_h-m_{h/2}=c h^p(1-2^{-p})+o(h^p). $$ Since $1-2^{-p}=(2^p-1)/2^p$, division by $2^p-1$ gives $c(h/2)^p+o(h^p)$, which is the leading bias $m_{h/2}-\theta$ on the finer grid. If $c=0$, the same calculation says that the displayed estimate is only $o(h^p)$; it does not identify a nonzero leading term. $\square$

The proposition concerns the means $m_h$, not single simulated outputs. In a numerical experiment those means are usually estimated from samples, so the difference between two sample averages also contains sampling noise. The formula is therefore a bias diagnostic based on an asymptotic model, not an exact replacement for the unknown bias at every step size.

Worked Example: Checking the Sign of a Two-Grid Bias Estimate

Suppose $m_h=\theta+2h$ exactly, with $\theta=5$, and take $h=1/10$. Then $m_h=5+2(1/10)=5.2$, while $m_{h/2}=5+2(1/20)=5.1$. For $p=1$, the finer-grid bias is $m_{h/2}-\theta=5.1-5=0.1$. The two-grid estimate gives $(m_h-m_{h/2})/(2^1-1)=(5.2-5.1)/1=0.1$.

Check: the sign is positive, as is the actual finer-grid bias. Reversing the difference would give $(5.1-5.2)/1=-0.1$, the wrong sign. The identity also checks exactly for any positive $h$ here: $m_h-m_{h/2}=2h-2(h/2)=h$, and division by $2^1-1=1$ gives $h$, which equals $2(h/2)$, the fine-grid bias.

A Practical Error Budget

Suppose the desired root-mean-square error is at most $\varepsilon>0$. Since the two terms in the exact decomposition are nonnegative, a practical plan is to assign a bias allowance $\varepsilon_b$ and a sampling allowance $\varepsilon_s$, chosen so that $\varepsilon_b^2+\varepsilon_s^2\le\varepsilon^2$. If $|b_h|\le\varepsilon_b$, then requiring $v_h/N\le\varepsilon_s^2$ ensures the total mean-square error is at most $\varepsilon^2$. When $\varepsilon_s>0$, this sampling requirement is met by any integer $N\ge v_h/\varepsilon_s^2$, so the least such positive integer is the ceiling of that ratio (or $1$ if the ratio is below $1$).

There is also an exact feasibility check if the bias and variance are known. For $v_h>0$, a finite positive sample size can achieve mean-square error at most $\varepsilon^2$ if and only if $|b_h|<\varepsilon$; in that case one needs $N\ge v_h/(\varepsilon^2-b_h^2)$. Indeed, rearranging the decomposition gives precisely $v_h/N\le\varepsilon^2-b_h^2$. If $|b_h|\ge\varepsilon$, the left side $v_h/N$ is positive for every finite $N$, so the target cannot be met. In the degenerate case $v_h=0$, the output is constant almost surely and the mean-square error is $b_h^2$ for every $N$; the target is met exactly when $|b_h|\le\varepsilon$.

Worked Example: Splitting a Total Tolerance

For a given step size, suppose the output has mean $m_h=2+3h$ and variance $v_h=1+h$, with target $\theta=2$. These moment formulas can, for example, define a simulation model $Y_h=2+3h+\sqrt{1+h}\,Z$, where $Z$ is standard normal. Take $h=1/10$ and total tolerance $\varepsilon=1/2$. The bias is $b_h=3(1/10)=3/10$, and the variance is $v_h=1+1/10=11/10$.

Use the full bias allowance $\varepsilon_b=3/10$. The remaining sampling allowance is $\varepsilon_s=\sqrt{(1/2)^2-(3/10)^2}=\sqrt{1/4-9/100}=\sqrt{16/100}=2/5$. Therefore $N$ must satisfy $(11/10)/N\le(2/5)^2=4/25$, or $N\ge(11/10)(25/4)=55/8=6.875$. Choose $N=7$. The resulting mean-square error is $9/100+(11/10)/7=9/100+11/70=0.2471$ to four significant digits, which is below $\varepsilon^2=1/4=0.25$.

Check: with $N=6$, the mean-square error would be $9/100+11/60=0.2733$, which exceeds $0.25$. Thus $N=7$ is the smallest integer satisfying this chosen error budget.

Common Mistake: More Paths Do Not Remove Bias

A frequent error is to report $v_h/N$ as the mean-square error relative to the continuous-time target. That quantity is only the sampling-variance term. In the Ornstein–Uhlenbeck example, it is approximately $0.006348$ for $N=81$, but the squared bias is approximately $0.002649$; the total is approximately $0.008998$, not $0.006348$. The calculation $v_h/N$ is appropriate for the variance of the sample average, or for its mean-square error relative to $m_h$, but not for its mean-square error relative to $\theta$ unless the bias is zero.

Refining the time step generally reduces discretization bias when the scheme has a positive weak order, but it also increases the number of time steps per path and may change $v_h$ and computational cost. Increasing $N$ reduces sampling error at the chosen step size but does not improve the discretized mean. A sound experiment therefore checks both components at a common $h$, states whether its interval is centered on $m_h$ or intended to approximate $\theta$, and reports the assumptions used to estimate the bias.

Check Your Understanding

Use the exact decomposition and distinguish the discretized mean from the continuous-time target.

If $v_h=2$, $b_h=1/5$, and $N=100$, compute the mean-square error relative to $\theta$.
At fixed $h$, what happens to the sampling standard error and the bias when $N$ is multiplied by four?
For $m_h=\theta+c h^2+o(h^2)$, write the two-grid estimate of the leading bias on the $h/2$ grid.
If $v_h>0$ and $|b_h|=\varepsilon$, can any finite $N$ attain mean-square error at most $\varepsilon^2$? Explain using the exact formula.
Why does a confidence interval centered at the sample average of $Y_h$ not, without further adjustment, quantify the time-discretization bias?