criptic1/ctrader-ml-systems

1"""2Stieltjes propagator cross-impact kernel calibration.3 4Implements the passivity-preserving reduced cross-impact model:5    G_r(t) = Σ_{k=1}^{r} A_k * exp(-ρ_k * t)6 7where each A_k ≽ 0 (PSD) and ρ_k > 0.8 9This ensures:101. Admissibility (no-manipulation) by construction — the kernel is a11   Stieltjes propagator, hence passive.122. Finite-dimensional state-space realization for fast execution optimization.133. Constructive error bounds: ||G - G_r||_{L^1} controls execution cost error.14 15Calibration pipeline:16    1. Empirical kernel estimation from trades/quotes data17    2. Scalar NNLS warmstart for decay rate initialization18    3. Matrix SDP fitting with PSD constraints via cvxpy19    4. Alternating optimization (poles + residues)20    5. Passivity verification test suite21 22References:23- Core theory: Passivity-preserving cross-impact reduction (uploaded study)24- Calibration: Inverse Laplace with positivity constraints25- SDP: cvxpy with SCS/MOSEK solvers26"""27 28import logging29from typing import Optional, List, Tuple, Dict30import numpy as np31from scipy.optimize import nnls, minimize32from scipy.linalg import eigh33 34from ..utils.types import CrossImpactKernel35 36logger = logging.getLogger(__name__)37 38 39class StieltjesKernelCalibrator:40    """41    Calibrate a reduced Stieltjes cross-impact kernel from market data.42    43    The model:44        G_r(t) = Σ_{k=1}^{K} A_k * exp(-ρ_k * t)45    46    where:47        - ρ_k > 0 are decay rates (poles)48        - A_k ≽ 0 are PSD residue matrices (d × d)49        - K is the model order (number of exponential modes)50    51    Admissibility is guaranteed by construction: any non-negative combination52    of decaying exponentials with PSD coefficients is a Stieltjes propagator.53    """54 55    def __init__(56        self,57        n_modes: int = 8,58        n_assets: int = 5,59        rho_range: Tuple[float, float] = (0.01, 100.0),60        regularization: float = 1e-4,61        solver: str = "SCS",62        max_alternating_iter: int = 20,63    ):64        """65        Args:66            n_modes: K — number of exponential decay modes67            n_assets: d — number of assets68            rho_range: (min, max) for decay rates69            regularization: nuclear norm penalty for low-rank bias70            solver: cvxpy solver ("SCS" or "MOSEK")71            max_alternating_iter: iterations for alternating optimization72        """73        self.n_modes = n_modes74        self.n_assets = n_assets75        self.rho_range = rho_range76        self.regularization = regularization77        self.solver = solver78        self.max_alternating_iter = max_alternating_iter79 80    def estimate_empirical_kernel(81        self,82        signed_flow: np.ndarray,83        price_changes: np.ndarray,84        n_lags: int = 50,85    ) -> Tuple[np.ndarray, List[np.ndarray]]:86        """87        Estimate empirical cross-impact kernel G(t) from trade data.88        89        Uses regression of price changes on lagged signed order flow:90            Δp_t = Σ_{s=0}^{L} G(s) * v_{t-s} + ε_t91        92        Args:93            signed_flow: (T, d) matrix of signed trade volume per asset94            price_changes: (T, d) matrix of price changes per asset95            n_lags: number of lag points to estimate96            97        Returns:98            t_obs: array of lag times99            G_obs: list of estimated d×d kernel matrices at each lag100        """101        T, d = signed_flow.shape102        assert price_changes.shape == (T, d), "Shape mismatch"103 104        t_obs = np.arange(n_lags, dtype=float)105        G_obs = []106 107        for lag in range(n_lags):108            # For each lag, regress price_changes on lagged flow109            if lag > 0:110                y = price_changes[lag:]        # (T-lag, d)111                X = signed_flow[:-lag]          # (T-lag, d)112            else:113                y = price_changes114                X = signed_flow115 116            # Cross-impact at this lag: G(lag) ≈ (X'X)^{-1} X' y117            # For d-dimensional: G is d×d118            XtX = X.T @ X119            XtY = X.T @ y120 121            try:122                G_lag = np.linalg.solve(XtX + 1e-8 * np.eye(d), XtY)123            except np.linalg.LinAlgError:124                G_lag = np.linalg.lstsq(X, y, rcond=None)[0]125 126            # Symmetrize (cross-impact should be symmetric for passivity)127            G_lag = 0.5 * (G_lag + G_lag.T)128            G_obs.append(G_lag)129 130        return t_obs, G_obs131 132    def fit_scalar_warmstart(133        self, t_obs: np.ndarray, G_obs: List[np.ndarray]134    ) -> Tuple[np.ndarray, np.ndarray]:135        """136        Step 1: Scalar NNLS warmstart.137        138        Fit trace(G(t)) = Σ a_k * exp(-ρ_k * t) with a_k ≥ 0.139        This gives good initial decay rates ρ_k.140        141        Returns:142            rho: initial decay rates, shape (K,)143            a_scalar: scalar weights, shape (K,)144        """145        # Use trace of G as scalar proxy146        g_scalar = np.array([np.trace(G) for G in G_obs])147 148        # Dense pole grid for NNLS149        rho_grid = np.logspace(150            np.log10(self.rho_range[0]),151            np.log10(self.rho_range[1]),152            100,153        )154 155        # Build basis: Φ[i,k] = exp(-ρ_k * t_i)156        Phi = np.exp(-np.outer(t_obs, rho_grid))157 158        # Non-negative least squares: a ≥ 0159        a_hat, residual = nnls(Phi, np.maximum(g_scalar, 0))160 161        # Select top-K modes by weight162        top_idx = np.argsort(a_hat)[-self.n_modes:]163        rho_init = rho_grid[top_idx]164        a_init = a_hat[top_idx]165 166        # Sort by decay rate167        sort_idx = np.argsort(rho_init)168        return rho_init[sort_idx], a_init[sort_idx]169 170    def fit_matrix_sdp(171        self,172        t_obs: np.ndarray,173        G_obs: List[np.ndarray],174        rho_fixed: np.ndarray,175    ) -> List[np.ndarray]:176        """177        Step 2: SDP-constrained matrix fitting.178        179        Fix decay rates ρ_k, solve for PSD residue matrices A_k:180            min Σ_t ||G_obs(t) - Σ_k A_k exp(-ρ_k t)||_F^2 + λ Σ_k ||A_k||_*181            s.t. A_k ≽ 0  ∀k182        183        Args:184            t_obs: time points185            G_obs: observed kernel matrices186            rho_fixed: fixed decay rates187            188        Returns:189            List of K PSD matrices A_k190        """191        try:192            import cvxpy as cp193        except ImportError:194            logger.error("cvxpy not installed. pip install cvxpy[scs]")195            raise196 197        K = len(rho_fixed)198        d = self.n_assets199 200        # Decision variables201        A_vars = [cp.Variable((d, d), symmetric=True) for _ in range(K)]202 203        # PSD constraints — THE passivity guarantee204        constraints = [A >> 0 for A in A_vars]205 206        # Reconstruction loss207        loss = 0208        for i, (t, G_true) in enumerate(zip(t_obs, G_obs)):209            G_pred = sum(210                A_vars[k] * float(np.exp(-rho_fixed[k] * t))211                for k in range(K)212            )213            loss += cp.sum_squares(G_pred - G_true)214 215        # Nuclear norm regularization (promotes low rank → fewer effective modes)216        reg = self.regularization * sum(cp.normNuc(A) for A in A_vars)217 218        objective = cp.Minimize(loss + reg)219        prob = cp.Problem(objective, constraints)220 221        # Solve222        solver_map = {"SCS": cp.SCS, "MOSEK": cp.MOSEK, "ECOS": cp.ECOS}223        solver = solver_map.get(self.solver, cp.SCS)224        225        try:226            prob.solve(solver=solver, verbose=False, eps=1e-4)227        except cp.SolverError:228            logger.warning(f"Solver {self.solver} failed, trying SCS fallback")229            prob.solve(solver=cp.SCS, verbose=False, eps=1e-3)230 231        if prob.status not in ("optimal", "optimal_inaccurate"):232            logger.warning(f"SDP solve status: {prob.status}")233 234        # Extract solutions (clip negative eigenvalues for numerical safety)235        A_fitted = []236        for A in A_vars:237            A_val = A.value238            if A_val is None:239                A_val = np.eye(d) * 1e-6240            # Project to PSD (numerical safety)241            A_val = self._project_psd(A_val)242            A_fitted.append(A_val)243 244        return A_fitted245 246    def fit_alternating(247        self,248        t_obs: np.ndarray,249        G_obs: List[np.ndarray],250    ) -> CrossImpactKernel:251        """252        Full calibration: alternating optimization of poles and residues.253        254        Phase 1: NNLS warmstart → initial ρ255        Phase 2: For each iteration:256            (a) Fix ρ → SDP for A_k (convex)257            (b) Fix A_k → L-BFGS-B for ρ (non-convex but low-dimensional)258        259        Returns:260            CrossImpactKernel with calibrated parameters261        """262        logger.info(f"Starting Stieltjes kernel calibration: "263                    f"K={self.n_modes}, d={self.n_assets}")264 265        # Phase 1: Warmstart266        rho, a_scalar = self.fit_scalar_warmstart(t_obs, G_obs)267        logger.info(f"Warmstart: ρ = {rho}")268 269        best_loss = float("inf")270        best_rho = rho.copy()271        best_A = None272 273        for it in range(self.max_alternating_iter):274            # Phase 2a: Fix ρ, solve SDP for A275            try:276                A_fitted = self.fit_matrix_sdp(t_obs, G_obs, rho)277            except Exception as e:278                logger.warning(f"SDP iteration {it} failed: {e}")279                break280 281            # Compute current loss282            loss = self._compute_loss(t_obs, G_obs, rho, A_fitted)283            logger.debug(f"Iteration {it}: loss = {loss:.6f}")284 285            if loss < best_loss:286                best_loss = loss287                best_rho = rho.copy()288                best_A = [A.copy() for A in A_fitted]289 290            # Phase 2b: Fix A, optimize ρ via L-BFGS-B291            rho = self._optimize_poles(t_obs, G_obs, A_fitted, rho)292 293        if best_A is None:294            raise RuntimeError("Calibration failed: no valid solution found")295 296        logger.info(f"Calibration complete: loss = {best_loss:.6f}")297 298        return CrossImpactKernel(299            rho=best_rho,300            A_residues=best_A,301            n_assets=self.n_assets,302            n_modes=self.n_modes,303        )304 305    def _optimize_poles(306        self,307        t_obs: np.ndarray,308        G_obs: List[np.ndarray],309        A_fixed: List[np.ndarray],310        rho_init: np.ndarray,311    ) -> np.ndarray:312        """Optimize decay rates with fixed residues via L-BFGS-B."""313 314        def objective(log_rho):315            rho = np.exp(log_rho)316            return self._compute_loss(t_obs, G_obs, rho, A_fixed)317 318        log_rho_init = np.log(rho_init)319        bounds = [320            (np.log(self.rho_range[0]), np.log(self.rho_range[1]))321        ] * self.n_modes322 323        result = minimize(324            objective, log_rho_init,325            method="L-BFGS-B", bounds=bounds,326            options={"maxiter": 50, "ftol": 1e-8},327        )328 329        return np.exp(result.x)330 331    def _compute_loss(332        self,333        t_obs: np.ndarray,334        G_obs: List[np.ndarray],335        rho: np.ndarray,336        A_list: List[np.ndarray],337    ) -> float:338        """Frobenius reconstruction loss."""339        loss = 0.0340        for t, G_true in zip(t_obs, G_obs):341            G_pred = sum(A * np.exp(-rho[k] * t) for k, A in enumerate(A_list))342            loss += np.sum((G_pred - G_true) ** 2)343        return loss344 345    @staticmethod346    def _project_psd(M: np.ndarray, tol: float = 0.0) -> np.ndarray:347        """Project symmetric matrix to nearest PSD via eigenvalue clipping."""348        eigvals, eigvecs = eigh(M)349        eigvals_clipped = np.maximum(eigvals, tol)350        return eigvecs @ np.diag(eigvals_clipped) @ eigvecs.T351 352 353class PassivityTestSuite:354    """355    Comprehensive passivity verification for cross-impact kernels.356    357    Tests that the kernel G satisfies the no-manipulation condition:358        ⟨u, K_G u⟩ ≥ 0  for all admissible trading programs u359    360    This is equivalent to passivity of the causal operator K_G.361    362    Tests:363    1. PSD residue check: eigenvalues of each A_k ≥ 0364    2. Kernel PSD check: eigenvalues of G(t) ≥ 0 at sample times365    3. Round-trip cost check: cost of round-trip programs ≥ 0366    4. Energy dissipation check: storage function is non-increasing367    5. L^1 approximation error bound tracking368    """369 370    def __init__(self, n_random_tests: int = 1000, tolerance: float = 1e-8):371        self.n_random_tests = n_random_tests372        self.tolerance = tolerance373 374    def run_all_tests(375        self,376        kernel: CrossImpactKernel,377        T: float = 10.0,378        dt: float = 0.01,379        reference_kernel: Optional[CrossImpactKernel] = None,380    ) -> Dict[str, dict]:381        """382        Run the complete passivity test suite.383        384        Args:385            kernel: the reduced kernel to test386            T: time horizon for tests387            dt: time discretization388            reference_kernel: full kernel for approximation error bounds389            390        Returns:391            Dict of test name -> {passed: bool, value: float, details: ...}392        """393        results = {}394 395        # Test 1: PSD residue check396        results["psd_residues"] = self._test_psd_residues(kernel)397 398        # Test 2: Kernel PSD at sample times399        results["kernel_psd"] = self._test_kernel_psd(kernel, T)400 401        # Test 3: Round-trip cost402        results["round_trip_cost"] = self._test_round_trip_cost(kernel, T, dt)403 404        # Test 4: Energy dissipation405        results["energy_dissipation"] = self._test_energy_dissipation(kernel, T, dt)406 407        # Test 5: L^1 error bound (if reference provided)408        if reference_kernel is not None:409            results["l1_error_bound"] = self._test_l1_error(410                kernel, reference_kernel, T411            )412 413        # Overall pass/fail414        all_passed = all(r["passed"] for r in results.values())415        results["overall"] = {"passed": all_passed}416 417        return results418 419    def _test_psd_residues(self, kernel: CrossImpactKernel) -> dict:420        """Test 1: Each A_k must have non-negative eigenvalues."""421        min_eigenvalues = []422        for k, A in enumerate(kernel.A_residues):423            eigvals = np.linalg.eigvalsh(A)424            min_ev = float(np.min(eigvals))425            min_eigenvalues.append(min_ev)426 427        all_psd = all(ev >= -self.tolerance for ev in min_eigenvalues)428        return {429            "passed": all_psd,430            "min_eigenvalues_per_mode": min_eigenvalues,431            "worst_eigenvalue": min(min_eigenvalues),432        }433 434    def _test_kernel_psd(435        self, kernel: CrossImpactKernel, T: float, n_points: int = 100436    ) -> dict:437        """Test 2: G(t) must be PSD at all time points."""438        t_grid = np.linspace(0, T, n_points)439        G_vals = kernel.evaluate_grid(t_grid)440 441        min_eigenvalues = []442        for i in range(n_points):443            eigvals = np.linalg.eigvalsh(G_vals[i])444            min_eigenvalues.append(float(np.min(eigvals)))445 446        worst = min(min_eigenvalues)447        return {448            "passed": worst >= -self.tolerance,449            "worst_eigenvalue": worst,450            "n_violations": sum(1 for ev in min_eigenvalues if ev < -self.tolerance),451        }452 453    def _test_round_trip_cost(454        self, kernel: CrossImpactKernel, T: float, dt: float455    ) -> dict:456        """457        Test 3: Cost of round-trip trading programs must be ≥ 0.458        459        A round-trip program satisfies: ∫_0^T u(t) dt = 0460        (net zero position change — bought and sold same amount).461        462        The no-manipulation condition requires:463            C(u) = ⟨u, K_G u⟩ ≥ 0  for all such u464        """465        n_steps = int(T / dt)466        d = kernel.n_assets467        t_grid = np.linspace(0, T, n_steps)468 469        violations = 0470        worst_cost = float("inf")471 472        for trial in range(self.n_random_tests):473            # Generate random round-trip program474            u = self._generate_round_trip(n_steps, d)475 476            # Compute impact cost via convolution477            cost = self._compute_impact_cost(kernel, u, t_grid, dt)478 479            if cost < -self.tolerance:480                violations += 1481            worst_cost = min(worst_cost, cost)482 483        return {484            "passed": violations == 0,485            "n_violations": violations,486            "worst_cost": float(worst_cost),487            "n_tests": self.n_random_tests,488            "violation_rate": violations / self.n_random_tests,489        }490 491    def _test_energy_dissipation(492        self, kernel: CrossImpactKernel, T: float, dt: float493    ) -> dict:494        """495        Test 4: Energy/storage function check.496        497        The state-space realization has storage function:498            S(t) = (1/2) Σ_k z_k(t)' A_k z_k(t)499        500        The dissipation identity ensures S(t) ≥ 0 and tracks total energy.501        Verify that the energy balance:502            ⟨u, y⟩_{[0,T]} = S(T) + dissipated_energy ≥ 0503        504        Uses exact exponential integration for stability:505            z_k(t+dt) = exp(-ρ_k dt) z_k(t) + (1-exp(-ρ_k dt))/ρ_k * u_t506        """507        # Use finer dt for numerical stability (max dt*rho < 0.1)508        max_rho = max(kernel.rho)509        stable_dt = min(dt, 0.05 / max_rho)510        n_steps = max(int(T / stable_dt), 200)511        d = kernel.n_assets512 513        np.random.seed(12345)514        # Random trading program515        u = np.random.randn(n_steps, d) * 0.1516 517        # Simulate state-space realization with EXACT exponential integrator518        z = np.zeros((kernel.n_modes, d))  # state per mode519        total_input_energy = 0.0520 521        for t_idx in range(n_steps):522            u_t = u[t_idx]523 524            # Save pre-update state525            z_old = [z[k].copy() for k in range(kernel.n_modes)]526 527            # Exact exponential state update for stability:528            # z_k(t+dt) = exp(-ρ_k dt) * z_k(t) + (1-exp(-ρ_k dt))/ρ_k * u_t529            for k in range(kernel.n_modes):530                decay = np.exp(-kernel.rho[k] * stable_dt)531                if kernel.rho[k] > 1e-10:532                    input_gain = (1 - decay) / kernel.rho[k]533                else:534                    input_gain = stable_dt535                z[k] = decay * z_old[k] + input_gain * u_t536 537            # Symmetric output: average of pre- and post-update (midpoint rule)538            # This eliminates the O(dt) energy bias from asymmetric evaluation539            y_t = sum(540                kernel.A_residues[k] @ (0.5 * (z_old[k] + z[k]))541                for k in range(kernel.n_modes)542            )543 544            # Input-output energy: u' y dt545            total_input_energy += np.dot(u_t, y_t) * stable_dt546 547        # Storage at final time548        S_T = 0.5 * sum(549            z[k] @ kernel.A_residues[k] @ z[k]550            for k in range(kernel.n_modes)551        )552 553        # Dissipation identity: ⟨u,y⟩ = S(T) + dissipated ≥ 0554        # Allow tolerance proportional to integration scale (discrete approx noise)555        energy_scale = max(abs(total_input_energy), S_T, 1e-6)556        relative_violation = -total_input_energy / energy_scale if total_input_energy < 0 else 0557        return {558            "passed": total_input_energy >= -max(self.tolerance, energy_scale * 0.01),559            "total_input_energy": float(total_input_energy),560            "final_storage": float(S_T),561            "dissipated": float(total_input_energy - S_T),562            "relative_violation": float(relative_violation),563        }564 565    def _test_l1_error(566        self,567        reduced: CrossImpactKernel,568        full: CrossImpactKernel,569        T: float,570        n_points: int = 1000,571    ) -> dict:572        """573        Test 5: L^1 approximation error bound.574        575        ||G - G_r||_{L^1(0,T)} controls:576        - Optimizer difference: ||u_r* - u*||_{L^2} ≤ C * ||G-G_r||_{L^1}577        - Value difference: |V(G_r) - V(G)| ≤ C * ||G-G_r||_{L^1}578        """579        t_grid = np.linspace(0, T, n_points)580        G_full = full.evaluate_grid(t_grid)581        G_red = reduced.evaluate_grid(t_grid)582 583        diff = G_full - G_red584        norms = np.linalg.norm(diff.reshape(n_points, -1), axis=1)585        l1_error = float(np.trapezoid(norms, t_grid))586 587        return {588            "passed": True,  # informational — no hard threshold589            "l1_error": l1_error,590            "max_pointwise_error": float(np.max(norms)),591            "mean_pointwise_error": float(np.mean(norms)),592        }593 594    def _generate_round_trip(self, n_steps: int, d: int) -> np.ndarray:595        """Generate a random round-trip trading program (net zero)."""596        u = np.random.randn(n_steps, d) * 0.1597        # Force round-trip: subtract mean to make sum ≈ 0598        u -= np.mean(u, axis=0, keepdims=True)599        return u600 601    def _compute_impact_cost(602        self,603        kernel: CrossImpactKernel,604        u: np.ndarray,605        t_grid: np.ndarray,606        dt: float,607    ) -> float:608        """609        Compute impact cost C(u) = ⟨u, K_G u⟩ via discrete convolution.610        611        (K_G u)(t) = ∫_0^t G(t-s) u(s) ds ≈ Σ_s G(t-s) u(s) dt612        C(u) = ∫_0^T u(t)' (K_G u)(t) dt613        """614        n = len(t_grid)615        d = u.shape[1]616        cost = 0.0617 618        for t_idx in range(n):619            # Compute (K_G u)(t) via summation620            impact = np.zeros(d)621            for s_idx in range(t_idx + 1):622                lag = t_grid[t_idx] - t_grid[s_idx]623                G_lag = kernel.evaluate(lag)624                impact += G_lag @ u[s_idx] * dt625 626            cost += np.dot(u[t_idx], impact) * dt627 628        return cost629