feat: add Ledoit-Wolf 2020 analytical nonlinear shrinkage estimator

pypfopt/risk_models.py

@@ -17,6 +17,7 @@
     - manual shrinkage
     - Ledoit Wolf shrinkage
     - Oracle Approximating shrinkage
+    - Analytical Nonlinear shrinkage (Ledoit and Wolf, 2020)
 
 - covariance to correlation matrix
 """
@@ -136,6 +137,7 @@ def risk_matrix(prices, method="sample_cov", **kwargs):
         - ``ledoit_wolf_single_factor``
         - ``ledoit_wolf_constant_correlation``
         - ``oracle_approximating``
+        - ``analytical_nonlinear_shrinkage``
 
     Raises
     ------
@@ -165,6 +167,8 @@ def risk_matrix(prices, method="sample_cov", **kwargs):
         )
     elif method == "oracle_approximating":
         return CovarianceShrinkage(prices, **kwargs).oracle_approximating()
+    elif method == "analytical_nonlinear_shrinkage":
+        return CovarianceShrinkage(prices, **kwargs).analytical_nonlinear_shrinkage()
     else:
         raise NotImplementedError("Risk model {} not implemented".format(method))
 
@@ -666,3 +670,80 @@ def oracle_approximating(self):
         X = np.nan_to_num(self.X.values)
         shrunk_cov, self.delta = self.covariance.oas(X)
         return self._format_and_annualize(shrunk_cov)
+
+    def analytical_nonlinear_shrinkage(self):
+        """
+        Nonlinear shrinkage estimator from Ledoit and Wolf (2020).
+
+        The standard sample covariance matrix is a bad estimator when you have
+        many assets relative to your history length -- it systematically blows up
+        the large eigenvalues and squashes the small ones. This method corrects
+        each eigenvalue individually using a closed-form formula derived from
+        random matrix theory, rather than applying one global shrinkage intensity
+        like ledoit_wolf() does.
+
+        Reference: Ledoit, O. and Wolf, M. (2020). Analytical Nonlinear Shrinkage
+        of Large-Dimensional Covariance Matrices. Annals of Statistics, 48(5).
+        The algorithm here follows the 7-step summary in Section 4.7.
+
+        :raises ValueError: if you have more assets than observations (p >= n),
+            since the sample covariance is singular and the formula breaks down.
+        :return: annualised shrunk covariance matrix
+        :rtype: pd.DataFrame
+        """
+        n, p = self.X.shape
+
+        if p >= n:
+            raise ValueError(
+                f"Need more observations than assets (p < n), but got "
+                f"p={p} and n={n}. Either use fewer assets or a longer lookback."
+            )
+
+        # Decompose the sample covariance: S = U * diag(lam) * U.T
+        # eigh is used instead of eig because S is symmetric -- gives real,
+        # sorted eigenvalues and is numerically more stable.
+        lam, u = np.linalg.eigh(self.S)
+
+        # Global bandwidth (Eq. 4.4). n^(-1/3) is the theoretically optimal
+        # choice for this kernel in the large-dimensional limit.
+        h = n ** (-1.0 / 3.0)
+
+        # Each eigenvalue gets its own bandwidth proportional to its size (Eq. 4.5).
+        # This is what makes the kernel locally adaptive -- small eigenvalues
+        # get a finer window, large ones a coarser window.
+        h_loc = lam * h
+
+        # Build the (p x p) matrix of standardised distances.
+        # u_mat[i, j] = (lam[i] - lam[j]) / h_loc[j]
+        u_mat = (lam[:, np.newaxis] - lam[np.newaxis, :]) / h_loc[np.newaxis, :]
+
+        # Epanechnikov kernel: k(u) = 0.75*(1 - u^2) for |u| <= 1, else 0.
+        epan = np.where(np.abs(u_mat) <= 1.0, 0.75 * (1.0 - u_mat ** 2), 0.0)
+        f_tilde = np.mean(epan / h_loc[np.newaxis, :], axis=1)
+
+        # Hilbert transform of the KDE (Proposition 4.1).
+        # The log blows up at the kernel boundary |u|=1, but the full expression
+        # has a finite limit there -- we just zero out any non-finite values.
+        with np.errstate(divide="ignore", invalid="ignore"):
+            log_term = 0.5 * (1.0 - u_mat ** 2) * np.log(
+                np.abs((1.0 + u_mat) / (1.0 - u_mat))
+            )
+        log_term = np.where(np.isfinite(log_term), log_term, 0.0)
+        hf_tilde = np.mean(
+            (3.0 / (4.0 * np.pi)) * (log_term + u_mat) / h_loc[np.newaxis, :], axis=1
+        )
+
+        # The shrinkage formula (Section 4.7).
+        # c = p/n is the concentration ratio -- the closer it is to 1, the harder
+        # the correction needs to be.
+        c = p / n
+        d_tilde = lam / (
+            (np.pi * c * lam * f_tilde) ** 2
+            + (1.0 - c - np.pi * c * lam * hf_tilde) ** 2
+        )
+
+        # Reconstruct and symmetrise to clean up any floating-point drift.
+        shrunk_cov = u @ np.diag(d_tilde) @ u.T
+        shrunk_cov = (shrunk_cov + shrunk_cov.T) / 2.0
+
+        return self._format_and_annualize(shrunk_cov)

tests/test_risk_models.py

@@ -310,6 +310,30 @@ def test_oracle_approximating():
     assert risk_models._is_positive_semidefinite(shrunk_cov)
 
 
+def test_analytical_nonlinear_shrinkage():
+    df = get_data()
+    cs = risk_models.CovarianceShrinkage(df)
+    shrunk_cov = cs.analytical_nonlinear_shrinkage()
+    assert shrunk_cov.shape == (20, 20)
+    assert list(shrunk_cov.index) == list(df.columns)
+    assert list(shrunk_cov.columns) == list(df.columns)
+    assert not shrunk_cov.isnull().any().any()
+    assert risk_models._is_positive_semidefinite(shrunk_cov)
+    # make sure the correction actually did something
+    assert not np.allclose(shrunk_cov.values, risk_models.sample_cov(df).values)
+
+
+def test_analytical_nonlinear_shrinkage_singular_raises():
+    # more assets than observations -- formula breaks down, should raise
+    np.random.seed(0)
+    prices = pd.DataFrame(
+        np.random.randn(5, 20).cumsum(axis=0) + 100,
+        columns=[f"Asset{i}" for i in range(20)],
+    )
+    with pytest.raises(ValueError, match="p < n"):
+        risk_models.CovarianceShrinkage(prices).analytical_nonlinear_shrinkage()
+
+
 def test_risk_matrix_and_returns_data():
     # Test the switcher method for simple calls
     df = get_data()
@@ -323,6 +347,7 @@ def test_risk_matrix_and_returns_data():
         "ledoit_wolf_single_factor",
         "ledoit_wolf_constant_correlation",
         "oracle_approximating",
+        "analytical_nonlinear_shrinkage",
     }:
         S = risk_models.risk_matrix(df, method=method)
         assert S.shape == (20, 20)