diff --git a/pypfopt/risk_models.py b/pypfopt/risk_models.py index 171ce0d0..7b4ae072 100644 --- a/pypfopt/risk_models.py +++ b/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) diff --git a/tests/test_risk_models.py b/tests/test_risk_models.py index 2b25072d..5a26c8ba 100644 --- a/tests/test_risk_models.py +++ b/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)