{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Lesson 6 — Discrete Distributions & Stochastic Simulation\n",
    "\n",
    "**Companion notebook** for the Lesson 6 slides. Because every simulation here uses a\n",
    "**seeded** generator, your numbers will match the checks exactly — as long as you make\n",
    "the same calls in the same order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-001",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check(label, got, expected, tol=1e-6):\n",
    "    \"\"\"Compare your answer with the expected value; prints ✅ or ❌ (never raises).\"\"\"\n",
    "    import math\n",
    "    try:\n",
    "        if isinstance(expected, float) or isinstance(got, float):\n",
    "            ok = math.isclose(float(got), float(expected), rel_tol=tol, abs_tol=tol)\n",
    "        elif isinstance(expected, (list, tuple)) and expected and isinstance(expected[0], str):\n",
    "            ok = list(got) == list(expected)\n",
    "        elif hasattr(expected, \"__len__\") and not isinstance(expected, str):\n",
    "            import numpy as np\n",
    "            ok = np.allclose(np.asarray(got, dtype=float), np.asarray(expected, dtype=float), rtol=tol, atol=tol)\n",
    "        else:\n",
    "            ok = got == expected\n",
    "    except Exception:\n",
    "        ok = False\n",
    "    print((\"✅ \" if ok else \"❌ \") + f\"{label}: got {got!r}, expected {expected!r}\")\n",
    "import math\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "## 1. Seeded generators are reproducible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(42)\n",
    "print(rng.random(3))\n",
    "rng = np.random.default_rng(42)\n",
    "print(rng.random(3))            # identical\n",
    "print(rng.integers(1, 7, size=5), rng.choice([\"TERM\", \"UL\", \"WL\"], size=4, p=[0.6, 0.3, 0.1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "## 2. Bernoulli, Binomial, Poisson — and checking against theory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-005",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(1)\n",
    "print((rng.random(10) < 0.3).astype(int))\n",
    "print(rng.binomial(1, 0.3, size=10))\n",
    "\n",
    "rng = np.random.default_rng(7)\n",
    "x = rng.binomial(n=1000, p=0.01, size=100_000)\n",
    "print(\"Binomial:\", x.mean(), x.var(), \"theory:\", 10.0, 9.9)\n",
    "y = rng.poisson(lam=2.0, size=100_000)\n",
    "print(\"Poisson:\", y.mean(), y.var(), np.quantile(y, [0.5, 0.95, 0.995]))\n",
    "pmf = [math.exp(-2) * 2**k / math.factorial(k) for k in range(4)]\n",
    "print(\"simulated P(Y=k):\", (np.bincount(y)[:4] / len(y)).round(4), \" theory:\", np.round(pmf, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "## 3. How many simulations are enough? The standard error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-007",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(3)\n",
    "for n in (100, 10_000, 1_000_000):\n",
    "    s = rng.poisson(2.0, size=n)\n",
    "    print(f\"{n:>9,}  mean {s.mean():.4f}   s.e. {math.sqrt(2 / n):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-008",
   "metadata": {},
   "source": [
    "---\n",
    "## Your turn — aggregate stop-loss on 10,000 lives\n",
    "\n",
    "Follow the steps **in this order** so your random numbers match the checks.\n",
    "\n",
    "1. `ages = rng.integers(30, 71, size=n_lives)`\n",
    "2. Mortality `q` from the Standard Ultimate Survival Model (Lesson 5)\n",
    "3. `benefits = rng.choice([50_000, 100_000, 250_000, 500_000], size=n_lives, p=[0.3, 0.4, 0.2, 0.1])`\n",
    "4. Exact `ES` $=\\sum q_k b_k$ and `SD` $=\\sqrt{\\sum q_k(1-q_k)b_k^2}$\n",
    "5. `deaths = rng.random((n_sims, n_lives)) < q`, then annual totals `S = deaths @ benefits`\n",
    "6. Stop-loss premium $E[(S-d)_+]$ with retention $d = 1.25 \\times ES$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-009",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(2026)\n",
    "n_lives, n_sims = 10_000, 5_000\n",
    "A_, B_, c_ = 0.00022, 2.7e-6, 1.124\n",
    "ES = SD = S = stop_loss = None\n",
    "...  # your code here\n",
    "\n",
    "check(\"E[S]\", ES, 3951407.01, tol=1e-8)\n",
    "check(\"SD[S]\", SD, 1055701.26, tol=1e-8)\n",
    "check(\"simulated mean of S\", None if S is None else S.mean(), 3952540.0)\n",
    "check(\"stop-loss premium at 1.25 × E[S]\", stop_loss, 110871.12, tol=1e-7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "**Extend it:** repeat for retentions of 1.0, 1.5, and 2.0 × E[S]; add P(S > d) and the\n",
    "standard error of each premium; and plot a histogram of `S`. Which premium would you\n",
    "trust least — and why?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}