{ "cells": [ { "cell_type": "markdown", "id": "signal-attribution-00", "metadata": { "part": "optimiser" }, "source": [ "# Reproduce the synthetic signal comparison\n", "\n", "Install `numpy scipy cvxpy clarabel`, then run the cells in order. The notebook is self-contained: it needs no checkout or download. NumPy handles the diagonal calculation; CVXPY with Clarabel handles a full covariance matrix. The final cell reproduces both published fixed-seed paths exactly. A separate small dense example checks what changes when assets are correlated.\n", "\n", "The formulas are those in section 3.3 and Bruno Durin's *Signal-wise performance attribution for constrained optimisation* (2014). Its equation numbers below make the code's conventions explicit. All values are synthetic. The fixed-seed examples are not the comparison tables' averages across 20 seeds." ] }, { "cell_type": "code", "execution_count": 1, "id": "signal-attribution-01", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:34.808085Z", "iopub.status.busy": "2026-10-29T10:22:34.808004Z", "iopub.status.idle": "2026-10-29T10:22:35.625287Z", "shell.execute_reply": "2026-10-29T10:22:35.624737Z" }, "part": "optimiser" }, "outputs": [], "source": [ "from dataclasses import dataclass\n", "from math import factorial\n", "import numpy as np\n", "from scipy.optimize import nnls\n", "import cvxpy as cp\n", "import warnings\n", "solver_events = []\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-02", "metadata": { "part": "optimiser" }, "source": [ "## 1. Optimiser\n", "\n", "`Desk` holds the signal coefficients, persistence, cost scales and position bounds. `covariance` returns either the identity matrix or a one-factor covariance with unit diagonal. The static problem uses $\\Lambda=\\lambda\\Sigma$ (section 3.3; Durin Eq. 7), so this compact example supports a general covariance matrix but does not expose an independent cost matrix. `gamma_ref`, `lam` and `lam0` keep the published scaling conventions." ] }, { "cell_type": "code", "execution_count": 2, "id": "signal-attribution-03", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.626648Z", "iopub.status.busy": "2026-10-29T10:22:35.626547Z", "iopub.status.idle": "2026-10-29T10:22:35.630476Z", "shell.execute_reply": "2026-10-29T10:22:35.629988Z" }, "part": "optimiser" }, "outputs": [], "source": [ "@dataclass(frozen=True)\n", "class Desk:\n", " n_assets: int = 50\n", " n_periods: int = 1000\n", " names: tuple = (\"fast A\", \"slow A\", \"fast B\", \"slow B\")\n", " ic: tuple = (0.06, 0.05, 0.04, 0.03)\n", " rho: tuple = (0.5, 0.95, 0.5, 0.95)\n", " corr: float = 0.3\n", " risk_target: float = 1.0\n", " box: float = 1.5\n", " lam_rel: float = 1.0\n", " spread_rel: float = 0.5\n", " cov: str = \"diag\"\n", " long_only: bool = False\n", " factor_share: float = 0.3\n", "\n", " def sigma2(self):\n", " return np.ones(self.n_assets)\n", "\n", " def innovation_corr(self):\n", " C = np.full((len(self.ic), len(self.ic)), float(self.corr))\n", " np.fill_diagonal(C, 1.0)\n", " return C\n", "\n", " def g_std(self):\n", " rho = np.asarray(self.rho, float)\n", " a = np.sqrt(1 - rho**2)\n", " S = self.innovation_corr() * np.outer(a, a) / (1 - np.outer(rho, rho))\n", " b = np.asarray(self.ic, float)\n", " return float(np.sqrt(b @ S @ b))\n", "\n", " def gamma_ref(self):\n", " return np.sqrt(self.n_assets) * self.g_std() / self.risk_target\n", "\n", " def lam(self):\n", " return self.lam_rel * self.gamma_ref()\n", "\n", " def lam0(self):\n", " return self.spread_rel * self.g_std()\n", "\n", " def bounds(self):\n", " hi = np.full(self.n_assets, self.box * self.risk_target / np.sqrt(self.n_assets))\n", " return (np.zeros(self.n_assets) if self.long_only else -hi), hi\n", "\n", " def factor_model(self):\n", " beta = np.linspace(0.8, 1.2, self.n_assets)\n", " f = np.sqrt(self.factor_share) * beta\n", " return f, 1.0 - f**2\n", "\n", " def covariance(self):\n", " if self.cov == \"diag\":\n", " return np.diag(self.sigma2())\n", " f, s2 = self.factor_model()\n", " return np.outer(f, f) + np.diag(s2)\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-04", "metadata": { "part": "optimiser" }, "source": [ "`simulate` generates all AR(1) innovations, then all return noise, in the published order. Returns equal the sum of forecast components plus noise (Durin Eqs. 1–4). The initial scores use the first innovations, leaving a small cross-correlation transient for mixed persistence. The dense example draws its factor noise before its idiosyncratic noise." ] }, { "cell_type": "code", "execution_count": 3, "id": "signal-attribution-05", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.631713Z", "iopub.status.busy": "2026-10-29T10:22:35.631636Z", "iopub.status.idle": "2026-10-29T10:22:35.634047Z", "shell.execute_reply": "2026-10-29T10:22:35.633704Z" }, "part": "optimiser" }, "outputs": [], "source": [ "def simulate(desk, seed):\n", " rng = np.random.default_rng(seed)\n", " K, T, N = len(desk.ic), desk.n_periods, desk.n_assets\n", " rho = np.asarray(desk.rho, float)[:, None]\n", " L = np.linalg.cholesky(desk.innovation_corr() + 1e-12 * np.eye(K))\n", " e = np.einsum(\"kl,ltn->ktn\", L, rng.standard_normal((K, T, N)))\n", " z = np.empty((K, T, N))\n", " z[:, 0] = e[:, 0]\n", " a = np.sqrt(1 - rho**2)\n", " for t in range(1, T):\n", " z[:, t] = rho * z[:, t - 1] + a * e[:, t]\n", " if desk.cov == \"diag\":\n", " u = rng.standard_normal((T, N))\n", " else:\n", " f, s2 = desk.factor_model()\n", " u = rng.standard_normal((T, 1)) * f + rng.standard_normal((T, N)) * np.sqrt(s2)\n", " g = np.asarray(desk.ic, float)[:, None, None] * z\n", " return g, g.sum(0) + u" ] }, { "cell_type": "markdown", "id": "signal-attribution-06", "metadata": { "part": "optimiser" }, "source": [ "`Solution` stores positions, trades, risk aversion, spread subgradients and signed box multipliers. For CVXPY runs, `dual_upper`, `dual_lower` and `primal_raw` retain the solver outputs. In particular, `nu` remains the raw upper-minus-lower dual; it is not overwritten to force attribution sums." ] }, { "cell_type": "code", "execution_count": 4, "id": "signal-attribution-07", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.635064Z", "iopub.status.busy": "2026-10-29T10:22:35.634996Z", "iopub.status.idle": "2026-10-29T10:22:35.637048Z", "shell.execute_reply": "2026-10-29T10:22:35.636726Z" }, "part": "optimiser" }, "outputs": [], "source": [ "@dataclass\n", "class Solution:\n", " \"\"\"Path of the algorithm for a batch of B runs.\"\"\"\n", " x: np.ndarray\n", " dx: np.ndarray\n", " gamma: np.ndarray\n", " s: np.ndarray\n", " nu: np.ndarray\n", " at_bound: np.ndarray\n", " dual_upper: np.ndarray | None = None\n", " dual_lower: np.ndarray | None = None\n", " primal_raw: np.ndarray | None = None\n", " diagnostics: dict | None = None\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-08", "metadata": { "part": "optimiser" }, "source": [ "`_bcast_bounds` repeats the same position bounds for each batched run (Durin Eq. 15)." ] }, { "cell_type": "code", "execution_count": 5, "id": "signal-attribution-09", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.638122Z", "iopub.status.busy": "2026-10-29T10:22:35.638053Z", "iopub.status.idle": "2026-10-29T10:22:35.639833Z", "shell.execute_reply": "2026-10-29T10:22:35.639373Z" }, "part": "optimiser" }, "outputs": [], "source": [ "def _bcast_bounds(lo, hi, B, N):\n", " lo = np.broadcast_to(np.asarray(lo, float), (B, N))\n", " hi = np.broadcast_to(np.asarray(hi, float), (B, N))\n", " return (lo, hi)" ] }, { "cell_type": "markdown", "id": "signal-attribution-10", "metadata": { "part": "optimiser" }, "source": [ "`_trade_diag` is the diagonal closed form: soft-threshold the trade and clip the position (section 3.3; Durin Eqs. 50–51)." ] }, { "cell_type": "code", "execution_count": 6, "id": "signal-attribution-11", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.640966Z", "iopub.status.busy": "2026-10-29T10:22:35.640900Z", "iopub.status.idle": "2026-10-29T10:22:35.642567Z", "shell.execute_reply": "2026-10-29T10:22:35.642240Z" }, "part": "optimiser" }, "outputs": [], "source": [ "def _trade_diag(G, xp, gamma, sig2, lam, lam0, lo, hi):\n", " \"\"\"Closed-form optimal position for given gamma (B,) (Durin Eqs. 50-51, p. 17).\"\"\"\n", " g = gamma[:, None]\n", " q = (lam + g) * sig2\n", " z = G - g * sig2 * xp\n", " dstar = np.sign(z) * np.maximum(np.abs(z) - lam0, 0.0) / q\n", " xu = xp + dstar\n", " return (np.clip(xu, lo, hi), xu, z, q)" ] }, { "cell_type": "markdown", "id": "signal-attribution-12", "metadata": { "part": "optimiser" }, "source": [ "`solve_diag` uses 32 bisection steps on log risk aversion. It retains the published floor $10^{-8}\\gamma_{ref}$ when the target is out of reach. A fixed risk-aversion path skips this search. This fast implementation is retained for exact parity with the published paths." ] }, { "cell_type": "code", "execution_count": 7, "id": "signal-attribution-13", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.643526Z", "iopub.status.busy": "2026-10-29T10:22:35.643467Z", "iopub.status.idle": "2026-10-29T10:22:35.647109Z", "shell.execute_reply": "2026-10-29T10:22:35.646693Z" }, "part": "optimiser" }, "outputs": [], "source": [ "def solve_diag(desk: Desk, G: np.ndarray, sigma_target: np.ndarray | None=None, gamma_path: np.ndarray | None=None, lo: np.ndarray | None=None, hi: np.ndarray | None=None, n_bisect: int=32) -> Solution:\n", " \"\"\"Run A on a batch of combined signals G (B, T, N), diagonal Sigma.\"\"\"\n", " (B, T, N) = G.shape\n", " sig2 = desk.sigma2()\n", " (lam, lam0) = (desk.lam(), desk.lam0())\n", " if lo is None or hi is None:\n", " (lo, hi) = desk.bounds()\n", " (lo, hi) = _bcast_bounds(lo, hi, B, N)\n", " x = np.zeros((B, T, N))\n", " gam = np.zeros((B, T))\n", " xp = np.zeros((B, N))\n", " g_ref = desk.gamma_ref()\n", " (g_min, g_max) = (1e-08 * g_ref, 100000000.0 * g_ref)\n", " g_prev = np.full(B, g_ref)\n", " if sigma_target is not None:\n", " tgt2 = np.broadcast_to(np.asarray(sigma_target, float), (B,)) ** 2\n", "\n", " def risk2(gm):\n", " xx = _trade_diag(G[:, t], xp, gm, sig2, lam, lam0, lo, hi)[0]\n", " return (sig2 * xx ** 2).sum(1)\n", " for t in range(T):\n", " if gamma_path is not None:\n", " gm = gamma_path[:, t]\n", " else:\n", " a = np.maximum(g_prev / 30, g_min)\n", " b = np.minimum(g_prev * 30, g_max)\n", " for _ in range(10):\n", " low_a = risk2(a) < tgt2\n", " a = np.where(low_a & (a > g_min), np.maximum(a / 1000.0, g_min), a)\n", " high_b = risk2(b) > tgt2\n", " b = np.where(high_b & (b < g_max), np.minimum(b * 1000.0, g_max), b)\n", " if not ((low_a & (a > g_min)).any() or (high_b & (b < g_max)).any()):\n", " break\n", " (la, lb) = (np.log(a), np.log(b))\n", " for _ in range(n_bisect):\n", " m = 0.5 * (la + lb)\n", " above = risk2(np.exp(m)) > tgt2\n", " la = np.where(above, m, la)\n", " lb = np.where(above, lb, m)\n", " gm = np.exp(0.5 * (la + lb))\n", " gm = np.where(risk2(a) < tgt2, a, gm)\n", " g_prev = gm\n", " x[:, t] = _trade_diag(G[:, t], xp, gm, sig2, lam, lam0, lo, hi)[0]\n", " gam[:, t] = gm\n", " xp = x[:, t]\n", " return _multipliers_diag(desk, G, x, gam, lo, hi)" ] }, { "cell_type": "markdown", "id": "signal-attribution-14", "metadata": { "part": "optimiser" }, "source": [ "`_multipliers_diag` recovers the diagonal spread subgradient and signed box multiplier from the optimality equation (section 3.3; Durin Eqs. 17, 21). These are closed-form multipliers, not outputs of CVXPY." ] }, { "cell_type": "code", "execution_count": 8, "id": "signal-attribution-15", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.648226Z", "iopub.status.busy": "2026-10-29T10:22:35.648163Z", "iopub.status.idle": "2026-10-29T10:22:35.650796Z", "shell.execute_reply": "2026-10-29T10:22:35.650251Z" }, "part": "optimiser" }, "outputs": [], "source": [ "def _multipliers_diag(desk, G, x, gam, lo, hi) -> Solution:\n", " \"\"\"Recover s and nu from the KKT equation given positions and gamma.\"\"\"\n", " sig2 = desk.sigma2()\n", " (lam, lam0) = (desk.lam(), desk.lam0())\n", " xp = np.concatenate([np.zeros_like(x[:, :1]), x[:, :-1]], axis=1)\n", " dx = x - xp\n", " g = gam[:, :, None]\n", " q = (lam + g) * sig2\n", " z = G - g * sig2 * xp\n", " dstar = np.sign(z) * np.maximum(np.abs(z) - lam0, 0.0) / q\n", " xu = xp + dstar\n", " at_bound = (xu > hi[:, None]) | (xu < lo[:, None])\n", " if lam0 > 0:\n", " s = np.where(dx != 0, np.sign(dx), np.clip(z / lam0, -1.0, 1.0))\n", " else:\n", " s = np.zeros_like(dx)\n", " nu = np.where(at_bound, z - q * dx - lam0 * s, 0.0)\n", " return Solution(x=x, dx=dx, gamma=gam, s=s, nu=nu, at_bound=at_bound)" ] }, { "cell_type": "markdown", "id": "signal-attribution-16", "metadata": { "part": "optimiser" }, "source": [ "`solve_cvxpy` solves the same static problem with a dense covariance (Durin Eqs. 5, 7, 15, 19). In risk-target mode the risk constraint's dual is $\\mu$, so $\\gamma=2\\mu$. In fixed-$\\gamma$ mode it enters the objective directly.\n", "\n", "Clarabel tolerances match the study. Trades and positions within $10^{-6}\\sigma^*/\\sqrt N$ of zero or a bound are snapped, as in the study; the raw primal values remain available. The solver's box duals remain untouched. Nonzero trades use their actual sign; at zero trades the KKT residual selects a spread subgradient in $[-1,1]$. Numerical KKT and attribution discrepancies are measured below, rather than removed by reconstructing box multipliers. This compact solver requires finite box bounds.\n", "\n", "If Clarabel fails while reusing solver state, retry the identical problem once from a cold start, with the same tolerances. Retry events and solver warnings are retained and printed; a second failure stops execution. Warning messages are recorded without local file paths.\n", "\n", "Finite bounds also provide a rigorous risk bound. Put $m_i=\\max(|lo_i|,|hi_i|)$. Every feasible position satisfies $x^\\top\\Sigma x\\leq m^\\top|\\Sigma|m$. If the requested squared risk target exceeds this bound, its constraint is redundant. We solve the same problem with $\\gamma=0$ and omit the redundant inequality, whose dual may be chosen as zero. This exact constraint elimination avoids a poorly scaled, enormous right-hand side; it does not cap targets or alter their updates.\n", "The redundant-risk check uses the bound $x^\\top\\Sigma x\\leq m^\\top|\\Sigma|m$,\n", "where $m_i=\\max(|l_i|,|u_i|)$. This inequality justifies removing the constraint;\n", "the huge-target versus zero-gamma test below checks branch equivalence, not an\n", "independent solver oracle.\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "signal-attribution-17", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.652069Z", "iopub.status.busy": "2026-10-29T10:22:35.651991Z", "iopub.status.idle": "2026-10-29T10:22:35.658855Z", "shell.execute_reply": "2026-10-29T10:22:35.658392Z" }, "part": "optimiser" }, "outputs": [], "source": [ "def solve_cvxpy(desk, G, sigma_target=None, gamma_path=None, lo=None, hi=None,\n", " snap_tolerance=1e-6):\n", " \"\"\"Dense static optimisation; retain raw primal values and box duals.\"\"\"\n", " B, T, N = G.shape\n", " Sigma = desk.covariance()\n", " lam, spread = desk.lam(), desk.lam0()\n", " if lo is None or hi is None:\n", " lo, hi = desk.bounds()\n", " lo, hi = _bcast_bounds(lo, hi, B, N)\n", " if not (np.isfinite(lo).all() and np.isfinite(hi).all()):\n", " raise ValueError(\"This compact CVXPY example requires finite position bounds.\")\n", " dxv, xv = cp.Variable(N), cp.Variable(N)\n", " linear, previous = cp.Parameter(N), cp.Parameter(N)\n", " qscale, risk_squared = cp.Parameter(nonneg=True), cp.Parameter(nonneg=True)\n", " lower, upper = cp.Parameter(N), cp.Parameter(N)\n", " upper_constraint, lower_constraint = xv <= upper, xv >= lower\n", " constraints = [xv - dxv == previous, upper_constraint, lower_constraint]\n", " risk_constraint = cp.quad_form(xv, Sigma) <= risk_squared\n", " objective = cp.Maximize(linear @ dxv - 0.5 * qscale * cp.quad_form(dxv, Sigma)\n", " - spread * cp.norm1(dxv))\n", " risk_problem = cp.Problem(objective, constraints + [risk_constraint])\n", " gamma_problem = cp.Problem(objective, constraints)\n", " shape = (B, T, N)\n", " x, raw, trades = np.zeros(shape), np.zeros(shape), np.zeros(shape)\n", " upper_dual, lower_dual = np.zeros(shape), np.zeros(shape)\n", " gamma = np.zeros((B, T))\n", " at_bound = np.zeros(shape, dtype=bool)\n", " statuses = np.empty((B, T), dtype=object)\n", " solver_iterations = np.zeros((B, T), dtype=int)\n", " equality_error = np.zeros((B, T))\n", " snap = snap_tolerance * desk.risk_target / np.sqrt(N)\n", " for b in range(B):\n", " maximum_position = np.maximum(np.abs(lo[b]), np.abs(hi[b]))\n", " variance_bound = maximum_position @ np.abs(Sigma) @ maximum_position\n", " risk_required = gamma_path is None and float(sigma_target[b])**2 < variance_bound\n", " xp = np.zeros(N)\n", " for t in range(T):\n", " previous.value, lower.value, upper.value = xp, lo[b], hi[b]\n", " if risk_required:\n", " linear.value, qscale.value = G[b, t], lam\n", " risk_squared.value = float(sigma_target[b])**2\n", " problem = risk_problem\n", " else:\n", " gm = 0.0 if gamma_path is None else float(gamma_path[b, t])\n", " linear.value = G[b, t] - gm * (Sigma @ xp)\n", " qscale.value = lam + gm\n", " problem = gamma_problem\n", " options = dict(solver=\"CLARABEL\", tol_gap_abs=1e-10,\n", " tol_gap_rel=1e-10, tol_feas=1e-10, tol_ktratio=1e-8)\n", " with warnings.catch_warnings(record=True) as caught:\n", " warnings.simplefilter(\"always\")\n", " try:\n", " problem.solve(**options)\n", " except cp.error.SolverError:\n", " solver_events.append(dict(event=\"cold retry\", run=b, period=t,\n", " risk_target=(None if gamma_path is not None\n", " else float(sigma_target[b]))))\n", " try:\n", " problem.solve(**options, warm_start=False)\n", " except cp.error.SolverError as error:\n", " target = None if gamma_path is not None else float(sigma_target[b])\n", " raise RuntimeError(\n", " f\"Clarabel failed after cold retry: run={b}, period={t}, \"\n", " f\"risk_target={target}; {error}\") from None\n", " for warning in caught:\n", " solver_events.append(dict(event=\"solver warning\", message=str(warning.message)))\n", " if problem.status not in (\"optimal\", \"optimal_inaccurate\"):\n", " raise RuntimeError(f\"Clarabel status: {problem.status}\")\n", " if risk_required:\n", " gm = 2 * float(np.asarray(risk_constraint.dual_value).item())\n", " statuses[b, t] = problem.status\n", " solver_iterations[b, t] = problem.solver_stats.num_iters\n", " equality_error[b, t] = np.max(np.abs(xv.value - dxv.value - xp))\n", " raw[b, t] = xv.value\n", " delta = np.where(np.abs(dxv.value) < snap, 0.0, dxv.value)\n", " xt = xp + delta\n", " on_hi, on_lo = np.abs(xt - hi[b]) < snap, np.abs(xt - lo[b]) < snap\n", " xt = np.where(on_hi, hi[b], np.where(on_lo, lo[b], xt))\n", " x[b, t], trades[b, t], gamma[b, t] = xt, xt - xp, gm\n", " upper_dual[b, t] = upper_constraint.dual_value\n", " lower_dual[b, t] = lower_constraint.dual_value\n", " at_bound[b, t] = on_hi | on_lo\n", " xp = xt\n", " nu = upper_dual - lower_dual\n", " # Use actual trade signs. At zero trades, choose a subgradient in [-1, 1]\n", " # from the KKT residual, leaving the solver's box duals unchanged.\n", " w = G - lam * trades @ Sigma - gamma[:, :, None] * (x @ Sigma) - nu\n", " s = (np.where(trades != 0, np.sign(trades), np.clip(w / spread, -1, 1))\n", " if spread > 0 else np.zeros_like(trades))\n", " risk = np.sqrt(np.einsum(\"bti,ij,btj->bt\", x, Sigma, x))\n", " raw_risk = np.sqrt(np.einsum(\"bti,ij,btj->bt\", raw, Sigma, raw))\n", " stationarity = w - spread * s\n", " diagnostics = dict(status=statuses, solver_iterations=solver_iterations,\n", " risk=risk, raw_risk=raw_risk, equality_error=equality_error,\n", " stationarity_error=np.max(np.abs(stationarity), axis=-1),\n", " box_violation=np.maximum(0, np.maximum(\n", " raw - hi[:, None], lo[:, None] - raw)).max(axis=-1),\n", " snap_change=np.max(np.abs(raw - x), axis=-1))\n", " return Solution(x, trades, gamma, s, nu, at_bound,\n", " upper_dual, lower_dual, raw, diagnostics)\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-18", "metadata": { "part": "optimiser" }, "source": [ "`solve` selects the diagonal closed form or CVXPY. All rerun methods below use this same entry point." ] }, { "cell_type": "code", "execution_count": 10, "id": "signal-attribution-19", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.660146Z", "iopub.status.busy": "2026-10-29T10:22:35.660080Z", "iopub.status.idle": "2026-10-29T10:22:35.661766Z", "shell.execute_reply": "2026-10-29T10:22:35.661400Z" }, "part": "optimiser" }, "outputs": [], "source": [ "def solve(desk, G, **options):\n", " return (solve_diag if desk.cov == \"diag\" else solve_cvxpy)(desk, G, **options)\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-20", "metadata": { "part": "optimiser" }, "source": [ "### Diagonal oracle\n", "\n", "This repeats the study's diagonal-oracle configuration: 8 assets, 30 periods, box 1, spread 0.8, seed 3. The declared study tolerances are positions $10^{-5}$ absolute, box multipliers $10^{-4}$ absolute, and active risk aversion $10^{-3}$ relative. CVXPY's inactive risk dual tends to zero; the closed form retains its floor, so inactive risk aversion is reported separately.\n", "\n", "The first configuration has no active risk periods. A second case keeps the same seed and dimensions but uses box 1.5 to exercise active risk duals. At a zero trade on a position bound, spread and box multipliers are not separately identified: the raw solver box dual can differ from the closed form while their combined force agrees. We print the full raw discrepancy, compare raw duals where the trade is nonzero, and separately compare the study's explicitly reconstructed spread-first convention. The raw duals used elsewhere remain unchanged.\n", "\n", "Risk aversion is also non-unique when the portfolio remains on the risk sphere without a resolved trade: the spread subgradient can offset a range of risk multipliers. We print that disagreement too. The relative gamma assertion uses active periods with a trade larger than the solver's declared snapping scale on an unbounded coordinate with nonzero position, and requires this mask to be nonempty in the box-1.5 case. This is an identifiability restriction, not a larger tolerance." ] }, { "cell_type": "code", "execution_count": 11, "id": "signal-attribution-21", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.662765Z", "iopub.status.busy": "2026-10-29T10:22:35.662699Z", "iopub.status.idle": "2026-10-29T10:22:35.772309Z", "shell.execute_reply": "2026-10-29T10:22:35.771800Z" }, "part": "optimiser" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Oracle box=1.0; active risk periods=0\n", "Unidentified box coordinates: 208\n", "Active risk periods with identifiable gamma: 0\n", "Maximum position error: 8.693970349193947e-09\n", "Raw box-dual error, all coordinates: 0.0938566251796392\n", "Raw box-dual error, identifiable coordinates: 4.006360382335927e-09\n", "Study-convention box-dual error: 4.005961229402999e-09\n", "Maximum fixed-gamma position error: 1.9367057402241983e-09\n", "Maximum raw-to-snapped position change: 8.329894583113052e-09\n", "Inactive gamma, CVXPY max / diagonal floor: 0.0 3.2756529752049163e-09\n", "Oracle box=1.5; active risk periods=26\n", "Unidentified box coordinates: 152\n", "Active risk periods with identifiable gamma: 22\n", "Raw active gamma max absolute difference (including periods not identified by this test): 0.07429351985869384\n", "Maximum position error: 6.956109848826886e-07\n", "Raw box-dual error, all coordinates: 0.10709899976812781\n", "Raw box-dual error, identifiable coordinates: 8.28780995165701e-07\n", "Study-convention box-dual error: 6.768218758224087e-07\n", "Maximum fixed-gamma position error: 7.985581536917863e-06\n", "Maximum raw-to-snapped position change: 1.2514183911349619e-08\n", "Identifiable active gamma max relative error: 7.426598902293691e-06\n", "Inactive gamma, CVXPY max / diagonal floor: 1.0608439388721236e-08 3.2756529752049163e-09\n", "PASS: positions, identifiable active risk duals, identifiable raw box duals and study-convention box duals\n", "Solver events: [{'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}, {'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}]\n" ] } ], "source": [ "for box in (1.0, 1.5):\n", " oracle = Desk(n_assets=8, n_periods=30, box=box, spread_rel=0.8)\n", " g, r = simulate(oracle, 3)\n", " G = g.sum(0)[None]\n", " diagonal = solve_diag(oracle, G, sigma_target=np.array([1.0]))\n", " convex = solve_cvxpy(oracle, G, sigma_target=np.array([1.0]))\n", " active = diagonal.gamma[0] > 1e-6 * oracle.gamma_ref()\n", " resolved_trade = np.abs(diagonal.dx) > 1e-6 * oracle.risk_target / np.sqrt(oracle.n_assets)\n", " identifiable = resolved_trade\n", " risk_identifiable = active & np.any(\n", " resolved_trade[0] & ~diagonal.at_bound[0]\n", " & (np.abs(diagonal.x[0]) > 1e-6 * oracle.risk_target / np.sqrt(oracle.n_assets)), axis=1)\n", " assert identifiable.any()\n", " np.testing.assert_allclose(convex.x, diagonal.x, atol=1e-5)\n", " np.testing.assert_allclose(convex.nu[identifiable], diagonal.nu[identifiable], atol=1e-4)\n", " np.testing.assert_allclose(convex.gamma[0, risk_identifiable], diagonal.gamma[0, risk_identifiable], rtol=1e-3)\n", " if box > 1:\n", " assert risk_identifiable.any(), \"Risk-dual comparison must be non-vacuous\"\n", " # Explicitly reproduce the study's spread-first convention for comparison.\n", " # This diagnostic does not overwrite convex.nu or its raw dual arrays.\n", " w = G - oracle.lam() * convex.dx - convex.gamma[:, :, None] * convex.x\n", " convention_s = np.where(convex.dx != 0, np.sign(convex.dx),\n", " np.clip(w / oracle.lam0(), -1, 1))\n", " convention_nu = np.where(convex.at_bound, w - oracle.lam0() * convention_s, 0)\n", " np.testing.assert_allclose(convention_nu, diagonal.nu, atol=1e-4)\n", " fixed = solve_cvxpy(oracle, G, gamma_path=diagonal.gamma)\n", " np.testing.assert_allclose(fixed.x, diagonal.x, atol=1e-5)\n", " print(f\"Oracle box={box}; active risk periods={active.sum()}\")\n", " print(\"Unidentified box coordinates:\", int((~identifiable).sum()))\n", " print(\"Active risk periods with identifiable gamma:\", int(risk_identifiable.sum()))\n", " if active.any():\n", " print(\"Raw active gamma max absolute difference (including periods not identified by this test):\",\n", " np.max(np.abs(convex.gamma[0, active] - diagonal.gamma[0, active])))\n", " print(\"Maximum position error:\", np.max(np.abs(convex.x - diagonal.x)))\n", " print(\"Raw box-dual error, all coordinates:\", np.max(np.abs(convex.nu - diagonal.nu)))\n", " print(\"Raw box-dual error, identifiable coordinates:\",\n", " np.max(np.abs(convex.nu[identifiable] - diagonal.nu[identifiable])))\n", " print(\"Study-convention box-dual error:\", np.max(np.abs(convention_nu - diagonal.nu)))\n", " print(\"Maximum fixed-gamma position error:\", np.max(np.abs(fixed.x - diagonal.x)))\n", " print(\"Maximum raw-to-snapped position change:\", np.max(np.abs(convex.primal_raw - convex.x)))\n", " if risk_identifiable.any():\n", " print(\"Identifiable active gamma max relative error:\", np.max(np.abs(\n", " convex.gamma[0, risk_identifiable] / diagonal.gamma[0, risk_identifiable] - 1)))\n", " if (~active).any():\n", " print(\"Inactive gamma, CVXPY max / diagonal floor:\",\n", " convex.gamma[0, ~active].max(), 1e-8 * oracle.gamma_ref())\n", "print(\"PASS: positions, identifiable active risk duals, identifiable raw box duals and study-convention box duals\")\n", "print(\"Solver events:\", solver_events)\n", "solver_events.clear()\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-22", "metadata": { "part": "durin" }, "source": [ "## 2. Durin attribution\n", "\n", "`effective_position_matrix` builds the general rank-one bound terms from vectors $v_c$. Write a position constraint as $\\epsilon_c v_c^\\top x\\leq\\epsilon_c M_c$, with signed bound $M_c$, orientation $\\epsilon_c\\in\\{-1,1\\}$ and non-negative solver dual $\\lambda^c$. Then\n", "\n", "$$\\eta_c=\\frac{\\epsilon_c\\lambda^c}{2M_c},\\qquad D_b=\\sum_c2\\eta_c v_cv_c^\\top.$$\n", "\n", "This is section 3.3's signed-bound form of Durin Eq. 37. For the coordinate box, $v_c$ is a coordinate vector and the oriented dual is `nu`. Zero bounds cannot use this division; their limiting condition is imposed in `durin_matrix`. The matrix helper accepts general constraint directions, while this notebook's optimiser supplies coordinate bounds." ] }, { "cell_type": "code", "execution_count": 12, "id": "signal-attribution-23", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.773617Z", "iopub.status.busy": "2026-10-29T10:22:35.773538Z", "iopub.status.idle": "2026-10-29T10:22:35.775335Z", "shell.execute_reply": "2026-10-29T10:22:35.774974Z" }, "part": "durin" }, "outputs": [], "source": [ "def effective_position_matrix(vectors, signed_bounds, oriented_duals):\n", " \"\"\"Finite bound penalties: 2 eta_c v_c v_c.T; zero bounds handled separately.\"\"\"\n", " eta = oriented_duals / (2 * signed_bounds)\n", " return vectors.T @ ((2 * eta)[:, None] * vectors)\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-24", "metadata": { "part": "durin" }, "source": [ "`durin_matrix` forms the full matrices $\\bar Q=(\\lambda+\\gamma)\\Sigma+D_b+D_s$ and $\\bar P=\\gamma\\Sigma+D_b$ (Durin Eqs. 30–31). On a nonzero joint trade, $\\eta_s=\\lambda_0/(2|\\Delta x|)$ and $D_s=\\operatorname{diag}(2\\eta_s)$ (Eq. 48).\n", "\n", "Each signal solves $\\bar Q\\Delta x^k=g^k-\\bar P x^k_{t-1}$. Zero joint trades impose $\\Delta x^k=0$ (Eq. 54); a binding zero position imposes $x^k=0$ (Eq. 44). The remaining coordinates solve the reduced matrix equation, including their covariance with the fixed coordinates. This is a matrix solve, not an asset-by-asset diagonal formula.\n", "\n", "Spread penalties act on trades and enter only $\\bar Q$; position penalties act on the new position and enter both $\\bar Q$ and $\\bar P$.\n", "\n", "The zero-trade limit applies only with positive spread. With zero spread, opposite component trades may cancel; the ordinary matrix recursion must retain them." ] }, { "cell_type": "code", "execution_count": 13, "id": "signal-attribution-25", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.776554Z", "iopub.status.busy": "2026-10-29T10:22:35.776489Z", "iopub.status.idle": "2026-10-29T10:22:35.779411Z", "shell.execute_reply": "2026-10-29T10:22:35.778983Z" }, "part": "durin" }, "outputs": [], "source": [ "def durin_matrix(desk, sol, g, b=0):\n", " K, T, N = g.shape\n", " Sigma, lam, spread = desk.covariance(), desk.lam(), desk.lam0()\n", " xk, xp = np.zeros((K, T, N)), np.zeros((K, N))\n", " coordinate_vectors = np.eye(N)\n", " for t in range(T):\n", " x, dx, nu = sol.x[b, t], sol.dx[b, t], sol.nu[b, t]\n", " bound, gamma = sol.at_bound[b, t], sol.gamma[b, t]\n", " zero_trade = (dx == 0) & (spread > 0)\n", " zero_position = bound & (x == 0) & ~zero_trade\n", " free = ~(zero_trade | zero_position)\n", " finite_bound = free & bound\n", " position_penalty = effective_position_matrix(\n", " coordinate_vectors[finite_bound], x[finite_bound], nu[finite_bound])\n", " eta_spread = np.zeros(N)\n", " moving = free & (dx != 0)\n", " eta_spread[moving] = spread / (2 * np.abs(dx[moving]))\n", " Qbar = (lam + gamma) * Sigma + position_penalty + np.diag(2 * eta_spread)\n", " Pbar = gamma * Sigma + position_penalty\n", " delta = np.where(zero_position, -xp, 0.0)\n", " indices = np.flatnonzero(free)\n", " rhs = g[:, t, indices] - xp @ Pbar[indices].T - delta @ Qbar[indices].T\n", " delta[:, indices] = np.linalg.solve(Qbar[np.ix_(indices, indices)], rhs.T).T\n", " xk[:, t] = xp + delta\n", " xp = xk[:, t]\n", " return xk, np.diff(xk, axis=1, prepend=0.0)\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-26", "metadata": { "part": "durin" }, "source": [ "`durin_positions` retains the diagonal specialisation for the published-path check and calls the full matrix recursion for dense covariance. Its spread and bound penalties are the same effective multipliers (section 3.3; Durin Eqs. 37, 48)." ] }, { "cell_type": "code", "execution_count": 14, "id": "signal-attribution-27", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.780728Z", "iopub.status.busy": "2026-10-29T10:22:35.780666Z", "iopub.status.idle": "2026-10-29T10:22:35.783682Z", "shell.execute_reply": "2026-10-29T10:22:35.783360Z" }, "part": "durin" }, "outputs": [], "source": [ "def durin_positions(desk: Desk, sol: Solution, g: np.ndarray, b: int=0):\n", " \"\"\"Signal-wise positions x^k and trades dx^k of run ``b`` of ``sol``.\"\"\"\n", " if desk.cov != \"diag\":\n", " return durin_matrix(desk, sol, g, b)\n", " (K, T, N) = g.shape\n", " sig2 = desk.sigma2()\n", " (lam, lam0) = (desk.lam(), desk.lam0())\n", " xk = np.zeros((K, T, N))\n", " xp = np.zeros((K, N))\n", " for t in range(T):\n", " (x, dx, nu, ab, gm) = (sol.x[b, t], sol.dx[b, t], sol.nu[b, t], sol.at_bound[b, t], sol.gamma[b, t])\n", " zero_trade = (dx == 0) & (lam0 > 0)\n", " zero_pos = ab & (x == 0) & ~zero_trade\n", " free = ~zero_trade & ~zero_pos\n", " s = sol.s[b, t]\n", " with np.errstate(divide='ignore', invalid='ignore'):\n", " Ds = np.divide(lam0 * s, dx, out=np.zeros_like(dx), where=free & (dx != 0))\n", " Db = np.where(free & ab, nu / x, 0.0)\n", " Qbar = (lam + gm) * sig2 + Ds + Db\n", " Pbar = gm * sig2 + Db\n", " dk = (g[:, t] - Pbar * xp) / Qbar\n", " dk = np.where(zero_trade, 0.0, np.where(zero_pos, -xp, dk))\n", " xk[:, t] = xp + dk\n", " xp = xk[:, t]\n", " dxk = np.diff(xk, axis=1, prepend=0.0)\n", " return (xk, dxk)" ] }, { "cell_type": "markdown", "id": "signal-attribution-28", "metadata": { "part": "durin" }, "source": [ "`allocations` returns gross P&L, quadratic cost, spread cost, net P&L and, when given the total position, the contribution $x_k^\\top\\Sigma x$ to portfolio **variance**. Costs use the joint trade: $\\tfrac12\\Delta x_k^\\top\\Lambda\\Delta x$ and $\\lambda_0\\Delta x_k^\\top\\operatorname{sgn}(\\Delta x)$ (Durin Eqs. 13–14, 22). These signed contributions can be negative; their sums are checked separately." ] }, { "cell_type": "code", "execution_count": 15, "id": "signal-attribution-29", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.784889Z", "iopub.status.busy": "2026-10-29T10:22:35.784827Z", "iopub.status.idle": "2026-10-29T10:22:35.786875Z", "shell.execute_reply": "2026-10-29T10:22:35.786536Z" }, "part": "durin" }, "outputs": [], "source": [ "def allocations(desk, xk, dxk, dx_total, returns, x_total=None):\n", " Sigma = desk.covariance()\n", " gross = (xk * returns).sum(-1)\n", " quadratic_cost = 0.5 * desk.lam() * (dxk * (dx_total @ Sigma)).sum(-1)\n", " spread_cost = desk.lam0() * (dxk * np.sign(dx_total)).sum(-1)\n", " result = dict(gross=gross, quadratic_cost=quadratic_cost, spread_cost=spread_cost,\n", " net=gross - quadratic_cost - spread_cost)\n", " if x_total is not None:\n", " result[\"variance\"] = (xk * (x_total @ Sigma)).sum(-1)\n", " return result\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-30", "metadata": { "part": "durin" }, "source": [ "`pnl_split` applies that cost convention to component positions. The diagonal arithmetic is preserved to reproduce the published endpoints exactly." ] }, { "cell_type": "code", "execution_count": 16, "id": "signal-attribution-31", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.788045Z", "iopub.status.busy": "2026-10-29T10:22:35.787970Z", "iopub.status.idle": "2026-10-29T10:22:35.790051Z", "shell.execute_reply": "2026-10-29T10:22:35.789659Z" }, "part": "durin" }, "outputs": [], "source": [ "def pnl_split(desk, xk, dxk, dx_total, r):\n", " \"\"\"P&L of component positions xk with costs split as Durin Eqs. 14 and 22.\"\"\"\n", " if desk.cov != \"diag\":\n", " return allocations(desk, xk, dxk, dx_total, r)[\"net\"]\n", " (lam, lam0) = (desk.lam(), desk.lam0())\n", " sig2 = desk.sigma2()\n", " gross = (xk * r).sum(-1)\n", " quad = 0.5 * lam * (dxk * sig2 * dx_total).sum(-1)\n", " spread = lam0 * (dxk * np.sign(dx_total)).sum(-1)\n", " return gross - quad - spread" ] }, { "cell_type": "markdown", "id": "signal-attribution-32", "metadata": { "part": "durin" }, "source": [ "`pnl` evaluates the full run with the same position, quadratic-cost and spread-cost ledger (Durin Eq. 5 and section 3.3)." ] }, { "cell_type": "code", "execution_count": 17, "id": "signal-attribution-33", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.791152Z", "iopub.status.busy": "2026-10-29T10:22:35.791087Z", "iopub.status.idle": "2026-10-29T10:22:35.792895Z", "shell.execute_reply": "2026-10-29T10:22:35.792536Z" }, "part": "durin" }, "outputs": [], "source": [ "def pnl(desk: Desk, sol: Solution, r: np.ndarray) -> np.ndarray:\n", " \"\"\"Net P&L per period (B, T): x_t . r_{t+1} minus quadratic and spread costs.\"\"\"\n", " return pnl_split(desk, sol.x, sol.dx, sol.dx, r)" ] }, { "cell_type": "markdown", "id": "signal-attribution-34", "metadata": { "part": "durin" }, "source": [ "`component_pnl` evaluates signal or constraint portfolios using the joint run's trades for the cost allocation (Durin Eqs. 14, 22)." ] }, { "cell_type": "code", "execution_count": 18, "id": "signal-attribution-35", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.793870Z", "iopub.status.busy": "2026-10-29T10:22:35.793819Z", "iopub.status.idle": "2026-10-29T10:22:35.795398Z", "shell.execute_reply": "2026-10-29T10:22:35.795044Z" }, "part": "durin" }, "outputs": [], "source": [ "def component_pnl(desk, xk, dxk, sol: Solution, r, b: int=0):\n", " \"\"\"P&L per component with Durin's cost split (Eqs. 14, 22).\"\"\"\n", " return pnl_split(desk, xk, dxk, sol.dx[b][None], r[None])" ] }, { "cell_type": "markdown", "id": "signal-attribution-36", "metadata": { "part": "other" }, "source": [ "## 3. Other methods\n", "\n", "`grinold_easton_positions` sends each signal, the spread force $-\\lambda_0s$ and the box force $-\\nu$ through $Q\\Delta x+Px_{t-1}=\\text{source}$ (section 3.3; Durin Eqs. 10–12, 17, 21). It keeps separate spread and position-limit portfolios and uses full matrices for dense covariance." ] }, { "cell_type": "code", "execution_count": 19, "id": "signal-attribution-37", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.796640Z", "iopub.status.busy": "2026-10-29T10:22:35.796578Z", "iopub.status.idle": "2026-10-29T10:22:35.798717Z", "shell.execute_reply": "2026-10-29T10:22:35.798390Z" }, "part": "other" }, "outputs": [], "source": [ "def grinold_easton_positions(desk, sol, g):\n", " K, T, N = g.shape\n", " lam, spread = desk.lam(), desk.lam0()\n", " sources = np.concatenate([g, -spread * sol.s[0][None], -sol.nu[0][None]])\n", " xk, xp = np.zeros((K + 2, T, N)), np.zeros((K + 2, N))\n", " Sigma = desk.covariance()\n", " for t in range(T):\n", " gamma = sol.gamma[0, t]\n", " if desk.cov == \"diag\":\n", " sig2 = desk.sigma2()\n", " delta = (sources[:, t] - gamma * sig2 * xp) / ((lam + gamma) * sig2)\n", " else:\n", " Q, P = (lam + gamma) * Sigma, gamma * Sigma\n", " delta = np.linalg.solve(Q, (sources[:, t] - xp @ P).T).T\n", " xk[:, t] = xp + delta\n", " xp = xk[:, t]\n", " return xk, np.diff(xk, axis=1, prepend=0.0)\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-38", "metadata": { "part": "other" }, "source": [ "`regress` fits non-negative least squares without an intercept to stand-alone net P&Ls (section 3.1; Durin's discussion on pp. 2–3). The unexplained book P&L remains a separate residual." ] }, { "cell_type": "code", "execution_count": 20, "id": "signal-attribution-39", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.799841Z", "iopub.status.busy": "2026-10-29T10:22:35.799780Z", "iopub.status.idle": "2026-10-29T10:22:35.801334Z", "shell.execute_reply": "2026-10-29T10:22:35.800967Z" }, "part": "other" }, "outputs": [], "source": [ "def regress(total, proxies):\n", " beta = nnls(proxies.T, total)[0]\n", " attributed = beta[:, None] * proxies\n", " return beta, attributed, total - attributed.sum(0)" ] }, { "cell_type": "markdown", "id": "signal-attribution-40", "metadata": { "part": "other" }, "source": [ "`iterated_proxies` reruns each signal at its updated risk target (section 3.1). Targets are multiplied by their NNLS coefficients without damping. Live zero coefficients drop proxies; zero-P&L proxies leave the convergence criterion. No live proxy gives `no_live_proxy`, not convergence. Convergence requires live coefficients within $10^{-3}$ of one, with at most 50 fits.\n", "With `trace=True`, each fit retains its targets, coefficients, realised risk paths,\n", "binding fractions (absolute tolerance $10^{-7}$, relative $10^{-6}$), proxy P&Ls\n", "and solver diagnostics. Dropped proxies have zero risk and no solver record.\n", "`update_power=1` is the method above; a power below one is a separately labelled\n", "experimental damped update, $\\sigma_k\\leftarrow\\sigma_k\\hat\\beta_k^p$.\n", "It does not change the coefficient convergence test.\n" ] }, { "cell_type": "code", "execution_count": 21, "id": "signal-attribution-41", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.802277Z", "iopub.status.busy": "2026-10-29T10:22:35.802212Z", "iopub.status.idle": "2026-10-29T10:22:35.805094Z", "shell.execute_reply": "2026-10-29T10:22:35.804794Z" }, "part": "other" }, "outputs": [], "source": [ "def iterated_proxies(desk, g, total, r, tol=1e-3, max_iter=50, *, trace=False, update_power=1.0):\n", " if not 0 < update_power <= 1:\n", " raise ValueError(\"update_power must be in (0, 1]; 1 is the undamped method\")\n", " history = []\n", " K, T, N = g.shape\n", " sigma = np.full(K, desk.risk_target)\n", " status = \"max_iter\"\n", " for iteration in range(1, max_iter + 1):\n", " active = np.flatnonzero(sigma > 0)\n", " proxies = np.zeros((K, T))\n", " if len(active):\n", " sol = solve(desk, g[active], sigma_target=sigma[active])\n", " proxies[active] = pnl(desk, sol, r)\n", " solved = active.copy()\n", " active = active[np.any(proxies[active] != 0, axis=1)]\n", " beta = np.zeros(K)\n", " if len(active):\n", " beta[active], _, _ = regress(total, proxies[active])\n", " if trace:\n", " risk = np.zeros((K, T))\n", " if len(solved):\n", " risk[solved] = np.sqrt(np.einsum(\n", " \"kti,ij,ktj->kt\", sol.x, desk.covariance(), sol.x))\n", " history.append(dict(iteration=iteration, beta=beta.copy(),\n", " risk_targets=sigma.copy(), realised_risk=risk,\n", " binding_fraction=np.mean(np.isclose(risk, sigma[:, None],\n", " atol=1e-7, rtol=1e-6) & (sigma[:, None] > 0), axis=1),\n", " solved=solved, proxies=proxies.copy(),\n", " solver=sol.diagnostics if len(solved) else None))\n", " if not len(active):\n", " status = \"no_live_proxy\"\n", " break\n", " if np.max(np.abs(beta[active] - 1)) < tol:\n", " status = \"converged\"\n", " break\n", " if iteration == max_iter:\n", " break\n", " sigma[active] *= beta[active] ** update_power\n", " return beta[:, None] * proxies, dict(\n", " status=status, n_iter=iteration, beta=beta, risk_targets=sigma, history=history)" ] }, { "cell_type": "markdown", "id": "signal-attribution-42", "metadata": { "part": "other" }, "source": [ "`coalition_signals` sums forecasts for every non-empty subset, indexed by bitmask (section 3.4's game). Every subset uses the same construction rules." ] }, { "cell_type": "code", "execution_count": 22, "id": "signal-attribution-43", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.806212Z", "iopub.status.busy": "2026-10-29T10:22:35.806143Z", "iopub.status.idle": "2026-10-29T10:22:35.807970Z", "shell.execute_reply": "2026-10-29T10:22:35.807529Z" }, "part": "other" }, "outputs": [], "source": [ "def coalition_signals(g: np.ndarray) -> np.ndarray:\n", " \"\"\"Combined signal of every coalition, indexed by bitmask 1..2^K-1.\"\"\"\n", " K = g.shape[0]\n", " masks = np.arange(1, 2 ** K)\n", " member = masks[:, None] >> np.arange(K) & 1\n", " return np.einsum('mk,ktn->mtn', member.astype(float), g)" ] }, { "cell_type": "markdown", "id": "signal-attribution-44", "metadata": { "part": "other" }, "source": [ "`shapley_weights` computes the probability weight of each coalition size in a random ordering of signals (section 3.4)." ] }, { "cell_type": "code", "execution_count": 23, "id": "signal-attribution-45", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.808884Z", "iopub.status.busy": "2026-10-29T10:22:35.808827Z", "iopub.status.idle": "2026-10-29T10:22:35.810500Z", "shell.execute_reply": "2026-10-29T10:22:35.810130Z" }, "part": "other" }, "outputs": [], "source": [ "def shapley_weights(K: int):\n", " return np.array([factorial(s) * factorial(K - s - 1) / factorial(K) for s in range(K)])" ] }, { "cell_type": "markdown", "id": "signal-attribution-46", "metadata": { "part": "other" }, "source": [ "`shapley` averages marginal P&L over all orders. The empty coalition's value is zero; this calculation is applied period by period (section 3.4)." ] }, { "cell_type": "code", "execution_count": 24, "id": "signal-attribution-47", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.811477Z", "iopub.status.busy": "2026-10-29T10:22:35.811418Z", "iopub.status.idle": "2026-10-29T10:22:35.813497Z", "shell.execute_reply": "2026-10-29T10:22:35.813030Z" }, "part": "other" }, "outputs": [], "source": [ "def shapley(v: np.ndarray, K: int) -> np.ndarray:\n", " \"\"\"Exact Shapley values per period.\"\"\"\n", " w = shapley_weights(K)\n", " masks = np.arange(2 ** K)\n", " size = np.array([bin(m).count('1') for m in masks])\n", " phi = np.zeros((K, v.shape[1]))\n", " for k in range(K):\n", " S = masks[masks >> k & 1 == 0]\n", " phi[k] = (w[size[S]][:, None] * (v[S | 1 << k] - v[S])).sum(0)\n", " return phi" ] }, { "cell_type": "markdown", "id": "signal-attribution-48", "metadata": { "part": "other" }, "source": [ "`singleton_index` identifies the stand-alone coalitions, so the first NNLS fit and the risk-target coalition game use exactly the same paths (sections 3.1 and 3.4)." ] }, { "cell_type": "code", "execution_count": 25, "id": "signal-attribution-49", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.814346Z", "iopub.status.busy": "2026-10-29T10:22:35.814289Z", "iopub.status.idle": "2026-10-29T10:22:35.815777Z", "shell.execute_reply": "2026-10-29T10:22:35.815474Z" }, "part": "other" }, "outputs": [], "source": [ "def singleton_index(K: int) -> np.ndarray:\n", " \"\"\"Bitmask of each singleton coalition.\"\"\"\n", " return 1 << np.arange(K)" ] }, { "cell_type": "markdown", "id": "signal-attribution-50", "metadata": { "part": "other" }, "source": [ "`compare` calculates all seven displayed methods. It evaluates both coalition risk conventions, computes leave-one-out as the full-coalition value minus each complementary coalition, and keeps proxy and leave-one-out residuals separate. The full-portfolio solver and the rerun solver are the same." ] }, { "cell_type": "code", "execution_count": 26, "id": "signal-attribution-51", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.816813Z", "iopub.status.busy": "2026-10-29T10:22:35.816755Z", "iopub.status.idle": "2026-10-29T10:22:35.819421Z", "shell.execute_reply": "2026-10-29T10:22:35.819125Z" }, "part": "other" }, "outputs": [], "source": [ "def compare(desk, seed=20261029, include_iterated=True):\n", " g, r = simulate(desk, seed)\n", " K, T, _ = g.shape\n", " full = solve(desk, g.sum(0)[None], sigma_target=np.array([desk.risk_target]))\n", " total = pnl(desk, full, r)[0]\n", " parts = {}\n", " for name, positions in [(\"durin\", durin_positions),\n", " (\"grinold_easton\", grinold_easton_positions)]:\n", " x, dx = positions(desk, full, g)\n", " parts[name] = component_pnl(desk, x, dx, full, r)\n", " coalitions = coalition_signals(g)\n", " for convention in (\"target\", \"gamma\"):\n", " options = (dict(sigma_target=np.full(len(coalitions), desk.risk_target))\n", " if convention == \"target\" else\n", " dict(gamma_path=np.broadcast_to(full.gamma, (len(coalitions), T))))\n", " sol = solve(desk, coalitions, **options)\n", " v = np.vstack([np.zeros(T), pnl(desk, sol, r)])\n", " parts[\"shapley_\" + convention] = shapley(v, K)\n", " if convention == \"target\":\n", " _, parts[\"regression\"], _ = regress(total, v[singleton_index(K)])\n", " parts[\"loo\"] = v[-1] - v[(2**K - 1) ^ singleton_index(K)]\n", " iteration = None\n", " if include_iterated:\n", " parts[\"emmanuel\"], iteration = iterated_proxies(desk, g, total, r)\n", " for method in (\"regression\", \"loo\", \"emmanuel\"):\n", " if method not in parts:\n", " continue\n", " parts[method] = np.vstack([parts[method], total - parts[method].sum(0)])\n", " return total, parts, iteration" ] }, { "cell_type": "markdown", "id": "signal-attribution-52", "metadata": { "part": "reproduction" }, "source": [ "## 4. Reproduction and checks\n", "\n", "First run a small one-factor example (8 assets, 30 periods, seed 0), using CVXPY only. Long-only bounds exercise the zero-position limit as well as the no-trade limit. Raw duals are used, so the declared absolute tolerance is $10^{-5}$ for dense attribution sums, consistent with the diagonal oracle's position tolerance. This is a numerical reconciliation check, not the exact endpoint check below. The output gives the errors in positions, trades, gross P&L, each cost, net P&L and variance, plus the KKT residual.\n", "\n", "The separate four-asset/four-period example also runs the methods through `compare`. Its iterated targets can exceed the maximum risk permitted by the box. The solver removes such redundant risk inequalities using the rigorous bound above; a separate check compares this branch with an explicit zero-risk-aversion run. Iteration status is printed as returned, including `max_iter` when applicable. Default `compare` raises any remaining solver failure." ] }, { "cell_type": "code", "execution_count": 27, "id": "signal-attribution-53", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:35.820551Z", "iopub.status.busy": "2026-10-29T10:22:35.820481Z", "iopub.status.idle": "2026-10-29T10:22:36.577219Z", "shell.execute_reply": "2026-10-29T10:22:36.576686Z" }, "part": "reproduction" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dense KKT max residual: 4.401681470917254e-09\n", "Raw box-dual identity error: 0.0\n", "Actual liquidations from positive position to zero: 8\n", "Zero joint trades / active zero bounds: 203 118\n", "Durin {'positions': 1.9557196417885336e-09, 'trades': 1.9557196417885336e-09, 'gross': 3.085415589687557e-09, 'quadratic_cost': 5.83463121567318e-11, 'spread_cost': 1.6965026605753053e-10, 'net': 3.085415589687557e-09, 'variance': 1.947450312123067e-09}\n", "Grinold–Easton {'positions': 2.2032357271939418e-08, 'trades': 1.762616990905186e-08, 'gross': 2.7894036414011225e-08, 'quadratic_cost': 9.632457158292862e-11, 'spread_cost': 1.2902979899809264e-09, 'net': 2.7888603204573315e-08, 'variance': 6.229143978409013e-09}\n", "PASS: dense allocations, variance sums and zero-trade/zero-bound limits\n", "PASS: zero-spread cancelling trades, diagonal and dense\n", "PASS: redundant risk target equals gamma=0; max position error 0.0\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Dense method check: durin max period reconciliation error 3.3075678536331754e-07\n", "Dense method check: grinold_easton max period reconciliation error 1.5779620150813134e-07\n", "Dense method check: shapley_target max period reconciliation error 2.6728619317850644e-14\n", "Dense method check: regression max period reconciliation error 1.3877787807814457e-17\n", "Dense method check: loo max period reconciliation error 0.0\n", "Dense method check: shapley_gamma max period reconciliation error 3.6393206892526564e-08\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Dense iteration status: max_iter fits: 50\n", "Solver events: [{'event': 'solver warning', 'message': 'Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.'}]\n" ] } ], "source": [ "dense = Desk(n_assets=8, n_periods=30, cov=\"factor\", long_only=True,\n", " box=1.0, spread_rel=0.8)\n", "g, r = simulate(dense, 0)\n", "full = solve_cvxpy(dense, g.sum(0)[None], sigma_target=np.array([1.0]))\n", "Sigma = dense.covariance()\n", "kkt = (g.sum(0) - dense.lam() * full.dx[0] @ Sigma\n", " - full.gamma[0, :, None] * (full.x[0] @ Sigma)\n", " - dense.lam0() * full.s[0] - full.nu[0])\n", "print(\"Dense KKT max residual:\", np.max(np.abs(kkt)))\n", "print(\"Raw box-dual identity error:\",\n", " np.max(np.abs(full.nu - (full.dual_upper - full.dual_lower))))\n", "zero_trades = full.dx[0] == 0\n", "zero_bounds = full.at_bound[0] & (full.x[0] == 0)\n", "previous = np.concatenate([np.zeros_like(full.x[0, :1]), full.x[0, :-1]])\n", "liquidations = zero_bounds & (previous > 0) & ~zero_trades\n", "assert zero_trades.any() and zero_bounds.any() and liquidations.any()\n", "print(\"Actual liquidations from positive position to zero:\", liquidations.sum())\n", "print(\"Zero joint trades / active zero bounds:\", zero_trades.sum(), zero_bounds.sum())\n", "ledger = allocations(dense, full.x[0], full.dx[0], full.dx[0], r, full.x[0])\n", "for name, positions in [(\"Durin\", durin_positions), (\"Grinold–Easton\", grinold_easton_positions)]:\n", " xk, dxk = positions(dense, full, g)\n", " np.testing.assert_allclose(component_pnl(dense, xk, dxk, full, r),\n", " allocations(dense, xk, dxk, full.dx[0], r)[\"net\"],\n", " atol=1e-12, rtol=0)\n", " a = allocations(dense, xk, dxk, full.dx[0], r, full.x[0])\n", " errors = dict(positions=np.max(np.abs(xk.sum(0) - full.x[0])),\n", " trades=np.max(np.abs(dxk.sum(0) - full.dx[0])))\n", " errors.update({key: np.max(np.abs(values.sum(0) - ledger[key]))\n", " for key, values in a.items()})\n", " if name == \"Durin\":\n", " np.testing.assert_array_equal(dxk[:, zero_trades], 0)\n", " np.testing.assert_array_equal(xk[:, zero_bounds], 0)\n", " for key, error in errors.items():\n", " assert error <= 1e-5, (name, key, error)\n", " print(name, {key: float(value) for key, value in errors.items()})\n", "print(\"PASS: dense allocations, variance sums and zero-trade/zero-bound limits\")\n", "\n", "# Zero-spread cancellation: zero joint trade must allow opposite component trades.\n", "for covariance in (\"diag\", \"factor\"):\n", " cancellation = Desk(n_assets=2, n_periods=2, cov=covariance, spread_rel=0.0)\n", " opposing = np.array([[[0.05, -0.03], [0.04, -0.02]],\n", " [[-0.05, 0.03], [-0.04, 0.02]]])\n", " joint = solve(cancellation, opposing.sum(0)[None], gamma_path=np.ones((1, 2)))\n", " xk, dxk = durin_positions(cancellation, joint, opposing)\n", " assert np.any(np.abs(dxk) > 1e-3)\n", " np.testing.assert_allclose(dxk.sum(0), joint.dx[0], atol=1e-8, rtol=0)\n", " np.testing.assert_allclose(xk.sum(0), joint.x[0], atol=1e-8, rtol=0)\n", "print(\"PASS: zero-spread cancelling trades, diagonal and dense\")\n", "\n", "# Check the redundant-risk branch against the same problem at gamma=0.\n", "small_dense = Desk(n_assets=4, n_periods=4, cov=\"factor\")\n", "g_small, r_small = simulate(small_dense, 0)\n", "large_target = solve_cvxpy(small_dense, g_small.sum(0)[None],\n", " sigma_target=np.array([11881.909026551642]))\n", "zero_gamma = solve_cvxpy(small_dense, g_small.sum(0)[None],\n", " gamma_path=np.zeros((1, small_dense.n_periods)))\n", "np.testing.assert_array_equal(large_target.gamma, 0)\n", "np.testing.assert_allclose(large_target.x, zero_gamma.x, atol=1e-8, rtol=0)\n", "np.testing.assert_allclose(large_target.nu, zero_gamma.nu, atol=1e-8, rtol=0)\n", "print(\"PASS: redundant risk target equals gamma=0; max position error\",\n", " float(np.max(np.abs(large_target.x - zero_gamma.x))))\n", "total_dense, methods_dense, _ = compare(small_dense, seed=0, include_iterated=False)\n", "for method, values in methods_dense.items():\n", " assert np.isfinite(values).all()\n", " error = np.max(np.abs(values.sum(0) - total_dense))\n", " assert error <= 1e-5, (method, error)\n", " print(\"Dense method check:\", method, \"max period reconciliation error\", float(error))\n", "g_small, r_small = simulate(small_dense, 0)\n", "try:\n", " attributed_dense, iteration_dense = iterated_proxies(small_dense, g_small, total_dense, r_small, trace=True)\n", " residual_dense = total_dense - attributed_dense.sum(0)\n", " assert np.isfinite(attributed_dense).all()\n", " np.testing.assert_allclose(attributed_dense.sum(0) + residual_dense,\n", " total_dense, atol=1e-5, rtol=0)\n", " print(\"Dense iteration status:\", iteration_dense[\"status\"], \"fits:\", iteration_dense[\"n_iter\"])\n", "except RuntimeError as error:\n", " print(\"DENSE ITERATION NOT VERIFIED:\", str(error))\n", "print(\"Solver events:\", solver_events)\n", "solver_events.clear()\n" ] }, { "cell_type": "markdown", "id": "8b7379bc", "metadata": {}, "source": [ "### Why the small dense iteration does not converge\n", "\n", "This is the four-asset, **four-period** seed-0 example used by `compare`,\n", "not the eight-asset, thirty-period allocation check above, nor the longer\n", "multi-seed study. Keep this short example as a failure case: its risk targets\n", "are upper bounds, not a promise that the optimiser will use that much risk.\n", "Costs and the individual forecasts leave the stand-alone solutions below\n", "those bounds even in the first fit. Two NNLS coefficients exceed one, so the\n", "undamped update raises already inactive bounds without enlarging the proxies.\n", "The other two coefficients are zero and their proxies drop out.\n", "\n", "The diagnostics below print every fit's coefficients, targets, maximum realised\n", "risks and binding fractions. `history` also holds each period's realised and\n", "raw risk, solver status and iteration count, primal equality and box errors,\n", "snapping changes and stationarity residuals. These residuals are numerical\n", "checks, not a claim that the solver's internal duality gap is available.\n", "The declared Clarabel tolerances remain unchanged.\n", "\n", "We check the cause by removing the risk constraint altogether, doubling the\n", "initial bounds, and repeating without snapping. We then compare a separate\n", "**optional** geometric damped update (power 0.5) with the original method.\n", "Finally, a proxy-stability diagnostic identifies redundant reruns; it is not\n", "the coefficient stopping rule and must not be called convergence.\n" ] }, { "cell_type": "code", "execution_count": 28, "id": "2e0d77f5", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:36.578619Z", "iopub.status.busy": "2026-10-29T10:22:36.578528Z", "iopub.status.idle": "2026-10-29T10:22:37.329552Z", "shell.execute_reply": "2026-10-29T10:22:37.329180Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 beta [ 3.3523 0. 10.440506 0. ] target [1. 1. 1. 1.] max risk [0.365357 0.002152 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "2 beta [ 3.3523 0. 10.440506 0. ] target [ 3.3523 0. 10.440506 1. ] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "3 beta [ 3.3523 0. 10.440506 0. ] target [ 11.237916 0. 109.004159 1. ] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "4 beta [ 3.3523 0. 10.440506 0. ] target [3.767287e+01 0.000000e+00 1.138059e+03 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "5 beta [ 3.3523 0. 10.440506 0. ] target [1.262907e+02 0.000000e+00 1.188191e+04 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "6 beta [ 3.3523 0. 10.440506 0. ] target [4.233645e+02 0.000000e+00 1.240531e+05 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "7 beta [ 3.3523 0. 10.440506 0. ] target [1.419245e+03 0.000000e+00 1.295177e+06 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "8 beta [ 3.3523 0. 10.440506 0. ] target [4.757734e+03 0.000000e+00 1.352231e+07 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "9 beta [ 3.3523 0. 10.440506 0. ] target [1.594935e+04 0.000000e+00 1.411797e+08 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "10 beta [ 3.3523 0. 10.440506 0. ] target [5.346702e+04 0.000000e+00 1.473987e+09 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "11 beta [ 3.3523 0. 10.440506 0. ] target [1.792375e+05 0.000000e+00 1.538917e+10 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "12 beta [ 3.3523 0. 10.440506 0. ] target [6.008578e+05 0.000000e+00 1.606708e+11 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "13 beta [ 3.3523 0. 10.440506 0. ] target [2.014256e+06 0.000000e+00 1.677484e+12 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "14 beta [ 3.3523 0. 10.440506 0. ] target [6.752390e+06 0.000000e+00 1.751378e+13 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "15 beta [ 3.3523 0. 10.440506 0. ] target [2.263604e+07 0.000000e+00 1.828527e+14 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "16 beta [ 3.3523 0. 10.440506 0. ] target [7.588279e+07 0.000000e+00 1.909075e+15 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "17 beta [ 3.3523 0. 10.440506 0. ] target [2.543819e+08 0.000000e+00 1.993171e+16 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "18 beta [ 3.3523 0. 10.440506 0. ] target [8.527643e+08 0.000000e+00 2.080971e+17 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "19 beta [ 3.3523 0. 10.440506 0. ] target [2.858722e+09 0.000000e+00 2.172639e+18 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "20 beta [ 3.3523 0. 10.440506 0. ] target [9.583294e+09 0.000000e+00 2.268345e+19 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "21 beta [ 3.3523 0. 10.440506 0. ] target [3.212608e+10 0.000000e+00 2.368267e+20 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "22 beta [ 3.3523 0. 10.440506 0. ] target [1.076962e+11 0.000000e+00 2.472590e+21 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "23 beta [ 3.3523 0. 10.440506 0. ] target [3.610301e+11 0.000000e+00 2.581509e+22 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "24 beta [ 3.3523 0. 10.440506 0. ] target [1.210281e+12 0.000000e+00 2.695226e+23 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "25 beta [ 3.3523 0. 10.440506 0. ] target [4.057226e+12 0.000000e+00 2.813953e+24 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "26 beta [ 3.3523 0. 10.440506 0. ] target [1.360104e+13 0.000000e+00 2.937909e+25 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "27 beta [ 3.3523 0. 10.440506 0. ] target [4.559477e+13 0.000000e+00 3.067325e+26 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "28 beta [ 3.3523 0. 10.440506 0. ] target [1.528473e+14 0.000000e+00 3.202443e+27 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "29 beta [ 3.3523 0. 10.440506 0. ] target [5.123901e+14 0.000000e+00 3.343512e+28 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "30 beta [ 3.3523 0. 10.440506 0. ] target [1.717685e+15 0.000000e+00 3.490796e+29 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "31 beta [ 3.3523 0. 10.440506 0. ] target [5.758197e+15 0.000000e+00 3.644567e+30 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "32 beta [ 3.3523 0. 10.440506 0. ] target [1.930320e+16 0.000000e+00 3.805113e+31 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "33 beta [ 3.3523 0. 10.440506 0. ] target [6.471013e+16 0.000000e+00 3.972730e+32 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "34 beta [ 3.3523 0. 10.440506 0. ] target [2.169278e+17 0.000000e+00 4.147731e+33 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "35 beta [ 3.3523 0. 10.440506 0. ] target [7.272070e+17 0.000000e+00 4.330441e+34 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "36 beta [ 3.3523 0. 10.440506 0. ] target [2.437816e+18 0.000000e+00 4.521199e+35 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "37 beta [ 3.3523 0. 10.440506 0. ] target [8.172291e+18 0.000000e+00 4.720361e+36 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "38 beta [ 3.3523 0. 10.440506 0. ] target [2.739597e+19 0.000000e+00 4.928295e+37 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "39 beta [ 3.3523 0. 10.440506 0. ] target [9.183952e+19 0.000000e+00 5.145389e+38 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "40 beta [ 3.3523 0. 10.440506 0. ] target [3.078736e+20 0.000000e+00 5.372047e+39 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "41 beta [ 3.3523 0. 10.440506 0. ] target [1.032085e+21 0.000000e+00 5.608688e+40 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "42 beta [ 3.3523 0. 10.440506 0. ] target [3.459858e+21 0.000000e+00 5.855754e+41 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "43 beta [ 3.3523 0. 10.440506 0. ] target [1.159848e+22 0.000000e+00 6.113703e+42 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "44 beta [ 3.3523 0. 10.440506 0. ] target [3.888159e+22 0.000000e+00 6.383015e+43 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "45 beta [ 3.3523 0. 10.440506 0. ] target [1.303428e+23 0.000000e+00 6.664191e+44 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "46 beta [ 3.3523 0. 10.440506 0. ] target [4.369480e+23 0.000000e+00 6.957752e+45 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "47 beta [ 3.3523 0. 10.440506 0. ] target [1.464781e+24 0.000000e+00 7.264245e+46 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "48 beta [ 3.3523 0. 10.440506 0. ] target [4.910385e+24 0.000000e+00 7.584239e+47 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "49 beta [ 3.3523 0. 10.440506 0. ] target [1.646108e+25 0.000000e+00 7.918329e+48 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "50 beta [ 3.3523 0. 10.440506 0. ] target [5.518249e+25 0.000000e+00 8.267136e+49 1.000000e+00] max risk [0.365357 0. 0.162856 0. ] binding fraction [0. 0. 0. 0.]\n", "Risk-unconstrained vs initial proxy P&L max error: 3.0635910858427956e-07\n", "Doubled targets vs initial P&L max error: 1.1062087527241582e-08\n", "Unsnapped full NNLS beta (includes numerical near-zero columns): [2.832514e+00 1.156902e+02 2.976582e+00 1.753866e+08]\n", "Unsnapped proxy P&L column norms: [2.886402e-01 6.154359e-03 1.788437e-01 4.912001e-09]\n", "Unsnapped NNLS restricted to the same surviving columns: [ 3.3523 10.440506]\n", "Risk-unconstrained beta: [ 3.3523 0. 10.440506 0. ]\n", "Unsnapped full-fit residual norm: 5.283137298900343e-16\n", "Snap tolerance 1e-08 full NNLS beta: [ 3.3523 0. 10.440506 0. ]\n", "Snap tolerance 1e-05 full NNLS beta: [ 3.3523 0. 10.440506 0. ]\n", "Snap tolerance 0.0001 full NNLS beta: [ 3.3523 0. 10.440506 0. ]\n", "Maximum equality_error 5.551115123125783e-17\n", "Maximum box_violation 0.0\n", "Maximum stationarity_error 1.0655092225198004e-09\n", "Maximum snap_change 9.609189828524719e-09\n", "Solver statuses: ['optimal']\n", "Maximum raw risk excess: -0.634643185121041\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Optional power-0.5 damping: max_iter fits: 50 beta: [ 3.3523 0. 10.440506 0. ] targets: [7.428492e+12 0.000000e+00 9.092379e+24 1.000000e+00]\n", "Optional absolute proxy-P&L stability diagnostic (1e-6): fit 3 — stalled proxies, NOT coefficient convergence\n" ] } ], "source": [ "np.set_printoptions(precision=6)\n", "for row in iteration_dense[\"history\"]:\n", " print(row[\"iteration\"], \"beta\", row[\"beta\"], \"target\", row[\"risk_targets\"],\n", " \"max risk\", row[\"realised_risk\"].max(axis=1),\n", " \"binding fraction\", row[\"binding_fraction\"])\n", "\n", "unconstrained = solve_cvxpy(small_dense, g_small,\n", " gamma_path=np.zeros((len(g_small), small_dense.n_periods)))\n", "doubled = solve_cvxpy(small_dense, g_small, sigma_target=np.full(len(g_small), 2.0))\n", "unsnapped = solve_cvxpy(small_dense, g_small,\n", " sigma_target=np.ones(len(g_small)), snap_tolerance=0)\n", "first = iteration_dense[\"history\"][0]\n", "free_pnl = pnl(small_dense, unconstrained, r_small)\n", "free_beta = regress(total_dense, free_pnl)[0]\n", "print(\"Risk-unconstrained vs initial proxy P&L max error:\", np.max(np.abs(free_pnl-first[\"proxies\"])))\n", "print(\"Doubled targets vs initial P&L max error:\",\n", " np.max(np.abs(pnl(small_dense, doubled, r_small)-first[\"proxies\"])))\n", "unsnapped_pnl = pnl(small_dense, unsnapped, r_small)\n", "print(\"Unsnapped full NNLS beta (includes numerical near-zero columns):\",\n", " regress(total_dense, unsnapped_pnl)[0])\n", "print(\"Unsnapped proxy P&L column norms:\", np.linalg.norm(unsnapped_pnl, axis=1))\n", "live = np.flatnonzero(first[\"beta\"] > 0)\n", "print(\"Unsnapped NNLS restricted to the same surviving columns:\",\n", " regress(total_dense, unsnapped_pnl[live])[0])\n", "print(\"Risk-unconstrained beta:\", free_beta)\n", "print(\"Unsnapped full-fit residual norm:\", np.linalg.norm(\n", " regress(total_dense, unsnapped_pnl)[2]))\n", "# Vary the numerical zero convention without preselecting the surviving columns.\n", "for snap_tol in (1e-8, 1e-5, 1e-4):\n", " sensitivity = solve_cvxpy(small_dense, g_small,\n", " sigma_target=np.ones(len(g_small)), snap_tolerance=snap_tol)\n", " sensitivity_beta = regress(total_dense, pnl(small_dense, sensitivity, r_small))[0]\n", " print(\"Snap tolerance\", snap_tol, \"full NNLS beta:\", sensitivity_beta)\n", " np.testing.assert_allclose(sensitivity_beta, first[\"beta\"], atol=1e-5, rtol=0)\n", "\n", "np.testing.assert_allclose(free_pnl, first[\"proxies\"], atol=1e-6, rtol=0)\n", "np.testing.assert_allclose(pnl(small_dense, doubled, r_small), first[\"proxies\"], atol=1e-6, rtol=0)\n", "np.testing.assert_allclose(unsnapped_pnl, first[\"proxies\"], atol=1e-6, rtol=0)\n", "np.testing.assert_allclose(regress(total_dense, unsnapped_pnl[live])[0],\n", " first[\"beta\"][live], atol=1e-5, rtol=0)\n", "assert free_beta[0] > 3 and free_beta[2] > 10\n", "assert all(np.all(row[\"binding_fraction\"] == 0) for row in iteration_dense[\"history\"])\n", "for key in (\"equality_error\", \"box_violation\", \"stationarity_error\", \"snap_change\"):\n", " print(\"Maximum\", key, max(np.max(row[\"solver\"][key]) for row in iteration_dense[\"history\"]))\n", "print(\"Solver statuses:\", sorted({s for row in iteration_dense[\"history\"]\n", " for s in row[\"solver\"][\"status\"].ravel()}))\n", "print(\"Maximum raw risk excess:\", max(np.max(row[\"solver\"][\"raw_risk\"] -\n", " row[\"risk_targets\"][row[\"solved\"], None]) for row in iteration_dense[\"history\"]))\n", "\n", "_, damped_dense = iterated_proxies(small_dense, g_small, total_dense, r_small,\n", " trace=True, update_power=0.5)\n", "print(\"Optional power-0.5 damping:\", damped_dense[\"status\"],\n", " \"fits:\", damped_dense[\"n_iter\"], \"beta:\", damped_dense[\"beta\"],\n", " \"targets:\", damped_dense[\"risk_targets\"])\n", "assert damped_dense[\"status\"] == iteration_dense[\"status\"] == \"max_iter\"\n", "# Compare only after dropping zero-coefficient proxies; otherwise a dropped\n", "# nonzero P&L column would itself look like movement.\n", "stationary_fit = None\n", "for previous, current in zip(iteration_dense[\"history\"], iteration_dense[\"history\"][1:]):\n", " if (np.array_equal(previous[\"solved\"], current[\"solved\"])\n", " and np.max(np.abs(current[\"proxies\"] - previous[\"proxies\"])) < 1e-6):\n", " stationary_fit = current[\"iteration\"]\n", " break\n", "print(\"Optional absolute proxy-P&L stability diagnostic (1e-6): fit\", stationary_fit,\n", " \"— stalled proxies, NOT coefficient convergence\")\n", "assert stationary_fit is not None\n" ] }, { "cell_type": "markdown", "id": "9ea7a0ad", "metadata": {}, "source": [ "The printed coefficients stay near 3.3523 and 10.4405 on the two surviving\n", "proxies. Their maximum realised risks stay near 0.36536 and 0.16286; no period's\n", "risk constraint binds. The risk-unconstrained counterfactual reproduces the same result. Removing snapping\n", "preserves proxy P&Ls within $10^{-6}$ but can change the full NNLS fit:\n", "a numerical near-zero column receives an enormous coefficient. Restricting\n", "that sensitivity check to the same surviving columns recovers their coefficients.\n", "The full unsnapped fit also moves the third coefficient from about 10.44 to\n", "2.98 and has a near-zero residual with four columns for four observations.\n", "Coefficient magnitudes therefore depend on the numerical zero convention and\n", "resulting active set. The separate full-fit checks vary the snapping tolerance\n", "without fixing that active set in advance. Thus NNLS is sensitive to numerical zero columns; the existing no-trade snapping\n", "convention matters and is retained. This is target saturation with growth\n", "in the targets, not a period-two coefficient cycle. Numerical warnings are\n", "retained, but numerical error does not explain coefficient errors of this size.\n", "\n", "Halving the logarithmic update only slows the target growth; it does not solve\n", "this case, so this case does not motivate damping as a remedy. A proxy-stability\n", "stop would save redundant work and report **stalled**, not converged. We leave\n", "the default method and its `max_iter` result intact. This trajectory does not\n", "prove that no fixed point exists at other initial targets. The study's 80%\n", "convergence across twenty longer one-factor runs is a different experiment.\n" ] }, { "cell_type": "markdown", "id": "signal-attribution-54", "metadata": { "part": "reproduction" }, "source": [ "## Check the two published paths\n", "\n", "The reference below contains the final cumulative values copied from the figure's `f13-pnl-paths.json`, including residuals and the two Grinold–Easton portfolios. It is embedded so this notebook works offline. The source JSON SHA-256 is `91b750176b9263260aa13a85063dbf2c1f3b12bb561139e752bf4f9cc8f53544`.\n", "\n", "The check uses absolute tolerance $10^{-8}$ simulation P&L units and zero relative tolerance. It compares every component endpoint for every method and the book endpoint. It also checks that each method, including any residual, adds back to the book in every period. Neither check changes the reference. The examples may take a few minutes." ] }, { "cell_type": "code", "execution_count": 29, "id": "signal-attribution-55", "metadata": { "execution": { "iopub.execute_input": "2026-10-29T10:22:37.330910Z", "iopub.status.busy": "2026-10-29T10:22:37.330833Z", "iopub.status.idle": "2026-10-29T10:22:43.683104Z", "shell.execute_reply": "2026-10-29T10:22:43.682579Z" }, "part": "reproduction" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "base: book=519.091093822072; iteration=converged (10 fits)\n", " regression: fast A=100.532420647, slow A=70.389548587, fast B=14.304523929, slow B=5.152252732, residual=328.712347927\n", " emmanuel: fast A=152.764529211, slow A=112.631195043, fast B=58.463719033, slow B=9.128433519, residual=186.103217017\n", " loo: fast A=111.606971613, slow A=119.726816978, fast B=45.547217806, slow B=46.361535410, residual=195.848552015\n", " grinold_easton: fast A=256.728484749, slow A=366.925527849, fast B=130.117258842, slow B=152.894025039, spread cost=-241.893975086, position limit=-145.680227571\n", " durin: fast A=175.062474127, slow A=179.264382658, fast B=101.859995050, slow B=62.904241988\n", " shapley_target: fast A=177.489940858, slow A=176.863660931, fast B=98.660627553, slow B=66.076864480\n", " shapley_gamma: fast A=165.110443254, slow A=201.341402322, fast B=84.365643592, slow B=68.273604654\n", "PASS: all component endpoints and period sums; max endpoint error=0\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "spread_2: book=104.328336610861; iteration=no_live_proxy (2 fits)\n", " regression: fast A=0.000000000, slow A=0.000000000, fast B=0.000000000, slow B=0.000000000, residual=104.328336611\n", " emmanuel: fast A=0.000000000, slow A=0.000000000, fast B=0.000000000, slow B=0.000000000, residual=104.328336611\n", " loo: fast A=76.941784382, slow A=79.718501859, fast B=83.241510123, slow B=38.198069495, residual=-173.771529248\n", " grinold_easton: fast A=584.582758852, slow A=1411.978222450, fast B=427.665383299, slow B=414.989056887, spread cost=-2717.716466795, position limit=-17.170618083\n", " durin: fast A=39.747235338, slow A=25.440018176, fast B=27.503588409, slow B=11.637494689\n", " shapley_target: fast A=30.384940844, slow A=30.230508341, fast B=29.720064717, slow B=13.992822709\n", " shapley_gamma: fast A=31.427851827, slow A=29.636409594, fast B=30.215339873, slow B=13.048735316\n", "PASS: all component endpoints and period sums; max endpoint error=0\n" ] } ], "source": [ "reference = {'base': {'total': 519.0910938220723,\n", " 'methods': {'regression': {'fast A': 100.53242064658039,\n", " 'slow A': 70.38954858667016,\n", " 'fast B': 14.304523929024574,\n", " 'slow B': 5.152252732430301,\n", " 'residual': 328.7123479273673},\n", " 'emmanuel': {'fast A': 152.764529210917,\n", " 'slow A': 112.63119504250166,\n", " 'fast B': 58.46371903282879,\n", " 'slow B': 9.128433518895262,\n", " 'residual': 186.10321701692942},\n", " 'loo': {'fast A': 111.60697161328982,\n", " 'slow A': 119.72681697759026,\n", " 'fast B': 45.54721780588971,\n", " 'slow B': 46.361535410303155,\n", " 'residual': 195.84855201499892},\n", " 'grinold_easton': {'fast A': 256.72848474852657,\n", " 'slow A': 366.92552784896986,\n", " 'fast B': 130.11725884186922,\n", " 'slow B': 152.89402503899424,\n", " 'spread cost': -241.89397508564036,\n", " 'position limit': -145.68022757064776},\n", " 'durin': {'fast A': 175.06247412657459,\n", " 'slow A': 179.2643826576884,\n", " 'fast B': 101.85999505010044,\n", " 'slow B': 62.904241987708474},\n", " 'shapley_target': {'fast A': 177.48994085845976,\n", " 'slow A': 176.86366093064652,\n", " 'fast B': 98.66062755250604,\n", " 'slow B': 66.07686448045925},\n", " 'shapley_gamma': {'fast A': 165.110443254377,\n", " 'slow A': 201.34140232216058,\n", " 'fast B': 84.36564359163518,\n", " 'slow B': 68.27360465389944}}},\n", " 'spread_2': {'total': 104.32833661086099,\n", " 'methods': {'regression': {'fast A': 0.0,\n", " 'slow A': 0.0,\n", " 'fast B': 0.0,\n", " 'slow B': 0.0,\n", " 'residual': 104.32833661086099},\n", " 'emmanuel': {'fast A': 0.0,\n", " 'slow A': 0.0,\n", " 'fast B': 0.0,\n", " 'slow B': 0.0,\n", " 'residual': 104.32833661086099},\n", " 'loo': {'fast A': 76.94178438231884,\n", " 'slow A': 79.71850185927482,\n", " 'fast B': 83.2415101226865,\n", " 'slow B': 38.198069494507166,\n", " 'residual': -173.77152924792617},\n", " 'grinold_easton': {'fast A': 584.582758851705,\n", " 'slow A': 1411.9782224503138,\n", " 'fast B': 427.665383299356,\n", " 'slow B': 414.98905688713205,\n", " 'spread cost': -2717.716466794506,\n", " 'position limit': -17.1706180831368},\n", " 'durin': {'fast A': 39.74723533758032,\n", " 'slow A': 25.44001817571785,\n", " 'fast B': 27.503588408741983,\n", " 'slow B': 11.637494688821034},\n", " 'shapley_target': {'fast A': 30.38494084386784,\n", " 'slow A': 30.230508340967734,\n", " 'fast B': 29.720064717153566,\n", " 'slow B': 13.992822708871934},\n", " 'shapley_gamma': {'fast A': 31.427851827138475,\n", " 'slow A': 29.636409594471246,\n", " 'fast B': 30.21533987280228,\n", " 'slow B': 13.048735316449124}}}}\n", "\n", "tolerance = 1e-8\n", "for scenario, spread in [(\"base\", 0.5), (\"spread_2\", 2.0)]:\n", " desk = Desk(spread_rel=spread)\n", " total, parts, iteration = compare(desk)\n", " expected = reference[scenario]\n", " book_endpoint = np.cumsum(total)[-1]\n", " np.testing.assert_allclose(book_endpoint, expected[\"total\"], atol=tolerance, rtol=0)\n", " print(f\"{scenario}: book={book_endpoint:.12f}; \"\n", " f\"iteration={iteration['status']} ({iteration['n_iter']} fits)\")\n", " largest_error = abs(book_endpoint - expected[\"total\"])\n", " for method, rows in expected[\"methods\"].items():\n", " endpoints = np.cumsum(parts[method], axis=1)[:, -1]\n", " np.testing.assert_allclose(endpoints, list(rows.values()), atol=tolerance, rtol=0)\n", " np.testing.assert_allclose(parts[method].sum(0), total, atol=tolerance, rtol=0)\n", " largest_error = max(largest_error, np.max(np.abs(endpoints - list(rows.values()))))\n", " print(f\" {method}: \" + \", \".join(\n", " f\"{label}={value:.9f}\" for label, value in zip(rows, endpoints)))\n", " print(f\"PASS: all component endpoints and period sums; max endpoint error={largest_error:.3g}\")" ] } ], "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.13.0" } }, "nbformat": 4, "nbformat_minor": 5 }