{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "91834d9c",
   "metadata": {},
   "source": [
    "# DCT Laboratory — Volume II, Chapter 11\n",
    "## Distributionally Robust Enterprise Optimization\n",
    "**Seed `26211`** · Companion to the chapter and AXIOM Module **AXIOM-11 (Vol. II)**\n",
    "\n",
    "Ambiguity with teeth — and everything in closed form. The **worst case over a\n",
    "total-variation ball**: mass slides from the best outcome to the worst,\n",
    "piecewise-linearly, kinking at $\\delta = 0.4$. The **decision flip**: nominal\n",
    "analysis picks bold project B (6.8 vs 4.8); DRO at $\\delta = 0.15$ picks safe\n",
    "A (4.5 vs 4.1), with the flip at exactly $\\delta^* = 0.125$. And the\n",
    "**data-driven radius** $\\delta_n = 0.55/\\sqrt{n}$: twenty samples buy back\n",
    "the bold choice. Mirrored in `DCT_V2_Ch11_Lab.xlsx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8701d3b1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:40:30.416703Z",
     "iopub.status.busy": "2026-07-14T00:40:30.416521Z",
     "iopub.status.idle": "2026-07-14T00:40:30.874113Z",
     "shell.execute_reply": "2026-07-14T00:40:30.872919Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.dpi']=110\n",
    "\n",
    "import numpy as np\n",
    "SEED = 26211\n",
    "Z = np.array([10.0, 6.0, 2.0, -4.0]); PHAT = np.array([0.4, 0.3, 0.2, 0.1])\n",
    "def wc_mean(z, phat, delta):\n",
    "    \"\"\"Worst-case E[z] over the TV ball {p : ||p - phat||_TV <= delta} (minimizing adversary):\n",
    "    slide mass from highest-z outcomes to the lowest-z outcome.\"\"\"\n",
    "    p = phat.astype(float).copy()\n",
    "    order = np.argsort(z)[::-1]           # highest payoff first\n",
    "    lo = int(np.argmin(z)); budget = delta\n",
    "    for i in order:\n",
    "        if i == lo or budget <= 0: continue\n",
    "        take = min(p[i], budget)\n",
    "        p[i] -= take; p[lo] += take; budget -= take\n",
    "    return float(p @ z)\n",
    "# Panel 2: two projects on the same support probabilities\n",
    "ZA = np.array([5.0, 5.0, 5.0, 3.0]); ZB = np.array([12.0, 8.0, 1.0, -6.0])\n",
    "def flip_delta():\n",
    "    # wc_A = 4.8 - 2d ; wc_B = 6.8 - 18d ; equal at d = 2/16\n",
    "    return 2.0/16.0\n",
    "# Panel 3: data-driven radius\n",
    "def delta_n(n): return 0.55/np.sqrt(n)\n",
    "def choice(n):\n",
    "    d = delta_n(n)\n",
    "    return \"B\" if wc_mean(ZB, PHAT, d) > wc_mean(ZA, PHAT, d) else \"A\"\n",
    "def n_switch():\n",
    "    n = 1\n",
    "    while choice(n) == \"A\": n += 1\n",
    "    return n\n",
    "\n",
    "def reference_values():\n",
    "    return {\n",
    "        \"nominal_mean\": round(float(PHAT @ Z),4),\n",
    "        \"wc_mean_d01\": round(wc_mean(Z, PHAT, 0.10),4),\n",
    "        \"wc_slope\": round((float(PHAT@Z)-wc_mean(Z,PHAT,0.10))/0.10,4),\n",
    "        \"kink_delta\": 0.4,\n",
    "        \"wc_mean_d05\": round(wc_mean(Z, PHAT, 0.50),4),\n",
    "        \"nom_A\": round(float(PHAT @ ZA),4), \"nom_B\": round(float(PHAT @ ZB),4),\n",
    "        \"wc_A_015\": round(wc_mean(ZA, PHAT, 0.15),4),\n",
    "        \"wc_B_015\": round(wc_mean(ZB, PHAT, 0.15),4),\n",
    "        \"dro_flip_delta\": round(flip_delta(),4),\n",
    "        \"delta_n20\": round(delta_n(20),4),\n",
    "        \"n_switch\": n_switch(),\n",
    "        \"delta_n100\": round(delta_n(100),4),\n",
    "        \"wc_B_n100\": round(wc_mean(ZB, PHAT, delta_n(100)),4),\n",
    "    }\n",
    "if __name__ == \"__main__\":\n",
    "    [print(f\"{k:16s} {v}\") for k,v in reference_values().items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cac17059",
   "metadata": {},
   "source": [
    "## Panel 1 — The worst-case distribution is extremal\n",
    "Support $z = (10, 6, 2, -4)$, nominal $\\hat{p} = (0.4, 0.3, 0.2, 0.1)$,\n",
    "ambiguity = the TV ball of radius $\\delta$ around $\\hat{p}$. The adversary's\n",
    "optimal move (Worst-Case Distribution Existence Theorem) is boundary-simple:\n",
    "slide mass $\\delta$ from the *highest* payoff to the *lowest*. So the\n",
    "worst-case mean falls **linearly at slope $z_{max} - z_{min} = 14$** until the\n",
    "best outcome's mass is exhausted at $\\delta = 0.4$ — then the slope kinks to\n",
    "the next spread. Ambiguity Sets Capture Distributional Uncertainty (Prop.),\n",
    "priced: each unit of radius costs 14 units of guaranteed mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "649bf54a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:40:30.876708Z",
     "iopub.status.busy": "2026-07-14T00:40:30.875868Z",
     "iopub.status.idle": "2026-07-14T00:40:31.107025Z",
     "shell.execute_reply": "2026-07-14T00:40:31.105988Z"
    }
   },
   "outputs": [],
   "source": [
    "ds = np.linspace(0, 0.6, 121)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.2))\n",
    "ax.plot(ds, [wc_mean(Z, PHAT, d) for d in ds], c=\"#C8A24B\", lw=2.4)\n",
    "ax.axvline(0.4, c=\"#8A8F8B\", ls=\":\", lw=1.4, label=\"kink: best outcome's mass exhausted\")\n",
    "ax.scatter([0.1],[wc_mean(Z,PHAT,0.1)], c=\"#0B3D2E\", s=60, zorder=5, label=f\"delta=0.1: {wc_mean(Z,PHAT,0.1):.2f}\")\n",
    "ax.set(xlabel=\"ambiguity radius delta (total variation)\", ylabel=\"worst-case E[z]\",\n",
    "       title=\"Nominal 5.80, sliding at slope 14 — seed 26211\")\n",
    "ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"nominal {PHAT@Z:.4f}   wc(0.1) = {wc_mean(Z,PHAT,0.1):.4f}   wc(0.5) = {wc_mean(Z,PHAT,0.5):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2e1fb18",
   "metadata": {},
   "source": [
    "## Panel 2 — The decision flip\n",
    "Two projects on the same states: safe $A = (5,5,5,3)$ and bold\n",
    "$B = (12,8,1,-6)$. Nominal means: $A = 4.8$, $B = 6.8$ — **nominal analysis\n",
    "funds B**. But B's payoff *spread* is 18 against A's 2, so the worst case\n",
    "punishes B nine times harder per unit of radius: the lines cross at\n",
    "$\\delta^* = 2/16 = 0.125$, and at $\\delta = 0.15$ DRO funds A (4.5 vs 4.1).\n",
    "DREO Generalizes Stochastic and Robust Optimization (Prop.): $\\delta = 0$\n",
    "recovers Chapter 9, $\\delta \\to$ everything recovers worst-case — and the\n",
    "interesting decisions live between."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9abbd5a2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:40:31.108961Z",
     "iopub.status.busy": "2026-07-14T00:40:31.108775Z",
     "iopub.status.idle": "2026-07-14T00:40:31.297286Z",
     "shell.execute_reply": "2026-07-14T00:40:31.296537Z"
    }
   },
   "outputs": [],
   "source": [
    "ds = np.linspace(0, 0.3, 61)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.2))\n",
    "ax.plot(ds, [wc_mean(ZA, PHAT, d) for d in ds], c=\"#1B6B52\", lw=2.2, label=\"safe A: 4.8 − 2·delta\")\n",
    "ax.plot(ds, [wc_mean(ZB, PHAT, d) for d in ds], c=\"#C8A24B\", lw=2.2, label=\"bold B: 6.8 − 18·delta\")\n",
    "ax.axvline(flip_delta(), c=\"#B0532F\", ls=\":\", lw=1.5, label=f\"flip at delta* = {flip_delta():.3f}\")\n",
    "ax.set(xlabel=\"ambiguity radius delta\", ylabel=\"worst-case value\",\n",
    "       title=\"Nominal funds B; robustness past 0.125 funds A — seed 26211\")\n",
    "ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"nominal: A {PHAT@ZA:.4f}, B {PHAT@ZB:.4f}   at delta=0.15: A {wc_mean(ZA,PHAT,0.15):.4f}, B {wc_mean(ZB,PHAT,0.15):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02fa217c",
   "metadata": {},
   "source": [
    "## Panel 3 — Data buys boldness\n",
    "The radius is not a mood — it is calibrated (Data-Driven Enterprise\n",
    "Optimization, Def.; Statistical Confidence Set, Def.): with $n$ observations,\n",
    "$\\delta_n = 0.55/\\sqrt{n}$, shrinking as evidence accumulates (the\n",
    "Finite-Sample Enterprise Performance Guarantee Theorem's shape). Few samples:\n",
    "the prudent radius exceeds 0.125 and DRO funds safe A. At **$n = 20$** the\n",
    "radius drops below the flip and the bold project earns its funding — with a\n",
    "guarantee attached. Data-Driven Ambiguity Sets Improve Enterprise Adaptability\n",
    "(Prop.): the same policy machinery, growing bolder exactly as fast as the\n",
    "evidence justifies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e7cce1f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:40:31.299748Z",
     "iopub.status.busy": "2026-07-14T00:40:31.299553Z",
     "iopub.status.idle": "2026-07-14T00:40:31.482551Z",
     "shell.execute_reply": "2026-07-14T00:40:31.480964Z"
    }
   },
   "outputs": [],
   "source": [
    "ns = np.arange(1, 41)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.2))\n",
    "ax.plot(ns, [delta_n(n) for n in ns], \"o-\", c=\"#C8A24B\", lw=2, ms=3.5, label=\"delta_n = 0.55/sqrt(n)\")\n",
    "ax.axhline(flip_delta(), c=\"#B0532F\", ls=\":\", lw=1.5, label=\"flip threshold 0.125\")\n",
    "ax.axvline(n_switch(), c=\"#0B3D2E\", ls=\"--\", lw=1.4, label=f\"n = {n_switch()}: bold B funded\")\n",
    "ax.set(xlabel=\"sample size n\", ylabel=\"calibrated radius\", title=\"The radius shrinks; the decision matures — seed 26211\")\n",
    "ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"delta(20) = {delta_n(20):.4f}   delta(100) = {delta_n(100):.4f}   wc_B at n=100: {wc_mean(ZB,PHAT,delta_n(100)):.4f}\")\n",
    "print(f\"choice(10) = {choice(10)}   choice(20) = {choice(20)}   choice(100) = {choice(100)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aeac5110",
   "metadata": {},
   "source": [
    "## Validation — agrees with `DCT_V2_Ch11_Lab.xlsx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb332868",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:40:31.484549Z",
     "iopub.status.busy": "2026-07-14T00:40:31.484364Z",
     "iopub.status.idle": "2026-07-14T00:40:31.491323Z",
     "shell.execute_reply": "2026-07-14T00:40:31.490347Z"
    }
   },
   "outputs": [],
   "source": [
    "ref = reference_values()\n",
    "expected = {\"nominal_mean\":5.8,\"wc_mean_d01\":4.4,\"wc_slope\":14.0,\"kink_delta\":0.4,\n",
    " \"wc_mean_d05\":-0.8,\"nom_A\":4.8,\"nom_B\":6.8,\"wc_A_015\":4.5,\"wc_B_015\":4.1,\n",
    " \"dro_flip_delta\":0.125,\"delta_n20\":0.123,\"n_switch\":20,\"delta_n100\":0.055,\"wc_B_n100\":5.81}\n",
    "for k,v in expected.items():\n",
    "    assert abs(ref[k]-v)<5e-4, f\"MISMATCH {k}\"\n",
    "    print(f\"PASS  {k:16s} {ref[k]}\")\n",
    "print(\"\\nAll checkpoints agree — seed 26211.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f0272f7",
   "metadata": {},
   "source": [
    "**Next**: Exercises 11.5–11.9 (Part C) resize the support and watch the kink structure multiply; AXIOM-11's adversary console lets you fight the worst-case distribution by hand. Chapter 12 turns from one contested objective to several: Pareto. Solutions: IM Vol. II, Ch. 11."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}