diff --git a/docs/source/algorithms/agent_based.rst b/docs/source/algorithms/agent_based.rst index 8414b90..80b5c62 100644 --- a/docs/source/algorithms/agent_based.rst +++ b/docs/source/algorithms/agent_based.rst @@ -1,9 +1,78 @@ Agent-based — bounded-confidence consensus ========================================== -Generic interacting-agent simulator (`consensus_dynamics`) — linear bounded-confidence rule $s_i^{k+1} = (1-α) s_i^k + α \bar s^k + ξ_i$. +Generic symmetric *interacting-agent* simulator implementing the linear bounded-confidence +update rule -.. note:: Companion executed notebook: `15_agent_based.ipynb <../../examples/notebooks/15_agent_based.ipynb>`_ +.. math:: + + s^{k+1}_i \;=\; (1 - \alpha)\, s^k_i \;+\; \alpha\, \bar s^k \;+\; \xi^k_i, + \qquad \bar s^k \;=\; \frac1N \sum_{j=1}^N s^k_j, + \qquad \xi^k_i \sim \mathcal{N}(0, \sigma^2), + +with $\alpha \in (0, 1]$ the *averaging weight* and $\sigma$ the noise scale. This is the +DeGroot–Friedkin–Johnsen baseline of opinion dynamics, and the *complete-graph* limit of the +Hegselmann–Krause and Vicsek flocking models. + +Mathematical background +----------------------- + +**Mean conservation.** Averaging the update over $i$ gives +$\bar s^{k+1} = \bar s^k + \bar\xi^k$ with $\mathbb{E}[\bar\xi^k] = 0$, so the empirical mean +is a *martingale* and is exactly preserved in expectation: + +.. math:: + + \mathbb{E}[\bar s^k] \;=\; \bar s^0 \quad \text{for all } k \ge 0. + +In the noiseless case $\sigma = 0$ the mean is preserved *path-by-path*. + +**Geometric contraction of the spread.** Define the deviation $d^k_i := s^k_i - \bar s^k$. +The update implies + +.. math:: + + d^{k+1}_i \;=\; (1 - \alpha)\, d^k_i \;+\; \bigl(\xi^k_i - \bar\xi^k\bigr) , + +so in the absence of noise $\| d^k \|_\infty \le (1 - \alpha)^k \| d^0 \|_\infty$ — the spread +*contracts geometrically* with rate $1 - \alpha$. The companion notebook plots +$\max_i s^k_i - \min_i s^k_i$ on a log scale across $\alpha \in \{0.05, \dots, 1\}$ and +recovers exactly this slope. + +**Stationary variance with noise.** Treating the deviation as an AR(1) process with input +variance $\sigma^2 (1 - 1/N)$, the steady-state variance of any single agent's deviation is + +.. math:: + + \mathrm{Var}_\infty(d_i) \;=\; \frac{\sigma^2 (1 - 1/N)}{1 - (1 - \alpha)^2} + \;\xrightarrow[\alpha \to 0]{}\; \frac{\sigma^2}{2\alpha}\,(1 - 1/N). + +**Continuous-time limit (linear Vlasov).** Sending $\alpha = \theta\, \Delta t$, +$\xi^k_i = \sigma \sqrt{\Delta t}\, W^i_k$ and $\Delta t \to 0$ recovers the McKean–Vlasov SDE +$dX^i_t = \theta(\bar X_t - X^i_t)\, dt + \sigma\, dW^i_t$ of :doc:`mckean_vlasov` — the +discrete consensus update is the prototype of mean-field interaction. + +**Spectral interpretation.** On a general weighted graph the update reads +$s^{k+1} = (I - \alpha L)\, s^k + \xi^k$, where $L$ is the normalised Laplacian. The +complete-graph case shipped here has $L = I - \tfrac1N \mathbf{1}\mathbf{1}^\top$ with +eigenvalue $1$ on the orthogonal complement of $\mathbf{1}$, hence the contraction rate +$1 - \alpha$ above. Replacing $\mathbf{1}\mathbf{1}^\top / N$ by an arbitrary stochastic +matrix produces the full DeGroot model and is a one-liner extension on the Rust side. + +Why it matters +-------------- + +* **Opinion dynamics & social learning.** Calibration of polarisation/consensus models + (Bayesian persuasion, social media echo chambers, voting-system stability). +* **Distributed estimation & federated learning.** Average-consensus protocols for sensor + networks, gossip algorithms, federated averaging — all reduce to the same contraction + argument with explicit convergence rate $1 - \alpha$. +* **Coupled-oscillator physics.** Linear approximation of the Kuramoto / Vicsek models near + the synchronised regime; direct comparison with the McKean–Vlasov continuous limit. + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 15 — Agent-based dynamics ========================= diff --git a/docs/source/algorithms/bsde.rst b/docs/source/algorithms/bsde.rst index 192007a..03cc535 100644 --- a/docs/source/algorithms/bsde.rst +++ b/docs/source/algorithms/bsde.rst @@ -1,10 +1,96 @@ BSDE — θ-scheme and deep-BSDE bridge ==================================== -This notebook exercises `optimizr.linear_bsde_constant_coeffs`, the Crank–Nicolson θ-scheme for the BSDE -`-dY = (a Y + b Z + c) dt - Z dW` with constant coefficients, and verifies the discrete trajectory against the analytic solution `Y_t = exp(-ρ (T - t))`. +A **backward stochastic differential equation** (BSDE) on $[0, T]$ is the inverse-time problem -.. note:: Companion executed notebook: `10_bsde.ipynb <../../examples/notebooks/10_bsde.ipynb>`_ +.. math:: + + Y_t \;=\; \xi \;+\; \int_t^T f(s, Y_s, Z_s)\, ds \;-\; \int_t^T Z_s\, dW_s, + \qquad Y_T = \xi, + +where $\xi \in L^2(\mathcal{F}_T)$ is the *terminal condition*, $f$ is the *driver* and the +unknowns are an adapted pair $(Y, Z) \in \mathcal{S}^2 \times \mathcal{H}^2$. The auxiliary +process $Z$ is a *non-anticipative hedge*: it makes the equation adapted despite the terminal +constraint. + +The primitive `linear_bsde_constant_coeffs` solves the constant-coefficient linear case + +.. math:: + + -dY_t \;=\; (a\, Y_t + b\, Z_t + c)\, dt \;-\; Z_t\, dW_t, + \qquad Y_T = \xi, + +by a **Crank–Nicolson θ-scheme** (θ = 0.5 → second-order in $\Delta t$). + +Mathematical background +----------------------- + +**Pardoux–Peng theorem (1990).** If $f$ is uniformly Lipschitz in $(y, z)$ and +$\mathbb{E}\!\int_0^T f(s, 0, 0)^2\, ds < \infty$, then the BSDE admits a unique solution +$(Y, Z) \in \mathcal{S}^2 \times \mathcal{H}^2$. The proof is a Banach–Picard fixed point on +$\Phi : (Y, Z) \mapsto (Y', Z')$ with +$Y'_t = \mathbb{E}\bigl[\xi + \int_t^T f(s, Y_s, Z_s)\, ds \bigm| \mathcal{F}_t\bigr]$ and $Z'$ +obtained by the martingale representation theorem. + +**Closed-form for the linear case.** For $a, b, c$ deterministic the solution is the +conditional expectation under a Girsanov-shifted measure: + +.. math:: + + Y_t \;=\; \mathbb{E}\!\left[\, \xi\, e^{\int_t^T a(s)\, ds} + \;+\; \int_t^T c(s)\, e^{\int_t^s a(r)\, dr}\, ds + \,\Big|\, \mathcal{F}_t \right], + +with the Girsanov density $\frac{d\mathbb{Q}}{d\mathbb{P}} = \mathcal{E}\bigl(\int_0^\cdot b(s)\,dW_s\bigr)$. +When $b = c = 0$, $a \equiv -\rho$ and $\xi = 1$ this collapses to the analytic ground truth +$Y_t = e^{-\rho(T-t)}$ used by the convergence test. + +**Feynman–Kac bridge.** Setting $f(s, y, z) = -r y$ and $\xi = g(X_T)$ for a forward SDE $X$ +recovers the discounted-payoff PDE: $Y_t = e^{-r(T-t)} \mathbb{E}[g(X_T) \mid \mathcal{F}_t]$. +More generally, the markovian BSDE + +.. math:: + + Y_t = g(X_T) + \int_t^T f(s, X_s, Y_s, Z_s)\, ds - \int_t^T Z_s\, dW_s, + +is the probabilistic representation of the semilinear PDE +$\partial_t u + \mathcal{L}u + f(t, x, u, \sigma^\top \nabla u) = 0$, $u(T, x) = g(x)$, with +$Y_t = u(t, X_t)$ and $Z_t = \sigma^\top(t, X_t)\nabla u(t, X_t)$. + +**Crank–Nicolson θ-scheme.** On a uniform grid $0 = t_0 < \cdots < t_N = T$ the scheme reads + +.. math:: + + Y^N_{t_i} \;=\; \mathbb{E}\!\bigl[\, Y^N_{t_{i+1}} \,\big|\, \mathcal{F}_{t_i}\bigr] + \;+\; \Delta t\,\bigl(\theta\, f(t_i, Y^N_{t_i}, Z^N_{t_i}) + + (1-\theta)\, f(t_{i+1}, Y^N_{t_{i+1}}, Z^N_{t_{i+1}})\bigr), + +with $Z^N_{t_i} = \Delta t^{-1}\,\mathbb{E}\bigl[Y^N_{t_{i+1}}(W_{t_{i+1}} - W_{t_i})\bigm|\mathcal{F}_{t_i}\bigr]$ +(discrete Clark–Ocone identity). For $\theta = 1/2$ the global truncation error is +$\sup_i \mathbb{E}|Y_{t_i} - Y^N_{t_i}|^2 = O(\Delta t^2)$ — the second-order rate verified +empirically by the convergence cell of the companion notebook. + +**Deep-BSDE bridge (E–Han–Jentzen, 2017).** In high dimension the conditional expectation +is intractable; one parametrises $Z_{t_i} = \zeta^i_\theta(X_{t_i})$ by a neural network and +minimises $\mathbb{E}\bigl[(Y^\theta_T - \xi)^2\bigr]$ over $(Y_0, \theta)$. The trait +`ConditionalExpectation` and the struct `DeepBsdeBridge` expose the same θ-scheme step so the +user can plug in any regression / neural-network conditional-expectation oracle. + +Why it matters +-------------- + +* **Pricing & hedging in incomplete markets.** $Y_t$ is the super-replication price of the + contingent claim $\xi$ and $Z_t$ is the instantaneous hedge ratio. Constraints (transaction + costs, portfolio caps, recursive utilities) are absorbed into the driver $f$. +* **Stochastic control.** Forward–backward SDEs are the probabilistic counterpart of the + Hamilton–Jacobi–Bellman PDE; deep-BSDE solves HJB up to $d \sim 100$ state variables, well + beyond grid-based PDE solvers. +* **Risk-sensitive optimisation.** Quadratic-driver BSDE + $-dY = \tfrac1{2\eta}|Z|^2 dt - Z\, dW$ encodes exponential utility hedging (Kramkov–Schachermayer 1999). + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 10 — BSDE θ-scheme ================== diff --git a/docs/source/algorithms/generative_calibration_hooks.rst b/docs/source/algorithms/generative_calibration_hooks.rst index bf47c2b..da311cf 100644 --- a/docs/source/algorithms/generative_calibration_hooks.rst +++ b/docs/source/algorithms/generative_calibration_hooks.rst @@ -1,9 +1,85 @@ -Generative calibration — Gaussian MMD loss +Generative calibration — Gaussian-MMD loss ========================================== -Maximum-Mean-Discrepancy distance with Gaussian kernel (`mmd_gaussian`). Self-distance is exactly zero; the metric grows monotonically with sample shift. +Kernel-based **Maximum Mean Discrepancy** distance (Gretton et al. 2012) — a closed-form, +differentiable, distribution-free metric between two empirical samples. Used as the loss +function of every generative-calibration loop in `optimiz-rs`. -.. note:: Companion executed notebook: `17_generative_calibration.ipynb <../../examples/notebooks/17_generative_calibration.ipynb>`_ +Mathematical background +----------------------- + +**Definition.** For a positive-definite kernel $k : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ +with reproducing-kernel Hilbert space (RKHS) $\mathcal{H}_k$, the *kernel mean embedding* of a +probability measure $P$ is $\mu_P := \mathbb{E}_{X \sim P}[k(X, \cdot)] \in \mathcal{H}_k$. +The **squared MMD** is the RKHS distance between embeddings: + +.. math:: + + \mathrm{MMD}^2(P, Q) \;:=\; \| \mu_P - \mu_Q \|_{\mathcal{H}_k}^2 + \;=\; \mathbb{E}\,[k(X, X')] \;-\; 2\, \mathbb{E}\,[k(X, Y)] \;+\; \mathbb{E}\,[k(Y, Y')] , + +where $X, X' \sim P$ and $Y, Y' \sim Q$ are independent. When $k$ is *characteristic* +(e.g. Gaussian RBF), $\mathrm{MMD}(P, Q) = 0 \iff P = Q$. + +**U-statistic estimator.** Given i.i.d. samples $\{x_i\}_{i=1}^n$ and $\{y_j\}_{j=1}^m$, the +unbiased estimator is + +.. math:: + + \widehat{\mathrm{MMD}}^2 \;=\; + \frac{1}{n(n-1)}\!\sum_{i \ne i'} k(x_i, x_{i'}) + \;-\; \frac{2}{n m}\!\sum_{i, j} k(x_i, y_j) + \;+\; \frac{1}{m(m-1)}\!\sum_{j \ne j'} k(y_j, y_{j'}) . + +It is unbiased, computable in $O((n + m)^2)$ for $d = 1$ (the case implemented), and asymptotically +normal under the alternative. Self-distance is **exactly zero**. + +**Kernel.** The shipped routine uses the Gaussian RBF +$k_\sigma(x, y) = \exp\!\bigl(-(x - y)^2 / (2\sigma^2)\bigr)$ with bandwidth $\sigma$. Standard +reproducing-kernel theory shows that this kernel is *characteristic*, hence MMD metrises weak +convergence on bounded subsets. + +**Closed forms for two notable cases.** + +* **Pure translation, equal samples.** If $Q$ is the law of $X + \Delta$ with $X \sim P$ on + $\mathbb{R}$ and $P = \delta$ atomic, the squared MMD is $2 - 2 e^{-\Delta^2 / (2\sigma^2)}$ — + smooth, monotone in $|\Delta|$, asymptote $2$ as $\Delta \to \infty$. This is the analytic + ground-truth verified by the *bandwidth dependence* cell of the companion notebook. +* **Two Gaussians.** For $P = \mathcal{N}(\mu_1, \sigma_1^2)$ and $Q = \mathcal{N}(\mu_2, \sigma_2^2)$, + + .. math:: + + \mathrm{MMD}^2_\sigma(P, Q) \;=\; + \frac{\sigma}{\sqrt{\sigma^2 + 2\sigma_1^2}} + \;-\; \frac{2\sigma}{\sqrt{\sigma^2 + \sigma_1^2 + \sigma_2^2}}\, e^{-\frac{(\mu_1 - \mu_2)^2}{2(\sigma^2 + \sigma_1^2 + \sigma_2^2)}} + \;+\; \frac{\sigma}{\sqrt{\sigma^2 + 2\sigma_2^2}} , + + giving an exact reference for unit tests. + +**Statistical guarantee.** Gretton et al. (2012, Thm. 12) give the deviation bound +$\Pr\!\bigl(\widehat{\mathrm{MMD}}^2 - \mathrm{MMD}^2 > \varepsilon\bigr) \le \exp\bigl(-\varepsilon^2 nm / (8 K^2 (n + m))\bigr)$ +for $|k| \le K$. Hence MMD detects fixed alternatives at the optimal $n^{-1/2}$ rate. + +**Connection with Wasserstein.** Both metrise weak convergence, but MMD is *quadratic in the +sample size* (no transport plan to solve) and admits unbiased low-variance gradient estimators — +the reason it is the loss of choice in implicit-generative-model training +(generator-loss / score-matching alternatives). + +Why it matters +-------------- + +* **Generative calibration.** Train an implicit sampler (neural SDE, copula generator, GAN-like + architecture) by minimising $\widehat{\mathrm{MMD}}^2$ between the simulator output and the + target distribution. The trait `GenerativeSampler` plus `calibration_step` is the abstract + glue. +* **Two-sample testing.** Distribution drift detection in streaming data, A/B-test signal + extraction, anomaly detection. +* **Model selection.** Replace likelihood ratios when likelihoods are intractable + (simulator-based inference, ABC). + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 17 — MMD calibration loss ========================= diff --git a/docs/source/algorithms/mckean_vlasov.rst b/docs/source/algorithms/mckean_vlasov.rst index a6bde62..2a2ffbf 100644 --- a/docs/source/algorithms/mckean_vlasov.rst +++ b/docs/source/algorithms/mckean_vlasov.rst @@ -1,9 +1,85 @@ McKean–Vlasov — propagation of chaos ==================================== -Interacting-particle Euler scheme for $dX_t = θ(\bar X_t - X_t) dt + σ dW_t$ (`mean_reverting_mckean_vlasov`). The empirical mean is preserved; the empirical variance approaches the diffusion-only equilibrium. +A **McKean–Vlasov SDE** is a stochastic differential equation whose drift and diffusion depend +on the *law* of the solution itself: -.. note:: Companion executed notebook: `14_mckean_vlasov.ipynb <../../examples/notebooks/14_mckean_vlasov.ipynb>`_ +.. math:: + + dX_t \;=\; b\bigl(t, X_t, \mathcal{L}(X_t)\bigr)\, dt \;+\; \sigma\bigl(t, X_t, \mathcal{L}(X_t)\bigr)\, dW_t, + \qquad X_0 \sim \mu_0 . + +It is the formal $N \to \infty$ limit of an exchangeable system of $N$ interacting diffusions + +.. math:: + + dX^{i,N}_t \;=\; b\!\Bigl(t, X^{i,N}_t, \tfrac1N\!\sum_{j=1}^N \delta_{X^{j,N}_t}\Bigr)\, dt + \;+\; \sigma\!\Bigl(t, X^{i,N}_t, \tfrac1N\!\sum_{j=1}^N \delta_{X^{j,N}_t}\Bigr)\, dW^i_t . + +The primitive shipped here, `mean_reverting_mckean_vlasov`, simulates the canonical example + +.. math:: + + dX_t \;=\; \theta\bigl(\bar X_t - X_t\bigr)\, dt \;+\; \sigma\, dW_t, + \qquad \bar X_t = \mathbb{E}[X_t], + +with the symmetric Euler particle scheme $X^{i,N}_{k+1} = X^{i,N}_k + \theta(\bar X^N_k - X^{i,N}_k)\Delta t + \sigma\sqrt{\Delta t}\,\xi^i_k$. + +Mathematical background +----------------------- + +**Sznitman's propagation of chaos (1991).** Under standard Lipschitz assumptions on $b, \sigma$ in +$(x, \mu)$ (the $\mu$ argument equipped with the Wasserstein distance $W_2$), the empirical +measure $\mu^N_t = \tfrac1N \sum_i \delta_{X^{i,N}_t}$ converges weakly to the deterministic flow +$\mathcal{L}(X_t)$, and any fixed sub-system of $k$ particles becomes asymptotically independent: + +.. math:: + + \sup_{0 \le t \le T} \, \mathbb{E}\bigl[\,W_2^2\!\bigl(\mu^N_t,\, \mathcal{L}(X_t)\bigr)\bigr] + \;\le\; \frac{C(T)}{N^{2/(d+4)}} . + +**Density flow (nonlinear Fokker–Planck).** The marginal density $\rho_t = \mathrm{law}(X_t)$ +satisfies the *nonlinear* PDE + +.. math:: + + \partial_t \rho_t \;+\; \nabla\!\cdot\!\bigl(b(t, x, \rho_t)\, \rho_t\bigr) + \;=\; \tfrac12\, \nabla^2\!:\!\bigl(\sigma\sigma^\top(t, x, \rho_t)\, \rho_t\bigr). + +**Closed-form for the mean-reverting case.** Taking expectation of the SDE gives +$\dot{\bar X}_t = 0$, so the population mean is *exactly preserved*: $\bar X_t \equiv \bar X_0$. +The deviation $\widetilde X^i_t := X^{i,N}_t - \bar X_0$ then solves a standard Ornstein–Uhlenbeck +SDE, so each marginal is Gaussian with + +.. math:: + + \mathbb{E}[X_t] \;=\; \bar X_0, + \qquad + \mathrm{Var}(X_t) \;=\; \mathrm{Var}(X_0)\, e^{-2\theta t} \;+\; \frac{\sigma^2}{2\theta}\bigl(1 - e^{-2\theta t}\bigr) + \;\xrightarrow[t\to\infty]{}\; \frac{\sigma^2}{2\theta}. + +The companion notebook checks both the mean conservation and the variance asymptote. + +**Connection with mean-field BSDEs.** Coupling the McKean–Vlasov forward SDE with a backward +equation $-dY_t = f(t, X_t, Y_t, Z_t, \mathcal{L}(X_t, Y_t))\, dt - Z_t\, dW_t$ produces the +*mean-field BSDE* of Carmona–Delarue (2018), itself the probabilistic representation of the +HJB side of mean-field games (cf. :doc:`stochastic_control`). + +Why it matters +-------------- + +* **Mean-field games.** At the Nash equilibrium of a symmetric $N$-player game, each player's + state follows a McKean–Vlasov SDE in which the population law $\mu_t$ is the consistent + fixed point of every player's best response. This is the master tool of Lasry–Lions theory + for systemic-risk modelling, optimal execution and price formation. +* **Statistical physics.** Vlasov, Boltzmann, and granular-media equations all arise as + density flows of mean-field particle systems; the same Euler scheme estimates their solutions. +* **Generative modelling.** Stein-variational gradient descent and score-based diffusion can + be analysed as McKean–Vlasov gradient flows on $W_2$. + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 14 — McKean–Vlasov mean-reverting dynamics ========================================== diff --git a/docs/source/algorithms/pde.rst b/docs/source/algorithms/pde.rst index 4edf4b3..00cca3b 100644 --- a/docs/source/algorithms/pde.rst +++ b/docs/source/algorithms/pde.rst @@ -1,9 +1,83 @@ PDE — Fokker–Planck, HJB, elliptic Poisson ========================================== -Three CPU-only finite-difference solvers: 1-D forward Fokker–Planck (`fokker_planck_constant`), 2-D explicit HJB (`hjb_quadratic_2d`) and 2-D Poisson SOR (`poisson_2d_zero_boundary`). Each routine is verified against an analytic ground truth. +Three CPU-only finite-difference solvers covering the two canonical PDE pillars of stochastic +analysis: the **forward** equation for the marginal density of a diffusion (Fokker–Planck), +the **backward** equation for an optimally controlled diffusion (Hamilton–Jacobi–Bellman), +and a static **elliptic** boundary-value problem (Poisson). -.. note:: Companion executed notebook: `11_pde.ipynb <../../examples/notebooks/11_pde.ipynb>`_ +Mathematical background +----------------------- + +**Fokker–Planck (Kolmogorov forward).** For a 1-D Itô diffusion +$dX_t = \mu(t, x)\, dt + \sigma(t, x)\, dW_t$, the marginal density $\rho(t, x)$ of $X_t$ +satisfies the parabolic PDE + +.. math:: + + \partial_t \rho \;+\; \partial_x\!\bigl(\mu(t,x)\, \rho\bigr) + \;=\; \tfrac12\, \partial^2_{xx}\!\bigl(\sigma^2(t,x)\, \rho\bigr), + \qquad \rho(0, \cdot) = \rho_0 . + +For the *pure-diffusion* test ($\mu \equiv 0$, $\sigma^2 \equiv 1$, $\rho_0 = \mathcal{N}(0, 1)$) +the analytic Gaussian heat kernel gives $\rho(t, x) = \frac{1}{\sqrt{2\pi(1+t)}}\exp\!\bigl(-\frac{x^2}{2(1+t)}\bigr)$, +so the variance grows linearly: $\mathrm{Var}(X_t) = 1 + t$. The conservative +Lax–Wendroff / centred-flux scheme implemented by `fokker_planck_constant` preserves total mass +(checked in the notebook to machine precision). + +**Hamilton–Jacobi–Bellman.** Consider the controlled diffusion +$dX_t = \mu(X_t, \alpha_t)\, dt + \sigma(X_t)\, dW_t$ and the value function +$v(t, x) = \sup_\alpha \mathbb{E}_{t,x}\!\bigl[\int_t^T r(X_s, \alpha_s)\, ds + g(X_T)\bigr]$. +Dynamic programming produces + +.. math:: + + \partial_t v \;+\; \sup_{a \in \mathcal{A}}\Bigl\{ \mu(x, a) \cdot \nabla v + \;+\; \tfrac12\, \mathrm{tr}\!\bigl(\sigma\sigma^\top(x)\, \nabla^2 v\bigr) + \;+\; r(x, a) \Bigr\} \;=\; 0, + \qquad v(T, x) = g(x). + +`hjb_quadratic_2d` discretises this in 2-D by an explicit finite-difference scheme; the simple +heat-only relaxation case (:math:`H \equiv 0`, :math:`\sigma^2 > 0`) preserves a constant value while a +quadratic terminal :math:`g(x) = \tfrac12 \lVert x \rVert^2` smooths into a Gaussian-shaped value surface. + +**Elliptic Poisson with zero Dirichlet boundary.** On the unit square $\Omega = (0,1)^2$, + +.. math:: + + -\Delta u(x, y) = f(x, y) \text{ in } \Omega, \qquad u\!\restriction_{\partial\Omega} = 0 . + +The Laplace eigenfunctions $\phi_{m,n}(x, y) = \sin(m\pi x)\sin(n\pi y)$ form an +orthonormal basis with eigenvalues $\lambda_{m,n} = (m^2 + n^2)\pi^2$, so for +$f = 2\pi^2 \sin(\pi x)\sin(\pi y)$ the *exact* solution is +$u(x, y) = \sin(\pi x)\sin(\pi y)$. `poisson_2d_zero_boundary` solves the 5-point stencil by +**Successive Over-Relaxation** with optimal relaxation parameter +$\omega^* = 2 / (1 + \sin(\pi h))$ for grid spacing $h = 1/(N-1)$, achieving spectral radius +$\rho \sim 1 - 2\pi h$ — i.e. $O(h^{-1})$ iterations to reach a fixed tolerance, against +$O(h^{-2})$ for plain Gauss–Seidel. + +**Probabilistic representation (Feynman–Kac).** Both the parabolic HJB and the elliptic +Poisson PDE admit stochastic representations: $u(x) = \mathbb{E}_x\!\bigl[\int_0^{\tau_\Omega} f(X_s)\, ds\bigr]$ +for the latter, where $\tau_\Omega$ is the first exit time of the diffusion from $\Omega$. +This links the PDE solvers above to the BSDE primitives of :doc:`bsde`. + +Why it matters +-------------- + +* **Density estimation under controlled noise.** Fokker–Planck is the workhorse of + non-equilibrium statistical physics, plasma transport, calibration of stochastic-volatility + models, and Langevin-based MCMC convergence diagnostics. +* **Optimal control & inverse problems.** HJB is the cornerstone of dynamic programming, + reinforcement learning (continuous-time policy iteration), and stochastic-control routing. +* **Mean-field games.** The MFG fixed point is exactly the coupled system + *(backward HJB + forward Fokker–Planck)* with cost depending on the density — building this + loop on top of the two solvers above is one of the v2.0 milestones. +* **Image processing & PDE-constrained optimisation.** Poisson editing, electric-potential + reconstruction, gravitational-potential inversion all reduce to the same elliptic stencil. + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 11 — PDE solvers ================ diff --git a/docs/source/algorithms/quadratic_impact_control.rst b/docs/source/algorithms/quadratic_impact_control.rst index e1fbe71..b7ffff3 100644 --- a/docs/source/algorithms/quadratic_impact_control.rst +++ b/docs/source/algorithms/quadratic_impact_control.rst @@ -1,9 +1,83 @@ Quadratic-impact control — closed-form Riccati ============================================== -Closed-form Riccati feedback for a controlled 1-D SDE with quadratic running cost (`quadratic_impact_control_py`). +Closed-form Riccati feedback for the canonical *single-state, quadratic-cost* linear control +problem with running quadratic *impact* penalty. -.. note:: Companion executed notebook: `13_quadratic_impact.ipynb <../../examples/notebooks/13_quadratic_impact.ipynb>`_ +Mathematical background +----------------------- + +Let $A_t$ be a controlled scalar state driven by an additive control $u_t$ and Gaussian noise. +The controller minimises the *finite-horizon quadratic objective* + +.. math:: + + J(u) \;=\; \mathbb{E}\!\left[\,\int_0^T \bigl(\,\tfrac{\gamma}{2}\, u_t^2 + \;+\; \tfrac{\phi}{2}\, A_t^2 \,\bigr)\, dt + \;+\; \tfrac{A_T}{2}\, A_T^2 \,\right] , + +where $\gamma > 0$ is the **impact / control cost**, $\phi \ge 0$ the **running risk weight** +and $A_T$ the **terminal penalty** (over-loaded notation: $A_T$ here is the *coefficient*). + +**Hamilton–Jacobi–Bellman.** With value function $v(t, A) = \tfrac12 h(t)\, A^2 + c(t)$, the +HJB equation collapses to a scalar Riccati ODE on $h$: + +.. math:: + + h'(t) \;=\; \frac{h(t)^2}{\gamma} \;-\; \phi, + \qquad + h(T) \;=\; A_T . + +The optimal feedback is the linear law + +.. math:: + + u^*(t, A) \;=\; -\, \frac{h(t)}{\gamma}\, A \;\equiv\; -\, k(t)\, A, + +with *feedback gain* $k(t) = h(t) / \gamma$. This is the structure returned by the primitive. + +**Closed-form solutions.** + +* **Symmetric fixed point** $\gamma = \phi = A_T = 1$: $h(t) \equiv 1$ is the unique solution + (RHS vanishes), so the feedback gain is constant $k \equiv 1$. The notebook checks this + to machine precision. +* **Generic $\phi > 0$.** Writing $\bar h = \sqrt{\gamma \phi}$ for the steady-state and + $\rho = \sqrt{\phi / \gamma}$, the Riccati ODE has the closed-form (separation of variables / + Bernoulli substitution) + + .. math:: + + h(t) \;=\; \bar h\, \frac{(\bar h + A_T)\, e^{2\rho(T-t)} \;-\; (\bar h - A_T)} + {(\bar h + A_T)\, e^{2\rho(T-t)} \;+\; (\bar h - A_T)} . + + In the limit $T - t \to \infty$ the trajectory relaxes to the stationary value $\bar h = \sqrt{\gamma\phi}$. +* **Free of running risk** $\phi = 0$. Then $h'(t) = h(t)^2/\gamma$ integrates explicitly to + + .. math:: + + h(t) \;=\; \frac{A_T}{1 + (A_T / \gamma)(T - t)} , + + recovering the Pontryagin LQR closed form $P(0) = 1/2$ of :doc:`stochastic_control`. + +**Connection with mean-field games.** Coupling this single-agent control with an interacting +population — the running cost depending on the *average* control $\bar u_t$ — yields the +Almgren–Chriss MFG (Lasry–Lions 2007); at the Nash equilibrium the optimal trajectory is the +uniform schedule $\dot A^*_t = -A_0 / T$ (cf. Sec. 3 of Carmona–Delarue 2018, Vol. I). + +Why it matters +-------------- + +* **Optimal execution.** Almgren–Chriss and its mean-field variants reduce to exactly this + Riccati ODE; the closed form means *real-time* feedback re-computation. +* **Stochastic regulators.** Temperature stabilisation, attitude control, queueing-network + smoothing all map to a quadratic-impact problem with a single state. +* **Building block for higher-dimensional MPC.** Vector generalisations of $h(t)$ are matrix + Riccati ODEs; this scalar primitive is the verification kernel against which the matrix + solver in :doc:`matrix_riccati` is tested. + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 13 — Quadratic-impact controlled SDE ==================================== diff --git a/docs/source/algorithms/robust_drift.rst b/docs/source/algorithms/robust_drift.rst index 852efc5..be7fd7c 100644 --- a/docs/source/algorithms/robust_drift.rst +++ b/docs/source/algorithms/robust_drift.rst @@ -1,9 +1,92 @@ -Inference — Huber-IRLS drift estimator -====================================== +Inference — Huber-IRLS robust drift estimator +============================================= -Robust drift estimator (`robust_drift`) for $x_{k+1} = x_k + (a + b x_k) Δt + σ ε_k$ via Huber IRLS — resists 5 % heavy-tailed innovations. +Heavy-tail-resistant maximum-likelihood estimator for the discrete Ornstein–Uhlenbeck-type model -.. note:: Companion executed notebook: `16_robust_drift.ipynb <../../examples/notebooks/16_robust_drift.ipynb>`_ +.. math:: + + x_{k+1} \;=\; x_k \;+\; (a + b\, x_k)\, \Delta t \;+\; \sigma\, \sqrt{\Delta t}\, \varepsilon_k, + \qquad \varepsilon_k \sim_{\text{i.i.d.}} P_\varepsilon , + +where $P_\varepsilon$ is *contaminated*: a fraction $1 - \eta$ of standard Gaussian innovations +plus a fraction $\eta$ of large outliers (jumps, fat tails, recording errors). + +Mathematical background +----------------------- + +**Naive OLS.** Setting $y_k := (x_{k+1} - x_k)/\Delta t$, the model is the linear regression +$y_k = a + b\, x_k + \sigma\, \Delta t^{-1/2}\, \varepsilon_k$. Ordinary least-squares +minimises $\sum_k (y_k - a - b x_k)^2$ but its breakdown point is $0$: a single outlier with +$|\varepsilon_k| \gg 1$ moves the estimate arbitrarily far. + +**Huber loss & IRLS.** Huber (1964) replaces the quadratic loss by the *piecewise* loss + +.. math:: + + \rho_\delta(r) \;=\; + \begin{cases} + \tfrac12\, r^2, & |r| \le \delta, \\[2pt] + \delta\,\bigl(|r| - \tfrac\delta2\bigr), & |r| > \delta, + \end{cases} + +which is *quadratic in the bulk* and *linear in the tails*. The first-order condition +$\sum_k \psi_\delta(r_k)\, \nabla_{a,b}\, r_k = 0$ with $\psi_\delta = \rho_\delta'$ rewrites +as a weighted least-squares problem with weights + +.. math:: + + w_k \;=\; \min\!\Bigl(1,\; \frac{\delta}{|r_k|}\Bigr) , + +so the **Iteratively Reweighted Least-Squares** algorithm reads + +.. math:: + + \widehat{(a, b)}^{(t+1)} \;=\; \arg\min_{a, b}\; \sum_k w^{(t)}_k\, (y_k - a - b\, x_k)^2, + \qquad w^{(t+1)}_k = \min\!\bigl(1, \delta / |r^{(t+1)}_k|\bigr). + +The sequence converges geometrically when the design matrix is well-conditioned +(Holland–Welsch 1977). `robust_drift` returns the limit pair $(\widehat a, \widehat b)$ and +the number of iterations. + +**Choice of the cut-off.** The default $\delta = 1.345 \cdot \hat\sigma$ delivers $95\%$ +asymptotic efficiency under Gaussian innovations while keeping the influence function bounded; +it is the Huber–Hampel value used as the standard reference in robust statistics. + +**Closed-form one-step (debiased OLS).** When the contamination is symmetric and the +innovations have finite variance $\sigma^2_\varepsilon$, the *consistent* one-step estimate at +the ordinary least-squares solution $(\hat a^0, \hat b^0)$ reads + +.. math:: + + \binom{\widehat a}{\widehat b} + \;=\; + \binom{\hat a^0}{\hat b^0} + \;+\; \bigl(X^\top W X\bigr)^{-1}\, X^\top \psi_\delta(r^0), + +where $X$ is the $(N - 1) \times 2$ design matrix and $W = \mathrm{diag}(w_k)$. Bahadur +linearisation shows $\widehat\theta - \theta^\star = O_P(N^{-1/2})$ even in the contaminated +model, with asymptotic variance $\sigma^2_\psi / I^2_\psi$ (Huber, *Robust Statistics*, 2004, +Thm. 7.7). + +**Connection with Malliavin calculus.** The driver $a + b\, x$ is exactly the linearised +drift of the Ornstein–Uhlenbeck process used in the Greeks formulae of +:doc:`stochastic_control` and the Vasicek interest-rate model; robust calibration is the +pre-requisite for any Monte-Carlo Greeks computation under noisy historical data. + +Why it matters +-------------- + +* **Heavy-tailed historical data.** Crypto returns, electricity prices, plasma confinement + signals, and bio-medical recordings all contain spikes that destroy OLS but leave Huber + estimates within statistical noise. +* **Online & streaming estimation.** IRLS with $\sim 10$ iterations is real-time on streaming + windows and exposes a stable derivative for downstream control loops. +* **Robust risk management.** Replacing raw OLS by IRLS in any volatility / mean-reversion + estimator dramatically reduces *parameter risk* in stress periods. + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 16 — Robust drift estimation ============================ diff --git a/docs/source/algorithms/stochastic_control.rst b/docs/source/algorithms/stochastic_control.rst index dc6d957..809bdaf 100644 --- a/docs/source/algorithms/stochastic_control.rst +++ b/docs/source/algorithms/stochastic_control.rst @@ -1,9 +1,84 @@ Stochastic control — switching, Pontryagin, two-sided intensities ================================================================= -Three primitives: discrete-time optimal switching (`optimal_switching_dp`), 1-D Pontryagin LQR shooting (`pontryagin_lqr`) and the bilateral intensity controller (`two_sided_intensities`). +Three complementary primitives covering the discrete and continuous worlds of stochastic +control: dynamic-programming **optimal switching** (Snell envelope), the continuous-time +**Pontryagin–Bismut maximum principle** for the linear-quadratic regulator, and a **two-sided +intensity controller** for jump processes. -.. note:: Companion executed notebook: `12_stochastic_control.ipynb <../../examples/notebooks/12_stochastic_control.ipynb>`_ +Mathematical background +----------------------- + +**1. Optimal switching as a Snell envelope.** Let $(Y^i_k)_{k, i}$ be the running rewards in +mode $i \in \{1, \dots, M\}$ and $c_{ij}$ the cost of switching from $i$ to $j$. The value +function $V_k(i)$ satisfies the backward dynamic-programming recursion + +.. math:: + + V_N(i) = g(i), + \qquad + V_k(i) \;=\; Y^i_k \;+\; \max_{j}\!\bigl( V_{k+1}(j) - c_{ij}\bigr). + +This is the *multi-mode Snell envelope* of El Karoui–Quenez (1995). When switching is free +($c_{ij} = 0$) and only mode 1 pays a unit reward at every period, $V_k(i) = N - k$ for +$i \neq 1$ and $V_k(1) = N - k + 1$ — reproduced exactly by `optimal_switching_dp`. + +**2. Pontryagin–Bismut maximum principle (LQR).** For the controlled SDE +$dX_t = (a X_t + b u_t)\, dt + \sigma\, dW_t$ with quadratic cost +$J(u) = \mathbb{E}\!\bigl[\int_0^T (q X_t^2 + r u_t^2)\, dt + s_T X_T^2\bigr]$, the +adjoint variable $P_t$ solves the **matrix Riccati ODE** + +.. math:: + + \dot P_t \;+\; 2 a\, P_t \;-\; \frac{b^2}{r}\, P_t^2 \;+\; q \;=\; 0, + \qquad P_T = s_T, + +and the optimal feedback is $u^*_t = -(b/r)\, P_t\, X_t$. In the canonical case +$a = q = 0$, $b = r = s_T = 1$, $T = 1$ the ODE simplifies to +$\dot P_t = P_t^2$, whose closed-form solution is + +.. math:: + + P_t \;=\; \frac{1}{1 + (T - t)} , + \qquad + P(0) = \tfrac12 . + +The primitive `pontryagin_lqr` reproduces this with relative error below $10^{-3}$ for +$N = 2000$ steps (the symmetric Strang splitting is second-order in $\Delta t$). + +**3. Two-sided intensity control.** For a jump-controller the agent picks the rates +:math:`\lambda_\pm \ge 0` at which up/down events fire. With *affine premia* +:math:`\delta_\pm(\lambda) = \alpha_\pm + \kappa_\pm \lambda` and value-function jumps +:math:`\Delta V_\pm`, the instantaneous Hamiltonian is + +.. math:: + + \sup_{\lambda_\pm \ge 0}\!\Bigl[\,\lambda_+\bigl(\delta_+(\lambda_+) - \Delta V_+\bigr) + \;+\; \lambda_-\bigl(\delta_-(\lambda_-) - \Delta V_-\bigr)\Bigr], + +and the first-order condition gives the closed-form maximiser + +.. math:: + + \lambda^*_\pm \;=\; \max\!\Bigl(0,\; \frac{\alpha_\pm - \Delta V_\pm}{2\, \kappa_\pm}\Bigr). + +The quantity :math:`\Delta V_\pm` is the (estimated) marginal value of an additional event; +`two_sided_intensities` returns :math:`(\lambda^*_+, \lambda^*_-)` in closed form, which is what +lets the broader optimal-execution loop run in real time. + +Why it matters +-------------- + +* **Optimal switching** powers production-mode selection (start/stop a power plant), regime + changes in algorithmic strategies, and American-style option pricing (Carmona–Touzi 2008). +* **Pontryagin LQR** is the linearised core of every continuous-control problem: target + tracking, Kalman-LQG, ground-up RL, robust $H_\infty$ design. +* **Two-sided intensity control** is the closed-form heart of optimal market making + (Avellaneda–Stoikov 2008, Cartea–Jaimungal–Penalva 2015) and limit-order placement. + +.. note:: + 📓 **Companion notebook** — `view on GitHub `_ + · `download .ipynb `_ 12 — Stochastic control ======================= @@ -102,7 +177,7 @@ Closed-form Riccati for $a=q=0$, $b=r=s_T=1$, $T=1$ is $P(t) = 1/(1 + (T - t))$, Two-sided intensity control --------------------------- -Affine premium $δ_±(λ) = α_± + κ_± λ$. First-order condition: $\lambda^*_\pm = \max(0, (α_\pm - ΔV_\pm) / (2 κ_\pm))$. +Affine premium :math:`\delta_\pm(\lambda) = \alpha_\pm + \kappa_\pm \lambda`. First-order condition: :math:`\lambda^*_\pm = \max(0, (\alpha_\pm - \Delta V_\pm) / (2 \kappa_\pm))`. .. code-block:: python