API reference

skepsis.evaluate.evaluate(returns, trials=None, params=None, freq='daily', chosen=None, pbo_blocks=16, n_resamples=5000, seed=0, thresholds=None)

Run every overfitting diagnostic the provided inputs allow. See README.

Source code in skepsis/evaluate.py

def evaluate(
    returns: Any,
    trials: Any | None = None,
    params: Any | None = None,
    freq: str | int | float = "daily",
    chosen: int | str | None = None,
    pbo_blocks: int = 16,
    n_resamples: int = 5000,
    seed: int = 0,
    thresholds: Thresholds | None = None,
) -> Result:
    """Run every overfitting diagnostic the provided inputs allow. See README."""
    periods = moments.periods_per_year(freq)
    validate_n_blocks(pbo_blocks)
    r = coerce_returns(returns)
    trials_arr: np.ndarray | None = None
    labels: list[str] = []
    params_arr: np.ndarray | None = None
    param_names: list[str] | None = None
    if trials is not None:
        trials_arr, labels = coerce_trials(trials)
    if params is not None:
        params_arr, param_names = coerce_params(params)
    validate_alignment(r, trials_arr, params_arr)

    skipped: dict[str, str] = {}
    with _warnings.catch_warnings(record=True) as caught:
        _warnings.simplefilter("always")

        if chosen is not None and params_arr is None:
            _warnings.warn(
                "`chosen=` was provided without `params=`; it only affects the "
                "sensitivity diagnostic and is ignored",
                SkepsisWarning,
                stacklevel=2,
            )

        n_obs = len(r)
        sr_p = moments.sharpe(r)
        skew = moments.skewness(r)
        kurt = moments.kurtosis(r)
        psr_value = probabilistic_sharpe_ratio(sr_p, 0.0, n_obs, skew, kurt)
        psr_res = PsrResult(
            value=psr_value,
            sharpe_periodic=sr_p,
            sharpe_annualized=moments.annualized_sharpe(r, periods),
            benchmark_sr=0.0,
            n_obs=n_obs,
            skewness=skew,
            kurtosis=kurt,
        )

        trial_metrics: np.ndarray | None = None
        chosen_label: str | None = None
        if trials_arr is not None:
            n_trials = trials_arr.shape[1]
            col_sharpes = _column_sharpes(trials_arr)
            finite = col_sharpes[np.isfinite(col_sharpes)]
            var_sr = float(np.var(finite, ddof=1)) if len(finite) >= 2 else 0.0
            if var_sr == 0.0:
                _warnings.warn(
                    "variance of trial Sharpe ratios is zero; DSR benchmark falls "
                    "back to 0 (no deflation)",
                    SkepsisWarning,
                    stacklevel=2,
                )
            single_trial = False
        else:
            n_trials, var_sr, single_trial = 1, 0.0, True
            _warnings.warn(
                "no trials provided — DSR uses trial count 1, which is almost "
                "certainly optimistic; pass `trials=` with every variant you tried",
                SkepsisWarning,
                stacklevel=2,
            )
        dsr_value = deflated_sharpe_ratio(sr_p, n_obs, skew, kurt, var_sr, n_trials)
        dsr_res = DsrResult(
            value=dsr_value,
            p_value=1.0 - dsr_value,
            benchmark_sr=expected_max_sharpe(var_sr, n_trials),
            n_trials=n_trials,
            var_trial_sr=var_sr,
            single_trial=single_trial,
        )

        pbo_res: PboResult | None = None
        if trials_arr is None:
            skipped["pbo"] = "trials not provided (pass the returns of every variant tried)"
        elif n_obs < 2 * pbo_blocks:
            skipped["pbo"] = (
                f"needs >= {2 * pbo_blocks} observations for n_blocks={pbo_blocks}, "
                f"got {n_obs}; pass a smaller pbo_blocks explicitly"
            )
        else:
            pbo_res = cscv(trials_arr, n_blocks=pbo_blocks)

        boot_res: BootstrapResult | None = None
        if n_obs < _MIN_BOOTSTRAP_OBS:
            skipped["bootstrap"] = f"needs >= {_MIN_BOOTSTRAP_OBS} observations, got {n_obs}"
        else:
            boot_res = bootstrap(r, periods, n_resamples=n_resamples, seed=seed)
            if boot_res.mean_block_length > _AUTOCORR_WARN_BLOCK_LENGTH:
                _warnings.warn(
                    f"estimated mean block length {boot_res.mean_block_length:.1f} "
                    f"exceeds {_AUTOCORR_WARN_BLOCK_LENGTH:.0f}: returns are heavily "
                    "autocorrelated, which strains the IID-ish assumptions behind "
                    "PSR/DSR — read those diagnostics with extra skepticism",
                    SkepsisWarning,
                    stacklevel=2,
                )

        sens_res: SensitivityResult | None = None
        if params_arr is None:
            skipped["sensitivity"] = "params not provided (one row per trial)"
        else:
            assert trials_arr is not None  # validate_alignment guarantees this
            metrics = _column_sharpes(trials_arr, warn=False) * float(np.sqrt(periods))
            chosen_idx, chosen_label = _find_chosen(r, trials_arr, labels, chosen, metrics)
            trial_metrics = metrics
            sens_res = sensitivity(params_arr, metrics, chosen_idx)

    warning_msgs = [str(w.message) for w in caught if issubclass(w.category, SkepsisWarning)]
    for w in caught:
        _warnings.warn_explicit(w.message, w.category, w.filename, w.lineno)

    verdict = decide(
        dsr=dsr_res.value,
        dsr_single_trial=dsr_res.single_trial,
        pbo=pbo_res.value if pbo_res else None,
        bootstrap_p=boot_res.p_value_no_skill if boot_res else None,
        stability_score=sens_res.stability_score if sens_res else None,
        thresholds=thresholds,
    )
    return Result(
        psr=psr_res,
        deflated_sharpe=dsr_res,
        pbo=pbo_res,
        bootstrap=boot_res,
        sensitivity=sens_res,
        verdict=verdict,
        skipped=skipped,
        warnings=warning_msgs,
        meta={
            "freq": freq,
            "periods_per_year": periods,
            "n_obs": n_obs,
            "n_trials": n_trials,
            "chosen_label": chosen_label,
            "skepsis_version": __version__,
        },
        returns=r,
        params=params_arr,
        param_names=param_names,
        trial_metrics=trial_metrics,
        thresholds=thresholds or Thresholds(),
    )

skepsis.evaluate.Result dataclass

Everything skepsis concluded, plus the data the HTML report needs.

Source code in skepsis/evaluate.py

@dataclass
class Result:
    """Everything skepsis concluded, plus the data the HTML report needs."""

    psr: PsrResult
    deflated_sharpe: DsrResult
    pbo: PboResult | None
    bootstrap: BootstrapResult | None
    sensitivity: SensitivityResult | None
    verdict: Verdict
    skipped: dict[str, str]
    warnings: list[str]
    meta: dict[str, Any]
    returns: np.ndarray
    params: np.ndarray | None = None
    param_names: list[str] | None = None
    trial_metrics: np.ndarray | None = None
    thresholds: Thresholds = field(default_factory=Thresholds)

    def to_dict(self) -> dict[str, Any]:
        """JSON-serializable summary (scalars only, no arrays)."""
        d: dict[str, Any] = {
            "psr": {
                "value": self.psr.value,
                "sharpe_periodic": self.psr.sharpe_periodic,
                "sharpe_annualized": self.psr.sharpe_annualized,
                "n_obs": self.psr.n_obs,
                "skewness": self.psr.skewness,
                "kurtosis": self.psr.kurtosis,
            },
            "deflated_sharpe": {
                "value": self.deflated_sharpe.value,
                "p_value": self.deflated_sharpe.p_value,
                "benchmark_sr": self.deflated_sharpe.benchmark_sr,
                "n_trials": self.deflated_sharpe.n_trials,
                "single_trial": self.deflated_sharpe.single_trial,
            },
            "verdict": {"level": self.verdict.level, "reasons": list(self.verdict.reasons)},
            "skipped": dict(self.skipped),
            "warnings": list(self.warnings),
            "meta": dict(self.meta),
        }
        if self.pbo is not None:
            d["pbo"] = {
                "value": self.pbo.value,
                "n_combinations": self.pbo.n_combinations,
                "n_blocks": self.pbo.n_blocks,
                "n_trials": self.pbo.n_trials,
            }
        if self.bootstrap is not None:
            d["bootstrap"] = {
                "sharpe_obs": self.bootstrap.sharpe_obs,
                "sharpe_ci": list(self.bootstrap.sharpe_ci),
                "drawdown_obs": self.bootstrap.drawdown_obs,
                "drawdown_ci": list(self.bootstrap.drawdown_ci),
                "p_value_no_skill": self.bootstrap.p_value_no_skill,
                "mean_block_length": self.bootstrap.mean_block_length,
                "n_degenerate_resamples": self.bootstrap.n_degenerate_resamples,
            }
        if self.sensitivity is not None:
            d["sensitivity"] = {
                "stability_score": _finite_or_none(self.sensitivity.stability_score),
                "neighbor_median": _finite_or_none(self.sensitivity.neighbor_median),
                "k": self.sensitivity.k,
                "flagged": self.sensitivity.flagged,
            }
        return d

    def summary(self) -> str:
        lines = [
            f"skepsis {self.meta['skepsis_version']} — verdict: {self.verdict.level}",
            f"  annualized Sharpe: {self.psr.sharpe_annualized:.3f}  "
            f"(PSR {self.psr.value:.3f}, DSR {self.deflated_sharpe.value:.3f} "
            f"over {count_trials(self.deflated_sharpe.n_trials)})",
        ]
        if self.pbo is not None:
            lines.append(f"  PBO: {self.pbo.value:.3f} ({self.pbo.n_combinations} combinations)")
        if self.bootstrap is not None:
            lines.append(
                f"  bootstrap no-skill p: {self.bootstrap.p_value_no_skill:.4f}, "
                f"Sharpe 95% CI [{self.bootstrap.sharpe_ci[0]:.2f}, "
                f"{self.bootstrap.sharpe_ci[1]:.2f}]"
            )
        if self.sensitivity is not None:
            lines.append(
                f"  stability score: {format_stability(self.sensitivity.stability_score)}"
            )
        for reason in self.verdict.reasons:
            lines.append(f"  - {reason}")
        for name, why in self.skipped.items():
            lines.append(f"  skipped {name}: {why}")
        return "\n".join(lines)

    def save_html(self, path: str | Path) -> Path:
        """Render the self-contained HTML report (report module lands in Task 10)."""
        from skepsis.report.render import render_html

        out = Path(path)
        out.write_text(render_html(self), encoding="utf-8")
        return out

save_html(path)

Render the self-contained HTML report (report module lands in Task 10).

Source code in skepsis/evaluate.py

def save_html(self, path: str | Path) -> Path:
    """Render the self-contained HTML report (report module lands in Task 10)."""
    from skepsis.report.render import render_html

    out = Path(path)
    out.write_text(render_html(self), encoding="utf-8")
    return out

to_dict()

JSON-serializable summary (scalars only, no arrays).

Source code in skepsis/evaluate.py

def to_dict(self) -> dict[str, Any]:
    """JSON-serializable summary (scalars only, no arrays)."""
    d: dict[str, Any] = {
        "psr": {
            "value": self.psr.value,
            "sharpe_periodic": self.psr.sharpe_periodic,
            "sharpe_annualized": self.psr.sharpe_annualized,
            "n_obs": self.psr.n_obs,
            "skewness": self.psr.skewness,
            "kurtosis": self.psr.kurtosis,
        },
        "deflated_sharpe": {
            "value": self.deflated_sharpe.value,
            "p_value": self.deflated_sharpe.p_value,
            "benchmark_sr": self.deflated_sharpe.benchmark_sr,
            "n_trials": self.deflated_sharpe.n_trials,
            "single_trial": self.deflated_sharpe.single_trial,
        },
        "verdict": {"level": self.verdict.level, "reasons": list(self.verdict.reasons)},
        "skipped": dict(self.skipped),
        "warnings": list(self.warnings),
        "meta": dict(self.meta),
    }
    if self.pbo is not None:
        d["pbo"] = {
            "value": self.pbo.value,
            "n_combinations": self.pbo.n_combinations,
            "n_blocks": self.pbo.n_blocks,
            "n_trials": self.pbo.n_trials,
        }
    if self.bootstrap is not None:
        d["bootstrap"] = {
            "sharpe_obs": self.bootstrap.sharpe_obs,
            "sharpe_ci": list(self.bootstrap.sharpe_ci),
            "drawdown_obs": self.bootstrap.drawdown_obs,
            "drawdown_ci": list(self.bootstrap.drawdown_ci),
            "p_value_no_skill": self.bootstrap.p_value_no_skill,
            "mean_block_length": self.bootstrap.mean_block_length,
            "n_degenerate_resamples": self.bootstrap.n_degenerate_resamples,
        }
    if self.sensitivity is not None:
        d["sensitivity"] = {
            "stability_score": _finite_or_none(self.sensitivity.stability_score),
            "neighbor_median": _finite_or_none(self.sensitivity.neighbor_median),
            "k": self.sensitivity.k,
            "flagged": self.sensitivity.flagged,
        }
    return d

skepsis.verdict.Thresholds dataclass

Default rule thresholds. fail => LIKELY_OVERFIT; warn counts toward WEAK/MODERATE.

Source code in skepsis/verdict.py

@dataclass(frozen=True)
class Thresholds:
    """Default rule thresholds. fail => LIKELY_OVERFIT; warn counts toward WEAK/MODERATE."""

    dsr_fail: float = 0.5    # DSR below this: worse than a coin flip after deflation
    dsr_warn: float = 0.95
    pbo_fail: float = 0.5    # in-sample winner lands bottom-half OOS more often than not
    pbo_warn: float = 0.2
    bootstrap_fail: float = 0.5   # no-skill p-value above this
    bootstrap_warn: float = 0.05
    sensitivity_fail: float = 2.0  # chosen config > 2x its neighbors' median
    sensitivity_warn: float = 1.5

skepsis.core.psr

Probabilistic and Deflated Sharpe Ratio.

References: - Bailey & Lopez de Prado (2012), "The Sharpe Ratio Efficient Frontier", Journal of Risk 15(2). [PSR] - Bailey & Lopez de Prado (2014), "The Deflated Sharpe Ratio: Correcting for Selection Bias, Backtest Overfitting and Non-Normality", Journal of Portfolio Management 40(5). [DSR]

Conventions: sr and sr_benchmark are PERIODIC (non-annualized) Sharpe ratios; kurt is NON-excess kurtosis (normal = 3.0).

deflated_sharpe_ratio(sr, n_obs, skew, kurt, var_trial_sr, n_trials)

PSR evaluated against the expected max Sharpe of n_trials null strategies.

Source code in skepsis/core/psr.py

def deflated_sharpe_ratio(
    sr: float, n_obs: int, skew: float, kurt: float, var_trial_sr: float, n_trials: int
) -> float:
    """PSR evaluated against the expected max Sharpe of n_trials null strategies."""
    sr_star = expected_max_sharpe(var_trial_sr, n_trials)
    return probabilistic_sharpe_ratio(sr, sr_star, n_obs, skew, kurt)

expected_max_sharpe(var_trial_sr, n_trials)

E[max periodic SR] across n_trials strategies under the no-skill null.

E[max SR] ~= sqrt(V) * ((1-gamma) * z(1 - 1/N) + gamma * z(1 - 1/(N*e))) where gamma is the Euler-Mascheroni constant. Returns 0.0 for n_trials == 1 (a single trial has no selection bias; the formula diverges at N=1).

Source code in skepsis/core/psr.py

def expected_max_sharpe(var_trial_sr: float, n_trials: int) -> float:
    """E[max periodic SR] across n_trials strategies under the no-skill null.

    E[max SR] ~= sqrt(V) * ((1-gamma) * z(1 - 1/N) + gamma * z(1 - 1/(N*e)))
    where gamma is the Euler-Mascheroni constant. Returns 0.0 for n_trials == 1
    (a single trial has no selection bias; the formula diverges at N=1).
    """
    if var_trial_sr < 0.0:
        raise InvalidInputError(f"var_trial_sr must be >= 0, got {var_trial_sr}")
    if n_trials < 1:
        raise InvalidInputError(f"n_trials must be >= 1, got {n_trials}")
    if n_trials == 1:
        return 0.0
    sd = math.sqrt(var_trial_sr)
    return sd * (
        (1.0 - _EULER_GAMMA) * float(stats.norm.ppf(1.0 - 1.0 / n_trials))
        + _EULER_GAMMA * float(stats.norm.ppf(1.0 - 1.0 / (n_trials * math.e)))
    )

probabilistic_sharpe_ratio(sr, sr_benchmark, n_obs, skew, kurt)

P[true SR > sr_benchmark], correcting for sample length, skew and kurtosis.

PSR = Phi( (sr - sr) * sqrt(T - 1) / sqrt(1 - skewsr + (kurt - 1)/4 * sr^2) )

Source code in skepsis/core/psr.py

def probabilistic_sharpe_ratio(
    sr: float, sr_benchmark: float, n_obs: int, skew: float, kurt: float
) -> float:
    """P[true SR > sr_benchmark], correcting for sample length, skew and kurtosis.

    PSR = Phi( (sr - sr*) * sqrt(T - 1) / sqrt(1 - skew*sr + (kurt - 1)/4 * sr^2) )
    """
    if n_obs < _MIN_OBS:
        raise InsufficientDataError(f"PSR needs >= {_MIN_OBS} observations, got {n_obs}")
    den_sq = 1.0 - skew * sr + (kurt - 1.0) / 4.0 * sr * sr
    if den_sq <= 0.0:
        raise InvalidInputError(
            f"PSR denominator is non-positive ({den_sq:.6f}); the skew/kurtosis "
            "estimates are too extreme relative to the Sharpe ratio for the "
            "PSR approximation to hold"
        )
    z = (sr - sr_benchmark) * math.sqrt(n_obs - 1.0) / math.sqrt(den_sq)
    return float(stats.norm.cdf(z))

skepsis.core.pbo

Probability of Backtest Overfitting via CSCV.

Reference: Bailey, Borwein, Lopez de Prado & Zhu (2015), "The Probability of Backtest Overfitting", Journal of Computational Finance.

CSCV: split the (T, N) trial-returns matrix into S equal time blocks. For every combination of S/2 blocks used as in-sample (IS): rank trials IS, take the IS winner, find its relative rank omega among out-of-sample (OOS) scores, and compute the logit lambda = ln(omega / (1 - omega)). PBO is the fraction of combinations with lambda <= 0 (the IS winner lands in the bottom half OOS).

The default metric (periodic column Sharpe) runs on a fast path that combines per-block sums and sums-of-squares instead of materializing submatrices. Columns with zero variance score -inf so they rank last. A custom metric callback forces the materializing path.

PboResult dataclass

CSCV output. value is the PBO in [0, 1]; logits has one entry per combination.

Source code in skepsis/core/pbo.py

@dataclass(frozen=True)
class PboResult:
    """CSCV output. `value` is the PBO in [0, 1]; `logits` has one entry per combination."""

    value: float
    logits: np.ndarray
    n_combinations: int
    n_blocks: int
    n_trials: int
    n_obs_used: int

cscv(trials, n_blocks=16, metric=None)

Run CSCV on a (T, N) matrix of per-period trial returns (columns = trials).

Source code in skepsis/core/pbo.py

def cscv(
    trials: np.ndarray,
    n_blocks: int = 16,
    metric: Callable[[np.ndarray], np.ndarray] | None = None,
) -> PboResult:
    """Run CSCV on a (T, N) matrix of per-period trial returns (columns = trials)."""
    if trials.ndim != 2 or trials.shape[1] < 2:
        raise InvalidInputError("trials must be a (T, N) matrix with N >= 2 columns")
    validate_n_blocks(n_blocks)
    n_obs, n_trials = trials.shape
    if n_obs < 2 * n_blocks:
        raise InsufficientDataError(
            f"CSCV needs at least 2 * n_blocks = {2 * n_blocks} observations "
            f"(2 rows per block), got {n_obs}"
        )
    trimmed = n_obs % n_blocks
    if trimmed:
        warnings.warn(
            f"trimmed last {trimmed} observation(s) so {n_obs} rows split into "
            f"{n_blocks} equal blocks",
            SkepsisWarning,
            stacklevel=2,
        )
    n_used = n_obs - trimmed
    rows = n_used // n_blocks
    blocks = [trials[i * rows : (i + 1) * rows] for i in range(n_blocks)]
    combos = list(itertools.combinations(range(n_blocks), n_blocks // 2))
    logits = np.empty(len(combos))

    if metric is None:
        s1 = np.array([b.sum(axis=0) for b in blocks])  # (S, N) block sums
        s2 = np.array([(b * b).sum(axis=0) for b in blocks])  # (S, N) block sums of squares
        all_blocks = frozenset(range(n_blocks))

        def perf(idx: frozenset[int]) -> np.ndarray:
            n = rows * len(idx)
            sel = list(idx)
            t1 = s1[sel].sum(axis=0)
            t2 = s2[sel].sum(axis=0)
            var = (t2 - t1 * t1 / n) / (n - 1)
            # A genuinely constant column should give var == 0 exactly, but
            # var is computed here from accumulated block sums/sums-of-squares
            # (t2 - t1^2/n), and that catastrophic-cancellation path can leave
            # a tiny float64 residual (~1e-12 relative) instead of an exact
            # zero -- `var > 0` would then treat a degenerate column as real.
            # Compare against a margin scaled by t2/n (>= 0, the column's
            # mean-square magnitude) instead: genuine return series have
            # var / (mean square) >> 1e-12, so only truly degenerate columns
            # (including the all-zero column, where t2/n == 0) are caught.
            degenerate = var <= (t2 / n) * 1e-12
            with np.errstate(divide="ignore", invalid="ignore"):
                out: np.ndarray = np.where(degenerate, -np.inf, (t1 / n) / np.sqrt(var))
            return out

        for i, combo in enumerate(combos):
            is_set = frozenset(combo)
            n_star = int(np.argmax(perf(is_set)))
            logits[i] = _rank_logit(perf(all_blocks - is_set), n_star, n_trials)
    else:
        for i, combo in enumerate(combos):
            is_set = frozenset(combo)
            m_is = np.vstack([blocks[j] for j in sorted(is_set)])
            m_oos = np.vstack([blocks[j] for j in range(n_blocks) if j not in is_set])
            n_star = int(np.argmax(metric(m_is)))
            logits[i] = _rank_logit(metric(m_oos), n_star, n_trials)

    return PboResult(
        value=float(np.mean(logits <= 0.0)),
        logits=logits,
        n_combinations=len(combos),
        n_blocks=n_blocks,
        n_trials=n_trials,
        n_obs_used=n_used,
    )

validate_n_blocks(n_blocks)

Raise InvalidInputError unless n_blocks is an integer, even, and within [4, 24].

Source code in skepsis/core/pbo.py

def validate_n_blocks(n_blocks: int) -> None:
    """Raise InvalidInputError unless n_blocks is an integer, even, and within [4, 24]."""
    if isinstance(n_blocks, bool):
        raise InvalidInputError(f"n_blocks must be an integer, got bool {n_blocks!r}")
    try:
        operator.index(n_blocks)
    except TypeError:
        raise InvalidInputError(f"n_blocks must be an integer, got {n_blocks!r}") from None
    if n_blocks % 2 != 0:
        raise InvalidInputError(f"n_blocks must be even, got {n_blocks}")
    if not _MIN_BLOCKS <= n_blocks <= _MAX_BLOCKS:
        raise InvalidInputError(
            f"n_blocks must be between {_MIN_BLOCKS} and {_MAX_BLOCKS}, got {n_blocks}"
        )

skepsis.core.bootstrap

Stationary block bootstrap for Sharpe and drawdown confidence intervals.

References: - Politis & Romano (1994), "The Stationary Bootstrap", JASA 89(428). - Politis & White (2004), "Automatic Block-Length Selection for the Dependent Bootstrap", Econometric Reviews 23(1), with the Patton-Politis-White (2009) correction.

The no-skill p-value resamples the DEMEANED series with the same index matrix and reports (1 + #{SR_null >= SR_obs}) / (n_resamples + 1). Degenerate (exactly-constant) resamples score a signed-infinite Sharpe -- except an exactly-zero-mean constant draw, which scores exactly 0.0 and so correctly TIES a zero observed Sharpe -- and count toward that exceedance total; they are excluded from the CI/report distributions via the exact constancy mask (not finiteness, since a zero-mean constant draw is finite). A bootstrap whose resamples are ALL degenerate raises rather than attempting to form a CI. See BootstrapResult for the exact convention.

BootstrapResult dataclass

Bootstrap distributions and the no-skill p-value. Sharpe values are annualized.

Degenerate (exactly constant) resamples score a signed-infinite Sharpe (sign(mean) * inf), EXCEPT an exactly-zero-mean constant resample (zero return, zero risk), which scores exactly 0.0 so it correctly TIES a zero observed Sharpe (0.0 >= 0.0 is True) instead of vanishing as sign(0) * inf == nan. Those rows count toward p_value_no_skill on the null side (+inf and an exact-tying 0.0 are exceedances; -inf is not), keeping the (1 + #exceedances) / (n_resamples + 1) denominator exact. They are excluded from sharpe_ci and sharpe_distribution via the exact constancy mask -- not finiteness, since a zero-mean constant row is finite (0.0) but still degenerate; n_degenerate_resamples reports how many raw-side resamples were degenerate. If every raw-side resample is degenerate, bootstrap() raises InvalidInputError instead of returning a result with no CI.

Source code in skepsis/core/bootstrap.py

@dataclass(frozen=True)
class BootstrapResult:
    """Bootstrap distributions and the no-skill p-value. Sharpe values are annualized.

    Degenerate (exactly constant) resamples score a signed-infinite Sharpe
    (`sign(mean) * inf`), EXCEPT an exactly-zero-mean constant resample
    (zero return, zero risk), which scores exactly `0.0` so it correctly
    TIES a zero observed Sharpe (`0.0 >= 0.0` is True) instead of vanishing
    as `sign(0) * inf == nan`. Those rows count toward `p_value_no_skill` on
    the null side (`+inf` and an exact-tying `0.0` are exceedances; `-inf`
    is not), keeping the `(1 + #exceedances) / (n_resamples + 1)`
    denominator exact. They are excluded from `sharpe_ci` and
    `sharpe_distribution` via the exact constancy mask -- not finiteness,
    since a zero-mean constant row is finite (`0.0`) but still degenerate;
    `n_degenerate_resamples` reports how many raw-side resamples were
    degenerate. If every raw-side resample is degenerate, `bootstrap()`
    raises `InvalidInputError` instead of returning a result with no CI.
    """

    sharpe_obs: float
    sharpe_ci: tuple[float, float]
    drawdown_obs: float
    drawdown_ci: tuple[float, float]
    p_value_no_skill: float
    mean_block_length: float
    n_resamples: int
    n_degenerate_resamples: int
    sharpe_distribution: np.ndarray
    drawdown_distribution: np.ndarray

bootstrap(returns, periods, n_resamples=5000, mean_block_length=None, seed=0, ci=0.95)

Stationary-bootstrap CIs for annualized Sharpe and max drawdown, plus a p-value against the no-skill (demeaned) null.

Source code in skepsis/core/bootstrap.py

def bootstrap(
    returns: np.ndarray,
    periods: float,
    n_resamples: int = 5000,
    mean_block_length: float | None = None,
    seed: int = 0,
    ci: float = 0.95,
) -> BootstrapResult:
    """Stationary-bootstrap CIs for annualized Sharpe and max drawdown, plus
    a p-value against the no-skill (demeaned) null."""
    if not math.isfinite(periods) or periods <= 0:
        raise InvalidInputError(f"periods must be finite and > 0, got {periods}")
    if len(returns) < MIN_OBS:
        raise InsufficientDataError(
            f"bootstrap needs >= {MIN_OBS} observations, got {len(returns)}"
        )
    if not 0.5 < ci < 1.0:
        raise InvalidInputError(f"ci must be in (0.5, 1), got {ci}")
    if mean_block_length is None:
        mean_block_length = politis_white_block_length(returns)
    n_obs = len(returns)
    _validate_index_knobs(n_obs, mean_block_length, n_resamples)
    # Fail fast on exactly-constant returns (before spending time resampling):
    # sharpe() rejects them, which is also the observed-Sharpe value we need.
    sharpe_obs = annualized_sharpe(returns, periods)
    rng = np.random.default_rng(seed)
    centered = returns - returns.mean()

    sharpe_chunks: list[np.ndarray] = []
    raw_mask_chunks: list[np.ndarray] = []
    dd_chunks: list[np.ndarray] = []
    null_chunks: list[np.ndarray] = []
    for idx_chunk in _iter_index_chunks(n_obs, mean_block_length, n_resamples, rng):
        raw_values, raw_mask = _sharpe_rows(returns[idx_chunk], periods)
        sharpe_chunks.append(raw_values)
        raw_mask_chunks.append(raw_mask)
        dd_chunks.append(_drawdown_rows(returns[idx_chunk]))
        null_values, _null_mask = _sharpe_rows(centered[idx_chunk], periods)
        null_chunks.append(null_values)

    sharpe_raw = np.concatenate(sharpe_chunks)
    raw_mask = np.concatenate(raw_mask_chunks)
    dd_dist = np.concatenate(dd_chunks)
    null_dist = np.concatenate(null_chunks)

    # Degenerate (exactly-constant) resamples score a signed-infinite Sharpe,
    # except an exactly-zero-mean constant draw which scores exactly 0.0 (see
    # _sharpe_rows). Degeneracy accounting and CI exclusion use the exact
    # constancy MASK, not finiteness -- a zero-mean constant row is finite
    # (0.0) but still degenerate and must not populate the CI/report-facing
    # sharpe_distribution. They remain in null_dist for the p-value
    # exceedance count below (unchanged formula: +inf and an exact 0.0 tying
    # a zero observed Sharpe both count; -inf never does), so the
    # (1 + #exceedances) / (n_resamples + 1) denominator stays exact.
    n_degenerate = int(raw_mask.sum())
    if raw_mask.all():
        raise InvalidInputError(
            f"all {n_resamples} resamples were constant (zero variance); the series "
            "is too close to constant to bootstrap confidence intervals"
        )
    if n_degenerate > 0:
        warnings.warn(
            f"{n_degenerate} of {n_resamples} resamples were constant (zero variance); "
            "they count toward the p-value as signed-infinite Sharpe (or 0.0 for an "
            "exactly-zero-mean constant draw) and are excluded from CI quantiles",
            SkepsisWarning,
            stacklevel=2,
        )
    sharpe_dist = sharpe_raw[~raw_mask]

    p_value = float(1 + np.sum(null_dist >= sharpe_obs)) / (n_resamples + 1)

    lo, hi = (1.0 - ci) / 2.0, 1.0 - (1.0 - ci) / 2.0
    dd_finite = dd_dist[np.isfinite(dd_dist)]  # always finite; belt-and-suspenders
    return BootstrapResult(
        sharpe_obs=sharpe_obs,
        sharpe_ci=(float(np.quantile(sharpe_dist, lo)), float(np.quantile(sharpe_dist, hi))),
        drawdown_obs=max_drawdown(returns),
        drawdown_ci=(float(np.quantile(dd_finite, lo)), float(np.quantile(dd_finite, hi))),
        p_value_no_skill=p_value,
        mean_block_length=float(mean_block_length),
        n_resamples=n_resamples,
        n_degenerate_resamples=n_degenerate,
        sharpe_distribution=sharpe_dist,
        drawdown_distribution=dd_dist,
    )

politis_white_block_length(x)

Automatic mean block length for the stationary bootstrap (Politis-White 2004).

Source code in skepsis/core/bootstrap.py

def politis_white_block_length(x: np.ndarray) -> float:
    """Automatic mean block length for the stationary bootstrap (Politis-White 2004)."""
    n = len(x)
    if n < MIN_OBS:
        raise InsufficientDataError(f"block-length selection needs >= {MIN_OBS} obs, got {n}")
    kn = max(5, int(math.sqrt(math.log10(n))))
    m_max = int(math.ceil(math.sqrt(n))) + kn
    b_max = math.ceil(min(3.0 * math.sqrt(n), n / 3.0))
    xc = x - x.mean()
    acov = np.array([float(np.dot(xc[: n - k], xc[k:])) / n for k in range(m_max + 1)])
    if acov[0] == 0.0:
        raise InvalidInputError("returns have zero variance; block length is undefined")
    acorr = acov / acov[0]
    band = 2.0 * math.sqrt(math.log10(n) / n)
    m_hat = None
    for m in range(1, m_max + 1):
        window = acorr[m + 1 : m + kn + 1]
        if len(window) < kn:
            break
        if bool(np.all(np.abs(window) < band)):
            m_hat = m
            break
    if m_hat is None:
        big = np.where(np.abs(acorr[1:]) > band)[0]
        m_hat = int(big.max()) + 1 if len(big) else 1
    m_top = min(2 * m_hat, m_max)
    lags = np.arange(1, m_top + 1)
    lam = np.where(lags / m_top <= 0.5, 1.0, 2.0 * (1.0 - lags / m_top))
    g_hat = 2.0 * float(np.sum(lam * lags * acov[1 : m_top + 1]))
    d_hat = 2.0 * float(acov[0] + 2.0 * np.sum(lam * acov[1 : m_top + 1])) ** 2
    if d_hat <= 0.0 or g_hat == 0.0:
        return 1.0
    b = ((2.0 * g_hat * g_hat) / d_hat) ** (1.0 / 3.0) * n ** (1.0 / 3.0)
    return float(np.clip(b, 1.0, b_max))

stationary_bootstrap_indices(n_obs, mean_block_length, n_resamples, rng)

(n_resamples, n_obs) index matrix: geometric block lengths, circular wrap.

Stationary-bootstrap formulation: each position independently starts a new block with probability p = 1/mean_block_length (this IS the stationary bootstrap of Politis & Romano — geometric block lengths emerge from the per-position Bernoulli trials); block starts are uniform on [0, n_obs); within a block, indices continue circularly from the block's start.

Generated in row chunks of _CHUNK_ROWS to bound peak memory instead of materializing several full (n_resamples, n_obs) temporaries at once; the resulting seeded stream is otherwise identical in distribution but differs bit-for-bit from earlier, fully-vectorized development builds.

Source code in skepsis/core/bootstrap.py

def stationary_bootstrap_indices(
    n_obs: int, mean_block_length: float, n_resamples: int, rng: np.random.Generator
) -> np.ndarray:
    """(n_resamples, n_obs) index matrix: geometric block lengths, circular wrap.

    Stationary-bootstrap formulation: each position independently starts a new
    block with probability p = 1/mean_block_length (this IS the stationary
    bootstrap of Politis & Romano — geometric block lengths emerge from the
    per-position Bernoulli trials); block starts are uniform on [0, n_obs);
    within a block, indices continue circularly from the block's start.

    Generated in row chunks of `_CHUNK_ROWS` to bound peak memory instead of
    materializing several full (n_resamples, n_obs) temporaries at once; the
    resulting seeded stream is otherwise identical in distribution but differs
    bit-for-bit from earlier, fully-vectorized development builds.
    """
    _validate_index_knobs(n_obs, mean_block_length, n_resamples)
    out = np.empty((n_resamples, n_obs), dtype=np.int64)
    row = 0
    for chunk in _iter_index_chunks(n_obs, mean_block_length, n_resamples, rng):
        out[row : row + chunk.shape[0]] = chunk
        row += chunk.shape[0]
    return out

skepsis.core.sensitivity

Parameter-neighborhood stability: is the chosen configuration a plateau or a spike?

Method: k = min(2*d, n-1) nearest neighbors of the chosen configuration in z-scored parameter space (for interior points of a regular grid these are the orthogonal grid neighbors). stability_score = chosen_metric / median(neighbor metrics): ~1.0 on a plateau, >> 1.0 on an isolated spike (fitted to noise).

Edge conventions: - chosen metric <= 0 -> score = nan, SkepsisWarning (a non-positive chosen metric is its own problem; the ratio is meaningless) - neighbor median <= 0 < chosen -> score = inf, flagged, SkepsisWarning

SensitivityResult dataclass

Neighborhood stability of the chosen parameter configuration.

Source code in skepsis/core/sensitivity.py

@dataclass(frozen=True)
class SensitivityResult:
    """Neighborhood stability of the chosen parameter configuration."""

    stability_score: float
    chosen_index: int
    k: int
    neighbor_indices: np.ndarray
    neighbor_median: float
    flagged: bool
    method: str = "knn"

sensitivity(params, metrics, chosen_index, spike_threshold=1.5)

Score the chosen configuration against its k nearest parameter neighbors.

Source code in skepsis/core/sensitivity.py

def sensitivity(
    params: np.ndarray,
    metrics: np.ndarray,
    chosen_index: int,
    spike_threshold: float = 1.5,
) -> SensitivityResult:
    """Score the chosen configuration against its k nearest parameter neighbors."""
    p = np.asarray(params, dtype=np.float64)
    if p.ndim == 1:
        p = p[:, None]
    m = np.asarray(metrics, dtype=np.float64)
    n, d = p.shape
    if n < _MIN_TRIALS:
        raise InsufficientDataError(f"sensitivity needs >= {_MIN_TRIALS} trials, got {n}")
    if m.shape != (n,):
        raise InvalidInputError(f"metrics length {m.shape} must match params rows ({n},)")
    if not 0 <= chosen_index < n:
        raise InvalidInputError(f"chosen_index {chosen_index} out of range [0, {n})")

    sd = p.std(axis=0)
    sd[sd == 0.0] = 1.0  # constant parameter column: distance contribution 0
    z = (p - p.mean(axis=0)) / sd
    dists = np.linalg.norm(z - z[chosen_index], axis=1)
    order = np.argsort(dists)
    k = min(2 * d, n - 1)
    neighbors = np.array([i for i in order if i != chosen_index][:k])
    neighbor_median = float(np.median(m[neighbors]))
    chosen = float(m[chosen_index])

    if chosen <= 0.0:
        warnings.warn(
            f"chosen configuration has non-positive metric ({chosen:.4f}); "
            "stability score is undefined",
            SkepsisWarning,
            stacklevel=2,
        )
        score, flagged = float("nan"), False
    elif neighbor_median <= 0.0:
        warnings.warn(
            f"neighbors have non-positive median metric ({neighbor_median:.4f}) while the "
            "chosen configuration is positive — an extreme isolated spike",
            SkepsisWarning,
            stacklevel=2,
        )
        score, flagged = float("inf"), True
    else:
        score = chosen / neighbor_median
        flagged = score > spike_threshold

    return SensitivityResult(
        stability_score=score,
        chosen_index=int(chosen_index),
        k=int(k),
        neighbor_indices=neighbors,
        neighbor_median=neighbor_median,
        flagged=flagged,
    )