From 6922ca23685bcaea0b3ce9f2416777e1303fc86c Mon Sep 17 00:00:00 2001 From: PureSaber Date: Wed, 14 Oct 2026 22:43:42 +0800 Subject: [PATCH 1/6] feat: solve joint portfolio constraints with auditable convex backends --- .github/workflows/ci.yml | 6 +- README.md | 77 ++- pyproject.toml | 7 +- requirements.lock | 91 +-- scripts/benchmark_optimizer.py | 96 ++++ src/quant_portfolio/__init__.py | 10 + src/quant_portfolio/allocation_v2.py | 386 +++++++++++++ src/quant_portfolio/executable.py | 225 ++++++++ src/quant_portfolio/objectives.py | 311 ++--------- src/quant_portfolio/optimization.py | 763 +++++++++----------------- src/quant_portfolio/risk_snapshot.py | 191 +++++++ src/quant_portfolio/solver.py | 589 ++++++++++++++++++++ tests/test_governance.py | 8 +- tests/test_multiperiod_constraints.py | 17 +- tests/test_optimization.py | 9 +- tests/test_solver_failure_contract.py | 107 ++++ tests/test_solver_v2.py | 362 ++++++++++++ 17 files changed, 2418 insertions(+), 837 deletions(-) create mode 100644 scripts/benchmark_optimizer.py create mode 100644 src/quant_portfolio/allocation_v2.py create mode 100644 src/quant_portfolio/executable.py create mode 100644 src/quant_portfolio/risk_snapshot.py create mode 100644 src/quant_portfolio/solver.py create mode 100644 tests/test_solver_failure_contract.py create mode 100644 tests/test_solver_v2.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f598495..779881c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,9 +11,10 @@ permissions: jobs: test: - runs-on: ubuntu-24.04 + runs-on: ${{ matrix.os }} strategy: matrix: + os: [ubuntu-24.04, windows-2025] python-version: ["3.10", "3.11", "3.12"] steps: @@ -30,7 +31,7 @@ jobs: run: | python -m pip install --no-deps --requirement requirements.lock python -m pip check - python -m pip install --no-deps --no-build-isolation --editable . + python -m pip install --no-deps --no-build-isolation . python -m pip check - name: Ruff check @@ -66,6 +67,7 @@ jobs: print(f"cross_asset.py pure branch coverage: {covered}/{total} ({covered / total:.2%})") - name: Run deterministic synthetic demo + shell: bash run: | python scripts/run_synthetic_demo.py --out "$RUNNER_TEMP/synthetic-demo" grep -q "SYNTHETIC DEMONSTRATION" "$RUNNER_TEMP/synthetic-demo/metrics.json" diff --git a/README.md b/README.md index 73dca9d..18995b5 100644 --- a/README.md +++ b/README.md @@ -173,9 +173,82 @@ it emits QExec `OrderIntent` suggestions and never alters a ledger, positions, o The module has no live-order, network, or credential capability. Missing, future, duplicate, or non-finite PIT inputs fail closed. -## M6 governance and reproducibility +## Optimizer protocol v2 and migration + +`allocate_research` returns an `AllocationResult` with continuous account weights, +independent constraint residuals, forecast/risk units, exact solver versions and +backend attempts. `research_allocation_weights` remains a Series convenience API. +New configurations explicitly specify `protocol_version: 2`, `mode`, +`risk_objective: absolute|active`, `forecast_horizon_sessions` and +`expected_return_model: ic_vol|provided|none`. For example: + +```yaml +allocation: + protocol_version: 2 + mode: cost_aware + risk_objective: absolute + forecast_horizon_sessions: 5 + expected_return_model: ic_vol + lookback: 120 + min_observations: 20 + ic_lookback: 120 + ic_min_periods: 20 + max_turnover: 2.0 + timeout_seconds: 30 +``` -Version `0.4.2` consumes only published annotated internal tags: +The expected-return and risk penalty use the declared decision horizon. Costs +are charged once per rebalance. Covariance inputs and tracking-error limits are +annualized. IC uses only labels available and matured by the decision, selects +non-overlapping periods and fails on insufficient history. `provided` forecasts +must cover every candidate. Black–Litterman remains an explicit lower-level +views API; v2 does not infer missing views or market weights. + +Mean–variance maximizes forecast return minus half the risk-aversion coefficient +times variance and trading costs. `active` centers risk on the full benchmark +and is currently supported by `cost_aware`. Target methods minimize squared +distance to equal, inverse-volatility, risk-parity or HRP targets plus declared +cost penalties; risk-contribution changes are reported. Ratio methods maximize +the positive excess-return/volatility or weighted-volatility/portfolio-volatility +ratio; transaction costs are reported and checked during executable repair, +but are not subtracted from these ratio definitions. CVaR/CDaR use exact empirical +tail-risk auxiliary variables; EVaR uses exponential cones. Their objectives are +forecast return minus risk aversion times scenario risk and trading costs. +CDaR uses cumulative simple-return paths, not compounded wealth drawdowns. + +All constraints enter the same problem, including trades out of the candidate +pool. Cash is retained in whole-account weights; benchmark weights are never +renormalized to the candidate pool or invested sleeve. `RiskSnapshot` carries +X/F/D or a dense covariance, availability time, horizon, annualization and model +identity. The factor path never constructs an asset-by-asset covariance. +HRP/risk-parity target construction currently requires explicit dense candidate +covariance; research configuration rejects their combination with factor-only +Barra inputs. The default exposure definitions are unchanged. + +CVXPY 1.7.5, OSQP 1.0.5 and Clarabel 0.11.1 are pinned. OSQP may retry the same +problem in Clarabel; an inaccurate result is never directly published. Failures +distinguish `infeasible`, `not_converged`, `numerical_failure`, `timeout`, +`no_positive_excess_return`, `undefined_ratio` and `repair_failed`. The legacy +mean–variance algorithm is available only through `solver_backend="legacy"` +for comparison; it is never an automatic fallback. Integer targets use exact +execution-layer cash quotes and deterministic local repair. They are verified +feasible targets, with no claim of a global integer optimum. + +Configurations without a protocol version retain their score interpretation +and are marked as compatibility mode. Solver correctness fixes can change their +weights. Exact historical reproduction requires the original frozen environment; +new runs receive new identities and must not overwrite previous evidence. + +Run `scripts/benchmark_optimizer.py --help` for the seeded 200/1,000/6,000 asset +factor benchmark. The regression suite contains the bounded-simplex KKT failure, +three six-asset projection cases, joint-turnover maximum Sharpe, pool-external +benchmark/cash tests and 200 independent SLSQP objective comparisons. + +## M6 governance and reproducibility (historical baseline) + +The original `0.4.2` governance baseline consumed these annotated internal tags. +The optimizer upgrade pins its reviewed protocol and execution revisions in +`pyproject.toml` and `requirements.lock` instead: - `quant-data-kit v0.8.1` (`b5621379e1a15562371be31c03f354a6acf512e4`) - `quant-execution v0.5.1` (`99ab9b1445d72164fa4e8c6d1ebde859b80dba1f`) diff --git a/pyproject.toml b/pyproject.toml index 40c86da..03a220b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,8 +10,11 @@ readme = "README.md" license = { text = "MIT" } requires-python = ">=3.10" dependencies = [ - "quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@f85128307dd0e3407d24350282fe7b679b591309", - "quant-execution @ git+https://github.com/PureSaber/quant-execution.git@49ef53d2c502ad6faad6818f853dadd1f8b8b182", + "cvxpy==1.7.5", + "osqp==1.0.5", + "clarabel==0.11.1", + "quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@ad3440700b226497c40093197f1df4321b49d1e8", + "quant-execution @ git+https://github.com/PureSaber/quant-execution.git@8423726a9d80b6509d2c016a7478fc0f40a29774", "pandas>=2.0", "pyarrow>=14.0", "pyyaml>=6.0", diff --git a/requirements.lock b/requirements.lock index 9c5b18a..5b74a97 100644 --- a/requirements.lock +++ b/requirements.lock @@ -1,105 +1,44 @@ -# -# This file is autogenerated by pip-compile with Python 3.10 -# by the following command: -# -# pip-compile --allow-unsafe --build-deps-for=editable --constraint=requirements-constraints.txt --extra=dev --index-url=https://pypi.org/simple --no-index --output-file=requirements.lock --pip-args='--timeout=15 --retries=1' --strip-extras pyproject.toml -# +# Exact dependency closure; optimizer additions resolved universally from solver-requirements.in. attrs==26.1.0 - # via - # jsonschema - # referencing +cffi==2.1.1 +clarabel==0.11.1 +cloudpickle==3.1.2 colorama==0.4.6 - # via pytest coverage==7.16.0 - # via pytest-cov +cvxpy==1.7.5 duckdb==1.5.5 - # via quant-data-kit exceptiongroup==1.3.1 - # via - # pytest - # quant-portfolio (pyproject.toml) iniconfig==2.3.0 - # via pytest +jinja2==3.1.6 +joblib==1.6.0 jsonschema==4.26.0 - # via - # quant-data-kit - # quant-execution jsonschema-rs==0.52.1 - # via quant-data-kit jsonschema-specifications==2025.9.1 - # via jsonschema +markupsafe==3.0.4 numpy==2.2.6 - # via - # pandas - # quant-data-kit - # quant-portfolio (pyproject.toml) orjson==3.12.0 - # via quant-data-kit +osqp==1.0.5 packaging==26.3 - # via pytest pandas==2.3.3 - # via - # quant-data-kit - # quant-portfolio (pyproject.toml) pluggy==1.6.0 - # via - # pytest - # pytest-cov pyarrow==25.0.1 - # via - # quant-data-kit - # quant-execution - # quant-portfolio (pyproject.toml) +pycparser==3.11 ; implementation_name != 'PyPy' pygments==2.21.0 - # via pytest pytest==9.1.1 - # via - # pytest-cov - # quant-portfolio (pyproject.toml) pytest-cov==7.1.0 - # via quant-portfolio (pyproject.toml) python-dateutil==2.9.0.post0 - # via pandas pytz==2026.3.post1 - # via pandas pyyaml==6.0.3 - # via - # quant-data-kit - # quant-portfolio (pyproject.toml) -quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@f85128307dd0e3407d24350282fe7b679b591309 - # via - # quant-execution - # quant-portfolio (pyproject.toml) -quant-execution @ git+https://github.com/PureSaber/quant-execution.git@49ef53d2c502ad6faad6818f853dadd1f8b8b182 - # via quant-portfolio (pyproject.toml) +quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@ad3440700b226497c40093197f1df4321b49d1e8 +quant-execution @ git+https://github.com/PureSaber/quant-execution.git@8423726a9d80b6509d2c016a7478fc0f40a29774 referencing==0.37.0 - # via - # jsonschema - # jsonschema-specifications rpds-py==0.30.0 - # via - # -c requirements-constraints.txt - # jsonschema - # referencing ruff==0.16.10 - # via quant-portfolio (pyproject.toml) +scipy==1.15.3 +scs==3.3.1 +setuptools==84.0.0 six==1.17.0 - # via python-dateutil tomli==2.4.1 ; python_version < "3.11" - # via - # coverage - # pytest - # quant-portfolio (pyproject.toml) typing-extensions==4.16.0 - # via - # exceptiongroup - # quant-data-kit - # referencing tzdata==2026.3 - # via pandas websockets==15.0.1 - # via quant-data-kit - -# The following packages are considered to be unsafe in a requirements file: -setuptools==84.0.0 - # via quant-portfolio (pyproject.toml::build-system.requires) diff --git a/scripts/benchmark_optimizer.py b/scripts/benchmark_optimizer.py new file mode 100644 index 0000000..cde130e --- /dev/null +++ b/scripts/benchmark_optimizer.py @@ -0,0 +1,96 @@ +"""Reproducible factor-only scalability probe (not a live trading latency claim).""" + +import argparse +import json +import os +from pathlib import Path +from time import perf_counter + +import numpy as np + +from quant_portfolio.optimization import OptimizationConstraints +from quant_portfolio.solver import RiskMatrix, solve_portfolio + + +def peak_memory(): + if os.name == "nt": + import ctypes + from ctypes import wintypes + + class Counters(ctypes.Structure): + _fields_ = [("cb", wintypes.DWORD), ("PageFaultCount", wintypes.DWORD)] + [ + (name, ctypes.c_size_t) + for name in ( + "PeakWorkingSetSize", + "WorkingSetSize", + "QuotaPeakPagedPoolUsage", + "QuotaPagedPoolUsage", + "QuotaPeakNonPagedPoolUsage", + "QuotaNonPagedPoolUsage", + "PagefileUsage", + "PeakPagefileUsage", + ) + ] + + counters = Counters() + counters.cb = ctypes.sizeof(counters) + api = ctypes.windll.psapi.GetProcessMemoryInfo + api.argtypes = [wintypes.HANDLE, ctypes.POINTER(Counters), wintypes.DWORD] + api.restype = wintypes.BOOL + if not api(wintypes.HANDLE(-1), ctypes.byref(counters), counters.cb): + raise ctypes.WinError() + return counters.PeakWorkingSetSize + import resource + + return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss * 1024 + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--assets", type=int, default=200) + parser.add_argument("--repeats", type=int, default=5) + parser.add_argument("--timeout", type=float, default=30) + parser.add_argument("--output", type=Path, required=True) + args = parser.parse_args() + n, k = args.assets, 25 + rng = np.random.default_rng(20261014 + n) + x = rng.normal(0, 0.4, (n, k)) + f = np.eye(k) * 0.003 + d = rng.uniform(0.01, 0.04, n) + risk = RiskMatrix(exposures=x, factor_covariance=f, specific_variance=d) + current = np.ones(n) / n + mu = rng.normal(0.001, 0.0001, n) + times, records = [], [] + for _ in range(args.repeats): + start = perf_counter() + got = solve_portfolio( + mu, + risk, + current, + OptimizationConstraints(max_weight=3 / n, max_turnover=0.4), + risk_scale=1 / 252, + risk_aversion=5, + linear_costs=np.ones(n) * 0.00005, + timeout_seconds=args.timeout, + ) + times.append(perf_counter() - start) + records.append(got.diagnostics) + payload = { + "assets": n, + "factors": k, + "seed": 20261014 + n, + "repeats": args.repeats, + "timeout_seconds": args.timeout, + "elapsed_seconds": times, + "p95_seconds": float(np.quantile(times, 0.95)), + "process_peak_memory_bytes": peak_memory(), + "asset_covariance_materialized": False, + "diagnostics": records, + } + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text(json.dumps(payload, indent=2), encoding="utf-8") + print(json.dumps({k: v for k, v in payload.items() if k != "diagnostics"})) + + +if __name__ == "__main__": + main() diff --git a/src/quant_portfolio/__init__.py b/src/quant_portfolio/__init__.py index 745eda3..dd362d2 100644 --- a/src/quant_portfolio/__init__.py +++ b/src/quant_portfolio/__init__.py @@ -1,5 +1,6 @@ """Sanitized quantitative research and portfolio engineering demos.""" +from quant_portfolio.allocation_v2 import AllocationResult, allocate_research from quant_portfolio.allocator import PortfolioSnapshot, allocate from quant_portfolio.cross_asset import ( CrossAssetConstraints, @@ -21,6 +22,7 @@ probability_of_backtest_overfitting, simulate_market, ) +from quant_portfolio.executable import ExecutableTarget, repair_integer_target from quant_portfolio.methods import ( black_litterman_expected_returns, constant_correlation_covariance, @@ -56,6 +58,8 @@ square_root_impact_cost, validate_research_allocation, ) +from quant_portfolio.risk_snapshot import RiskSnapshot +from quant_portfolio.solver import OptimizationError from quant_portfolio.synthetic_spread import ( BacktestMetrics, SyntheticSpreadConfig, @@ -64,24 +68,29 @@ ) __all__ = [ + "AllocationResult", "BacktestMetrics", "CrossAssetConstraints", "CrossAssetInput", "CrossAssetOptimizationResult", "DeflatedSharpe", + "ExecutableTarget", "MultiPeriodResult", "OptimizationConstraints", + "OptimizationError", "OptimizationFailure", "OptimizationResult", "OverfitResult", "PITFixedPoint", "PITMarketSnapshot", "PortfolioSnapshot", + "RiskSnapshot", "SimulationResult", "SyntheticSpreadConfig", "TargetPortfolio", "allocate", "allocate_integer_lots", + "allocate_research", "backtest_synthetic_spread", "black_litterman_expected_returns", "conditional_drawdown_at_risk", @@ -110,6 +119,7 @@ "optimize_multiperiod", "probability_of_backtest_overfitting", "random_matrix_covariance", + "repair_integer_target", "research_allocation_weights", "simulate_market", "square_root_impact_cost", diff --git a/src/quant_portfolio/allocation_v2.py b/src/quant_portfolio/allocation_v2.py new file mode 100644 index 0000000..b0762c5 --- /dev/null +++ b/src/quant_portfolio/allocation_v2.py @@ -0,0 +1,386 @@ +"""Unified research allocation in whole-account weights, including cash.""" + +from __future__ import annotations + +import importlib.metadata +from dataclasses import dataclass, field + +import numpy as np +import pandas as pd +from quant_data_kit.allocation_protocol import METHODS, validate_v2 + +from quant_portfolio.risk_snapshot import RiskSnapshot, mature_ic +from quant_portfolio.solver import RiskMatrix, solve_portfolio + + +@dataclass(frozen=True) +class AllocationResult: + weights: pd.Series + diagnostics: dict = field(default_factory=dict) + validation_context: dict = field(default_factory=dict, repr=False) + + def to_dict(self): + return {"continuous_weights": self.weights.to_dict(), "diagnostics": self.diagnostics} + + +def allocate_research( + scores, + trailing_returns, + current_weights, + linear_costs, + *, + invested_limit, + max_weight, + config, + factor_exposures=None, + factor_bounds=None, + covariance_override=None, + benchmark_weights=None, + factor_loadings=None, + factor_covariance=None, + specific_variance=None, + risk_snapshot=None, + decision_at=None, + expected_returns=None, + ic_labels=None, + lower_bounds=None, + upper_bounds=None, + buy_costs=None, + sell_costs=None, + impact_adv=None, + impact_volatility=None, + impact_coefficient=0.0, + portfolio_nav=None, + market_weights=None, + signal_cross_section=None, + max_tracking_error_override=None, + factor_bounds_absolute=False, +): + from quant_portfolio.optimization import ( + OptimizationConstraints, + _research_covariance, + _research_expected_returns, + validate_research_allocation, + ) + + v2 = config.get("protocol_version") == 2 + cfg = validate_v2(config) if v2 else validate_research_allocation(config) + if not isinstance(scores, pd.Series) or scores.empty or not scores.index.is_unique: + raise ValueError("scores must be a non-empty Series with unique assets") + scores = scores.astype(float) + if not np.isfinite(scores).all(): + raise ValueError("scores must be finite") + candidates = list(scores.index) + current = pd.Series(dtype=float) if current_weights is None else current_weights.astype(float) + if ( + not current.index.is_unique + or not np.isfinite(current).all() + or (current < 0).any() + or current.sum() > 1 + 1e-10 + ): + raise ValueError("current weights must be unique finite nonnegative account weights") + if not 0 < invested_limit <= 1 or not 0 < max_weight <= 1: + raise ValueError("invalid invested_limit or max_weight") + if risk_snapshot is not None: + if not isinstance(risk_snapshot, RiskSnapshot): + raise TypeError("risk_snapshot must be RiskSnapshot") + if decision_at is None: + raise ValueError("risk snapshot requires explicit decision_at") + if benchmark_weights is not None: + raise ValueError("benchmark must be provided only once in risk_snapshot") + benchmark_weights = risk_snapshot.benchmark + benchmark = None if benchmark_weights is None else benchmark_weights.astype(float) + if benchmark is not None and not benchmark.index.is_unique: + raise ValueError("benchmark labels must be unique") + assets = sorted( + set(candidates) + | set(current.index) + | (set() if benchmark is None else set(benchmark.index)) + ) + current = current.reindex(assets, fill_value=0.0) + wb = None if benchmark is None else benchmark.reindex(assets, fill_value=0.0).to_numpy() + lo = np.zeros(len(assets)) + hi = np.array([max_weight if a in candidates else 0.0 for a in assets]) + for bounds, destination in ((lower_bounds, lo), (upper_bounds, hi)): + if bounds is not None: + if not bounds.index.is_unique or set(bounds.index) - set(assets): + raise ValueError("invalid asset bound labels") + for i, asset in enumerate(assets): + if asset in bounds.index: + destination[i] = float(bounds.loc[asset]) + mode = cfg["mode"] + risk = None + history = None + if trailing_returns is not None: + history = trailing_returns.reindex(columns=assets).tail(cfg["lookback"]) + if v2: + from quant_portfolio.risk_snapshot import timestamp + + if decision_at is None or not isinstance(history.index, pd.DatetimeIndex): + raise ValueError("v2 history requires timestamped returns and decision_at") + # Index must denote availability, not a label's start date. + if any(timestamp(t) > timestamp(decision_at) for t in history.index): + raise ValueError("return history contains unavailable observations") + if risk_snapshot is not None: + risk = risk_snapshot.matrix(assets, decision_at) + elif covariance_override is not None: + if set(assets) - set(covariance_override.index): + raise ValueError("covariance does not cover the complete asset union") + risk = RiskMatrix(covariance=covariance_override.loc[assets, assets].to_numpy()) + elif factor_loadings is not None: + risk = RiskMatrix( + exposures=factor_loadings.loc[assets].to_numpy(), + factor_covariance=factor_covariance.loc[ + factor_loadings.columns, factor_loadings.columns + ].to_numpy(), + specific_variance=specific_variance.loc[assets].to_numpy(), + ) + elif mode != "equal" or cfg.get("max_tracking_error") is not None: + if ( + history is None + or len(history) < cfg["min_observations"] + or not np.isfinite(history).all().all() + ): + raise ValueError("complete return history for the asset union is insufficient") + if v2: + from quant_portfolio.optimization import estimate_covariance + + covariance = estimate_covariance(history, shrinkage=cfg["covariance_shrinkage"]) + else: + covariance = _research_covariance(cfg, history, None, assets, None, None, None) + risk = RiskMatrix(covariance=covariance.to_numpy()) + if ( + v2 + and risk_snapshot is not None + and risk_snapshot.horizon_sessions != cfg["forecast_horizon_sessions"] + ): + raise ValueError("risk prediction horizon differs from the forecast horizon") + horizon_scale = ( + ( + cfg["forecast_horizon_sessions"] + / (252 if risk_snapshot is None else risk_snapshot.annualization) + ) + if v2 + else 1.0 + ) + diag = { + "protocol_version": 2 if v2 else 1, + "compatibility_mode": not v2, + "risk_objective": cfg.get("risk_objective", "active" if wb is not None else "absolute"), + "risk_usage": METHODS[mode]["risk_usage"], + "asset_union": assets, + "invested_budget": invested_limit, + "cash_target": 1 - invested_limit, + "risk_horizon_scale": horizon_scale, + } + if risk_snapshot is not None: + diag["risk_snapshot"] = risk_snapshot.metadata() + mu = np.zeros(len(assets)) + if mode in {"cost_aware", "max_sharpe", "cvar", "evar", "cdar"}: + if v2: + if cfg["expected_return_model"] == "provided": + if expected_returns is None or not expected_returns.index.is_unique: + raise ValueError("provided return forecasts are required") + forecasts = expected_returns.reindex(candidates).to_numpy(dtype=float) + else: + if risk is None: + raise ValueError("ic_vol requires risk inputs") + ic, evidence = mature_ic( + ic_labels, + decision_at=decision_at, + horizon_sessions=cfg["forecast_horizon_sessions"], + lookback=cfg["ic_lookback"], + min_periods=cfg["ic_min_periods"], + ) + diag.update(evidence, information_coefficient=ic) + cross_section = ( + scores if signal_cross_section is None else signal_cross_section.astype(float) + ) + if ( + not cross_section.index.is_unique + or not np.isfinite(cross_section).all() + or set(candidates) - set(cross_section.index) + ): + raise ValueError( + "standardization universe must cover candidates with finite signals" + ) + raw = scores.to_numpy() + scale = float(cross_section.std(ddof=0)) + z = ( + np.zeros(len(scores)) + if scale == 0 + else (raw - float(cross_section.mean())) / scale + ) + diag["signal_standardization_assets"] = len(cross_section) + diag["constant_signal_cross_section"] = scale == 0 + vol = pd.Series(np.sqrt(risk.diagonal() * horizon_scale), index=assets) + forecasts = ic * vol.loc[candidates].to_numpy() * z + if not np.isfinite(forecasts).all(): + raise ValueError("return forecasts must cover all candidates") + mu = pd.Series(forecasts, index=candidates).reindex(assets, fill_value=0).to_numpy() + elif mode == "cost_aware" and cfg["expected_return_model"] != "score": + if cfg["expected_return_model"] == "black_litterman": + raise ValueError( + "Black-Litterman requires an explicit views API and full market weights" + ) + covariance = pd.DataFrame(np.diag(risk.diagonal()), index=assets, columns=assets) + mu = _research_expected_returns( + cfg, scores.reindex(assets, fill_value=0), covariance, market_weights + ).to_numpy() + else: + mu = scores.reindex(assets, fill_value=0).to_numpy() + x, bounds = None, dict(factor_bounds or {}) + if bounds: + if factor_exposures is None: + raise ValueError("factor bounds require exposures") + x = factor_exposures.reindex(index=assets, columns=list(bounds)).to_numpy() + if not np.isfinite(x).all(): + raise ValueError("constraint exposures must cover the asset union") + if cfg.get("factor_bound_reference") == "active" and not factor_bounds_absolute: + if wb is None: + raise ValueError("active factor bounds require full benchmark") + baseline = x.T @ wb + bounds = { + key: (pair[0] + baseline[i], pair[1] + baseline[i]) + for i, (key, pair) in enumerate(bounds.items()) + } + target = None + if METHODS[mode]["objective"] == "target": + if mode == "equal": + raw = pd.Series(1 / len(candidates), index=candidates) + elif mode == "inverse_vol": + vol = pd.Series(np.sqrt(risk.diagonal()), index=assets).loc[candidates] + if (vol <= 0).any(): + raise ValueError("inverse_vol requires positive variance") + raw = 1 / vol + raw /= raw.sum() + else: + from quant_portfolio.methods import ( + equal_risk_contribution_weights, + hierarchical_risk_parity_weights, + ) + + if risk.covariance is None: + raise ValueError( + f"{mode} target construction requires explicit dense candidate covariance" + ) + cov = pd.DataFrame(risk.covariance, index=assets, columns=assets).loc[ + candidates, candidates + ] + raw = ( + equal_risk_contribution_weights(cov) + if mode == "risk_parity" + else hierarchical_risk_parity_weights(cov) + ) + target = raw.reindex(assets, fill_value=0).to_numpy() * invested_limit + + def cost_vector(value): + if value is None: + return None + vector = value.reindex(assets) + # Liquidations need a cost estimate too; never silently assign zero. + needed = (current > 0) | pd.Series([a in candidates for a in assets], index=assets) + if vector[needed].isna().any(): + raise ValueError("cost rates must cover candidates and all held assets") + return vector.fillna(0).to_numpy(dtype=float) + + impact = None + if impact_coefficient or impact_adv is not None or impact_volatility is not None: + from quant_portfolio.optimization import _validated_impact + + impact = _validated_impact( + assets, impact_adv, impact_volatility, impact_coefficient, portfolio_nav + ) + tracking_limits = [ + value + for value in (cfg.get("max_tracking_error"), max_tracking_error_override) + if value is not None + ] + settings = OptimizationConstraints( + max_weight=max_weight, + max_turnover=cfg["max_turnover"], + max_tracking_error=min(tracking_limits) if tracking_limits else None, + ) + scenarios = None if history is None else history.to_numpy() + if v2 and mode in {"cvar", "evar", "cdar"}: + horizon = cfg["forecast_horizon_sessions"] + blocks = len(history) // horizon + if blocks < cfg["min_observations"]: + raise ValueError("insufficient non-overlapping horizon scenarios") + scenarios = ( + np.prod( + 1 + scenarios[-blocks * horizon :].reshape(blocks, horizon, len(assets)), axis=1 + ) + - 1 + ) + solved = solve_portfolio( + mu, + risk, + current.to_numpy(), + settings, + assets=assets, + budget=invested_limit, + lower=lo, + upper=hi, + factor_exposures=x, + factor_bounds=bounds, + benchmark=wb, + active_risk=diag["risk_objective"] == "active", + kind=METHODS[mode]["objective"], + target=target, + linear_costs=cost_vector(linear_costs), + buy_costs=cost_vector(buy_costs), + sell_costs=cost_vector(sell_costs), + risk_aversion=cfg["risk_aversion"], + risk_scale=horizon_scale, + turnover_penalty=cfg["turnover_penalty"], + impact=impact, + scenarios=scenarios, + beta=cfg.get(mode + "_beta", 0.95), + risk_free_rate=cfg.get("risk_free_rate", 0), + timeout_seconds=cfg.get("timeout_seconds", 30), + ) + w = solved.weights + diag.update( + solved.diagnostics, expected_return=float(mu @ w), turnover=float(np.abs(w - current).sum()) + ) + diag["versions"] = { + name: importlib.metadata.version(name) for name in ("cvxpy", "osqp", "clarabel", "numpy") + } + if risk is not None: + diag["annualized_volatility"] = float(np.sqrt(max(risk.variance(w), 0))) + if risk_snapshot is not None and risk_snapshot.exposures is not None: + names = list(risk_snapshot.exposures.columns) + diag["factor_exposures"] = dict(zip(names, (risk.exposures.T @ w).tolist())) + if wb is not None: + diag["active_exposures"] = dict(zip(names, (risk.exposures.T @ (w - wb)).tolist())) + if wb is not None: + diag["annualized_tracking_error"] = float(np.sqrt(max(risk.variance(w - wb), 0))) + if target is not None: + + def contributions(point): + marginal = ( + risk.covariance @ point + if risk.covariance is not None + else risk.exposures @ (risk.factor_covariance @ (risk.exposures.T @ point)) + + risk.specific_variance * point + ) + return (point * marginal).tolist() + + diag["target_risk_contributions"] = contributions(target) + diag["repaired_risk_contributions"] = contributions(w) + result = pd.Series(w, index=assets, name="weight") + result.attrs["allocation_diagnostics"] = diag + return AllocationResult( + result, + diag, + { + "assets": assets, + "risk": risk, + "benchmark": wb, + "settings": settings, + "factor_exposures": x, + "factor_bounds": bounds, + "current": current.to_numpy(), + }, + ) diff --git a/src/quant_portfolio/executable.py b/src/quant_portfolio/executable.py new file mode 100644 index 0000000..28ede05 --- /dev/null +++ b/src/quant_portfolio/executable.py @@ -0,0 +1,225 @@ +"""Deterministic local repair of continuous targets; never an integer optimality claim.""" + +from __future__ import annotations + +from dataclasses import dataclass +from decimal import Decimal +from time import perf_counter + +import numpy as np +import pandas as pd + +from quant_portfolio.solver import OptimizationError + + +@dataclass(frozen=True) +class ExecutableTarget: + quantities: dict[str, int] + weights: pd.Series + cash: Decimal + costs: Decimal + diagnostics: dict + + +def repair_integer_target( + weights, + *, + current, + prices, + lots, + cash, + nav, + quote_trade, + constraint_residuals, + max_invested, + max_buy=None, + max_sell=None, + timeout_seconds=30, + max_iterations=2000, +): + """quote_trade(asset, delta) returns exact (cash change, total friction). + + It is supplied by the execution adapter so rounding and fee rules have one + owner. Residuals must be independently computed in post-fee NAV weights. + Current quantities can contain odd-lot remainders; only board-lot trades + are proposed. No unobserved future market status is used. + """ + started = perf_counter() + if not np.isfinite(timeout_seconds) or timeout_seconds <= 0: + raise ValueError("repair timeout must be positive and finite") + if not np.isfinite(max_invested) or not 0 <= max_invested <= 1: + raise ValueError("max_invested must be in [0, 1]") + if type(max_iterations) is not int or max_iterations <= 0: + raise ValueError("max_iterations must be a positive integer") + assets = sorted(set(weights.index) | set(current)) + if not weights.index.is_unique or not np.isfinite(weights).all() or (weights < -1e-7).any(): + raise ValueError("invalid continuous weights") + nav, cash = Decimal(str(nav)), Decimal(str(cash)) + if not nav.is_finite() or nav <= 0 or not cash.is_finite() or cash < 0: + raise ValueError("invalid account cash or NAV") + n = len(assets) + if n == 0: + raise ValueError("at least one executable asset is required") + for mapping in (current, lots, max_buy or {}, max_sell or {}): + if any(isinstance(v, bool) or not np.isfinite(v) or int(v) != v for v in mapping.values()): + raise ValueError("quantities and board lots must be integers") + quantities = np.array([current.get(a, 0) for a in assets], dtype=np.int64) + steps = np.array([lots[a] for a in assets], dtype=np.int64) + if (quantities < 0).any() or (steps <= 0).any(): + raise ValueError("quantities and board lots must be non-negative integers") + marks = [Decimal(str(prices[a])) for a in assets] + if any(not p.is_finite() or p <= 0 for p in marks): + raise ValueError("positive exact prices are required for every asset") + target = weights.reindex(assets, fill_value=0).to_numpy() + sell = np.array( + [ + min(quantities[i], (max_sell or {}).get(a, int(quantities[i]))) + for i, a in enumerate(assets) + ] + ) + buy = np.array([(max_buy or {}).get(a, int(nav / marks[i])) for i, a in enumerate(assets)]) + if (sell < 0).any() or (buy < 0).any(): + raise ValueError("invalid tradable quantity") + minimum = quantities - (sell // steps) * steps + maximum = quantities + (buy // steps) * steps + quote_cache = {} + + def evaluate(q): + balance, costs = cash, Decimal(0) + for i, asset in enumerate(assets): + delta = int(q[i] - quantities[i]) + if delta: + key = (asset, delta) + if key not in quote_cache: + change, friction = quote_trade(asset, delta) + if not isinstance(change, Decimal) or not isinstance(friction, Decimal): + raise TypeError("execution quote must use exact Decimal amounts") + if not change.is_finite() or not friction.is_finite() or friction < 0: + raise ValueError("invalid execution quote") + quote_cache[key] = change, friction + change, friction = quote_cache[key] + balance += change + costs += friction + post_nav = nav - costs + if post_nav <= 0: + return (float("inf"), float("inf")), None + invested = sum((Decimal(int(q[i])) * marks[i] for i in range(n)), Decimal(0)) + # Ledger conservation is exact at the supplied marks and quote precision. + if balance + invested != post_nav: + raise ValueError("execution quote and account marks do not reconcile") + w = pd.Series( + [float(Decimal(int(q[i])) * marks[i] / post_nav) for i in range(n)], index=assets + ) + violations = constraint_residuals(w) + if not np.isfinite(list(violations.values())).all(): + raise ValueError("constraint residuals must be finite") + bound_error = max( + float(np.max(minimum - q, initial=0)), float(np.max(q - maximum, initial=0)) + ) + residual = max( + [ + 0.0, + max(1e-6, -float(balance / nav)) if balance < 0 else 0.0, + float(invested / nav) - max_invested, + bound_error, + *(max(float(v), 0) for v in violations.values()), + ] + ) + objective = float(np.sum((w.to_numpy() - target) ** 2)) + return (residual if residual > 1e-7 else 0.0, objective), (w, balance, costs, violations) + + desired = np.array([int(nav * Decimal(str(max(target[i], 0))) / marks[i]) for i in range(n)]) + # Preserve the account's odd-lot remainder and enforce known quantity limits. + floor = quantities + np.floor_divide(desired - quantities, steps) * steps + floor = np.minimum(np.maximum(floor, minimum), maximum) + starts = [floor, quantities.copy()] + best = None + attempts = [] + for seed in starts: + actions = [] + q = seed.copy() + score, details = evaluate(q) + seen = set() + for iteration in range(max_iterations): + if perf_counter() - started > timeout_seconds: + break + key = tuple(q) + if key in seen: + break + seen.add(key) + winner, winner_score, winner_details, move = q, score, details, None + + # Single reductions, additions and then single-lot exchanges. Sorting + # asset identifiers and strict improvement make ties deterministic. + def candidates(q): + for i in range(n): + for direction in (-1, 1): + candidate = q.copy() + candidate[i] += direction * steps[i] + if minimum[i] <= candidate[i] <= maximum[i]: + yield candidate, [(assets[i], direction * int(steps[i]))] + for i in range(n): + if q[i] - steps[i] < minimum[i]: + continue + for j in range(n): + if i == j or q[j] + steps[j] > maximum[j]: + continue + candidate = q.copy() + candidate[i] -= steps[i] + candidate[j] += steps[j] + yield candidate, [(assets[i], -int(steps[i])), (assets[j], int(steps[j]))] + + for candidate, changes in candidates(q): + if perf_counter() - started > timeout_seconds: + break + candidate_score, candidate_details = evaluate(candidate) + if candidate_score < winner_score: + winner, winner_score, winner_details, move = ( + candidate, + candidate_score, + candidate_details, + changes, + ) + if move is None: + break + q, score, details = winner, winner_score, winner_details + actions.append(move) + if score[0] == 0 and (best is None or score < best[0]): + best = (score, q.copy(), details, list(actions)) + attempts.append({"seed": seed.tolist(), "actions": actions, "residual": score[0]}) + if perf_counter() - started > timeout_seconds: + raise OptimizationError( + "timeout", "integer repair deadline exceeded", {"repair_attempts": attempts} + ) + if best is None: + raise OptimizationError( + "repair_failed", + "local repair found no verified executable target", + {"repair_attempts": attempts, "seconds": perf_counter() - started}, + ) + score, q, details, actions = best + final_score, details = evaluate(q) + w, balance, costs, residuals = details + if final_score[0] != 0 or balance < 0: + raise OptimizationError("repair_failed", "final account/risk verification failed") + return ExecutableTarget( + {asset: int(q[i]) for i, asset in enumerate(assets) if q[i]}, + w, + balance, + costs, + { + "status": "repaired", + "integer_global_optimum": False, + "repair_actions": actions, + "repair_attempts": attempts, + "constraint_residuals": residuals, + "seconds": perf_counter() - started, + "cash": str(balance), + "costs": str(costs), + "post_fee_nav": str(nav - costs), + "budget_shortfall": float(max_invested) + - float(sum(Decimal(int(q[i])) * marks[i] for i in range(n)) / nav), + "shortfall_reasons": ["fees", "board_lot_granularity", "known_quantity_constraints"], + "integer_weights": w.to_dict(), + }, + ) diff --git a/src/quant_portfolio/objectives.py b/src/quant_portfolio/objectives.py index 215a3fc..c8b6004 100644 --- a/src/quant_portfolio/objectives.py +++ b/src/quant_portfolio/objectives.py @@ -1,10 +1,8 @@ """Independent portfolio objectives that share the long-only feasible set. -Mean-variance stays in ``optimization.py``. The solvers here are projected -gradient methods on that same set: maximum Sharpe, maximum diversification, -entropic value at risk, and conditional drawdown at risk. Integer share lots -are a separate exact local search, not a continuous relaxation that is rounded -only when an order is created. +Public objective adapters use the shared CVXPY solver. Ratio feasibility and +path-risk epigraphs include the real constraints. Integer lots use deterministic +local repair and carry no global integer optimality claim. """ from __future__ import annotations @@ -23,13 +21,12 @@ def maximum_drawdown(portfolio_returns: np.ndarray) -> float: def conditional_drawdown_at_risk(portfolio_returns: np.ndarray, beta: float) -> float: - """Mean drawdown at and above the empirical ``beta`` quantile.""" + """Empirical RU tail drawdown, including a fractional boundary observation.""" _validate_beta(beta) drawdowns = _drawdown_series(portfolio_returns) - level = float(np.quantile(drawdowns, float(beta))) - tail = drawdowns[drawdowns >= level - 1e-15] - return float(tail.mean()) if len(tail) else float(drawdowns.max()) + level = float(np.quantile(drawdowns, float(beta), method="higher")) + return level + float(np.maximum(drawdowns - level, 0).mean()) / (1 - float(beta)) def entropic_value_at_risk(losses: np.ndarray, beta: float) -> float: @@ -110,7 +107,7 @@ def allocate_integer_lots( counts = np.maximum(np.floor(point / grid + 1e-10).astype(int), 0) counts = _shed_infeasible_lots(counts, grid, objective, is_feasible) if not is_feasible(counts * grid): - raise ValueError("no feasible integer-lot portfolio exists under the current constraints") + raise ValueError("repair_failed: no verified integer-lot target was found") return _improve_lots(counts, grid, objective, is_feasible) @@ -136,25 +133,10 @@ def optimize_max_sharpe( ) risk_free = _finite("risk_free_rate", risk_free_rate, positive=False) - def ratio(weights: np.ndarray) -> float: - volatility = math.sqrt(max(float(weights @ problem.covariance @ weights), 0.0)) - if volatility <= 1e-18: - raise ValueError("Sharpe ratio is undefined for a zero-volatility portfolio") - return (float(problem.expected @ weights) - risk_free) / volatility - - def gradient(weights: np.ndarray) -> np.ndarray: - excess = float(problem.expected @ weights) - risk_free - variance = max(float(weights @ problem.covariance @ weights), 1e-18) - volatility = math.sqrt(variance) - risk_gradient = problem.covariance @ weights - return problem.expected / volatility - excess * risk_gradient / variance**1.5 - - start = _positive_projection(problem.covariance, problem.expected - risk_free) - return _maximize_ratio( + return _convex_result( problem, - ratio, - gradient, - start, + "max_sharpe", + risk_free_rate=risk_free, max_iterations=max_iterations, tolerance=tolerance, ) @@ -184,26 +166,8 @@ def optimize_max_diversification( if np.any(volatility <= 0): raise ValueError("maximum diversification requires positive asset variance") - def ratio(weights: np.ndarray) -> float: - denominator = math.sqrt(max(float(weights @ problem.covariance @ weights), 0.0)) - if denominator <= 1e-18: - raise ValueError("diversification ratio is undefined for a zero-volatility portfolio") - return float(volatility @ weights) / denominator - - def gradient(weights: np.ndarray) -> np.ndarray: - numerator = float(volatility @ weights) - variance = max(float(weights @ problem.covariance @ weights), 1e-18) - denominator = math.sqrt(variance) - return volatility / denominator - numerator * (problem.covariance @ weights) / variance**1.5 - - start = _positive_projection(problem.covariance, volatility) - return _maximize_ratio( - problem, - ratio, - gradient, - start, - max_iterations=max_iterations, - tolerance=tolerance, + return _convex_result( + problem, "max_diversification", max_iterations=max_iterations, tolerance=tolerance ) @@ -221,36 +185,18 @@ def optimize_evar( ): """Maximize expected return minus entropic value at risk.""" - scenarios, expected, problem = _path_problem( + scenarios, _expected, problem = _path_problem( scenario_returns, expected_returns, current_weights, constraints ) aversion = _finite("risk_aversion", risk_aversion, positive=True) _validate_beta(beta) - anchor = _drawdown_anchor(scenarios, problem.project, max_drawdown) - - def objective_at(weights: np.ndarray) -> float: - losses = -(scenarios @ weights) - return float(expected @ weights) - aversion * entropic_value_at_risk(losses, beta) - - def gradient(weights: np.ndarray) -> np.ndarray: - losses = -(scenarios @ weights) - _value, temperature = _evar_temperature(losses, beta) - shifted = losses / temperature - shifted -= float(np.max(shifted)) - probabilities = np.exp(shifted) - probabilities /= probabilities.sum() - return expected + aversion * (scenarios.T @ probabilities) - - return _maximize_path( + return _convex_result( problem, - scenarios, - objective_at, - gradient, - anchor, - max_drawdown, - risk_name="evar", - risk_value=lambda weights: entropic_value_at_risk(-(scenarios @ weights), beta), - aversion=aversion, + "evar", + scenarios=scenarios, + beta=beta, + risk_aversion=aversion, + max_drawdown=max_drawdown, max_iterations=max_iterations, tolerance=tolerance, ) @@ -274,44 +220,57 @@ def optimize_cdar( Scenario order is the path order. """ - scenarios, expected, problem = _path_problem( + scenarios, _expected, problem = _path_problem( scenario_returns, expected_returns, current_weights, constraints ) aversion = _finite("risk_aversion", risk_aversion, positive=True) _validate_beta(beta) - anchor = _drawdown_anchor(scenarios, problem.project, max_drawdown) - - def objective_at(weights: np.ndarray) -> float: - return float(expected @ weights) - aversion * conditional_drawdown_at_risk( - scenarios @ weights, beta - ) - - def gradient(weights: np.ndarray) -> np.ndarray: - path = scenarios @ weights - drawdowns = _drawdown_series(path) - level = float(np.quantile(drawdowns, float(beta))) - tail = drawdowns >= level - 1e-15 - if not tail.any(): - tail = np.zeros(len(drawdowns), dtype=bool) - tail[int(np.argmax(drawdowns))] = True - risk_gradient = _tail_drawdown_gradient(path, scenarios, tail) - return expected - aversion * risk_gradient - - return _maximize_path( + return _convex_result( problem, - scenarios, - objective_at, - gradient, - anchor, - max_drawdown, - risk_name="cdar", - risk_value=lambda weights: conditional_drawdown_at_risk(scenarios @ weights, beta), - aversion=aversion, + "cdar", + scenarios=scenarios, + beta=beta, + risk_aversion=aversion, + max_drawdown=max_drawdown, max_iterations=max_iterations, tolerance=tolerance, ) +def _convex_result(problem, kind, *, tolerance, **kwargs): + from quant_portfolio.solver import RiskMatrix, solve_portfolio + + _finite("tolerance", tolerance, positive=True) + solved = solve_portfolio( + problem.expected, + RiskMatrix(covariance=problem.covariance), + problem.current, + problem.settings, + assets=problem.assets, + benchmark=problem.benchmark, + kind=kind, + **kwargs, + ) + w = solved.weights + variance = float(w @ problem.covariance @ w) + volatility = ( + float(np.std(kwargs["scenarios"] @ w, ddof=1)) + if "scenarios" in kwargs + else math.sqrt(max(variance, 0)) + ) + return _solver_api()[1]( + weights=pd.Series(w, index=problem.assets, name="weight"), + expected_return=float(problem.expected @ w), + volatility=volatility, + turnover=float(np.abs(w - problem.current).sum()), + objective=solved.objective, + group_weights=_group_weights(problem.assets, w, problem.settings), + converged=True, + iterations=solved.diagnostics["iterations"], + diagnostics=solved.diagnostics, + ) + + class _LongOnlyInputs: def __init__( self, @@ -430,73 +389,6 @@ def _aligned_weights( return numeric -def _positive_projection(covariance: np.ndarray, direction: np.ndarray) -> np.ndarray: - stabilized = covariance + np.eye(covariance.shape[0]) * 1e-10 - raw = np.linalg.solve(stabilized, np.maximum(direction, 0.0)) - if not np.isfinite(raw).all() or float(raw.sum()) <= 0: - raw = np.maximum(direction, 0.0) - if float(raw.sum()) <= 0: - raw = np.ones(covariance.shape[0]) - return raw / raw.sum() - - -def _maximize_ratio(problem, ratio, gradient, start, *, max_iterations: int, tolerance: float): - _result_type = _solver_api()[1] - _iteration_limit(max_iterations) - tolerance = _finite("tolerance", tolerance, positive=True) - candidates = ( - problem.project(np.full(len(problem.assets), 1.0 / len(problem.assets))), - problem.project(start), - ) - best = candidates[0] - best_value = ratio(best) - for candidate in candidates[1:]: - value = ratio(candidate) - if value > best_value: - best = candidate - best_value = value - weights = best.copy() - converged = False - iteration = max_iterations - for iteration in range(1, max_iterations + 1): - step = 0.25 - current_value = ratio(weights) - slope = gradient(weights) - candidate = weights - for _ in range(20): - proposal = problem.project(weights + step * slope) - proposal_value = ratio(proposal) - if proposal_value >= current_value - 1e-15: - candidate = proposal - break - step *= 0.5 - else: - converged = True - break - if proposal_value > best_value + tolerance: - best = candidate.copy() - best_value = proposal_value - if float(np.max(np.abs(candidate - weights))) <= tolerance: - weights = candidate - converged = True - break - weights = candidate - weights = best - if not converged and iteration >= max_iterations: - raise RuntimeError(f"ratio optimization did not converge after {max_iterations} iterations") - variance = max(float(weights @ problem.covariance @ weights), 0.0) - return _result_type( - weights=pd.Series(weights, index=problem.assets, name="weight"), - expected_return=float(problem.expected @ weights), - volatility=float(math.sqrt(variance)), - turnover=float(np.abs(weights - problem.current).sum()), - objective=float(ratio(weights)), - group_weights=_group_weights(problem.assets, weights, problem.settings), - converged=True, - iterations=iteration, - ) - - def _path_problem(scenario_returns, expected_returns, current_weights, constraints): if getattr(constraints, "max_tracking_error", None) is not None: raise ValueError( @@ -526,89 +418,6 @@ def _path_problem(scenario_returns, expected_returns, current_weights, constrain return values, expected, problem -def _drawdown_anchor(scenarios, project, limit: float | None): - if limit is None: - return None - if ( - isinstance(limit, bool) - or not isinstance(limit, (int, float)) - or not np.isfinite(limit) - or limit < 0 - ): - raise ValueError("max_drawdown must be finite and non-negative") - anchor = minimum_drawdown_weights(scenarios, project) - if maximum_drawdown(scenarios @ anchor) > float(limit) + 1e-8: - raise ValueError("max_drawdown is infeasible for these scenarios and constraints") - return anchor - - -def _maximize_path( - problem, - scenarios, - objective_at, - gradient, - anchor, - max_drawdown_limit, - *, - risk_name: str, - risk_value, - aversion: float, - max_iterations: int, - tolerance: float, -): - del risk_name - result_type = _solver_api()[1] - _iteration_limit(max_iterations) - tolerance = _finite("tolerance", tolerance, positive=True) - weights = problem.project(problem.current) - best = weights.copy() - best_value = -np.inf - converged = False - iteration = max_iterations - quiet = 0 - for iteration in range(1, max_iterations + 1): - slope = gradient(weights) - scale = max(float(np.max(np.abs(slope))), 1e-8) - step = 0.25 / scale / (1.0 + 0.02 * iteration) - candidate = problem.project(weights + step * slope) - if anchor is not None: - candidate = enforce_drawdown_limit( - candidate, scenarios, float(max_drawdown_limit), anchor - ) - value = objective_at(candidate) - if value > best_value + tolerance: - best = candidate.copy() - best_value = value - quiet = 0 - else: - quiet += 1 - if float(np.max(np.abs(candidate - weights))) <= tolerance or quiet >= 40: - weights = best - converged = True - break - weights = candidate - if not converged: - raise RuntimeError(f"path optimization did not converge after {max_iterations} iterations") - weights = best - if ( - anchor is not None - and maximum_drawdown(scenarios @ weights) > float(max_drawdown_limit) + 1e-8 - ): - raise RuntimeError("optimizer produced a portfolio outside the drawdown limit") - risk = float(risk_value(weights)) - expected = float(problem.expected @ weights) - return result_type( - weights=pd.Series(weights, index=problem.assets, name="weight"), - expected_return=expected, - volatility=float(np.std(scenarios @ weights, ddof=1)), - turnover=float(np.abs(weights - problem.current).sum()), - objective=expected - aversion * risk, - group_weights=_group_weights(problem.assets, weights, problem.settings), - converged=True, - iterations=iteration, - ) - - def _shed_infeasible_lots(counts, steps, objective, is_feasible): guard = int(counts.sum()) + 1 while not is_feasible(counts * steps) and counts.sum() > 0 and guard > 0: diff --git a/src/quant_portfolio/optimization.py b/src/quant_portfolio/optimization.py index b74d66c..6423b29 100644 --- a/src/quant_portfolio/optimization.py +++ b/src/quant_portfolio/optimization.py @@ -24,15 +24,30 @@ ) from quant_portfolio.objectives import ( allocate_integer_lots, - enforce_drawdown_limit, - maximum_drawdown, - minimum_drawdown_weights, optimize_cdar, optimize_evar, optimize_max_diversification, optimize_max_sharpe, ) +__all__ = [ + "MultiPeriodResult", + "OptimizationConstraints", + "OptimizationResult", + "estimate_capacity", + "estimate_covariance", + "optimize_cdar", + "optimize_cvar", + "optimize_evar", + "optimize_max_diversification", + "optimize_max_sharpe", + "optimize_mean_variance", + "optimize_multiperiod", + "research_allocation_weights", + "square_root_impact_cost", + "validate_research_allocation", +] + @dataclass(frozen=True) class OptimizationConstraints: @@ -54,6 +69,7 @@ class OptimizationResult: group_weights: dict[str, float] converged: bool iterations: int + diagnostics: dict[str, object] = field(default_factory=dict) _RESEARCH_ALLOCATION_FIELDS = { @@ -108,148 +124,9 @@ class OptimizationResult: def validate_research_allocation(value: Mapping[str, object]) -> dict[str, object]: - """Validate and normalize the closed research-allocation recipe fragment.""" - - if not isinstance(value, Mapping): - raise TypeError("allocation must be a mapping") - unknown = set(value) - _RESEARCH_ALLOCATION_FIELDS - if unknown: - raise ValueError(f"allocation contains unknown fields: {sorted(unknown)}") - mode = value.get("mode") - if mode not in _ALLOCATION_MODES: - raise ValueError( - "allocation.mode must be equal, inverse_vol, cost_aware, risk_parity, hrp or cvar, " - "max_sharpe, max_diversification, evar or cdar" - ) - - def integer(name: str, default: int) -> int: - result = value.get(name, default) - if isinstance(result, bool) or not isinstance(result, int): - raise TypeError(f"allocation.{name} must be an integer") - return result + from quant_data_kit.allocation_protocol import validate_research_allocation as validate - def number(name: str, default: float) -> float: - result = value.get(name, default) - if isinstance(result, bool) or not isinstance(result, (int, float)): - raise TypeError(f"allocation.{name} must be numeric") - result = float(result) - if not np.isfinite(result): - raise ValueError(f"allocation.{name} must be finite") - return result - - lookback = integer("lookback", 20) - min_observations = integer("min_observations", 10) - shrinkage = number("covariance_shrinkage", 0.2) - risk_aversion = number("risk_aversion", 5.0) - turnover_penalty = number("turnover_penalty", 0.0) - max_turnover = number("max_turnover", 1.0) - if lookback < 2: - raise ValueError("allocation.lookback must be at least 2") - if not 2 <= min_observations <= lookback: - raise ValueError("allocation.min_observations must be between 2 and lookback") - if not 0 <= shrinkage <= 1: - raise ValueError("allocation.covariance_shrinkage must be in [0, 1]") - if risk_aversion <= 0: - raise ValueError("allocation.risk_aversion must be positive") - if turnover_penalty < 0: - raise ValueError("allocation.turnover_penalty must be non-negative") - if not 0 <= max_turnover <= 2: - raise ValueError("allocation.max_turnover must be in [0, 2]") - estimator = value.get("covariance_estimator", "diagonal") - if estimator not in _COVARIANCE_ESTIMATORS: - raise ValueError( - "allocation.covariance_estimator must be diagonal, ledoit_wolf, factor, ewma, " - "constant_correlation, oas or random_matrix" - ) - if mode not in _COVARIANCE_MODES and estimator != "diagonal": - raise ValueError( - "allocation.covariance_estimator is only used by cost_aware, risk_parity, hrp, " - "max_sharpe and max_diversification" - ) - if "ewma_decay" in value and estimator != "ewma": - raise ValueError("allocation.ewma_decay is only supported for the ewma estimator") - decay = number("ewma_decay", 0.94) - if estimator == "ewma" and not 0 < decay < 1: - raise ValueError("allocation.ewma_decay must be in (0, 1)") - if "risk_free_rate" in value and mode != "max_sharpe": - raise ValueError("allocation.risk_free_rate is only supported for max_sharpe") - risk_free_rate = number("risk_free_rate", 0.0) - if "max_tracking_error" in value and mode != "cost_aware": - raise ValueError("allocation.max_tracking_error is only supported for cost_aware") - raw_tracking_error = value.get("max_tracking_error") - if raw_tracking_error is None: - max_tracking_error = None - else: - max_tracking_error = number("max_tracking_error", 0.0) - if max_tracking_error < 0: - raise ValueError("allocation.max_tracking_error must be non-negative") - if "evar_beta" in value and mode != "evar": - raise ValueError("allocation.evar_beta is only supported for evar") - if "cdar_beta" in value and mode != "cdar": - raise ValueError("allocation.cdar_beta is only supported for cdar") - expected_return_model = value.get("expected_return_model", "score") - if expected_return_model not in {"score", "ic_vol", "black_litterman"}: - raise ValueError( - "allocation.expected_return_model must be score, ic_vol or black_litterman" - ) - if mode != "cost_aware" and ( - expected_return_model != "score" - or "information_coefficient" in value - or "black_litterman_tau" in value - ): - raise ValueError("expected-return fields are only supported for cost_aware") - information_coefficient = number("information_coefficient", 0.05) - tau = number("black_litterman_tau", 0.05) - if ( - expected_return_model in {"ic_vol", "black_litterman"} - and not 0 < information_coefficient <= 1 - ): - raise ValueError("allocation.information_coefficient must be in (0, 1]") - if expected_return_model == "black_litterman" and not 0 < tau <= 1: - raise ValueError("allocation.black_litterman_tau must be in (0, 1]") - if mode != "cvar" and "cvar_beta" in value: - raise ValueError("allocation.cvar_beta is only supported for cvar") - cvar_beta = number("cvar_beta", 0.95) - if mode == "cvar" and not 0 < cvar_beta < 1: - raise ValueError("allocation.cvar_beta must be in (0, 1)") - evar_beta = number("evar_beta", 0.95) - cdar_beta = number("cdar_beta", 0.95) - if mode == "evar" and not 0 < evar_beta < 1: - raise ValueError("allocation.evar_beta must be in (0, 1)") - if mode == "cdar" and not 0 < cdar_beta < 1: - raise ValueError("allocation.cdar_beta must be in (0, 1)") - factor_bound_reference = value.get("factor_bound_reference", "absolute") - if factor_bound_reference not in {"absolute", "active"}: - raise ValueError("allocation.factor_bound_reference must be absolute or active") - if mode != "cost_aware" and factor_bound_reference != "absolute": - raise ValueError("allocation.factor_bound_reference is only supported for cost_aware") - normalized = { - "mode": mode, - "lookback": lookback, - "min_observations": min_observations, - "covariance_shrinkage": shrinkage, - "covariance_estimator": estimator, - "risk_aversion": risk_aversion, - "turnover_penalty": turnover_penalty, - "max_turnover": max_turnover, - "expected_return_model": expected_return_model, - "factor_bound_reference": factor_bound_reference, - } - # Persist only applicable optional fields. The execution adapter validates this - # normalized fragment again after serializing its plan. - if estimator == "ewma": - normalized["ewma_decay"] = decay - if mode == "cost_aware": - normalized.update( - information_coefficient=information_coefficient, - black_litterman_tau=tau, - max_tracking_error=max_tracking_error, - ) - if mode == "max_sharpe": - normalized["risk_free_rate"] = risk_free_rate - if mode in _PATH_MODES: - normalized[f"{mode}_beta"] = {"cvar": cvar_beta, "evar": evar_beta, "cdar": cdar_beta}[mode] - return normalized + return validate(value) def _validated_covariance(covariance: pd.DataFrame, assets: list[str]) -> np.ndarray: @@ -364,6 +241,15 @@ def research_allocation_weights( impact_volatility: pd.Series | None = None, impact_coefficient: float = 0.0, portfolio_nav: float | None = None, + risk_snapshot=None, + decision_at=None, + expected_returns=None, + ic_labels=None, + lower_bounds=None, + upper_bounds=None, + buy_costs=None, + sell_costs=None, + max_tracking_error_override=None, ) -> pd.Series: """Build long-only research weights with one explicit, closed allocation policy. @@ -381,6 +267,45 @@ def research_allocation_weights( liquidation of current assets absent from ``scores``. """ + if ( + config.get("protocol_version") == 2 + or risk_snapshot is not None + or config.get("mode") in {"max_sharpe", "max_diversification", "cvar", "evar", "cdar"} + or lower_bounds is not None + or upper_bounds is not None + ): + from quant_portfolio.allocation_v2 import allocate_research + + return allocate_research( + scores, + trailing_returns, + current_weights, + linear_costs, + invested_limit=invested_limit, + max_weight=max_weight, + config=config, + factor_exposures=factor_exposures, + factor_bounds=factor_bounds, + covariance_override=covariance_override, + benchmark_weights=benchmark_weights, + factor_loadings=factor_loadings, + factor_covariance=factor_covariance, + specific_variance=specific_variance, + risk_snapshot=risk_snapshot, + decision_at=decision_at, + expected_returns=expected_returns, + ic_labels=ic_labels, + lower_bounds=lower_bounds, + upper_bounds=upper_bounds, + buy_costs=buy_costs, + sell_costs=sell_costs, + impact_adv=impact_adv, + impact_volatility=impact_volatility, + impact_coefficient=impact_coefficient, + portfolio_nav=portfolio_nav, + market_weights=market_weights, + max_tracking_error_override=max_tracking_error_override, + ).weights settings = validate_research_allocation(config) if not isinstance(scores, pd.Series) or scores.empty or not scores.index.is_unique: raise ValueError("scores must be a non-empty Series with unique assets") @@ -440,30 +365,22 @@ def research_allocation_weights( exposures, bounds = _validated_factor_constraints(assets, factor_exposures, factor_bounds) desired = np.repeat(total / len(assets), len(assets)) if exposures is not None: - # The joint projector uses a unit-budget sleeve. All absolute - # portfolio limits, including liquidation outside the sleeve, scale - # by the same invested fraction; cash has zero factor exposure. - constraints = _validated_constraints( - OptimizationConstraints( - max_weight=min(1.0, cap / total), - max_turnover=turnover_limit / total, - ), - assets, - ) - weights = _project_joint_feasible_set( - desired / total, - assets, - current.to_numpy(dtype=float) / total, - outside_turnover / total, - constraints, - exposures, - { - factor: (lower / total, upper / total) - for factor, (lower, upper) in bounds.items() - }, - tolerance=1e-10, + from quant_portfolio.solver import solve_portfolio + + solved = solve_portfolio( + np.zeros(len(assets)), + None, + current.to_numpy(dtype=float), + OptimizationConstraints(max_weight=cap, max_turnover=turnover_limit), + assets=assets, + budget=total, + kind="target", + target=desired, + factor_exposures=exposures, + factor_bounds=bounds, + turnover_offset=outside_turnover, ) - return pd.Series(weights * total, index=assets, name="weight") + return pd.Series(solved.weights, index=assets, name="weight") weights = _constrain_turnover( desired, current.to_numpy(dtype=float), @@ -515,73 +432,22 @@ def research_allocation_weights( ) return pd.Series(weights, index=assets, name="weight") - if mode in (_COVARIANCE_MODES - {"cost_aware"}) | _PATH_MODES: - sleeve_cap = min(1.0, cap / total) - if mode in _PATH_MODES: - assert window is not None - path_constraints = OptimizationConstraints(max_weight=sleeve_cap, max_turnover=2.0) - if mode == "cvar": - solved = optimize_cvar( - window, - numeric_scores, - beta=float(settings["cvar_beta"]), - risk_aversion=float(settings["risk_aversion"]), - constraints=path_constraints, - ) - elif mode == "evar": - solved = optimize_evar( - window, - numeric_scores, - beta=float(settings["evar_beta"]), - risk_aversion=float(settings["risk_aversion"]), - constraints=path_constraints, - ) - else: - solved = optimize_cdar( - window, - numeric_scores, - beta=float(settings["cdar_beta"]), - risk_aversion=float(settings["risk_aversion"]), - constraints=path_constraints, - ) - desired = solved.weights.to_numpy(dtype=float) * total - elif mode in {"max_sharpe", "max_diversification"}: - covariance = _research_covariance( - settings, - window, - covariance_override, - assets, - factor_loadings, - factor_covariance, - specific_variance, - ) - ratio_constraints = OptimizationConstraints(max_weight=sleeve_cap, max_turnover=2.0) - if mode == "max_sharpe": - solved = optimize_max_sharpe( - numeric_scores, - covariance, - risk_free_rate=float(settings["risk_free_rate"]), - constraints=ratio_constraints, - ) - else: - solved = optimize_max_diversification(covariance, constraints=ratio_constraints) - desired = solved.weights.to_numpy(dtype=float) * total - else: - covariance = _research_covariance( - settings, - window, - covariance_override, - assets, - factor_loadings, - factor_covariance, - specific_variance, - ) - raw = ( - equal_risk_contribution_weights(covariance) - if mode == "risk_parity" - else hierarchical_risk_parity_weights(covariance) - ) - desired = raw.reindex(assets).to_numpy(dtype=float) * total + if mode in {"risk_parity", "hrp"}: + covariance = _research_covariance( + settings, + window, + covariance_override, + assets, + factor_loadings, + factor_covariance, + specific_variance, + ) + raw = ( + equal_risk_contribution_weights(covariance) + if mode == "risk_parity" + else hierarchical_risk_parity_weights(covariance) + ) + desired = raw.reindex(assets).to_numpy(dtype=float) * total desired = _project_box_simplex(desired, 0.0, cap, total=total) weights = _constrain_turnover( desired, @@ -781,25 +647,6 @@ def _project_box_simplex( return np.repeat(lower, count) if abs(upper * count - total) <= 1e-14: return np.repeat(upper, count) - projected = np.zeros(count, dtype=float) - free = np.ones(count, dtype=bool) - for _ in range(count + 1): - remaining = total - float(projected[~free].sum()) - shift = (float(values[free].sum()) - remaining) / int(free.sum()) - candidates = values[free] - shift - below = candidates < lower - above = candidates > upper - if not below.any() and not above.any(): - projected[free] = candidates - return projected - free_indices = np.flatnonzero(free) - projected[free_indices[below]] = lower - projected[free_indices[above]] = upper - free[free_indices[below | above]] = False - if not free.any(): - if abs(float(projected.sum()) - total) <= 1e-12: - return projected - break low = float(np.min(values - upper)) high = float(np.max(values - lower)) for _ in range(80): @@ -826,28 +673,27 @@ def _constrain_turnover( if max_turnover < 0 or turnover_offset < 0: raise ValueError("turnover limits and offsets must be non-negative") - anchor = _project_box_simplex(current, lower, upper, total=total) + from quant_portfolio.solver import OptimizationError, solve_portfolio - def turnover(target: np.ndarray) -> float: - return float(np.abs(target - current).sum()) + turnover_offset - - minimum = turnover(anchor) - if minimum > max_turnover + 1e-10: - raise ValueError( - "allocation.max_turnover is infeasible for the current portfolio and target sleeve" - ) - if turnover(desired) <= max_turnover + 1e-12: - return desired.copy() - low = 0.0 - high = 1.0 - for _ in range(100): - midpoint = (low + high) / 2 - candidate = anchor + midpoint * (desired - anchor) - if turnover(candidate) <= max_turnover: - low = midpoint - else: - high = midpoint - return anchor + low * (desired - anchor) + try: + return solve_portfolio( + np.zeros(len(desired)), + None, + current, + OptimizationConstraints(min_weight=lower, max_weight=upper, max_turnover=max_turnover), + budget=total, + kind="target", + target=desired, + turnover_offset=turnover_offset, + ).weights + except OptimizationError as exc: + if exc.status == "infeasible": + raise OptimizationError( + "infeasible", + "allocation.max_turnover is infeasible for the current portfolio and budget", + exc.diagnostics, + ) from exc + raise def _project_l1_ball(values: np.ndarray, center: np.ndarray, radius: float) -> np.ndarray: @@ -1333,7 +1179,9 @@ def project_halfspace( detail = ", ".join(f"{name}={value:.6g}" for name, value in violations.items()) if not detail: detail = "projection did not converge" - raise ValueError(f"joint constraints are infeasible: {detail}") + from quant_portfolio.solver import OptimizationError + + raise OptimizationError("not_converged", f"joint projection did not converge: {detail}") def optimize_mean_variance( @@ -1357,8 +1205,10 @@ def optimize_mean_variance( impact_coefficient: float = 0.0, portfolio_nav: float | None = None, weight_steps: pd.Series | None = None, + solver_backend: str = "cvxpy", + timeout_seconds: float = 30.0, ) -> OptimizationResult: - """Proximal-gradient long-only optimizer. + """Joint convex long-only optimizer with an explicit legacy reference backend. Maximizes expected return minus quadratic risk, linear trading costs and an additional turnover penalty. Pass ``benchmark_weights`` to charge active @@ -1368,8 +1218,7 @@ def optimize_mean_variance( are all supplied. ``constraints.max_tracking_error`` is a hard cap on active volatility and requires ``benchmark_weights``. ``weight_steps`` snaps the continuous solution to integer lots and may leave residual cash. - The non-smooth trading objective and all constraints share one Dykstra - proximal step on every iteration. + All constraints and trading costs enter the same convex problem. """ if not isinstance(expected_returns, pd.Series): raise TypeError("expected_returns must be a pandas Series") @@ -1494,45 +1343,46 @@ def optimize_mean_variance( if not np.isfinite(costs).all() or (costs < 0).any(): raise ValueError("linear_costs must be finite and non-negative for every asset") mu = numeric_expected.to_numpy(dtype=float) - largest_eigenvalue = max(float(np.linalg.eigvalsh(cov).max()), 1e-12) - impact_lipschitz = 0.0 - if impact is not None: - adv, vol, coefficient, nav = impact - curvature = float(np.max(coefficient * vol * np.sqrt(nav / adv)) * 1.5 * (1e-8**-0.25)) - impact_lipschitz = curvature - step = 0.5 / (risk_aversion * largest_eigenvalue + impact_lipschitz + 1.0) - projection_tolerance = max(min(float(tolerance), 1e-10), 1e-12) + from quant_portfolio.solver import RiskMatrix, solve_portfolio + + diagnostics = {} tracking_covariance = cov if settings.max_tracking_error is not None else None tracking_benchmark = benchmark if settings.max_tracking_error is not None else None - weights = _project_joint_feasible_set( - current, - assets, - current, - float(turnover_offset), - settings, - validated_exposures, - validated_factor_bounds, - tolerance=projection_tolerance, - covariance=tracking_covariance, - benchmark=tracking_benchmark, - max_tracking_error=settings.max_tracking_error, - ) - converged = False - - for iteration in range(1, max_iterations + 1): - impact_gradient = np.zeros(len(assets)) + if solver_backend == "cvxpy": + solved = solve_portfolio( + mu, + RiskMatrix(covariance=cov), + current, + settings, + assets=assets, + factor_exposures=validated_exposures, + factor_bounds=validated_factor_bounds, + benchmark=benchmark if active_risk else None, + active_risk=active_risk, + linear_costs=costs, + turnover_penalty=turnover_penalty, + turnover_offset=turnover_offset, + risk_aversion=risk_aversion, + impact=impact, + timeout_seconds=timeout_seconds, + max_iterations=max_iterations, + ) + weights = solved.weights + diagnostics = solved.diagnostics + converged, iteration = True, diagnostics["iterations"] + elif solver_backend == "legacy": + largest_eigenvalue = max(float(np.linalg.eigvalsh(cov).max()), 1e-12) + impact_lipschitz = 0.0 if impact is not None: adv, vol, coefficient, nav = impact - _, impact_gradient = square_root_impact_penalty( - weights - current, - adv, - vol, - impact_coefficient=coefficient, - portfolio_nav=nav, - ) - gradient = mu - risk_aversion * (cov @ (weights - risk_center)) - impact_gradient - candidate = _project_joint_feasible_set( - weights + step * gradient, + curvature = float(np.max(coefficient * vol * np.sqrt(nav / adv)) * 1.5 * (1e-8**-0.25)) + impact_lipschitz = curvature + step = 0.5 / (risk_aversion * largest_eigenvalue + impact_lipschitz + 1.0) + projection_tolerance = max(min(float(tolerance), 1e-10), 1e-12) + tracking_covariance = cov if settings.max_tracking_error is not None else None + tracking_benchmark = benchmark if settings.max_tracking_error is not None else None + weights = _project_joint_feasible_set( + current, assets, current, float(turnover_offset), @@ -1540,19 +1390,49 @@ def optimize_mean_variance( validated_exposures, validated_factor_bounds, tolerance=projection_tolerance, - trade_penalties=costs + turnover_penalty, - penalty_scale=step, covariance=tracking_covariance, benchmark=tracking_benchmark, max_tracking_error=settings.max_tracking_error, ) - if float(np.max(np.abs(candidate - weights))) <= tolerance: + converged = False + + for iteration in range(1, max_iterations + 1): + impact_gradient = np.zeros(len(assets)) + if impact is not None: + adv, vol, coefficient, nav = impact + _, impact_gradient = square_root_impact_penalty( + weights - current, + adv, + vol, + impact_coefficient=coefficient, + portfolio_nav=nav, + ) + gradient = mu - risk_aversion * (cov @ (weights - risk_center)) - impact_gradient + candidate = _project_joint_feasible_set( + weights + step * gradient, + assets, + current, + float(turnover_offset), + settings, + validated_exposures, + validated_factor_bounds, + tolerance=projection_tolerance, + trade_penalties=costs + turnover_penalty, + penalty_scale=step, + covariance=tracking_covariance, + benchmark=tracking_benchmark, + max_tracking_error=settings.max_tracking_error, + ) + if float(np.max(np.abs(candidate - weights))) <= tolerance: + weights = candidate + converged = True + break weights = candidate - converged = True - break - weights = candidate + else: + iteration = max_iterations + else: - iteration = max_iterations + raise ValueError("solver_backend must be cvxpy or legacy") violations = _constraint_violations( weights, @@ -1562,7 +1442,7 @@ def optimize_mean_variance( settings, validated_exposures, validated_factor_bounds, - tolerance=max(float(tolerance) * 10, 1e-9), + tolerance=max(float(tolerance) * 10, 1e-7), covariance=tracking_covariance, benchmark=tracking_benchmark, max_tracking_error=settings.max_tracking_error, @@ -1625,6 +1505,7 @@ def optimize_mean_variance( group_weights=group_weights, converged=converged, iterations=iteration, + diagnostics=diagnostics, ) @@ -1801,117 +1682,32 @@ def optimize_cvar( if numeric_current.isna().any() or (numeric_current < 0).any(): raise ValueError("current_weights must be finite and non-negative") current = numeric_current.to_numpy(dtype=float) - scenario_values = scenarios.to_numpy(dtype=float) - drawdown_anchor = None - if max_drawdown is not None: - if ( - isinstance(max_drawdown, bool) - or not isinstance(max_drawdown, (int, float)) - or not np.isfinite(max_drawdown) - or max_drawdown < 0 - ): - raise ValueError("max_drawdown must be finite and non-negative") + from quant_portfolio.solver import solve_portfolio - def project_drawdown(point: np.ndarray) -> np.ndarray: - return _project_joint_feasible_set( - point, - assets, - current, - 0.0, - settings, - None, - {}, - tolerance=1e-10, - ) - - drawdown_anchor = minimum_drawdown_weights(scenario_values, project_drawdown) - if maximum_drawdown(scenario_values @ drawdown_anchor) > float(max_drawdown) + 1e-8: - raise ValueError("max_drawdown is infeasible for these scenarios and constraints") - mu_values = mu.to_numpy(dtype=float) - scale = max(float(np.max(np.abs(scenario_values))), 1e-6) - base_step = 0.25 / (float(risk_aversion) * scale + 1.0) - weights = _project_joint_feasible_set( - current, - assets, + solved = solve_portfolio( + mu.to_numpy(dtype=float), + None, current, - 0.0, settings, - None, - {}, - tolerance=1e-10, + assets=assets, + kind="cvar", + scenarios=scenarios.to_numpy(dtype=float), + beta=beta, + risk_aversion=risk_aversion, + max_drawdown=max_drawdown, + max_iterations=max_iterations, ) - converged = False - iteration = max_iterations - best_weights = weights.copy() - best_objective = -np.inf - quiet = 0 - - def objective_at(point: np.ndarray) -> float: - point_losses = -(scenario_values @ point) - point_var = float(np.quantile(point_losses, float(beta))) - tail_loss = point_losses[point_losses >= point_var - 1e-15] - cvar_value = float(tail_loss.mean()) if len(tail_loss) else float(point_losses.max()) - return float(mu_values @ point) - float(risk_aversion) * cvar_value - - for iteration in range(1, max_iterations + 1): - losses = -(scenario_values @ weights) - var = float(np.quantile(losses, float(beta))) - tail = losses >= var - 1e-15 - if not tail.any(): - tail = np.zeros(len(losses), dtype=bool) - tail[int(np.argmax(losses))] = True - cvar_gradient = -scenario_values[tail].mean(axis=0) - gradient = mu_values - float(risk_aversion) * cvar_gradient - step = base_step / (1.0 + 0.05 * iteration) - candidate = _project_joint_feasible_set( - weights + step * gradient, - assets, - current, - 0.0, - settings, - None, - {}, - tolerance=1e-10, - ) - if drawdown_anchor is not None: - candidate = enforce_drawdown_limit( - candidate, scenario_values, float(max_drawdown), drawdown_anchor - ) - candidate_objective = objective_at(candidate) - if candidate_objective > best_objective + tolerance: - best_objective = candidate_objective - best_weights = candidate.copy() - quiet = 0 - else: - quiet += 1 - if float(np.max(np.abs(candidate - weights))) <= tolerance or quiet >= 30 or step < 1e-4: - weights = best_weights - converged = True - break - weights = candidate - if not converged: - raise RuntimeError(f"CVaR optimization did not converge after {max_iterations} iterations") - if ( - drawdown_anchor is not None - and maximum_drawdown(scenario_values @ weights) > float(max_drawdown) + 1e-8 - ): - raise RuntimeError("optimizer produced a portfolio outside the drawdown limit") - portfolio_returns = scenario_values @ weights - volatility = float(np.std(portfolio_returns, ddof=1)) - expected = float(mu_values @ weights) - losses = -(scenario_values @ weights) - var = float(np.quantile(losses, float(beta))) - tail_losses = losses[losses >= var - 1e-15] - cvar = float(tail_losses.mean()) + w = solved.weights return OptimizationResult( - weights=pd.Series(weights, index=assets, name="weight"), - expected_return=expected, - volatility=volatility, - turnover=float(np.abs(weights - current).sum()), - objective=expected - float(risk_aversion) * cvar, + weights=pd.Series(w, index=assets, name="weight"), + expected_return=float(mu @ w), + volatility=float(np.std(scenarios.to_numpy() @ w, ddof=1)), + turnover=float(np.abs(w - current).sum()), + objective=solved.objective, group_weights={}, converged=True, - iterations=iteration, + iterations=solved.diagnostics["iterations"], + diagnostics=solved.diagnostics, ) @@ -2010,77 +1806,62 @@ def optimize_multiperiod( if not np.isfinite(numeric_costs.to_numpy(dtype=float)).all() or (numeric_costs < 0).any(): raise ValueError("linear_costs must be finite and non-negative for every asset") costs = numeric_costs.to_numpy(dtype=float) - eigenvalues = np.linalg.eigvalsh(cov) - step = 0.5 / (float(risk_aversion) * max(float(eigenvalues.max()), 1e-12) + 1.0) - projection_tolerance = max(min(float(tolerance), 1e-10), 1e-12) - anchor = _project_joint_feasible_set( - current, - assets, - current, - 0.0, - settings, - None, - {}, - tolerance=projection_tolerance, - ) - plan = np.tile(anchor, (len(mu_path), 1)) - mu = mu_path.to_numpy(dtype=float) - converged = False - iteration = max_iterations - for iteration in range(1, max_iterations + 1): - updated = np.empty_like(plan) - previous = current - for period in range(len(mu)): - gradient = mu[period] - float(risk_aversion) * (cov @ plan[period]) - proposal = plan[period] + step * gradient - proposal = _project_joint_feasible_set( - proposal, - assets, - previous, - 0.0, - settings, - None, - {}, - tolerance=projection_tolerance, - trade_penalties=costs, - penalty_scale=step, - ) - updated[period] = proposal - previous = proposal - if float(np.max(np.abs(updated - plan))) <= tolerance: - plan = updated - converged = True - break - plan = updated - if not converged: - raise RuntimeError( - f"multi-period optimization did not converge after {max_iterations} iterations" + from time import perf_counter + + import cvxpy as cp + + from quant_portfolio.solver import OptimizationError, RiskMatrix, _run + + risk = RiskMatrix(covariance=cov) + plan = cp.Variable(mu_path.shape) + rules, terms = [], [] + previous = current + for period, mu in enumerate(mu_path.to_numpy(dtype=float)): + w = plan[period] + trade = w - previous + rules.extend( + [ + cp.sum(w) == 1, + w >= settings.min_weight, + w <= settings.max_weight, + cp.norm1(trade) <= settings.max_turnover, + ] + ) + for group, cap in settings.group_caps.items(): + mask = [i for i, a in enumerate(assets) if settings.group_by_asset.get(a) == group] + rules.append(cp.sum(w[mask]) <= cap) + terms.append( + mu @ w + - 0.5 * risk_aversion * cp.sum_squares(risk.root_expression(cp, w)) + - costs @ cp.abs(trade) + ) + previous = w + problem = cp.Problem(cp.Maximize(sum(terms)), rules) + diagnostics = {} + if not _run( + cp, + problem, + conic=False, + deadline=perf_counter() + 30, + diagnostics=diagnostics, + max_iterations=max_iterations, + ): + raise OptimizationError( + "infeasible", "multi-period joint constraints are infeasible", diagnostics ) - objective = 0.0 previous = current - for period in range(len(mu)): + for row in plan.value: violations = _constraint_violations( - plan[period], - assets, - previous, - 0.0, - settings, - None, - {}, - tolerance=max(projection_tolerance * 10, 1e-9), + row, assets, previous, 0.0, settings, None, {}, tolerance=1e-7 ) - if not np.isfinite(plan[period]).all() or violations: - raise RuntimeError( - f"multi-period plan violates constraints at period {period}: {violations}" + if not np.isfinite(row).all() or violations: + raise OptimizationError( + "numerical_failure", f"multi-period plan violates constraints: {violations}" ) - trade = plan[period] - previous - objective += float(mu[period] @ plan[period]) - objective -= 0.5 * float(risk_aversion) * float(plan[period] @ cov @ plan[period]) - objective -= float(costs @ np.abs(trade)) - previous = plan[period] + previous = row return MultiPeriodResult( - weights=pd.DataFrame(plan, index=expected_returns.index, columns=assets), - objective=objective, + weights=pd.DataFrame(plan.value, index=expected_returns.index, columns=assets), + objective=float(problem.value), converged=True, - iterations=iteration, + iterations=diagnostics["iterations"], ) diff --git a/src/quant_portfolio/risk_snapshot.py b/src/quant_portfolio/risk_snapshot.py new file mode 100644 index 0000000..2b34f90 --- /dev/null +++ b/src/quant_portfolio/risk_snapshot.py @@ -0,0 +1,191 @@ +"""Point-in-time, labeled, unit-bearing optimizer risk inputs.""" + +from __future__ import annotations + +import hashlib +from dataclasses import dataclass, field + +import numpy as np +import pandas as pd + +from quant_portfolio.solver import RiskMatrix + + +def timestamp(value): + result = pd.Timestamp(value) + if pd.isna(result) or result.tzinfo is None: + raise ValueError("risk timestamps must be valid and timezone-aware") + return result.tz_convert("UTC") + + +@dataclass(frozen=True, kw_only=True) +class RiskSnapshot: + estimated_at: object + available_at: object + horizon_sessions: int + annualization: int + model_version: str + input_digest: str + covariance: pd.DataFrame | None = None + exposures: pd.DataFrame | None = None + factor_covariance: pd.DataFrame | None = None + specific_variance: pd.Series | None = None + benchmark: pd.Series | None = None + _fingerprint: str | None = field(default=None, init=False, repr=False) + + def __post_init__(self): + for key in ("estimated_at", "available_at"): + object.__setattr__(self, key, timestamp(getattr(self, key))) + if self.estimated_at > self.available_at: + raise ValueError("risk cannot be available before estimation") + for key in ("horizon_sessions", "annualization"): + if type(getattr(self, key)) is not int or getattr(self, key) < 1: + raise ValueError(f"invalid risk {key}") + if not self.model_version or not self.input_digest: + raise ValueError("risk model version and input digest are required") + for key in ( + "covariance", + "exposures", + "factor_covariance", + "specific_variance", + "benchmark", + ): + value = getattr(self, key) + if value is not None: + if not value.index.is_unique or ( + isinstance(value, pd.DataFrame) and not value.columns.is_unique + ): + raise ValueError("risk labels must be unique") + object.__setattr__(self, key, value.copy(deep=True)) + self.matrix(self.assets, self.available_at) + object.__setattr__(self, "_fingerprint", self._digest()) + + @property + def assets(self): + source = self.covariance if self.covariance is not None else self.exposures + if source is None: + raise ValueError("risk representation is missing") + return list(source.index) + + def matrix(self, assets, decision_at): + self._check_unchanged() + if self.available_at > timestamp(decision_at): + raise ValueError("risk snapshot was unavailable at decision time") + if set(assets) - set(self.assets): + raise ValueError( + "risk snapshot does not cover candidates, holdings and benchmark union" + ) + if self.benchmark is not None: + b = self.benchmark + if ( + not np.isfinite(b).all() + or (b < 0).any() + or b.sum() > 1 + 1e-10 + or set(b.index) - set(self.assets) + ): + raise ValueError("benchmark must be complete and covered by risk snapshot") + if self.covariance is not None: + if any( + x is not None + for x in (self.exposures, self.factor_covariance, self.specific_variance) + ): + raise ValueError("provide exactly one risk representation") + if list(self.covariance.index) != list(self.covariance.columns): + raise ValueError("covariance labels must match in order") + return RiskMatrix(covariance=self.covariance.loc[assets, assets].to_numpy()) + if self.factor_covariance is None or self.specific_variance is None: + raise ValueError( + "factor risk requires exposures, factor covariance and specific variance" + ) + factors = list(self.exposures.columns) + if ( + set(self.factor_covariance.index) != set(factors) + or set(self.factor_covariance.columns) != set(factors) + or set(self.specific_variance.index) != set(self.assets) + ): + raise ValueError("factor risk labels must match exactly") + return RiskMatrix( + exposures=self.exposures.loc[assets].to_numpy(), + factor_covariance=self.factor_covariance.loc[factors, factors].to_numpy(), + specific_variance=self.specific_variance.loc[assets].to_numpy(), + ) + + def _check_unchanged(self): + if self._fingerprint is not None and self._digest() != self._fingerprint: + raise ValueError("risk snapshot was mutated after validation") + + def _digest(self): + digest = hashlib.sha256() + for value in ( + self.covariance, + self.exposures, + self.factor_covariance, + self.specific_variance, + self.benchmark, + ): + if value is not None: + digest.update(pd.util.hash_pandas_object(value, index=True).to_numpy().tobytes()) + if isinstance(value, pd.DataFrame): + digest.update(repr(tuple(value.columns)).encode()) + return digest.hexdigest() + + def metadata(self): + self._check_unchanged() + return { + "estimated_at": self.estimated_at.isoformat(), + "available_at": self.available_at.isoformat(), + "horizon_sessions": self.horizon_sessions, + "annualization": self.annualization, + "units": "annualized_variance", + "model_version": self.model_version, + "input_digest": self.input_digest, + "snapshot_digest": self._fingerprint, + } + + +def mature_ic(labels, *, decision_at, horizon_sessions, lookback, min_periods): + """Cross-sectional IC from available, non-overlapping, completed labels only.""" + required = { + "signal_at", + "label_end", + "available_at", + "asset", + "score", + "return", + "horizon_sessions", + } + if not isinstance(labels, pd.DataFrame) or not required <= set(labels.columns): + raise ValueError("ic_vol requires timestamped mature label history") + data = labels.copy() + for key in ("signal_at", "label_end", "available_at"): + data[key] = data[key].map(timestamp) + cutoff = timestamp(decision_at) + data = data.loc[ + (data.available_at <= cutoff) + & (data.label_end <= cutoff) + & (data.signal_at < cutoff) + & (data.horizon_sessions == horizon_sessions) + ] + if ((data.label_end <= data.signal_at) | (data.available_at < data.label_end)).any(): + raise ValueError("invalid label availability") + # Keep the last revision available at the decision, never a future revision. + data = data.sort_values("available_at").drop_duplicates(["signal_at", "asset"], keep="last") + periods, next_start = [], cutoff + for start, rows in reversed(list(data.groupby("signal_at", sort=True))): + if rows.label_end.nunique() != 1 or rows.label_end.max() > next_start: + continue + values = rows[["score", "return"]].to_numpy(dtype=float) + if len(rows) < 3 or not np.isfinite(values).all() or (np.std(values, axis=0) == 0).any(): + continue + periods.append((start, float(np.corrcoef(values.T)[0, 1]))) + next_start = start + if len(periods) == lookback: + break + if len(periods) < min_periods: + raise ValueError("insufficient mature non-overlapping IC periods") + return float(np.mean([x[1] for x in periods])), { + "ic_periods": len(periods), + "ic_first_signal": periods[-1][0].isoformat(), + "ic_last_signal": periods[0][0].isoformat(), + "ic_horizon_sessions": horizon_sessions, + } diff --git a/src/quant_portfolio/solver.py b/src/quant_portfolio/solver.py new file mode 100644 index 0000000..ec751c1 --- /dev/null +++ b/src/quant_portfolio/solver.py @@ -0,0 +1,589 @@ +"""Auditable convex portfolio problems; no strategy-changing solver fallback.""" + +from __future__ import annotations + +import warnings +from dataclasses import dataclass, field +from time import perf_counter +from typing import Any + +import numpy as np + + +class OptimizationError(ValueError, RuntimeError): + """A classified failure. A numerical failure is not proof of infeasibility.""" + + def __init__(self, status: str, message: str, diagnostics=None): + super().__init__(f"{status}: {message}") + self.status = status + self.diagnostics = diagnostics or {} + + +@dataclass(frozen=True) +class RiskMatrix: + """Covariance or X F X' + diag(D), without forming the asset matrix.""" + + covariance: np.ndarray | None = None + exposures: np.ndarray | None = None + factor_covariance: np.ndarray | None = None + specific_variance: np.ndarray | None = None + + def __post_init__(self): + dense = self.covariance is not None + factor = any( + x is not None for x in (self.exposures, self.factor_covariance, self.specific_variance) + ) + if dense == factor: + raise ValueError("provide exactly one dense or factor risk representation") + names = ( + ("covariance",) if dense else ("exposures", "factor_covariance", "specific_variance") + ) + for name in names: + value = getattr(self, name) + if value is None: + raise ValueError("factor risk requires X, F and D") + value = np.array(value, dtype=float, copy=True) + if not np.isfinite(value).all(): + raise ValueError("risk matrices must be finite") + value.setflags(write=False) + object.__setattr__(self, name, value) + matrix = self.covariance if dense else self.factor_covariance + if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1] or not len(matrix): + raise ValueError("covariance must be square") + if not np.allclose(matrix, matrix.T, rtol=1e-10, atol=1e-12): + raise ValueError("covariance must be symmetric") + if np.linalg.eigvalsh(matrix).min() < -1e-12: + raise ValueError("covariance must be positive semidefinite") + if not dense: + x, d = self.exposures, self.specific_variance + if x.ndim != 2 or x.shape[1] != len(matrix) or d.shape != (len(x),): + raise ValueError("factor risk dimensions differ") + if (d < 0).any(): + raise ValueError("specific variance must be non-negative") + + def variance(self, weights): + w = np.asarray(weights) + if self.covariance is not None: + return float(w @ self.covariance @ w) + exposure = self.exposures.T @ w + return float(exposure @ self.factor_covariance @ exposure + self.specific_variance @ w**2) + + def diagonal(self): + if self.covariance is not None: + return np.diag(self.covariance).copy() + return ( + np.einsum("ij,ij->i", self.exposures @ self.factor_covariance, self.exposures) + + self.specific_variance + ) + + def root_expression(self, cp, weights): + matrix = self.covariance if self.covariance is not None else self.factor_covariance + eigenvalues, vectors = np.linalg.eigh(matrix) + root = np.sqrt(np.maximum(eigenvalues, 0))[:, None] * vectors.T + if self.covariance is not None: + return root @ weights + return cp.hstack( + ( + root @ (self.exposures.T @ weights), + cp.multiply(np.sqrt(self.specific_variance), weights), + ) + ) + + +@dataclass +class SolveResult: + weights: np.ndarray + objective: float + diagnostics: dict[str, Any] = field(default_factory=dict) + + +def constraint_residuals( + weights, + current, + settings, + *, + assets, + budget=1.0, + lower=None, + upper=None, + factor_exposures=None, + factor_bounds=None, + risk=None, + benchmark=None, + turnover_offset=0.0, + budget_equality=True, +): + """Independent NumPy acceptance check; positive values denote violations.""" + w = np.asarray(weights, float) + if not np.isfinite(w).all(): + return {"finite": float("inf")} + lo = settings.min_weight if lower is None else np.asarray(lower) + hi = settings.max_weight if upper is None else np.asarray(upper) + result = { + "budget": abs(float(w.sum()) - budget) if budget_equality else float(w.sum()) - budget, + "asset_lower": float(np.max(lo - w)), + "asset_upper": float(np.max(w - hi)), + "turnover": float(np.abs(w - current).sum()) + turnover_offset - settings.max_turnover, + } + for group, cap in settings.group_caps.items(): + result["group:" + group] = ( + sum(w[i] for i, a in enumerate(assets) if settings.group_by_asset.get(a) == group) - cap + ) + if factor_bounds: + exposure = np.asarray(factor_exposures).T @ w + for i, (name, (left, right)) in enumerate(factor_bounds.items()): + result["factor_lower:" + name] = float(left - exposure[i]) + result["factor_upper:" + name] = float(exposure[i] - right) + if settings.max_tracking_error is not None: + result["tracking_error"] = ( + np.sqrt(max(risk.variance(w - benchmark), 0)) - settings.max_tracking_error + ) + return result + + +def _run(cp, problem, *, conic, deadline, diagnostics, max_iterations): + attempts = diagnostics.setdefault("attempts", []) + solvers = ("CLARABEL", "CLARABEL") if conic else ("OSQP", "CLARABEL") + for attempt, solver in enumerate(solvers): + remaining = deadline - perf_counter() + if remaining <= 0: + raise OptimizationError("timeout", "rebalance solve deadline exceeded", diagnostics) + started = perf_counter() + try: + if solver == "OSQP": + options = { + "eps_abs": 1e-10, + "eps_rel": 1e-10, + "max_iter": max_iterations, + "time_limit": remaining, + "polishing": True, + } + else: + options = { + "tol_gap_abs": 1e-9 if attempt == 0 else 1e-8, + "tol_gap_rel": 1e-9 if attempt == 0 else 1e-8, + "tol_feas": 1e-9, + "max_iter": max_iterations, + "time_limit": remaining, + } + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", message="Solution may be inaccurate.*") + problem.solve(solver=solver, warm_start=attempt == 0, **options) + except cp.error.SolverError as exc: + attempts.append({"solver": solver, "status": "numerical_failure", "error": str(exc)}) + continue + elapsed = perf_counter() - started + status = problem.status + residual = ( + max((float(np.max(np.abs(c.violation()))) for c in problem.constraints), default=0.0) + if status == cp.OPTIMAL + else None + ) + stats = problem.solver_stats + attempts.append( + { + "solver": solver, + "status": status, + "seconds": elapsed, + "compilation_seconds": problem.compilation_time, + "solver_seconds": stats.solve_time, + "iterations": stats.num_iters, + "max_constraint_residual": residual, + "options": options, + } + ) + if perf_counter() > deadline: + raise OptimizationError("timeout", "rebalance solve deadline exceeded", diagnostics) + if status == cp.OPTIMAL and residual is not None and residual <= 1e-7: + diagnostics.update( + status="optimal", solver=solver, residual=residual, iterations=stats.num_iters or 0 + ) + return True + if status == cp.INFEASIBLE: + diagnostics.update(status="infeasible", solver=solver) + return False + if status == cp.UNBOUNDED: + raise OptimizationError("unbounded", "objective is unbounded", diagnostics) + status = ( + "timeout" + if perf_counter() >= deadline + else ( + "not_converged" + if any(a["status"] == cp.USER_LIMIT for a in attempts) + else "numerical_failure" + ) + ) + raise OptimizationError( + status, "no accurate solver result; constraints were not relaxed", diagnostics + ) + + +def solve_portfolio( + expected, + risk, + current, + settings, + *, + assets=None, + budget=1.0, + lower=None, + upper=None, + factor_exposures=None, + factor_bounds=None, + benchmark=None, + active_risk=False, + linear_costs=None, + buy_costs=None, + sell_costs=None, + turnover_penalty=0.0, + turnover_offset=0.0, + risk_aversion=5.0, + risk_scale=1.0, + impact=None, + kind="mean_variance", + target=None, + scenarios=None, + beta=0.95, + max_drawdown=None, + risk_free_rate=0.0, + timeout_seconds=30.0, + max_iterations=2000, +): + """Solve with all constraints present, including real pre-trade holdings.""" + import cvxpy as cp + + started = perf_counter() + if ( + isinstance(timeout_seconds, bool) + or not np.isfinite(timeout_seconds) + or timeout_seconds <= 0 + ): + raise ValueError("timeout_seconds must be positive and finite") + if type(max_iterations) is not int or max_iterations < 1: + raise ValueError("max_iterations must be a positive integer") + deadline = started + timeout_seconds + mu, current = np.asarray(expected, float), np.asarray(current, float) + n = len(mu) + if n == 0 or mu.shape != (n,) or current.shape != (n,): + raise ValueError("expected returns and current weights must be aligned vectors") + if not np.isfinite(mu).all() or not np.isfinite(current).all() or (current < 0).any(): + raise ValueError( + "expected returns and current weights must be finite; holdings non-negative" + ) + if not np.isfinite( + [budget, risk_scale, risk_aversion, turnover_offset, turnover_penalty] + ).all(): + raise ValueError("problem parameters must be finite") + if ( + not 0 <= budget <= 1 + or risk_scale <= 0 + or risk_aversion <= 0 + or turnover_offset < 0 + or turnover_penalty < 0 + ): + raise ValueError("invalid budget, risk scale, aversion or turnover offset") + assets = list(range(n)) if assets is None else list(assets) + if len(assets) != n or len(set(assets)) != n: + raise ValueError("assets must be unique and aligned") + if not np.isfinite(settings.max_turnover) or settings.max_turnover < 0: + raise ValueError("max_turnover must be finite and non-negative") + if settings.max_tracking_error is not None and ( + not np.isfinite(settings.max_tracking_error) or settings.max_tracking_error < 0 + ): + raise ValueError("max_tracking_error must be finite and non-negative") + if factor_bounds and factor_exposures is None: + raise ValueError("factor bounds require factor exposures") + lo = np.broadcast_to(settings.min_weight if lower is None else lower, (n,)).copy() + hi = np.broadcast_to(settings.max_weight if upper is None else upper, (n,)).copy() + if ( + not np.isfinite(lo).all() + or not np.isfinite(hi).all() + or (lo < 0).any() + or (hi > 1).any() + or (lo > hi).any() + ): + raise ValueError("invalid per-asset bounds") + if lo.sum() > budget + 1e-10 or hi.sum() < budget - 1e-10: + raise OptimizationError("infeasible", "weight bounds cannot satisfy the budget constraint") + wb = np.zeros(n) if benchmark is None else np.asarray(benchmark, float) + if wb.shape != (n,) or not np.isfinite(wb).all() or (wb < 0).any() or wb.sum() > 1 + 1e-8: + raise ValueError("benchmark must be a complete non-negative weight vector") + if (active_risk or settings.max_tracking_error is not None) and benchmark is None: + raise ValueError("active risk and max_tracking_error require benchmark_weights") + w = cp.Variable(n, name="weights") + trade = w - current + constraints = [ + cp.sum(w) == budget, + w >= lo, + w <= hi, + cp.norm1(trade) + turnover_offset <= settings.max_turnover, + ] + labels = ["budget", "asset_lower", "asset_upper", "turnover"] + for group, cap in settings.group_caps.items(): + if not np.isfinite(cap) or not 0 <= cap <= 1: + raise ValueError("group caps must be in [0, 1]") + mask = np.array([settings.group_by_asset.get(a, "__ungrouped__") == group for a in assets]) + constraints.append(cp.sum(w[mask]) <= cap) + labels.append("group:" + group) + x = None if factor_exposures is None else np.asarray(factor_exposures, float) + if x is not None: + bounds = list((factor_bounds or {}).values()) + if x.shape != (n, len(bounds)) or not np.isfinite(x).all(): + raise ValueError("factor exposures and bounds must align") + for index, (name, (left, right)) in enumerate(factor_bounds.items()): + if left > right or not np.isfinite([left, right]).all(): + raise ValueError("factor bounds must be finite ordered pairs") + constraints.extend([x[:, index] @ w >= left, x[:, index] @ w <= right]) + labels.extend(["factor_lower:" + name, "factor_upper:" + name]) + if risk is None: + if settings.max_tracking_error is not None or kind in { + "mean_variance", + "max_sharpe", + "max_diversification", + }: + raise ValueError("this objective requires risk inputs") + risk_root = None + else: + risk_root = risk.root_expression(cp, w) + conic = settings.max_tracking_error is not None or impact is not None + if settings.max_tracking_error is not None: + constraints.append( + cp.norm(risk.root_expression(cp, w - wb), 2) <= settings.max_tracking_error + ) + labels.append("tracking_error") + cost = np.zeros(n) if linear_costs is None else np.asarray(linear_costs, float) + buy = cost if buy_costs is None else np.asarray(buy_costs, float) + sell = cost if sell_costs is None else np.asarray(sell_costs, float) + if any(c.shape != (n,) or not np.isfinite(c).all() or (c < 0).any() for c in (buy, sell)): + raise ValueError("cost rates must be complete, finite and non-negative") + penalty = ( + buy @ cp.pos(trade) + + sell @ cp.pos(-trade) + + turnover_penalty * (cp.norm1(trade) + turnover_offset) + ) + if impact is not None: + adv, vol, coefficient, nav = impact + adv, vol = np.asarray(adv, float), np.asarray(vol, float) + if ( + adv.shape != (n,) + or vol.shape != (n,) + or not np.isfinite(adv).all() + or not np.isfinite(vol).all() + or (adv <= 0).any() + or (vol < 0).any() + or not np.isfinite([coefficient, nav]).all() + or coefficient < 0 + or nav <= 0 + ): + raise ValueError( + "impact requires positive ADV/NAV and non-negative volatility/coefficient" + ) + coefficients = coefficient * np.asarray(vol) * np.sqrt(nav / np.asarray(adv)) + penalty += coefficients @ cp.power(cp.abs(trade), 1.5) + diagnostics = { + "backend": "cvxpy", + "objective_kind": kind, + "constraints": labels, + "risk_representation": "none" + if risk is None + else ("dense" if risk.covariance is not None else "factor"), + "attempts": [], + } + + if kind in {"max_sharpe", "max_diversification"}: + numerator = ( + mu @ w - risk_free_rate + if kind == "max_sharpe" + else np.sqrt(risk.diagonal() * risk_scale) @ w + ) + scale = ( + max(float(np.max(np.abs(mu))), abs(risk_free_rate), 1e-4) + if kind == "max_sharpe" + else max(float(np.max(np.sqrt(risk.diagonal() * risk_scale))), 1e-4) + ) + positive = cp.Problem(cp.Maximize(numerator / scale), constraints) + if not _run( + cp, + positive, + conic=True, + deadline=deadline, + diagnostics=diagnostics, + max_iterations=max_iterations, + ): + raise OptimizationError("infeasible", "joint constraints are infeasible", diagnostics) + if float(numerator.value) <= 1e-12: + raise OptimizationError( + "no_positive_excess_return", "ratio has no positive numerator", diagnostics + ) + if risk.variance(np.asarray(w.value)) <= 1e-18: + raise OptimizationError( + "undefined_ratio", "positive numerator with zero risk", diagnostics + ) + ratio = cp.Parameter(nonneg=True) + # Maximize the feasibility margin. This avoids an ill-conditioned + # zero-objective feasibility problem exactly at the ratio boundary. + problem = cp.Problem( + cp.Maximize((numerator - ratio * np.sqrt(risk_scale) * cp.norm(risk_root, 2)) / scale), + constraints, + ) + + def feasible_ratio(value): + ratio.value = value + if not _run( + cp, + problem, + conic=True, + deadline=deadline, + diagnostics=diagnostics, + max_iterations=max_iterations, + ): + raise OptimizationError("numerical_failure", "ratio domain changed", diagnostics) + return float(problem.value) >= 0.0 + + low, high = 0.0, 1.0 + best = np.asarray(w.value).copy() + for _ in range(40): + feasible = feasible_ratio(high) + if not feasible: + break + best = np.asarray(w.value).copy() + low, high = high, 2 * high + else: + raise OptimizationError( + "undefined_ratio", "zero risk or unbounded positive ratio", diagnostics + ) + for _ in range(60): + if high - low <= 1e-9 * max(1.0, low): + break + if feasible_ratio((low + high) / 2): + low, best = ratio.value, np.asarray(w.value).copy() + else: + high = ratio.value + w.value = best + variance = risk.variance(best) + if variance <= 1e-18: + raise OptimizationError("undefined_ratio", "zero-volatility portfolio", diagnostics) + objective_value = float(numerator.value) / np.sqrt(variance * risk_scale) + diagnostics["ratio_bracket"] = [low, high] + else: + if kind == "mean_variance": + centered = w - wb if active_risk else w + objective = ( + mu @ w + - risk_aversion + * 0.5 + * risk_scale + * cp.sum_squares(risk.root_expression(cp, centered)) + - penalty + ) + elif kind == "target": + target = np.asarray(target, float) + if target.shape != (n,) or not np.isfinite(target).all(): + raise ValueError("target must be finite and aligned") + objective = -cp.sum_squares(w - target) - penalty + elif kind in {"cvar", "evar", "cdar"}: + rows = np.asarray(scenarios, float) + if rows.ndim != 2 or rows.shape[1] != n or len(rows) < 2 or not np.isfinite(rows).all(): + raise ValueError("scenarios must be finite and complete with at least two rows") + if not 0 < beta < 1: + raise ValueError("beta must be in (0, 1)") + losses = -(rows @ w) + if kind == "cdar" or max_drawdown is not None: + cumulative = np.cumsum(rows, axis=0) @ w + peaks = cp.Variable(len(rows)) + constraints.extend([peaks >= 0, peaks >= cumulative, peaks[1:] >= peaks[:-1]]) + drawdowns = peaks - cumulative + if max_drawdown is not None: + if not np.isfinite(max_drawdown) or max_drawdown < 0: + raise ValueError("max_drawdown must be finite and non-negative") + constraints.append(drawdowns <= max_drawdown) + if kind == "cdar": + losses = drawdowns + threshold = cp.Variable() + if kind == "evar": + temperature = cp.Variable(nonneg=True) + aux = cp.Variable(len(rows)) + constraints.extend( + [ + cp.ExpCone( + losses - threshold, cp.multiply(np.ones(len(rows)), temperature), aux + ), + cp.sum(aux) <= temperature * len(rows) * (1 - beta), + ] + ) + tail_risk = threshold + else: + excess = cp.Variable(len(rows), nonneg=True) + constraints.append(excess >= losses - threshold) + tail_risk = threshold + cp.sum(excess) / (len(rows) * (1 - beta)) + objective = mu @ w - risk_aversion * tail_risk - penalty + conic = True + else: + raise ValueError(f"unsupported objective {kind}") + problem = cp.Problem(cp.Maximize(objective), constraints) + if not _run( + cp, + problem, + conic=conic, + deadline=deadline, + diagnostics=diagnostics, + max_iterations=max_iterations, + ): + raise OptimizationError("infeasible", "joint constraints are infeasible", diagnostics) + objective_value = float(objective.value) + weights = np.asarray(w.value, float) + checks = constraint_residuals( + weights, + current, + settings, + assets=assets, + budget=budget, + lower=lo, + upper=hi, + factor_exposures=x, + factor_bounds=factor_bounds, + risk=risk, + benchmark=wb, + turnover_offset=turnover_offset, + ) + residual = max(0.0, *checks.values()) + if kind in {"cvar", "evar", "cdar"} and max_drawdown is not None: + path = np.r_[0.0, np.cumsum(np.asarray(scenarios) @ weights)] + drawdown = float(np.max(np.maximum.accumulate(path) - path)) + residual = max(residual, drawdown - max_drawdown) + if not np.isfinite(weights).all() or residual > 1e-7: + raise OptimizationError( + "numerical_failure", "independent constraint recheck failed", diagnostics + ) + diagnostics.update( + status="optimal", + max_constraint_residual=residual, + constraint_residuals=checks, + total_seconds=perf_counter() - started, + iterations=sum(a.get("iterations") or 0 for a in diagnostics["attempts"]), + ) + delta = weights - current + estimated_buy = float(buy @ np.maximum(delta, 0)) + estimated_sell = float(sell @ np.maximum(-delta, 0)) + estimated_impact = 0.0 if impact is None else float(coefficients @ np.abs(delta) ** 1.5) + diagnostics["estimated_costs"] = { + "buy": estimated_buy, + "sell": estimated_sell, + "impact": estimated_impact, + "total": estimated_buy + estimated_sell + estimated_impact, + "units": "fraction_of_pre_trade_nav", + } + diagnostics["costs_in_objective"] = kind not in {"max_sharpe", "max_diversification"} + diagnostics["constraint_values"] = { + "budget": {"value": float(weights.sum()), "target": budget}, + "turnover": { + "value": float(np.abs(delta).sum()) + turnover_offset, + "upper": settings.max_turnover, + }, + "asset_weights": dict(zip(assets, weights.tolist())), + "asset_lower": dict(zip(assets, lo.tolist())), + "asset_upper": dict(zip(assets, hi.tolist())), + "factor_exposures": {} if x is None else dict(zip(factor_bounds, (x.T @ weights).tolist())), + "factor_bounds": factor_bounds or {}, + } + return SolveResult(weights.copy(), objective_value, diagnostics) diff --git a/tests/test_governance.py b/tests/test_governance.py index 2ae9615..8e954ff 100644 --- a/tests/test_governance.py +++ b/tests/test_governance.py @@ -26,21 +26,21 @@ def test_workspace_declaration_and_internal_release_tags_are_locked() -> None: dependencies = project["project"]["dependencies"] assert ( - "quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@f85128307dd0e3407d24350282fe7b679b591309" + "quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@ad3440700b226497c40093197f1df4321b49d1e8" in dependencies ) assert ( - "quant-execution @ git+https://github.com/PureSaber/quant-execution.git@49ef53d2c502ad6faad6818f853dadd1f8b8b182" + "quant-execution @ git+https://github.com/PureSaber/quant-execution.git@8423726a9d80b6509d2c016a7478fc0f40a29774" in dependencies ) lock = (ROOT / "requirements.lock").read_text(encoding="utf-8") assert ( - "quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@f85128307dd0e3407d24350282fe7b679b591309" + "quant-data-kit @ git+https://github.com/PureSaber/quant-data-kit.git@ad3440700b226497c40093197f1df4321b49d1e8" in lock ) assert ( - "quant-execution @ git+https://github.com/PureSaber/quant-execution.git@49ef53d2c502ad6faad6818f853dadd1f8b8b182" + "quant-execution @ git+https://github.com/PureSaber/quant-execution.git@8423726a9d80b6509d2c016a7478fc0f40a29774" in lock ) assert 'tomli==2.4.1 ; python_version < "3.11"' in lock diff --git a/tests/test_multiperiod_constraints.py b/tests/test_multiperiod_constraints.py index 7c83c6f..313376b 100644 --- a/tests/test_multiperiod_constraints.py +++ b/tests/test_multiperiod_constraints.py @@ -28,7 +28,12 @@ def test_multiperiod_enforces_group_caps_on_every_date(): assert weights[["A", "B"]].sum() <= 0.4 + 1e-8 assert (weights - previous).abs().sum() <= 0.2 + 1e-8 previous = weights - assert result.weights.iloc[0]["A"] > current["A"] + # Joint optimum anticipates tomorrow's B signal; the old greedy path + # bought A today and could not move enough into B tomorrow. + assert result.weights.to_numpy() == pytest.approx( + np.array([[0.1, 0.3, 0.6], [0, 0.4, 0.6]]), abs=1e-7 + ) + assert result.objective == pytest.approx(0.2051, abs=1e-7) def test_multiperiod_rejects_group_cap_and_turnover_conflict(): @@ -135,10 +140,12 @@ def test_multiperiod_rejects_invalid_solver_inputs(kwargs, message): def test_multiperiod_independently_checks_returned_plan(monkeypatch): # A regressed solver must not publish an infeasible plan as converged. - monkeypatch.setattr( - "quant_portfolio.optimization._project_joint_feasible_set", - lambda *args, **kwargs: np.array([0.8, 0.2]), - ) + def faulty_run(cp, problem, **kwargs): + problem.variables()[0].value = np.array([[0.8, 0.2]]) + kwargs["diagnostics"]["iterations"] = 1 + return True + + monkeypatch.setattr("quant_portfolio.solver._run", faulty_run) with pytest.raises(RuntimeError, match="violates constraints"): optimize_multiperiod( pd.DataFrame([[0.1, 0.0]], columns=["A", "B"]), diff --git a/tests/test_optimization.py b/tests/test_optimization.py index 0f2360c..aeb6289 100644 --- a/tests/test_optimization.py +++ b/tests/test_optimization.py @@ -161,7 +161,7 @@ def test_equal_allocation_moves_toward_target_within_turnover_limit() -> None: config={"mode": "equal", "max_turnover": 0.2}, ) assert weights.to_dict() == pytest.approx({"A": 0.4, "B": 0.2}) - assert (weights - pd.Series({"A": 0.5, "B": 0.1})).abs().sum() <= 0.2 + 1e-12 + assert (weights - pd.Series({"A": 0.5, "B": 0.1})).abs().sum() <= 0.2 + 1e-9 def test_research_allocation_rejects_open_or_unsupported_inputs() -> None: @@ -195,7 +195,7 @@ def test_optimizer_rejects_group_turnover_conflict_instead_of_false_convergence( with pytest.raises( ValueError, - match="minimum required turnover 0.6 exceeds max_turnover 0.1", + match="infeasible: joint constraints", ): optimize_mean_variance( expected, @@ -302,7 +302,7 @@ def test_factor_constraint_validation_and_infeasibility_are_explicit() -> None: with pytest.raises(ValueError, match="provided together"): optimize_mean_variance(expected, covariance, factor_exposures=exposures) - with pytest.raises(ValueError, match="factor_bounds.*cannot be satisfied"): + with pytest.raises(ValueError, match="infeasible: joint constraints"): optimize_mean_variance( expected, covariance, @@ -355,5 +355,6 @@ def test_linear_cost_kink_uses_stable_proximal_convergence() -> None: ) assert result.converged - assert result.iterations == 1 + assert result.diagnostics["status"] == "optimal" + assert result.diagnostics["max_constraint_residual"] <= 1e-7 assert result.weights.to_dict() == pytest.approx(current.to_dict()) diff --git a/tests/test_solver_failure_contract.py b/tests/test_solver_failure_contract.py new file mode 100644 index 0000000..ddfced4 --- /dev/null +++ b/tests/test_solver_failure_contract.py @@ -0,0 +1,107 @@ +from decimal import Decimal + +import cvxpy as cp +import numpy as np +import pandas as pd +import pytest + +from quant_portfolio import OptimizationConstraints, RiskSnapshot +from quant_portfolio.executable import repair_integer_target +from quant_portfolio.solver import OptimizationError, RiskMatrix, solve_portfolio + + +def problem(**kwargs): + return solve_portfolio( + [0.05, 0.02], + RiskMatrix(covariance=np.diag([0.04, 0.09])), + [0.5, 0.5], + OptimizationConstraints(max_weight=0.8, max_turnover=0.2), + **kwargs, + ) + + +@pytest.mark.parametrize("failure", ["exception", "inaccurate"]) +def test_osqp_failure_retries_same_problem_without_publishing_inaccurate(monkeypatch, failure): + reference = problem() + original = cp.Problem.solve + calls = [] + + def solve(self, *args, **kwargs): + calls.append((id(self), kwargs["solver"])) + if kwargs["solver"] == "OSQP": + if failure == "exception": + raise cp.error.SolverError("injected numerical failure") + original(self, *args, **kwargs) + self._status = cp.OPTIMAL_INACCURATE + return self.value + return original(self, *args, **kwargs) + + monkeypatch.setattr(cp.Problem, "solve", solve) + actual = problem() + assert [name for _, name in calls] == ["OSQP", "CLARABEL"] + assert calls[0][0] == calls[1][0] + np.testing.assert_allclose(actual.weights, reference.weights, atol=1e-7) + assert actual.diagnostics["solver"] == "CLARABEL" + + +def test_two_inaccurate_results_are_not_publishable(monkeypatch): + original = cp.Problem.solve + + def solve(self, *args, **kwargs): + original(self, *args, **kwargs) + self._status = cp.OPTIMAL_INACCURATE + return self.value + + monkeypatch.setattr(cp.Problem, "solve", solve) + with pytest.raises(OptimizationError) as caught: + problem() + assert caught.value.status == "numerical_failure" + assert len(caught.value.diagnostics["attempts"]) == 2 + + +def test_zero_risk_positive_return_has_no_finite_sharpe(): + with pytest.raises(OptimizationError) as caught: + solve_portfolio( + [1, 2], + RiskMatrix(covariance=np.zeros((2, 2))), + [0.5, 0.5], + OptimizationConstraints(max_weight=1), + kind="max_sharpe", + ) + assert caught.value.status == "undefined_ratio" + + +def test_factor_snapshot_cannot_change_after_validation(): + snapshot = RiskSnapshot( + estimated_at="2026-01-01T00:00Z", + available_at="2026-01-01T00:00Z", + horizon_sessions=1, + annualization=252, + model_version="fixed", + input_digest="fixture", + covariance=pd.DataFrame(np.eye(2), index=list("AB"), columns=list("AB")), + ) + snapshot.covariance.iloc[0, 0] = 2 + with pytest.raises(ValueError, match="mutated"): + snapshot.matrix(list("AB"), "2026-01-02T00:00Z") + + +def test_known_unsellable_holdings_and_fractional_quantity_rejection(): + args = { + "weights": pd.Series({"A": 0.2, "B": 0.8}), + "current": {"A": 500}, + "prices": {"A": 10, "B": 10}, + "lots": {"A": 100, "B": 100}, + "cash": Decimal(5000), + "nav": Decimal(10000), + "quote_trade": lambda a, q: (-Decimal(q) * 10, Decimal(0)), + "constraint_residuals": lambda w: {}, + "max_invested": 1, + "max_sell": {"A": 0}, + "max_buy": {"A": 0, "B": 1000}, + } + result = repair_integer_target(**args) + assert result.quantities == {"A": 500, "B": 500} + assert result.cash == 0 + with pytest.raises(ValueError, match="integers"): + repair_integer_target(**{**args, "current": {"A": 500.5}}) diff --git a/tests/test_solver_v2.py b/tests/test_solver_v2.py new file mode 100644 index 0000000..b3b79cb --- /dev/null +++ b/tests/test_solver_v2.py @@ -0,0 +1,362 @@ +from decimal import Decimal + +import numpy as np +import pandas as pd +import pytest +from scipy.optimize import minimize + +from quant_portfolio.allocation_v2 import allocate_research +from quant_portfolio.executable import repair_integer_target +from quant_portfolio.optimization import ( + OptimizationConstraints, + _project_box_simplex, + research_allocation_weights, +) +from quant_portfolio.risk_snapshot import RiskSnapshot, mature_ic +from quant_portfolio.solver import OptimizationError, RiskMatrix, solve_portfolio + + +def test_projection_kkt_regression(): + y = np.array([-0.470931, 0.461829, -0.810164, 1.860959, -1.46164, -0.256543]) + x = _project_box_simplex(y, 0, 0.35) + assert x == pytest.approx([0.042806, 0.35, 0, 0.35, 0, 0.257194], abs=1e-6) + assert x.sum() == pytest.approx(1, abs=1e-12) + free = (x > 1e-12) & (x < 0.35 - 1e-12) + theta = (y - x)[free][0] + assert x == pytest.approx(np.clip(y - theta, 0, 0.35), abs=1e-12) + + +@pytest.mark.parametrize("case", [0, 6, 7]) +def test_legacy_reference_retains_fixed_six_asset_cases(case): + from quant_portfolio.optimization import optimize_mean_variance + + rng = np.random.default_rng(20261014) + for _ in range(case + 1): + a = rng.normal(size=(6, 6)) + cov = 0.02 * np.eye(6) + a @ a.T * 0.005 / 6 + mu, cost = rng.normal(0.05, 0.06, 6), rng.uniform(0.001, 0.005, 6) + assets = list("ABCDEF") + common = { + "current_weights": pd.Series(1 / 6, index=assets), + "linear_costs": pd.Series(cost, index=assets), + "constraints": OptimizationConstraints( + max_weight=0.35, + max_turnover=0.6, + group_by_asset={a: "G" if i < 3 else "H" for i, a in enumerate(assets)}, + group_caps={"G": 0.55}, + ), + "factor_exposures": pd.DataFrame({"style": np.linspace(-1, 1, 6)}, index=assets), + "factor_bounds": {"style": (-0.15, 0.15)}, + } + expected, covariance = ( + pd.Series(mu, index=assets), + pd.DataFrame(cov, index=assets, columns=assets), + ) + reference = optimize_mean_variance(expected, covariance, solver_backend="legacy", **common) + production = optimize_mean_variance(expected, covariance, **common) + assert reference.converged + assert reference.objective == pytest.approx(production.objective, abs=1e-8) + + +@pytest.mark.parametrize( + "mode", + [ + "equal", + "inverse_vol", + "risk_parity", + "hrp", + "cost_aware", + "max_sharpe", + "max_diversification", + "cvar", + "evar", + "cdar", + ], +) +def test_v2_all_methods_solve_real_current_joint_constraints(mode): + rng = np.random.default_rng(9) + assets = list("ABCD") + history = pd.DataFrame( + rng.normal(0.001, 0.01, (80, 4)), + index=pd.date_range("2025-01-01", periods=80, tz="UTC"), + columns=assets, + ) + uses = mode in {"cost_aware", "max_sharpe", "cvar", "evar", "cdar"} + result = allocate_research( + pd.Series([1.0, 2.0, 3.0, 4.0], index=assets), + history, + pd.Series(0.2, index=assets), + None, + invested_limit=0.8, + max_weight=0.35, + config={ + "protocol_version": 2, + "mode": mode, + "risk_objective": "absolute", + "forecast_horizon_sessions": 1, + "expected_return_model": "provided" if uses else "none", + "max_turnover": 0.2, + }, + expected_returns=pd.Series(0.001, index=assets), + decision_at=history.index[-1], + ) + assert result.weights.sum() == pytest.approx(0.8) + assert abs(result.weights - 0.2).sum() <= 0.2 + 1e-7 + assert result.weights.max() <= 0.35 + 1e-7 + assert result.diagnostics["compatibility_mode"] is False + if not uses and mode != "max_diversification": + assert "repaired_risk_contributions" in result.diagnostics or mode == "equal" + + +@pytest.mark.parametrize("case", range(200)) +def test_200_small_problems_against_independent_slsqp(case): + # Cases 0, 6 and 7 include the original false-infeasibility regressions. + rng = np.random.default_rng(20261014) + for _ in range(case + 1): + matrix = rng.normal(size=(6, 6)) + cov = 0.02 * np.eye(6) + matrix @ matrix.T * 0.005 / 6 + mu = rng.normal(0.05, 0.06, 6) + cost = rng.uniform(0.001, 0.005, 6) + current, factor = np.ones(6) / 6, np.linspace(-1, 1, 6) + settings = OptimizationConstraints( + max_weight=0.35, + max_turnover=0.6, + group_by_asset={i: "G" if i < 3 else "H" for i in range(6)}, + group_caps={"G": 0.55}, + ) + got = solve_portfolio( + mu, + RiskMatrix(covariance=cov), + current, + settings, + linear_costs=cost, + factor_exposures=factor[:, None], + factor_bounds={"style": (-0.15, 0.15)}, + ) + + def loss(z): + w, t = z[:6], z[6:] + return -mu @ w + 2.5 * w @ cov @ w + cost @ t + + def inequalities(z): + w, t = z[:6], z[6:] + return np.r_[ + t - w + current, + t + w - current, + 0.6 - t.sum(), + 0.55 - w[:3].sum(), + 0.15 - factor @ w, + 0.15 + factor @ w, + ] + + oracle = minimize( + loss, + np.r_[current, np.zeros(6)], + method="SLSQP", + jac=lambda z: np.r_[-mu + 5 * cov @ z[:6], cost], + bounds=[(0, 0.35)] * 6 + [(0, None)] * 6, + constraints=[ + {"type": "eq", "fun": lambda z: z[:6].sum() - 1}, + {"type": "ineq", "fun": inequalities}, + ], + options={"maxiter": 1000, "ftol": 1e-12}, + ) + assert oracle.success, oracle.message + assert inequalities(oracle.x).min() >= -1e-7 + assert abs(got.objective + oracle.fun) <= 1e-6 * max(abs(oracle.fun), 1e-6) + assert got.diagnostics["max_constraint_residual"] <= 1e-7 + + +def test_turnover_limited_sharpe_matches_grid_reference(): + assets = list("ABC") + mu = pd.Series([0.12, 0.08, 0.04], index=assets) + cov = pd.DataFrame( + [[0.04, 0.008, 0.002], [0.008, 0.01, 0.001], [0.002, 0.001, 0.09]], + index=assets, + columns=assets, + ) + w = research_allocation_weights( + mu, + None, + pd.Series([0.0, 0.0, 1.0], index=assets), + None, + invested_limit=1.0, + max_weight=1.0, + config={"mode": "max_sharpe", "max_turnover": 0.4}, + covariance_override=cov, + ) + assert w.to_numpy() == pytest.approx([0.2, 0, 0.8], abs=2e-6) + grid = np.array([[a, 0.2 - a, 0.8] for a in np.linspace(0, 0.2, 2001)]) + ratios = (grid @ mu.to_numpy()) / np.sqrt(np.einsum("ij,ij->i", grid @ cov.to_numpy(), grid)) + actual = float(mu @ w) / np.sqrt(float(w @ cov @ w)) + assert actual >= ratios.max() - 1e-8 + + +def risk_snapshot(assets, *, benchmark=None, exposures=None, available="2026-01-02T07:00Z"): + x = ( + pd.DataFrame({"market": np.ones(len(assets))}, index=assets) + if exposures is None + else exposures + ) + return RiskSnapshot( + estimated_at="2026-01-02T07:00Z", + available_at=available, + horizon_sessions=1, + annualization=252, + model_version="test-v1", + input_digest="test-input", + exposures=x, + factor_covariance=pd.DataFrame( + np.eye(x.shape[1]) * 0.01, index=x.columns, columns=x.columns + ), + specific_variance=pd.Series(0.04, index=assets), + benchmark=benchmark, + ) + + +def test_factor_dense_risk_and_gradient_agree_without_materialization(): + rng = np.random.default_rng(4) + x, a = rng.normal(size=(10, 3)), rng.normal(size=(3, 3)) + f, d = a @ a.T, rng.uniform(0.001, 0.1, 10) + factor = RiskMatrix(exposures=x, factor_covariance=f, specific_variance=d) + dense = RiskMatrix(covariance=x @ f @ x.T + np.diag(d)) + w = rng.uniform(0, 1, 10) + assert factor.variance(w) == pytest.approx(dense.variance(w), rel=1e-12) + assert factor.diagonal() == pytest.approx(dense.diagonal()) + eps = 1e-5 + gradient = np.array( + [ + (factor.variance(w + np.eye(10)[i] * eps) - factor.variance(w - np.eye(10)[i] * eps)) + / (2 * eps) + for i in range(10) + ] + ) + assert gradient == pytest.approx(2 * dense.covariance @ w, rel=1e-8) + + +@pytest.mark.parametrize("budget", [1.0, 0.8]) +def test_partial_cash_and_pool_external_benchmark_survive(budget): + assets = list("ABCD") + benchmark = pd.Series([0.25, 0.25, 0.25, 0.25], index=assets) + snapshot = risk_snapshot(assets, benchmark=benchmark) + result = allocate_research( + pd.Series({"A": 1.0, "B": 2.0, "C": 3.0}), + None, + pd.Series({"D": 0.1}), + None, + invested_limit=budget, + max_weight=0.7, + risk_snapshot=snapshot, + decision_at="2026-01-02T07:00Z", + expected_returns=pd.Series({"A": 0.001, "B": 0.002, "C": 0.003}), + config={ + "protocol_version": 2, + "mode": "cost_aware", + "risk_objective": "active", + "forecast_horizon_sessions": 1, + "expected_return_model": "provided", + "max_turnover": 2.0, + "max_tracking_error": 0.09, + }, + ) + w = result.weights + assert list(w.index) == assets + assert w.D == pytest.approx(0, abs=1e-8) + assert w.sum() == pytest.approx(budget) + risk = snapshot.matrix(assets, "2026-01-02T07:00Z") + assert result.diagnostics["annualized_tracking_error"] == pytest.approx( + np.sqrt(risk.variance(w - benchmark)) + ) + assert result.diagnostics["annualized_tracking_error"] <= 0.09 + 1e-7 + assert result.diagnostics["risk_snapshot"]["model_version"] == "test-v1" + + +def test_future_risk_rejected_and_no_positive_sharpe_explicit(): + snapshot = risk_snapshot(list("AB"), available="2026-01-03T07:00Z") + with pytest.raises(ValueError, match="unavailable"): + snapshot.matrix(list("AB"), "2026-01-02T07:00Z") + with pytest.raises(OptimizationError) as error: + solve_portfolio( + [-1.0, -2.0], + RiskMatrix(covariance=np.eye(2)), + [0.5, 0.5], + OptimizationConstraints(max_weight=1), + kind="max_sharpe", + ) + assert error.value.status == "no_positive_excess_return" + + +def test_mature_ic_ignores_future_labels_and_revisions(): + rows = [] + dates = pd.date_range("2026-01-01", periods=8, tz="UTC") + for i in range(6): + for j in range(3): + rows.append( + { + "signal_at": dates[i], + "label_end": dates[i + 2], + "available_at": dates[i + 2], + "asset": str(j), + "score": j, + "return_": j * 0.01, + "horizon_sessions": 2, + } + ) + labels = pd.DataFrame(rows).rename(columns={"return_": "return"}) + before = mature_ic(labels, decision_at=dates[6], horizon_sessions=2, lookback=10, min_periods=3) + labels.loc[labels.available_at > dates[6], "return"] *= -100 + revisions = labels.loc[labels.signal_at == dates[0]].copy() + revisions["available_at"], revisions["return"] = dates[7], -999 + after = mature_ic( + pd.concat([labels, revisions]), + decision_at=dates[6], + horizon_sessions=2, + lookback=10, + min_periods=3, + ) + assert before == after + assert before[0] == pytest.approx(1) + assert before[1]["ic_periods"] == 3 + + +def test_integer_repair_accounts_for_floor_fee_and_rechecks_risk(): + def quote(asset, delta): + fee = Decimal(5) + return -Decimal(delta) * 10 - fee, fee + + def residuals(w): + return {"A_cap": w.A - 0.51, "B_floor": 0.45 - w.B} + + got = repair_integer_target( + pd.Series({"A": 0.5, "B": 0.5}), + current={}, + prices={"A": Decimal(10), "B": Decimal(10)}, + lots={"A": 100, "B": 100}, + cash=Decimal(10000), + nav=Decimal(10000), + quote_trade=quote, + constraint_residuals=residuals, + max_invested=1.0, + ) + assert got.cash >= 0 + assert got.quantities == {"A": 400, "B": 500} + assert got.costs == Decimal(10) + assert got.cash == Decimal(990) + assert max(residuals(got.weights).values()) <= 1e-7 + assert got.diagnostics["integer_global_optimum"] is False + + +def test_unsuccessful_integer_search_does_not_claim_infeasibility(): + with pytest.raises(OptimizationError) as exc: + repair_integer_target( + pd.Series({"A": 1.0}), + current={}, + prices={"A": Decimal(10)}, + lots={"A": 100}, + cash=Decimal(100), + nav=Decimal(100), + quote_trade=lambda a, q: (-Decimal(q) * 10, Decimal(0)), + constraint_residuals=lambda w: {"lower": 0.9 - w.A}, + max_invested=1.0, + ) + assert exc.value.status == "repair_failed" From 6ac45f4cefe640e72d57e7d7dd33d0c3d83a6054 Mon Sep 17 00:00:00 2001 From: PureSaber Date: Wed, 14 Oct 2026 23:10:48 +0800 Subject: [PATCH 2/6] Reconcile unpaid dividends without spending receivables --- src/quant_portfolio/executable.py | 9 ++++++++- tests/test_solver_v2.py | 21 +++++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/src/quant_portfolio/executable.py b/src/quant_portfolio/executable.py index 28ede05..90a8b41 100644 --- a/src/quant_portfolio/executable.py +++ b/src/quant_portfolio/executable.py @@ -32,6 +32,7 @@ def repair_integer_target( quote_trade, constraint_residuals, max_invested, + non_cash_assets=Decimal(0), max_buy=None, max_sell=None, timeout_seconds=30, @@ -41,6 +42,8 @@ def repair_integer_target( It is supplied by the execution adapter so rounding and fee rules have one owner. Residuals must be independently computed in post-fee NAV weights. + non_cash_assets is independently evidenced account value (e.g. unpaid + dividends): it contributes to NAV, but can never fund a purchase. Current quantities can contain odd-lot remainders; only board-lot trades are proposed. No unobserved future market status is used. """ @@ -55,6 +58,9 @@ def repair_integer_target( if not weights.index.is_unique or not np.isfinite(weights).all() or (weights < -1e-7).any(): raise ValueError("invalid continuous weights") nav, cash = Decimal(str(nav)), Decimal(str(cash)) + non_cash_assets = Decimal(str(non_cash_assets)) + if not non_cash_assets.is_finite() or non_cash_assets < 0: + raise ValueError("invalid non-cash account assets") if not nav.is_finite() or nav <= 0 or not cash.is_finite() or cash < 0: raise ValueError("invalid account cash or NAV") n = len(assets) @@ -105,7 +111,7 @@ def evaluate(q): return (float("inf"), float("inf")), None invested = sum((Decimal(int(q[i])) * marks[i] for i in range(n)), Decimal(0)) # Ledger conservation is exact at the supplied marks and quote precision. - if balance + invested != post_nav: + if balance + invested + non_cash_assets != post_nav: raise ValueError("execution quote and account marks do not reconcile") w = pd.Series( [float(Decimal(int(q[i])) * marks[i] / post_nav) for i in range(n)], index=assets @@ -215,6 +221,7 @@ def candidates(q): "constraint_residuals": residuals, "seconds": perf_counter() - started, "cash": str(balance), + "non_cash_assets": str(non_cash_assets), "costs": str(costs), "post_fee_nav": str(nav - costs), "budget_shortfall": float(max_invested) diff --git a/tests/test_solver_v2.py b/tests/test_solver_v2.py index b3b79cb..2bfd127 100644 --- a/tests/test_solver_v2.py +++ b/tests/test_solver_v2.py @@ -360,3 +360,24 @@ def test_unsuccessful_integer_search_does_not_claim_infeasibility(): max_invested=1.0, ) assert exc.value.status == "repair_failed" + + +def test_unpaid_dividend_counts_in_nav_but_cannot_fund_purchases(): + kwargs = { + "current": {}, + "prices": {"A": Decimal(10)}, + "lots": {"A": 100}, + "cash": Decimal(900), + "nav": Decimal(1100), + "non_cash_assets": Decimal(200), + "quote_trade": lambda a, q: (-Decimal(q) * 10, Decimal(0)), + "constraint_residuals": lambda w: {}, + "max_invested": 1.0, + } + result = repair_integer_target(pd.Series({"A": 1.0}), **kwargs) + assert result.quantities == {} + assert result.cash == Decimal(900) + assert result.diagnostics["non_cash_assets"] == "200" + kwargs["non_cash_assets"] = Decimal(199) + with pytest.raises(ValueError, match="do not reconcile"): + repair_integer_target(pd.Series({"A": 1.0}), **kwargs) From 67cae5523f208fb835caa4d79ff358c1c140dd67 Mon Sep 17 00:00:00 2001 From: PureSaber Date: Wed, 14 Oct 2026 23:29:54 +0800 Subject: [PATCH 3/6] Preserve fixed integer targets throughout local repair --- src/quant_portfolio/executable.py | 12 +++++++++++- tests/test_solver_v2.py | 17 +++++++++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/src/quant_portfolio/executable.py b/src/quant_portfolio/executable.py index 90a8b41..193f748 100644 --- a/src/quant_portfolio/executable.py +++ b/src/quant_portfolio/executable.py @@ -33,6 +33,7 @@ def repair_integer_target( constraint_residuals, max_invested, non_cash_assets=Decimal(0), + fixed_quantities=None, max_buy=None, max_sell=None, timeout_seconds=30, @@ -66,7 +67,10 @@ def repair_integer_target( n = len(assets) if n == 0: raise ValueError("at least one executable asset is required") - for mapping in (current, lots, max_buy or {}, max_sell or {}): + fixed_quantities = fixed_quantities or {} + if set(fixed_quantities) - set(assets): + raise ValueError("fixed quantities contain an unknown asset") + for mapping in (current, lots, max_buy or {}, max_sell or {}, fixed_quantities): if any(isinstance(v, bool) or not np.isfinite(v) or int(v) != v for v in mapping.values()): raise ValueError("quantities and board lots must be integers") quantities = np.array([current.get(a, 0) for a in assets], dtype=np.int64) @@ -88,6 +92,12 @@ def repair_integer_target( raise ValueError("invalid tradable quantity") minimum = quantities - (sell // steps) * steps maximum = quantities + (buy // steps) * steps + for i, asset in enumerate(assets): + if asset in fixed_quantities: + fixed = int(fixed_quantities[asset]) + if not minimum[i] <= fixed <= maximum[i] or (fixed - quantities[i]) % steps[i]: + raise ValueError("fixed target violates known quantity or board-lot constraints") + minimum[i] = maximum[i] = fixed quote_cache = {} def evaluate(q): diff --git a/tests/test_solver_v2.py b/tests/test_solver_v2.py index 2bfd127..ed92a64 100644 --- a/tests/test_solver_v2.py +++ b/tests/test_solver_v2.py @@ -381,3 +381,20 @@ def test_unpaid_dividend_counts_in_nav_but_cannot_fund_purchases(): kwargs["non_cash_assets"] = Decimal(199) with pytest.raises(ValueError, match="do not reconcile"): repair_integer_target(pd.Series({"A": 1.0}), **kwargs) + + +def test_fixed_exit_target_cannot_be_undone_by_integer_distance_repair(): + result = repair_integer_target( + pd.Series({"A": 0.0, "B": 0.5}), + current={"A": 100}, + prices={"A": Decimal(10), "B": Decimal(100)}, + lots={"A": 100, "B": 100}, + cash=Decimal(1000), + nav=Decimal(2000), + quote_trade=lambda a, q: (-Decimal(q) * (10 if a == "A" else 100), Decimal(0)), + constraint_residuals=lambda w: {}, + max_invested=0.5, + fixed_quantities={"A": 0}, + ) + assert result.quantities == {} + assert result.cash == Decimal(2000) From 8db176e6f92210cc43e0c6043ac05e1f53ceb6d1 Mon Sep 17 00:00:00 2001 From: PureSaber Date: Wed, 14 Oct 2026 23:31:04 +0800 Subject: [PATCH 4/6] Start repair from positions clipped to known quantity constraints --- src/quant_portfolio/executable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/quant_portfolio/executable.py b/src/quant_portfolio/executable.py index 193f748..758b0d9 100644 --- a/src/quant_portfolio/executable.py +++ b/src/quant_portfolio/executable.py @@ -148,7 +148,7 @@ def evaluate(q): # Preserve the account's odd-lot remainder and enforce known quantity limits. floor = quantities + np.floor_divide(desired - quantities, steps) * steps floor = np.minimum(np.maximum(floor, minimum), maximum) - starts = [floor, quantities.copy()] + starts = [floor, np.minimum(np.maximum(quantities, minimum), maximum)] best = None attempts = [] for seed in starts: From 7cba8409559fb1625e88b15b983a446e2e6e61d6 Mon Sep 17 00:00:00 2001 From: PureSaber Date: Wed, 14 Oct 2026 23:39:46 +0800 Subject: [PATCH 5/6] Preserve absolute risk semantics when adapting legacy benchmark configurations --- src/quant_portfolio/allocation_v2.py | 11 ++++++----- tests/test_solver_v2.py | 22 ++++++++++++++++++++++ 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/src/quant_portfolio/allocation_v2.py b/src/quant_portfolio/allocation_v2.py index b0762c5..aa61832 100644 --- a/src/quant_portfolio/allocation_v2.py +++ b/src/quant_portfolio/allocation_v2.py @@ -59,7 +59,7 @@ def allocate_research( from quant_portfolio.optimization import ( OptimizationConstraints, _research_covariance, - _research_expected_returns, + ic_vol_expected_returns, validate_research_allocation, ) @@ -166,7 +166,7 @@ def allocate_research( diag = { "protocol_version": 2 if v2 else 1, "compatibility_mode": not v2, - "risk_objective": cfg.get("risk_objective", "active" if wb is not None else "absolute"), + "risk_objective": cfg.get("risk_objective", "absolute"), "risk_usage": METHODS[mode]["risk_usage"], "asset_union": assets, "invested_budget": invested_limit, @@ -223,9 +223,10 @@ def allocate_research( raise ValueError( "Black-Litterman requires an explicit views API and full market weights" ) - covariance = pd.DataFrame(np.diag(risk.diagonal()), index=assets, columns=assets) - mu = _research_expected_returns( - cfg, scores.reindex(assets, fill_value=0), covariance, market_weights + mu = ic_vol_expected_returns( + scores.reindex(assets, fill_value=0), + pd.Series(np.sqrt(risk.diagonal()), index=assets), + float(cfg["information_coefficient"]), ).to_numpy() else: mu = scores.reindex(assets, fill_value=0).to_numpy() diff --git a/tests/test_solver_v2.py b/tests/test_solver_v2.py index ed92a64..617f964 100644 --- a/tests/test_solver_v2.py +++ b/tests/test_solver_v2.py @@ -234,6 +234,28 @@ def test_factor_dense_risk_and_gradient_agree_without_materialization(): assert gradient == pytest.approx(2 * dense.covariance @ w, rel=1e-8) +@pytest.mark.parametrize("return_model", ["score", "ic_vol"]) +def test_legacy_benchmark_does_not_silently_change_absolute_risk_objective(return_model): + assets = list("AB") + common = { + "scores": pd.Series([0.001, 0.002], index=assets), + "trailing_returns": None, + "current_weights": pd.Series([0.5, 0.5], index=assets), + "linear_costs": None, + "invested_limit": 1.0, + "max_weight": 1.0, + "decision_at": "2026-01-02T07:00Z", + "config": {"mode": "cost_aware", "expected_return_model": return_model}, + } + without = allocate_research(**common, risk_snapshot=risk_snapshot(assets)) + with_benchmark = allocate_research( + **common, risk_snapshot=risk_snapshot(assets, benchmark=pd.Series([0.9, 0.1], index=assets)) + ) + assert with_benchmark.diagnostics["risk_objective"] == "absolute" + assert with_benchmark.diagnostics["compatibility_mode"] is True + assert with_benchmark.weights.to_numpy() == pytest.approx(without.weights.to_numpy(), abs=1e-7) + + @pytest.mark.parametrize("budget", [1.0, 0.8]) def test_partial_cash_and_pool_external_benchmark_survive(budget): assets = list("ABCD") From da7a6203b7df42d9e92697986cb43dad180e3696 Mon Sep 17 00:00:00 2001 From: PureSaber Date: Thu, 15 Oct 2026 00:05:41 +0800 Subject: [PATCH 6/6] Preserve required Linux check names alongside Windows matrix --- .github/workflows/ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 779881c..a212ffd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,6 +11,7 @@ permissions: jobs: test: + name: test (${{ matrix.python-version }})${{ matrix.os == 'windows-2025' && ' Windows' || '' }} runs-on: ${{ matrix.os }} strategy: matrix: