Portfolio_Risk_Management/scripts/experiments/is_noise.py

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