{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d8af038e",
   "metadata": {},
   "source": [
    "# DCT Laboratory — Volume II, Chapter 8\n",
    "## Hamilton–Jacobi–Bellman Enterprise Framework\n",
    "**Seed `26208`** · Companion to the chapter and AXIOM Module **AXIOM-08 (Vol. II)**\n",
    "\n",
    "HJB with closed forms. The **stationary consumption problem**: value function\n",
    "$V(x) = A + \\ln(x)/\\rho$ verifies the HJB equation with **residual identically\n",
    "zero** — the Verification Theorem executed — and the feedback law is one line:\n",
    "$c^*(x) = \\rho x$. The **closed loop** simulated Euler-vs-exact. And the\n",
    "**finite-horizon LQ Riccati equation** solved twice: closed form and backward\n",
    "Euler, with the discretization error measured honestly.\n",
    "Mirrored in `DCT_V2_Ch08_Lab.xlsx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f68729c7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:11:59.480761Z",
     "iopub.status.busy": "2026-07-14T00:11:59.480580Z",
     "iopub.status.idle": "2026-07-14T00:11:59.985208Z",
     "shell.execute_reply": "2026-07-14T00:11:59.983395Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import log, exp\n",
    "plt.rcParams['figure.dpi']=110\n",
    "\n",
    "import numpy as np\n",
    "from math import log, exp\n",
    "SEED = 26208\n",
    "# --- stationary HJB: max int e^{-rho t} ln(c), dx = (r x - c) dt ---\n",
    "RHO, R = 0.10, 0.05\n",
    "A_CONST = (log(RHO) + R/RHO - 1)/RHO\n",
    "def V(x): return A_CONST + log(x)/RHO\n",
    "def c_star(x): return RHO*x\n",
    "def hjb_residual(x):\n",
    "    Vp = 1/(RHO*x)\n",
    "    return RHO*V(x) - (log(c_star(x)) + Vp*(R*x - c_star(x)))\n",
    "# closed loop: dx = (R - RHO) x dt\n",
    "X0, T, H = 10.0, 10.0, 0.1\n",
    "def euler_x():\n",
    "    x = X0\n",
    "    for _ in range(int(T/H)): x += H*(R-RHO)*x\n",
    "    return x\n",
    "# --- finite-horizon LQ: min int u^2/2 + q x(T)^2/2, dx = u dt ---\n",
    "QT, TL = 2.0, 4.0\n",
    "def p_exact(t): return QT/(1+QT*(TL-t))\n",
    "def p_euler(h=0.05):\n",
    "    p = QT; err = 0.0\n",
    "    t = TL\n",
    "    while t > 1e-9:\n",
    "        err = max(err, abs(p - p_exact(t)))\n",
    "        p = p - h*(p*p)    # dp/dt = p^2; stepping t -> t-h subtracts h*p^2\n",
    "        t -= h\n",
    "    return p, max(err, abs(p - p_exact(0.0)))\n",
    "\n",
    "def reference_values():\n",
    "    pe, perr = p_euler()\n",
    "    return {\n",
    "        \"c_star_x10\": round(c_star(10.0),4),\n",
    "        \"V_x10\": round(V(10.0),4),\n",
    "        \"hjb_residual_x10\": round(hjb_residual(10.0),4),\n",
    "        \"hjb_residual_x25\": round(hjb_residual(25.0),4),\n",
    "        \"x_T_exact\": round(X0*exp((R-RHO)*T),4),\n",
    "        \"x_T_euler\": round(euler_x(),4),\n",
    "        \"euler_gap\": round(abs(X0*exp((R-RHO)*T)-euler_x()),4),\n",
    "        \"p0_exact\": round(p_exact(0.0),4),\n",
    "        \"p0_euler\": round(pe,4),\n",
    "        \"riccati_maxerr\": round(perr,4),\n",
    "        \"V_LQ_x3\": round(p_exact(0.0)*9/2,4),\n",
    "        \"xT_LQ_x3\": round(3.0/(1+QT*TL),4),\n",
    "    }\n",
    "if __name__ == \"__main__\":\n",
    "    [print(f\"{k:18s} {v}\") for k,v in reference_values().items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6081005",
   "metadata": {},
   "source": [
    "## Panel 1 — Verification: the residual is zero, everywhere\n",
    "$\\max_c \\int e^{-\\rho t}\\ln c\\,dt$, $dx = (rx - c)dt$, $\\rho = 0.10$,\n",
    "$r = 0.05$. Candidate: $V(x) = A + \\ln(x)/\\rho$. The stationary HJB\n",
    "$\\rho V = \\max_c[\\ln c + V'(x)(rx - c)]$: the inner maximization gives\n",
    "$c^* = \\rho x$, and substituting back the equation holds **identically** —\n",
    "residual $0.0000$ at every $x$. That is the Verification Theorem's logic: a\n",
    "smooth solution of HJB that is attained IS the value function (The HJB\n",
    "Framework Provides Sufficient Conditions, Prop.) — sufficiency, where\n",
    "Pontryagin gave necessity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e58b0b1d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:11:59.988178Z",
     "iopub.status.busy": "2026-07-14T00:11:59.987263Z",
     "iopub.status.idle": "2026-07-14T00:12:00.290599Z",
     "shell.execute_reply": "2026-07-14T00:12:00.289470Z"
    }
   },
   "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(x) for x in xs], c=\"#0B3D2E\", lw=2.4)\n",
    "axes[0].set(xlabel=\"x\", ylabel=\"V(x)\", title=\"The value surface: V = A + ln(x)/ρ\")\n",
    "axes[0].grid(alpha=.25)\n",
    "axes[1].plot(xs, [hjb_residual(x) for x in xs], c=\"#C8A24B\", lw=2.4)\n",
    "axes[1].set(xlabel=\"x\", ylabel=\"HJB residual\", ylim=(-0.5,0.5), title=\"ρV − max_c[...] ≡ 0\")\n",
    "axes[1].grid(alpha=.25)\n",
    "plt.tight_layout(); plt.show()\n",
    "print(f\"c*(10) = {c_star(10.0):.4f}   V(10) = {V(10.0):.4f}\")\n",
    "print(f\"residual at x=10: {hjb_residual(10.0):.6f}   at x=25: {hjb_residual(25.0):.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31df1ac8",
   "metadata": {},
   "source": [
    "## Panel 2 — The feedback law, in closed loop\n",
    "$c^*(x) = \\rho x$: consume a fixed fraction of capital, whatever happens — the\n",
    "Optimal Enterprise Feedback Law (Def.) as one line of policy. The closed loop\n",
    "$dx = (r-\\rho)x\\,dt$ decays at 5%/yr: exact $x(10) = 6.0653$; Euler with\n",
    "$h = 0.1$ gives $6.0577$ — a 0.0076 gap that is pure discretization, measured\n",
    "because the closed form exists to measure against."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5e67156",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:12:00.292512Z",
     "iopub.status.busy": "2026-07-14T00:12:00.292311Z",
     "iopub.status.idle": "2026-07-14T00:12:00.475036Z",
     "shell.execute_reply": "2026-07-14T00:12:00.473800Z"
    }
   },
   "outputs": [],
   "source": [
    "ts = np.linspace(0, T, 200)\n",
    "exact = X0*np.exp((R-RHO)*ts)\n",
    "xe, path = X0, [X0]\n",
    "for _ in range(int(T/H)): xe += H*(R-RHO)*xe; path.append(xe)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.0))\n",
    "ax.plot(ts, exact, c=\"#0B3D2E\", lw=2.2, label=f\"exact → x(10) = {X0*exp((R-RHO)*T):.4f}\")\n",
    "ax.plot(np.arange(len(path))*H, path, \"o\", c=\"#C8A24B\", ms=3, label=f\"Euler h=0.1 → {euler_x():.4f}\")\n",
    "ax.set(xlabel=\"t\", ylabel=\"x(t)\", title=\"Closed loop under c* = ρx (seed 26208)\")\n",
    "ax.legend(frameon=False); ax.grid(alpha=.25); plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6d04737",
   "metadata": {},
   "source": [
    "## Panel 3 — Time-dependent HJB: the Riccati equation\n",
    "Finite-horizon LQ: $\\min \\int_0^4 u^2/2\\,dt + q\\,x(4)^2/2$, $\\dot{x} = u$,\n",
    "$q = 2$. The quadratic ansatz $V(x,t) = p(t)x^2/2$ reduces HJB to the Riccati\n",
    "ODE $\\dot{p} = p^2$, $p(4) = 2$ — closed form $p(t) = 2/(1+2(4-t))$, so\n",
    "$p(0) = 0.2222$ and $V(3, 0) = 1.0$. Backward Euler ($h = 0.05$) lands at\n",
    "$0.2168$: max error $0.0397$, concentrated where $\\dot{p}$ is steep — Numerical\n",
    "Approximation Becomes Necessary (Prop.), and this is what its error looks like\n",
    "when a closed form is standing by to grade it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15ad9c1f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:12:00.477177Z",
     "iopub.status.busy": "2026-07-14T00:12:00.476972Z",
     "iopub.status.idle": "2026-07-14T00:12:00.677063Z",
     "shell.execute_reply": "2026-07-14T00:12:00.676046Z"
    }
   },
   "outputs": [],
   "source": [
    "ts = np.linspace(0, TL, 200)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.0))\n",
    "ax.plot(ts, [p_exact(t) for t in ts], c=\"#0B3D2E\", lw=2.2, label=\"p(t) = 2/(1+2(4−t)) exact\")\n",
    "# backward Euler trace\n",
    "h=0.05; p=QT; tt=TL; te=[TL]; pe=[QT]\n",
    "while tt > 1e-9:\n",
    "    p = p - h*p*p; tt -= h; te.append(tt); pe.append(p)\n",
    "ax.plot(te, pe, \"o\", c=\"#C8A24B\", ms=2.5, label=f\"backward Euler h=0.05 → p(0) = {pe[-1]:.4f}\")\n",
    "ax.set(xlabel=\"t\", ylabel=\"p(t)\", title=\"The Riccati gain: exact vs numerical\")\n",
    "ax.legend(frameon=False); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "pe0, perr = p_euler()\n",
    "print(f\"p(0) exact {p_exact(0.0):.4f}   Euler {pe0:.4f}   max error {perr:.4f}\")\n",
    "print(f\"V(x=3, t=0) = p(0)*9/2 = {p_exact(0.0)*4.5:.4f}   closed-loop x(4) from x0=3: {3.0/(1+QT*TL):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5826fc0",
   "metadata": {},
   "source": [
    "## Validation — agrees with `DCT_V2_Ch08_Lab.xlsx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a8ba89b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-14T00:12:00.679059Z",
     "iopub.status.busy": "2026-07-14T00:12:00.678845Z",
     "iopub.status.idle": "2026-07-14T00:12:00.685569Z",
     "shell.execute_reply": "2026-07-14T00:12:00.684503Z"
    }
   },
   "outputs": [],
   "source": [
    "ref = reference_values()\n",
    "expected = {\"c_star_x10\":1.0,\"V_x10\":-5.0,\"hjb_residual_x10\":0.0,\"hjb_residual_x25\":0.0,\n",
    " \"x_T_exact\":6.0653,\"x_T_euler\":6.0577,\"euler_gap\":0.0076,\n",
    " \"p0_exact\":0.2222,\"p0_euler\":0.2168,\"riccati_maxerr\":0.0397,\n",
    " \"V_LQ_x3\":1.0,\"xT_LQ_x3\":0.3333}\n",
    "for k,v in expected.items():\n",
    "    assert abs(ref[k]-v)<5e-4, f\"MISMATCH {k}\"\n",
    "    print(f\"PASS  {k:18s} {ref[k]}\")\n",
    "print(\"\\nAll checkpoints agree — seed 26208.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8894c86",
   "metadata": {},
   "source": [
    "**Next**: Exercises 8.5–8.9 (Part C) shrink h and watch the Riccati error fall at first order; AXIOM-08's value-surface explorer renders V(x,t) live. Chapter 9 adds Volume I's randomness: stochastic optimization. Solutions: IM Vol. II, Ch. 8."
   ]
  }
 ],
 "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
}