Insurance Claims Simulator

/**
 * Fintech Math Bootcamp · Module 07 demo · Insurance Claims Simulator
 * A small insurer covers 1,000 synthetic home policies. How many claims arrive in a year, when, how big
 * are they, and how much should it hold for a bad year? Three claim-size models share the same average
 * claim, and a seeded Monte Carlo shows how different their 1-in-100 year losses are.
 * Lessons 061–070: PMF/CDF/survival/quantile, Bernoulli and binomial, Poisson, uniform sampling, the normal,
 * lognormal, Student-t, exponential and gamma waiting times, mixtures, Monte Carlo.
 * Synthetic parameters; not a pricing or reserving model for any real insurer.
 */

// 064 · a seeded uniform generator (mulberry32): the same seed always gives the same draws
export function mulberry32(seed: number): () => number {
  let a = seed >>> 0;
  return () => {
    a = (a + 0x6d2b79f5) >>> 0;
    let t = a;
    t = Math.imul(t ^ (t >>> 15), t | 1);
    t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
    return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
  };
}
// 064 · scale a uniform draw to any interval
export const uniformBetween = (a: number, b: number, u: number) => a + (b - a) * u;

// 065 · standard normal draw by Box-Muller, and the standard normal CDF (Abramowitz-Stegun erf)
export function normalDraw(rand: () => number): number {
  const u1 = 1 - rand(), u2 = rand();
  return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
}
export function phi(z: number): number {
  const x = Math.abs(z) / Math.SQRT2, t = 1 / (1 + 0.3275911 * x);
  const erf = 1 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * Math.exp(-x * x);
  return z >= 0 ? 0.5 * (1 + erf) : 0.5 * (1 - erf);
}

// 062 · binomial PMF for n comparable policies, each claiming with probability p
export function binomialPmf(n: number, p: number, kMax: number): number[] {
  const out = [(1 - p) ** n];
  for (let k = 0; k < kMax; k++) out.push(out[k] * (n - k) / (k + 1) * p / (1 - p));
  return out;
}
// 063 · Poisson PMF for claim counts over one year of exposure
export function poissonPmf(lambda: number, kMax: number): number[] {
  const out = [Math.exp(-lambda)];
  for (let k = 0; k < kMax; k++) out.push(out[k] * lambda / (k + 1));
  return out;
}
export function poissonDraw(lambda: number, rand: () => number): number {
  const limit = Math.exp(-lambda);
  let k = 0, prod = rand();
  while (prod > limit) { k++; prod *= rand(); }
  return k;
}
// 061 · CDF and survival read from a PMF
export const cdfAt = (pmf: number[], k: number) => pmf.slice(0, k + 1).reduce((s, v) => s + v, 0);

// 068 · exponential waiting time by inverse CDF, and the gamma wait for the k-th claim
export const exponentialWait = (scale: number, u: number) => -scale * Math.log(1 - u);
export const exponentialSurvival = (scale: number, t: number) => Math.exp(-t / scale);
export const gammaMoments = (shape: number, scale: number) => ({mean: shape * scale, sd: Math.sqrt(shape) * scale});

// 066 · lognormal claim sizes: log(size) is normal(mu, sigma)
export const lognormal = (mu: number, sigma: number) => ({
  median: Math.exp(mu), mean: Math.exp(mu + sigma * sigma / 2), mode: Math.exp(mu - sigma * sigma),
  sd: Math.exp(mu + sigma * sigma / 2) * Math.sqrt(Math.exp(sigma * sigma) - 1),
});
export const lognormalSurvival = (mu: number, sigma: number, x: number) => 1 - phi((Math.log(x) - mu) / sigma);

// 067 · Student-t with 4 degrees of freedom: closed-form CDF, variance ν/(ν−2) = 2
export function t4Cdf(t: number): number {
  const s = t * t / 4;
  return 0.5 + 0.375 * (t / Math.sqrt(1 + s)) * (1 - s / (3 * (1 + s)));
}
const t4Pdf = (t: number) => 0.375 * (1 + t * t / 4) ** -2.5;
// draw: Z / sqrt(chi-square(4) / 4), where chi-square(4) = −2 ln(U1 U2)
export const t4Draw = (rand: () => number) => normalDraw(rand) / Math.sqrt(-2 * Math.log((1 - rand()) * (1 - rand())) / 4);
// expected claim when log(size) = mu + s·T and the policy pays at most `limit` (Simpson's rule)
export function logTCappedMean(mu: number, s: number, limit: number): number {
  const xc = (Math.log(limit) - mu) / s, a = -60, n = 6000, h = (xc - a) / n;
  let acc = 0;
  for (let i = 0; i <= n; i++) {
    const x = a + i * h, w = i === 0 || i === n ? 1 : i % 2 ? 4 : 2;
    acc += w * Math.exp(mu + s * x) * t4Pdf(x);
  }
  return acc * h / 3 + limit * (1 - t4Cdf(xc));
}

// 069 · a two-group mixture: mean and variance from the components
export function mixture(weights: number[], means: number[], variances: number[]) {
  const mean = weights.reduce((s, w, i) => s + w * means[i], 0);
  const variance = weights.reduce((s, w, i) => s + w * (variances[i] + (means[i] - mean) ** 2), 0);
  return {mean, variance, sd: Math.sqrt(variance)};
}

// find mu so that f(mu) hits a target (f increasing)
function solve(f: (mu: number) => number, target: number): number {
  let lo = -5, hi = 20;
  for (let i = 0; i < 80; i++) { const mid = (lo + hi) / 2; if (f(mid) < target) lo = mid; else hi = mid; }
  return (lo + hi) / 2;
}

// 061 · empirical quantile Q(p) = smallest simulated value with F(x) >= p
export const quantile = (sorted: number[], p: number) => sorted[Math.min(sorted.length - 1, Math.max(0, Math.ceil(p * sorted.length) - 1))];

// 070 · Monte Carlo: simulate many years of total claims for one claim-size model
export function simulateYears(years: number, lambda: number, claim: (rand: () => number) => number, seed: number) {
  const rand = mulberry32(seed), totals: number[] = [], firstYear = {count: 0, claims: [] as number[], total: 0};
  for (let y = 0; y < years; y++) {
    const k = poissonDraw(lambda, rand);
    let s = 0;
    for (let i = 0; i < k; i++) { const x = claim(rand); s += x; if (y === 0) firstYear.claims.push(x); }
    if (y === 0) { firstYear.count = k; firstYear.total = s; }
    totals.push(s);
  }
  const sorted = [...totals].sort((a, b) => a - b), n = years;
  const mean = totals.reduce((s, v) => s + v, 0) / n;
  const sd = Math.sqrt(totals.reduce((s, v) => s + (v - mean) ** 2, 0) / (n - 1));
  const bins = 24, top = 2_400_000, hist = Array(bins).fill(0);
  for (const v of totals) hist[Math.min(bins - 1, Math.floor(v / top * bins))]++;
  return {
    mean, sd, standardError: sd / Math.sqrt(n),
    median: quantile(sorted, 0.5), q90: quantile(sorted, 0.9), q99: quantile(sorted, 0.99), q995: quantile(sorted, 0.995),
    survivalAt1M: totals.filter(v => v > 1_000_000).length / n, max: sorted[n - 1],
    hist: hist.map(c => c / n), histTop: top, firstYear,
  };
}

export function runDemo() {
  const policies = 1000, pClaim = 0.05, lambda = policies * pClaim, meanClaim = 10000, limit = 1_000_000, years = 10000;

  // 062 / 063 · claim counts
  const binom = binomialPmf(policies, pClaim, 90), pois = poissonPmf(lambda, 90);
  const counts = {
    policies, pClaim, lambda,
    binomialMean: policies * pClaim, binomialVariance: policies * pClaim * (1 - pClaim), poissonVariance: lambda,
    atMean: {binomial: binom[lambda], poisson: pois[lambda]},
    cdf60: cdfAt(pois, 60), survival60: 1 - cdfAt(pois, 60),
    q99: pois.findIndex((_, k) => cdfAt(pois, k) >= 0.99),
    poissonPmf: pois.slice(25, 80), binomialPmf: binom.slice(25, 80), pmfFrom: 25,
  };

  // 064 / 068 · when claims arrive: seeded uniforms turned into exponential waits (days)
  const scale = 365 / lambda, r = mulberry32(7);
  const us = Array.from({length: 8}, () => r());
  const waits = us.map(u => exponentialWait(scale, u));
  const arrival = {
    scale, uniforms: us, waits, claimDays: waits.reduce((acc: number[], w) => [...acc, (acc.at(-1) ?? 0) + w], []),
    dayOfYear: us.map(u => uniformBetween(0, 365, u)),
    survival14: exponentialSurvival(scale, 14), fifthClaim: gammaMoments(5, scale),
  };

  // 065 / 066 · lognormal sizes with a 10,000 average claim
  const sigma = 0.8, mu = Math.log(meanClaim) - sigma * sigma / 2, ln = lognormal(mu, sigma);
  const z30k = (Math.log(30000) - mu) / sigma;
  const logN = {mu, sigma, ...ln, z30k, survival30k: 1 - phi(z30k), survival100k: lognormalSurvival(mu, sigma, 100000)};

  // 069 · mixture: 96% ordinary customers, 4% catastrophe-exposed customers with 12 times the average claim
  const w = [0.96, 0.04], sig = [0.7, 1.0], ratio = 12;
  const meanO = meanClaim / (w[0] + w[1] * ratio), means = [meanO, meanO * ratio];
  const mus = means.map((m, i) => Math.log(m) - sig[i] ** 2 / 2);
  const vars = means.map((m, i) => (Math.exp(sig[i] ** 2) - 1) * m * m);
  const mix = {weights: w, sigmas: sig, means, medians: mus.map(Math.exp), ...mixture(w, means, vars),
    survival100k: w[0] * lognormalSurvival(mus[0], sig[0], 100000) + w[1] * lognormalSurvival(mus[1], sig[1], 100000)};

  // 067 · Student-t in log space (ν = 4, scale 0.8): the average exists only because of the policy limit
  const s = 0.8, muT = solve(m => logTCappedMean(m, s, limit), meanClaim);
  const logT = {nu: 4, scale: s, mu: muT, median: Math.exp(muT), limit, tVariance: 4 / (4 - 2), cappedMean: logTCappedMean(muT, s, limit),
    survival100k: 1 - t4Cdf((Math.log(100000) - muT) / s), survivalLimit: 1 - t4Cdf((Math.log(limit) - muT) / s)};

  // 070 · ten thousand simulated years for each claim-size model, each with its own fixed seed
  const sims = {
    lognormal: simulateYears(years, lambda, rr => Math.exp(mu + sigma * normalDraw(rr)), 101),
    mixture: simulateYears(years, lambda, rr => { const g = rr() < w[1] ? 1 : 0; return Math.exp(mus[g] + sig[g] * normalDraw(rr)); }, 202),
    studentT: simulateYears(years, lambda, rr => Math.min(limit, Math.exp(muT + s * t4Draw(rr))), 303),
  };

  return {meanClaim, expectedAnnualLoss: lambda * meanClaim, years, counts, arrival, logN, mix, logT, sims};
}

export const checkedResult = {"meanClaim":10000,"expectedAnnualLoss":500000,"years":10000,"counts":{"policies":1000,"pClaim":0.05,"lambda":50,"binomialMean":50,"binomialVariance":47.5,"poissonVariance":50,"atMean":{"binomial":0.05778798371410466,"poisson":0.056325006325190816},"cdf60":0.9278398201867428,"survival60":0.07216017981325717,"q99":67,"poissonPmf":[0.00003705785984050479,0.00007126511507789382,0.00013197243532943301,0.00023566506308827324,0.0004063190742901263,0.0006771984571502104,0.0010922555760487266,0.0017066493375761354,0.0025858323296608114,0.00380269460244237,0.0054324208606319575,0.007545028973099941,0.010195985098783705,0.013415769866820663,0.01719970495746239,0.021499631196827986,0.02621906243515608,0.031213169565661995,0.03629438321588604,0.041243617290779584,0.04582624143419954,0.049811131993695155,0.052990565950739525,0.055198506198687,0.05632500632519082,0.056325006325190816,0.05522059443646159,0.0530967254196746,0.05009125039591943,0.0463807874036291,0.04216435218511736,0.03764674302242622,0.03302345879160195,0.028468498958277542,0.02412584657481148,0.020104872145676234,0.016479403398095276,0.013289841450076835,0.010547493214346695,0.008240229073708355,0.006338637749006427,0.004801998294701838,0.00358358081694167,0.0026349858948100514,0.0019094100687029358,0.0013638643347878114,0.0009604678413998673,0.0006669915565276856,0.0004568435318682778,0.0003086780620731607,0.00020578537471544045,0.0001353851149443687,0.00008791241230153812,0.000056354110449703924,0.00003566715851247084],"binomialPmf":[0.000027086568961362176,0.000053460333476372716,0.00010150168578165114,0.00018564124110065144,0.00032748327831185695,0.0005578706372645843,0.0009187343262252068,0.0014642328324214234,0.002260569986896233,0.0033838563116542686,0.004915496536929359,0.006934874500200046,0.00950955763612069,0.01268380055898092,0.016466688444992775,0.020821694204786922,0.025659597479583372,0.03083653381318353,0.03615838358999978,0.041391833846447124,0.04628139550550111,0.050570632388734056,0.054025065284269085,0.05645382370165399,0.05772721822124019,0.05778798371410466,0.056654885994220264,0.054418508915501054,0.0512301355033714,0.04728551493342371,0.04280583457130988,0.038018339915308126,0.033138793056372004,0.028357424548238478,0.023829343375950624,0.019669659751552226,0.01595295959142286,0.012716323477373571,0.009964838280514964,0.007678497918456022,0.00581949315925088,0.004339095776634428,0.0031835942304607673,0.0022989887128636965,0.0016343687874820485,0.0011440581512374341,0.0007887131806158738,0.0005356100473626806,0.0003583605796341512,0.00023627329823674122,0.00015353619239805082,0.00009835247781731097,0.00006211735441093323,0.00003868712423838825,0.000023763843136438355],"pmfFrom":25},"arrival":{"scale":7.3,"uniforms":[0.011704753153026104,0.06195825757458806,0.97690763277933,0.6990287057124078,0.5214452685322613,0.4055216880515218,0.4662326325196773,0.23992518591694534],"waits":[0.08594868915435112,0.46691405495906085,27.508247917647967,8.765404822153492,5.379988261070922,3.7965186261636585,4.582904787456116,2.0026704000346554],"claimDays":[0.08594868915435112,0.552862744113412,28.06111066176138,36.826515483914875,42.206503744985795,46.003022371149456,50.58592715860557,52.588597558640224],"dayOfYear":[4.272234900854528,22.614764014724642,356.57128596445546,255.14547758502886,190.32752301427536,148.01541613880545,170.1749108696822,87.57269285968505],"survival14":0.14692864485854654,"fifthClaim":{"mean":36.5,"sd":16.323296235748465}},"logN":{"mu":8.890340371976183,"sigma":0.8,"median":7261.490370736914,"mean":10000.00000000001,"mode":3828.9288597511204,"sd":9468.267419675858,"z30k":1.7732653608351368,"survival30k":0.03809232599744261,"survival100k":0.000522358476235385},"mix":{"weights":[0.96,0.04],"sigmas":[0.7,1],"means":[6944.444444444444,83333.33333333333],"medians":[5435.448182235202,50544.22164271945],"mean":10000,"variance":730648481.0513726,"sd":27030.51018851425,"survival100k":0.009915973323981303},"logT":{"nu":4,"scale":0.8,"mu":8.382111628496077,"median":4368.223258865045,"limit":1000000,"tVariance":2,"cappedMean":9999.999999999922,"survival100k":0.008671129956251633,"survivalLimit":0.001227148334422945},"sims":{"lognormal":{"mean":499654.25360023405,"sd":97881.54304530079,"standardError":978.8154304530078,"median":493309.5327436938,"q90":628026.853164651,"q99":749692.5829425494,"q995":785409.7335921137,"survivalAt1M":0,"max":931010.314957701,"hist":[0,0.0001,0.01,0.1424,0.3743,0.3227,0.1215,0.0252,0.0033,0.0005,0,0,0,0,0,0,0,0,0,0,0,0,0,0],"histTop":2400000,"firstYear":{"count":37,"claims":[16542.84785219984,11708.121507531934,1159.5004424593126,7678.779460426996,7349.18066152039,4189.2534380909965,22632.905765773292,4353.448860250661,1198.1712298732793,7701.389836256539,8858.920735613981,2661.263086278998,13572.01581659995,8394.50470573574,5686.447449369472,3551.020470730281,9533.694820702734,19371.652078102095,3597.105444105638,7298.103770820749,4901.83760233824,13463.14169196204,10728.609309679767,2577.5244919217157,2561.3911041984256,5326.614270065563,11546.853010500363,19226.203474246355,1794.328423745996,9088.509234671004,23764.641237973097,5726.8243872079565,2322.3950663674086,9234.960818745572,9760.528600961494,11290.222533696084,3746.9827005764546],"total":314099.89539130044}},"mixture":{"mean":499806.0269839704,"sd":196538.92457437803,"standardError":1965.3892457437803,"median":454029.5272963864,"q90":734033.2761909914,"q99":1224570.3041031733,"q995":1353569.6691777022,"survivalAt1M":0.0246,"max":2543044.242792573,"hist":[0,0.0017,0.0713,0.2754,0.2633,0.1672,0.0974,0.053,0.0299,0.0162,0.0076,0.0058,0.0048,0.0024,0.001,0.0007,0.001,0.0002,0.0004,0.0001,0.0002,0.0001,0.0001,0.0002],"histTop":2400000,"firstYear":{"count":52,"claims":[11917.074420495963,1956.6532563078802,7092.045767687936,1631.9290852068805,7979.51033576637,11087.555746188715,12556.861695148038,12012.542558559278,3427.783762045298,5659.594821628765,15424.34629691286,7061.418312473099,13491.008338689591,3257.3605097888785,5707.976390005595,2713.6622696026434,2822.1315386563415,4027.1632653402467,3081.304106742231,7551.982220006485,11311.607729697424,11721.110247967867,6855.581886511196,2658.1676498385023,1091.752156809583,9721.917241329971,2055.2084081672233,2812.3273188880326,4110.878376281802,23667.869975210982,2416.8833070279607,4532.129935838706,22819.128722248675,33202.94853863699,248873.0580879698,4503.226778632149,5081.29359435715,9955.93403279336,15561.626654284657,6716.560357182989,4417.799507880968,3985.4459282269036,3290.5074225683484,9000.116410053206,5732.354011938669,13622.692964700615,4235.452097595659,3344.611532226504,3804.2487368802895,18286.870663038957,12411.131025387222,1605.5591882579058],"total":653865.9051856834}},"studentT":{"mean":498422.8754359541,"sd":322738.30742791994,"standardError":3227.3830742791993,"median":391907.3777397704,"q90":890637.4398487129,"q99":1642166.4356619008,"q995":1899191.6245694186,"survivalAt1M":0.0873,"max":3400536.269275151,"hist":[0,0.013,0.1904,0.3159,0.1968,0.0926,0.0477,0.0271,0.0175,0.0117,0.0083,0.0072,0.0188,0.0211,0.0141,0.006,0.0038,0.0018,0.0012,0.0011,0.0007,0.0009,0.0006,0.0017],"histTop":2400000,"firstYear":{"count":60,"claims":[12563.403136770383,2090.178479713406,12971.604231462557,3479.942253047861,4392.528678659179,4309.9221833596175,69752.21646410748,6304.154103902798,13389.597838273596,4043.962185908174,898.0545068321697,2258.7899835172307,2295.643463339357,7941.835184433551,12060.330630114255,5745.749405925179,1256.1253434983462,5685.500977598058,2962.998715494391,1122.9991039314323,7325.513065906155,5949.80129219574,13018.892316118763,6780.54181540207,7919.600804929286,9602.474366338405,1024.6344211550036,2909.0432896357634,9485.993394037265,2358.1361863242837,4657.568535022494,9675.987847934886,4139.117201232215,6350.9560409554,9930.135180848953,13085.008517853825,4801.757158247268,8468.062342517982,6474.470933545774,6045.25418163941,31217.858782048228,1631.4882995143025,1776.105926727912,15480.944340356627,7702.020963524603,3931.057058421665,824.3557086771999,4344.852159156454,5888.397759317364,2407.62434283574,9729.104397761334,583.7266357496864,6580.4315952381185,542.517377880413,724.4224670253213,2128.7087597306713,9964.419934808267,6064.549660522049,2680.564695388044,9242.73377603425],"total":438974.3704024482}}}};

// Run this file directly: npx tsx lessons/07-probability-distributions-and-simulation-basics/demo-insurance-claims-simulator.ts
if (process.argv[1] && import.meta.url.endsWith(process.argv[1].replace(/\\/g, "/").split("/").pop()!)) {
  console.log(JSON.stringify(runDemo(), null, 2));
}