criptic1/ctrader-ml-systems

1"""2Calibrated uncertainty quantification pipeline for the Aurea allocator.3 4Produces σ̂ (the uncertainty penalty term) from three sources:51. Ensemble disagreement across TimesFM, Chronos-2, Kronos, LightGBM/XGBoost62. Probabilistic model outputs (Chronos-2 calibrated quantiles, Kronos MC samples)73. Conformal prediction intervals (SFOCP: scale-free online conformal prediction)8 9The combined σ̂ feeds directly into the allocator's scoring:10    net_utility = edge - κ_u * σ̂ - friction11 12References:13- arxiv:2302.07869 — SAOCP / Scale-Free Online Conformal Prediction14- arxiv:2212.03463 — SPCI (Sequential Predictive Conformal Inference)15- arxiv:2508.02686 — MoE volatility-sensitive routing16"""17 18import logging19from typing import Optional, Dict, List, Tuple20import numpy as np21from dataclasses import dataclass, field22 23from ..utils.types import ForecastResult, UncertaintyEstimate24 25logger = logging.getLogger(__name__)26 27 28class ScaleFreeOnlineConformal:29    """30    Adaptive Online Conformal Prediction (ACI variant).31    32    Inspired by arxiv:2302.07869 (SAOCP) but uses the practical ACI33    algorithm from Gibbs & Candes (2021) which has better finite-sample34    convergence. The key idea:35    36        α̂_{t+1} = α̂_t + γ * (α - 𝟙[y_t ∉ C_t(x_t)])37        38    where α̂_t is the adaptive miscoverage level and C_t is the prediction39    set at the α̂_t quantile of recent residuals. This drives empirical40    coverage toward (1-α) without the gradient-magnitude asymmetry of41    raw pinball-loss OGD.42    43    For the radius-based formulation:44        s_hat tracks the (1-α̂_t) quantile of a rolling window of residuals,45        with α̂ adapted online to match the target coverage.46    """47 48    def __init__(self, alpha: float = 0.1, eta: float = 0.05, window: int = 200):49        """50        Args:51            alpha: miscoverage rate (0.1 = 90% coverage target)52            eta: step size for α̂ adaptation (0.01-0.1 typical)53            window: rolling window for residual quantile estimation54        """55        self.alpha = alpha56        self.eta = eta57        self.window = window58        self.s_hat = 0.0               # predicted radius59        self.alpha_hat = alpha          # adaptive miscoverage level60        self.n_updates = 061        self._residuals: List[float] = []62        self._coverage_history: List[bool] = []63        self._width_history: List[float] = []64 65    def predict_interval(self, y_hat: float) -> Tuple[float, float]:66        """Return (lower, upper) prediction interval centered on y_hat."""67        return y_hat - self.s_hat, y_hat + self.s_hat68 69    def predict_interval_array(self, y_hat: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:70        """Return (lower, upper) arrays for a forecast vector."""71        return y_hat - self.s_hat, y_hat + self.s_hat72 73    def update(self, y_true: float, y_hat: float):74        """75        Online update after observing true value.76        77        Two-step process:78        1. Record residual and whether current interval covered y_true79        2. Adapt α̂ and recompute s_hat from residual quantile80        """81        S_t = abs(y_true - y_hat)  # absolute residual82 83        # Track coverage BEFORE update84        covered = S_t <= self.s_hat85        self._coverage_history.append(covered)86        self._width_history.append(2 * self.s_hat)87 88        # Store residual89        self._residuals.append(S_t)90        if len(self._residuals) > self.window:91            self._residuals = self._residuals[-self.window:]92 93        # ACI update (Gibbs & Candes 2021):94        # α̂_{t+1} = α̂_t + γ * (α - miss_t)95        # miss_t = 1 if y_t NOT covered, 0 if covered96        # When missing too often: α̂ decreases → quantile level (1-α̂) increases → wider interval97        # When covering too much: α̂ increases → quantile level (1-α̂) decreases → narrower interval98        miss_t = 0.0 if covered else 1.099        self.alpha_hat = self.alpha_hat + self.eta * (self.alpha - miss_t)100        self.alpha_hat = np.clip(self.alpha_hat, 0.001, 0.999)101 102        # Recompute s_hat as (1 - α̂) quantile of recent residuals103        if len(self._residuals) >= 5:104            self.s_hat = float(np.quantile(105                self._residuals, 1 - self.alpha_hat106            ))107        else:108            # Warmup: use the max residual seen so far109            self.s_hat = max(self._residuals) if self._residuals else 0.0110 111        self.n_updates += 1112 113    def batch_update(self, y_true_arr: np.ndarray, y_hat_arr: np.ndarray):114        """Update with a batch of observations."""115        for yt, yh in zip(y_true_arr, y_hat_arr):116            self.update(float(yt), float(yh))117 118    @property119    def empirical_coverage(self) -> float:120        """Empirical coverage rate over update history."""121        if not self._coverage_history:122            return 1.0 - self.alpha123        return float(np.mean(self._coverage_history[-100:]))124 125    @property126    def avg_width(self) -> float:127        """Average interval width over recent history."""128        if not self._width_history:129            return 0.0130        return float(np.mean(self._width_history[-100:]))131 132    @property133    def calibration_error(self) -> float:134        """|empirical_coverage - (1-alpha)| — 0 is perfect."""135        return abs(self.empirical_coverage - (1 - self.alpha))136 137 138class EnsembleDisagreement:139    """140    Compute uncertainty from disagreement across multiple forecast models.141    142    Uses normalized standard deviation of forecasts, weighted by143    recent model performance (inverse MSE weighting).144    """145 146    def __init__(self, normalize: bool = True):147        self.normalize = normalize148 149    def compute(150        self,151        forecasts: Dict[str, np.ndarray],152        weights: Optional[Dict[str, float]] = None,153    ) -> Tuple[float, np.ndarray]:154        """155        Compute ensemble disagreement.156        157        Args:158            forecasts: model_name -> point_forecast array159            weights: model_name -> weight (optional, uniform if None)160            161        Returns:162            (scalar_disagreement, per_step_disagreement)163        """164        if len(forecasts) < 2:165            return 0.0, np.zeros(1)166 167        # Align lengths168        min_len = min(len(f) for f in forecasts.values())169        aligned = {k: v[:min_len] for k, v in forecasts.items()}170 171        if weights is None:172            weights = {k: 1.0 / len(forecasts) for k in forecasts}173 174        # Normalize weights175        total_w = sum(weights.get(k, 1.0) for k in aligned)176        norm_weights = {k: weights.get(k, 1.0) / total_w for k in aligned}177 178        # Weighted mean179        stacked = np.array(list(aligned.values()))180        w_arr = np.array([norm_weights[k] for k in aligned])181        weighted_mean = np.average(stacked, axis=0, weights=w_arr)182 183        # Weighted std (disagreement)184        sq_diffs = w_arr[:, None] * (stacked - weighted_mean[None, :]) ** 2185        per_step_var = np.sum(sq_diffs, axis=0)186        per_step_std = np.sqrt(per_step_var)187 188        if self.normalize and np.mean(np.abs(weighted_mean)) > 1e-10:189            # Coefficient of variation (normalized by level)190            per_step_std = per_step_std / (np.abs(weighted_mean) + 1e-10)191 192        scalar = float(np.mean(per_step_std))193        return scalar, per_step_std194 195 196class QuantileCalibrator:197    """198    Online calibration of model-produced quantile forecasts.199    200    Tracks empirical coverage for each quantile level and adjusts201    the prediction intervals via isotonic regression / histogram binning.202    """203 204    def __init__(self, quantile_levels: List[float] = None):205        self.quantile_levels = quantile_levels or [0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95]206        self._hits: Dict[float, List[bool]] = {q: [] for q in self.quantile_levels}207        self._adjustments: Dict[float, float] = {q: 0.0 for q in self.quantile_levels}208 209    def update(self, y_true: float, quantile_forecasts: Dict[float, float]):210        """211        Update calibration with observed value.212        213        Args:214            y_true: realized value215            quantile_forecasts: q -> forecast_at_q216        """217        for q in self.quantile_levels:218            if q in quantile_forecasts:219                hit = y_true <= quantile_forecasts[q]220                self._hits[q].append(hit)221                # Keep rolling window222                if len(self._hits[q]) > 200:223                    self._hits[q] = self._hits[q][-200:]224 225                # Compute adjustment (simple additive correction)226                empirical_q = np.mean(self._hits[q])227                self._adjustments[q] = q - empirical_q  # positive = need wider228 229    def adjust_quantiles(230        self, quantile_forecasts: Dict[float, np.ndarray]231    ) -> Dict[float, np.ndarray]:232        """Apply calibration adjustments to quantile forecasts."""233        adjusted = {}234        for q, forecast in quantile_forecasts.items():235            if q in self._adjustments:236                # Scale adjustment by the median level237                median = quantile_forecasts.get(0.5, forecast)238                scale = np.abs(median) * 0.01 + 1e-6239                adjusted[q] = forecast + self._adjustments[q] * scale240            else:241                adjusted[q] = forecast242        return adjusted243 244    def get_calibration_scores(self) -> Dict[float, float]:245        """Return calibration error for each quantile level."""246        scores = {}247        for q in self.quantile_levels:248            if self._hits[q]:249                empirical = np.mean(self._hits[q])250                scores[q] = abs(empirical - q)251            else:252                scores[q] = 0.0253        return scores254 255    @property256    def mean_calibration_error(self) -> float:257        """Mean absolute calibration error across all quantile levels."""258        scores = self.get_calibration_scores()259        if not scores:260            return 0.0261        return float(np.mean(list(scores.values())))262 263 264class CalibratedUncertaintyPipeline:265    """266    Full uncertainty quantification pipeline.267    268    Combines:269    1. Ensemble disagreement (across TimesFM/Chronos-2/Kronos/LightGBM)270    2. Probabilistic model quantiles (Chronos-2, Kronos MC)271    3. Conformal prediction (SFOCP for online adaptation)272    4. Quantile calibration (online isotonic correction)273    274    Produces a single σ̂ per symbol for the Aurea allocator, plus275    detailed prediction intervals for risk management.276    """277 278    def __init__(279        self,280        alpha: float = 0.1,281        ensemble_weight: float = 0.3,282        quantile_weight: float = 0.35,283        conformal_weight: float = 0.35,284    ):285        """286        Args:287            alpha: miscoverage rate for prediction intervals288            ensemble_weight: weight for ensemble disagreement in σ̂289            quantile_weight: weight for model quantile width in σ̂290            conformal_weight: weight for conformal interval width in σ̂291        """292        self.alpha = alpha293        self.w_ensemble = ensemble_weight294        self.w_quantile = quantile_weight295        self.w_conformal = conformal_weight296 297        # Per-symbol conformal predictors298        self._conformal: Dict[str, ScaleFreeOnlineConformal] = {}299        self._calibrator = QuantileCalibrator()300        self._disagreement = EnsembleDisagreement(normalize=True)301 302    def _get_conformal(self, symbol: str) -> ScaleFreeOnlineConformal:303        """Get or create conformal predictor for a symbol."""304        if symbol not in self._conformal:305            self._conformal[symbol] = ScaleFreeOnlineConformal(306                alpha=self.alpha, eta=1.0307            )308        return self._conformal[symbol]309 310    def compute_uncertainty(311        self,312        symbol: str,313        primary_forecast: ForecastResult,314        ensemble_forecasts: Optional[List[ForecastResult]] = None,315        model_weights: Optional[Dict[str, float]] = None,316    ) -> UncertaintyEstimate:317        """318        Compute calibrated uncertainty estimate for a symbol.319        320        This is the main entry point. Call once per symbol per cycle.321        322        Args:323            symbol: ticker324            primary_forecast: the forecast from the regime-gated model325            ensemble_forecasts: forecasts from all models (for disagreement)326            model_weights: adaptive weights from performance tracking327            328        Returns:329            UncertaintyEstimate with σ̂ and prediction intervals330        """331        components = {}332        horizon = primary_forecast.horizon333 334        # ─── 1. Ensemble disagreement ──────────────────────────────335        ensemble_disagreement = 0.0336        if ensemble_forecasts and len(ensemble_forecasts) > 0:337            forecasts_dict = {338                primary_forecast.model.value: primary_forecast.point_forecast339            }340            for ef in ensemble_forecasts:341                forecasts_dict[ef.model.value] = ef.point_forecast342 343            ensemble_disagreement, per_step_disagree = self._disagreement.compute(344                forecasts_dict, model_weights345            )346        components["ensemble_disagreement"] = ensemble_disagreement347 348        # ─── 2. Model quantile width ──────────────────────────────349        quantile_width = 0.0350        if primary_forecast.quantile_forecasts:351            # Apply calibration adjustment352            adjusted_q = self._calibrator.adjust_quantiles(353                primary_forecast.quantile_forecasts354            )355 356            # Compute interval width at the target alpha357            lo_key = min(adjusted_q.keys(), key=lambda q: abs(q - self.alpha / 2))358            hi_key = min(adjusted_q.keys(), key=lambda q: abs(q - (1 - self.alpha / 2)))359 360            lo_arr = adjusted_q[lo_key]361            hi_arr = adjusted_q[hi_key]362 363            width_arr = hi_arr - lo_arr364            quantile_width = float(np.mean(np.abs(width_arr)))365 366            # Normalize by price level367            price_level = np.abs(np.mean(primary_forecast.point_forecast))368            if price_level > 1e-10:369                quantile_width = quantile_width / price_level370 371        components["quantile_width"] = quantile_width372 373        # ─── 3. Conformal prediction interval ─────────────────────374        conformal = self._get_conformal(symbol)375        conformal_width = conformal.avg_width376 377        # Normalize by price level378        price_level = np.abs(np.mean(primary_forecast.point_forecast))379        if price_level > 1e-10:380            conformal_width_norm = conformal_width / price_level381        else:382            conformal_width_norm = conformal_width383 384        components["conformal_width"] = conformal_width_norm385        components["conformal_coverage"] = conformal.empirical_coverage386        components["conformal_calibration_error"] = conformal.calibration_error387 388        # ─── 4. Combined σ̂ ──────────────────────────────────────389        sigma_hat = (390            self.w_ensemble * ensemble_disagreement391            + self.w_quantile * quantile_width392            + self.w_conformal * conformal_width_norm393        )394 395        # ─── 5. Build prediction intervals ────────────────────────396        prediction_intervals = {}397        for test_alpha in [0.05, 0.1, 0.2]:398            lo, hi = conformal.predict_interval_array(primary_forecast.point_forecast)399            # Scale by alpha ratio400            scale = (self.alpha / test_alpha) if test_alpha < self.alpha else 1.0401            center = primary_forecast.point_forecast402            prediction_intervals[test_alpha] = (403                center - (center - lo) * scale,404                center + (hi - center) * scale,405            )406 407        # ─── 6. Overall calibration score ─────────────────────────408        calibration_score = (409            conformal.calibration_error * 0.5410            + self._calibrator.mean_calibration_error * 0.5411        )412 413        return UncertaintyEstimate(414            symbol=symbol,415            sigma_hat=float(sigma_hat),416            prediction_intervals=prediction_intervals,417            ensemble_disagreement=float(ensemble_disagreement),418            conformal_width=float(conformal_width),419            calibration_score=float(calibration_score),420            components=components,421        )422 423    def update_with_realized(424        self,425        symbol: str,426        y_true: float,427        y_predicted: float,428        quantile_forecasts: Optional[Dict[float, float]] = None,429    ):430        """431        Update all calibration components with realized values.432        Call this after each cycle when actual prices are known.433        434        Args:435            symbol: ticker436            y_true: realized price437            y_predicted: point forecast that was made438            quantile_forecasts: quantile forecasts that were made (for calibrator)439        """440        # Update conformal predictor441        conformal = self._get_conformal(symbol)442        conformal.update(y_true, y_predicted)443 444        # Update quantile calibrator445        if quantile_forecasts:446            self._calibrator.update(y_true, quantile_forecasts)447 448    def get_diagnostics(self) -> Dict:449        """Return diagnostic information for monitoring."""450        diag = {451            "n_symbols_tracked": len(self._conformal),452            "quantile_calibration": self._calibrator.get_calibration_scores(),453            "mean_calibration_error": self._calibrator.mean_calibration_error,454        }455        for symbol, conf in self._conformal.items():456            diag[f"conformal_{symbol}"] = {457                "coverage": conf.empirical_coverage,458                "avg_width": conf.avg_width,459                "calibration_error": conf.calibration_error,460                "n_updates": conf.n_updates,461            }462        return diag463