Race three predictability measures; the literature's fix loses

docs/EVAL.md

@@ -1120,3 +1120,79 @@ LightGBM result — needs narrowing to a statement about decisions rather than a
 No claim was contradicted outright. The Q13 negative result is independently predicted by
 the forecastability literature, which is the strongest single piece of evidence that the
 measurement apparatus here is working correctly.
+
+---
+
+# Racing the predictability measures (2026-08-08)
+
+Q13 showed the shipped diagnostic does not generalise, and the forecastability literature
+offered a specific reason and a specific fix: a single lagged correlation reads one
+frequency, spectral entropy reads all of them. **The fix was tried and it did not work.**
+
+## What was raced
+
+Three training-free measures, scored against the same target — did forecasting pay, by total
+economic cost, at a 6 h and 12 h commitment — on the same 43 serverless workloads from Q13
+plus the 5 non-GPU fleet workloads:
+
+1. **Daily autocorrelation** — what shipped.
+2. **Spectral predictability**, `1 − normalised spectral entropy` over a Welch-averaged
+   periodogram — the field's standard forecastability feature.
+3. **Low-frequency power fraction** — the share of spectral power at periods *longer than
+   the commitment window*. Not from the literature; motivated by this project's own Q1
+   result, that a decision fixed for N hours can only exploit structure slower than N hours.
+   Alone among the three it is a function of the decision horizon, not of the series only.
+
+Comparison is by **AUC**, which needs no cutoff. Reporting each measure at its own best
+threshold would let a measure win by fitting a cutoff to 43 points, so best-threshold
+accuracy is shown only as a clearly-labelled optimistic bound.
+
+## Result
+
+Serverless workloads, n = 43, with a 20,000-draw permutation test on each AUC:
+
+| measure | AUC @ 6 h | p | AUC @ 12 h | p |
+|---|---:|---:|---:|---:|
+| daily autocorrelation | 0.539 | 0.67 | 0.438 | 0.50 |
+| spectral predictability | 0.396 | 0.26 | 0.456 | 0.63 |
+| low-frequency power | 0.595 | 0.31 | **0.700** | **0.027** |
+
+**Spectral entropy failed.** At 0.396 and 0.456 it is at or below a coin flip, and on the
+fleet workloads it is actively anti-predictive (AUC 0.167 at 6 h, though n = 5 there makes
+that number nearly meaningless). The literature's recommended instrument did not rescue the
+diagnostic on this population, and predicting that it would was wrong.
+
+**The horizon-relative measure is the only one showing signal**, and it is not from the
+literature — it comes from this project's own Q1 finding. At a twelve-hour commitment it
+reaches AUC 0.700 with a nominal p = 0.027.
+
+## Why that 0.700 is not being shipped
+
+Six tests were run (3 measures x 2 windows). A Bonferroni-corrected threshold is p < 0.0083,
+and 0.027 does not clear it. Picking the largest of six AUCs and quoting its uncorrected
+p-value is precisely how a result gets manufactured, and this project's own doctrine says an
+implausibly clean number is an artefact until proven otherwise.
+
+So the honest status is: **a promising lead on one cohort, at one horizon, that does not
+survive correction for the number of things tried.** It is not a finding, and the shipped
+diagnostic is unchanged.
+
+## What this means for the diagnostic
+
+The product keeps daily autocorrelation and keeps the scope caveat Q13 forced onto it. That
+is not because autocorrelation is good — on serverless it is a coin flip — but because
+nothing tested beat it well enough to justify swapping the claim a visitor reads.
+
+The measures are implemented, unit-tested against signals with known answers, and available
+in `delphi.forecast.predictability` for the follow-up, which is now well defined: test
+low-frequency power on a much larger cohort with the horizon fixed in advance, as a single
+pre-registered hypothesis rather than one of six. On a fresh cohort with one test, p < 0.05
+would mean something.
+
+## Added to the negatives ledger
+
+- **Spectral entropy did not improve on daily autocorrelation** for predicting whether
+  forecasting pays, on either population — contradicting the expectation drawn from the
+  forecastability literature.
+- **The one measure that showed signal does not survive multiple-comparison correction**, and
+  is reported as a lead rather than promoted into the product.

scripts/compare_predictability.py

@@ -0,0 +1,164 @@
+"""Which training-free measure actually predicts whether forecasting pays?
+
+Q13 established that the shipped diagnostic — daily autocorrelation — does not generalise
+past the fleet-aggregate traces it was derived on. The forecastability literature says why:
+a single lagged correlation reads one frequency, while spectral entropy reads all of them,
+and spectral measures are the established indicator for exactly this question.
+
+This races three candidates against the same target on the same workloads:
+
+* ``daily_autocorrelation`` — what shipped;
+* ``spectral_predictability`` (1 - spectral entropy) — the field's standard;
+* ``low_frequency_power_fraction`` — the share of power at periods longer than the
+  commitment window, motivated by this project's own Q1 result that what a fixed-for-N-hours
+  decision can exploit is structure slower than N hours.
+
+**Comparison is threshold-free.** Reporting each measure at its own best cutoff would let a
+measure win by overfitting a cutoff to 43 points. AUC — the probability that a randomly
+chosen paying workload scores above a randomly chosen non-paying one — needs no cutoff and
+is the honest way to ask which measure carries more signal. Best-threshold accuracy is
+printed alongside as an optimistic upper bound, clearly labelled as one.
+
+Run: `python scripts/compare_predictability.py --cohort 400 --sample 60`
+"""
+
+from __future__ import annotations
+
+import argparse
+from dataclasses import dataclass
+
+import numpy as np
+from evaluate_diagnostic import build_workloads  # type: ignore[import-not-found]
+from evaluate_threshold import (  # type: ignore[import-not-found]
+    BIN_SECONDS,
+    build_candidates,
+    stratify,
+    strict_dominance,
+)
+
+from delphi.forecast.predictability import (
+    daily_autocorrelation,
+    low_frequency_power_fraction,
+    spectral_predictability,
+)
+
+
+@dataclass(frozen=True)
+class Scored:
+    name: str
+    population: str
+    measures: dict[str, float]
+    paid: bool
+
+
+def auc(scores: list[float], labels: list[bool]) -> float:
+    """Mann-Whitney U / ROC AUC. 0.5 is a coin flip; below 0.5 is anti-predictive."""
+    positive = [s for s, y in zip(scores, labels, strict=True) if y]
+    negative = [s for s, y in zip(scores, labels, strict=True) if not y]
+    if not positive or not negative:
+        return float("nan")
+    order = np.argsort(np.asarray(scores, dtype=np.float64), kind="mergesort")
+    ranks = np.empty(len(scores), dtype=np.float64)
+    ranks[order] = np.arange(1, len(scores) + 1, dtype=np.float64)
+    # Average ranks within ties so a constant measure scores exactly 0.5, not 1.0.
+    values = np.asarray(scores, dtype=np.float64)
+    for value in np.unique(values):
+        mask = values == value
+        if mask.sum() > 1:
+            ranks[mask] = ranks[mask].mean()
+    positive_rank_sum = float(ranks[np.asarray(labels, dtype=bool)].sum())
+    n_pos, n_neg = len(positive), len(negative)
+    return (positive_rank_sum - n_pos * (n_pos + 1) / 2.0) / (n_pos * n_neg)
+
+
+def best_threshold_accuracy(scores: list[float], labels: list[bool]) -> tuple[float, float]:
+    """Optimistic upper bound: the best cutoff chosen *on this same data*."""
+    best = (float("nan"), 0.0)
+    for cut in sorted(set(scores)):
+        predicted = [s >= cut for s in scores]
+        acc = sum(p == y for p, y in zip(predicted, labels, strict=True)) / len(labels)
+        if acc > best[1]:
+            best = (cut, acc)
+    return best
+
+
+def measures_for(values: np.ndarray, day_steps: int, window_steps: int) -> dict[str, float]:
+    return {
+        "daily autocorr": daily_autocorrelation(values, day_steps),
+        "spectral predictability": spectral_predictability(values),
+        "low-freq power": low_frequency_power_fraction(values, window_steps),
+    }
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--cohort", type=int, default=400)
+    parser.add_argument("--sample", type=int, default=60)
+    parser.add_argument("--windows", nargs="+", type=int, default=[6, 12])
+    parser.add_argument("--seed", type=int, default=20260808)
+    args = parser.parse_args()
+
+    print("# Which predictability measure actually works?\n")
+    print(
+        "Threshold-free comparison by AUC — the probability a paying workload scores above "
+        "a non-paying one. 0.50 is a coin flip. Best-threshold accuracy is an optimistic "
+        "upper bound, since the cutoff is fitted on the same data it is scored on.\n"
+    )
+
+    serverless = stratify(build_candidates(args.cohort, args.seed), args.sample, args.seed)
+    fleet = build_workloads(include_gpu=False)
+
+    for window_hours in args.windows:
+        rows: list[Scored] = []
+
+        window_steps = window_hours * 3600 // BIN_SECONDS
+        for candidate in serverless:
+            workload = candidate.workload
+            wins, _ = strict_dominance(workload, window_steps)
+            rows.append(
+                Scored(
+                    workload.name,
+                    "serverless",
+                    measures_for(workload.series.values, workload.day_steps, window_steps),
+                    wins > 0,
+                )
+            )
+
+        for workload in fleet:
+            steps = window_hours * 3600 // workload.series.step_seconds
+            wins, _ = strict_dominance(workload, steps)
+            rows.append(
+                Scored(
+                    workload.name,
+                    "fleet",
+                    measures_for(workload.series.values, workload.day_steps, steps),
+                    wins > 0,
+                )
+            )
+
+        print(f"\n## {window_hours}-hour commitment\n")
+        for population in ("serverless", "fleet", "both"):
+            subset = [r for r in rows if population in (r.population, "both")]
+            labels = [r.paid for r in subset]
+            if len(set(labels)) < 2:
+                print(
+                    f"### {population}: all outcomes identical, AUC undefined (n={len(subset)})\n"
+                )
+                continue
+            print(
+                f"### {population} (n={len(subset)}, "
+                f"{sum(labels)} paid / {len(labels) - sum(labels)} did not)\n"
+            )
+            print("| measure | AUC | best-threshold accuracy (optimistic) |")
+            print("|---|---:|---:|")
+            for key in ("daily autocorr", "spectral predictability", "low-freq power"):
+                scores = [r.measures[key] for r in subset]
+                area = auc(scores, labels)
+                cut, acc = best_threshold_accuracy(scores, labels)
+                print(f"| {key} | {area:.3f} | {acc:.0%} at >= {cut:.3f} |")
+            baseline = max(sum(labels), len(labels) - sum(labels)) / len(labels)
+            print(f"\nAlways-guess-the-majority baseline: {baseline:.0%}.\n")
+
+
+if __name__ == "__main__":
+    main()

src/delphi/api/app.py

@@ -32,6 +32,7 @@
 )
 from delphi.config import get_settings
 from delphi.control.newsvendor import CostRatio
+from delphi.forecast.predictability import optional_lag_autocorrelation
 
 SNAPSHOT_ENV = "DELPHI_SNAPSHOT_PATH"
 DEFAULT_SNAPSHOT = Path("data/snapshot.json")
@@ -190,12 +191,8 @@ def score(request: Annotated[DiagnosticRequest, Body()]) -> DiagnosticResponse:
             status_code=422, detail="a constant series has no structure to diagnose"
         )
 
-    centred = values - values.mean()
-
     def autocorrelation(lag: int) -> float | None:
-        if lag < 1 or len(values) <= lag + 8:
-            return None
-        return float(np.corrcoef(centred[:-lag], centred[lag:])[0, 1])
+        return optional_lag_autocorrelation(values, lag)
 
     per_day = max(round(86400 / request.step_seconds), 1)
     daily = autocorrelation(per_day)