{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Lesson 8A — Production Pipelines & Model Comparison\n",
    "\n",
    "**Companion notebook** for the Lesson 8A slides. You'll build scikit-learn Pipelines for\n",
    "a GLM and a gradient-boosted model, tune a penalty, wrap frequency and severity into one\n",
    "pure-premium object, save it with `joblib`, and run a champion-vs-challenger test on\n",
    "the 2025 holdout."
   ]
  },
  {
   "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 tempfile\n",
    "from dataclasses import dataclass, field\n",
    "from pathlib import Path\n",
    "\n",
    "import joblib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from sklearn.compose import ColumnTransformer\n",
    "from sklearn.ensemble import HistGradientBoostingRegressor\n",
    "from sklearn.linear_model import GammaRegressor, PoissonRegressor\n",
    "from sklearn.metrics import mean_poisson_deviance\n",
    "from sklearn.model_selection import GroupKFold\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import FunctionTransformer, OneHotEncoder, OrdinalEncoder\n",
    "\n",
    "df = pd.read_parquet(\"../data/processed/motor_model.parquet\")\n",
    "edges = [16, 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\"], edges, right=False, labels=labels).astype(str)\n",
    "cov = pd.read_csv(\"../data/motor/coverages.csv\")\n",
    "df[\"liab_premium\"] = df[\"vehicle_id\"].map(cov[cov[\"coverage\"] == \"LIAB\"].set_index(\"vehicle_id\")[\"annual_premium\"])\n",
    "train = df[df[\"policy_year\"] <= 2024].copy()\n",
    "hold = df[df[\"policy_year\"] == 2025].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "## 1. A GLM pipeline\n",
    "\n",
    "A `ColumnTransformer` prepares each column its own way; the `Pipeline` chains it to the\n",
    "model. Use a named function (not a lambda) so the pipeline can be saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_miles(miles):\n",
    "    return np.log(miles / 12_000)\n",
    "\n",
    "\n",
    "def log_miles_name(transformer, names):\n",
    "    return [\"log_miles\"]\n",
    "\n",
    "\n",
    "def glm_prep(columns):\n",
    "    cats = [c for c in (\"vehicle_style\", \"age_band\") if c in columns]\n",
    "    base = {\"vehicle_style\": \"Sedan\", \"age_band\": \"40-49\"}\n",
    "    steps = [(\"cat\", OneHotEncoder(drop=[base[c] for c in cats], handle_unknown=\"ignore\",\n",
    "                                   sparse_output=False), cats), (\"urban\", \"passthrough\", [\"urban\"])]\n",
    "    if \"annual_mileage\" in columns:\n",
    "        steps.append((\"log_miles\", FunctionTransformer(log_miles, feature_names_out=log_miles_name),\n",
    "                      [\"annual_mileage\"]))\n",
    "    return ColumnTransformer(steps, verbose_feature_names_out=False)\n",
    "\n",
    "\n",
    "def glm_frequency(alpha):\n",
    "    return Pipeline([(\"prep\", glm_prep([\"vehicle_style\", \"age_band\", \"urban\", \"annual_mileage\"])),\n",
    "                     (\"glm\", PoissonRegressor(alpha=alpha, solver=\"newton-cholesky\", tol=1e-8, max_iter=1000))])\n",
    "\n",
    "\n",
    "def glm_severity(alpha):\n",
    "    return Pipeline([(\"prep\", glm_prep([\"vehicle_style\", \"urban\"])),\n",
    "                     (\"glm\", GammaRegressor(alpha=alpha, solver=\"newton-cholesky\", tol=1e-8, max_iter=1000))])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "---\n",
    "## Your turn 1 — the same model in two libraries\n",
    "\n",
    "Fit `glm_frequency(alpha=0.0)` on the claim **rate** with exposure as the sample weight\n",
    "(`glm__sample_weight=...`). Fit the equivalent statsmodels GLM with an offset. Store the\n",
    "largest absolute difference between the two sets of coefficients (intercept included)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-005",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_diff = None\n",
    "...  # your code here\n",
    "check(\"coefficients agree to 1e-4\", None if max_diff is None else max_diff < 1e-4, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "## Your turn 2 — tune the penalty\n",
    "\n",
    "For α in [0, 1e-4, 1e-3, 1e-2, 1e-1], compute the mean Poisson deviance over\n",
    "`GroupKFold(5, shuffle=True, random_state=0)` folds grouped by `policyholder_id`.\n",
    "Predicted claims on the left-out fold = predicted rate × exposure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-007",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_alpha = None\n",
    "...  # your code here\n",
    "check(\"best alpha\", best_alpha, 0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-008",
   "metadata": {},
   "source": [
    "## 2. A boosted challenger, and one object for pure premium\n",
    "\n",
    "`PurePremiumModel` holds a frequency pipeline and a severity pipeline, and rescales each\n",
    "so that predicted totals equal actual totals on the training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-009",
   "metadata": {},
   "outputs": [],
   "source": [
    "GBM_FEATURES = [\"vehicle_style\", \"driver_age\", \"years_licensed\", \"n_drivers\", \"vehicle_age\",\n",
    "                \"vehicle_value\", \"annual_mileage\", \"urban\", \"coast_km\"]\n",
    "\n",
    "\n",
    "def gbm(loss):\n",
    "    prep = ColumnTransformer([(\"style\", OrdinalEncoder(), [\"vehicle_style\"]),\n",
    "                              (\"num\", \"passthrough\", GBM_FEATURES[1:])],\n",
    "                             verbose_feature_names_out=False).set_output(transform=\"pandas\")\n",
    "    return Pipeline([(\"prep\", prep), (\"gbm\", HistGradientBoostingRegressor(\n",
    "        loss=loss, categorical_features=[\"vehicle_style\"], learning_rate=0.05, max_iter=400,\n",
    "        max_leaf_nodes=8, min_samples_leaf=400, l2_regularization=1.0, random_state=0))])\n",
    "\n",
    "\n",
    "@dataclass\n",
    "class PurePremiumModel:\n",
    "    name: str\n",
    "    frequency: Pipeline\n",
    "    severity: Pipeline\n",
    "    metadata: dict = field(default_factory=dict)\n",
    "    freq_scale: float = 1.0\n",
    "    sev_scale: float = 1.0\n",
    "\n",
    "    def fit(self, frame):\n",
    "        self.frequency.fit(frame, frame[\"n_liab\"] / frame[\"exposure\"],\n",
    "                           **{f\"{self.frequency.steps[-1][0]}__sample_weight\": frame[\"exposure\"]})\n",
    "        claims = frame[frame[\"n_liab\"] > 0]\n",
    "        self.severity.fit(claims, claims[\"inc_liab\"] / claims[\"n_liab\"],\n",
    "                          **{f\"{self.severity.steps[-1][0]}__sample_weight\": claims[\"n_liab\"]})\n",
    "        self.freq_scale = self.sev_scale = 1.0\n",
    "        self.freq_scale = frame[\"n_liab\"].sum() / (self.predict_frequency(frame) * frame[\"exposure\"]).sum()\n",
    "        self.sev_scale = claims[\"inc_liab\"].sum() / (self.predict_severity(claims) * claims[\"n_liab\"]).sum()\n",
    "        self.metadata.update(trained_rows=len(frame), sklearn_version=sklearn.__version__)\n",
    "        return self\n",
    "\n",
    "    def predict_frequency(self, frame):\n",
    "        return self.frequency.predict(frame) * self.freq_scale\n",
    "\n",
    "    def predict_severity(self, frame):\n",
    "        return self.severity.predict(frame) * self.sev_scale\n",
    "\n",
    "    def predict_pure_premium(self, frame):\n",
    "        return self.predict_frequency(frame) * self.predict_severity(frame)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "## Your turn 3 — fit the champion and the challenger\n",
    "\n",
    "Build `models = {\"glm\": PurePremiumModel(\"glm\", glm_frequency(1e-4), glm_severity(1e-4)),\n",
    "\"gbm\": PurePremiumModel(\"gbm\", gbm(\"poisson\"), gbm(\"gamma\"))}` and fit both on `train`.\n",
    "Look at the balance factors: which model needed correcting, and by how much?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-011",
   "metadata": {},
   "outputs": [],
   "source": [
    "models = None\n",
    "...  # your code here\n",
    "check(\"GBM severity balance factor\", None if models is None else round(models[\"gbm\"].sev_scale, 3), 1.051)\n",
    "check(\"GLM frequency balance factor\", None if models is None else round(models[\"glm\"].freq_scale, 3), 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-012",
   "metadata": {},
   "source": [
    "## Your turn 4 — save, reload, verify\n",
    "\n",
    "Save the GLM model with `joblib.dump` to a temporary folder, load it back, and check that\n",
    "the reloaded model gives identical pure premiums on the holdout."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "identical = None\n",
    "if models is not None:\n",
    "    ...  # your code here\n",
    "check(\"reloaded predictions identical\", identical, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-014",
   "metadata": {},
   "source": [
    "## Your turn 5 — the A/B test on holdout loss ratios\n",
    "\n",
    "For each model on the 2025 holdout:\n",
    "1. actual ÷ expected losses (expected = pure premium × exposure);\n",
    "2. the Gini (see `gini` below);\n",
    "3. actual loss ratio (`inc_liab` ÷ `ep_liab`) by quintile of expected loss ÷ current\n",
    "   earned premium (`ep_liab`). Use `pd.qcut(ratio.rank(method=\"first\"), 5, labels=range(1, 6))`.\n",
    "\n",
    "Store results in a dict `ab` like `{\"glm\": {\"ae\": ..., \"gini\": ..., \"lr\": Series}, ...}`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-015",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gini(losses, predicted, exposure):\n",
    "    order = np.argsort(predicted)\n",
    "    cum_e = np.r_[0, np.cumsum(exposure[order]) / exposure.sum()]\n",
    "    cum_l = np.r_[0, np.cumsum(losses[order]) / losses.sum()]\n",
    "    return float(1 - 2 * np.trapezoid(cum_l, cum_e))\n",
    "\n",
    "\n",
    "ab = None\n",
    "if models is not None:\n",
    "    ...  # your code here\n",
    "check(\"GLM A/E\", None if ab is None else round(ab[\"glm\"][\"ae\"], 3), 1.014)\n",
    "check(\"GLM Gini\", None if ab is None else round(ab[\"glm\"][\"gini\"], 3), 0.275)\n",
    "check(\"GBM Gini\", None if ab is None else round(ab[\"gbm\"][\"gini\"], 3), 0.235)\n",
    "check(\"GLM quintile-5 loss ratio\", None if ab is None else round(ab[\"glm\"][\"lr\"].iloc[-1], 3), 1.204)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}