Monte Carlo Standard Errors for Empirical Quantiles

Why Quantile Precision Depends on Local Density

A simulation can produce thousands of outcomes and still estimate a high quantile imprecisely. The key is not just how many outcomes were generated, but how densely they cluster near the quantile being estimated. If outcomes are crowded near the target, a modest change in cumulative probability corresponds to a small change in the quantile. If outcomes are sparse there, the same probability change can correspond to a large change in the reported value.

This tutorial develops a standard error approximation for an empirical quantile based on independent simulated outcomes. It is a plug-in estimate: the target probability and sample size determine the uncertainty on the probability scale, while an estimate of the local density translates that uncertainty into the units of the simulated quantity.

Definition: Let $X_1,\ldots,X_n$ be independent outcomes from a continuous distribution with cumulative distribution function $F_X$ and density $f_X$. For a probability $p$ strictly between 0 and 1, the population $p$-quantile is $q_p=F_X^{-1}(p)$. The empirical distribution function is $\widehat F_n(x)=n^{-1}\sum_{j=1}^n I(X_j\le x)$, and the empirical $p$-quantile is $\widehat q_p=\inf\{x:\widehat F_n(x)\ge p\}$. With this convention, $\widehat q_p$ is the order statistic $X_{(\lceil np\rceil)}$, where the simulated outcomes are sorted from smallest to largest.

The empirical quantile is a data-based estimate of $q_p$, not an exact value of it. Its sampling uncertainty refers to variation across repeated independent simulation runs of the same size from the same distribution. The usual approximation requires that $f_X(q_p)$ exist, be positive, and be reasonably smooth near $q_p$. It is not an exact finite-sample formula.

The Monte Carlo Standard Error Approximation

Under the stated conditions, the standard error of the empirical quantile is approximated by

\operatorname{SE}(\widehat q_p) \approx \frac{\sqrt{p(1-p)/n}}{f_X(q_p)}.

In practice, the density at the unknown population quantile is replaced by an estimated local density, often evaluated near $\widehat q_p$. Thus the usable plug-in estimate is

\widehat{\operatorname{SE}}(\widehat q_p) = \frac{\sqrt{p(1-p)/n}}{\widehat f_X(\widehat q_p)}.

The numerator is the standard deviation of the empirical cumulative probability at the target: for each simulated outcome, the indicator $I(X_j\le q_p)$ is a Bernoulli variable with success probability $p$, so the sample proportion has variance $p(1-p)/n$. The denominator converts a probability-scale fluctuation into a fluctuation in the units of $X$.

Interpretation: The formula applies to independent, identically distributed simulation outcomes from a continuous model, with a positive, well-behaved density at the target quantile. If the density is estimated, the result is a plug-in approximation and does not separately account for uncertainty in that density estimate. It does not directly apply to a quantile at a probability atom, to a point where the density is zero or irregular, or to correlated outcomes treated as if they were independent.

Why the Density Appears: A Short Derivation

The empirical quantile is defined so that the empirical cumulative probability at its location is approximately $p$. For a large sample from a continuous distribution, this gives $\widehat F_n(\widehat q_p)\approx p=F_X(q_p)$. Near $q_p$, a first-order expansion of the population distribution function gives

F_X(\widehat q_p)-F_X(q_p) \approx f_X(q_p)(\widehat q_p-q_p).

The empirical fluctuation at $q_p$, namely $\widehat F_n(q_p)-p$, has approximate variance $p(1-p)/n$, because it is an average of independent Bernoulli indicators. Equating this fluctuation to the corresponding first-order change in cumulative probability yields

\widehat q_p-q_p \approx -\frac{\widehat F_n(q_p)-p}{f_X(q_p)}.

The sign indicates that if the empirical cumulative probability at $q_p$ is too high, the empirical quantile tends to fall below $q_p$. For the standard error, only the variance matters. Dividing the probability-scale standard deviation by $f_X(q_p)$ gives the stated approximation. This argument also explains the sparse-tail effect: a small $f_X(q_p)$ makes the quantile more sensitive to the same empirical probability fluctuation.

Worked Example: Quantile Precision in a Simulated Claim Model

A hypothetical insurer simulates $n=3{,}600$ independent claim amounts, measured in thousands of dollars, to estimate the 90th percentile. The empirical quantile is $42$ thousand dollars, and a local density estimate there is $0.01$ per thousand dollars. Estimate the Monte Carlo standard error.

Method. Use the plug-in empirical-quantile standard error. The target probability is $p=0.90$, so the probability-scale standard error is

\sqrt{\frac{p(1-p)}{n}} = \sqrt{\frac{0.90(0.10)}{3{,}600}} = \sqrt{\frac{0.09}{3{,}600}} = \sqrt{0.000025} = 0.005.

The local density is $0.01$ per thousand dollars. Therefore,

\widehat{\operatorname{SE}}(\widehat q_{0.90}) = \frac{0.005}{0.01\text{ per thousand dollars}} = 0.5\text{ thousand dollars}.

The estimated standard error is 0.5 thousand dollars, or $500. The empirical quantile is estimated as 42 thousand dollars, but this calculation concerns its simulation error; it does not imply that the true population quantile is exactly 42 thousand dollars. The units check is important: dividing a probability-scale quantity by a density measured per thousand dollars produces a result in thousands of dollars.

Check. Squaring the probability-scale standard error gives $0.005^2=0.000025=0.09/3{,}600$, and multiplying $0.5$ thousand dollars by $0.01$ per thousand dollars recovers $0.005$. The common mistake avoided is reporting the numerator alone as the quantile standard error; it is only the cumulative-probability standard error.

Estimating the Local Density

The density is usually unknown, so it must be estimated from observations near the empirical quantile. One simple local estimate counts outcomes in a window centered at the quantile and divides by both the sample size and the window width. If the window is $[\widehat q_p-h,\widehat q_p+h]$, and $m$ simulated outcomes fall in that window, the estimate is

\widehat f_X(\widehat q_p) = \frac{m/n}{2h} = \frac{m}{2hn}.

This is a local frequency per unit of $X$, not merely the fraction of observations in the window. A narrower window focuses more closely on the density at the quantile but may contain few observations and be noisy. A wider window uses more observations but can average over regions where the density differs. Kernel density estimates use weighted nearby observations rather than a single window; whichever method is used, its units and its behavior near the target should be checked.

Worked Example: Estimating Density From a Local Window

A hypothetical simulation generates $n=5{,}000$ independent annual loss outcomes, in thousands of dollars. The estimated 90th percentile is $\widehat q_{0.90}=70$. A window extending 2 thousand dollars on each side, from 68 to 72, contains 100 simulated outcomes. Use this window to estimate the density and then the quantile standard error.

Step 1: Find the window width. The full width is $72-68=4$ thousand dollars. The fraction of outcomes in the window is $100/5{,}000=0.02$.

Step 2: Estimate the density. Divide that fraction by the width:

\widehat f_X(70) = \frac{0.02}{4} = 0.005\text{ per thousand dollars}.

Equivalently, $100/(5{,}000\times4)=0.005$ per thousand dollars.

Step 3: Calculate the probability-scale standard error.

\sqrt{\frac{0.90(0.10)}{5{,}000}} = \sqrt{0.000018} \approx 0.0042426.

Step 4: Convert to the quantile scale.

\widehat{\operatorname{SE}}(\widehat q_{0.90}) = \frac{0.0042426}{0.005} \approx 0.8485\text{ thousand dollars}.

To four decimal places, the estimated standard error is 0.8485 thousand dollars. This is a plug-in result: it treats the window-based density estimate as fixed. The window contains 100 observations, so it supplies some local information, but a different window width could give a different density estimate. The common mistake avoided is dividing by the half-width of 2; the density estimate uses the full interval width of 4.

Check. The estimated density times the reported standard error is approximately $0.005(0.8485)=0.0042425$, which agrees with the probability-scale standard error up to rounding.

What Sparse Tail Data Do to Precision

At a high target probability, only a small fraction of outcomes lie above the quantile. If the distribution also has a low density near that quantile, the nearby order statistics can be widely separated. In the formula, the probability-scale component may be small, yet division by a small density can still produce a large standard error. A large simulation size alone therefore does not guarantee a stable tail quantile.

A density estimate based on only a few nearby outcomes is itself fragile. Moving the window boundary slightly, or changing one simulated outcome, can materially change the estimated local density and hence the estimated standard error. It is useful to report the density method and the amount of local data used, particularly when the target lies in a sparse region. Repeating the calculation with reasonable alternative window widths can reveal sensitivity, although it does not remove the underlying scarcity of tail information.

Worked Example: A Sparse High-Quantile Region

A hypothetical portfolio model produces $n=2{,}500$ independent simulated annual losses, measured in thousands of dollars. For the 99th percentile, a local window centered at the empirical quantile has half-width 5 thousand dollars and contains 18 outcomes. Estimate the density and the plug-in standard error, and interpret the result.

Step 1: Calculate the full window width and local fraction. The full width is $2(5)=10$ thousand dollars, and the fraction of outcomes in the window is $18/2{,}500=0.0072$.

Step 2: Estimate the local density.

\widehat f_X(\widehat q_{0.99}) = \frac{0.0072}{10} = 0.00072\text{ per thousand dollars}.

Step 3: Calculate the probability-scale standard error.

\sqrt{\frac{0.99(0.01)}{2{,}500}} = \sqrt{\frac{0.0099}{2{,}500}} = \sqrt{0.00000396} \approx 0.0019900.

Step 4: Convert to the loss scale.

\widehat{\operatorname{SE}}(\widehat q_{0.99}) = \frac{0.0019900}{0.00072} \approx 2.7639\text{ thousand dollars}.

To four decimal places, the estimated standard error is 2.7639 thousand dollars. The window contains only 18 outcomes, so the density estimate is based on limited local information and should not be treated as highly stable. The large standard error is consistent with the small estimated density: a small probability-scale fluctuation translates into a comparatively large change in the loss quantile.

Check. The estimate $0.00072$ per thousand dollars multiplied by $2.7639$ thousand dollars is approximately $0.0019900$. The method avoids a common tail trap: assuming that a 99th-percentile estimate must be precise simply because 2,500 simulations were run.

Common Mistakes and Scope

Stopping at the probability-scale standard error. The numerator $\sqrt{p(1-p)/n}$ is not yet in the units of $X$. Divide it by the local density.
Using the density in the wrong direction. A smaller density means a larger quantile standard error. The formula divides by $f_X(q_p)$; it does not multiply by it.
Ignoring units. If the density is reported per thousand dollars, the resulting standard error is in thousands of dollars. Convert units only after the calculation.
Treating a noisy density estimate as exact. A plug-in standard error does not include the additional uncertainty from estimating the density. Sparse local data can make this omission important.
Applying the independent-sample formula to dependent output without adjustment. The derivation uses independent outcomes. If simulated output is correlated, its dependence must be addressed rather than treating the number of recorded observations as an independent sample size.
Using the formula at a point where its conditions fail. Discrete outcomes, a probability atom, or a zero or irregular density near the target require different treatment; the smooth positive-density approximation is not justified there.

Key takeaway: For independent simulation outcomes with a smooth, positive density at the target quantile, estimate quantile Monte Carlo error by dividing $\sqrt{p(1-p)/n}$ by the local density. Low density and sparse nearby outcomes signal that a tail quantile may be unstable even when the simulation count appears substantial.

Check Your Understanding

For each problem, use the stated independent-simulation assumptions and round numerical standard errors to four decimal places where needed.

A simulation has $n=6{,}400$ independent outcomes. At $p=0.75$, the estimated local density is $0.012$ per unit. Find the plug-in standard error in units of the simulated outcome. (Answer: $0.4511$ units.)
A hypothetical loss simulation has $n=10{,}000$, targets $p=0.99$, and estimates the local density as $0.0025$ per thousand dollars. Find the estimated standard error. (Answer: $0.3980$ thousand dollars.)
Two independent simulations each use $n=4{,}000$ outcomes and target $p=0.95$. Their estimated local densities are $0.01$ and $0.0025$ per unit. Find both standard errors and identify the less precise estimate. (Answer: $0.3446$ units and $1.3784$ units, respectively; the second is less precise.)
For $n=8{,}000$ outcomes, a window of half-width 1.5 units around the estimated 90th percentile contains 96 outcomes. Estimate the local density and then the quantile standard error. (Answer: density $0.004$ per unit; standard error $0.8385$ units.)
With $p$ and the local density unchanged, a simulation's sample size is increased from $n$ to $9n$. By what factor does the standard error change? (Answer: it becomes one-third as large.)

What You'll Learn