356017e689
strategy='auto' now routes every scenario-mode CVaR optimisation through quantcore.optimization.large_scale (LP relaxation on the tail rows, screening, round-and-repair, certified against the relaxation bound, escalation to the integer MILP otherwise); strategy='milp' is the explicit opt-out, and row_generation applies to that path only. When the escalated MILP proves the rounded portfolio optimal within its gap, the result is Optimal even where the relaxation bound is loose (small portfolios with a real integrality gap). UI: a Strategy selector and the driver's bound, gap, screened size, rows and solver calls in the provenance card. Deploy env examples document the new cuOpt settings (presolve off, dual simplex for LPs, 50 ms polling). The service-equivalence test asks for strategy='milp' where it checks the integer row-generation path. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
284 lines
18 KiB
Python
284 lines
18 KiB
Python
"""Service-level checks: the scenario optimiser reproduces the exact Gaussian (mean-variance-equivalent)
|
|
optimum when the scenarios are moment-matched Sobol' points, the drift convention is arithmetic, and
|
|
importance sampling / row generation leave the optimum unchanged."""
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
import pytest
|
|
|
|
from quantcore.optimization.problem import MeanCVaRProblem
|
|
from quantcore.optimization.service import OptimizationRequest, OptimizationService
|
|
from tests.conftest import SYMBOLS
|
|
|
|
BASE = dict(symbols=SYMBOLS, cash=500_000.0, cvar_limit=12_000.0, solver="highs", horizon="week",
|
|
costs={"proportional_buy": 0.0, "proportional_sell": 0.0}, importance_sampling="off",
|
|
constraints={"max_weight": 0.15, "max_sector_weight": 0.4})
|
|
|
|
|
|
@pytest.fixture(scope="module")
|
|
def svc(settings, md):
|
|
return OptimizationService(settings, md)
|
|
|
|
|
|
def gaussian_cvar_of(svc, req, result):
|
|
"""Closed-form CVaR of a proposal under the Gaussian model of the request."""
|
|
symbols = svc.candidate_symbols(req)
|
|
mu_ann, _, _ = svc.expected_returns(req, symbols)
|
|
spec, _ = svc.build_gaussian_spec(req, symbols, mu_ann)
|
|
shares = {p["symbol"]: p["shares"] for p in result.positions}
|
|
value = np.array([shares.get(s, 0) for s in symbols]) * spec.prices
|
|
return MeanCVaRProblem(spec).gaussian_risk(value, result.total_cost)
|
|
|
|
|
|
def test_sobol_moment_matched_scenarios_reproduce_the_exact_optimum(svc):
|
|
exact = svc.optimize(OptimizationRequest(scenario_model="gaussian_exact", **BASE))
|
|
assert exact.status == "Optimal" and exact.reproducible and exact.solver_info["outer_approximation"]["converged"]
|
|
scen = svc.optimize(OptimizationRequest(scenario_model="gbm", n_scenarios=4096, **BASE))
|
|
assert scen.status == "Optimal" and scen.model["scenarios"]["moment_matched"]
|
|
# The GBM scenarios are lognormal in simple returns while the exact mode is normal: at weekly
|
|
# volatilities of 4-5 % the lognormal left tail is ~2-3 % thinner, so the scenario optimum takes
|
|
# ~3 % more exposure. Within that model difference the two agree (the sharp test of the sampling
|
|
# itself, with Gaussian scenarios used linearly, is in test_rowgen_gaussian.py).
|
|
assert abs(scen.expected_pnl - exact.expected_pnl) <= 0.05 * exact.expected_pnl
|
|
_, cvar_true, _ = gaussian_cvar_of(svc, OptimizationRequest(scenario_model="gaussian_exact", **BASE), scen)
|
|
assert cvar_true <= BASE["cvar_limit"] * 1.05
|
|
# the expected horizon return reported per position is the model's (arithmetic) expectation
|
|
mu_ann, _, _ = svc.expected_returns(OptimizationRequest(**BASE), SYMBOLS)
|
|
for p in scen.positions:
|
|
assert abs(p["expected_return_horizon"] - mu_ann[p["symbol"]] / 252 * 5) < 1e-12
|
|
# plain Monte Carlo without moment matching: the sample mean disagrees with the model mean
|
|
mc = svc.optimize(OptimizationRequest(scenario_model="gbm", n_scenarios=4096, sampling="mc", moment_match=False, **BASE))
|
|
assert mc.status == "Optimal"
|
|
dev_mc = abs(mc.solver_info["objective"] - mc.expected_pnl) # objective uses the model mean either way
|
|
assert dev_mc < 1e-6 # (the objective is always the model expectation) ...
|
|
scen_mean_mc = svc.build_spec(OptimizationRequest(scenario_model="gbm", n_scenarios=4096, sampling="mc",
|
|
moment_match=False, **BASE), SYMBOLS,
|
|
svc.scenarios(OptimizationRequest(scenario_model="gbm", n_scenarios=4096, sampling="mc",
|
|
moment_match=False, **BASE), SYMBOLS, mu_ann))
|
|
# ... but the sample mean of a plain MC set is visibly off the model mean, a moment-matched one is not
|
|
err_mc = np.abs(scen_mean_mc.scenarios.mean(axis=0) - scen_mean_mc.expected_returns).max()
|
|
sobol_spec = svc.build_spec(OptimizationRequest(scenario_model="gbm", n_scenarios=4096, **BASE), SYMBOLS,
|
|
svc.scenarios(OptimizationRequest(scenario_model="gbm", n_scenarios=4096, **BASE), SYMBOLS, mu_ann))
|
|
err_sobol = np.abs(sobol_spec.scenarios.mean(axis=0) - sobol_spec.expected_returns).max()
|
|
assert err_sobol < 0.1 * err_mc
|
|
|
|
|
|
def test_importance_sampling_and_row_generation_keep_the_optimum(svc):
|
|
plain = svc.optimize(OptimizationRequest(scenario_model="gbm", n_scenarios=4096, **BASE))
|
|
tilted = svc.optimize(OptimizationRequest(scenario_model="gbm", n_scenarios=4096,
|
|
**{**BASE, "importance_sampling": "market"}))
|
|
assert tilted.status == "Optimal" and tilted.model["scenarios"]["weighted"]
|
|
assert tilted.model["scenarios"]["effective_sample_size"] > 2000
|
|
assert abs(tilted.expected_pnl - plain.expected_pnl) <= 0.02 * plain.expected_pnl
|
|
exact_req = OptimizationRequest(scenario_model="gaussian_exact", **BASE)
|
|
assert gaussian_cvar_of(svc, exact_req, tilted)[1] <= BASE["cvar_limit"] * 1.05
|
|
# the integer MILP path with tail-row generation (the default strategy is the large-scale driver)
|
|
rg = svc.optimize(OptimizationRequest(scenario_model="gbm", n_scenarios=4096,
|
|
**{**BASE, "strategy": "milp", "row_generation": "on"}))
|
|
assert rg.status == "Optimal" and rg.solver_info["row_generation"]["converged"]
|
|
assert rg.scenario_hash == plain.scenario_hash # same scenario set ...
|
|
assert abs(rg.expected_pnl - plain.expected_pnl) <= 1e-4 * plain.expected_pnl # ... same optimum
|
|
assert max(rg.solver_info["row_generation"]["working_set_sizes"]) < 4096
|
|
|
|
|
|
def test_empirical_model_options(svc):
|
|
req = OptimizationRequest(scenario_model="empirical", n_scenarios=2048, **{**BASE, "importance_sampling": "market"})
|
|
res = svc.optimize(req)
|
|
assert res.status == "Optimal" and res.model["scenarios"]["model"] == "empirical"
|
|
assert all(t > 0 for t in res.model["scenarios"]["importance_sampling"]["theta"])
|
|
assert res.model["scenarios"]["importance_sampling"]["direction"] == "market"
|
|
assert res.cvar <= BASE["cvar_limit"] * (1 + 1e-6)
|
|
bal = svc.optimize(OptimizationRequest(scenario_model="empirical", n_scenarios=1000, sampling="balanced", **BASE))
|
|
assert bal.status == "Optimal" and bal.model["scenarios"]["n_scenarios"] == 1000
|
|
with pytest.raises(ValueError):
|
|
OptimizationRequest(scenario_model="bogus", **BASE)
|
|
with pytest.raises(ValueError):
|
|
OptimizationRequest(**{**BASE, "importance_sampling": "maybe"})
|
|
|
|
|
|
def test_one_scenario_set_for_a_family_of_runs(svc):
|
|
"""Risk level, constraints, costs, objective, cash and holdings must not change the scenarios;
|
|
a different expected-return model must change only their drift, not their shocks."""
|
|
svc.clear_scenario_cache()
|
|
base = dict(BASE, importance_sampling="market", scenario_model="gbm", n_scenarios=2048)
|
|
ref = svc.optimize(OptimizationRequest(**base))
|
|
variants = {
|
|
"other CVaR limit": dict(base, cvar_limit=6_000.0),
|
|
"limit as a fraction": dict(base, cvar_limit=None, cvar_limit_fraction=0.02),
|
|
"confidence 99 %": dict(base, confidence=0.99, cvar_limit=20_000.0),
|
|
"tighter caps": dict(base, constraints={"max_weight": 0.10, "max_sector_weight": 0.3}),
|
|
"with costs": dict(base, costs={"proportional_buy": 0.002, "proportional_sell": 0.002}),
|
|
"min-CVaR objective": dict(base, objective="min_cvar", cvar_limit=None, target_return=0.001),
|
|
"utility objective": dict(base, objective="utility", cvar_limit=None, risk_aversion=0.2),
|
|
"more cash": dict(base, cash=2_000_000.0),
|
|
"rebalancing": dict(base, mode="rebalance", cash=50_000.0,
|
|
holdings={p["symbol"]: p["shares"] for p in ref.positions[:3]}),
|
|
}
|
|
for label, kwargs in variants.items():
|
|
res = svc.optimize(OptimizationRequest(**kwargs))
|
|
assert res.status == "Optimal", label
|
|
assert res.scenario_key["id"] == ref.scenario_key["id"], label
|
|
assert res.scenario_hash == ref.scenario_hash, label
|
|
assert res.solver_info["scenario_cache"]["reused"], label # and it was not regenerated
|
|
# a different expected-return model shifts the scenarios but leaves the shocks alone
|
|
symbols = svc.candidate_symbols(OptimizationRequest(**base))
|
|
req = OptimizationRequest(**base)
|
|
mu_a, _, _ = svc.expected_returns(req, symbols)
|
|
scen_a = svc.scenarios(req, symbols, mu_a)
|
|
req_views = OptimizationRequest(**dict(base, views=[{"symbol": SYMBOLS[0], "sentiment": 0.9, "confidence": 0.8}]))
|
|
mu_b, _, _ = svc.expected_returns(req_views, symbols)
|
|
scen_b = svc.scenarios(req_views, symbols, mu_b)
|
|
assert scen_b.hash != scen_a.hash and scen_b.shock_hash == scen_a.shock_hash
|
|
np.testing.assert_allclose(scen_b.shocks, scen_a.shocks, atol=1e-15)
|
|
assert (scen_b.expected_simple[0] > scen_a.expected_simple[0]) # the view moved the drift
|
|
# what *must* change the set: symbols, horizon, covariance model, seed, sampling, scenario count
|
|
for label, kwargs in {
|
|
"fewer symbols": dict(base, symbols=SYMBOLS[:8]),
|
|
"monthly horizon": dict(base, horizon="month"),
|
|
"other covariance": dict(base, covariance={"method": "rolling", "window": 250}),
|
|
"other seed": dict(base, seed=7),
|
|
"plain Monte Carlo": dict(base, sampling="mc"),
|
|
"more scenarios": dict(base, n_scenarios=4096),
|
|
"no tilt": dict(base, importance_sampling="off"),
|
|
}.items():
|
|
other = OptimizationRequest(**kwargs)
|
|
assert svc.scenario_key(other, svc.candidate_symbols(other))["id"] != ref.scenario_key["id"], label
|
|
# the adaptive tilt is honest about depending on the risk parameters
|
|
ad = OptimizationRequest(**dict(base, importance_sampling="adaptive"))
|
|
ad_other = OptimizationRequest(**dict(base, importance_sampling="adaptive", cvar_limit=6_000.0))
|
|
assert svc.scenario_key(ad, symbols)["id"] != svc.scenario_key(ad_other, symbols)["id"]
|
|
assert "adaptive_tilt" in svc.scenario_key(ad, symbols)
|
|
|
|
|
|
def test_tilt_is_driven_by_the_covariance_and_reports_its_alignment(svc):
|
|
"""The tilt shifts each asset by its covariance with the reference portfolio (not by its own
|
|
volatility), leaves the model correlations exactly intact, and the result says how well the
|
|
solution's loss lines up with the tilted direction."""
|
|
from scipy.stats import norm
|
|
|
|
from quantcore.models import sampling as smp
|
|
from quantcore.models.covariance import CovarianceModel
|
|
from quantcore.models.simulation import GBMSimulator, TiltSpec
|
|
from quantcore.optimization.problem import OptimizationConstraints, ProblemSpec, TransactionCosts
|
|
|
|
req = OptimizationRequest(**dict(BASE, scenario_model="gbm", n_scenarios=2048, importance_sampling="market"))
|
|
symbols = svc.candidate_symbols(req)
|
|
r, cov_model = svc._returns_and_cov_model(req, symbols)
|
|
h = req.horizon_periods
|
|
cov_h = cov_model.estimate(r).values * h
|
|
root = smp.pca_root(cov_h)
|
|
reference = svc.market_exposure(symbols) * 1e6
|
|
|
|
# 1. the shift of a tilted block is exactly -z_beta * Sigma v / sigma: proportional to each
|
|
# asset's covariance with the reference, i.e. to its beta, not to its own volatility
|
|
z = norm.ppf(0.95)
|
|
direction = smp.loss_direction(root, reference)
|
|
shift = z * (root @ direction)
|
|
sigma = np.sqrt(reference @ cov_h @ reference)
|
|
np.testing.assert_allclose(shift, -z * (cov_h @ reference) / sigma, atol=1e-12)
|
|
beta_i = (cov_h @ reference) / (reference @ cov_h @ reference)
|
|
assert np.corrcoef(shift, beta_i)[0, 1] < -0.999 # perfectly ordered by beta
|
|
if np.std(np.sqrt(np.diag(cov_h))) > 0: # ... and only loosely by own volatility
|
|
assert np.corrcoef(shift, np.sqrt(np.diag(cov_h)))[0, 1] > -0.999
|
|
|
|
# 2. the weighted tilted sample has exactly the model covariance (and correlations)
|
|
mu_ann, _, _ = svc.expected_returns(req, symbols)
|
|
sim = GBMSimulator(cov_model)
|
|
tilted = sim.simulate(r, n_scenarios=2048, horizon_periods=h, mu=mu_ann / 252, seed=1,
|
|
tilt=TiltSpec(exposure=reference))
|
|
w = tilted.weights
|
|
centred = tilted.log_returns - w @ tilted.log_returns
|
|
cov_t = (centred * w[:, None]).T @ centred
|
|
np.testing.assert_allclose(cov_t, cov_h, rtol=1e-8, atol=1e-14)
|
|
unweighted = np.cov(tilted.log_returns, rowvar=False)
|
|
assert np.abs(unweighted - cov_h).max() > 1e-6 # the raw tilted sample is *not* the model
|
|
|
|
# 3. alignment: 1 for the reference itself, negative for a portfolio built to oppose it
|
|
spec = ProblemSpec(scenarios=tilted.simple_returns(), probabilities=tilted.probabilities,
|
|
expected_returns=tilted.expected_simple, **svc._spec_common(req, symbols))
|
|
assert abs(svc.tilt_alignment(spec, reference, svc.market_exposure(symbols)) - 1) < 1e-9 # itself
|
|
opposed = np.zeros(len(symbols))
|
|
opposed[np.argsort(beta_i)[:2]] = 1e5 # the two assets least exposed to the reference
|
|
assert svc.tilt_alignment(spec, opposed, svc.market_exposure(symbols)) < 0.9
|
|
# a real run reports the figure and stays well above the warning threshold
|
|
res = svc.optimize(req)
|
|
assert res.solver_info["tilt_alignment"] > 0.5 and not res.warnings
|
|
|
|
|
|
def test_risk_drift_zero_keeps_the_limit_binding_when_the_forecast_is_large(md, settings):
|
|
"""Over a long horizon a big expected return shifts the whole tail into profit, and a CVaR limit
|
|
measured on the drifted scenarios stops constraining anything. The zero-drift convention measures
|
|
risk on the centred scenarios -- so the risk figure does not depend on the return forecast at all
|
|
-- while the objective still uses it."""
|
|
svc = OptimizationService(settings, md)
|
|
base = dict(symbols=SYMBOLS[:8], cash=1_000_000.0, horizon=252, freq="daily", confidence=0.99,
|
|
cvar_limit_fraction=0.12, scenario_model="gbm", n_scenarios=2048, solver="highs",
|
|
expected_returns="zero", row_generation="off", importance_sampling="off")
|
|
symbols = svc.candidate_symbols(OptimizationRequest(**base))
|
|
v = np.full(len(symbols), 125_000.0) # an equally weighted portfolio of the whole wealth
|
|
|
|
risk = {}
|
|
for mu_val in (0.8, 1.5): # 80 % and 150 % p.a. forecasts
|
|
for mode in ("model", "zero"):
|
|
req = OptimizationRequest(**{**base, "risk_drift": mode})
|
|
spec = svc.build_spec(req, symbols, svc.scenarios(req, symbols, pd.Series(mu_val, index=symbols)))
|
|
risk[(mu_val, mode)] = OptimizationService.portfolio_risk_on_scenarios(spec, v / spec.prices, 0.99)
|
|
if mode == "zero":
|
|
np.testing.assert_allclose(spec.mean_returns, svc.build_spec(
|
|
OptimizationRequest(**{**base, "risk_drift": "model"}), symbols,
|
|
svc.scenarios(OptimizationRequest(**base), symbols, pd.Series(mu_val, index=symbols))
|
|
).mean_returns) # the objective sees the forecast either way
|
|
|
|
# the zero-drift risk is the same whatever the forecast; the drifted one is not
|
|
assert risk[(0.8, "zero")]["cvar"] == pytest.approx(risk[(1.5, "zero")]["cvar"], rel=1e-12)
|
|
assert risk[(1.5, "model")]["cvar"] < 0 < risk[(1.5, "zero")]["cvar"] # the tail became a gain
|
|
assert risk[(0.8, "model")]["cvar"] < risk[(0.8, "zero")]["cvar"]
|
|
# the drift enters multiplicatively (e^{r+d} - 1), so the gap is well below mu * V
|
|
gap = risk[(1.5, "zero")]["cvar"] - risk[(1.5, "model")]["cvar"]
|
|
assert 0.3 * 1.5e6 < gap < 1.5e6
|
|
|
|
res = svc.optimize(OptimizationRequest(**{**base, "expected_returns": "historical"}))
|
|
if res.cvar <= 0: # the historical forecast is large enough to make the limit vacuous
|
|
assert any("constrains nothing" in w for w in res.warnings)
|
|
|
|
|
|
def test_an_all_cash_utility_solution_says_why(md, settings):
|
|
"""objective='utility' with a large lambda is optimal in cash -- which must not look like a bug."""
|
|
from quantcore.optimization.service import OptimizationRequest, OptimizationService
|
|
|
|
svc = OptimizationService(settings, md)
|
|
common = dict(cash=1_000_000.0, max_assets=6, horizon="week", scenario_model="gaussian_exact",
|
|
solver="highs", time_limit=120, objective="utility")
|
|
big = svc.optimize(OptimizationRequest(risk_aversion=50.0, **common))
|
|
assert big.status in ("Optimal", "FeasibleFound") and big.invested == 0
|
|
assert any("nothing was invested" in w and "risk_aversion" in w for w in big.warnings)
|
|
small = svc.optimize(OptimizationRequest(risk_aversion=0.0, **common))
|
|
assert small.invested > 0 and not any("nothing was invested" in w for w in small.warnings)
|
|
|
|
|
|
def test_a_weighted_cvar_objective_is_routed_to_highs_and_says_so(md, settings):
|
|
"""cuOpt aborts when the CVaR is the objective *and* the scenarios carry probabilities.
|
|
|
|
The routing is that narrow on purpose: the same objective with uniform weights solves on the GPU,
|
|
so only the weighted case may be diverted, and only under solver='auto'.
|
|
"""
|
|
from quantcore.optimization.service import OptimizationRequest, OptimizationService
|
|
|
|
svc = OptimizationService(settings, md)
|
|
common = dict(cash=1_000_000.0, max_assets=6, horizon="week", scenario_model="gbm",
|
|
n_scenarios=1024, solver="auto", time_limit=120, risk_drift="zero")
|
|
weighted = svc.optimize(OptimizationRequest(objective="min_cvar", target_return=0.0005,
|
|
importance_sampling="market", **common))
|
|
assert weighted.solver == "highs"
|
|
assert any("probability-weighted" in w for w in weighted.warnings)
|
|
|
|
# uniform weights: no diversion (cuOpt is unreachable in the tests, so it still lands on HiGHS --
|
|
# through the fallback, which says something different)
|
|
plain = svc.optimize(OptimizationRequest(objective="min_cvar", target_return=0.0005,
|
|
importance_sampling="off", **common))
|
|
assert not any("probability-weighted" in w for w in plain.warnings)
|
|
# the constrained form is never diverted either
|
|
other = svc.optimize(OptimizationRequest(objective="max_return", cvar_limit_fraction=0.05,
|
|
importance_sampling="market", **common))
|
|
assert not any("probability-weighted" in w for w in other.warnings)
|