Handle missing return observations explicitly in shrinkage and semicovariance

docs/RiskModels.rst

@@ -163,3 +163,53 @@ 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.
+
+.. warning::
+
+    This changes results for inputs containing missing returns. Previous
+    versions treated missing returns as zero downside observations and divided
+    all semicovariance entries by the total number of rows. Saved portfolio
+    weights, reported risk and Sharpe ratios can therefore change materially;
+    recompute them and review the chosen missing-data policy. Fully observed
+    inputs retain the same semicovariance formula.
+
+    On the repository stock-price fixture with ``benchmark=0`` and a 2% risk-free
+    rate, the long-only ``EfficientFrontier.max_sharpe`` example's reported risk
+    changes from approximately 0.0737 to 0.1733 and its Sharpe ratio from 3.48 to
+    1.81. These are in-sample outputs for different optimized portfolios, using
+    a semicovariance risk matrix; they are not realized volatility, a universal
+    multiplier, or evidence of better out-of-sample performance.
+
+Pairwise deletion does not guarantee a positive-semidefinite matrix. For example,
+the three return rows ``(-.1, -.1)``, ``(0, NaN)``, ``(NaN, 0)`` at benchmark zero
+and frequency one produce raw entries ``[[.005, .01], [.01, .005]]``. Its
+eigenvalues are -.005 and .015. The default spectral repair warns and sets the
+negative eigenvalue to zero, producing .0075 in every entry: each diagonal
+increases by .0025 (50%). Thus a repaired diagonal is no longer necessarily the
+observed marginal downside moment. ``fix_method="diag"`` instead adds .0055 to
+each diagonal in this example. Neither repair replaces a considered missing-data
+model or restores observations that were never recorded.
+
+There is no configurable minimum joint-observation count: a pair with one
+shared return is allowed, while a pair with no shared returns raises an error.
+Different entries can therefore be estimated from very different sample sizes.
+In the example above, the off-diagonal uses one observation and each diagonal
+uses two; dividing the off-diagonal by the marginal scales gives 2, not a valid
+correlation. Spectral repair restores positive semidefiniteness but does not
+make that one-observation estimate reliable or add missing information. Check
+joint-observation counts and choose a missing-data policy before using the
+matrix for portfolio optimization.

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,120 @@
+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))
+
+
+@pytest.mark.parametrize("frequency", [1, 252])
+@pytest.mark.parametrize(
+    "fix_method, expected",
+    [
+        ("spectral", [[0.0075, 0.0075], [0.0075, 0.0075]]),
+        ("diag", [[0.0105, 0.01], [0.01, 0.0105]]),
+    ],
+)
+def test_semicovariance_repairs_indefinite_pairwise_estimate(
+    fix_method, expected, frequency
+):
+    returns = pd.DataFrame({"A": [-0.1, 0.0, np.nan], "B": [-0.1, np.nan, 0.0]})
+    original = returns.copy()
+    # Each variance uses two observations; the cross-product uses one.
+    # Raw matrix [[.005, .01], [.01, .005]] has eigenvalues -.005 and .015.
+    # Spectral clipping gives .015/2 in every entry (diagonal increases 50%).
+    # Diagonal repair instead adds 1.1 * .005 to both diagonal entries.
+    with pytest.warns(UserWarning, match="non positive semidefinite") as caught:
+        actual = risk_models.semicovariance(
+            returns,
+            returns_data=True,
+            benchmark=0,
+            frequency=frequency,
+            fix_method=fix_method,
+        )
+    assert len(caught) == 1
+    expected = pd.DataFrame(
+        np.array(expected) * frequency, index=returns.columns, columns=returns.columns
+    )
+    pd.testing.assert_frame_equal(actual, expected)
+    assert np.linalg.eigvalsh(actual).min() >= -1e-12
+    pd.testing.assert_frame_equal(returns, original)
+
+
+def test_semicovariance_complete_data_retains_full_sample_formula():
+    returns = pd.DataFrame({"A": [-0.1, 0.02, 0.0], "B": [0.03, -0.05, -0.02]})
+    benchmark = 0.01
+    downside = np.minimum(returns.to_numpy() - benchmark, 0)
+    expected = downside.T @ downside / len(returns) * 252
+    actual = risk_models.semicovariance(
+        returns, returns_data=True, benchmark=benchmark, frequency=252
+    )
+    np.testing.assert_allclose(actual, expected, rtol=1e-12, atol=1e-12)