{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1278bc87",
   "metadata": {},
   "source": [
    "# DCT Laboratory — Volume I, Chapter 14\n",
    "## Mathematical Analysis of the Unified Architecture\n",
    "**Seed `26114`** · Companion to the chapter and AXIOM Module **AXIOM-14**\n",
    "\n",
    "Analysis instruments on Chapter 13's unified $(K, p)$ system: the **Lipschitz\n",
    "constant** behind existence and uniqueness, **convergence** to equilibrium,\n",
    "**perturbation propagation**, and the Sensitivity Theorem's warning made\n",
    "numerical — steady-state sensitivity **explodes** as coupling nears criticality.\n",
    "Mirrored in `DCT_V1_Ch14_Lab.xlsx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d76c07f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T23:17:23.624583Z",
     "iopub.status.busy": "2026-07-13T23:17:23.624329Z",
     "iopub.status.idle": "2026-07-13T23:17:24.695563Z",
     "shell.execute_reply": "2026-07-13T23:17:24.694179Z"
    }
   },
   "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 = 26114\n",
    "AKK, AKP, BK = 0.90, 0.20, 4.0\n",
    "APP, BP = 0.80, 2.0\n",
    "def A_of(g): return np.array([[AKK, AKP],[g, APP]])\n",
    "def K_ss(g):\n",
    "    det = (1-AKK)*(1-APP) - AKP*g\n",
    "    return ((1-APP)*BK + AKP*BP)/det\n",
    "def rho(g): return float(max(abs(np.linalg.eigvals(A_of(g)))))\n",
    "def lip(g):   # spectral norm = sqrt(lambda_max(A^T A))\n",
    "    return float(np.sqrt(max(np.linalg.eigvals(A_of(g).T @ A_of(g)).real)))\n",
    "\n",
    "def perturb_norm(k=12, g=0.06, d0=(0.0, 5.0)):\n",
    "    d = np.array(d0)\n",
    "    for _ in range(k): d = A_of(g) @ d\n",
    "    return float(np.linalg.norm(d))\n",
    "\n",
    "def sens_fd(g, h=0.01):\n",
    "    return (K_ss(g+h) - K_ss(g))/h\n",
    "\n",
    "def reference_values():\n",
    "    return {\n",
    "        \"lip_pre\": round(lip(0.06), 4),\n",
    "        \"rho_pre\": round(rho(0.06), 4),\n",
    "        \"halve_quarters\": round(float(np.log(2)/np.log(1/rho(0.06))), 4),\n",
    "        \"K_ss_006\": round(K_ss(0.06), 4), \"K_ss_007\": round(K_ss(0.07), 4),\n",
    "        \"K_ss_008\": round(K_ss(0.08), 4), \"K_ss_009\": round(K_ss(0.09), 4),\n",
    "        \"sens_fd_low\":  round(sens_fd(0.06), 4),\n",
    "        \"sens_fd_high\": round(sens_fd(0.08), 4),\n",
    "        \"sens_ratio\":   round(sens_fd(0.08)/sens_fd(0.06), 4),\n",
    "        \"pert_norm_12\": round(perturb_norm(), 4),\n",
    "    }\n",
    "if __name__ == \"__main__\":\n",
    "    [print(f\"{k:16s} {v}\") for k,v in reference_values().items()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6ff7744",
   "metadata": {},
   "source": [
    "## Panel 1 — Existence, uniqueness, convergence\n",
    "The Existence and Uniqueness Theorems ask for a Lipschitz bound: for the linear\n",
    "system it is the spectral norm, $L = 0.9922 < 1$ — a contraction, so one solution\n",
    "exists per start and iterates converge. The rate is $\\rho = 0.9704$: the error\n",
    "halves every **23.08 quarters**. Well-posed and *slow* — both facts matter, and\n",
    "they are different facts ($L$ vs. $\\rho$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "866f1b1d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T23:17:24.698608Z",
     "iopub.status.busy": "2026-07-13T23:17:24.697660Z",
     "iopub.status.idle": "2026-07-13T23:17:24.722679Z",
     "shell.execute_reply": "2026-07-13T23:17:24.721393Z"
    }
   },
   "outputs": [],
   "source": [
    "g = 0.06\n",
    "print(f\"Lipschitz (spectral norm) L = {lip(g):.4f}   spectral radius ρ = {rho(g):.4f}\")\n",
    "print(f\"error half-life: ln2/ln(1/ρ) = {np.log(2)/np.log(1/rho(g)):.4f} quarters\")\n",
    "z = np.array([40.0, 50.0]); ss = np.array([K_ss(g), 0]);  # track K error only for display\n",
    "errs = []\n",
    "x = np.array([40.0, 50.0])\n",
    "target = None\n",
    "# iterate far to find the fixed point numerically, then measure error decay\n",
    "for _ in range(4000): x = A_of(g)@x + np.array([BK, BP])\n",
    "fp = x.copy()\n",
    "x = np.array([40.0, 50.0])\n",
    "for k in range(60):\n",
    "    errs.append(float(np.linalg.norm(x-fp))); x = A_of(g)@x + np.array([BK, BP])\n",
    "r = errs[40]/errs[39]\n",
    "print(f\"measured per-step error ratio at k=40: {r:.4f}  (→ ρ)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3f09689",
   "metadata": {},
   "source": [
    "## Panel 2 — Sensitivity explodes near criticality\n",
    "Steady-state capital as a function of the coupling $g$: 150 at 0.06, 200 at\n",
    "0.07, 300 at 0.08, **600 at 0.09**. The finite-difference sensitivity rises from\n",
    "5,000 to 30,000 over the same interval — a factor of 6. Sensitivity Depends on\n",
    "Coupling Strength (Prop.): near criticality, small parameter errors become\n",
    "large forecast errors, which is the Sensitivity Theorem's operational content."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "913a93d8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T23:17:24.724742Z",
     "iopub.status.busy": "2026-07-13T23:17:24.724277Z",
     "iopub.status.idle": "2026-07-13T23:17:25.028211Z",
     "shell.execute_reply": "2026-07-13T23:17:25.026461Z"
    }
   },
   "outputs": [],
   "source": [
    "gs = np.linspace(0.05, 0.092, 200)\n",
    "Ks = [K_ss(g) for g in gs]\n",
    "fig, ax = plt.subplots(figsize=(8.2,4.4))\n",
    "ax.plot(gs, Ks, c=\"#C8A24B\", lw=2.5)\n",
    "for g in (0.06,0.07,0.08,0.09):\n",
    "    ax.scatter([g],[K_ss(g)], c=\"#0B3D2E\", zorder=5)\n",
    "    ax.annotate(f\"{K_ss(g):.0f}\", (g, K_ss(g)), textcoords=\"offset points\", xytext=(6,6), fontsize=9)\n",
    "ax.set(xlabel=\"coupling g\", ylabel=\"steady-state capital K*\", \n",
    "       title=\"K*(g) = ((1−a_pp)b_K + a_Kp b_p) / det(g) — seed 26114\")\n",
    "ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"sensitivity (fd, h=0.01): at g=0.06 → {sens_fd(0.06):.1f}   at g=0.08 → {sens_fd(0.08):.1f}   ratio {sens_fd(0.08)/sens_fd(0.06):.1f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e30498e",
   "metadata": {},
   "source": [
    "## Panel 3 — Perturbations propagate\n",
    "A +5 measurement error in performance at $t=0$. The difference system\n",
    "$\\delta_{k+1} = A\\delta_k$ carries it through the coupling: after 12 quarters\n",
    "the error has migrated into capital and decayed to norm **3.01** — Perturbations\n",
    "Propagate Through Architectural Dependencies (Prop.), traced coordinate by\n",
    "coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2317679b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T23:17:25.030658Z",
     "iopub.status.busy": "2026-07-13T23:17:25.030314Z",
     "iopub.status.idle": "2026-07-13T23:17:25.235076Z",
     "shell.execute_reply": "2026-07-13T23:17:25.233668Z"
    }
   },
   "outputs": [],
   "source": [
    "d = np.array([0.0, 5.0]); path = [d.copy()]\n",
    "for _ in range(12): d = A_of(0.06)@d; path.append(d.copy())\n",
    "path = np.array(path)\n",
    "fig, ax = plt.subplots(figsize=(7.8,4.0))\n",
    "ax.plot(path[:,0], \"o-\", c=\"#0B3D2E\", lw=2, ms=4, label=\"δK (was 0)\")\n",
    "ax.plot(path[:,1], \"s-\", c=\"#C8A24B\", lw=2, ms=4, label=\"δp (was 5)\")\n",
    "ax.set(xlabel=\"quarter\", ylabel=\"perturbation\", title=\"A +5 error in p migrates into K\")\n",
    "ax.legend(frameon=False); ax.grid(alpha=.25); plt.tight_layout(); plt.show()\n",
    "print(f\"‖δ‖ at k=12: {np.linalg.norm(path[12]):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98035680",
   "metadata": {},
   "source": [
    "## Validation — agrees with `DCT_V1_Ch14_Lab.xlsx`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df3b97c9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T23:17:25.237539Z",
     "iopub.status.busy": "2026-07-13T23:17:25.237314Z",
     "iopub.status.idle": "2026-07-13T23:17:25.244523Z",
     "shell.execute_reply": "2026-07-13T23:17:25.243431Z"
    }
   },
   "outputs": [],
   "source": [
    "ref = reference_values()\n",
    "expected = {\"lip_pre\":0.9922,\"rho_pre\":0.9704,\"halve_quarters\":23.0814,\n",
    " \"K_ss_006\":150.0,\"K_ss_007\":200.0,\"K_ss_008\":300.0,\"K_ss_009\":600.0,\n",
    " \"sens_fd_low\":5000.0,\"sens_fd_high\":30000.0,\"sens_ratio\":6.0,\"pert_norm_12\":3.0097}\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 26114.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd554db8",
   "metadata": {},
   "source": [
    "**Next**: Exercises 14.9–14.12 (Part C) push g toward det = 0; AXIOM-14's analysis bench plots every instrument against the coupling live. Solutions: IM Ch. 14."
   ]
  }
 ],
 "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
}