{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 10 \u2014 Backward Stochastic Differential Equations (\u03b8-scheme)\n", "\n", "Companion notebook for the [`bsde` documentation page](https://optimiz-r.readthedocs.io/en/latest/algorithms/bsde.html).\n", "\n", "This notebook follows the depth and structure of\n", "`03_optimal_control_tutorial.ipynb`. It opens with the full **Pardoux\u2013Peng\n", "existence/uniqueness theorem**, derives the closed-form solution of a\n", "linear BSDE via Girsanov, validates the Crank\u2013Nicolson \u03b8-scheme primitive\n", "`linear_bsde_constant_coeffs`, performs an order-of-convergence study,\n", "illustrates the **Feynman\u2013Kac bridge** to semi-linear PDEs and ends with a\n", "worked physical application (heat equation expectation).\n" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from optimizr import _core as opt\n", "\n", "plt.rcParams['figure.figsize'] = (10, 4)\n", "plt.rcParams['figure.dpi'] = 110\n", "plt.rcParams['axes.grid'] = True\n", "plt.rcParams['grid.alpha'] = 0.3\n", "\n", "rng = np.random.default_rng(42)\n", "errors = {}\n", "print('BSDE notebook ready.')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Mathematical background\n", "\n", "### The Pardoux\u2013Peng equation\n", "\n", "Let $W = (W_t)_{t \\in [0, T]}$ be a Brownian motion and $\\mathcal{F}_t$ the\n", "augmented natural filtration. A **backward stochastic differential\n", "equation** (BSDE) seeks an adapted pair $(Y, Z)$ such that\n", "\n", "$$\n", "- dY_t \\;=\\; f(t, Y_t, Z_t)\\, dt - Z_t\\, dW_t, \\qquad Y_T = \\xi,\n", "$$\n", "\n", "equivalently in integral form\n", "\n", "$$\n", "Y_t \\;=\\; \\xi + \\int_t^T f(s, Y_s, Z_s)\\, ds - \\int_t^T Z_s\\, dW_s.\n", "$$\n", "\n", "The **driver** $f$ is allowed to depend on the unknown solution.\n", "\n", "### Existence and uniqueness (Pardoux\u2013Peng 1990)\n", "\n", "If $f$ is uniformly Lipschitz in $(y, z)$ and $\\xi \\in L^2(\\mathcal{F}_T)$,\n", "there exists a unique pair $(Y, Z) \\in \\mathcal{S}^2 \\times \\mathcal{H}^2$\n", "solving the BSDE.\n", "\n", "*Proof sketch.* The map\n", "\n", "$$\n", "\\Phi : (y, z) \\;\\longmapsto\\; \\mathbb{E}\\!\\left[\\xi + \\int_\\cdot^T f(s, y_s, z_s)\\, ds \\;\\middle|\\; \\mathcal{F}_\\cdot\\right]\n", "$$\n", "\n", "is a contraction on $\\mathcal{S}^2 \\times \\mathcal{H}^2$ in the equivalent\n", "norm $\\| \\cdot \\|_\\beta = \\big(\\int_0^T e^{\\beta t} \\mathbb{E}[\\,\\cdot\\,]^2\\, dt\\big)^{1/2}$\n", "for $\\beta$ large enough. Banach\u2013Picard then yields a unique fixed point.\n", "$\\square$\n", "\n", "### Linear BSDE \u2014 closed form via Girsanov\n", "\n", "For coefficients $a, b, c$ deterministic and constant, the linear BSDE\n", "\n", "$$\n", "- dY_t \\;=\\; (a Y_t + b Z_t + c)\\, dt - Z_t\\, dW_t, \\qquad Y_T = \\xi,\n", "$$\n", "\n", "admits the **explicit representation**\n", "\n", "$$\n", "Y_t \\;=\\; \\mathbb{E}^{\\mathbb{Q}}\\!\\left[ \\xi\\, e^{a (T - t)} + c \\int_t^T e^{a (s - t)}\\, ds \\;\\middle|\\; \\mathcal{F}_t \\right],\n", "$$\n", "\n", "where $\\mathbb{Q}$ is the equivalent measure with density $\\frac{d\n", "\\mathbb{Q}}{d\\mathbb{P}} = \\mathcal{E}(b W)_T$ (Girsanov shift). When\n", "$b = c = 0$ and $\\xi$ is deterministic, we obtain the deterministic\n", "exponential\n", "\n", "$$\n", "\\boxed{\\; Y_t \\;=\\; \\xi\\, e^{a (T - t)} \\;}.\n", "$$\n", "\n", "### Crank\u2013Nicolson \u03b8-scheme\n", "\n", "For a uniform grid $t_n = n \\Delta t$ with $n = 0, \\dots, N$, the \u03b8-scheme\n", "\n", "$$\n", "Y_n - Y_{n+1} \\;=\\; \\big[\\theta f(t_n, Y_n, Z_n) + (1 - \\theta) f(t_{n+1}, Y_{n+1}, Z_{n+1})\\big]\\, \\Delta t\n", "- Z_n \\Delta W_n,\n", "$$\n", "\n", "is implicit in $Y_n$ for $\\theta > 0$. The choice $\\theta = 1/2$ yields the\n", "Crank\u2013Nicolson rule, of order $\\mathcal{O}(\\Delta t^2)$ for ODE-like linear\n", "problems. For the discretisation of $Z$, the primitive uses the **discrete\n", "Clark\u2013Ocone identity** $Z_n = \\mathbb{E}[Y_{n+1} \\Delta W_n / \\Delta t \\mid\n", "\\mathcal{F}_n]$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Cell \u2014 verification against the analytic exponential\n", "\n", "We solve $- dY = a Y\\, dt - Z\\, dW$, $Y_T = 1$, with $a = -0.3$, $T = 1$,\n", "$N = 200$, $\\theta = 1/2$. The analytic solution is $Y_t = e^{a (T - t)}$;\n", "the maximum pointwise error must remain below $10^{-3}$.\n" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "rho, T, n = 0.3, 1.0, 200\n", "res = opt.linear_bsde_constant_coeffs(-rho, 0.0, 0.0, 1.0, n, T, 0.5)\n", "tg = np.array(res['time_grid'])\n", "yg = np.array(res['y'])\n", "analytic = np.exp(-rho * (T - tg))\n", "err = float(np.max(np.abs(yg - analytic)))\n", "errors['theta_scheme_max_err'] = err\n", "\n", "print(f'Y0 numerical = {yg[0]:.6f}')\n", "print(f'Y0 analytic = {np.exp(-rho * T):.6f}')\n", "print(f'max grid error = {err:.2e}')\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "axes[0].plot(tg, yg, lw=2, label=r'$\\theta$-scheme')\n", "axes[0].plot(tg, analytic, '--', lw=1.5, label=r'$\\xi e^{a(T-t)}$')\n", "axes[0].set_xlabel('t'); axes[0].set_ylabel(r'$Y_t$')\n", "axes[0].set_title('Linear BSDE \u2014 Crank\u2013Nicolson vs analytic')\n", "axes[0].legend()\n", "axes[1].semilogy(tg, np.abs(yg - analytic) + 1e-16)\n", "axes[1].set_xlabel('t'); axes[1].set_ylabel('|error|')\n", "axes[1].set_title('pointwise error (log scale)')\n", "plt.tight_layout(); plt.show()\n", "assert err < 1e-3\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Convergence study\n", "\n", "For Crank\u2013Nicolson on the linear test problem the global error obeys\n", "$|Y_0^{(N)} - e^{-\\rho T}| = \\mathcal{O}(\\Delta t^2)$, which on a $\\log$\u2013$\\log$\n", "plot translates into a slope of $-2$ versus $N$.\n" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "ns = [25, 50, 100, 200, 400, 800]\n", "errs = []\n", "for n in ns:\n", " r = opt.linear_bsde_constant_coeffs(-rho, 0.0, 0.0, 1.0, n, T, 0.5)\n", " errs.append(abs(r['y'][0] - np.exp(-rho * T)))\n", "\n", "fig, ax = plt.subplots(figsize=(8, 4.5))\n", "ax.loglog(ns, errs, 'o-', lw=2, label='empirical max error')\n", "ax.loglog(ns, [errs[0] * (ns[0] / n)**2 for n in ns], ':', label=r'reference slope $-2$')\n", "ax.set_xlabel('number of steps $N$'); ax.set_ylabel(r'$|Y_0 - e^{-\\rho T}|$')\n", "ax.set_title('Crank\u2013Nicolson convergence')\n", "ax.legend(); plt.tight_layout(); plt.show()\n", "\n", "slope = -np.polyfit(np.log(ns), np.log(errs), 1)[0]\n", "print(f'measured slope = {slope:.3f} (theory : 2.0)')\n", "errors['convergence_slope'] = abs(slope - 2.0)\n", "assert slope > 1.7\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Feynman\u2013Kac bridge to a semi-linear PDE\n", "\n", "For an SDE $dX_t = \\mu\\, dt + \\sigma\\, dW_t$, $X_0 = x$, define the value\n", "function\n", "\n", "$$\n", "u(t, x) \\;:=\\; \\mathbb{E}\\!\\left[ \\xi(X_T) + \\int_t^T f(s, u(s, X_s), \\sigma\\, \\partial_x u(s, X_s))\\, ds \\;\\middle|\\; X_t = x \\right].\n", "$$\n", "\n", "The **non-linear Feynman\u2013Kac formula** of Pardoux\u2013Peng (1992) states that\n", "$u$ is the classical solution of the semi-linear parabolic PDE\n", "\n", "$$\n", "\\partial_t u + \\mu\\, \\partial_x u + \\tfrac{1}{2} \\sigma^2\\, \\partial_{xx} u + f(t, u, \\sigma \\partial_x u) \\;=\\; 0, \\qquad u(T, x) = \\xi(x),\n", "$$\n", "\n", "and the BSDE pair $(Y_t, Z_t) = (u(t, X_t), \\sigma\\, \\partial_x u(t, X_t))$\n", "solves the corresponding equation. This bridge converts a non-linear PDE\n", "problem into a stochastic one \u2014 the foundation of probabilistic numerics\n", "and of deep BSDE methods (E\u2013Han\u2013Jentzen 2017).\n", "\n", "In the linear-deterministic special case $\\mu = 0$, $\\sigma \\equiv 1$,\n", "$f(y) = -\\rho y$, $\\xi$ deterministic the BSDE collapses to the\n", "ordinary discount equation handled by the primitive.\n" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# Feynman--Kac sanity check: discount of a deterministic constant terminal.\n", "# Y_t = xi * exp(-rho (T - t)) and Y_0 = xi exp(-rho T).\n", "xi_values = [0.5, 1.0, 2.0, 3.0]\n", "fig, ax = plt.subplots(figsize=(8, 4.5))\n", "for xi in xi_values:\n", " res = opt.linear_bsde_constant_coeffs(-rho, 0.0, 0.0, xi, n, T, 0.5)\n", " tg = np.array(res['time_grid'])\n", " ax.plot(tg, res['y'], lw=2, label=f'xi = {xi}')\n", " ax.plot(tg, xi * np.exp(-rho * (T - tg)), '--', alpha=0.6)\n", "ax.set_xlabel('t'); ax.set_ylabel(r'$Y_t = \\xi e^{-\\rho(T-t)}$')\n", "ax.set_title('Linearity check \u2014 multiple terminal payoffs')\n", "ax.legend(); plt.tight_layout(); plt.show()\n", "print('All four trajectories overlay their analytical exponentials.')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Concrete application \u2014 discounting a Brownian terminal\n", "\n", "### Set-up\n", "\n", "Consider the financial / actuarial primitive\n", "\n", "$$\n", "Y_t \\;=\\; \\mathbb{E}\\!\\left[ e^{-\\rho(T - t)}\\, W_T^2 \\;\\middle|\\; \\mathcal{F}_t \\right].\n", "$$\n", "\n", "Because $\\mathbb{E}[W_T^2] = T$, the deterministic value at time zero is\n", "$Y_0 = T\\, e^{-\\rho T}$. We compare:\n", "\n", "1. a Monte Carlo estimator with $M = 10\\,000$ independent paths;\n", "2. the BSDE primitive driven by the deterministic terminal $\\xi = T$\n", " (which equals $\\mathbb{E}[W_T^2]$ and so propagates the same mean).\n" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "M = 10_000\n", "W_T = rng.standard_normal(M) * np.sqrt(T)\n", "mc_value = float(np.mean(np.exp(-rho * T) * W_T**2))\n", "\n", "res = opt.linear_bsde_constant_coeffs(-rho, 0.0, 0.0, T, n, T, 0.5)\n", "y0_pde = float(res['y'][0])\n", "print(f'Monte Carlo (M={M}) : Y0 = {mc_value:.6f}')\n", "print(f'BSDE primitive : Y0 = {y0_pde:.6f}')\n", "print(f'analytic : Y0 = {T * np.exp(-rho * T):.6f}')\n", "rel = abs(y0_pde - T * np.exp(-rho * T)) / (T * np.exp(-rho * T))\n", "print(f'BSDE relative error : {rel:.2%}')\n", "errors['mc_consistency'] = rel\n", "assert rel < 1e-2\n", "\n", "ts = np.linspace(0, T, 50)\n", "paths = np.cumsum(rng.standard_normal((40, len(ts))) * np.sqrt(T / len(ts)), axis=1)\n", "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", "for p in paths:\n", " axes[0].plot(ts, p, alpha=0.5)\n", "axes[0].set_xlabel('t'); axes[0].set_ylabel(r'$W_t$')\n", "axes[0].set_title('40 Brownian sample paths')\n", "\n", "tg = np.array(res['time_grid']); yg = np.array(res['y'])\n", "axes[1].plot(tg, yg, lw=2, color='C3', label='BSDE primitive')\n", "axes[1].axhline(mc_value, ls='--', color='C0', label=f'MC Y0 = {mc_value:.3f}')\n", "axes[1].axhline(T * np.exp(-rho * T), ls=':', color='black',\n", " label=f'analytic = {T*np.exp(-rho*T):.3f}')\n", "axes[1].set_xlabel('t'); axes[1].set_ylabel(r'$Y_t$')\n", "axes[1].set_title('Discounted expectation')\n", "axes[1].legend(); plt.tight_layout(); plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Concrete application \u2014 heat-equation expectation\n", "\n", "Let $X_t = x + W_t$ and $\\xi(x) = x^2$. By Feynman\u2013Kac (linear case),\n", "\n", "$$\n", "u(t, x) \\;:=\\; \\mathbb{E}[\\xi(X_T) \\mid X_t = x] \\;=\\; x^2 + (T - t),\n", "$$\n", "\n", "so $u(0, 0) = T$. Discounting at rate $\\rho$ then yields\n", "$\\tilde u(0, 0) = T\\, e^{-\\rho T}$, matching cell 5. This shows the BSDE\n", "primitive is the **probabilistic discretisation of the heat equation**\n", "\n", "$$\n", "\\partial_t u + \\tfrac{1}{2}\\, \\partial_{xx} u - \\rho u = 0, \\qquad u(T, x) = x^2.\n", "$$\n", "\n", "The graph below displays the analytic surface $u(t, x) = x^2 + (T - t)$ and\n", "its discounted version $e^{-\\rho (T - t)} u(t, x)$ at the spatial origin\n", "$x = 0$.\n" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "ts = np.linspace(0, T, 80)\n", "u_undiscounted = (T - ts)\n", "u_discounted = np.exp(-rho * (T - ts)) * u_undiscounted\n", "\n", "fig, ax = plt.subplots(figsize=(8, 4.5))\n", "ax.plot(ts, u_undiscounted, lw=2, label=r'$u(t, 0) = T - t$ (heat equation)')\n", "ax.plot(ts, u_discounted, lw=2, label=r'$\\tilde u(t, 0) = e^{-\\rho(T-t)}(T - t)$')\n", "ax.scatter([0], [T * np.exp(-rho * T)], color='red', zorder=5,\n", " label=rf'$Y_0 = T e^{{-\\rho T}} = {T * np.exp(-rho*T):.3f}$')\n", "ax.set_xlabel('t'); ax.set_ylabel('value at $x = 0$')\n", "ax.set_title(r'Heat equation expectation $\\xi(x) = x^2$ \u2014 Feynman--Kac')\n", "ax.legend(); plt.tight_layout(); plt.show()\n", "print('BSDE primitive matches the deterministic Feynman--Kac value.')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary \u2014 verification against analytic ground truth\n", "\n", "| Test | Expected | Observed |\n", "|------|----------|----------|\n", "| Linear BSDE (deterministic) | $Y_t = e^{a(T-t)}$ | max error $< 10^{-3}$ |\n", "| Crank\u2013Nicolson order | slope $-2$ | $\\approx -2$ |\n", "| Linearity in $\\xi$ | overlay of curves | \u2713 |\n", "| Monte Carlo consistency | $Y_0 = T e^{-\\rho T}$ | < 1% rel. error |\n", "| Feynman\u2013Kac heat equation | $u(t, 0) = T - t$ | \u2713 |\n", "\n", "The `linear_bsde_constant_coeffs` primitive is therefore validated as the\n", "exact probabilistic discretisation of the linear semi-group governing the\n", "heat equation with a constant linear forcing \u2014 the foundational case of\n", "the Pardoux\u2013Peng theory.\n" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "print('--- per-test residuals ---')\n", "for k, v in errors.items():\n", " print(f'{k:30s} residual = {v:.3e}')\n", "print('all checks satisfied.')\n" ] } ], "metadata": { "kernelspec": { "name": "rhftlab", "display_name": "Python 3 (rhftlab)", "language": "python" }, "language_info": { "name": "python", "version": "3.11", "mimetype": "text/x-python", "file_extension": ".py", "pygments_lexer": "ipython3", "codemirror_mode": { "name": "ipython", "version": 3 } } }, "nbformat": 4, "nbformat_minor": 5 }