Reproducible randomness, claim-count models, and a 10,000-life portfolio simulated five thousand times.
np.random.default_rngCompanion notebook: lesson-06-discrete-simulation.ipynb — worked examples and self-checking exercises.
Deterministic sequences that behave like chance.
import numpy as np
rng = np.random.default_rng(42) # create a Generator
rng.random(3) # uniform on [0, 1)
# [0.77395605 0.43887844 0.85859792]
rng = np.random.default_rng(42) # same seed …
rng.random(3) # … same numbers
# [0.77395605 0.43887844 0.85859792]
rng.integers(1, 7, size=5) # dice
rng.choice(["TERM", "UL", "WL"], size=4,
p=[0.6, 0.3, 0.1]) # weighted pick
rng.normal(0, 1, size=2) # standard normalsize= — ask for a million draws at once, not one at a time in a loopdefault_rng() with no seed gives different numbers every runYou’ll see old code using np.random.seed() and np.random.rand(). That global-state API still works but is legacy — use a Generator.
42Record the seed with every simulation result. “Rerun the capital model with seed 2026” must reproduce last quarter’s number to the cent — that’s what an audit trail means for stochastic models.
Bernoulli, Binomial, and Poisson.
| Distribution | Models | Mean | Variance | NumPy |
|---|---|---|---|---|
| Bernoulli(p) | does one life die / one policy claim? | p | p(1−p) | rng.random(n) < p |
| Binomial(n, p) | deaths among n identical lives | np | np(1−p) | rng.binomial(n, p, size) |
| Poisson(λ) | claim counts with no fixed upper limit | λ | λ | rng.poisson(lam, size) |
rng = np.random.default_rng(1)
(rng.random(10) < 0.3).astype(int)
# [0 0 1 0 0 0 0 0 0 1]
rng.binomial(1, 0.3, size=10) # Bernoulli = Binomial(1, p)
# [1 0 0 1 0 0 0 0 0 0]rng.random(n) < p is the most flexible Bernoulli: p can be an array — a different probability for every life. You’ll need that in the lab.
rng = np.random.default_rng(7)
x = rng.binomial(n=1000, p=0.01, size=100_000)
x.mean(), x.var() # (10.0038, 9.8718...)
# theory: (10.0, 9.9)
np.quantile(x, [0.5, 0.95, 0.995]) # [10. 15. 19.]
y = rng.poisson(lam=2.0, size=100_000)
y.mean(), y.var() # (1.99862, 1.99215...)
np.quantile(y, [0.5, 0.95, 0.995]) # [2. 5. 6.]
np.bincount(y)[:4] / len(y)
# [0.13556 0.27077 0.27009 0.18024]| k | Simulated P(Y = k) | Theory e−22k/k! |
|---|---|---|
| 0 | 0.1356 | 0.1353 |
| 1 | 0.2708 | 0.2707 |
| 2 | 0.2701 | 0.2707 |
| 3 | 0.1802 | 0.1804 |
Theoretical Poisson(2) quantiles at 50%, 95%, 99.5%: 2, 5, 6 — matched exactly.
Always check a simulation against something you can calculate. If the easy case is wrong, the hard case is too.
rng = np.random.default_rng(3)
for n in (100, 10_000, 1_000_000):
s = rng.poisson(2.0, size=n)
print(n, s.mean(), np.sqrt(2 / n))| n | Sample mean | Std. error σ/√n |
|---|---|---|
| 100 | 2.1700 | 0.1414 |
| 10,000 | 1.9942 | 0.0141 |
| 1,000,000 | 2.0037 | 0.0014 |
Report a simulated result with its standard error: “stop-loss premium 110,871 ± 4,819”. It tells the reader which digits mean anything.
A reinsurer offers stop-loss cover: it pays the amount by which a year’s total death claims S exceed a retention d.
$$\text{Stop-loss premium} = E\big[(S - d)_+\big]$$
E[S] has a simple formula. But E[(S − d)+] depends on the whole shape of S’s distribution — a sum of 10,000 differently weighted Bernoullis. There’s no neat closed form; simulation gets it in seconds.
rng = np.random.default_rng(2026)
n_lives = 10_000
ages = rng.integers(30, 71, size=n_lives) # 30..70
A, B, c = 0.00022, 2.7e-6, 1.124 # Lesson 5's SUSM
q = 1 - np.exp(-A - B * c**ages * (c - 1) / np.log(c))
benefits = rng.choice([50_000, 100_000, 250_000, 500_000],
size=n_lives, p=[0.3, 0.4, 0.2, 0.1])
ES = q @ benefits # E[S]
VS = (q * (1 - q)) @ benefits.astype(float) ** 2 # Var[S]Each life k contributes bk · Bernoulli(qk), independently:
$$E[S] = \sum_k q_k b_k \qquad \operatorname{Var}[S] = \sum_k q_k(1-q_k)\,b_k^2$$
| Expected deaths Σqk | 24.657 |
| E[S] | 3,951,407.01 |
| SD[S] | 1,055,701.26 |
Dot products (Lesson 5) do the sums. .astype(float) before squaring: NumPy integers have a fixed size and silently wrap around on overflow; floats don’t.
n_sims = 5_000
deaths = rng.random((n_sims, n_lives)) < q # (5000, 10000)
S = deaths @ benefits # (5000,)
S.mean(), S.std() # (3952540.0, 1057259.45)
d = 1.25 * ES
payout = np.maximum(S - d, 0)
payout.mean() # 110,871.12
(S > d).mean() # 0.1742
payout.std() / np.sqrt(n_sims) # 4,818.86 (s.e.)Broadcasting compares each row of uniforms against the (10000,) vector q. The boolean matrix is 50 MB — fine for a laptop.
| Retention | d | Stop-loss premium | P(S > d) |
|---|---|---|---|
| 1.00 × E[S] | 3,951,407 | 421,951.52 | 0.4680 |
| 1.25 × E[S] | 4,939,259 | 110,871.12 | 0.1742 |
| 1.50 × E[S] | 5,927,111 | 20,299.80 | 0.0380 |
| 2.00 × E[S] | 7,902,814 | 446.62 | 0.0012 |
Replace per-life Bernoullis with a Poisson(Σq) death count per year and randomly chosen benefits. How close is the stop-loss premium?
Give 20 lives a 5,000,000 benefit. How much do the tail probabilities change? Why?
Simulate in chunks of 1,000 years so 100,000 simulations fit in memory. Confirm the answers agree within standard error.
This is a genuine reinsurance pricing method. Real models add expenses, a risk margin, and a better mortality basis — the simulation engine stays the same.
Generator and pass it aroundscipy.stats: PDFs, CDFs, and inverse CDFs — and generating realistic, long-tailed claim severities.