Point72 Cubist 机器学习管线：机器学习交易策略

组件 1：特征工程与 alpha 因子

特征工程把原始的价量数据转化为有预测力的信号。WorldQuant 的研究显示， 他们 101 个公式化 alpha 中有 80 个仍在生产环境运行，两两相关性平均只有 15.9%。这种低相关性让 ensemble 模型能够捕捉多样化的信号。

三类特征：

第 1 类：技术指标（动量、趋势、波动率）

技术指标捕捉的是机器学习模型可以利用的价格形态。研究显示，SVM 的准确率能达到 65-85% ，前提是喂给它经过工程处理的技术特征来预测趋势。

import pandas as pd
import numpy as np
import yfinance as yf
from ta.momentum import RSIIndicator, StochasticOscillator
from ta.trend import MACD, SMAIndicator, EMAIndicator
from ta.volatility import BollingerBands, AverageTrueRange

def calculate_technical_features(df):
    """
    Calculate technical indicators for ML feature engineering.

    Args:
        df: DataFrame with OHLCV data (columns: Open, High, Low, Close, Volume)

    Returns:
        DataFrame with additional technical indicator columns
    """
    close = df['Close']
    high = df['High']
    low = df['Low']
    volume = df['Volume']

    # Momentum Indicators
    rsi = RSIIndicator(close=close, window=14)
    df['rsi'] = rsi.rsi()

    stoch = StochasticOscillator(high=high, low=low, close=close, window=14, smooth_window=3)
    df['stoch_k'] = stoch.stoch()
    df['stoch_d'] = stoch.stoch_signal()

    # Trend Indicators
    macd = MACD(close=close, window_slow=26, window_fast=12, window_sign=9)
    df['macd'] = macd.macd()
    df['macd_signal'] = macd.macd_signal()
    df['macd_diff'] = macd.macd_diff()

    df['sma_20'] = SMAIndicator(close=close, window=20).sma_indicator()
    df['sma_50'] = SMAIndicator(close=close, window=50).sma_indicator()
    df['sma_200'] = SMAIndicator(close=close, window=200).sma_indicator()

    df['ema_12'] = EMAIndicator(close=close, window=12).ema_indicator()
    df['ema_26'] = EMAIndicator(close=close, window=26).ema_indicator()

    # Volatility Indicators
    bb = BollingerBands(close=close, window=20, window_dev=2)
    df['bb_high'] = bb.bollinger_hband()
    df['bb_mid'] = bb.bollinger_mavg()
    df['bb_low'] = bb.bollinger_lband()
    df['bb_width'] = (df['bb_high'] - df['bb_low']) / df['bb_mid']  # Normalized width

    atr = AverageTrueRange(high=high, low=low, close=close, window=14)
    df['atr'] = atr.average_true_range()
    df['atr_pct'] = df['atr'] / close  # ATR as % of price

    # Volume Indicators
    df['volume_sma_20'] = df['Volume'].rolling(window=20).mean()
    df['volume_ratio'] = df['Volume'] / df['volume_sma_20']

    return df

# Example usage
ticker = 'AAPL'
df = yf.download(ticker, start='2020-01-01', end='2025-01-01')
df = calculate_technical_features(df)
print(df[['Close', 'rsi', 'macd', 'bb_width', 'atr_pct']].tail())

这些特征为什么有效： RSI 识别超买超卖状态，MACD 捕捉动量转向，布林带辨别波动率状态，ATR 度量风险。机器学习模型学到的是它们的非线性组合——例如， 高动量（MACD 高于 0）+ 低波动率（bb_width 低于 0.05）= 强烈买入信号.

第 2 类：统计变换（Z 分数、收益率、滞后项）

原始价格是非平稳的（带趋势）。机器学习模型需要平稳的特征。Z 分数和对数收益率能解决这个问题：

def calculate_statistical_features(df):
    """
    Calculate statistical transformations for stationarity.

    Returns z-scores, log returns, and lagged features.
    """
    # Log Returns (stationary)
    df['returns'] = np.log(df['Close'] / df['Close'].shift(1))
    df['returns_5d'] = np.log(df['Close'] / df['Close'].shift(5))
    df['returns_20d'] = np.log(df['Close'] / df['Close'].shift(20))

    # Z-Scores (normalized features)
    for feature in ['Close', 'Volume', 'rsi', 'macd']:
        if feature in df.columns:
            rolling_mean = df[feature].rolling(window=60).mean()
            rolling_std = df[feature].rolling(window=60).std()
            df[f'{feature}_zscore'] = (df[feature] - rolling_mean) / rolling_std

    # Lagged Features (past values as predictors)
    df['close_lag_1'] = df['Close'].shift(1)
    df['close_lag_5'] = df['Close'].shift(5)
    df['volume_lag_1'] = df['Volume'].shift(1)

    # Rolling Statistics
    df['close_std_20'] = df['Close'].rolling(window=20).std()
    df['returns_std_20'] = df['returns'].rolling(window=20).std()  # Realized volatility

    # Skewness & Kurtosis (tail risk indicators)
    df['returns_skew_60'] = df['returns'].rolling(window=60).skew()
    df['returns_kurt_60'] = df['returns'].rolling(window=60).kurt()

    return df

df = calculate_statistical_features(df)
print(df[['returns', 'Close_zscore', 'returns_std_20']].tail())

关键洞见： Z 分数之所以能避免前视偏差，是因为它用的是 滚动窗口 （60 天的均值和标准差），而不是全样本统计量。这确保了特征在实时环境下是算得出来的。

第 3 类：多因子 alpha（Fama-French、WorldQuant 风格）

机构量化用因子模型来捕捉系统性风险溢价。下面是一个简化的实现：

def calculate_alpha_factors(df_dict):
    """
    Calculate cross-sectional alpha factors across multiple stocks.
    Inspired by WorldQuant 101 Alphas and Fama-French factors.

    Args:
        df_dict: Dictionary {ticker: DataFrame with OHLCV + features}

    Returns:
        DataFrame with alpha factors for each stock
    """
    # Combine all stocks into single DataFrame with MultiIndex
    dfs = []
    for ticker, df in df_dict.items():
        df = df.copy()
        df['ticker'] = ticker
        dfs.append(df)

    combined = pd.concat(dfs)
    combined = combined.set_index(['ticker', combined.index])

    # Calculate cross-sectional rankings (key to WorldQuant approach)
    def cross_sectional_rank(group):
        """Rank stocks from 0 (worst) to 1 (best) within each date."""
        return group.rank(pct=True)

    # Momentum Factor (12-month return, skip last month to avoid reversal)
    combined['momentum_12m'] = combined.groupby(level=0)['Close'].pct_change(252)
    combined['momentum_rank'] = combined.groupby(level=1)['momentum_12m'].transform(cross_sectional_rank)

    # Short-Term Reversal (1-month return, negative predictor)
    combined['reversal_1m'] = combined.groupby(level=0)['Close'].pct_change(21)
    combined['reversal_rank'] = combined.groupby(level=1)['reversal_1m'].transform(cross_sectional_rank)

    # Volatility Factor (lower vol = higher rank)
    combined['volatility_60d'] = combined.groupby(level=0)['returns'].transform(lambda x: x.rolling(60).std())
    combined['volatility_rank'] = combined.groupby(level=1)['volatility_60d'].transform(lambda x: 1 - cross_sectional_rank(x))  # Invert: low vol = high rank

    # Volume Factor (abnormal volume)
    combined['volume_20d_avg'] = combined.groupby(level=0)['Volume'].transform(lambda x: x.rolling(20).mean())
    combined['volume_shock'] = combined['Volume'] / combined['volume_20d_avg']
    combined['volume_rank'] = combined.groupby(level=1)['volume_shock'].transform(cross_sectional_rank)

    # Quality Factor (proxied by price stability)
    combined['quality'] = -combined['volatility_60d']  # Simple proxy: stable stocks = high quality
    combined['quality_rank'] = combined.groupby(level=1)['quality'].transform(cross_sectional_rank)

    # Composite Alpha Score (equal-weighted combination)
    alpha_factors = ['momentum_rank', 'reversal_rank', 'volatility_rank', 'volume_rank', 'quality_rank']
    combined['alpha_composite'] = combined[alpha_factors].mean(axis=1)

    return combined

# Example: Download S&P 500 stocks (using top 20 for demo)
tickers = ['AAPL', 'MSFT', 'GOOGL', 'AMZN', 'NVDA', 'META', 'TSLA', 'BRK-B', 'UNH', 'JNJ',
           'V', 'XOM', 'WMT', 'JPM', 'PG', 'MA', 'CVX', 'HD', 'MRK', 'ABBV']

df_dict = {}
for ticker in tickers:
    try:
        df = yf.download(ticker, start='2020-01-01', end='2025-01-01', progress=False)
        df = calculate_technical_features(df)
        df = calculate_statistical_features(df)
        df_dict[ticker] = df
    except:
        print(f"Failed to download {ticker}")

alpha_df = calculate_alpha_factors(df_dict)
print(alpha_df[['momentum_rank', 'volatility_rank', 'alpha_composite']].tail(20))

横截面排序为什么重要： WorldQuant 的 101 Alphas 用的是排序而不是原始数值。这让因子变得 市场中性 （看的是相对表现），并且 对市场状态切换有韧性 （排序在牛市和熊市中都依然成立）。

特征筛选：别把厨房水槽也扔进去

研究提醒我们："把所有能想到的信号一股脑塞进模型，只会稀释预测力。"下面是系统筛选特征的方法：

from sklearn.ensemble import RandomForestRegressor
from sklearn.feature_selection import SelectKBest, f_regression
import matplotlib.pyplot as plt

def select_top_features(X, y, n_features=20, method='random_forest'):
    """
    Select top N features using either RandomForest importance or F-statistic.

    Args:
        X: Feature matrix (DataFrame)
        y: Target variable (returns, rankings, etc.)
        n_features: Number of features to select
        method: 'random_forest' or 'f_statistic'

    Returns:
        List of top feature names
    """
    # Remove NaN values
    valid_idx = ~(X.isna().any(axis=1) | y.isna())
    X_clean = X[valid_idx]
    y_clean = y[valid_idx]

    if method == 'random_forest':
        # Train Random Forest and extract feature importances
        rf = RandomForestRegressor(n_estimators=100, max_depth=10, random_state=42, n_jobs=-1)
        rf.fit(X_clean, y_clean)

        # Get feature importances
        importances = pd.Series(rf.feature_importances_, index=X.columns).sort_values(ascending=False)
        top_features = importances.head(n_features).index.tolist()

        # Plot feature importances
        plt.figure(figsize=(10, 6))
        importances.head(n_features).plot(kind='barh')
        plt.xlabel('Feature Importance')
        plt.title(f'Top {n_features} Features (Random Forest)')
        plt.tight_layout()
        plt.savefig('feature_importance.png')

    elif method == 'f_statistic':
        # Use F-statistic (linear correlation with target)
        selector = SelectKBest(score_func=f_regression, k=n_features)
        selector.fit(X_clean, y_clean)

        # Get selected features
        scores = pd.Series(selector.scores_, index=X.columns).sort_values(ascending=False)
        top_features = scores.head(n_features).index.tolist()

        print(f"Top {n_features} features by F-statistic:")
        print(scores.head(n_features))

    return top_features

# Example: Select top 20 features predicting 5-day forward returns
# Prepare features and target
feature_cols = [col for col in df.columns if col not in ['Open', 'High', 'Low', 'Close', 'Volume', 'ticker']]
X = df[feature_cols]
y = df['Close'].pct_change(5).shift(-5)  # 5-day forward return (target)

top_features = select_top_features(X, y, n_features=20, method='random_forest')
print(f"\nSelected features: {top_features}")

最佳实践： 把特征筛选放在 每一个 walk-forward 折内部 执行，以避免前视偏差。在全样本上选出的特征会泄漏未来信息。

# 错误：用了未来数据
df['close_zscore'] = (df['Close'] - df['Close'].mean()) / df['Close'].std()

# 正确：只用过去的数据
df['close_zscore'] = (df['Close'] - df['Close'].rolling(60).mean()) / df['Close'].rolling(60).std()

# 错误：目标变量用了当天收盘价（下午 4 点才拿得到）
# 但特征用的是上午 9:30 的开盘价 → 若在开盘时交易就产生前视偏差
y = df['Close'].pct_change().shift(-1)  # 明天的回报
X = df[['Open', 'rsi', 'macd']]  # 今天的开盘价（上午 9:30 已知）

# 正确：把时间对齐
y = df['Close'].pct_change().shift(-1)  # 明天的回报
X = df[['Open', 'rsi', 'macd']].shift(1)  # 基于昨天收盘价算出的特征

特征工程检查清单

要点回顾： 特征工程正是散户交易者能够达到 机构效率 70-80%的地方。WorldQuant 的 101 Alphas 是公开的，TA-Lib 是免费的，数据由 yfinance 提供。接下来就是机器学习算法。

组件 2：机器学习模型管线（XGBoost、LightGBM、CatBoost）

特征做好之后，我们来训练梯度提升模型。研究显示，在金融预测任务上，XGBoost、LightGBM 和 Random Forest 持续跑赢 AdaBoost、Bagging 和 ExtraTrees。2025 年一项混合模型研究通过 LightGBM + CatBoost 的 ensemble，取得了 准确率提升 10-15% 的成绩。

算法选择：梯度提升三剑客

训练管线：walk-forward 验证

经典的 k 折交叉验证会在时间序列数据上造成 灾难性的前视偏差 。walk-forward 才是行业黄金标准：

import xgboost as xgb
import lightgbm as lgb
from catboost import CatBoostRegressor
from sklearn.metrics import mean_squared_error, r2_score
import warnings
warnings.filterwarnings('ignore')

class WalkForwardValidator:
    """
    Walk-forward validation for time-series ML models.
    Prevents look-ahead bias by training on past data, testing on future data.
    """

    def __init__(self, train_window_days=504, test_window_days=63, step_days=21):
        """
        Args:
            train_window_days: Training window size (504 days ≈ 2 years)
            test_window_days: Test window size (63 days ≈ 3 months)
            step_days: How many days to shift forward each iteration (21 days ≈ 1 month)
        """
        self.train_window_days = train_window_days
        self.test_window_days = test_window_days
        self.step_days = step_days

    def split(self, df, date_col='Date'):
        """
        Generate train/test splits for walk-forward validation.

        Yields:
            (train_indices, test_indices) for each fold
        """
        df = df.sort_values(date_col).reset_index(drop=True)
        dates = df[date_col].unique()

        # Start after we have enough data for first training window
        start_idx = self.train_window_days

        while start_idx + self.test_window_days < len(dates):
            # Training window: [start - train_window, start)
            train_start = start_idx - self.train_window_days
            train_end = start_idx

            # Test window: [start, start + test_window)
            test_start = start_idx
            test_end = start_idx + self.test_window_days

            # Get date ranges
            train_dates = dates[train_start:train_end]
            test_dates = dates[test_start:test_end]

            # Get corresponding row indices
            train_idx = df[df[date_col].isin(train_dates)].index
            test_idx = df[df[date_col].isin(test_dates)].index

            yield train_idx, test_idx

            # Step forward
            start_idx += self.step_days

    def validate(self, X, y, model_class, model_params, feature_cols):
        """
        Perform walk-forward validation and collect predictions.

        Returns:
            DataFrame with predictions, actuals, and metrics for each fold
        """
        results = []
        fold_num = 0

        for train_idx, test_idx in self.split(X):
            fold_num += 1
            print(f"Fold {fold_num}: Train {len(train_idx)} samples, Test {len(test_idx)} samples")

            # Split data
            X_train = X.loc[train_idx, feature_cols]
            y_train = y.loc[train_idx]
            X_test = X.loc[test_idx, feature_cols]
            y_test = y.loc[test_idx]

            # Remove NaN values
            train_valid = ~(X_train.isna().any(axis=1) | y_train.isna())
            test_valid = ~(X_test.isna().any(axis=1) | y_test.isna())

            X_train_clean = X_train[train_valid]
            y_train_clean = y_train[train_valid]
            X_test_clean = X_test[test_valid]
            y_test_clean = y_test[test_valid]

            if len(X_train_clean) < 100 or len(X_test_clean) < 10:
                print(f"Fold {fold_num}: Insufficient data, skipping")
                continue

            # Train model
            if model_class == 'xgboost':
                model = xgb.XGBRegressor(**model_params)
            elif model_class == 'lightgbm':
                model = lgb.LGBMRegressor(**model_params)
            elif model_class == 'catboost':
                model = CatBoostRegressor(**model_params, verbose=0)

            model.fit(X_train_clean, y_train_clean)

            # Predict
            y_pred = model.predict(X_test_clean)

            # Calculate metrics
            mse = mean_squared_error(y_test_clean, y_pred)
            r2 = r2_score(y_test_clean, y_pred)

            # Store results
            fold_results = pd.DataFrame({
                'fold': fold_num,
                'date': X.loc[test_idx[test_valid], 'Date'].values if 'Date' in X.columns else test_idx[test_valid],
                'actual': y_test_clean.values,
                'predicted': y_pred,
                'mse': mse,
                'r2': r2
            })
            results.append(fold_results)

            print(f"Fold {fold_num}: MSE={mse:.6f}, R²={r2:.4f}")

        return pd.concat(results, ignore_index=True)

# Example usage
# Prepare data (assuming 'alpha_df' from previous feature engineering section)
alpha_df = alpha_df.reset_index()
alpha_df['Date'] = alpha_df['level_1']  # Date is in level_1 from MultiIndex

# Define target: 5-day forward return
alpha_df = alpha_df.sort_values(['ticker', 'Date'])
alpha_df['target'] = alpha_df.groupby('ticker')['Close'].pct_change(5).shift(-5)

# Define features
feature_cols = ['momentum_rank', 'reversal_rank', 'volatility_rank', 'volume_rank',
                'quality_rank', 'rsi_zscore', 'macd', 'bb_width', 'atr_pct',
                'returns_std_20', 'returns_skew_60']

# XGBoost parameters
xgb_params = {
    'n_estimators': 100,
    'max_depth': 5,
    'learning_rate': 0.05,
    'subsample': 0.8,
    'colsample_bytree': 0.8,
    'reg_alpha': 0.1,  # L1 regularization
    'reg_lambda': 1.0,  # L2 regularization
    'random_state': 42,
    'n_jobs': -1
}

# Run walk-forward validation
validator = WalkForwardValidator(train_window_days=504, test_window_days=63, step_days=21)
results = validator.validate(alpha_df, alpha_df['target'], 'xgboost', xgb_params, feature_cols)

print(f"\nOverall Performance:")
print(f"Mean MSE: {results['mse'].mean():.6f}")
print(f"Mean R²: {results['r2'].mean():.4f}")

为什么这样做有效： 每一折都用过去 2 年的数据训练，在接下来 3 个月上测试，然后向前推 1 个月。模型永远看不到未来数据，这模拟的正是你每月重训练一次的真实交易场景。

用 Optuna 做超参数调优（贝叶斯优化）

GridSearch 要搜索 810 种组合，RandomSearch 采样 100 次。Optuna 的贝叶斯方法找到最优 hyperparameter 只需 区区 67 次迭代 （对比基准测试）：

import optuna
from optuna.samplers import TPESampler

def objective_xgboost(trial, X_train, y_train, X_val, y_val, feature_cols):
    """
    Optuna objective function for XGBoost hyperparameter tuning.
    Uses Tree-structured Parzen Estimator (TPE) for Bayesian optimization.
    """
    # Define hyperparameter search space
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 50, 300),
        'max_depth': trial.suggest_int('max_depth', 3, 10),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
        'subsample': trial.suggest_float('subsample', 0.6, 1.0),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
        'reg_alpha': trial.suggest_float('reg_alpha', 0.0, 1.0),
        'reg_lambda': trial.suggest_float('reg_lambda', 0.0, 2.0),
        'min_child_weight': trial.suggest_int('min_child_weight', 1, 10),
        'gamma': trial.suggest_float('gamma', 0.0, 1.0),
        'random_state': 42,
        'n_jobs': -1
    }

    # Train model
    model = xgb.XGBRegressor(**params)

    # Handle NaN values
    train_valid = ~(X_train.isna().any(axis=1) | y_train.isna())
    val_valid = ~(X_val.isna().any(axis=1) | y_val.isna())

    model.fit(X_train[train_valid][feature_cols], y_train[train_valid])

    # Predict on validation set
    y_pred = model.predict(X_val[val_valid][feature_cols])

    # Calculate MSE (minimize)
    mse = mean_squared_error(y_val[val_valid], y_pred)

    return mse

# Run Optuna optimization
def optimize_hyperparameters(X, y, feature_cols, n_trials=50):
    """
    Optimize hyperparameters using Optuna with single train/val split.
    In production, run this within each walk-forward fold.
    """
    # Single train/val split (80/20) for hyperparameter search
    split_idx = int(len(X) * 0.8)
    X_train = X.iloc[:split_idx]
    y_train = y.iloc[:split_idx]
    X_val = X.iloc[split_idx:]
    y_val = y.iloc[split_idx:]

    # Create Optuna study
    study = optuna.create_study(
        direction='minimize',  # Minimize MSE
        sampler=TPESampler(seed=42)
    )

    # Run optimization
    study.optimize(
        lambda trial: objective_xgboost(trial, X_train, y_train, X_val, y_val, feature_cols),
        n_trials=n_trials,
        show_progress_bar=True
    )

    print(f"\nBest hyperparameters:")
    print(study.best_params)
    print(f"Best MSE: {study.best_value:.6f}")

    return study.best_params

# Example: Optimize XGBoost for 50 trials
best_params = optimize_hyperparameters(alpha_df, alpha_df['target'], feature_cols, n_trials=50)

生产环境提示： 运行 Optuna 每季度一次 ，用最近的数据（过去 2 年）。hyperparameter 在月与月之间不会剧烈变化，所以每月重训练可以沿用上一季度的最优参数。

组件 3：模型可解释性（SHAP 值）

"建立对 AI 的信任，是加速数据科学与机器学习在金融服务业落地的关键"——XAI in Finance 系统性综述。SHAP（SHapley Additive exPlanations）提供的正是机构量化用来理解模型 原因 有效的那套框架。

可解释性为什么重要（不只是为了合规）

SHAP vs LIME：为什么金融领域是 SHAP 胜出

为 XGBoost 交易模型接入 SHAP

import shap
import matplotlib.pyplot as plt
import numpy as np

def analyze_shap_values(model, X_train, X_test, feature_names):
    """
    Calculate and visualize SHAP values for XGBoost/LightGBM/CatBoost models.

    Args:
        model: Trained gradient boosting model
        X_train: Training data (for background distribution)
        X_test: Test data (predictions to explain)
        feature_names: List of feature names

    Returns:
        shap_values: SHAP values array
        explainer: SHAP explainer object (reusable)
    """
    # Create SHAP explainer (TreeExplainer for gradient boosting models)
    # Uses TreeSHAP algorithm: polynomial time instead of exponential
    explainer = shap.TreeExplainer(model)

    # Calculate SHAP values for test set
    # This explains each prediction: how much did each feature contribute?
    shap_values = explainer.shap_values(X_test)

    # Expected value: average model output (baseline prediction)
    expected_value = explainer.expected_value
    print(f"Baseline prediction (expected value): {expected_value:.4f}")

    return shap_values, explainer

def plot_shap_summary(shap_values, X_test, feature_names, max_display=20):
    """
    Create summary plot showing global feature importance.

    Each dot is a stock-date prediction. Color = feature value (red high, blue low).
    X-axis = SHAP value (impact on prediction).
    """
    plt.figure(figsize=(10, 8))
    shap.summary_plot(shap_values, X_test, feature_names=feature_names,
                      max_display=max_display, show=False)
    plt.title('Feature Importance (SHAP Values)')
    plt.tight_layout()
    plt.savefig('shap_summary.png', dpi=300)
    print("Saved SHAP summary plot to shap_summary.png")

def plot_shap_waterfall(shap_values, X_test, feature_names, prediction_idx=0):
    """
    Create waterfall plot for a single prediction.
    Shows step-by-step how each feature pushes prediction up or down.

    Args:
        prediction_idx: Which test sample to explain (default: first sample)
    """
    # Create explanation object for single prediction
    shap_exp = shap.Explanation(
        values=shap_values[prediction_idx],
        base_values=explainer.expected_value,
        data=X_test.iloc[prediction_idx],
        feature_names=feature_names
    )

    plt.figure(figsize=(10, 6))
    shap.plots.waterfall(shap_exp, max_display=15, show=False)
    plt.title(f'SHAP Waterfall - Prediction {prediction_idx}')
    plt.tight_layout()
    plt.savefig(f'shap_waterfall_{prediction_idx}.png', dpi=300)
    print(f"Saved SHAP waterfall plot to shap_waterfall_{prediction_idx}.png")

def get_feature_importance_ranking(shap_values, feature_names):
    """
    Calculate global feature importance ranking.

    Returns:
        DataFrame with feature names and importance scores (mean absolute SHAP value)
    """
    # Mean absolute SHAP value = average impact on predictions
    importance = np.abs(shap_values).mean(axis=0)

    # Create DataFrame and sort
    importance_df = pd.DataFrame({
        'feature': feature_names,
        'importance': importance
    }).sort_values('importance', ascending=False)

    return importance_df

def detect_feature_interactions(shap_values, X_test, feature_names,
                                feature_1='momentum_rank', feature_2='volatility_rank'):
    """
    Visualize interaction between two features using SHAP dependence plot.

    Shows how feature_1's effect on prediction changes based on feature_2's value.
    """
    # Find feature indices
    idx_1 = feature_names.index(feature_1)
    idx_2 = feature_names.index(feature_2)

    plt.figure(figsize=(10, 6))
    shap.dependence_plot(
        idx_1, shap_values, X_test, feature_names=feature_names,
        interaction_index=idx_2, show=False
    )
    plt.title(f'Feature Interaction: {feature_1} vs {feature_2}')
    plt.tight_layout()
    plt.savefig(f'shap_interaction_{feature_1}_{feature_2}.png', dpi=300)
    print(f"Saved interaction plot to shap_interaction_{feature_1}_{feature_2}.png")

# Example usage (continuing from walk-forward validation)
# Assume we have trained model and test data from previous section

# Train final model on full training set
X_train_full = alpha_df[feature_cols].iloc[:int(len(alpha_df)*0.8)]
y_train_full = alpha_df['target'].iloc[:int(len(alpha_df)*0.8)]
X_test_full = alpha_df[feature_cols].iloc[int(len(alpha_df)*0.8):]

# Remove NaN
train_valid = ~(X_train_full.isna().any(axis=1) | y_train_full.isna())
test_valid = ~(X_test_full.isna().any(axis=1))

final_model = xgb.XGBRegressor(**xgb_params)
final_model.fit(X_train_full[train_valid], y_train_full[train_valid])

# Calculate SHAP values
shap_values, explainer = analyze_shap_values(
    final_model,
    X_train_full[train_valid],
    X_test_full[test_valid],
    feature_cols
)

# Plot summary (global feature importance)
plot_shap_summary(shap_values, X_test_full[test_valid], feature_cols, max_display=20)

# Get feature importance ranking
importance_df = get_feature_importance_ranking(shap_values, feature_cols)
print("\nTop 10 Features by SHAP Importance:")
print(importance_df.head(10))

# Explain single prediction (e.g., strongest buy signal)
predictions = final_model.predict(X_test_full[test_valid])
strongest_buy_idx = predictions.argmax()
print(f"\nExplaining strongest buy signal (index {strongest_buy_idx}, predicted return: {predictions[strongest_buy_idx]:.4f})")
plot_shap_waterfall(shap_values, X_test_full[test_valid], feature_cols, prediction_idx=strongest_buy_idx)

# Detect feature interactions
detect_feature_interactions(shap_values, X_test_full[test_valid], feature_cols,
                            feature_1='momentum_rank', feature_2='volatility_rank')

解读 SHAP 输出：真实案例

情景： 你的模型预测 AAPL 在接下来 5 天将上涨 3.5%（强烈买入）。SHAP 揭示的内容如下：

Baseline prediction (expected value): 0.0012 (0.12% average return)

SHAP Waterfall for Prediction 147 (AAPL, 2024-08-15):
Feature                     SHAP Value    Cumulative
----------------------------------------
Expected Value                           0.0012
+ momentum_rank (0.92)      +0.0185      0.0197
+ quality_rank (0.88)       +0.0095      0.0292
+ volatility_rank (0.15)    -0.0032      0.0260
+ reversal_rank (0.23)      -0.0018      0.0242
+ volume_rank (0.78)        +0.0048      0.0290
+ rsi_zscore (1.2)          +0.0035      0.0325
+ macd (0.05)               +0.0012      0.0337
+ bb_width (0.02)           +0.0008      0.0345
----------------------------------------
Final Prediction                         0.0345 (3.45%)

解读：

• momentum_rank（0.92）： AAPL 的 12 个月动量排在所有股票前 8%。这为预测回报加上了 +1.85%。 ✓ 有效信号
• quality_rank（0.88）： 低波动、走势稳定的股票，加上 +0.95%。 ✓ 防守型的质量溢价
• volatility_rank（0.15）： 近期波动率飙升（排名靠后），减去 0.32%。 ⚠ 检测到风险因子
• reversal_rank（0.23）： 一个月内的强势上涨意味着均值回归风险，减去 0.18%。 ⚠ 短期涨过头了

行动： 买入信号是成立的（由动量 + 质量驱动），但由于短期反转风险和波动率飙升，仓位应缩减 25%。

生产环境监控：跟踪特征重要性随时间的变化

def monitor_feature_importance_drift(shap_values_history, feature_names, window_size=3):
    """
    Track how feature importance changes across walk-forward folds.
    Detects regime changes when importance rankings shift dramatically.

    Args:
        shap_values_history: List of SHAP value arrays from each fold
        feature_names: List of feature names
        window_size: Number of recent folds to average

    Returns:
        DataFrame with importance trends
    """
    importance_by_fold = []

    for fold_idx, shap_vals in enumerate(shap_values_history):
        importance = np.abs(shap_vals).mean(axis=0)
        importance_by_fold.append({
            'fold': fold_idx,
            **{f: imp for f, imp in zip(feature_names, importance)}
        })

    importance_df = pd.DataFrame(importance_by_fold)

    # Calculate rolling average importance
    for feature in feature_names:
        importance_df[f'{feature}_ma'] = importance_df[feature].rolling(window_size).mean()

    # Detect sudden changes (>50% importance shift)
    importance_df['regime_change'] = False
    for feature in feature_names:
        pct_change = importance_df[feature].pct_change().abs()
        importance_df.loc[pct_change > 0.5, 'regime_change'] = True

    return importance_df

def alert_feature_importance_anomaly(importance_df, feature_names, threshold=0.5):
    """
    Send alert when feature importance shifts dramatically.
    Indicates potential regime change or model degradation.
    """
    latest_fold = importance_df.iloc[-1]

    if latest_fold['regime_change']:
        print("⚠️ REGIME CHANGE DETECTED")
        print(f"Fold {latest_fold['fold']}: Feature importance shifted >50%")

        # Find which features changed
        prev_fold = importance_df.iloc[-2]
        for feature in feature_names:
            pct_change = (latest_fold[feature] - prev_fold[feature]) / prev_fold[feature]
            if abs(pct_change) > threshold:
                direction = "↑" if pct_change > 0 else "↓"
                print(f"  {direction} {feature}: {pct_change:.1%} change")

        print("\nAction: Review model assumptions, consider retraining with different features")

# Example: Track importance across 10 walk-forward folds
# (Assume we've stored SHAP values from each fold in shap_values_history list)
importance_trend = monitor_feature_importance_drift(shap_values_history, feature_cols, window_size=3)
alert_feature_importance_anomaly(importance_trend, feature_cols, threshold=0.5)

要点回顾： SHAP 把黑箱式的机器学习模型变成可解释的系统。对散户交易者来说，这意味着能在回撤期守住纪律（理解模型 原因 有效的原因），也意味着能避免灾难性失败（在伪相关炸掉你的账户之前就发现它）。

组件 4：风险管理与生产部署

拥有成熟 MLOps 管线的机构，把模型失效率降低了 60% 并且部署更新的速度快了 5 倍 ，优于人工监控的做法。本节把机构级的风险管理和 drift 检测搬到散户体量上落地。

模型 drift 的三种类型

用 Evidently AI 实现 drift 检测

from evidently.report import Report
from evidently.metric_preset import DataDriftPreset, RegressionPreset
from evidently import ColumnMapping
import pandas as pd

class DriftMonitor:
    """
    Monitor ML model drift using Evidently AI (free open-source library).
    Tracks data drift, prediction drift, and model performance degradation.
    """

    def __init__(self, reference_data, feature_cols, target_col, prediction_col):
        """
        Args:
            reference_data: Training data (baseline distribution)
            feature_cols: List of feature column names
            target_col: Target variable column name
            prediction_col: Model prediction column name
        """
        self.reference_data = reference_data
        self.feature_cols = feature_cols
        self.target_col = target_col
        self.prediction_col = prediction_col

        # Define column mapping for Evidently
        self.column_mapping = ColumnMapping(
            target=target_col,
            prediction=prediction_col,
            numerical_features=feature_cols
        )

    def check_data_drift(self, current_data, save_report=True):
        """
        Check if current data distribution differs from reference (training) data.

        Returns:
            drift_detected: Boolean
            drift_summary: Dictionary with drift details
        """
        # Create drift report
        data_drift_report = Report(metrics=[DataDriftPreset()])

        data_drift_report.run(
            reference_data=self.reference_data,
            current_data=current_data,
            column_mapping=self.column_mapping
        )

        # Save HTML report
        if save_report:
            data_drift_report.save_html('drift_report.html')
            print("Drift report saved to drift_report.html")

        # Extract drift results
        drift_results = data_drift_report.as_dict()
        drift_detected = drift_results['metrics'][0]['result']['dataset_drift']

        # Get per-feature drift scores
        drift_summary = {
            'dataset_drift': drift_detected,
            'n_drifted_features': drift_results['metrics'][0]['result']['number_of_drifted_columns'],
            'drift_share': drift_results['metrics'][0]['result']['share_of_drifted_columns']
        }

        return drift_detected, drift_summary

    def check_prediction_drift(self, current_data, save_report=True):
        """
        Check if model predictions have drifted from reference period.
        """
        prediction_drift_report = Report(metrics=[RegressionPreset()])

        prediction_drift_report.run(
            reference_data=self.reference_data,
            current_data=current_data,
            column_mapping=self.column_mapping
        )

        if save_report:
            prediction_drift_report.save_html('prediction_drift_report.html')
            print("Prediction drift report saved to prediction_drift_report.html")

        # Extract performance metrics
        results = prediction_drift_report.as_dict()

        # Compare current vs reference performance
        ref_metrics = results['metrics'][0]['result']['reference']
        curr_metrics = results['metrics'][0]['result']['current']

        drift_summary = {
            'ref_mae': ref_metrics['mean_abs_error'],
            'curr_mae': curr_metrics['mean_abs_error'],
            'mae_change': (curr_metrics['mean_abs_error'] - ref_metrics['mean_abs_error']) / ref_metrics['mean_abs_error'],
            'ref_r2': ref_metrics['r2_score'],
            'curr_r2': curr_metrics['r2_score']
        }

        return drift_summary

    def automated_retraining_decision(self, drift_summary, pred_drift_summary,
                                     drift_threshold=0.5, performance_threshold=0.3):
        """
        Decide whether to trigger automated retraining based on drift metrics.

        Args:
            drift_threshold: If >50% of features drift, retrain
            performance_threshold: If MAE increases >30%, retrain

        Returns:
            retrain: Boolean decision
            reason: String explaining decision
        """
        reasons = []

        # Check data drift
        if drift_summary['drift_share'] > drift_threshold:
            reasons.append(f"Data drift: {drift_summary['drift_share']:.1%} of features drifted (threshold: {drift_threshold:.1%})")

        # Check performance degradation
        if pred_drift_summary['mae_change'] > performance_threshold:
            reasons.append(f"Performance degradation: MAE increased {pred_drift_summary['mae_change']:.1%} (threshold: {performance_threshold:.1%})")

        # Decision
        if reasons:
            return True, " | ".join(reasons)
        else:
            return False, "No significant drift detected"

# Example usage (monthly drift monitoring)
# Assume we have reference data (last 2 years training) and current data (last month)

# Prepare reference data (training period)
train_end_idx = int(len(alpha_df) * 0.8)
reference_df = alpha_df.iloc[:train_end_idx].copy()
reference_df['prediction'] = final_model.predict(reference_df[feature_cols].fillna(0))

# Prepare current data (last month of data)
current_df = alpha_df.iloc[train_end_idx:].copy()
current_df['prediction'] = final_model.predict(current_df[feature_cols].fillna(0))

# Initialize drift monitor
monitor = DriftMonitor(
    reference_data=reference_df,
    feature_cols=feature_cols,
    target_col='target',
    prediction_col='prediction'
)

# Check data drift
drift_detected, drift_summary = monitor.check_data_drift(current_df, save_report=True)
print(f"\nData Drift Detected: {drift_detected}")
print(f"Drifted Features: {drift_summary['n_drifted_features']} ({drift_summary['drift_share']:.1%})")

# Check prediction drift
pred_drift_summary = monitor.check_prediction_drift(current_df, save_report=True)
print(f"\nPerformance Change:")
print(f"  Reference MAE: {pred_drift_summary['ref_mae']:.6f}")
print(f"  Current MAE: {pred_drift_summary['curr_mae']:.6f}")
print(f"  Change: {pred_drift_summary['mae_change']:+.1%}")

# Automated retraining decision
retrain, reason = monitor.automated_retraining_decision(drift_summary, pred_drift_summary)
print(f"\nRetrain Model: {retrain}")
print(f"Reason: {reason}")

仓位规模与风险管理

class RiskManager:
    """
    Production risk management for ML trading strategy.
    Implements position sizing, portfolio limits, and automated circuit breakers.
    """

    def __init__(self, initial_capital=50000, max_position_pct=0.02, max_portfolio_risk=0.15):
        """
        Args:
            initial_capital: Starting capital ($)
            max_position_pct: Max risk per trade (2% = $1,000 risk on $50k account)
            max_portfolio_risk: Max portfolio drawdown before circuit breaker (15%)
        """
        self.initial_capital = initial_capital
        self.current_capital = initial_capital
        self.max_position_pct = max_position_pct
        self.max_portfolio_risk = max_portfolio_risk
        self.positions = {}
        self.peak_capital = initial_capital

    def calculate_position_size(self, ticker, predicted_return, atr_pct, confidence=1.0):
        """
        Calculate position size using Kelly Criterion variant.

        Args:
            ticker: Stock ticker
            predicted_return: Model's predicted return (e.g., 0.03 for +3%)
            atr_pct: Average True Range as % of price (volatility measure)
            confidence: Model confidence (0-1), reduce if SHAP shows weak features

        Returns:
            position_size: Dollar amount to invest
        """
        # Base position: 2% risk per trade
        base_risk = self.current_capital * self.max_position_pct

        # Adjust for predicted return magnitude (higher prediction = larger size)
        # But cap at 5% of portfolio to avoid concentration
        predicted_magnitude = min(abs(predicted_return), 0.10)  # Cap at 10% predicted return
        size_multiplier = predicted_magnitude / 0.03  # Normalize to 3% baseline

        # Adjust for volatility (lower vol = larger size)
        volatility_adjustment = 0.02 / max(atr_pct, 0.01)  # 2% baseline ATR

        # Adjust for model confidence (from SHAP analysis)
        confidence_adjustment = confidence

        # Calculate position size
        position_size = base_risk * size_multiplier * volatility_adjustment * confidence_adjustment

        # Cap at 5% of portfolio
        max_position = self.current_capital * 0.05
        position_size = min(position_size, max_position)

        return position_size

    def check_portfolio_limits(self):
        """
        Check if portfolio drawdown exceeds limit.
        Returns:
            circuit_breaker_triggered: Boolean
            current_drawdown: Current DD as decimal (e.g., -0.12 for -12%)
        """
        self.peak_capital = max(self.peak_capital, self.current_capital)
        current_drawdown = (self.current_capital - self.peak_capital) / self.peak_capital

        circuit_breaker = current_drawdown < -self.max_portfolio_risk

        return circuit_breaker, current_drawdown

    def update_capital(self, realized_pnl):
        """
        Update capital after closing position.
        """
        self.current_capital += realized_pnl

        # Check circuit breaker
        circuit_breaker, drawdown = self.check_portfolio_limits()

        if circuit_breaker:
            print(f"⚠️ CIRCUIT BREAKER TRIGGERED")
            print(f"Current Drawdown: {drawdown:.2%} (Limit: {-self.max_portfolio_risk:.2%})")
            print("Action: Close all positions, halt new trades, review strategy")

        return circuit_breaker

# Example usage
risk_mgr = RiskManager(initial_capital=50000, max_position_pct=0.02, max_portfolio_risk=0.15)

# Calculate position size for AAPL prediction
predicted_return = 0.0345  # +3.45% from SHAP example
atr_pct = 0.025  # 2.5% ATR
confidence = 0.9  # High confidence (SHAP showed strong features)

position_size = risk_mgr.calculate_position_size('AAPL', predicted_return, atr_pct, confidence)
print(f"Recommended position size for AAPL: ${position_size:,.0f}")

# Simulate trade outcome
realized_pnl = position_size * predicted_return * 0.6  # Realized 60% of predicted return
circuit_breaker = risk_mgr.update_capital(realized_pnl)

if not circuit_breaker:
    print(f"New capital: ${risk_mgr.current_capital:,.0f}")

生产部署检查清单

要点回顾： 生产环境的机器学习交易要求 时刻保持警觉。市场会变，模型会漂移，昨天管用的东西明天就失灵。drift 检测 + 自动重训练 + SHAP 监控，就是 Point72 长期成功与散户爆仓之间的分水岭。