Handle missing return observations explicitly in shrinkage and semicovariance

docs/RiskModels.rst

@@ -163,3 +163,17 @@ References
 .. [3] Ledoit, O., & Wolf, M. (2003). `Honey, I Shrunk the Sample Covariance Matrix <http://www.ledoit.net/honey.pdf>`_ 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 <http://www.ledoit.net/ole2.pdf>`_, 10, 603–621.
 .. [5] Chen et al. (2010),  `Shrinkage Algorithms for MMSE Covariance Estimation <https://arxiv.org/pdf/0907.4698.pdf>`_, 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.

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)

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,
             ]
         ),
     )

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,
     )

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),
     )
 
 

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))