Statistics of the Sharpe ratio

"""M04 exercise: statistics of the Sharpe ratio (notes §4.5 and §4.6).

Implement every function, then run from the repo root:

    cd exercises && .venv/bin/python -m pytest modules/M04-systematic-research-swing/tests/test_m04_sharpe_stats.py

All Sharpe ratios here are PER PERIOD (not annualised) unless a function says
otherwise. If returns are daily, S is a daily Sharpe ratio and T is a number of
days. Kurtosis `kurt` is the ordinary (non-excess) kurtosis gamma_4 (normal = 3).
Use numpy and scipy.stats.norm only.
"""

from __future__ import annotations

import numpy as np
from scipy.stats import norm

EULER_GAMMA = 0.5772156649015329


def sharpe(returns, rf: float = 0.0) -> float:
    """Sample Sharpe ratio S_hat = (mean(r) - rf) / std(r, ddof=1), per period."""
    raise NotImplementedError


def sharpe_se(sr: float, T: int, skew: float = 0.0, kurt: float = 3.0) -> float:
    """Asymptotic standard error of the sample Sharpe ratio (notes eq. 4.12, Prop. 4.11).

    SE = sqrt( (1 - skew*S + (kurt - 1)/4 * S^2) / T ).
    With skew = 0, kurt = 3 this is Lo's (2002) iid-normal result sqrt((1 + S^2/2)/T).
    """
    raise NotImplementedError


def annualise_sharpe(sr: float, q: int, rho=None) -> float:
    """Annualise a per-period Sharpe ratio with q periods per year (notes eq. 4.14).

    rho=None: iid scaling sqrt(q) * S.
    rho = sequence rho_1, rho_2, ... of autocorrelations (lags >= q ignored, missing lags = 0):
    eta(q) = q / sqrt(q + 2 * sum_{k=1}^{q-1} (q-k) rho_k), return eta(q) * S.
    """
    raise NotImplementedError


def probabilistic_sharpe(sr: float, T: int, skew: float = 0.0, kurt: float = 3.0, sr_benchmark: float = 0.0) -> float:
    """PSR(S*) = Phi( (S_hat - S*) sqrt(T - 1) / sqrt(1 - skew S_hat + (kurt-1)/4 S_hat^2) ) (notes eq. 4.15).

    Reference value: PSR(0.1, 252, 0, 3) = 0.942987.
    """
    raise NotImplementedError


def min_track_record_length(sr: float, sr_benchmark: float = 0.0, skew: float = 0.0, kurt: float = 3.0, prob: float = 0.95) -> float:
    """Smallest T with PSR(S*) >= prob, holding S_hat fixed (notes eq. 4.16).

    MinTRL = 1 + (1 - skew S + (kurt-1)/4 S^2) * (Phi^{-1}(prob) / (S - S*))^2.
    Return float('inf') if sr <= sr_benchmark.
    """
    raise NotImplementedError


def expected_max_sharpe(n_trials: int, var_sr: float = 1.0, mean_sr: float = 0.0) -> float:
    """Gumbel approximation to E[max of N iid N(mean_sr, var_sr)] (notes Prop. 4.13, eq. 4.17).

    E[max] ~ mean + sqrt(V) * ((1 - g) * Phi^{-1}(1 - 1/N) + g * Phi^{-1}(1 - 1/(N e))),
    g = EULER_GAMMA. With V = 1: N=10 -> 1.574598, N=100 -> 2.530603, N=1000 -> 3.255122.
    Raise ValueError if n_trials < 2.
    """
    raise NotImplementedError


def exact_expected_max_normal(n_trials: int, n_grid: int = 20001) -> float:
    """Exact E[max of N iid N(0,1)] = int x N phi(x) Phi(x)^{N-1} dx, by the trapezoid rule on [-10, 10]."""
    raise NotImplementedError


def deflated_sharpe(sr: float, T: int, n_trials: int, var_sr_trials: float, skew: float = 0.0, kurt: float = 3.0) -> float:
    """Deflated Sharpe ratio (notes eq. 4.18): PSR with benchmark S0 = expected_max_sharpe(n_trials, var_sr_trials).

    var_sr_trials is the cross-trial variance of the per-period Sharpe estimates.
    """
    raise NotImplementedError