diff --git a/docs/RiskModels.rst b/docs/RiskModels.rst index 3fc8b152..5c56bf88 100644 --- a/docs/RiskModels.rst +++ b/docs/RiskModels.rst @@ -163,3 +163,17 @@ References .. [3] Ledoit, O., & Wolf, M. (2003). `Honey, I Shrunk the Sample Covariance Matrix `_ The Journal of Portfolio Management, 30(4), 110–119. https://doi.org/10.3905/jpm.2004.110 .. [4] Ledoit, O., & Wolf, M. (2001). `Improved estimation of the covariance matrix of stock returns with an application to portfolio selection `_, 10, 603–621. .. [5] Chen et al. (2010), `Shrinkage Algorithms for MMSE Covariance Estimation `_, IEEE Transactions on Signals Processing, 58(10), 5016-5029. + +Missing returns +--------------- + +Automatic Ledoit-Wolf and OAS estimates use the common window of fully observed +returns. They warn when incomplete rows are discarded and require at least two +complete observations. Missing values are not zero returns; explicitly choose +an imputation model before calling these estimators if dropping rows would lose +too much information. Manual shrinkage retains the pairwise sample covariance. + +``semicovariance`` averages each pair over jointly observed returns, including +observed zero returns. It raises an error if any pair has no joint observation. +As with pairwise sample covariance, the resulting matrix may need the existing +positive-semidefinite repair, which can alter individual entries. diff --git a/pypfopt/risk_models.py b/pypfopt/risk_models.py index 171ce0d0..6db87ae9 100644 --- a/pypfopt/risk_models.py +++ b/pypfopt/risk_models.py @@ -246,10 +246,13 @@ def semicovariance( returns = prices else: returns = returns_from_prices(prices, log_returns) + observed = returns.notna().astype(float) + counts = observed.T @ observed + if (counts == 0).any().any(): + raise ValueError("Every asset pair must have at least one joint observation") drops = np.fmin(returns - benchmark, 0) - T = drops.shape[0] return fix_nonpositive_semidefinite( - (drops.T @ drops) / T * frequency, kwargs.get("fix_method", "spectral") + (drops.T @ drops) / counts * frequency, kwargs.get("fix_method", "spectral") ) @@ -462,6 +465,24 @@ def __init__(self, prices, returns_data=False, frequency=252, log_returns=False) self.S = self.X.cov().values self.delta = None # shrinkage constant + def _complete_returns(self): + """Select common observations for estimators requiring a complete sample.""" + complete = self.X.dropna(how="any") + if len(complete) < 2: + raise ValueError( + "Shrinkage estimation requires at least two complete observations" + ) + dropped = len(self.X) - len(complete) + if dropped: + warnings.warn( + "Dropping {} incomplete return observations for shrinkage estimation".format( + dropped + ), + UserWarning, + stacklevel=3, + ) + return complete.values + def _format_and_annualize(self, raw_cov_array): """ Helper method which annualises the output of shrinkage calculations, @@ -529,7 +550,7 @@ def ledoit_wolf(self, shrinkage_target="constant_variance"): shrunk sample covariance matrix """ if shrinkage_target == "constant_variance": - X = np.nan_to_num(self.X.values) + X = self._complete_returns() shrunk_cov, self.delta = self.covariance.ledoit_wolf(X) elif shrinkage_target == "single_factor": shrunk_cov, self.delta = self._ledoit_wolf_single_factor() @@ -553,7 +574,7 @@ def _ledoit_wolf_single_factor(self): np.ndarray, float shrunk sample covariance matrix, shrinkage constant """ - X = np.nan_to_num(self.X.values) + X = self._complete_returns() # De-mean returns t, n = np.shape(X) @@ -610,10 +631,11 @@ def _ledoit_wolf_constant_correlation(self): np.ndarray, float shrunk sample covariance matrix, shrinkage constant """ - X = np.nan_to_num(self.X.values) + X = self._complete_returns() t, n = np.shape(X) - S = self.S # sample cov matrix + # Use the same complete observations for covariance and higher moments. + S = np.atleast_2d(np.cov(X, rowvar=False)) # Constant correlation target var = np.diag(S).reshape(-1, 1) @@ -663,6 +685,6 @@ def oracle_approximating(self): np.ndarray shrunk sample covariance matrix """ - X = np.nan_to_num(self.X.values) + X = self._complete_returns() shrunk_cov, self.delta = self.covariance.oas(X) return self._format_and_annualize(shrunk_cov) diff --git a/tests/test_black_litterman.py b/tests/test_black_litterman.py index 25bdff39..5e241048 100644 --- a/tests/test_black_litterman.py +++ b/tests/test_black_litterman.py @@ -180,26 +180,26 @@ def test_bl_returns_all_views(): posterior_rets, np.array( [ - 0.11168648, - 0.16782938, - 0.12516799, - 0.24067997, - 0.32848296, - -0.22789895, - 0.16311297, - 0.11928542, - 0.25414308, - 0.11007738, - 0.06282615, - -0.03140218, - -0.16977172, - 0.05254821, - -0.10463884, - 0.32173375, - 0.26399864, - 0.1118594, - 0.22999558, - 0.08977448, + 0.0659959179376, + 0.142395737648, + 0.0897272713478, + 0.203540995698, + 0.329077173837, + -0.194678572889, + 0.140774034267, + 0.106320649238, + 0.206079787474, + 0.110853248026, + 0.0604675037668, + -0.0334756812694, + -0.127710177297, + 0.0364737295266, + -0.0939289294804, + 0.294362055441, + 0.239905888459, + 0.0852597727498, + 0.210724156209, + 0.061193023013, ] ), ) diff --git a/tests/test_cla.py b/tests/test_cla.py index c5f69e23..a96cff78 100644 --- a/tests/test_cla.py +++ b/tests/test_cla.py @@ -106,7 +106,7 @@ def test_cla_max_sharpe_semicovariance(): np.testing.assert_almost_equal(cla.weights.sum(), 1) np.testing.assert_allclose( cla.portfolio_performance(risk_free_rate=0.02), - (0.2721798377099145, 0.07258537193305141, 3.474251505420551), + (0.328899626820958, 0.1707193786308789, 1.809399901160873), atol=1e-4, rtol=1e-4, ) diff --git a/tests/test_efficient_frontier.py b/tests/test_efficient_frontier.py index 038788e0..119eba3d 100644 --- a/tests/test_efficient_frontier.py +++ b/tests/test_efficient_frontier.py @@ -1257,7 +1257,7 @@ def test_max_sharpe_semicovariance(): assert all([i >= 0 for i in w.values()]) np.testing.assert_allclose( ef.portfolio_performance(risk_free_rate=0.02), - (0.2762965426962885, 0.07372667096108301, 3.476307004714425), + (0.3337156666549726, 0.1733001108784157, 1.810245042919043), ) @@ -1271,7 +1271,7 @@ def test_max_sharpe_short_semicovariance(): np.testing.assert_almost_equal(ef.weights.sum(), 1) np.testing.assert_allclose( ef.portfolio_performance(risk_free_rate=0.02), - (0.42444834528495234, 0.0898263632679403, 4.50255727350929), + (0.6050331853361551, 0.1756415312028965, 3.330836285299403), ) @@ -1291,7 +1291,7 @@ def test_min_volatility_shrunk_L2_reg(): assert all([i >= 0 for i in w.values()]) np.testing.assert_allclose( ef.portfolio_performance(risk_free_rate=0.02), - (0.17358178582309983, 0.19563960638632416, 0.7850239972361532), + (0.1811698577736233, 0.1511922398792176, 1.065992923329772), ) @@ -1308,7 +1308,7 @@ def test_efficient_return_shrunk(): assert all([i >= 0 for i in w.values()]) np.testing.assert_allclose( ef.portfolio_performance(risk_free_rate=0.02), - (0.22, 0.08892163937693882, 2.2491713085967766), + (0.22, 0.1263908584841322, 1.582392923022269), ) diff --git a/tests/test_missing_return_estimation.py b/tests/test_missing_return_estimation.py new file mode 100644 index 00000000..748b2d90 --- /dev/null +++ b/tests/test_missing_return_estimation.py @@ -0,0 +1,75 @@ +import numpy as np +import pandas as pd +import pytest + +from pypfopt import risk_models + + +def estimate(data, method): + model = risk_models.CovarianceShrinkage(data, returns_data=True, frequency=1) + return ( + model.oracle_approximating() if method == "oas" else model.ledoit_wolf(method) + ) + + +@pytest.mark.parametrize( + "method", ["constant_variance", "single_factor", "constant_correlation", "oas"] +) +def test_shrinkage_uses_complete_observations(method): + rng = np.random.default_rng(52) + returns = pd.DataFrame(rng.normal(size=(120, 3)) / 100, columns=list("ABC")) + returns.loc[:79, "B"] = np.nan + returns.loc[95, "C"] = np.nan + original = returns.copy() + with pytest.warns(UserWarning, match="81.*incomplete"): + actual = estimate(returns, method) + expected = estimate(returns.dropna(), method) + pd.testing.assert_frame_equal(actual, expected) + pd.testing.assert_frame_equal(returns, original) + + +@pytest.mark.parametrize( + "method", ["constant_variance", "single_factor", "constant_correlation", "oas"] +) +@pytest.mark.parametrize("rows", [0, 1]) +def test_shrinkage_requires_two_complete_observations(method, rows): + returns = pd.DataFrame({"A": [0.1, 0.2, 0.3], "B": [np.nan, np.nan, np.nan]}) + if rows: + returns.loc[2, "B"] = 0.1 + with pytest.raises(ValueError, match="two complete"): + estimate(returns, method) + + +def test_semicovariance_uses_pairwise_observation_counts(): + returns = pd.DataFrame( + {"A": [-0.1, 0.2, np.nan, -0.2], "B": [np.nan, -0.1, -0.1, 0.4]} + ) + actual = risk_models.semicovariance( + returns, returns_data=True, benchmark=0, frequency=1 + ) + expected = pd.DataFrame( + [[1 / 60, 0], [0, 1 / 150]], index=returns.columns, columns=returns.columns + ) + pd.testing.assert_frame_equal(actual, expected) + + +def test_semicovariance_keeps_observed_zero_returns(): + returns = pd.DataFrame({"A": [-0.1, 0.0, -0.2], "B": [0.0, -0.1, 0.0]}) + actual = risk_models.semicovariance( + returns, returns_data=True, benchmark=0, frequency=1 + ) + np.testing.assert_allclose(actual, [[1 / 60, 0], [0, 1 / 300]]) + + +def test_semicovariance_rejects_pairs_without_joint_observations(): + returns = pd.DataFrame({"A": [0.1, np.nan], "B": [np.nan, 0.2]}) + with pytest.raises(ValueError, match="joint observation"): + risk_models.semicovariance(returns, returns_data=True) + + +def test_semicovariance_uses_joint_counts_off_diagonal(): + returns = pd.DataFrame({"A": [-0.1, -0.1, np.nan], "B": [-0.1, np.nan, -0.1]}) + actual = risk_models.semicovariance( + returns, returns_data=True, benchmark=0, frequency=1 + ) + np.testing.assert_allclose(actual, np.full((2, 2), 0.01))