Coverage for src/jsharpe/sharpe/corrections.py: 100%

1"""Multiple-testing corrections and false-discovery-rate control. 

2 

3This module groups the family-wise error rate corrections (Bonferroni, 

4Šidák, Holm) and the false-discovery-rate routines (critical value, FDR 

5control, predictive pFDR and observed oFDR) used when screening many 

6candidate strategies. 

7""" 

8# ruff: noqa: N802, N803, N806, TRY003 

9 

10import math 

11import warnings 

12 

13import numpy as np 

14import scipy 

15 

16from .psr import probabilistic_sharpe_ratio, sharpe_ratio_variance 

17 

18 

19def _valid_fdr_inputs(q: float, SR0: float, SR1: float, sigma0: float, sigma1: float, p_H1: float) -> bool: 

20 """Return whether the FDR-critical-value inputs are within their valid ranges. 

21 

22 Args: 

23 q: Desired False Discovery Rate; must be in (0, 1). 

24 SR0: Mean under the null hypothesis; must be < SR1. 

25 SR1: Mean under the alternative hypothesis. 

26 sigma0: Standard deviation under the null hypothesis; must be > 0. 

27 sigma1: Standard deviation under the alternative hypothesis; must be > 0. 

28 p_H1: Prior probability of the alternative hypothesis; must be in (0, 1). 

29 

30 Returns: 

31 ``True`` if every parameter is within its valid range. 

32 """ 

33 return SR0 < SR1 and 0 < q < 1 and 0 < p_H1 < 1 and sigma0 > 0 and sigma1 > 0 

34 

35 

36def _fdr_posterior(c: float, SR0: float, SR1: float, sigma0: float, sigma1: float, p_H1: float) -> float: 

37 """Posterior false-discovery probability P[H=0 | X > c] under the two-normal mixture. 

38 

39 Args: 

40 c: Candidate critical value. 

41 SR0: Mean under the null hypothesis. 

42 SR1: Mean under the alternative hypothesis. 

43 sigma0: Standard deviation under the null hypothesis. 

44 sigma1: Standard deviation under the alternative hypothesis. 

45 p_H1: Prior probability of the alternative hypothesis. 

46 

47 Returns: 

48 Posterior false discovery probability at threshold ``c`` (0 where non-finite). 

49 """ 

50 a = 1 / (1 + scipy.stats.norm.sf((c - SR1) / sigma1) / scipy.stats.norm.sf((c - SR0) / sigma0) * p_H1 / (1 - p_H1)) 

51 return float(np.where(np.isfinite(a), a, 0)) 

52 

53 

54def adjusted_p_values_bonferroni(ps: np.ndarray) -> np.ndarray: 

55 """Adjust p-values using Bonferroni correction for FWER control. 

56 

57 Multiplies each p-value by the number of tests M, capping at 1. 

58 This is the most conservative multiple testing correction. 

59 

60 Args: 

61 ps: Array of unadjusted p-values. 

62 

63 Returns: 

64 Array of Bonferroni-adjusted p-values, each in [0, 1]. 

65 

66 Example: 

67 >>> p_vals = np.array([0.01, 0.03, 0.05]) 

68 >>> adj_p = adjusted_p_values_bonferroni(p_vals) 

69 >>> np.allclose(adj_p, [0.03, 0.09, 0.15]) 

70 True 

71 """ 

72 M = len(ps) 

73 result: np.ndarray = np.minimum(1, M * ps) 

74 return result 

75 

76 

77def adjusted_p_values_sidak(ps: np.ndarray) -> np.ndarray: 

78 """Adjust p-values using Šidák correction for FWER control. 

79 

80 Uses the formula: 1 - (1 - p)^M, which is slightly less conservative 

81 than Bonferroni when tests are independent. 

82 

83 Args: 

84 ps: Array of unadjusted p-values. 

85 

86 Returns: 

87 Array of Šidák-adjusted p-values, each in [0, 1]. 

88 

89 Example: 

90 >>> p_vals = np.array([0.01, 0.03, 0.05]) 

91 >>> adj_p = adjusted_p_values_sidak(p_vals) 

92 >>> np.all(adj_p <= adjusted_p_values_bonferroni(p_vals)) 

93 np.True_ 

94 """ 

95 M = len(ps) 

96 return 1 - (1 - ps) ** M 

97 

98 

99def adjusted_p_values_holm(ps: np.ndarray, *, variant: str = "bonferroni") -> np.ndarray: 

100 """Adjust p-values using Holm's step-down procedure for FWER control. 

101 

102 A step-down procedure that is uniformly more powerful than Bonferroni 

103 while still controlling the Family-Wise Error Rate. 

104 

105 Args: 

106 ps: Array of unadjusted p-values. 

107 variant: Correction method for each step. Either "bonferroni" 

108 (default) or "sidak". 

109 

110 Returns: 

111 Array of Holm-adjusted p-values, each in [0, 1]. 

112 

113 Raises: 

114 ValueError: If variant is not "bonferroni" or "sidak". 

115 

116 Example: 

117 >>> p_vals = np.array([0.01, 0.04, 0.03]) 

118 >>> adj_p = adjusted_p_values_holm(p_vals) 

119 >>> float(adj_p[0]) # Smallest p-value adjusted most 

120 0.03 

121 """ 

122 if variant not in ("bonferroni", "sidak"): 

123 raise ValueError(f"Unknown Holm variant {variant!r}; expected 'bonferroni' or 'sidak'.") 

124 i = np.argsort(ps) 

125 M = len(ps) 

126 p_adjusted = np.zeros(M) 

127 previous = 0 

128 for j, idx in enumerate(i): 

129 candidate = min(1, ps[idx] * (M - j)) if variant == "bonferroni" else 1 - (1 - ps[idx]) ** (M - j) 

130 p_adjusted[idx] = max(previous, candidate) 

131 previous = p_adjusted[idx] 

132 return p_adjusted 

133 

134 

135def FDR_critical_value(q: float, SR0: float, SR1: float, sigma0: float, sigma1: float, p_H1: float) -> float: 

136 """Compute critical value for FDR control in hypothesis testing. 

137 

138 Given a mixture model where H ~ Bernoulli(p_H1) determines whether 

139 X follows N(SR0, sigma0^2) or N(SR1, sigma1^2), finds the critical 

140 value c such that P[H=0 | X > c] = q. 

141 

142 Args: 

143 q: Desired False Discovery Rate in (0, 1). 

144 SR0: Mean under null hypothesis (must be < SR1). 

145 SR1: Mean under alternative hypothesis. 

146 sigma0: Standard deviation under null hypothesis (must be > 0). 

147 sigma1: Standard deviation under alternative hypothesis (must be > 0). 

148 p_H1: Prior probability of alternative hypothesis in (0, 1). 

149 

150 Returns: 

151 Critical value c. Returns -inf if solution is outside [-10, 10], 

152 or nan if no solution exists. 

153 

154 Raises: 

155 ValueError: If parameters are out of valid ranges. 

156 

157 Example: 

158 >>> c = FDR_critical_value(q=0.2, SR0=0, SR1=0.5, sigma0=0.2, sigma1=0.3, p_H1=0.1) 

159 >>> c > 0 # Critical value should be positive 

160 True 

161 """ 

162 if not _valid_fdr_inputs(q, SR0, SR1, sigma0, sigma1, p_H1): 

163 raise ValueError( 

164 "Invalid FDR inputs: require SR0 < SR1, 0 < q < 1, 0 < p_H1 < 1, sigma0 > 0, and sigma1 > 0; " 

165 f"got q={q}, SR0={SR0}, SR1={SR1}, sigma0={sigma0}, sigma1={sigma1}, p_H1={p_H1}." 

166 ) 

167 

168 with warnings.catch_warnings(): 

169 warnings.filterwarnings("ignore", message="invalid value encountered in scalar divide") 

170 warnings.filterwarnings("ignore", message="divide by zero encountered in scalar divide") 

171 

172 f_lo = _fdr_posterior(-10, SR0, SR1, sigma0, sigma1, p_H1) 

173 f_hi = _fdr_posterior(10, SR0, SR1, sigma0, sigma1, p_H1) 

174 

175 if f_lo < q: # Solution outside of the search interval 

176 return float(-np.inf) 

177 

178 if (f_lo - q) * (f_hi - q) > 0: # No solution, for instance if σ₀≫σ₁ and q small 

179 return float(np.nan) 

180 

181 return float(scipy.optimize.brentq(lambda c: _fdr_posterior(c, SR0, SR1, sigma0, sigma1, p_H1) - q, -10, 10)) 

182 

183 

184def control_for_FDR( 

185 q: float, 

186 *, 

187 SR0: float = 0, 

188 SR1: float = 0.5, 

189 p_H1: float = 0.05, 

190 T: int = 24, 

191 gamma3: float = 0.0, 

192 gamma4: float = 3.0, 

193 rho: float = 0.0, 

194 K: int = 1, 

195) -> tuple[float, float, float, float]: 

196 """Compute critical value to test multiple Sharpe ratios controlling FDR. 

197 

198 Determines the critical Sharpe ratio threshold and associated error rates 

199 to control the False Discovery Rate at level q when testing multiple 

200 strategies. 

201 

202 Args: 

203 q: Desired False Discovery Rate level in (0, 1). 

204 SR0: Sharpe ratio under null hypothesis H0. Default 0. 

205 SR1: Sharpe ratio under alternative hypothesis H1. Default 0.5. 

206 p_H1: Prior probability that H1 is true. Default 0.05. 

207 T: Number of observations. Default 24. 

208 gamma3: Skewness of returns. Default 0. 

209 gamma4: Kurtosis of returns (non-excess). Default 3 (Gaussian). 

210 rho: Autocorrelation of returns. Default 0. 

211 K: Number of strategies (K=1 for FDR; K>1 for FWER-FDR). Default 1. 

212 

213 Returns: 

214 Tuple of (alpha, beta, SR_c, q_hat): 

215 - alpha: Significance level P[SR > SR_c | H0] 

216 - beta: Type II error P[SR <= SR_c | H1]; power is 1 - beta 

217 - SR_c: Critical Sharpe ratio threshold 

218 - q_hat: Estimated FDR (should be close to q) 

219 

220 Example: 

221 >>> alpha, beta, SR_c, q_hat = control_for_FDR(q=0.25, T=24) 

222 >>> bool(0 < alpha < 1) 

223 True 

224 >>> bool(SR_c > 0) # Critical value is positive 

225 True 

226 """ 

227 Z = scipy.stats.norm.cdf 

228 

229 s0 = math.sqrt(sharpe_ratio_variance(SR0, T, gamma3=gamma3, gamma4=gamma4, rho=rho, K=K)) 

230 s1 = math.sqrt(sharpe_ratio_variance(SR1, T, gamma3=gamma3, gamma4=gamma4, rho=rho, K=K)) 

231 SRc = FDR_critical_value(q, SR0, SR1, s0, s1, p_H1) 

232 

233 beta = Z((SRc - SR1) / s1) 

234 alpha = q / (1 - q) * p_H1 / (1 - p_H1) * (1 - beta) 

235 q_hat = 1 / (1 + (1 - beta) / alpha * p_H1 / (1 - p_H1)) 

236 

237 return alpha, beta, SRc, q_hat 

238 

239 

240def pFDR( 

241 p_H1: float, 

242 alpha: float, 

243 beta: float, 

244) -> float: 

245 """Compute posterior FDR given test outcome exceeds critical value. 

246 

247 Calculates P[H0 | SR > SR_c], the probability that the null hypothesis 

248 is true given that the observed Sharpe ratio exceeds the critical value. 

249 This is the "predictive" FDR based on the critical value, not the 

250 observed value. 

251 

252 Args: 

253 p_H1: Prior probability that H1 is true. 

254 alpha: Significance level (Type I error rate). 

255 beta: Type II error rate (1 - power). 

256 

257 Returns: 

258 Posterior probability of H0 given rejection. 

259 

260 Example: 

261 >>> # With 5% prior on H1 and 5% significance 

262 >>> fdr = pFDR(p_H1=0.05, alpha=0.05, beta=0.3) 

263 >>> 0 < fdr < 1 

264 True 

265 """ 

266 p_H0 = 1 - p_H1 

267 return 1 / (1 + (1 - beta) * p_H1 / alpha / p_H0) 

268 

269 

270def oFDR( 

271 SR: float, 

272 SR0: float, 

273 SR1: float, 

274 T: int, 

275 p_H1: float, 

276 *, 

277 gamma3: float = 0.0, 

278 gamma4: float = 3.0, 

279 rho: float = 0.0, 

280 K: int = 1, 

281) -> float: 

282 """Compute observed FDR given the observed Sharpe ratio. 

283 

284 Calculates P[H0 | SR > SR_obs], the probability that the null hypothesis 

285 is true given the observed Sharpe ratio value. This is the "observed" 

286 FDR which conditions on the actual observation rather than just the 

287 critical value. 

288 

289 Args: 

290 SR: Observed Sharpe ratio. 

291 SR0: Sharpe ratio under null hypothesis. 

292 SR1: Sharpe ratio under alternative hypothesis. 

293 T: Number of observations. 

294 p_H1: Prior probability that H1 is true. 

295 gamma3: Skewness of returns. Default 0. 

296 gamma4: Kurtosis of returns (non-excess). Default 3 (Gaussian). 

297 rho: Autocorrelation of returns. Default 0. 

298 K: Number of strategies for variance adjustment. Default 1. 

299 

300 Returns: 

301 Posterior probability of H0 given the observed SR. 

302 

303 Example: 

304 >>> # Higher observed SR should give lower probability of H0 

305 >>> fdr_low = oFDR(SR=0.3, SR0=0, SR1=0.5, T=24, p_H1=0.1) 

306 >>> fdr_high = oFDR(SR=0.8, SR0=0, SR1=0.5, T=24, p_H1=0.1) 

307 >>> bool(fdr_high < fdr_low) 

308 True 

309 """ 

310 p0 = 1 - probabilistic_sharpe_ratio(SR, SR0, T=T, gamma3=gamma3, gamma4=gamma4, rho=rho, K=K) 

311 p1 = 1 - probabilistic_sharpe_ratio(SR, SR1, T=T, gamma3=gamma3, gamma4=gamma4, rho=rho, K=K) 

312 p_H0 = 1 - p_H1 

313 return p0 * p_H0 / (p0 * p_H0 + p1 * p_H1)