"""Sampling noise of the mean-CVaR MILP under different scenario generators (GBM model). Reference: the exact Gaussian problem max mu'x s.t. -mu'x + k_beta * sqrt(x'Sigma x) <= Lambda, caps, budget, integer shares -- solved with an outer-approximation (Kelley cutting planes) on the same MILP (no scenarios at all). Under a normal model this is the mean-variance solution. """ import sys, time, json import numpy as np, pandas as pd from scipy.stats import norm, qmc from scipy.optimize import milp, LinearConstraint, Bounds import scipy.sparse as sps from quantcore.config import get_settings from quantcore.data.market_data import MarketDataService from quantcore.optimization.service import OptimizationService, OptimizationRequest from quantcore.optimization.problem import MeanCVaRProblem, ProblemSpec, OptimizationConstraints, TransactionCosts from quantcore.optimization.solvers import HighsSolver from quantcore.models.covariance import CovarianceModel, symmetric_inv_sqrt N_ASSETS = int(sys.argv[1]) if len(sys.argv) > 1 else 30 S = int(sys.argv[2]) if len(sys.argv) > 2 else 2048 N_SEEDS = int(sys.argv[3]) if len(sys.argv) > 3 else 8 INTEGER = (sys.argv[4] == "int") if len(sys.argv) > 4 else False KINDS = sys.argv[5].split(",") if len(sys.argv) > 5 else ["mc", "antithetic", "mc+mm", "lhs", "sobol", "sobol+mm"] BETA = float(sys.argv[6]) if len(sys.argv) > 6 else 0.95 W = 1_000_000.0 LAMBDA = 0.05 * W H = 5 settings = get_settings() md = MarketDataService(settings) svc = OptimizationService(settings, md) req = OptimizationRequest(max_assets=N_ASSETS, scenario_model="gbm", n_scenarios=S, solver="highs", cash=W) symbols = svc.candidate_symbols(req) mu_ann, _, _ = svc.expected_returns(req, symbols) r = md.returns(symbols, currency="base", freq="daily").dropna() cov_model = CovarianceModel.from_dict(req.covariance) cov_d = cov_model.estimate(r).values mu_h = (mu_ann / 252).reindex(symbols).values * H cov_h = cov_d * H prices = md.latest_prices(symbols, currency="base").reindex(symbols).values sectors = md.sectors(symbols).fillna("Unknown").tolist() N = len(symbols) k_beta = norm.pdf(norm.ppf(BETA)) / (1 - BETA) print(f"N={N} assets, S={S} scenarios, horizon {H} days, beta={BETA}, k_beta={k_beta:.4f}") print("weekly vol range (%%): %.2f .. %.2f" % (100 * np.sqrt(np.diag(cov_h)).min(), 100 * np.sqrt(np.diag(cov_h)).max())) # eigen root (PCA ordering) and Cholesky root vals, vecs = np.linalg.eigh(cov_h) order = np.argsort(vals)[::-1] A_pca = vecs[:, order] * np.sqrt(vals[order]) # columns = principal directions, largest first A_chol = np.linalg.cholesky(cov_h) def moment_match(z): z = z - z.mean(axis=0) c = z.T @ z / len(z) return z @ symmetric_inv_sqrt(c) def base_sample(kind, seed): rng = np.random.default_rng(seed) if kind == "mc": return rng.standard_normal((S, N)), A_chol if kind == "antithetic": z = rng.standard_normal((S // 2, N)) return np.vstack([z, -z]), A_chol if kind == "mc+mm": return moment_match(rng.standard_normal((S, N))), A_chol if kind == "lhs": u = qmc.LatinHypercube(d=N, scramble=True, seed=seed).random(S) return norm.ppf(u), A_pca if kind == "sobol": u = qmc.Sobol(d=N, scramble=True, seed=seed).random(S) return norm.ppf(u), A_pca if kind == "sobol+mm": u = qmc.Sobol(d=N, scramble=True, seed=seed).random(S) return moment_match(norm.ppf(u)), A_pca raise ValueError(kind) def scenarios(kind, seed): z, a = base_sample(kind, seed) return mu_h + z @ a.T # S x N horizon log-returns (used linearly: exact Gaussian P&L) def solve_milp(scen, integer): spec = ProblemSpec(symbols=symbols, prices=prices, scenarios=scen, sectors=sectors, cash=W, holdings=np.zeros(N), confidence=BETA, cvar_limit=LAMBDA, objective="max_return", constraints=OptimizationConstraints(), costs=TransactionCosts(proportional_buy=0.0, proportional_sell=0.0)) prob = MeanCVaRProblem(spec) model = prob.build() if not integer: model.var_types = ["C"] * len(model.var_types) sol = HighsSolver().solve(model, time_limit=300, mip_gap=1e-5) n = np.rint(sol.x[prob.i_n]) if integer else sol.x[prob.i_n] pnl = scen @ (n * prices) from quantcore.measurement.risk import scenario_cvar _, cvar_in = scenario_cvar(-pnl, BETA) return n, float(pnl.mean()), cvar_in, sol def true_stats(n): v = n * prices mean = v @ mu_h sd = np.sqrt(v @ cov_h @ v) return mean, -mean + k_beta * sd, sd # ---- exact Gaussian reference by outer approximation on the same integer feasible set ---------------- def exact_reference(tol=1e-6, max_iter=200, integer=True): ub = np.floor(np.minimum(0.10, 1.0) * W / prices) c = -(prices * mu_h) # minimise -E[P&L] rows = [prices] # budget lb = [-np.inf]; ubs = [W] for sec in sorted(set(sectors)): # sector caps (defaults of OptimizationConstraints) mask = np.array([s == sec for s in sectors], dtype=float) rows.append(prices * mask); lb.append(-np.inf); ubs.append(OptimizationConstraints().max_sector_weight * W) cuts = [] n_prev = None for it in range(max_iter): A = np.vstack(rows + [cut[0] for cut in cuts]) lo = np.array(lb + [-np.inf] * len(cuts)); hi = np.array(ubs + [cut[1] for cut in cuts]) res = milp(c, constraints=LinearConstraint(A, lo, hi), integrality=np.ones(N) if integer else np.zeros(N), bounds=Bounds(np.zeros(N), ub), options={"mip_rel_gap": 1e-9}) n = np.round(res.x) if integer else res.x mean, cvar, sd = true_stats(n) if cvar <= LAMBDA + tol * LAMBDA: return n, mean, cvar, it + 1 v = n * prices grad = -prices * mu_h + k_beta * prices * (cov_h @ v) / sd # d CVaR / d n # cut: cvar(n_k) + grad'(n - n_k) <= Lambda -> grad'n <= Lambda - cvar(n_k) + grad'n_k cuts.append((grad, LAMBDA - cvar + grad @ n)) raise RuntimeError("no convergence") t0 = time.time() n_ref, mean_ref, cvar_ref, iters = exact_reference() x_rel, mean_rel, cvar_rel, iters_rel = exact_reference(integer=False) print(f"\nexact Gaussian reference (integer): E[P&L]={mean_ref:.2f} CHF, CVaR={cvar_ref:.2f}, " f"{int((n_ref > 0).sum())} positions, {iters} cuts, {time.time()-t0:.1f}s") print(f"continuous relaxation: E[P&L]={mean_rel:.2f} CHF, CVaR={cvar_rel:.2f}") results = {} for kind in KINDS: rows = [] t0 = time.time() for seed in range(N_SEEDS): n, epnl, cvar_in, sol = solve_milp(scenarios(kind, seed), INTEGER) mean_t, cvar_t, _ = true_stats(n) rows.append({"seed": seed, "epnl_in": epnl, "cvar_in": cvar_in, "epnl_true": mean_t, "cvar_true": cvar_t, "n_pos": int((n > 0).sum()), "shares": n, "status": sol.status}) df = pd.DataFrame(rows) shares = np.vstack(df["shares"].values) # dissimilarity between seeds: mean over pairs of (value traded to go from one to the other) / W pair_turn = [] for i in range(len(shares)): for j in range(i + 1, len(shares)): pair_turn.append(np.abs(shares[i] - shares[j]) @ prices / W) dist_ref = np.abs(shares - (n_ref if INTEGER else x_rel)) @ prices / W results[kind] = { "E[P&L] in-sample mean": df.epnl_in.mean(), "E[P&L] in-sample sd": df.epnl_in.std(), "E[P&L] true mean": df.epnl_true.mean(), "E[P&L] true sd": df.epnl_true.std(), "loss vs ref (CHF)": (mean_ref if INTEGER else mean_rel) - df.epnl_true.mean(), "true CVaR mean": df.cvar_true.mean(), "true CVaR sd": df.cvar_true.std(), "CVaR excess over limit (%)": 100 * (df.cvar_true.mean() / LAMBDA - 1), "positions": df.n_pos.mean(), "turnover between seeds (%)": 100 * np.mean(pair_turn), "turnover to ref (%)": 100 * dist_ref.mean(), "seconds": time.time() - t0, } print(f"{kind:10s} done in {time.time()-t0:.0f}s statuses={df.status.unique().tolist()}") out = pd.DataFrame(results).T pd.set_option("display.width", 250); pd.set_option("display.float_format", lambda v: f"{v:10.2f}") print("\n%s problem; reference E[P&L] = %.2f CHF, CVaR limit = %.0f CHF (5%% of wealth), weekly" % ("INTEGER" if INTEGER else "CONTINUOUS (weights)", mean_ref if INTEGER else mean_rel, LAMBDA)) print(out.T.to_string()) out.to_csv(f"{'int' if INTEGER else 'cont'}_scenario_noise_N{N}_S{S}.csv")