[kadyrbekovhamit-cyber]

Automatic shrinkage estimators currently replace missing returns with zero. semicovariance also converts NaNs to zero via fmin and divides every pair by the full row count. An unobserved return is therefore treated as an observed flat day.

Proposed policy:

• Ledoit–Wolf (all three targets) and OAS use the common complete-return sample, warn with the number of dropped rows, and require at least two complete observations. Constant-correlation covariance and higher moments use that same sample, retaining the existing ddof convention. Manual shrunk_covariance(delta) remains unchanged.
• Semicovariance divides each downside-product sum by that asset pair's joint-observation count, rejects pairs with no overlap, and retains the existing PSD-repair option. Actual zero returns still count as observations.

This makes a statistical choice explicit: complete-case shrinkage can discard much of the older history; pairwise estimates can need PSD repair. The docs state these limitations. I would welcome feedback on the policy/API before merge.

Six old integration goldens change because they included the zero-imputation result. I cross-checked the updated moments with pairwise means and sklearn, the BL posterior with its precision-form expression, and five optimized objective values with an independent scipy SLSQP solve (largest objective gap below 5e-13), before updating those goldens.

Validation: original initial regressions 14 failed / 1 passed; patched risk_models + 16 regressions + Black–Litterman/CLA/efficient-frontier modules: 146 passed / 4 skipped. Tests include internal/leading gaps, input immutability, insufficient overlap, observed zeros and off-diagonal joint counts. On the repository stock-price fixture at benchmark=0, BABA downside variance is .00560896 with the old denominator versus .04465234 using observed days; the long-only minimum-risk weight changes from .532475 to .081924. This is a fixture example, not observed investment losses or a universal multiplier.

Related prior discussion: #270 (common history recommendation), #294 (zero-fill criticism in a different estimator), #500/#514 (new-stock overweight symptoms). I do not claim that the general missing-data issue is new. No matching current patch found in the reviewed open PRs. Tested main a6638d2, Python 3.12 / numpy 2.2.6 / pandas 3.0.6. AI-assisted investigation and patch, one numerical worker.

Review follow-up / compatibility impact — 24 October 2026

Thanks to @arhancanli for flagging both the compatibility impact and the missing PSD-repair regression. Independently rerunning the repository stock-price fixture with benchmark=0 and risk-free rate=0.02 gives:

Example	Previous reported risk	Proposed reported risk	Previous Sharpe	Proposed Sharpe
EfficientFrontier max_sharpe	0.0737267	0.1733001	3.476307	1.810245
CLA max_sharpe	0.0725854	0.1707194	3.474252	1.809400

These are in-sample risk measures from a semicovariance matrix for different optimized portfolios, not realized volatility or a universal effect. The fixture has 7,125 return rows but only 895 observed BABA returns: dividing BABA's downside-product sum by all rows suppresses that entry by 7,125/895 = 7.96089. A separate pairwise-mean calculation matches the proposed matrix. Both old and new fixture matrices are already PSD, so this particular large shift is not caused by spectral repair.

New end-to-end tests use return rows (-.1,-.1), (0,NaN), (NaN,0). At benchmark=0 and frequency=1, raw pairwise semicovariance is [[.005,.01],[.01,.005]], with eigenvalues -.005 and .015. Default spectral repair must warn and return .0075 in every entry, increasing each diagonal by .0025 (50%). The alternative diagonal repair returns diagonal .0105 and off-diagonal .01. Tests cover both methods and frequencies 1/252, PSD, labels, input immutability, and a fully observed compatibility control. The four repair cases fail against the unmodified base estimator; the compatibility control passes. All 21 missing-data regressions now pass; affected risk-model, BL, CLA and efficient-frontier modules: 151 passed, 4 skipped (existing warning summary retained).

Proposed release note (for the next release, if merged)

Behavior change for missing-return inputs: automatic Ledoit–Wolf/OAS estimation uses complete observations instead of implicit zero filling, with a warning when rows are discarded. Semicovariance uses joint observation counts instead of the full-row denominator. Existing portfolio weights, risk and Sharpe results can change materially; recompute saved results and review the missing-data policy. Fully observed semicovariance retains its prior formula. Pairwise estimates can be indefinite; PSD repair can change even marginal downside moments, as the documented 50% diagonal-change example demonstrates. This is an explicit estimation policy, not a claim of universally better out-of-sample performance.

The project publishes release notes on GitHub Releases; the paragraph above is supplied for maintainers to include upon release. No release is being announced by this PR.

[arhancanli]

Read the diff (did not run the suite). The semicovariance change divides each pair by its joint observation count instead of T, and the updated expectations show how much that moves results on the repo's test prices: test_max_sharpe_semicovariance goes from volatility 0.0737 / Sharpe 3.48 to 0.1733 / 1.81, and test_cla_max_sharpe_semicovariance from 0.0726 / 3.47 to 0.1707 / 1.81. That is a 2x change in reported risk for existing users of risk_models.semicovariance on data with gaps, so it deserves a line in the PR description and the release notes, not only a docs section. The old divisor-T version effectively treated missing returns as zero drops and shrank the matrix; the new one is closer to right, but anyone with a saved max-Sharpe result will see it change.

One question the tests could answer: with the pairwise counts, entries are averaged over different windows, so the matrix is not guaranteed PSD (the docs mention the spectral fix). Is there a test showing fix_nonpositive_semidefinite actually fires on a gappy input and by how much it moves the diagonal?

[kadyrbekovhamit-cyber]

Thanks @arhancanli — agreed on both points. I reran the actual fixture: EfficientFrontier reported risk changes 0.0737267 → 0.1733001 (Sharpe 3.476307 → 1.810245); CLA changes 0.0725854 → 0.1707194 (Sharpe 3.474252 → 1.809400). The PR body now highlights that compatibility change and includes a proposed release-note paragraph, and the docs carry an explicit warning. These are in-sample semicovariance-based outputs for different portfolios, not realized risk or a universal multiplier. The fixture itself needs no PSD repair; its change comes from observation counts.

The requested end-to-end regression is now in test_missing_return_estimation.py: return rows (-.1,-.1), (0,NaN), (NaN,0) give raw [[.005,.01],[.01,.005]] at benchmark=0, frequency=1. Its eigenvalues are -.005 and .015. Spectral repair must actually warn and produce .0075 in every entry: each diagonal rises by .0025, or 50%. The test checks that exact matrix, PSD, labels and unchanged inputs; it also covers fix_method="diag" (diagonal .0105), annualization, and fully observed compatibility.

The four new repair cases fail against the unmodified base estimator and pass on the proposal. Current results: 21 missing-data regressions passed; affected risk-model, BL, CLA and efficient-frontier modules 151 passed / 4 skipped. Pairwise-moment values on the stock fixture were also independently recomputed. I have left the policy trade-off explicit: pairwise deletion and spectral repair do not recover missing data or establish better out-of-sample performance.

[arhancanli]

I checked the new end-to-end case by hand and it holds. For rows (-.1,-.1), (0,NaN), (NaN,0) at benchmark 0, the off-diagonal has a single joint observation (.1 × .1 = .01) and each diagonal averages over its own two observations (.01 / 2 = .005), so the raw matrix is [[.005,.01],[.01,.005]]. Its eigenvalues are -.005 and .015. Clipping the negative one leaves .015 · ½[[1,1],[1,1]], so every entry is .0075 and each diagonal rises by .0025, which is the 50% you describe. Thanks for stating that the 0.0737 → 0.1733 change comes from the observation counts and not from any repair.

This example also shows where the non-PSD matrix comes from: a pair with one joint observation yields a correlation of .01 / .005 = 2. pandas.DataFrame.cov has min_periods for this, and entries below it come back NaN. Does the PR have a minimum joint-observation count for a pair, or is a pair estimated from a single shared return allowed? If it is allowed, I would document that next to the spectral-repair warning, because the repair then moves the numbers most for exactly those pairs.