Monte Carlo Primer – When, Why, Sampling Pitfalls

A single-point risk estimate (“our expected annual loss is $2.4M”) invites false precision. When people hear a single number, the number tends to become a budget line. Then the budget line becomes a commitment. But the estimate was uncertain to begin with so the real outcome might be anywhere within a (sometimes very) wide range.

Monte Carlo simulation replaces the single point with a distribution of outcomes. Instead of “$2.1M,” it becomes “$1.1M to $5.8M as the 90% range of outcomes, with a median of $2.1M.” This framing:

• Communicates uncertainty honestly,

• Supports threshold-based decisions (“is there more than a 10% chance we exceed our risk appetite?”), and

• Naturally accommodates the skewed, heavy-tailed distributions that dominate security loss data.

This notebook builds the core Monte Carlo workflow from scratch: compound Poisson loss simulation, convergence diagnostics, risk band extraction, VaR/ES tail measures, baseline-vs-control sensitivity analysis, and reproducibility practices.

Setup

We use decision-security’s Monte Carlo engine (simulate_aggregate_losses), risk band functions, and visualization helpers.

import numpy as np
import matplotlib.pyplot as plt

from decision_security.synth import make_rng, sample
from decision_security.montecarlo import (
    simulate_aggregate_losses,
    make_lognormal_severity,
    risk_bands,
    var_es,
)
from decision_security.viz import plot_risk_bands, plot_loss_distribution

plt.rcParams.update({
    "font.family": "serif",
    "font.size": 10,
    "axes.labelsize": 11,
    "axes.titlesize": 12,
    "xtick.labelsize": 9,
    "ytick.labelsize": 9,
    "legend.fontsize": 9,
    "figure.dpi": 150,
    "axes.spines.top": False,
    "axes.spines.right": False,
})

PRIMARY = "#1A1A1A"
ACCENT = "#E74C3C"
DARK_BG = "#34495E"
LIGHT_GRAY = "#95A5A6"
MED_GRAY = "#7F8C8D"
VERY_LIGHT = "#BDC3C7"

Standard Error \(\sim 1/\sqrt{N}\)

Every Monte Carlo estimate is a sample statistic with a standard error (SE) that shrinks as \(1/\sqrt{N}\). Doubling precision requires quadrupling samples. But the constant in front of \(1/\sqrt{N}\) is the coefficient of variation (CV) of what you are sampling: the relative SE of the mean is \(\text{CV}/\sqrt{N}\). Heavy-tailed severities have large CVs, so “N = 10,000 is enough” is not a universal rule.

The example below samples a lognormal(\(\mu = 10\), \(\sigma = 2\)) severity. Its CV is \(\sqrt{e^{\sigma^2} - 1} \approx 7.3\) – so at \(N = 10{,}000\) the relative SE of the mean is about 7% (\(7.3/\sqrt{10{,}000}\)), not below 1%. The sample estimate below shows less at that \(N\), because the sample SD underestimates the true one until the rare huge draws turn up; once a few of them are in the sample (here by \(N = 100{,}000\)), it overshoots instead. Driving it down to a 1% relative SE would require \(N = (\text{CV}/0.01)^2 \approx 535{,}000\) draws. Tail quantiles (p95, p99) converge slower than the median, because fewer samples populate the tail – but with tails this heavy, the mean can be the noisiest of all.

The practical takeaway: pick \(N\) from the CV of your data and the precision you need, not from a fixed habit. For board-level reporting, run twice with different seeds and check stability directly.

rng = make_rng(42)
sdlog = 2
big = sample("lognormal", 100_000, rng=rng, meanlog=10, sdlog=sdlog)

# The sample CV is driven by the few largest draws, so anchor on the closed form.
cv = np.sqrt(np.exp(sdlog**2) - 1)
sample_cv = big.std(ddof=1) / big.mean()
print(f"CV of lognormal(10, {sdlog}) = sqrt(exp({sdlog}^2) - 1) = {cv:.1f}")
print(f"Sample CV at N=100,000 = {sample_cv:.1f} ({'above' if sample_cv > cv else 'below'} the true value)")
print()
for n in [100, 1_000, 10_000, 100_000]:
    prefix = big[:n]
    se = prefix.std(ddof=1) / np.sqrt(n)
    rel_se = se / prefix.mean()
    print(f"N={n:>7,}  mean={prefix.mean():>14,.2f}  SE={se:>12,.2f}  "
          f"relative SE={rel_se:>6.1%}  theory CV/sqrt(N)={cv / np.sqrt(n):>6.1%}")

target = 0.01
n_needed = (cv / target) ** 2
print(f"\nN for a 1% relative SE = (CV / 0.01)^2 = {n_needed:,.0f}")

CV of lognormal(10, 2) = sqrt(exp(2^2) - 1) = 7.3
Sample CV at N=100,000 = 11.6 (above the true value)

N=    100  mean=     61,573.99  SE=   16,758.39  relative SE= 27.2%  theory CV/sqrt(N)= 73.2%
N=  1,000  mean=    140,299.18  SE=   19,166.36  relative SE= 13.7%  theory CV/sqrt(N)= 23.2%
N= 10,000  mean=    164,858.50  SE=    9,967.10  relative SE=  6.0%  theory CV/sqrt(N)=  7.3%
N=100,000  mean=    171,837.42  SE=    6,284.71  relative SE=  3.7%  theory CV/sqrt(N)=  2.3%

N for a 1% relative SE = (CV / 0.01)^2 = 535,982

rng = make_rng(99)
sizes = np.logspace(1.5, 4.5, 30, dtype=int)
means = []
ses = []
for n in sizes:
    draws = sample("lognormal", n, rng=rng, meanlog=10, sdlog=2)
    means.append(draws.mean())
    ses.append(draws.std(ddof=1) / np.sqrt(n))

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(sizes, means, "o-", markersize=3, color=DARK_BG)
axes[0].axhline(np.exp(10 + 2**2/2), color=ACCENT, linestyle="--", label="True E[X]")
axes[0].set_xscale("log")
axes[0].set_xlabel("Sample Size (N)")
axes[0].set_ylabel("Sample Mean")
axes[0].set_title("Mean Estimate vs Sample Size")
axes[0].legend()

sigma_true = np.exp(10 + 2**2/2) * np.sqrt(np.exp(2**2) - 1)
axes[1].plot(sizes, ses, "o-", markersize=3, color=ACCENT, label="Sample SE")
axes[1].plot(sizes, sigma_true / np.sqrt(sizes), "--", color=PRIMARY, label="Theoretical SE = sigma/sqrt(N)")
axes[1].set_xscale("log")
axes[1].set_yscale("log")
axes[1].set_xlabel("Sample Size (N)")
axes[1].set_ylabel("Standard Error")
axes[1].set_title("SE Shrinks as 1/sqrt(N)")
axes[1].legend()
plt.tight_layout()
plt.show()

As N grows by 10x, the theoretical standard error drops by \(\sqrt{10} \approx 3.2\text{x}\). The left plot shows the sample mean converging toward the true expectation. The right plot shows SE on a log-log scale. The dashed line is the theoretical \(\sigma/\sqrt{N}\), with slope \(-0.5\). The sample SEs scatter widely and mostly sit below it, because the sample SD misses the rare huge draws that drive the true \(\sigma\).

Practical rule of thumb: 10,000 runs gives stable p50 and p90 estimates. For p95 and beyond, you may need 50,000–100,000 runs.

Note: At \(N = 10{,}000\), your p50 estimate is rock-solid, your p90 is decent, and your p99 is still noisy. Always report tail estimates with a caveat about simulation precision.

Compound Poisson Loss Model

A common workhorse of FAIR-style risk quantification is the compound Poisson model: draw a frequency \(n\) from \(\text{Poisson}(\lambda)\), then draw \(n\) independent severity values from a loss distribution (typically lognormal), and sum them. One iteration gives one simulated annual loss. Repeat \(N\) times to build the loss exceedance curve.

In pseudocode:

for each simulation run:
    n_events = Poisson(lambda)
    losses   = [Lognormal(mu, sigma) for _ in range(n_events)]
    annual_loss = sum(losses)

Everything else – controls modeling, scenario comparison, portfolio aggregation – is layered on top of this basic structure.

rng = make_rng(42)
lam = 0.6
sev = make_lognormal_severity(10, 2)
losses = simulate_aggregate_losses(10_000, lam, sev, rng=rng)

print(f"Simulated {len(losses):,} periods")
print(f"Zero-loss periods: {(losses == 0).sum():,} ({(losses == 0).mean():.1%})")

# Closed-form cross-check: for a compound Poisson, E[loss] = lam * E[severity].
mean_sev = np.exp(10 + 2**2 / 2)
print(f"\nClosed-form E[loss] = lam * E[sev] = {lam} * ${mean_sev:,.0f} = ${lam * mean_sev:,.0f}")
print(f"Simulated mean aggregate loss:                 ${losses.mean():,.0f}")

Simulated 10,000 periods
Zero-loss periods: 5,405 (54.0%)

Closed-form E[loss] = lam * E[sev] = 0.6 * $162,755 = $97,653
Simulated mean aggregate loss:                 $100,899

nonzero = losses[losses > 0]
frac_zero = (losses == 0).mean()

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].bar(["No loss", "Loss > 0"], [frac_zero, 1 - frac_zero],
            color=[LIGHT_GRAY, DARK_BG], edgecolor="white")
axes[0].set_ylabel("Fraction of simulated years")
axes[0].set_title(f"{frac_zero:.0%} of years see no loss at all")
axes[0].set_ylim(0, 1)

log_bins = np.logspace(np.log10(nonzero.min()), np.log10(nonzero.max()), 60)
axes[1].hist(nonzero, bins=log_bins, color=DARK_BG, edgecolor="white", alpha=0.7)
for pct, col in [(50, PRIMARY), (90, MED_GRAY), (95, ACCENT)]:
    val = np.percentile(nonzero, pct)
    axes[1].axvline(val, color=col, ls="--", lw=1.5,
                    label=f"p{pct} of losses = ${val:,.0f}")
axes[1].set_xscale("log")
axes[1].set_xlabel("Loss (USD, log scale)")
axes[1].set_ylabel("Frequency")
axes[1].set_title(f"When a loss occurs ({1 - frac_zero:.0%} of years)")
axes[1].legend(fontsize=9)

plt.tight_layout()
plt.show()

Risk Bands

The simulated loss distribution gives you a full menu of quantiles. Each maps to a different decision context:

• p50 (median): The “typical year” outcome. Use for resource planning, not for expected budgeting: the expected loss is the mean, which heavy tails pull well above the median.

• p90: The stress case. Use for setting risk appetite thresholds and insurance retention levels.

• p95: The planning ceiling. Use for capital reserves, worst-case staffing, and board-level risk reporting.

When communicating results, the p50–p90 band is the most useful single summary: “In a typical year we expect $1–3M in losses; in a bad year, up to $7M.”

bands = risk_bands(losses)
for label, value in bands.items():
    print(f"{label}: ${value:>14,.0f}")

frac_zero = (losses == 0).mean()
median = bands["p50"]
print(f"\nMedian (p50) = ${median:,.0f}; {frac_zero:.0%} of years see no loss.")
if median == 0:
    print("A $0 median says nothing about the bad years; the tail measures (VaR/ES) carry the decision.")
else:
    print(f"p95 = {bands['p95'] / median:.1f}x the median; the tail measures (VaR/ES) carry the decision.")

p50: $             0
p90: $       158,928
p95: $       385,922

Median (p50) = $0; 54% of years see no loss.
A $0 median says nothing about the bad years; the tail measures (VaR/ES) carry the decision.

# Horizontal bars read well even when p50 is zero
band_labels = list(bands.keys())
band_values = list(bands.values())

fig, ax = plt.subplots(figsize=(10, 4))
colors = [LIGHT_GRAY if v == 0 else DARK_BG for v in band_values]
bars = ax.barh(band_labels, band_values, color=colors, edgecolor="white")
for bar, val in zip(bars, band_values):
    if val > 0:
        ax.text(val + max(band_values) * 0.02, bar.get_y() + bar.get_height() / 2,
                f"${val:,.0f}", va="center", fontsize=9)
    else:
        ax.text(max(band_values) * 0.02, bar.get_y() + bar.get_height() / 2,
                "$0", va="center", fontsize=9, color=MED_GRAY)
ax.set_xlabel("Annual Aggregate Loss (USD)")
ax.set_title("Loss Distribution Risk Bands")
ax.invert_yaxis()
plt.tight_layout()
plt.show()

Value-at-Risk and Expected Shortfall

Value at Risk (VaR) at confidence level \(\alpha\) is the \(\alpha\)-quantile of the loss distribution: “there is only a \((1-\alpha)\) probability that losses exceed VaR.” It answers “how bad could it get?” but says nothing about how much worse beyond that threshold.

Expected Shortfall (ES), also called Conditional VaR or Tail VaR, is the average loss given that losses exceed the VaR threshold:

\[ES_\alpha = \mathbb{E}[L \mid L > \text{VaR}_\alpha]\]

ES is strictly more informative for heavy-tailed distributions because it captures tail severity, not just tail probability.

Use VaR for threshold-based decisions (SLAs, appetite statements). Use ES when you need to understand the magnitude of tail outcomes – for instance, when sizing cyber insurance limits or evaluating whether a catastrophic scenario is survivable.

Tip: In some regulatory frameworks (Basel’s FRTB for market risk, the Swiss Solvency Test), ES has replaced VaR precisely because VaR ignores what happens beyond the threshold. Solvency II still uses 99.5% VaR.

v, es = var_es(losses, 0.95)
print(f"VaR(95%): ${v:>14,.0f}")
print(f" ES(95%): ${es:>14,.0f}")
print(f"\nES / VaR = {es / v:.1f}x")

VaR(95%): $       385,860
 ES(95%): $     1,513,947

ES / VaR = 3.9x

nonzero_sorted = np.sort(losses[losses > 0])
log_bins = np.logspace(np.log10(nonzero_sorted.min()), np.log10(nonzero_sorted.max()), 60)

fig, ax = plt.subplots(figsize=(10, 5))
ax.hist(nonzero_sorted, bins=log_bins, color=DARK_BG, edgecolor="none", alpha=0.7)

tail = nonzero_sorted[nonzero_sorted >= v]
ax.hist(tail, bins=log_bins[log_bins >= v * 0.9], color=ACCENT, edgecolor="none",
        alpha=0.8, label=f"Tail (>{v:,.0f})")

ax.axvline(v, color=ACCENT, ls="--", lw=1.5, label=f"VaR(95%) = ${v:,.0f}")
ax.axvline(es, color=PRIMARY, ls=":", lw=1.5, label=f"ES(95%) = ${es:,.0f}")
ax.set_xscale("log")
ax.set_xlabel("Aggregate Loss (USD, log scale)")
ax.set_ylabel("Frequency")
ax.set_title("Tail Region: VaR and Expected Shortfall (nonzero losses)")
ax.legend()
plt.tight_layout()
plt.show()

Sensitivity: Baseline vs. Control Scenario

Monte Carlo makes scenario comparison trivial. Run the simulation twice: once with current parameters (baseline) and once with modified parameters. The difference in output distributions is the control’s risk reduction.

A control might reduce the mean loss modestly but dramatically compress the tail, or vice versa.

For formal sensitivity analysis, vary one input at a time across its plausible range and record the change in output p90. A tornado chart of these sensitivities reveals which input assumptions matter most – and therefore where better data collection or expert calibration has the highest payoff.

Below we deploy an EDR solution that halves the expected event frequency from \(\lambda=0.6\) to \(\lambda=0.3\) and overlay the two loss distributions.

rng_base = make_rng(99)
rng_edr = make_rng(99)

sev = make_lognormal_severity(10, 2)
losses_base = simulate_aggregate_losses(10_000, 0.6, sev, rng=rng_base)
losses_edr  = simulate_aggregate_losses(10_000, 0.3, sev, rng=rng_edr)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

nz_base = losses_base[losses_base > 0]
nz_edr = losses_edr[losses_edr > 0]
all_nz = np.concatenate([nz_base, nz_edr])
log_bins = np.logspace(np.log10(all_nz.min()), np.log10(all_nz.max()), 60)

axes[0].hist(nz_base, bins=log_bins, alpha=0.5,
             label=f"Baseline (lam=0.6, {(losses_base==0).mean():.0%} zero)", color=ACCENT)
axes[0].hist(nz_edr, bins=log_bins, alpha=0.5,
             label=f"With EDR (lam=0.3, {(losses_edr==0).mean():.0%} zero)", color=DARK_BG)
for data, col, lbl in [(losses_base, ACCENT, "Baseline"), (losses_edr, DARK_BG, "EDR")]:
    p90 = np.percentile(data, 90)
    if p90 > 0:
        axes[0].axvline(p90, color=col, ls=":", lw=1.5, alpha=0.8,
                        label=f"{lbl} p90 = ${p90:,.0f}")
axes[0].set_xscale("log")
axes[0].set_xlabel("Aggregate Loss (USD, log scale)")
axes[0].set_ylabel("Frequency")
axes[0].set_title("Distribution Shift (nonzero losses)")
axes[0].legend(fontsize=8)

# Exceedance curves include the zero-loss years
for data, col, lbl in [(losses_base, ACCENT, "Baseline"), (losses_edr, DARK_BG, "With EDR")]:
    sorted_l = np.sort(data)
    exceed = 1 - np.arange(1, len(sorted_l) + 1) / len(sorted_l)
    axes[1].plot(sorted_l, exceed, color=col, lw=1.5, label=lbl)
axes[1].set_xscale("log")
axes[1].set_yscale("log")
axes[1].set_xlabel("Loss Threshold (USD)")
axes[1].set_ylabel("P(Loss > threshold)")
axes[1].set_title("Loss Exceedance Curve")
axes[1].axhline(0.05, color=LIGHT_GRAY, ls="--", lw=1, label="5% exceedance (p95)")
axes[1].axhline(0.10, color=VERY_LIGHT, ls="--", lw=1, label="10% exceedance (p90)")
axes[1].legend(fontsize=8)
axes[1].set_ylim(1e-4, 1)

plt.tight_layout()
plt.show()

print("Risk bands comparison:")
print(f"{'':>8} {'Baseline':>14} {'With EDR':>14} {'Reduction':>12}")
qs = (0.5, 0.9, 0.95, 0.99)
rb_base = risk_bands(losses_base, qs=qs)
rb_edr  = risk_bands(losses_edr, qs=qs)
for k in rb_base:
    b = rb_base[k]
    e = rb_edr[k]
    pct = (b - e) / b * 100 if b > 0 else 0
    print(f"{k:>8} ${b:>13,.0f} ${e:>13,.0f} {pct:>10.1f}%")

Risk bands comparison:
               Baseline       With EDR    Reduction
     p50 $            0 $            0        0.0%
     p90 $      159,257 $       54,501       65.8%
     p95 $      385,307 $      162,323       57.9%
     p99 $    1,832,680 $      969,807       47.1%

Reproducibility

Monte Carlo results are random – but they should be reproducibly random. Setting a random seed before simulation ensures identical results across runs, which is essential for debugging, peer review, and audit trails. In Python: rng = np.random.default_rng(seed=42). Pass this generator to all sampling calls.

When presenting results, note “N = 50,000, seed = 42” so any reviewer can reproduce the exact output. For production risk reports, run with 3–5 different seeds and confirm that key quantiles are stable.

Always pass an explicit seed via make_rng(seed) in any analysis that will be reported, reviewed, or compared across runs.

sev = make_lognormal_severity(10, 2)

run_a = simulate_aggregate_losses(1_000, 0.6, sev, rng=make_rng(7))
run_b = simulate_aggregate_losses(1_000, 0.6, sev, rng=make_rng(7))

print(f"Identical? {np.array_equal(run_a, run_b)}")
print(f"Run A mean: ${run_a.mean():,.2f}")
print(f"Run B mean: ${run_b.mean():,.2f}")

Identical? True
Run A mean: $95,155.16
Run B mean: $95,155.16

Pitfalls to Watch For

1. Small-N tail noise. With \(N = 1{,}000\), your p95 estimate is based on roughly 50 observations in the tail. That is not enough for a stable estimate, especially with heavy-tailed severity. Increase \(N\) or acknowledge the imprecision explicitly. The demo below shows the instability at N=500.

2. Ignored correlations. The basic compound Poisson model assumes frequency and severity are independent, and that multiple risk scenarios are uncorrelated. In reality, a threat campaign may simultaneously increase frequency and severity, and correlated scenarios (e.g., supply-chain compromise affecting multiple systems) can produce aggregate losses far beyond what independent simulation predicts. Model correlations explicitly when they matter – Part 2 covers copulas and correlated simulation.

3. Unit and time-horizon mismatch. If the model is annual but the event frequency was counted over the last quarter, the resulting ALE will be wrong. Align all inputs to the same time horizon and unit basis before simulating. This sounds obvious but is a common error in practice.

# Pitfall 1: small-N tail instability
sev = make_lognormal_severity(10, 2)
n_trials = 300
p95_estimates = []

for trial in range(n_trials):
    small = simulate_aggregate_losses(500, 0.6, sev, rng=make_rng(trial))
    p95_estimates.append(np.quantile(small, 0.95))

p95_estimates = np.array(p95_estimates)
cv = p95_estimates.std(ddof=1) / p95_estimates.mean()

print(f"P95 across {n_trials} trials (N=500 each):")
print(f"  min  = ${p95_estimates.min():>12,.0f}")
print(f"  max  = ${p95_estimates.max():>12,.0f}")
print(f"  mean = ${p95_estimates.mean():>12,.0f}")
print(f"  coefficient of variation = {cv:.0%}   (max/min = {p95_estimates.max()/max(p95_estimates.min(),1):.1f}x)")

P95 across 300 trials (N=500 each):
  min  = $     189,194
  max  = $     737,405
  mean = $     368,857
  coefficient of variation = 21%   (max/min = 3.9x)

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(p95_estimates, "o", markersize=4, alpha=0.6, color=DARK_BG)
ax.axhline(np.median(p95_estimates), color=ACCENT, linestyle="--", label=f"Median p95 = ${np.median(p95_estimates):,.0f}")
ax.fill_between(range(len(p95_estimates)),
                np.percentile(p95_estimates, 10),
                np.percentile(p95_estimates, 90),
                alpha=0.15, color=ACCENT, label="10th-90th range")
ax.set_xlabel("Trial")
ax.set_ylabel("p95 Estimate ($)")
ax.set_title("p95 Instability at N=500: Each Dot is One Simulation")
ax.legend()
plt.tight_layout()
plt.show()