Coverage for src / jsharpe / sharpe.py: 99%

1"""Sharpe-related utilities for statistical analysis and hypothesis testing. 

2 

3This module provides comprehensive tools for Sharpe ratio analysis, including: 

4 - Variance estimation under non-Gaussian returns 

5 - Statistical significance testing 

6 - Multiple testing corrections (FDR, FWER) 

7 - Portfolio optimization utilities 

8""" 

9# ruff: noqa: N802, N803, N806, S101, TRY003 

10 

11import math 

12import warnings 

13from collections.abc import Callable 

14 

15import numpy as np 

16import scipy 

17 

18 

19def ppoints(n: int, a: float | None = None) -> np.ndarray: 

20 """Generate probability points for Q-Q plots and distribution fitting. 

21 

22 Creates n equidistant points in the interval (0, 1), suitable for 

23 probability plotting positions. The boundaries 0 and 1 are excluded. 

24 This function mirrors the behavior of R's `ppoints()` function. 

25 

26 Args: 

27 n: Number of probability points to generate. 

28 a: Plotting position parameter in [0, 1]. Defaults to 0.5 if n > 10, 

29 otherwise 3/8 (following R's convention). 

30 

31 Returns: 

32 Array of n equidistant probability points in (0, 1). 

33 

34 Raises: 

35 ValueError: If a is not in [0, 1]. 

36 

37 Example: 

38 >>> import numpy as np 

39 >>> pts = ppoints(5) 

40 >>> len(pts) 

41 5 

42 >>> np.all((pts > 0) & (pts < 1)) 

43 np.True_ 

44 >>> ppoints(5, a=0.5) 

45 array([0.1, 0.3, 0.5, 0.7, 0.9]) 

46 """ 

47 if a is None: 

48 a = 0.5 if n > 10 else 3 / 8 

49 if not (0 <= a <= 1): 

50 raise ValueError(f"the offset should be in [0,1], got {a}") 

51 return np.linspace(1 - a, n - a, n) / (n + 1 - 2 * a) 

52 

53 

54def robust_covariance_inverse(V: np.ndarray) -> np.ndarray: 

55 r"""Compute inverse of a constant-correlation covariance matrix. 

56 

57 Uses the Sherman-Morrison formula for efficient computation. 

58 Assumes the variance matrix has the form: 

59 $V = \rho \sigma \sigma' + (1-\rho) \text{diag}(\sigma^2)$ 

60 (constant correlations across all pairs). 

61 

62 Its inverse is computed as: 

63 $V^{-1} = A^{-1} - \dfrac{ A^{-1} \rho \sigma \sigma' A^{-1} } 

64 { 1 + \rho \sigma' A^{-1} \sigma }$ 

65 

66 Args: 

67 V: Covariance matrix with constant off-diagonal correlations. 

68 Shape (n, n). 

69 

70 Returns: 

71 Inverse of the covariance matrix. Shape (n, n). 

72 

73 Example: 

74 >>> import numpy as np 

75 >>> # Create a constant-correlation covariance matrix 

76 >>> rho = 0.5 

77 >>> sigma = np.array([1.0, 2.0, 1.5]) 

78 >>> C = rho * np.ones((3, 3)) 

79 >>> np.fill_diagonal(C, 1) 

80 >>> V = (C * sigma.reshape(-1, 1)).T * sigma.reshape(-1, 1) 

81 >>> V_inv = robust_covariance_inverse(V) 

82 >>> np.allclose(V @ V_inv, np.eye(3), atol=1e-10) 

83 True 

84 """ 

85 sigma = np.sqrt(np.diag(V)) 

86 C = (V.T / sigma).T / sigma 

87 rho = np.mean(C[np.triu_indices_from(C, 1)]) 

88 A = np.diag(1 / sigma**2) / (1 - rho) 

89 sigma = sigma.reshape(-1, 1) 

90 result: np.ndarray = A - (rho * A @ sigma @ sigma.T @ A) / (1 + rho * sigma.T @ A @ sigma) 

91 return result 

92 

93 

94def minimum_variance_weights_for_correlated_assets(V: np.ndarray) -> np.ndarray: 

95 """Compute weights of the minimum variance portfolio for correlated assets. 

96 

97 Computes the portfolio weights that minimize portfolio variance subject 

98 to the constraint that weights sum to 1. Assumes a constant-correlation 

99 covariance structure for efficient computation. 

100 

101 Args: 

102 V: Covariance matrix of asset returns. Shape (n, n). 

103 

104 Returns: 

105 Portfolio weights that minimize variance. Shape (n,). 

106 Weights sum to 1. 

107 

108 Example: 

109 >>> import numpy as np 

110 >>> # Create a simple covariance matrix 

111 >>> V = np.array([[0.04, 0.01], [0.01, 0.09]]) 

112 >>> w = minimum_variance_weights_for_correlated_assets(V) 

113 >>> np.isclose(w.sum(), 1.0) 

114 np.True_ 

115 """ 

116 ones = np.ones(shape=V.shape[0]) 

117 S = robust_covariance_inverse(V) 

118 w = S @ ones 

119 w = w / np.sum(w) 

120 return w 

121 

122 

123def effective_rank(C: np.ndarray) -> float: 

124 """Compute the effective rank of a positive semi-definite matrix. 

125 

126 The effective rank measures the "effective dimensionality" of a matrix 

127 by computing the exponential of the entropy of its normalized eigenvalues. 

128 This provides a continuous measure between 1 (perfectly correlated) and 

129 n (perfectly uncorrelated/identity matrix). 

130 

131 Algorithm: 

132 1. Compute eigenvalues (non-negative for PSD matrices) 

133 2. Discard zero eigenvalues 

134 3. Normalize to form a probability distribution 

135 4. Compute entropy H = -Σ p_i log(p_i) 

136 5. Return exp(H) 

137 

138 Args: 

139 C: Positive semi-definite matrix (e.g., correlation matrix). 

140 Shape (n, n). 

141 

142 Returns: 

143 Effective rank, a value in [1, n] where n is the matrix dimension. 

144 

145 References: 

146 Roy, O. and Vetterli, M. (2007). "The effective rank: a measure of 

147 effective dimensionality." EURASIP Journal on Advances in Signal 

148 Processing. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.177.2721 

149 

150 Example: 

151 >>> import numpy as np 

152 >>> # Identity matrix has effective rank equal to its dimension 

153 >>> abs(effective_rank(np.eye(3)) - 3.0) < 1e-10 

154 True 

155 >>> # Perfectly correlated matrix has effective rank 1 

156 >>> C = np.ones((3, 3)) 

157 >>> abs(effective_rank(C) - 1.0) < 1e-10 

158 True 

159 """ 

160 p = np.linalg.eigvalsh(C) 

161 p = p[p > 0] 

162 p = p / sum(p) 

163 H = np.sum(-p * np.log(p)) 

164 return math.exp(H) 

165 

166 

167def _silhouette_samples(X: np.ndarray, labels: np.ndarray) -> np.ndarray: 

168 """Compute silhouette coefficients for each sample. 

169 

170 For each sample i: 

171 a(i) = mean distance from i to all other samples in the same cluster 

172 b(i) = min over other clusters c of mean distance from i to samples in c 

173 s(i) = (b(i) - a(i)) / max(a(i), b(i)) 

174 

175 Args: 

176 X: Feature matrix of shape (n_samples, n_features). 

177 labels: Cluster labels of shape (n_samples,). 

178 

179 Returns: 

180 Silhouette coefficients of shape (n_samples,). 

181 """ 

182 n = X.shape[0] 

183 unique_labels = np.unique(labels) 

184 n_clusters = len(unique_labels) 

185 

186 # Pairwise Euclidean distances: dist[i,j] = ||X[i] - X[j]|| 

187 diff = X[:, np.newaxis, :] - X[np.newaxis, :, :] 

188 dist = np.sqrt(np.einsum("ijk,ijk->ij", diff, diff)) 

189 

190 # For each cluster c, compute mean distance from every sample to cluster c. 

191 # For samples IN cluster c, the self-distance (zero) is excluded. 

192 label_to_idx = {c: idx for idx, c in enumerate(unique_labels)} 

193 cluster_mean_dist = np.zeros((n, n_clusters)) 

194 for c_idx, c in enumerate(unique_labels): 

195 mask = labels == c 

196 count = int(mask.sum()) 

197 dist_to_cluster = dist[:, mask].sum(axis=1) 

198 cluster_mean_dist[:, c_idx] = dist_to_cluster / max(count, 1) 

199 # For members of cluster c, exclude self-distance (which is 0) 

200 in_cluster = np.where(mask)[0] 

201 if count > 1: 

202 cluster_mean_dist[in_cluster, c_idx] = dist_to_cluster[in_cluster] / (count - 1) 

203 else: 

204 cluster_mean_dist[in_cluster, c_idx] = 0.0 

205 

206 # a(i): mean distance to same cluster 

207 own_cluster_idx = np.array([label_to_idx[lbl] for lbl in labels]) 

208 a = cluster_mean_dist[np.arange(n), own_cluster_idx] 

209 

210 # b(i): min mean distance to any other cluster 

211 b_mat = cluster_mean_dist.copy() 

212 b_mat[np.arange(n), own_cluster_idx] = np.inf 

213 b = b_mat.min(axis=1) 

214 b = np.where(b == np.inf, 0.0, b) 

215 

216 max_ab = np.maximum(a, b) 

217 return np.where(max_ab > 0, (b - a) / max_ab, 0.0) 

218 

219 

220def number_of_clusters( 

221 C: np.ndarray, 

222 *, 

223 retries: int = 10, 

224 max_clusters: int = 100, 

225) -> tuple[int, dict[int, float], np.ndarray]: 

226 """Compute the optimal number of clusters from a correlation matrix. 

227 

228 Implements the algorithm from section 8.1 of Lopez de Prado (2018): 

229 1. Convert the correlation matrix into a distance matrix. 

230 2. Using the columns of the distance matrix as features, run the 

231 k-means algorithm for each k and compute the quality of the 

232 clustering. 

233 3. Return the clustering with the highest quality. 

234 

235 Quality is defined as the mean of the silhouette scores divided by 

236 their standard deviation. 

237 

238 Args: 

239 C: Correlation matrix. Must be square, symmetric, finite, with ones 

240 on the diagonal and all entries in [-1, 1]. 

241 retries: Number of times to run k-means for each k to reduce the 

242 impact of random initialisation. Default 10. 

243 max_clusters: Maximum number of clusters to evaluate. Capped at 

244 ``C.shape[0] - 1``. Default 100. 

245 

246 Returns: 

247 Tuple of (n_clusters, qualities, labels): 

248 - n_clusters: Optimal number of clusters. 

249 - qualities: Dict mapping k to its quality score. 

250 - labels: Cluster assignment for each observation (shape (n,)). 

251 

252 References: 

253 Lopez de Prado, M. (2018). "Detection of false investment strategies 

254 using unsupervised learning methods." SSRN 3167017. 

255 https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3167017 

256 

257 Example: 

258 >>> np.random.seed(42) 

259 >>> C, _, _ = get_random_correlation_matrix( 

260 ... number_of_trials=20, effective_number_of_trials=4 

261 ... ) 

262 >>> n, qualities, labels = number_of_clusters(C, retries=3, max_clusters=8) 

263 >>> 2 <= n <= 8 

264 True 

265 >>> labels.shape 

266 (20,) 

267 """ 

268 assert isinstance(C, np.ndarray) 

269 assert np.all(C >= -1) 

270 assert np.all(C <= 1) 

271 assert np.all(np.diag(C) == 1) 

272 assert np.all(np.isfinite(C)) 

273 

274 max_clusters = min(max_clusters, C.shape[0] - 1) 

275 

276 # Convert correlations to distances 

277 D = np.sqrt((1 - C) / 2) 

278 assert np.all(np.isfinite(D)) 

279 

280 qualities: dict[int, float] = {} 

281 best_labels: dict[int, np.ndarray] = {} 

282 for k in range(2, max_clusters + 1): 

283 qualities[k] = -np.inf 

284 for _ in range(retries): 

285 with warnings.catch_warnings(): 

286 warnings.simplefilter("ignore") 

287 _centroids, labels = scipy.cluster.vq.kmeans2(D, k, minit="points", iter=300) 

288 # Skip degenerate solutions with empty clusters 

289 if len(np.unique(labels)) < k: 

290 continue 

291 silhouette_vals = _silhouette_samples(D, labels) 

292 std = silhouette_vals.std() 

293 if std == 0: 

294 continue 

295 q = float(silhouette_vals.mean() / std) 

296 if q > qualities[k]: 

297 qualities[k] = q 

298 best_labels[k] = labels.copy() 

299 

300 # Select the best k among those for which a valid solution was found 

301 valid_k = {k: q for k, q in qualities.items() if k in best_labels} 

302 if not valid_k: 

303 raise RuntimeError("No valid clustering solution found; try increasing retries or reducing max_clusters.") 

304 best_k = max(valid_k, key=lambda x: valid_k[x]) 

305 return best_k, qualities, best_labels[best_k] 

306 

307 

308def moments_Mk(k: int, *, rho: float = 0) -> tuple[float, float, float]: 

309 """Compute moments of the maximum of k standard normal random variables. 

310 

311 Computes E[M_k], E[M_k^2], and Var[M_k] where M_k = max(Z_1, ..., Z_k) 

312 and Z_i are standard normal. Supports both independent (rho=0) and 

313 equi-correlated (rho > 0) cases. 

314 

315 For the correlated case with equi-correlation rho: 

316 Z_i = sqrt(rho) * X + sqrt(1-rho) * Y_i 

317 M = sqrt(rho) * X + sqrt(1-rho) * max(Y_i) 

318 

319 Args: 

320 k: Number of standard normal random variables. 

321 rho: Equi-correlation coefficient in [0, 1). Default is 0 (independent). 

322 

323 Returns: 

324 Tuple of (Ez, Ez2, var): 

325 - Ez: Expected value E[M_k] 

326 - Ez2: Second moment E[M_k^2] 

327 - var: Variance Var[M_k] 

328 

329 Example: 

330 >>> Ez, Ez2, var = moments_Mk(1) 

331 >>> bool(abs(Ez) < 1e-10) # E[max of one N(0,1)] = 0 

332 True 

333 >>> bool(abs(var - 1.0) < 1e-10) # Var[max of one N(0,1)] = 1 

334 True 

335 """ 

336 Phi = scipy.stats.norm.cdf 

337 Ez = E_under_normal(lambda z: k * z * (Phi(z) ** (k - 1))) 

338 Ez2 = E_under_normal(lambda z: k * z**2 * (Phi(z) ** (k - 1))) 

339 var = Ez2 - Ez**2 

340 

341 Ez = (1 - rho) * Ez 

342 var = rho + (1 - rho) * var 

343 Ez2 = var + Ez**2 

344 

345 return Ez, Ez2, var 

346 

347 

348def sharpe_ratio_variance( 

349 SR: float, 

350 T: int, 

351 *, 

352 gamma3: float = 0.0, 

353 gamma4: float = 3.0, 

354 rho: float = 0.0, 

355 K: int = 1, 

356) -> float: 

357 """Compute the asymptotic variance of the Sharpe ratio estimator. 

358 

359 Accounts for non-Gaussian returns (skewness, kurtosis) and autocorrelation. 

360 Under Gaussian iid returns, reduces to (1 + SR^2/2) / T. 

361 

362 Args: 

363 SR: Sharpe ratio value (annualized or per-period). 

364 T: Number of observations (time periods). 

365 gamma3: Skewness of returns. Default 0 (symmetric). 

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

367 rho: First-order autocorrelation of returns. Default 0 (iid). 

368 K: Number of strategies for multiple testing adjustment. Default 1. 

369 

370 Returns: 

371 Variance of the Sharpe ratio estimator. 

372 

373 Example: 

374 >>> # Variance under Gaussian iid assumptions 

375 >>> var_gaussian = sharpe_ratio_variance(SR=0.5, T=24) 

376 >>> # Variance with non-Gaussian returns (higher kurtosis) 

377 >>> var_nongauss = sharpe_ratio_variance(SR=0.5, T=24, gamma4=6.0) 

378 >>> bool(var_nongauss > var_gaussian) # Higher kurtosis increases variance 

379 True 

380 """ 

381 A = 1 

382 B = rho / (1 - rho) 

383 C = rho**2 / (1 - rho**2) 

384 a = A + 2 * B 

385 b = A + B + C 

386 c = A + 2 * C 

387 V = (a * 1 - b * gamma3 * SR + c * (gamma4 - 1) / 4 * SR**2) / T 

388 return float(V * moments_Mk(K)[2]) 

389 

390 

391def variance_of_the_maximum_of_k_Sharpe_ratios(number_of_trials: int, variance: float) -> float: 

392 """Compute the variance of the maximum Sharpe ratio across K strategies. 

393 

394 Selection across a larger pool increases the uncertainty of the selected 

395 (maximum) Sharpe ratio estimate. This function applies a logarithmic 

396 variance inflation factor that increases monotonically with K. 

397 

398 Args: 

399 number_of_trials: Number of strategies (K >= 1). 

400 variance: Base variance of individual Sharpe ratio estimates. 

401 

402 Returns: 

403 Inflated variance accounting for selection from K strategies. 

404 

405 Example: 

406 >>> # More trials increases variance due to selection bias 

407 >>> v1 = variance_of_the_maximum_of_k_Sharpe_ratios(1, 0.1) 

408 >>> v10 = variance_of_the_maximum_of_k_Sharpe_ratios(10, 0.1) 

409 >>> bool(v10 > v1) 

410 True 

411 """ 

412 # Monotone increasing variance inflation with K (K>=1) 

413 inflation = 1.0 + np.log(max(1, int(number_of_trials))) 

414 return float(variance * inflation) 

415 

416 

417def control_for_FDR( 

418 q: float, 

419 *, 

420 SR0: float = 0, 

421 SR1: float = 0.5, 

422 p_H1: float = 0.05, 

423 T: int = 24, 

424 gamma3: float = 0.0, 

425 gamma4: float = 3.0, 

426 rho: float = 0.0, 

427 K: int = 1, 

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

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

430 

431 Determines the critical Sharpe ratio threshold and associated error rates 

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

433 strategies. 

434 

435 Args: 

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

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

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

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

440 T: Number of observations. Default 24. 

441 gamma3: Skewness of returns. Default 0. 

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

443 rho: Autocorrelation of returns. Default 0. 

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

445 

446 Returns: 

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

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

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

450 - SR_c: Critical Sharpe ratio threshold 

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

452 

453 Example: 

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

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

456 True 

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

458 True 

459 """ 

460 Z = scipy.stats.norm.cdf 

461 

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

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

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

465 

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

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

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

469 

470 return alpha, beta, SRc, q_hat 

471 

472 

473def expected_maximum_sharpe_ratio(number_of_trials: int, variance: float, SR0: float = 0) -> float: 

474 """Compute expected maximum Sharpe ratio under multiple testing. 

475 

476 Estimates E[max(SR_1, ..., SR_K)] assuming K independent Sharpe ratio 

477 estimates, each with the same variance. Uses the Gumbel approximation 

478 for the expected maximum of normals. 

479 

480 Args: 

481 number_of_trials: Number of strategies (K). 

482 variance: Variance of individual Sharpe ratio estimates. 

483 SR0: Baseline Sharpe ratio to add. Default 0. 

484 

485 Returns: 

486 Expected value of the maximum Sharpe ratio across K strategies. 

487 

488 Example: 

489 >>> # Expected maximum increases with number of trials 

490 >>> e1 = expected_maximum_sharpe_ratio(1, 0.1) 

491 >>> e10 = expected_maximum_sharpe_ratio(10, 0.1) 

492 >>> bool(e10 > e1) 

493 True 

494 """ 

495 return float( 

496 SR0 

497 + ( 

498 np.sqrt(variance) 

499 * ( 

500 (1 - np.euler_gamma) * scipy.stats.norm.ppf(1 - 1 / number_of_trials) 

501 + np.euler_gamma * scipy.stats.norm.ppf(1 - 1 / number_of_trials / np.exp(1)) 

502 ) 

503 ) 

504 ) 

505 

506 

507def minimum_track_record_length( 

508 SR: float, 

509 SR0: float, 

510 *, 

511 gamma3: float = 0.0, 

512 gamma4: float = 3.0, 

513 rho: float = 0.0, 

514 alpha: float = 0.05, 

515) -> float: 

516 """Compute minimum track record length for statistical significance. 

517 

518 Determines the minimum number of observations T required for the observed 

519 Sharpe ratio SR to be significantly greater than SR0 at confidence level 

520 1 - alpha. 

521 

522 Args: 

523 SR: Observed Sharpe ratio. 

524 SR0: Sharpe ratio under null hypothesis H0. 

525 gamma3: Skewness of returns. Default 0. 

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

527 rho: Autocorrelation of returns. Default 0. 

528 alpha: Significance level. Default 0.05. 

529 

530 Returns: 

531 Minimum track record length (number of periods) required. 

532 

533 Example: 

534 >>> # Higher Sharpe ratio needs shorter track record 

535 >>> mtrl_high = minimum_track_record_length(SR=1.0, SR0=0) 

536 >>> mtrl_low = minimum_track_record_length(SR=0.5, SR0=0) 

537 >>> bool(mtrl_high < mtrl_low) 

538 True 

539 """ 

540 var = sharpe_ratio_variance(SR0, T=1, gamma3=gamma3, gamma4=gamma4, rho=rho, K=1) 

541 return float(var * (scipy.stats.norm.ppf(1 - alpha) / (SR - SR0)) ** 2) 

542 

543 

544def make_expectation_gh( 

545 n_nodes: int = 200, 

546) -> Callable[[Callable[[np.ndarray], np.ndarray]], float]: 

547 """Create an expectation function using Gauss-Hermite quadrature. 

548 

549 Returns a function that computes E[g(Z)] where Z ~ N(0,1) using 

550 Gauss-Hermite quadrature with the specified number of nodes. 

551 

552 The approximation is: E[g(Z)] ≈ (1/√π) Σ w_i g(√2 t_i) 

553 where (t_i, w_i) are Gauss-Hermite nodes and weights. 

554 

555 Args: 

556 n_nodes: Number of quadrature nodes. Higher values increase accuracy. 

557 Default 200. 

558 

559 Returns: 

560 Function E(g) that computes E[g(Z)] for Z ~ N(0,1). 

561 

562 Example: 

563 >>> E = make_expectation_gh(n_nodes=100) 

564 >>> # E[Z^2] should be 1 for standard normal 

565 >>> bool(abs(E(lambda z: z**2) - 1.0) < 1e-6) 

566 True 

567 >>> # E[Z] should be 0 for standard normal 

568 >>> bool(abs(E(lambda z: z)) < 1e-10) 

569 True 

570 """ 

571 nodes, weights = np.polynomial.hermite.hermgauss(n_nodes) 

572 scale = np.sqrt(2.0) 

573 norm = 1.0 / np.sqrt(np.pi) 

574 x = scale * nodes 

575 

576 def E(g: Callable[[np.ndarray], np.ndarray]) -> float: 

577 """Compute E[g(Z)] for Z ~ N(0,1) via Gauss-Hermite quadrature. 

578 

579 Args: 

580 g: Function to integrate against the standard normal density. 

581 

582 Returns: 

583 Approximation of E[g(Z)]. 

584 """ 

585 vals = g(x) 

586 return float(norm * np.dot(weights, vals)) 

587 

588 return E 

589 

590 

591E_under_normal = make_expectation_gh(n_nodes=200) 

592 

593 

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

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

596 

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

598 This is the most conservative multiple testing correction. 

599 

600 Args: 

601 ps: Array of unadjusted p-values. 

602 

603 Returns: 

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

605 

606 Example: 

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

608 >>> adj_p = adjusted_p_values_bonferroni(p_vals) 

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

610 True 

611 """ 

612 M = len(ps) 

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

614 return result 

615 

616 

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

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

619 

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

621 than Bonferroni when tests are independent. 

622 

623 Args: 

624 ps: Array of unadjusted p-values. 

625 

626 Returns: 

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

628 

629 Example: 

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

631 >>> adj_p = adjusted_p_values_sidak(p_vals) 

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

633 np.True_ 

634 """ 

635 M = len(ps) 

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

637 

638 

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

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

641 

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

643 while still controlling the Family-Wise Error Rate. 

644 

645 Args: 

646 ps: Array of unadjusted p-values. 

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

648 (default) or "sidak". 

649 

650 Returns: 

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

652 

653 Raises: 

654 AssertionError: If variant is not "bonferroni" or "sidak". 

655 

656 Example: 

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

658 >>> adj_p = adjusted_p_values_holm(p_vals) 

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

660 0.03 

661 """ 

662 assert variant in ["bonferroni", "sidak"] 

663 i = np.argsort(ps) 

664 M = len(ps) 

665 p_adjusted = np.zeros(M) 

666 previous = 0 

667 for j, idx in enumerate(i): 

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

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

670 previous = p_adjusted[idx] 

671 return p_adjusted 

672 

673 

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

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

676 

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

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

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

680 

681 Args: 

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

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

684 SR1: Mean under alternative hypothesis. 

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

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

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

688 

689 Returns: 

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

691 or nan if no solution exists. 

692 

693 Raises: 

694 AssertionError: If parameters are out of valid ranges. 

695 

696 Example: 

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

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

699 True 

700 """ 

701 assert SR0 < SR1 

702 assert 0 < q < 1 

703 assert 0 < p_H1 < 1 

704 assert sigma0 > 0 

705 assert sigma1 > 0 

706 

707 with warnings.catch_warnings(): 

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

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

710 

711 def f(c: float) -> float: 

712 """Compute the posterior probability P[H=0 | X > c]. 

713 

714 Args: 

715 c: Candidate critical value. 

716 

717 Returns: 

718 Posterior false discovery probability at threshold c. 

719 """ 

720 a = 1 / ( 

721 1 

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

723 ) 

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

725 

726 if f(-10) < q: # Solution outside of the search interval 

727 return float(-np.inf) 

728 

729 if (f(-10) - q) * (f(10) - q) > 0: # No solution, for instance if σ₀≫σ₁ and q small 

730 return float(np.nan) 

731 

732 return float(scipy.optimize.brentq(lambda c: f(c) - q, -10, 10)) 

733 

734 

735def critical_sharpe_ratio( 

736 SR0: float, 

737 T: int, 

738 *, 

739 gamma3: float = 0.0, 

740 gamma4: float = 3.0, 

741 rho: float = 0.0, 

742 alpha: float = 0.05, 

743 K: int = 1, 

744) -> float: 

745 """Compute critical Sharpe ratio for hypothesis testing. 

746 

747 Determines the threshold SR_c for the one-sided test: 

748 H0: SR = SR0 vs H1: SR > SR0 

749 at significance level alpha. 

750 

751 Args: 

752 SR0: Sharpe ratio under null hypothesis. 

753 T: Number of observations. 

754 gamma3: Skewness of returns. Default 0. 

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

756 rho: Autocorrelation of returns. Default 0. 

757 alpha: Significance level. Default 0.05. 

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

759 

760 Returns: 

761 Critical Sharpe ratio threshold. Reject H0 if observed SR > SR_c. 

762 

763 Example: 

764 >>> SR_c = critical_sharpe_ratio(SR0=0, T=24, alpha=0.05) 

765 >>> bool(SR_c > 0) # Need positive SR to reject H0: SR=0 

766 True 

767 """ 

768 variance = sharpe_ratio_variance(SR0, T, gamma3=gamma3, gamma4=gamma4, rho=rho, K=K) 

769 return float(SR0 + scipy.stats.norm.ppf(1 - alpha) * math.sqrt(variance)) 

770 

771 

772def probabilistic_sharpe_ratio( 

773 SR: float, 

774 SR0: float, 

775 *, 

776 variance: float | None = None, 

777 T: int | None = None, 

778 gamma3: float = 0.0, 

779 gamma4: float = 3.0, 

780 rho: float = 0.0, 

781 K: int = 1, 

782) -> float: 

783 """Compute the Probabilistic Sharpe Ratio (PSR). 

784 

785 The PSR is 1 - p, where p is the p-value of testing H0: SR = SR0 vs 

786 H1: SR > SR0. It can be interpreted as a "Sharpe ratio on a probability 

787 scale", i.e., mapping the SR to [0, 1]. 

788 

789 Args: 

790 SR: Observed Sharpe ratio. 

791 SR0: Sharpe ratio under null hypothesis. 

792 variance: Variance of SR estimator. Provide this OR (T, gamma3, ...). 

793 T: Number of observations (if variance not provided). 

794 gamma3: Skewness of returns. Default 0. 

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

796 rho: Autocorrelation of returns. Default 0. 

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

798 

799 Returns: 

800 Probabilistic Sharpe Ratio in [0, 1]. Values near 1 indicate 

801 strong evidence that the true SR exceeds SR0. 

802 

803 Raises: 

804 AssertionError: If both variance and T are provided. 

805 

806 Example: 

807 >>> psr = probabilistic_sharpe_ratio(SR=0.5, SR0=0, T=24) 

808 >>> bool(0 < psr < 1) 

809 True 

810 >>> # Higher observed SR gives higher PSR 

811 >>> psr_high = probabilistic_sharpe_ratio(SR=1.0, SR0=0, T=24) 

812 >>> bool(psr_high > psr) 

813 True 

814 """ 

815 if variance is None: 

816 assert T is not None, "T must be provided if variance is not provided" 

817 variance = sharpe_ratio_variance(SR0, T, gamma3=gamma3, gamma4=gamma4, rho=rho, K=K) 

818 else: 

819 assert T is None, "Provide either the variance or (T, gamma3, gamma4, rho)" 

820 return float(scipy.stats.norm.cdf((SR - SR0) / math.sqrt(variance))) 

821 

822 

823def sharpe_ratio_power( 

824 SR0: float, 

825 SR1: float, 

826 T: int, 

827 *, 

828 gamma3: float = 0.0, 

829 gamma4: float = 3.0, 

830 rho: float = 0.0, 

831 alpha: float = 0.05, 

832 K: int = 1, 

833) -> float: 

834 """Compute statistical power for Sharpe ratio hypothesis test. 

835 

836 Calculates the power (1 - β) of the test H0: SR = SR0 vs H1: SR = SR1, 

837 which is the probability of correctly rejecting H0 when the true 

838 Sharpe ratio is SR1. 

839 

840 Note: Power is equivalent to recall in classification: 

841 Power = P[reject H0 | H1] = TP / (TP + FN) 

842 

843 Args: 

844 SR0: Sharpe ratio under null hypothesis. 

845 SR1: Sharpe ratio under alternative hypothesis (should be > SR0). 

846 T: Number of observations. 

847 gamma3: Skewness of returns. Default 0. 

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

849 rho: Autocorrelation of returns. Default 0. 

850 alpha: Significance level. Default 0.05. 

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

852 

853 Returns: 

854 Statistical power in [0, 1]. 

855 

856 Example: 

857 >>> # More observations increases power 

858 >>> power_short = sharpe_ratio_power(SR0=0, SR1=0.5, T=12) 

859 >>> power_long = sharpe_ratio_power(SR0=0, SR1=0.5, T=48) 

860 >>> bool(power_long > power_short) 

861 True 

862 """ 

863 critical_SR = critical_sharpe_ratio(SR0, T, gamma3=gamma3, gamma4=gamma4, rho=rho, alpha=alpha) 

864 variance = sharpe_ratio_variance(SR1, T, gamma3=gamma3, gamma4=gamma4, rho=rho, K=K) 

865 beta = scipy.stats.norm.cdf((critical_SR - SR1) / math.sqrt(variance)) 

866 return float(1 - beta) 

867 

868 

869def generate_autocorrelated_non_gaussian_data( 

870 N: int, 

871 n: int, 

872 SR0: float = 0, 

873 name: str = "gaussian", 

874 rho: float | None = None, 

875 gaussian_autocorrelation: float = 0, 

876) -> np.ndarray: 

877 """Generate autocorrelated non-Gaussian return data for simulation. 

878 

879 Creates a matrix of simulated returns with specified autocorrelation 

880 and marginal distribution characteristics (skewness/kurtosis). 

881 

882 Uses a copula-like approach: 

883 1. Generate AR(1) Gaussian processes 

884 2. Transform to uniform via Gaussian CDF 

885 3. Transform to target marginals via inverse CDF 

886 

887 Args: 

888 N: Number of time periods (rows). 

889 n: Number of assets/strategies (columns). 

890 SR0: Target Sharpe ratio. Default 0. 

891 name: Distribution type. One of "gaussian", "mild", "moderate", 

892 "severe". Default "gaussian". 

893 rho: Autocorrelation coefficient. If None, uses gaussian_autocorrelation. 

894 gaussian_autocorrelation: Autocorrelation for Gaussian case. Default 0. 

895 

896 Returns: 

897 Array of shape (N, n) containing simulated returns. 

898 

899 Example: 

900 >>> np.random.seed(42) 

901 >>> X = generate_autocorrelated_non_gaussian_data(100, 2, SR0=0.1, name="mild") 

902 >>> X.shape 

903 (100, 2) 

904 """ 

905 if rho is None: 

906 # With the distributions we consider the autocorrelation is almost the same. 

907 rho = gaussian_autocorrelation 

908 

909 shape = (N, n) 

910 

911 # Marginal distribution: ppf 

912 R = 10_000 

913 marginal = generate_non_gaussian_data(R, 1, SR0=SR0, name=name)[:, 0] 

914 ppf = scipy.interpolate.interp1d(ppoints(R), sorted(marginal), fill_value="extrapolate") 

915 

916 # AR(1) processes 

917 X = np.random.normal(size=shape) 

918 for i in range(1, shape[0]): 

919 X[i, :] = rho * X[i - 1, :] + np.sqrt(1 - rho**2) * X[i, :] 

920 

921 # Convert the margins to uniform, with the Gaussian cdf 

922 X = scipy.stats.norm.cdf(X) 

923 

924 # Convert the uniforms to the target margins, using the ppf 

925 result: np.ndarray = ppf(X) 

926 

927 return result 

928 

929 

930def get_random_correlation_matrix( 

931 number_of_trials: int = 100, 

932 effective_number_of_trials: int = 10, 

933 number_of_observations: int = 200, 

934 noise: float = 0.1, 

935) -> tuple[np.ndarray, np.ndarray, np.ndarray]: 

936 """Generate a random correlation matrix with block structure. 

937 

938 Creates a correlation matrix representing clustered strategies, where 

939 strategies within the same cluster are highly correlated and strategies 

940 across clusters have lower correlation. 

941