{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1fe07d7b",
   "metadata": {},
   "source": [
    "# Model Monitoring & Drift Lab\n",
    "\n",
    "Companion notebook to [Data Science Module 6](https://rodoslay.com/data-science/06-production-ml-for-risk) at **rodoslay.com**.\n",
    "\n",
    "A monitoring system built from scratch and stress-tested: PSI/CSI with decile decomposition,\n",
    "a **drift injector** (mean shift / new segment / pipeline null-bug) so you learn each signature,\n",
    "and the **retraining trap** — a model retrained mid-drift on lagged labels that loses to a simple\n",
    "recalibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "703baa71",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "rng = np.random.default_rng(21)\n",
    "plt.rcParams[\"figure.figsize\"] = (10, 4.5)\n",
    "plt.rcParams[\"axes.grid\"] = True\n",
    "plt.rcParams[\"grid.alpha\"] = 0.3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a99f4fa",
   "metadata": {},
   "source": [
    "## 1. PSI from scratch, with the decomposition that makes it useful"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84e26ca4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(expected, actual, bins=10):\n",
    "    \"\"\"PSI of `actual` vs `expected`, using expected's decile edges.\n",
    "    Returns (psi_total, per-bin DataFrame).\"\"\"\n",
    "    edges = np.quantile(expected, np.linspace(0, 1, bins + 1))\n",
    "    edges[0], edges[-1] = -np.inf, np.inf\n",
    "    e = np.histogram(expected, edges)[0] / len(expected)\n",
    "    a = np.histogram(actual, edges)[0] / len(actual)\n",
    "    a = np.maximum(a, 1e-6); e = np.maximum(e, 1e-6)\n",
    "    contrib = (a - e) * np.log(a / e)\n",
    "    table = pd.DataFrame({\"bin\": range(1, bins + 1), \"expected %\": e.round(4),\n",
    "                          \"actual %\": a.round(4), \"PSI contribution\": contrib.round(4)})\n",
    "    return contrib.sum(), table\n",
    "\n",
    "baseline = rng.normal(0, 1, 50_000)      # development-time score distribution\n",
    "same = rng.normal(0, 1, 20_000)\n",
    "value, _ = psi(baseline, same)\n",
    "print(f\"PSI, same population resampled: {value:.4f}   (noise floor — anything near this is nothing)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "475e153b",
   "metadata": {},
   "source": [
    "## 2. The drift injector — learn each signature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f583910",
   "metadata": {},
   "outputs": [],
   "source": [
    "def inject(kind, n=20_000):\n",
    "    if kind == \"mean shift\":       return rng.normal(0.30, 1, n)\n",
    "    if kind == \"wider\":            return rng.normal(0, 1.35, n)\n",
    "    if kind == \"new segment\":      return np.r_[rng.normal(0, 1, int(n*0.85)), rng.normal(-1.8, 0.5, int(n*0.15))]\n",
    "    if kind == \"pipeline bug\":     # a feature nulls out -> scores collapse to a constant band\n",
    "        x = rng.normal(0, 1, n); x[rng.random(n) < 0.25] = 0.0; return x\n",
    "\n",
    "rows = []\n",
    "for kind in [\"mean shift\", \"wider\", \"new segment\", \"pipeline bug\"]:\n",
    "    v, table = psi(baseline, inject(kind))\n",
    "    worst = table.iloc[table[\"PSI contribution\"].idxmax()]\n",
    "    rows.append([kind, round(v, 3),\n",
    "                 \"stable\" if v < 0.10 else \"monitor\" if v < 0.25 else \"INVESTIGATE\",\n",
    "                 f\"bin {int(worst['bin'])}: expected {worst['expected %']:.1%}, actual {worst['actual %']:.1%}\"])\n",
    "pd.DataFrame(rows, columns=[\"injected drift\", \"PSI\", \"status\", \"worst bin (the pointer)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "593fab84",
   "metadata": {},
   "source": [
    "Each failure mode has a fingerprint: a mean shift leaks mass across every boundary; a new\n",
    "segment slams one tail bin; a pipeline bug spikes the middle decile where the constant lands.\n",
    "The *decomposition*, not the total, is what tells you where to look.\n",
    "\n",
    "## 3. A two-year production run, with the outcome lag made explicit\n",
    "\n",
    "A model ships at month 0. Concept drift begins at month 8 (a coefficient rotates). Labels take\n",
    "**6 months** to resolve — so at any month *m* you can only measure real performance for loans\n",
    "booked before *m − 6*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be9231b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_PER_MONTH, N_MONTHS, LAG = 4_000, 24, 6\n",
    "\n",
    "def make_month(m, n=N_PER_MONTH):\n",
    "    X = rng.standard_normal((n, 3))\n",
    "    w1 = 1.3 if m < 8 else 1.3 - 1.5 * min((m - 8) / 10, 1)     # drifting coefficient\n",
    "    logit = -2.4 + w1 * X[:, 0] + 0.8 * X[:, 1] + 0.5 * X[:, 2]\n",
    "    y = (rng.random(n) < 1 / (1 + np.exp(-logit))).astype(int)\n",
    "    return X, y\n",
    "\n",
    "months = [make_month(m) for m in range(N_MONTHS)]\n",
    "X0 = np.vstack([months[m][0] for m in range(6)])\n",
    "y0 = np.hstack([months[m][1] for m in range(6)])\n",
    "model = LogisticRegression(max_iter=1000).fit(X0, y0)          # trained on months 0-5\n",
    "\n",
    "rows = []\n",
    "for m in range(6, N_MONTHS):\n",
    "    Xm, ym = months[m]\n",
    "    p = model.predict_proba(Xm)[:, 1]\n",
    "    score_psi, _ = psi(model.predict_proba(X0)[:, 1], p)\n",
    "    known = m - LAG      # latest month whose labels have resolved\n",
    "    if known >= 6:\n",
    "        Xk, yk = months[known]\n",
    "        pk = model.predict_proba(Xk)[:, 1]\n",
    "        auc_known = roc_auc_score(yk, pk)\n",
    "        gap_known = yk.mean() - pk.mean()\n",
    "    else:\n",
    "        auc_known, gap_known = np.nan, np.nan\n",
    "    auc_true = roc_auc_score(ym, p)                             # what you CAN'T see yet\n",
    "    rows.append([m, score_psi, auc_known, gap_known, auc_true])\n",
    "\n",
    "log = pd.DataFrame(rows, columns=[\"month\", \"score PSI (today)\", \"AUC (lagged view)\",\n",
    "                                  \"calib gap (lagged)\", \"true AUC (invisible today)\"])\n",
    "log.round(3).head(18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b3c254a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\n",
    "ax1.plot(log[\"month\"], log[\"score PSI (today)\"], \"o-\", color=\"#1f4d3a\")\n",
    "ax1.axhline(0.10, color=\"#8a6d1a\", ls=\"--\", lw=0.8); ax1.axhline(0.25, color=\"#b3341e\", ls=\"--\", lw=0.8)\n",
    "ax1.axvline(8, color=\"k\", ls=\":\", lw=0.8); ax1.set_title(\"Score PSI — available same-day\")\n",
    "ax2.plot(log[\"month\"], log[\"true AUC (invisible today)\"], \"o-\", color=\"#b3341e\", label=\"true AUC (unknowable)\")\n",
    "ax2.plot(log[\"month\"], log[\"AUC (lagged view)\"], \"o-\", color=\"#5c5c5c\", label=\"AUC as visible (6m lag)\")\n",
    "ax2.axvline(8, color=\"k\", ls=\":\", lw=0.8); ax2.legend(); ax2.set_title(\"Performance — what you see arrives late\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e2e3d03",
   "metadata": {},
   "source": [
    "The picture is the module's core claim: drift starts at month 8; **PSI reacts within a month or\n",
    "two**, while the lagged performance view stays serene long after the true AUC has sagged. The input\n",
    "monitor is your smoke detector; the outcome monitor is the fire report.\n",
    "\n",
    "## 4. The retraining trap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "446d4b60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# stand at month 14 (drift began at month 8). Your choices:\n",
    "DECIDE = 14\n",
    "labeled = range(0, DECIDE - LAG)              # only months 0..7 have resolved labels (pre-drift!)\n",
    "\n",
    "# option A: keep the original model\n",
    "# option B: \"retrain on the latest data\" -> latest LABELED data = pre-drift months\n",
    "Xr = np.vstack([months[m][0] for m in labeled]); yr = np.hstack([months[m][1] for m in labeled])\n",
    "retrained = LogisticRegression(max_iter=1000).fit(Xr, yr)\n",
    "\n",
    "# option C: recalibrate the original model's intercept on EARLY outcomes (30dpd-style proxy at 2m lag)\n",
    "# proxy: use months DECIDE-2 with partially-resolved (noisy) labels\n",
    "Xe, ye = months[DECIDE - 2]\n",
    "ye_noisy = np.where(rng.random(len(ye)) < 0.75, ye, rng.random(len(ye)) < ye.mean())  # imperfect early labels\n",
    "p_e = model.predict_proba(Xe)[:, 1]\n",
    "shift = np.log(ye_noisy.mean() / (1 - ye_noisy.mean())) - np.log(p_e.mean() / (1 - p_e.mean()))\n",
    "\n",
    "rows = []\n",
    "for m_eval in range(DECIDE, DECIDE + 4):\n",
    "    Xm, ym = months[m_eval]\n",
    "    pA = model.predict_proba(Xm)[:, 1]\n",
    "    pB = retrained.predict_proba(Xm)[:, 1]\n",
    "    zC = np.log(pA / (1 - pA)) + shift\n",
    "    pC = 1 / (1 + np.exp(-zC))\n",
    "    rows.append([m_eval, roc_auc_score(ym, pA), roc_auc_score(ym, pB),\n",
    "                 abs(ym.mean() - pA.mean()), abs(ym.mean() - pC.mean())])\n",
    "pd.DataFrame(rows, columns=[\"month\", \"AUC: original\", \"AUC: 'retrained on latest'\",\n",
    "                            \"|calib gap|: original\", \"|calib gap|: recalibrated\"]).round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65c08c0b",
   "metadata": {},
   "source": [
    "The retrained model is no better on ranking — its \"latest data\" was all pre-drift, so it learned\n",
    "the same stale physics with fresh confidence. The cheap intercept recalibration, fed by *early*\n",
    "outcome proxies, meaningfully fixes the calibration gap. Moral: **retraining is not a reflex; check\n",
    "what the retraining data actually contains.**\n",
    "\n",
    "## Exercises\n",
    "\n",
    "1. **CSI triage.** Add per-feature PSIs to the month loop. When score PSI first breaches 0.10,\n",
    "   which feature's CSI points at the cause?\n",
    "2. **Alert budget.** Tune thresholds so the two-year run produces at most three alerts — do you\n",
    "   still catch month 8? What did you trade?\n",
    "3. **Override simulation.** Add underwriters who overturn 5% of decisions at random pre-drift,\n",
    "   but 15% of the model's approvals post-drift. Build the override-rate monitor; when does it fire\n",
    "   relative to PSI?\n",
    "4. **The declined blind spot.** Track feature drift separately for applicants below the cutoff.\n",
    "   Inject drift *only* there. Which of your monitors sees it? (This motivates the randomized\n",
    "   approval slice.)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}