决策树因子可解释性：用规则集替代黑盒给出合规理由

import numpy as np

rng = np.random.default_rng(20260828)
N = 500

# 4 个因子
momentum = rng.uniform(-0.5, 0.5, N)
volatility = rng.uniform(0.01, 0.05, N)
rsi = rng.uniform(20, 90, N)
vol_change = rng.uniform(0.5, 2.5, N)

X = np.column_stack([momentum, volatility, rsi, vol_change])
feature_names = ["Mom_20d", "Vol_20d", "RSI_14", "VolChg"]

# 真实收益由规则驱动（含交互项）
returns = np.where(
    momentum > 0.15,
    np.where(rsi > 70, -0.02 + rng.normal(0, 0.01, N),  # 高动量+高RSI=超买→做空
             0.005 + rng.normal(0, 0.01, N)),            # 高动量+正常RSI→温和做多
    np.where(volatility > 0.025,
             np.where(vol_change > 1.5, -0.015 + rng.normal(0, 0.01, N),  # 高波动+放量→做空
                      0.01 + rng.normal(0, 0.01, N)),                     # 高波动+缩量→做多
             0.012 + rng.normal(0, 0.01, N))                               # 低波动→做多
)

# 标签：正收益=1（做多），负收益=-1（做空），接近0=0（中性）
labels = np.where(returns > 0.008, 1, np.where(returns < -0.008, -1, 0))

# 划分训练/测试
split = 350
X_train, X_test = X[:split], X[split:]
y_train, y_test = returns[:split], returns[split:]
labels_train, labels_test = labels[:split], labels[split:]

class SimpleCART:
    """从零实现的 CART 决策树（回归模式）"""
    def __init__(self, max_depth=3, min_samples_leaf=30):
        self.max_depth = max_depth
        self.min_samples_leaf = min_samples_leaf
        self.tree = None
        self.feature_names = None

    def _best_split(self, X, y):
        """找最优分裂点：遍历所有特征和切分点，选方差减少最大的。"""
        best_gain, best_feat, best_thresh = -np.inf, None, None
        n = len(y)
        parent_var = np.var(y)
        for feat in range(X.shape[1]):
            thresholds = np.percentile(X[:, feat], [25, 50, 75])  # 3 个候选切分点
            for thresh in thresholds:
                left = X[:, feat] <= thresh
                right = ~left
                if left.sum() < self.min_samples_leaf or right.sum() < self.min_samples_leaf:
                    continue
                # 方差减少量 = 信息增益
                gain = parent_var - (
                    left.sum()/n * np.var(y[left]) + right.sum()/n * np.var(y[right])
                )
                if gain > best_gain:
                    best_gain, best_feat, best_thresh = gain, feat, thresh
        return best_feat, best_thresh, best_gain

    def _build(self, X, y, depth=0):
        """递归构建树"""
        # 终止条件
        if depth >= self.max_depth or len(y) < self.min_samples_leaf * 2:
            return {"leaf": True, "value": y.mean(), "n": len(y), "ic": None}

        feat, thresh, gain = self._best_split(X, y)
        if feat is None or gain < 1e-6:
            return {"leaf": True, "value": y.mean(), "n": len(y), "ic": None}

        left_mask = X[:, feat] <= thresh
        right_mask = ~left_mask

        node = {
            "leaf": False,
            "feat": feat,
            "feat_name": self.feature_names[feat] if self.feature_names else f"feat_{feat}",
            "thresh": thresh,
            "gain": gain,
            "left": self._build(X[left_mask], y[left_mask], depth+1),
            "right": self._build(X[right_mask], y[right_mask], depth+1),
        }
        return node

    def fit(self, X, y, feature_names=None):
        self.feature_names = feature_names
        self.tree = self._build(X, y)

    def predict_one(self, x):
        node = self.tree
        while not node["leaf"]:
            if x[node["feat"]] <= node["thresh"]:
                node = node["left"]
            else:
                node = node["right"]
        return node["value"]

    def predict(self, X):
        return np.array([self.predict_one(x) for x in X])

    def extract_rules(self):
        """提取所有根到叶子的路径作为规则集"""
        rules = []
        def _traverse(node, conditions):
            if node["leaf"]:
                rules.append({
                    "conditions": list(conditions),
                    "prediction": node["value"],
                    "n_samples": node["n"],
                })
                return
            _traverse(node["left"], conditions + [(node["feat_name"], "<=", node["thresh"])])
            _traverse(node["right"], conditions + [(node["feat_name"], ">", node["thresh"])])
        _traverse(self.tree, [])
        return rules

# 训练
tree = SimpleCART(max_depth=3, min_samples_leaf=30)
tree.fit(X_train, y_train, feature_names)

# 预测与 IC
pred_test = tree.predict(X_test)
ic_tree = np.corrcoef(y_test, pred_test)[0, 1]
print(f"决策树 OOS IC = {ic_tree:.3f}")

# 提取规则
rules = tree.extract_rules()
for i, rule in enumerate(rules):
    conds = " & ".join([f"{n} {op} {t:.4f}" for n, op, t in rule["conditions"]])
    print(f"R{i+1}: IF {conds} THEN pred={rule['prediction']:.4f} (n={rule['n_samples']})")

# 对照组1：线性回归
from numpy.linalg import lstsq
beta, *_ = lstsq(X_train, y_train, rcond=None)
pred_linear = X_test @ beta
ic_linear = np.corrcoef(y_test, pred_linear)[0, 1]

# 对照组2：两层 MLP（简化版，纯 numpy）
def relu(x): return np.maximum(0, x)
W1 = rng.normal(0, 0.1, (4, 16))
b1 = rng.normal(0, 0.01, 16)
W2 = rng.normal(0, 0.1, (16, 1))
b2 = rng.normal(0, 0.01, 1)
# 训练（简化版，仅做前向 + 梯度下降 1000 步）
lr = 0.001
for _ in range(1000):
    h = relu(X_train @ W1 + b1)
    out = (h @ W2 + b2).flatten()
    err = out - y_train
    grad_out = err.reshape(-1, 1)
    W2 -= lr * (h.T @ grad_out)
    b2 -= lr * grad_out.mean(axis=0)
    grad_h = (grad_out @ W2.T) * (h > 0)
    W1 -= lr * (X_train.T @ grad_h)
    b1 -= lr * grad_h.mean(axis=0)
pred_mlp = (relu(X_test @ W1 + b1) @ W2 + b2).flatten()
ic_mlp = np.corrcoef(y_test, pred_mlp)[0, 1]

print(f"Linear OLS IC  = {ic_linear:.3f}")
print(f"Decision Tree  = {ic_tree:.3f}")
print(f"MLP (2-layer)  = {ic_mlp:.3f}")

# 逐规则评估：覆盖率 + IC + 边际贡献
rule_metrics = []
for i, rule in enumerate(rules):
    # 计算该规则在测试集上的覆盖率
    mask = np.ones(len(X_test), dtype=bool)
    for feat_name, op, thresh in rule["conditions"]:
        feat_idx = feature_names.index(feat_name)
        if op == "<=":
            mask &= (X_test[:, feat_idx] <= thresh)
        else:
            mask &= (X_test[:, feat_idx] > thresh)

    if mask.sum() > 0:
        r_true = y_test[mask]
        r_pred = np.full(mask.sum(), rule["prediction"])
        if np.std(r_pred) > 0 and np.std(r_true) > 0:
            r_ic = np.corrcoef(r_true, r_pred)[0, 1]
        else:
            r_ic = 0
        coverage = mask.sum() / len(X_test) * 100
        contribution = coverage / 100 * r_ic  # 加权贡献
    else:
        r_ic, coverage, contribution = 0, 0, 0

    rule_metrics.append({
        "rule_id": f"R{i+1}",
        "coverage": coverage,
        "ic": r_ic,
        "contribution": contribution,
    })

# 打印
for rm in rule_metrics:
    print(f"{rm['rule_id']}: coverage={rm['coverage']:.1f}%  IC={rm['ic']:.3f}  contribution={rm['contribution']:.4f}")

approved_rules = {
    "R1": {"conditions": [("Mom_20d", "<=", 0.15), ("Vol_20d", "<=", 0.025)],
            "direction": "long", "max_position": 0.05, "approved_by": "compliance", "date": "2026-08-28"},
    "R5": {"conditions": [("Mom_20d", ">", 0.15), ("RSI_14", ">", 70)],
            "direction": "short", "max_position": 0.03, "approved_by": "compliance", "date": "2026-08-28"},
}

def check_compliance(x, feature_names):
    """每笔交易必须命中已批准规则"""
    for rule_id, rule in approved_rules.items():
        match = True
        for feat_name, op, thresh in rule["conditions"]:
            idx = feature_names.index(feat_name)
            if op == "<=" and x[idx] > thresh:
                match = False
            if op == ">" and x[idx] <= thresh:
                match = False
        if match:
            return rule_id, rule["direction"], rule["max_position"]
    return None, None, 0  # 未命中任何规则 → 不交易