Coverage for src/fast_minimum_variance/shrinkage/util.py: 100%

1"""Shrinkage intensity and target utilities (LW, OAS, RMT, constant-correlation).""" 

2 

3import numpy as np 

4from sklearn.covariance import ledoit_wolf, oas 

5 

6 

7def lw_alpha_and_target(X: np.ndarray) -> tuple[float, np.ndarray]: # noqa: N803 

8 """Return (alpha_lw, target) for LW scaled-identity shrinkage via sklearn. 

9 

10 X must be column-demeaned. The shrinkage target is bar_lambda * I where 

11 bar_lambda = ||X||_F^2 / (n * T) is the average per-asset variance 

12 (equivalently, the mean eigenvalue of the sample covariance X^T X / T). 

13 """ 

14 _, alpha = ledoit_wolf(X, assume_centered=True) 

15 n = X.shape[1] 

16 T = X.shape[0] # noqa: N806 

17 bar_lam = float(np.linalg.norm(X, "fro") ** 2) / (n * T) 

18 cov = (X.transpose() @ X) / T 

19 bar_lam_cov = np.linalg.trace(cov) / n 

20 assert np.isclose(bar_lam, bar_lam_cov) # noqa: S101 

21 

22 delta2 = np.linalg.norm(cov - bar_lam * np.eye(n), "fro") ** 2 

23 norms_sq = np.sum(X**2, axis=1) 

24 beta2 = (np.sum(norms_sq**2) - 2 * np.sum(X * (X @ cov)) + T * np.linalg.norm(cov, "fro") ** 2) / T**2 

25 alpha_manual = min(1.0, beta2 / delta2) 

26 assert np.isclose(alpha_manual, alpha) # noqa: S101 

27 

28 return float(alpha), bar_lam * np.eye(n) 

29 

30 

31def lw_alpha_and_target_hard(X: np.ndarray, alpha: float = 0.5) -> tuple[float, np.ndarray]: # noqa: N803 

32 """Return (alpha, target) for scaled-identity shrinkage with a fixed alpha. 

33 

34 X must be column-demeaned. The shrinkage target is bar_lambda * I where 

35 bar_lambda = ||X||_F^2 / (n * T) is the average per-asset variance 

36 (equivalently, the mean eigenvalue of the sample covariance X^T X / T). 

37 """ 

38 # _, alpha = ledoit_wolf(X, assume_centered=True) 

39 n = X.shape[1] 

40 T = X.shape[0] # noqa: N806 

41 bar_lam = float(np.linalg.norm(X, "fro") ** 2) / (n * T) 

42 

43 return alpha, bar_lam * np.eye(n) 

44 

45 

46def oas_alpha_and_target(X: np.ndarray) -> tuple[float, np.ndarray]: # noqa: N803 

47 """Return (alpha_oas, target) for OAS scaled-identity shrinkage via sklearn. 

48 

49 Uses the same bar_lambda * I target as LW but the Oracle Approximating 

50 Shrinkage formula (Chen et al. 2010), which has lower MSE when n/T is 

51 non-negligible. 

52 """ 

53 _, alpha = oas(X, assume_centered=True) 

54 n = X.shape[1] 

55 bar_lam = float(np.linalg.norm(X, "fro") ** 2) / (n * X.shape[0]) 

56 return float(alpha), bar_lam * np.eye(n) 

57 

58 

59def cc_target(X: np.ndarray) -> tuple[np.ndarray, float]: # noqa: N803 

60 """Constant-correlation shrinkage target (Ledoit-Wolf 2004 JoPM). 

61 

62 T0_ij = rho_bar * sigma_i * sigma_j (i != j) 

63 T0_ii = sigma_i^2 

64 

65 where sigma_i = sqrt(Sigma_ii) is the per-asset sample standard deviation 

66 and rho_bar is the mean off-diagonal sample correlation coefficient. 

67 Always PSD for 0 < rho_bar < 1. 

68 

69 Returns (target, rho_bar). 

70 """ 

71 T, n = X.shape # noqa: N806 

72 cov = (X.T @ X) / T 

73 std = np.sqrt(np.diag(cov)) 

74 corr = cov / np.outer(std, std) 

75 np.fill_diagonal(corr, 0.0) 

76 rho_bar = float(corr.sum()) / (n * (n - 1)) 

77 target = rho_bar * np.outer(std, std) 

78 np.fill_diagonal(target, std**2) 

79 return target, rho_bar 

80 

81 

82def lw_alpha_for_target(X: np.ndarray, target: np.ndarray) -> float: # noqa: N803 

83 """LW oracle alpha for an arbitrary SPD shrinkage target T0. 

84 

85 alpha* = min(1, beta2 / delta2) where 

86 delta2 = ||S - T0||_F^2 (distance of sample cov from target) 

87 beta2 = (1/T^2) sum_t ||x_t x_t' - S||_F^2 (noise in S) 

88 """ 

89 T, _n = X.shape # noqa: N806 

90 cov = (X.T @ X) / T 

91 delta2 = float(np.linalg.norm(cov - target, "fro") ** 2) 

92 norms_sq = np.sum(X**2, axis=1) 

93 beta2 = float((np.sum(norms_sq**2) - 2 * np.sum(X * (X @ cov)) + T * np.linalg.norm(cov, "fro") ** 2) / T**2) 

94 return min(1.0, beta2 / delta2) 

95 

96 

97def rmt_target_and_alpha( 

98 X: np.ndarray, # noqa: N803 

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

100 """RMT-clipped shrinkage target with alpha=1. 

101 

102 Eigenvalues of the sample covariance above the Marchenko-Pastur bulk edge 

103 are kept as-is (signal); all others are clipped to bar_lambda (noise floor). 

104 The resulting target has lambda_min = bar_lambda (same as the scaled-identity 

105 target) so it provides equally effective lambda_min lifting at any alpha, while 

106 being 9x closer to the sample covariance in Frobenius norm than bar_lambda * I. 

107 

108 T0 = bar_lambda * I + U_k @ diag(lambda_k - bar_lambda) @ U_k^T 

109 

110 where (U_k, lambda_k) are the k eigenpairs of S = X^T X / T whose eigenvalues 

111 exceed the MP upper edge bar_lambda * (1 + sqrt(n/T))^2. 

112 

113 Returns (target, lr_factors, k, 1.0). alpha=1 means the system matrix is 

114 T0^RMT directly; the _kkt_step Woodbury path applies this in O(n_a k + k^3). 

115 """ 

116 T, n = X.shape # noqa: N806 

117 cov = (X.T @ X) / T 

118 bar_lam = np.trace(cov) / n 

119 mp_upper = bar_lam * (1.0 + np.sqrt(n / T)) ** 2 

120 

121 eigs, vecs = np.linalg.eigh(cov) # ascending order 

122 signal = eigs > mp_upper 

123 k = int(signal.sum()) 

124 eigs_k = eigs[signal] 

125 vecs_k = vecs[:, signal] 

126 

127 delta_k = eigs_k - bar_lam # (k,) eigenvalue excesses 

128 target = bar_lam * np.eye(n) + vecs_k @ np.diag(delta_k) @ vecs_k.T 

129 lr_factors = (float(bar_lam), vecs_k, delta_k) # for O(nk) matvec 

130 return target, lr_factors, k, 1.0