Portfolio_Risk_Management/scripts/experiments/scenario_noise.py

"""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")