Portfolio_Risk_Management/tests/test_large_scale.py

"""The large-scale driver (relaxation, tail rows, screening, round and repair) against the full MILP."""

import numpy as np
import pytest

from quantcore.config import Settings
from quantcore.models import sampling as smp
from quantcore.optimization.large_scale import (
    exposure_and_cost,
    relative_gap,
    relaxed,
    solve_large_scale,
)
from quantcore.optimization.problem import MeanCVaRProblem, OptimizationConstraints, ProblemSpec
from quantcore.optimization.rowgen import solve_by_row_generation
from quantcore.optimization.service import OptimizationRequest, OptimizationService
from quantcore.optimization.solvers import HighsSolver

SETTINGS = Settings(cuopt_url="http://127.0.0.1:9", mip_rel_gap=1e-6)


def highs(model):
    return HighsSolver(SETTINGS).solve(model, time_limit=120)


def inputs(seed=0, n=12, S=2048):
    """A synthetic universe with moment-matched Sobol' scenarios and a fat-tailed scale mixture."""
    rng = np.random.default_rng(seed)
    symbols = [f"S{i}" for i in range(n)]
    prices = rng.uniform(20, 400, n)
    vol = rng.uniform(0.02, 0.06, n)
    corr = 0.4 * np.ones((n, n)) + 0.6 * np.eye(n)
    cov = corr * np.outer(vol, vol)
    mean = rng.uniform(0.0, 0.006, n)
    z = smp.moment_match(smp.standard_normal_sample(S, n, "sobol", seed))
    w = np.sqrt(4.0 / rng.chisquare(4.0, size=(S, 1))) * np.sqrt(0.5)  # Student-t(4) scale mixture
    scen = mean + (z @ smp.pca_root(cov).T) * w
    sectors = ["Tech", "Fin", "Health"] * (n // 3) + ["Tech"] * (n % 3)
    return symbols, prices, scen, sectors


def check_feasible(spec, problem, x, tol=1e-6):
    """The returned vector satisfies the full model, and the decoded portfolio the business limits."""
    m = problem.build()
    Ax = m.A @ x
    assert (Ax >= m.row_lb - tol).all() and (Ax <= m.row_ub + tol).all()
    assert (x >= m.var_lb - tol).all() and (x <= m.var_ub + tol).all()
    ints = [i for i, t in enumerate(m.var_types) if t == "I"]
    assert np.abs(x[ints] - np.rint(x[ints])).max() < 1e-9
    d = problem.decode(x)
    assert d["cash_end"] >= -tol
    if spec.cvar_limit is not None:
        assert d["cvar"] <= spec.cvar_limit * (1 + 1e-7)
    return d


def test_matches_full_milp_and_certifies_the_gap():
    # a portfolio of realistic size: share counts in the thousands, so rounding costs far less than the gap
    symbols, prices, scen, sectors = inputs(seed=2)
    spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=5_000_000.0,
                       cvar_limit=150_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
    full = MeanCVaRProblem(spec)
    res_full = highs(full.build())
    res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-4)
    assert res.ok and res.status == "Optimal" and info.certified and not info.escalated
    assert info.relaxation_converged and info.working_set < 0.5 * spec.n_scenarios
    assert info.integer_stage == "repair"
    d = check_feasible(spec, problem, res.x)
    d_full = full.decode(res_full.x)
    # the relaxation bounds the MILP optimum, and the rounded portfolio is within the gap of it
    assert info.bound >= d_full["expected_pnl"] - 1e-6
    assert d["expected_pnl"] >= d_full["expected_pnl"] - 1e-4 * abs(info.bound) - 1e-6
    assert d["expected_pnl"] <= d_full["expected_pnl"] * (1 + 1e-6) + 1e-6
    assert res.mip_gap == info.gap and info.gap <= 1e-4
    assert res.raw["large_scale"]["solver_calls"] == len(info.calls)


def test_relaxation_is_a_valid_bound_at_the_continuous_point():
    symbols, prices, scen, sectors = inputs(seed=5, n=9, S=1024)
    spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=50_000.0,
                       cvar_limit=2_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
    full = MeanCVaRProblem(spec)
    lp = highs(relaxed(full.build()))  # the LP relaxation on all rows
    res, _, info = solve_large_scale(spec, highs, mip_gap=1e-4)
    assert res.ok and info.relaxation_converged
    assert abs(info.bound - lp.objective) <= 1e-6 * max(1.0, abs(lp.objective))
    value, cost = exposure_and_cost(full, lp.x)
    assert value.shape == (9,) and cost > 0


def test_cardinality_screening_is_bounded():
    symbols, prices, scen, sectors = inputs(seed=7, n=30, S=1024)
    cons = OptimizationConstraints(max_weight=0.2, max_sector_weight=0.6, max_positions=5)
    spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=200_000.0,
                       cvar_limit=6_000.0, constraints=cons)
    full = MeanCVaRProblem(spec)
    res_full = highs(full.build())
    d_full = full.decode(res_full.x)
    res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-4, screen_min=10, screen_multiple=2)
    assert res.ok and info.screened_to == 10
    d = check_feasible(spec, problem, res.x)
    assert int((d["shares"] > 0).sum()) <= 5
    # the heuristic never beats the true optimum, and the reported gap covers the real one
    assert d["expected_pnl"] <= d_full["expected_pnl"] * (1 + 1e-6) + 1e-6
    true_gap = relative_gap(info.bound, d["expected_pnl"], maximize=True)
    assert info.gap >= true_gap - 1e-9
    assert info.gap >= relative_gap(d_full["expected_pnl"], d["expected_pnl"], maximize=True) - 1e-9


def test_min_cvar_objective_and_rebalancing():
    symbols, prices, scen, sectors = inputs(seed=3, n=9, S=1024)
    holdings = np.array([40, 0, 10, 0, 5, 0, 0, 8, 0], dtype=float)
    target = 0.001 * (20_000.0 + holdings @ prices)  # a third of what a fully invested portfolio expects
    spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=20_000.0,
                       holdings=holdings, objective="min_cvar", target_pnl=target,
                       constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
    full = MeanCVaRProblem(spec)
    res_full = highs(full.build())
    d_full = full.decode(res_full.x)
    res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-4)
    assert res.ok and info.relaxation_converged
    d = check_feasible(spec, problem, res.x)
    np.testing.assert_array_equal(d["shares"], holdings + d["buys"] - d["sells"])
    assert d["expected_pnl"] >= target - 1e-6
    assert info.bound <= d_full["cvar"] + 1e-6  # a lower bound on the minimal CVaR
    assert d["cvar"] <= d_full["cvar"] + 1e-4 * max(1.0, abs(info.bound)) + 1e-6


def test_escalates_to_the_integer_milp_when_not_certified():
    symbols, prices, scen, sectors = inputs(seed=11, n=12, S=1024)
    spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=100_000.0,
                       cvar_limit=3_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
    res_rg, problem_rg, _ = solve_by_row_generation(spec, highs)
    d_rg = problem_rg.decode(res_rg.x)
    res, problem, info = solve_large_scale(spec, highs, mip_gap=1e-12)  # nothing rounds this well
    assert res.ok and info.escalated and info.integer_stage == "milp"
    d = check_feasible(spec, problem, res.x)
    assert abs(d["expected_pnl"] - d_rg["expected_pnl"]) <= 1e-6 * abs(d_rg["expected_pnl"]) + 1e-6
    res_no, _, info_no = solve_large_scale(spec, highs, mip_gap=1e-12, escalate=False)
    assert res_no.ok and not info_no.escalated and res_no.status == "FeasibleFound" and not info_no.certified


def test_solver_selection_by_size():
    symbols, prices, scen, sectors = inputs(seed=4, n=9, S=512)
    spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=50_000.0,
                       cvar_limit=2_000.0, constraints=OptimizationConstraints(max_weight=0.3, max_sector_weight=0.6))
    seen = []

    def small(model):
        seen.append(("small", model.n_rows))
        return highs(model)

    def large(model):
        seen.append(("large", model.n_rows))
        return highs(model)

    res, _, info = solve_large_scale(spec, small, large, large_rows=100)
    assert res.ok
    assert all(w == ("large" if rows >= 100 else "small") for w, rows in seen)
    assert {c["stage"] for c in info.calls} <= {"relaxation", "repair", "milp"}
    with pytest.raises(ValueError):
        gspec = ProblemSpec(symbols=symbols[:3], prices=prices[:3], scenarios=None, sectors=sectors[:3], cash=1e4,
                            cvar_limit=500.0, gaussian_mean=np.zeros(3), gaussian_cov=np.eye(3) * 1e-4)
        solve_large_scale(gspec, highs)


BASE = {"cash": 500_000.0, "cvar_limit": 12_000.0, "solver": "highs", "horizon": "week", "scenario_model": "gbm",
        "n_scenarios": 2048, "importance_sampling": "off", "costs": {"proportional_buy": 0.0, "proportional_sell": 0.0},
        "constraints": {"max_weight": 0.15, "max_sector_weight": 0.4}}


def test_service_routes_the_strategy(settings, md):
    from tests.conftest import SYMBOLS

    svc = OptimizationService(settings, md)
    milp = svc.optimize(OptimizationRequest(symbols=SYMBOLS, strategy="milp", **BASE))
    ls = svc.optimize(OptimizationRequest(symbols=SYMBOLS, strategy="large_scale", **BASE))
    assert "large_scale" not in milp.solver_info and ls.solver_info["large_scale"]["certified"]
    assert ls.status == "Optimal" and ls.reproducible
    assert abs(ls.expected_pnl - milp.expected_pnl) <= 1e-4 * abs(milp.expected_pnl) + 1e-6
    assert ls.cvar <= 12_000.0 * (1 + 1e-7)
    assert OptimizationRequest(symbols=SYMBOLS, **BASE).uses_large_scale(12, 2048)  # the default path
    assert not OptimizationRequest(symbols=SYMBOLS, strategy="milp", **BASE).uses_large_scale(750, 65536)
    with pytest.raises(ValueError):
        OptimizationRequest(symbols=SYMBOLS, strategy="bogus", **BASE)