{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2e7969ae",
   "metadata": {},
   "source": [
    "# DCT Laboratory — Volume II, Chapter 9\n",
    "## Stochastic Enterprise Optimization\n",
    "**Seed `26209`** · Companion to the chapter and AXIOM Module **AXIOM-09 (Vol. II)**\n",
    "\n",
    "Volume I built the randomness; this chapter finally optimizes on it — and the\n",
    "lab shows uncertainty doing three different things. It **flips a policy**\n",
    "(Chapter 7's machine with drift risk: gentle stops being optimal at\n",
    "$p = 0.55$). It **reprices without changing the policy** (log-utility Merton:\n",
    "$c^* = \\rho x$ survives $\\sigma$; the value drops by exactly\n",
    "$\\sigma^2/2\\rho^2 = 2.0$). And it is **purchasable as a guarantee** (chance\n",
    "constraints: 95% certainty costs 3.49 extra units). Mirrored in\n",
    "`DCT_V2_Ch09_Lab.xlsx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc7fc9ff",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:28:12.018083Z",
     "iopub.status.busy": "2026-07-14T00:28:12.017896Z",
     "iopub.status.idle": "2026-07-14T00:28:13.499066Z",
     "shell.execute_reply": "2026-07-14T00:28:13.497702Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "plt.rcParams['figure.dpi']=110\n",
    "\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "SEED = 26209\n",
    "BETA = 0.9\n",
    "# --- Panel 1: the two-state machine, gentle now risky ---\n",
    "# G: gentle (r=7): stays G w.p. 1-p, drifts to B w.p. p ; hard (r=10): -> B surely\n",
    "# B: repair (r=-4): -> G ; rundown (r=2): -> B\n",
    "def gentle_fp(p):\n",
    "    VG = (7 - 4*BETA*p)/(1 - BETA*(1-p) - BETA**2*p)\n",
    "    return VG, -4 + BETA*VG\n",
    "def hard_fp():\n",
    "    VG = (10 - 4*BETA)/(1 - BETA**2)   # hard -> B surely; repair -> G\n",
    "    return VG, -4 + BETA*VG\n",
    "def policy_in_G(p):\n",
    "    \"\"\"1 = hard optimal, 0 = gentle optimal (optimal V = max over stationary policies).\"\"\"\n",
    "    return int(hard_fp()[0] > gentle_fp(p)[0])\n",
    "def p_flip(grid=None):\n",
    "    grid = grid if grid is not None else [round(0.05*i,2) for i in range(21)]\n",
    "    for p in grid:\n",
    "        if policy_in_G(p): return p\n",
    "    return None\n",
    "# --- Panel 2: stochastic HJB, log-utility Merton ---\n",
    "RHO, R, SIG = 0.10, 0.05, 0.20\n",
    "def A_sig(sig): return (np.log(RHO) + R/RHO - 1 - sig**2/(2*RHO))/RHO\n",
    "def V_sig(x, sig=SIG): return A_sig(sig) + np.log(x)/RHO\n",
    "def resid_sig(x, sig=SIG):\n",
    "    Vp, Vpp = 1/(RHO*x), -1/(RHO*x*x)\n",
    "    c = RHO*x\n",
    "    return RHO*V_sig(x,sig) - (np.log(c) + Vp*(R*x - c) + 0.5*sig**2*x*x*Vpp)\n",
    "# --- Panel 3: chance constraint ---\n",
    "MU, SM, L = 2.0, 0.5, 10.0\n",
    "def i_min(conf): return L/(MU - norm.ppf(conf)*SM)\n",
    "\n",
    "def reference_values():\n",
    "    VG02, VB02 = gentle_fp(0.2)\n",
    "    VGh, _ = hard_fp()\n",
    "    return {\n",
    "        \"VG_p02\": round(VG02,4), \"VB_p02\": round(VB02,4),\n",
    "        \"policy_G_p02\": policy_in_G(0.2),\n",
    "        \"VG_hardFP\": round(VGh,4),\n",
    "        \"policy_G_p06\": policy_in_G(0.6),\n",
    "        \"p_flip\": p_flip(),\n",
    "        \"c_star_sigma\": round(RHO*10,4),\n",
    "        \"value_drop\": round(SIG**2/(2*RHO**2),4),\n",
    "        \"V_sigma_x10\": round(V_sig(10.0),4),\n",
    "        \"resid_sigma_x10\": round(resid_sig(10.0),4),\n",
    "        \"resid_sigma_x25\": round(resid_sig(25.0),4),\n",
    "        \"i_min_95\": round(i_min(0.95),4),\n",
    "        \"i_det\": round(L/MU,4),\n",
    "        \"guarantee_premium\": round(i_min(0.95)-L/MU,4),\n",
    "        \"i_min_99\": round(i_min(0.99),4),\n",
    "    }\n",
    "if __name__ == \"__main__\":\n",
    "    [print(f\"{k:20s} {v}\") for k,v in reference_values().items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43612d59",
   "metadata": {},
   "source": [
    "## Panel 1 — Uncertainty flips a policy\n",
    "Chapter 7's machine, with gentle operation now risky: it preserves the Good\n",
    "state only with probability $1-p$. Both candidate policies have closed-form\n",
    "values; the optimal value is their maximum (Stochastic Enterprise Dynamic\n",
    "Programming Theorem, expectations inside the Bellman backup). At $p = 0.2$\n",
    "gentle still wins ($V_G = 53.22$, down from the deterministic 70 — Enterprise\n",
    "Uncertainty Modifies Optimal Enterprise Policies, Prop.). Sweep $p$: at\n",
    "**$p = 0.55$ the policy flips to hard** — when care no longer reliably\n",
    "preserves the asset, harvesting becomes optimal. Volatility is not just a\n",
    "haircut on value; past a threshold it **rewrites the rule**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45351932",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:28:13.501586Z",
     "iopub.status.busy": "2026-07-14T00:28:13.501200Z",
     "iopub.status.idle": "2026-07-14T00:28:13.734679Z",
     "shell.execute_reply": "2026-07-14T00:28:13.733017Z"
    }
   },
   "outputs": [],
   "source": [
    "ps = [round(0.05*i,2) for i in range(21)]\n",
    "Vg = [gentle_fp(p)[0] for p in ps]\n",
    "Vh = [hard_fp()[0]]*len(ps)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.2))\n",
    "ax.plot(ps, Vg, \"o-\", c=\"#1B6B52\", lw=2.2, ms=4, label=\"gentle policy value (drift risk p)\")\n",
    "ax.plot(ps, Vh, \"--\", c=\"#C8A24B\", lw=2.2, label=\"hard policy value (p-independent)\")\n",
    "ax.axvline(p_flip(), c=\"#B0532F\", ls=\":\", lw=1.5, label=f\"policy flips at p = {p_flip()}\")\n",
    "ax.set(xlabel=\"p — probability gentle operation still drifts to Bad\", ylabel=\"V_G\",\n",
    "       title=\"Care wins until it stops working — seed 26209\")\n",
    "ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"p=0.2: V=({gentle_fp(0.2)[0]:.4f}, {gentle_fp(0.2)[1]:.4f}), gentle optimal: {policy_in_G(0.2)==0}\")\n",
    "print(f\"p=0.6: hard fixed point V_G = {hard_fp()[0]:.4f}, hard optimal: {policy_in_G(0.6)==1}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2fd7039",
   "metadata": {},
   "source": [
    "## Panel 2 — Uncertainty reprices; this policy survives\n",
    "The stochastic HJB (Thm.) adds curvature: $\\rho V = \\max_c[\\ln c +\n",
    "V'(rx-c) + \\tfrac{1}{2}\\sigma^2 x^2 V'']$. With log utility the ansatz\n",
    "$V = A_\\sigma + \\ln(x)/\\rho$ still verifies — **residual identically zero**\n",
    "— and the FOC still gives $c^* = \\rho x$: the feedback law is untouched by\n",
    "$\\sigma$. What moves is the level: $A_\\sigma$ falls by\n",
    "$\\sigma^2/2\\rho^2 = 2.0$ exactly, so $V(10) = -5 - 2 = -7$. The sharpest\n",
    "contrast with Panel 1: there uncertainty rewrote the rule; here it only\n",
    "reprices the enterprise — which of the two happens is a property of the\n",
    "problem, not of uncertainty itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "130a448f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:28:13.736659Z",
     "iopub.status.busy": "2026-07-14T00:28:13.736436Z",
     "iopub.status.idle": "2026-07-14T00:28:14.062666Z",
     "shell.execute_reply": "2026-07-14T00:28:14.061356Z"
    }
   },
   "outputs": [],
   "source": [
    "xs = np.linspace(2, 40, 200)\n",
    "fig, axes = plt.subplots(1,2, figsize=(10,3.9))\n",
    "axes[0].plot(xs, [V_sig(x,0.0) for x in xs], c=\"#8A8F8B\", lw=2, label=\"σ = 0 (Ch. 8)\")\n",
    "axes[0].plot(xs, [V_sig(x,SIG) for x in xs], c=\"#0B3D2E\", lw=2.2, label=\"σ = 0.2\")\n",
    "axes[0].set(xlabel=\"x\", ylabel=\"V(x)\", title=\"Same shape, lower level: the σ²/2ρ² haircut\")\n",
    "axes[0].legend(frameon=False, fontsize=9); axes[0].grid(alpha=.25)\n",
    "axes[1].plot(xs, [resid_sig(x) for x in xs], c=\"#C8A24B\", lw=2.4)\n",
    "axes[1].set(xlabel=\"x\", ylabel=\"stochastic HJB residual\", ylim=(-0.5,0.5), title=\"Verification, with ½σ²x²V″ inside: ≡ 0\")\n",
    "axes[1].grid(alpha=.25)\n",
    "plt.tight_layout(); plt.show()\n",
    "print(f\"c*(10) = {RHO*10:.4f} (unchanged)   V_sigma(10) = {V_sig(10.0):.4f}   drop = {SIG**2/(2*RHO**2):.4f}\")\n",
    "print(f\"residuals: x=10 -> {resid_sig(10.0):.6f}   x=25 -> {resid_sig(25.0):.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f767276",
   "metadata": {},
   "source": [
    "## Panel 3 — Uncertainty as a purchasable guarantee\n",
    "Revenue $R = m \\cdot i$ with $m \\sim N(2, 0.5^2)$; require\n",
    "$P(R \\geq 10) \\geq 1-\\alpha$. The Chance Constraint (Def.) has a\n",
    "deterministic equivalent through the normal quantile:\n",
    "$i \\geq L/(\\mu - z_{1-\\alpha}\\sigma_m)$. Deterministic plan: 5 units. The\n",
    "95% guarantee: **8.49 units — a 3.49 premium**, and 99% costs 11.95. A hard\n",
    "(worst-case) guarantee is impossible under unbounded noise; the chance\n",
    "constraint makes the problem feasible at a price you can read off a quantile\n",
    "(Chance Constraints Enlarge the Class of Feasible Problems, Prop.)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4be1d8ce",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:28:14.064771Z",
     "iopub.status.busy": "2026-07-14T00:28:14.064548Z",
     "iopub.status.idle": "2026-07-14T00:28:14.286645Z",
     "shell.execute_reply": "2026-07-14T00:28:14.285428Z"
    }
   },
   "outputs": [],
   "source": [
    "confs = np.linspace(0.50, 0.995, 200)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.0))\n",
    "ax.plot(confs*100, [i_min(c) for c in confs], c=\"#C8A24B\", lw=2.4)\n",
    "for c,lab in ((0.95,\"95%\"),(0.99,\"99%\")):\n",
    "    ax.scatter([c*100],[i_min(c)], c=\"#0B3D2E\", s=60, zorder=5)\n",
    "    ax.annotate(f\"{lab}: {i_min(c):.2f}\", (c*100, i_min(c)), textcoords=\"offset points\", xytext=(-64,6), fontsize=10, color=\"#0B3D2E\")\n",
    "ax.axhline(L/MU, c=\"#8A8F8B\", ls=\"--\", lw=1.5, label=\"deterministic plan: 5.0\")\n",
    "ax.set(xlabel=\"required confidence (%)\", ylabel=\"minimum investment i\",\n",
    "       title=\"The price of certainty is convex — seed 26209\")\n",
    "ax.legend(frameon=False, fontsize=9); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"i(95%) = {i_min(0.95):.4f}   premium over deterministic: {i_min(0.95)-L/MU:.4f}   i(99%) = {i_min(0.99):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fd10d36",
   "metadata": {},
   "source": [
    "## Validation — agrees with `DCT_V2_Ch09_Lab.xlsx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93299694",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:28:14.288751Z",
     "iopub.status.busy": "2026-07-14T00:28:14.288533Z",
     "iopub.status.idle": "2026-07-14T00:28:14.296742Z",
     "shell.execute_reply": "2026-07-14T00:28:14.295537Z"
    }
   },
   "outputs": [],
   "source": [
    "ref = reference_values()\n",
    "expected = {\"VG_p02\":53.2203,\"VB_p02\":43.8983,\"policy_G_p02\":0,\"VG_hardFP\":33.6842,\n",
    " \"policy_G_p06\":1,\"p_flip\":0.55,\"c_star_sigma\":1.0,\"value_drop\":2.0,\"V_sigma_x10\":-7.0,\n",
    " \"resid_sigma_x10\":0.0,\"resid_sigma_x25\":0.0,\"i_min_95\":8.492,\"i_det\":5.0,\n",
    " \"guarantee_premium\":3.492,\"i_min_99\":11.9499}\n",
    "for k,v in expected.items():\n",
    "    assert abs(ref[k]-v)<5e-4, f\"MISMATCH {k}\"\n",
    "    print(f\"PASS  {k:20s} {ref[k]}\")\n",
    "print(\"\\nAll checkpoints agree — seed 26209.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "656f7b93",
   "metadata": {},
   "source": [
    "**Next**: Exercises 9.5–9.9 (Part C) sweep β against the flip threshold and add drift risk to the Bad state; AXIOM-09's uncertainty console animates the fan of random trajectories. Chapter 10 asks what happens when even the distribution is unknown. Solutions: IM Vol. II, Ch. 9."
   ]
  }
 ],
 "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
}