Python for Actuarial Students · Course 1

Lesson 6

Discrete Distributions & Stochastic Simulation

Reproducible randomness, claim-count models, and a 10,000-life portfolio simulated five thousand times.

Lesson 6 · OverviewBy the end of this lesson you can…

Generate random numbers with np.random.default_rng
Explain seeds and why a “random” result can be reproduced exactly
Simulate Bernoulli, Binomial, and Poisson claim counts
Compare empirical and theoretical means, variances, and quantiles
Judge simulation accuracy with the standard error
Price aggregate stop-loss cover by simulation

75 min

Concepts & worked examples

75 min

Lab: stop-loss on 10,000 lives

Companion notebook: lesson-06-discrete-simulation.ipynb — worked examples and self-checking exercises.

Part 1

Pseudo-random Numbers

Deterministic sequences that behave like chance.

GeneratorsThe modern NumPy way

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 normal

Create one generator per analysis and pass it around
Every method takes size= — ask for a million draws at once, not one at a time in a loop
default_rng() with no seed gives different numbers every run

You’ll see old code using np.random.seed() and np.random.rand(). That global-state API still works but is legacy — use a Generator.

MechanicsHow a seed produces “randomness”

Seed42

→

Internal statea large number, set from the seed

→step

Next statea fixed mathematical update

→

Outputbits → 0.77395605…

The algorithm (NumPy’s default is PCG64) is completely deterministic
Its outputs pass statistical tests for randomness, and the sequence doesn’t repeat for 2¹²⁸ draws
Same seed + same sequence of calls ⇒ identical results, on any machine

Record 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.

Part 2

Discrete Claim-Count Models

Bernoulli, Binomial, and Poisson.

DistributionsThree workhorses

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.

Empirical vs. theoryDoes the simulation match the math?

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⁻²2^k/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.

AccuracyHow many simulations are enough?

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

The law of large numbers: sample means converge to the true mean
The standard error σ/√n measures how far off a sample mean typically is
100× more simulations buys only 10× more accuracy

Report a simulated result with its standard error: “stop-loss premium 110,871 ± 4,819”. It tells the reader which digits mean anything.

Actuarial Lab · 75 minAggregate stop-loss on 10,000 lives

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]$$

Build a portfolio: ages, mortality rates, sums assured
Compute E[S] and Var[S] exactly
Simulate 5,000 years of deaths for every life
Price the stop-loss at several retentions

Why simulate?

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.

Lab · Steps 1–2The portfolio and exact moments

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 b_k · Bernoulli(q_k), 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 Σq_k	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.

Lab · Steps 3–4Simulate, then price

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

Tasks

Reproduce both tables — then change the seed. Which digits move?
Plot a histogram of S (matplotlib). Is it symmetric?
Add a standard-error column. Is the 2.0× premium estimated reliably?

Lab · ExtensionsGoing further

Poisson shortcut

Replace per-life Bernoullis with a Poisson(Σq) death count per year and randomly chosen benefits. How close is the stop-loss premium?

Concentration risk

Give 20 lives a 5,000,000 benefit. How much do the tail probabilities change? Why?

Memory

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.

RecapLesson 6 checklist

I create one seeded Generator and pass it around
I can explain why a seeded result is reproducible
I can simulate Bernoulli, Binomial, and Poisson counts
I compare simulated moments and quantiles with theory

I quote simulation results with a standard error
I computed E[S] and Var[S] with dot products
I simulated 5,000 × 10,000 lives without a loop
I priced stop-loss cover at several retentions

Up next

Lesson 7

Continuous Distributions & Loss Generation

scipy.stats: PDFs, CDFs, and inverse CDFs — and generating realistic, long-tailed claim severities.

← Course home Next: Lesson 7 →