简介：这份资源是华泰证券2019年6月发布的金工深度研究报告，聚焦遗传规划在选股因子挖掘中的应用，面向量化投资研究者、因子开发人员及金融工程方向的学习者。报告系统梳理了遗传规划的原理与完整流程，涵盖公式的树形结构表示、适应度计算、选择、交叉、变异与终止条件等核心环节，并深入讲解gplearn程序包的定制改进思路，包括扩充函数集、引入单因子测试中性化及并行运算加速。资源包内含1个PDF文件，大小约3.32MB，内容完整呈现从理论到测试的全过程，并给出6个具有增量信息且可解释性良好的选股因子案例。目前已有929人学习下载。读者可借此掌握一套“先有公式、后有逻辑”的因子研究方法，理解遗传规划突破人脑思维局限、从量价数据中挖掘隐藏因子的优势，同时认识其因子复杂、可解释性降低等局限，为构建个性化选股框架提供参考。
1. 遗传规划选股因子挖掘：为什么它比手工挖因子更值得投入
做量化选股的人迟早会撞上一堵墙：手工构造的因子——动量、反转、波动率、换手率——翻来覆去就那么几十个，IC 衰减越来越快，同质化严重到令人绝望。你需要的不是再拍脑袋想一个公式，而是一套能自动搜索因子表达式的系统。遗传规划（Genetic Programming，GP）就是干这个的：它把因子公式当作可进化的个体，通过选择、交叉、变异不断迭代，从海量算子组合中“进化”出对收益率有预测力的表达式。华泰这份 2019 年的研报把这个思路系统化地落到了 A 股选股场景里，核心工具是 Python 的 gplearn 库。这篇文章不讲虚的，我会把 GP 因子挖掘的完整链路拆开——从算子集设计、适应度函数定义，到 gplearn 的参数配置、因子有效性检验，再到实盘落地时那些让人翻车的细节。适合已经会用 Python 做基础因子回测、想往自动化因子挖掘方向走的从业者。
2. 遗传规划挖因子的底层逻辑与算子集设计
2.1 为什么是遗传规划而不是穷举或随机搜索
因子表达式的搜索空间是组合爆炸的。假设你有 10 个基础变量、8 个算子、表达式最大深度为 5，可能的表达式数量在亿级别以上。穷举不现实，随机搜索效率极低——大部分随机组合出来的公式没有经济学含义，纯粹是噪声拟合。
遗传规划的优势在于它利用了“积累信息”。每一代种群中表现好的个体，其子结构会被保留和重组。一个在早期被发现有预测力的子表达式（比如close/volume这种量价结构），会在后续世代中作为构建块被反复复用和组合。这比从零开始随机试要高效得多。
另一个关键点是 GP 天然支持非线性组合。手工因子往往是线性的或者简单单调的，但市场中的定价偏差经常是非线性的——比如动量因子在小市值股票中的效果和大市值股票中完全不同。GP 通过算子嵌套可以自动捕捉这类交互效应。
常见做法是：先用 gplearn 的SymbolicTransformer做一轮因子生成，把产出的因子当作新特征喂给 LightGBM 或线性模型做最终打分。这样既保留了 GP 的搜索能力，又用树模型兜住了非线性拟合的底。
2.2 算子集和终端集的选型：不是越多越好
算子集决定了 GP 能表达什么。gplearn 内置了基础算术算子（add、sub、mul、div）、超越函数（sin、cos、exp、log）和比较算子。但在选股场景里，不是什么都能用。
我一般会把算子集限制在以下几类：
- 算术算子：add、sub、mul、div。div 必须开启保护模式（
protected=True），否则除零会产生 inf 污染整个种群。 - 一元变换：log、sqrt、abs、neg。log 和 sqrt 要求输入为正，gplearn 会自动做保护处理，但你要清楚这会在负值区域产生截断。
- 时序算子：这是 gplearn 原生不支持的，需要自己扩展。常见的做法是用
ts_mean、ts_std、ts_rank、delay等自定义函数注册进去。
终端集就是基础变量。不要一上来就扔几百个变量进去，维度太高会导致种群收敛慢且容易过拟合。我通常从 10-20 个核心量价变量起步：开盘价、收盘价、最高价、最低价、成交量、成交额、换手率、市值、行业哑变量等。
from gplearn.genetic import SymbolicTransformer from gplearn.functions import make_function import numpy as np # 自定义时序算子：过去 n 日均值 def _ts_mean(x, n): """x: 一维数组，n: 窗口长度""" result = np.full_like(x, np.nan) for i in range(n - 1, len(x)): result[i] = np.nanmean(x[i - n + 1:i + 1]) return result ts_mean = make_function(function=_ts_mean, name='ts_mean', arity=2) # 自定义时序算子：过去 n 日标准差 def _ts_std(x, n): result = np.full_like(x, np.nan) for i in range(n - 1, len(x)): result[i] = np.nanstd(x[i - n + 1:i + 1]) return result ts_std = make_function(function=_ts_std, name='ts_std', arity=2) # 注册到 function_set function_set = ['add', 'sub', 'mul', 'div', 'log', 'sqrt', 'abs', 'neg', ts_mean, ts_std]上面这段代码做了两件事：定义了两个时序算子并注册为 gplearn 可识别的函数对象。arity=2表示接受两个参数——第一个是数据序列，第二个是窗口长度。注意 gplearn 的自定义函数要求输入输出都是 numpy 数组，且要处理好 NaN。
参数说明：ts_mean的第二个参数 n 在 GP 进化过程中会被当作常数处理，但 gplearn 默认的常数是浮点数，你需要做取整处理，否则x[i - 2.7 + 1:i + 1]这种索引会直接报错。我一般会在函数内部加n = int(round(n))并限制 n 的范围在 2 到 60 之间。
2.3 适应度函数：IC 还是 Rank IC
gplearn 默认用make_fitness自定义适应度。选股场景下，最常用的适应度是因子值与未来收益的 Rank IC（Spearman 相关系数）。用 Rank IC 而不是 Pearson IC 的原因是：A 股收益分布厚尾严重，Pearson 相关系数容易被极端值主导。
from gplearn.fitness import make_fitness from scipy.stats import spearmanr def _rank_ic(y_pred, y_true, w): """计算 Rank IC，返回负值因为 gplearn 默认最小化""" ic, _ = spearmanr(y_pred, y_true) if np.isnan(ic): return 0.0 return -abs(ic) # 取负绝对值，最小化负值等价于最大化绝对值 rank_ic = make_fitness(function=_rank_ic, greater_is_better=False)这里有个容易翻车的点：gplearn 的适应度函数签名是(y, y_pred, sample_weight)，顺序不能搞反。另外返回负值是因为 gplearn 默认做最小化，你如果想最大化 IC，要么取负，要么设greater_is_better=True。我习惯取负绝对值，这样正 IC 和负 IC 的因子都能被保留——负 IC 因子取反就是正 IC 因子。
3. 用 gplearn 跑通因子挖掘的完整流程
3.1 数据准备与面板数据处理
GP 挖因子需要的是面板数据：每只股票在每个时间截面上有一组特征值和对应的未来收益。gplearn 本身不处理面板结构，你需要把数据展平或者按截面循环。
我一般会这样做：先把所有股票在所有交易日的数据拼成一个大矩阵，每一行是一个“股票-日期”样本，列是特征。然后按日期分组，在每个截面上分别计算 IC，最后取均值作为适应度。但 gplearn 的适应度函数是全局的，不支持分组计算。折中方案是：用整个面板的 Rank IC 作为适应度，虽然忽略了口截面的差异，但在因子挖掘阶段够用了。
import pandas as pd import numpy as np # 假设 df 是面板数据，包含 stock_id, trade_date, close, volume, future_ret 等列 # 构造特征矩阵 feature_cols = ['close', 'volume', 'turnover', 'market_cap', 'ret_1d', 'ret_5d'] X = df[feature_cols].values y = df['future_ret'].values # 处理缺失值：GP 不能处理 NaN，用行业中位数填充或直接删除 mask = ~(np.isnan(X).any(axis=1) | np.isnan(y)) X, y = X[mask], y[mask] # 标准化：GP 对量纲敏感，不同特征的尺度差异会导致某些算子失效 from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X = scaler.fit_transform(X)数据准备阶段有三个坑：第一，未来收益的计算必须严格对齐，不能用未来数据；第二，标准化要在每个截面内做还是全局做？我建议全局做，因为 GP 进化过程中种群是固定的，截面标准化会导致同一表达式在不同截面的输出不可比；第三，缺失值必须处理干净，gplearn 遇到 NaN 会直接报错或者产出全 NaN 的因子。
3.2 SymbolicTransformer 的关键参数配置
gplearn 提供了两个类：SymbolicRegressor做回归，SymbolicTransformer做特征生成。因子挖掘用后者，因为它输出的是变换后的特征矩阵，可以直接喂给下游模型。
from gplearn.genetic import SymbolicTransformer gp = SymbolicTransformer( generations=20, # 进化代数 population_size=2000, # 种群大小 hall_of_fame=100, # 名人堂保留的最优个体数 n_components=20, # 最终输出的因子数量 function_set=function_set, metric=rank_ic, # 自定义适应度 parsimony_coefficient=0.001, # 简约系数，控制表达式复杂度 max_samples=0.8, # 每代抽样比例 crossover_prob=0.7, # 交叉概率 mutation_prob=0.1, # 变异概率 p_point_replace=0.05, # 变异时替换节点的概率 random_state=42, n_jobs=-1, # 并行 verbose=1 ) gp.fit(X, y) X_new = gp.transform(X) # 生成的新因子矩阵参数逐个说清楚：
generations和population_size是最影响计算时间的。2000 个种群跑 20 代，在 10 万样本量下大概需要 30-60 分钟（取决于表达式深度和算子复杂度）。如果只是验证思路，可以先跑 500 种群、10 代看看效果。
hall_of_fame和n_components的关系：名人堂是所有代中最优个体的集合，n_components是从名人堂中选出的最终因子数。一般hall_of_fame设为n_components的 3-5 倍，给后续去相关留空间。
parsimony_coefficient是控制过拟合的关键。它惩罚表达式长度，值越大越倾向于短表达式。0.001 是个保守的起点，如果发现生成的因子全是深度 5 以上的复杂公式，可以调到 0.005 甚至 0.01。
max_samples=0.8是 bagging 的思路，每代只用 80% 的样本进化，剩下的做验证。这在样本量充足时能有效降低过拟合。
crossover_prob和mutation_prob之和通常控制在 0.8-0.9，剩下的概率是复制。变异概率不宜过高，否则退化成随机搜索。
3.3 因子有效性检验：IC 衰减、换手率和相关性
GP 跑出来的因子不能直接用，必须过三道检验。
第一道是 IC 衰减。计算因子在 T+1、T+5、T+10、T+20 的 Rank IC，看衰减速度。好的因子应该在 T+5 还有一半以上的 IC。
from scipy.stats import spearmanr def ic_decay(factor_values, forward_returns_dict): """forward_returns_dict: {period: returns_array}""" results = {} for period, ret in forward_returns_dict.items(): ic, _ = spearmanr(factor_values, ret) results[period] = ic return results # 假设 X_new 是 GP 生成的因子矩阵，每列一个因子 for i in range(X_new.shape[1]): decay = ic_decay(X_new[:, i], {1: ret_1d, 5: ret_5d, 10: ret_10d, 20: ret_20d}) print(f"Factor {i}: {decay}")第二道是换手率。GP 因子容易在相邻截面上剧烈变化，导致换手率过高。计算因子值的自相关系数，如果 T 日和 T+1 日的秩相关低于 0.5，说明换手率偏高，实盘中交易成本会吃掉大部分收益。
第三道是因子间相关性。n_components=20输出的 20 个因子如果两两相关都在 0.8 以上，那实际上只有 2-3 个独立因子。用pandas.DataFrame.corr()算一下相关矩阵，把高相关的因子合并或剔除。
4. 避坑与常见问题排查
4.1 生成的因子全是 NaN 或 inf
现象：gp.transform(X)输出的矩阵中大量 NaN 或 inf。
原因：最常见的是除零和 log 负数。gplearn 的div默认是保护除法，但如果你自定义了算子或者用了sqrt、log，在负值区域会产生 NaN。另一个原因是输入数据本身有 NaN，GP 在进化过程中把 NaN 传播到了整个表达式。
解决：在fit之前用np.nan_to_num或中位数填充把所有 NaN 处理掉。对于log和sqrt，在自定义函数里加np.abs保护。检查function_set里是否有未保护的除法。
4.2 适应度很高但实盘 IC 为负
现象：训练集上 Rank IC 达到 0.08，实盘跑出来 IC 是 -0.02。
原因：过拟合。GP 的搜索能力太强，在样本内能找到纯粹拟合噪声的表达式。特别是当parsimony_coefficient设得太小、generations设得太大时，这个问题非常严重。
解决：把数据按时间切分，用前 70% 做训练，后 30% 做验证。在适应度函数里加入验证集 IC 的惩罚项。或者用max_samples做 bagging，每代只抽样部分数据。另外，把parsimony_coefficient调大，强制表达式变短。
4.3 种群多样性崩溃，所有个体长得一样
现象：跑了 5 代之后，种群中大部分个体的表达式结构完全相同，适应度不再提升。
原因：选择压力过大，少数高适应度个体迅速占领整个种群。交叉和变异概率太低，无法产生新的结构。
解决：提高mutation_prob到 0.15-0.2，降低选择压力（gplearn 内部用锦标赛选择，可以通过tournament_size参数调整，默认是 20，调小到 5-10 可以增加多样性）。另外，增大population_size也能缓解这个问题。
4.4 计算时间过长，跑一次要几个小时
现象：population_size=5000、generations=50，跑了一天还没结束。
原因：适应度函数计算太慢。如果自定义的时序算子用了 Python 循环，每个个体每代都要遍历整个面板数据，计算量是 O(种群大小 × 代数 × 样本量 × 窗口长度)。
解决：把时序算子用 numpy 的滑动窗口函数重写，避免 Python 循环。或者用 numba 做 JIT 加速。另一个思路是减少样本量——GP 不需要全量数据，随机抽样 30% 的股票做训练就够了。
4.5 因子在大小盘股票上表现差异巨大
现象：因子在全市场 IC 是 0.05，但在沪深 300 成分股内 IC 只有 0.01。
原因：GP 在进化过程中偏向于拟合小市值股票的收益特征，因为小市值股票波动大、噪声多，更容易找到“看起来有效”的表达式。
解决：在适应度函数里做市值中性化。具体做法是：先对因子值和未来收益分别做市值回归，取残差再算 IC。或者在数据预处理阶段就把市值因子从特征中剔除，只用量价数据。
5. 进阶：把 GP 因子接入多因子模型的实战技巧
5.1 因子去相关与正交化
GP 输出的 20 个因子之间往往高度相关。直接全部塞进模型会导致多重共线性，线性模型的系数会变得极不稳定。我一般会做两步处理：先算相关矩阵，把相关系数大于 0.7 的因子配对，保留 IC 更高的那个；然后对剩下的因子做对称正交化（Symmetric Orthogonalization），确保两两正交。
from scipy.linalg import sqrtm import numpy as np def symmetric_orthogonalize(factor_matrix): """对称正交化，保持因子与原始因子的相似性""" cov = np.cov(factor_matrix.T) inv_sqrt = np.linalg.inv(sqrtm(cov)) return factor_matrix @ inv_sqrt X_orth = symmetric_orthogonalize(X_new)对称正交化比施密特正交化的好处是：它不依赖于因子的排列顺序，且正交化后的因子与原始因子的相关性最大。施密特正交化会把第一个因子保留原样，后面的因子被过度改变。
5.2 滚动窗口重新挖掘
市场结构在变，2019 年有效的因子到 2021 年可能就失效了。我习惯每季度重新跑一次 GP，用最近 3 年的数据做训练，生成新的因子池。但要注意：每次重新挖掘后，因子的经济含义可能完全不同，不能简单地把新旧因子混在一起用。
一个实用的做法是：维护一个因子池，每季度新增一批 GP 因子，同时根据最近 6 个月的 IC 表现淘汰表现最差的 20%。这样因子池始终保持更新，又不会因为频繁换血导致策略不稳定。
5.3 用 GP 因子做增强而非替代
GP 因子最大的价值不是替代手工因子，而是提供增量信息。我一般会把 GP 因子和手工因子放在一起，用 LightGBM 做特征重要性排序。如果 GP 因子的重要性排在前 30%，说明它确实带来了手工因子没有捕捉到的信息。如果全部排在末尾，那这轮 GP 挖掘基本是白跑了。
最后说一个我踩过的坑：GP 生成的因子表达式一定要打印出来看。gplearn 的gp._programs属性可以拿到每个因子的表达式字符串。如果看到sin(cos(exp(close)))这种完全没有经济学含义的公式，即使 IC 很高也要警惕——它大概率是过拟合了。好的 GP 因子通常长得像ts_mean(close/volume, 20)或者rank(turnover) * ret_5d这种能讲出逻辑的结构。希望帮到你。
本文还有配套的精品资源，点击获取
