{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Lesson 3 — GLMs for Severity & Pure Premium\n",
    "\n",
    "**Companion notebook** for the Lesson 3 slides. You'll fit Gamma and inverse Gaussian\n",
    "severity models, combine frequency and severity into a pure premium, write a rating\n",
    "function, check it on the 2025 holdout, and fit a Tweedie model in one step."
   ]
  },
  {
   "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 numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as st\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "df = pd.read_parquet(\"../data/processed/motor_model.parquet\")\n",
    "bins = [17, 21, 25, 30, 40, 50, 65, 75, 120]\n",
    "labels = [\"17-20\", \"21-24\", \"25-29\", \"30-39\", \"40-49\", \"50-64\", \"65-74\", \"75+\"]\n",
    "df[\"age_band\"] = pd.cut(df[\"driver_age\"], bins, right=False, labels=labels)\n",
    "df[\"log_miles\"] = np.log(df[\"annual_mileage\"] / 12_000)\n",
    "df[\"young\"] = df[\"driver_age\"] < 25\n",
    "train = df[df[\"policy_year\"] <= 2024].copy()\n",
    "hold = df[df[\"policy_year\"] == 2025].copy()\n",
    "\n",
    "claims = pd.read_csv(\"../data/motor/claims.csv\")\n",
    "claims[\"incurred\"] = claims[\"paid\"] + claims[\"outstanding\"]\n",
    "liab = claims[claims[\"coverage\"] == \"LIAB\"].merge(\n",
    "    train[[\"vehicle_id\", \"vehicle_style\", \"age_band\", \"young\", \"urban\", \"log_miles\", \"channel\"]],\n",
    "    on=\"vehicle_id\", how=\"inner\", validate=\"many_to_one\")\n",
    "print(f\"{len(liab):,} liability claims; mean {liab['incurred'].mean():,.0f}, \"\n",
    "      f\"median {liab['incurred'].median():,.0f}, CV {liab['incurred'].std() / liab['incurred'].mean():.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "---\n",
    "## Your turn 1 — a Gamma severity GLM\n",
    "\n",
    "Fit `incurred ~ C(vehicle_style, Treatment('Sedan')) + young + urban` with a Gamma\n",
    "family and an explicit **log** link. Store the model as `sev`. What are the Sports\n",
    "relativity and the dispersion φ (`.scale`)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "SEV = \"C(vehicle_style, Treatment('Sedan')) + young + urban\"\n",
    "gamma_log = sm.families.Gamma(sm.families.links.Log())\n",
    "sev = None\n",
    "...  # your code here\n",
    "check(\"Sports severity relativity\",\n",
    "      None if sev is None else round(np.exp(sev.params[\"C(vehicle_style, Treatment('Sedan'))[T.Sports]\"]), 3), 1.204)\n",
    "check(\"dispersion\", None if sev is None else round(sev.scale, 3), 1.832)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "## 1. Do the extra factors help?\n",
    "\n",
    "An F-test compares the full severity model with the parsimonious one (Gamma models have\n",
    "an estimated dispersion, so use F rather than χ²)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-005",
   "metadata": {},
   "outputs": [],
   "source": [
    "FULL = (\"C(vehicle_style, Treatment('Sedan')) + C(age_band, Treatment('40-49')) + urban\"\n",
    "        \" + log_miles + C(channel, Treatment('Agent'))\")\n",
    "sev_full = smf.glm(\"incurred ~ \" + FULL, data=liab, family=gamma_log).fit()\n",
    "if sev is not None:\n",
    "    df_num = sev.df_resid - sev_full.df_resid\n",
    "    F = ((sev.deviance - sev_full.deviance) / df_num) / sev_full.scale\n",
    "    print(f\"F = {F:.3f} on {df_num:.0f} df, p = {st.f.sf(F, df_num, sev_full.df_resid):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "## Your turn 2 — which family fits best?\n",
    "\n",
    "Fit the same formula with an inverse Gaussian family (log link) and compare AIC with the\n",
    "Gamma. Then compare with a lognormal: fit `np.log(incurred)` by OLS and convert its\n",
    "log-likelihood to the dollar scale by subtracting `np.log(liab[\"incurred\"]).sum()`.\n",
    "Store the three AICs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-007",
   "metadata": {},
   "outputs": [],
   "source": [
    "aic_gamma = aic_ig = aic_lognormal = None\n",
    "if sev is not None:\n",
    "    ...  # your code here\n",
    "check(\"Gamma AIC\", None if aic_gamma is None else round(aic_gamma, 1), 85001.1)\n",
    "check(\"inverse Gaussian AIC\", None if aic_ig is None else round(aic_ig, 1), 84687.6)\n",
    "check(\"lognormal AIC\", None if aic_lognormal is None else round(aic_lognormal, 1), 83865.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-008",
   "metadata": {},
   "source": [
    "## 2. The frequency model from Lesson 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-009",
   "metadata": {},
   "outputs": [],
   "source": [
    "FREQ = (\"C(vehicle_style, Treatment('Sedan')) + C(age_band, Treatment('40-49'))\"\n",
    "        \" + urban + log_miles\")\n",
    "freq = smf.glm(\"n_liab ~ \" + FREQ, data=train, family=sm.families.Poisson(),\n",
    "               offset=np.log(train[\"exposure\"])).fit()\n",
    "np.exp(freq.params).round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "## Your turn 3 — pure premium and a rating function\n",
    "\n",
    "1. Add `pp_hat` = predicted frequency (per car-year: use `offset=np.zeros(len(d))`)\n",
    "   × predicted severity to `train`, and check expected vs actual losses.\n",
    "2. Write `rate(style, age_band, young, urban, miles)` returning the pure premium for one\n",
    "   risk. Price a Sedan · 40–49 · rural · 12,000 miles and a Sports · 17–20 · young ·\n",
    "   urban · 15,000 miles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-011",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rate(style, age_band, young, urban, miles):\n",
    "    ...  # your code here\n",
    "\n",
    "\n",
    "base = sports = None\n",
    "if sev is not None:\n",
    "    train[\"pp_hat\"] = freq.predict(train, offset=np.zeros(len(train))) * sev.predict(train)\n",
    "    print(f\"expected {(train['pp_hat'] * train['exposure']).sum():,.0f} vs actual {train['inc_liab'].sum():,.0f}\")\n",
    "    base = rate(\"Sedan\", \"40-49\", False, False, 12_000)\n",
    "    sports = rate(\"Sports\", \"17-20\", True, True, 15_000)\n",
    "check(\"base pure premium\", None if base is None else round(base, 2), 447.03)\n",
    "check(\"young urban sports pure premium\", None if sports is None else round(sports, 2), 3883.87)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-012",
   "metadata": {},
   "source": [
    "## Your turn 4 — the holdout year\n",
    "\n",
    "Predict pure premium for the 2025 holdout and compute actual ÷ expected losses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "ae = None\n",
    "if sev is not None:\n",
    "    ...  # your code here\n",
    "check(\"holdout A/E\", None if ae is None else round(ae, 3), 1.022)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-014",
   "metadata": {},
   "source": [
    "## 3. Where the current rates go wrong"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-015",
   "metadata": {},
   "outputs": [],
   "source": [
    "if sev is not None:\n",
    "    display(train.groupby(\"urban\").apply(lambda g: pd.Series({\n",
    "        \"premium_per_car_year\": g[\"ep_liab\"].sum() / g[\"exposure\"].sum(),\n",
    "        \"actual_pp\": g[\"inc_liab\"].sum() / g[\"exposure\"].sum(),\n",
    "        \"model_pp\": (g[\"pp_hat\"] * g[\"exposure\"]).sum() / g[\"exposure\"].sum(),\n",
    "        \"loss_ratio\": g[\"inc_liab\"].sum() / g[\"ep_liab\"].sum()})).round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-016",
   "metadata": {},
   "source": [
    "## Your turn 5 — a Tweedie GLM on pure premium\n",
    "\n",
    "Model `pp = inc_liab / exposure` with `var_weights=exposure` and a Tweedie family with a\n",
    "log link, for power p in 1.3, 1.4, …, 1.8. Keep the power with the highest\n",
    "log-likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-017",
   "metadata": {},
   "outputs": [],
   "source": [
    "train[\"pp\"] = train[\"inc_liab\"] / train[\"exposure\"]\n",
    "best_p = None\n",
    "...  # your code here\n",
    "check(\"best Tweedie power\", best_p, 1.6)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}