"""Importance sampling for the mean-CVaR LP (weights): two-stage exponential tilt along the loss direction of a first-pass solution, defensive mixture with the nominal density, self-normalised weights, weighted moment matching, exact mean in the objective. Compared with Sobol+mm at equal S against the exact Gaussian optimum.""" import sys, time import numpy as np, pandas as pd import scipy.sparse as sps from scipy.optimize import milp, LinearConstraint, Bounds from scipy.stats import norm, qmc 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 OptimizationConstraints from quantcore.models.covariance import CovarianceModel, symmetric_inv_sqrt N_ASSETS = int(sys.argv[1]); SIZES = [int(x) for x in sys.argv[2].split(",")]; N_SEEDS = int(sys.argv[3]) BETAS = [float(x) for x in sys.argv[4].split(",")] W = 1_000_000.0; LAMBDA = 0.05 * W; H = 5; WMAX = 0.10; SMAX = 0.25 settings = get_settings(); md = MarketDataService(settings); svc = OptimizationService(settings, md) req = OptimizationRequest(max_assets=N_ASSETS, 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_h = CovarianceModel.from_dict(req.covariance).estimate(r).values * H mu_h = (mu_ann / 252).reindex(symbols).values * H prices = md.latest_prices(symbols, currency="base").reindex(symbols).values sectors = md.sectors(symbols).fillna("Unknown").tolist() N = len(symbols) vals, vecs = np.linalg.eigh(cov_h); o = np.argsort(vals)[::-1]; A = vecs[:, o] * np.sqrt(vals[o]) sector_rows = [(np.array([s == sec for s in sectors], dtype=float)) for sec in sorted(set(sectors))] def k_of(beta): return norm.pdf(norm.ppf(beta)) / (1 - beta) def true_stats(n, beta): v = n * prices; mean = v @ mu_h; sd = np.sqrt(v @ cov_h @ v) return mean, -mean + k_of(beta) * sd, sd def weighted_mm(z, p): m = p @ z zc = z - m c = (zc * p[:, None]).T @ zc return zc @ symmetric_inv_sqrt(c) def solve_lp(scen, p, beta): """max mu_h'(n*P) s.t. RU-CVaR with scenario probabilities p <= Lambda, budget, caps, sectors; n continuous.""" S = len(scen) nv = N + 1 + S c = np.zeros(nv); c[:N] = -(prices * mu_h) # exact model mean in the objective G = scen * prices[None, :] rows = [np.repeat(np.arange(S), N), np.arange(S), np.arange(S)] cols = [np.tile(np.arange(N), S), np.full(S, N), N + 1 + np.arange(S)] vals_ = [G.ravel(), np.ones(S), np.ones(S)] lb = [np.zeros(S)]; ub = [np.full(S, np.inf)] ri = S rows.append(np.concatenate([[ri], np.full(S, ri)])); cols.append(np.concatenate([[N], N + 1 + np.arange(S)])) vals_.append(np.concatenate([[1.0], p / (1 - beta)])); lb.append([-np.inf]); ub.append([LAMBDA]); ri += 1 rows.append(np.full(N, ri)); cols.append(np.arange(N)); vals_.append(prices); lb.append([-np.inf]); ub.append([W]); ri += 1 for mask in sector_rows: rows.append(np.full(N, ri)); cols.append(np.arange(N)); vals_.append(prices * mask); lb.append([-np.inf]); ub.append([SMAX * W]); ri += 1 Am = sps.csr_matrix((np.concatenate(vals_), (np.concatenate(rows), np.concatenate(cols))), shape=(ri, nv)) vlb = np.zeros(nv); vub = np.full(nv, np.inf); vub[:N] = WMAX * W / prices; vlb[N] = -np.inf res = milp(c, constraints=LinearConstraint(Am, np.concatenate(lb), np.concatenate(ub)), integrality=np.zeros(nv), bounds=Bounds(vlb, vub)) n = res.x[:N]; alpha = res.x[N] return n, alpha def scenarios_sobol(S, seed): z = norm.ppf(qmc.Sobol(d=N, scramble=True, seed=seed).random(S)) p = np.full(S, 1.0 / S) return mu_h + weighted_mm(z, p) @ A.T, p def scenarios_is(S, seed, beta, n1, eps=0.5, shift=None): """Mixture proposal: first eps*S Sobol points shifted by theta along the loss direction of n1, rest nominal.""" v = n1 * prices a = A.T @ v # loss = -v'mu - a'Z d = -a / np.linalg.norm(a) # unit direction of increasing loss in Z-space theta = (shift if shift is not None else norm.ppf(beta)) * d z = norm.ppf(qmc.Sobol(d=N, scramble=True, seed=seed).random(S)) m = int(round(eps * S)) z[:m] += theta # likelihood ratio f/g with g = (1-eps) f + eps f(.-theta): 1 / ((1-eps) + eps exp(theta'z - |theta|^2/2)) lr = 1.0 / ((1 - eps) + eps * np.exp(z @ theta - theta @ theta / 2)) p = lr / lr.sum() ess = 1.0 / (p ** 2).sum() return mu_h + weighted_mm(z, p) @ A.T, p, ess def scenarios_is_multi(S, seed, n1, shifts=(1.6449, 2.3263), eps=(0.25, 0.25)): """Multi-component defensive mixture: nominal + one shifted block per shift, all along the loss direction of n1.""" v = n1 * prices a = A.T @ v d = -a / np.linalg.norm(a) z = norm.ppf(qmc.Sobol(d=N, scramble=True, seed=seed).random(S)) thetas = [sh * d for sh in shifts] start = 0 for th, e in zip(thetas, eps): m = int(round(e * S)); z[start:start + m] += th; start += m dens = (1 - sum(eps)) * np.ones(S) for th, e in zip(thetas, eps): dens = dens + e * np.exp(z @ th - th @ th / 2) lr = 1.0 / dens p = lr / lr.sum() return mu_h + weighted_mm(z, p) @ A.T, p, 1.0 / (p ** 2).sum() def exact_reference(beta, tol=1e-7): k = k_of(beta) cuts = [] ub = WMAX * W / prices c = -(prices * mu_h) base_rows = [prices] + [prices * m for m in sector_rows] base_ub = [W] + [SMAX * W] * len(sector_rows) for it in range(500): Am = np.vstack(base_rows + [g for g, _ in cuts]); hi = np.array(base_ub + [b for _, b in cuts]) res = milp(c, constraints=LinearConstraint(Am, np.full(len(hi), -np.inf), hi), integrality=np.zeros(N), bounds=Bounds(np.zeros(N), ub)) x = res.x; mean, cvar, sd = true_stats(x, beta) if cvar <= LAMBDA * (1 + tol): return x, mean v = x * prices grad = -prices * mu_h + k * prices * (cov_h @ v) / sd cuts.append((grad, LAMBDA - cvar + grad @ x)) raise RuntimeError market_w = md.market_weights(symbols).reindex(symbols).fillna(0.0).values market_n = market_w * W / prices # cap-weighted candidate portfolio (continuous shares) KINDS = sys.argv[5].split(",") if len(sys.argv) > 5 else ["sobol+mm", "sobol+mm+IS(eps=.5)", "IS-market(z_b)", "IS-market(2.33)"] pd.set_option("display.width", 250); pd.set_option("display.float_format", lambda v: f"{v:10.2f}") for beta in BETAS: x_ref, mean_ref = exact_reference(beta) print(f"\n===== beta={beta}: exact optimum E[P&L]={mean_ref:.2f} CHF, CVaR limit {LAMBDA:.0f}") for S in SIZES: results = {} for kind in KINDS: rows = [] t0 = time.time() for seed in range(N_SEEDS): if kind == "sobol+mm": scen, p = scenarios_sobol(S, seed); ess = S n, alpha = solve_lp(scen, p, beta) elif kind == "IS-market(mix)": scen, p, ess = scenarios_is_multi(S, seed + 1000, market_n) n, alpha = solve_lp(scen, p, beta) elif kind.startswith("IS-market"): shift = norm.ppf(beta) if "z_b" in kind else 2.33 scen, p, ess = scenarios_is(S, seed + 1000, beta, market_n, eps=0.5, shift=shift) n, alpha = solve_lp(scen, p, beta) else: eps = 0.5 if ".5" in kind else 0.3 scen1, p1 = scenarios_sobol(S, seed) # stage 1: direction from a plain solve n1, _ = solve_lp(scen1, p1, beta) scen, p, ess = scenarios_is(S, seed + 1000, beta, n1, eps=eps) n, alpha = solve_lp(scen, p, beta) mean_t, cvar_t, _ = true_stats(n, beta) L = -(scen * prices) @ n rows.append({"epnl_true": mean_t, "cvar_true": cvar_t, "w": n * prices / W, "ess": ess, "tail_pts": int((L > alpha).sum())}) df = pd.DataFrame(rows); ws = np.vstack(df["w"].values) pair = [np.abs(ws[i] - ws[j]).sum() for i in range(len(ws)) for j in range(i + 1, len(ws))] results[kind] = { "true E[P&L] mean": df.epnl_true.mean(), "loss vs optimum (CHF)": mean_ref - df.epnl_true.mean(), "true CVaR excess (%)": 100 * (df.cvar_true.mean() / LAMBDA - 1), "true CVaR sd (% of limit)": 100 * df.cvar_true.std() / LAMBDA, "turnover between seeds (%)": 100 * np.mean(pair), "turnover to optimum (%)": 100 * np.mean(np.abs(ws - x_ref * prices / W).sum(axis=1)), "scenarios beyond VaR": df.tail_pts.mean(), "ESS": df.ess.mean(), "seconds": time.time() - t0, } print(f"--- S={S}") print(pd.DataFrame(results).to_string())