风险模型与因子归因

r = X·f + u

r：N 只股票的收益向量
X：N×K 的因子暴露矩阵
f：K 维因子收益
u：特质收益（假设互不相关）

Σ = X·Σf·Xᵀ + Δ

Σf：K×K 的因子协方差（K 通常 30~60）
Δ：N×N 的对角矩阵（特质风险）

def standardize_exposure(raw, mktcap):
    # 市值加权均值与标准差
    w = mktcap / mktcap.sum()
    mean = (raw * w).sum()
    std = np.sqrt((w * (raw - mean) ** 2).sum())
    return (raw - mean) / std          # 市值加权的 z-score

import statsmodels.api as sm

def estimate_factor_returns(returns, X, weights=None):
    # 每期横截面回归 r = X·f + u，加权最小二乘
    Xw = sm.add_constant(X)                  # 加截距（可吸收市场均值）
    model = sm.WLS(returns, Xw, weights=weights).fit()
    return model.params, model.resid         # 因子收益与特质收益

def volatility_regime_adjust(cov, current_vol, long_run_vol, beta=0.5):
    # 把当前协方差按波动率偏离长期均值的程度做修正
    scale = (long_run_vol / current_vol) ** beta
    return cov * scale

def shrink_specific_var(sample_var, prior_var, n_obs, prior_weight=50):
    # 向先验收缩，样本量越少收缩越强
    lam = prior_weight / (prior_weight + n_obs)
    return lam * prior_var + (1 - lam) * sample_var

先验可以按「行业 × 市值组」计算：
小盘股的特质风险天然高于大盘股，先验应该分组而非全局

def decompose_risk(w, X, Sigma_f, Delta):
    total_var = w @ (X @ Sigma_f @ X.T + Delta) @ w
    factor_var = w @ X @ Sigma_f @ X.T @ w       # 因子贡献
    specific_var = w @ Delta @ w                 # 特质贡献
    return {
        'total_vol': np.sqrt(total_var) * np.sqrt(252),
        'factor_vol': np.sqrt(factor_var) * np.sqrt(252),
        'specific_vol': np.sqrt(specific_var) * np.sqrt(252),
        'factor_ratio': factor_var / total_var,   # 因子风险占比
    }

def factor_risk_contribution(w, X, Sigma_f):
    expo = X.T @ w                       # 组合的因子暴露（K 维）
    factor_var = expo @ Sigma_f @ expo
    contrib = expo * (Sigma_f @ expo) / factor_var   # 各因子的风险贡献，和为 1
    return dict(zip(FACTOR_NAMES, contrib))

主动收益 = 配置效应 + 选股效应 + 交互效应

配置效应：行业权重的偏离 × 行业基准收益
选股效应：行业内的选股 × 行业权重
交互效应：配置与选股的交叉项

def brinson_attribution(w_p, w_b, r_p, r_b, sectors):
    alloc = ((w_p - w_b) * r_b).sum()        # 配置效应
    select = (w_b * (r_p - r_b)).sum()       # 选股效应
    interact = ((w_p - w_b) * (r_p - r_b)).sum()   # 交互效应
    return {'allocation': alloc, 'selection': select, 'interaction': interact}

def factor_attribution(portfolio_returns, exposures, factor_returns):
    explained = (exposures * factor_returns).sum(axis=1)   # 因子解释的收益
    alpha = portfolio_returns - explained                  # 特质 alpha
    return {
        'factor_return': explained.sum(),
        'alpha': alpha.sum(),
        'alpha_ir': alpha.mean() / alpha.std() * np.sqrt(252),   # alpha 的信息比率
    }

MRC_i = ∂σ_p / ∂w_i = (Σw)_i / σ_p

风险贡献（RC）= w_i · MRC_i
组合波动 = Σ RC_i     （欧拉分解，贡献之和等于总风险）

def risk_contributions(w, cov):
    port_var = w @ cov @ w
    port_vol = np.sqrt(port_var)
    mrc = (cov @ w) / port_vol          # 边际贡献
    rc = w * mrc                        # 总贡献
    return rc, rc / port_vol            # 绝对贡献与占比

from scipy.optimize import minimize

def risk_budget(cov, budget):
    # budget: 各资产的目标风险贡献占比，和为 1
    n = len(budget)
    def objective(w):
        rc = w * (cov @ w)
        rc_frac = rc / rc.sum()
        return ((rc_frac - budget) ** 2).sum()
    cons = [{'type': 'eq', 'fun': lambda w: w.sum() - 1}]
    res = minimize(objective, np.ones(n) / n, bounds=[(1e-6, 1)] * n, constraints=cons)
    return res.x

def scenario_loss(w, X, factor_shock, specific_shock=0.0):
    # factor_shock: K 维的因子冲击（如「小盘 -10%、成长 -8%」）
    factor_pnl = (X @ factor_shock) @ w
    specific_pnl = specific_shock * (w ** 2).sum() ** 0.5
    return factor_pnl + specific_pnl

def reverse_stress(w, X, max_loss):
    # 找到让组合损失达到 max_loss 的最小因子冲击方向
    exposure = X.T @ w                       # 组合的因子暴露
    # 最脆弱的方向 = 与暴露方向相反的冲击
    direction = -exposure / np.linalg.norm(exposure)
    magnitude = max_loss / np.linalg.norm(exposure)
    return direction, magnitude

import cvxpy as cp

def optimize_with_risk_model(mu, w, X, Sigma_f, Delta, max_factor_expo=0.05):
    n, k = X.shape
    cov = X @ Sigma_f @ X.T + Delta          # 由风险模型构造协方差
    objective = cp.Maximize(mu @ w - 5.0 / 2 * cp.quad_form(w, cov))
    constraints = [
        cp.sum(w) == 1,
        w <= 0.03,
        cp.abs(X.T @ w) <= max_factor_expo,  # 限制因子暴露
    ]
    prob = cp.Problem(objective, constraints)
    prob.solve(solver=cp.OSQP)
    return w.value