
简介这份资源是华泰证券2019年6月发布的金工深度研究报告聚焦遗传规划在选股因子挖掘中的应用面向量化投资研究者、因子开发人员及金融工程方向的学习者。报告系统梳理了遗传规划的原理与完整流程涵盖公式的树形结构表示、适应度计算、选择、交叉、变异与终止条件等核心环节并深入讲解gplearn程序包的定制改进思路包括扩充函数集、引入单因子测试中性化及并行运算加速。资源包内含1个PDF文件大小约3.32MB内容完整呈现从理论到测试的全过程并给出6个具有增量信息且可解释性良好的选股因子案例。目前已有929人学习下载。读者可借此掌握一套“先有公式、后有逻辑”的因子研究方法理解遗传规划突破人脑思维局限、从量价数据中挖掘隐藏因子的优势同时认识其因子复杂、可解释性降低等局限为构建个性化选股框架提供参考。1. 遗传规划选股因子挖掘为什么它比手工挖因子更值得投入做量化选股的人迟早会撞上一堵墙手工构造的因子——动量、反转、波动率、换手率——翻来覆去就那么几十个IC 衰减越来越快同质化严重到令人绝望。你需要的不是再拍脑袋想一个公式而是一套能自动搜索因子表达式的系统。遗传规划Genetic ProgrammingGP就是干这个的它把因子公式当作可进化的个体通过选择、交叉、变异不断迭代从海量算子组合中“进化”出对收益率有预测力的表达式。华泰这份 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 必须开启保护模式protectedTrue否则除零会产生 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, namets_mean, arity2) # 自定义时序算子过去 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, namets_std, arity2) # 注册到 function_set function_set [add, sub, mul, div, log, sqrt, abs, neg, ts_mean, ts_std]上面这段代码做了两件事定义了两个时序算子并注册为 gplearn 可识别的函数对象。arity2表示接受两个参数——第一个是数据序列第二个是窗口长度。注意 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 ICgplearn 默认用make_fitness自定义适应度。选股场景下最常用的适应度是因子值与未来收益的 Rank ICSpearman 相关系数。用 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_betterFalse)这里有个容易翻车的点gplearn 的适应度函数签名是(y, y_pred, sample_weight)顺序不能搞反。另外返回负值是因为 gplearn 默认做最小化你如果想最大化 IC要么取负要么设greater_is_betterTrue。我习惯取负绝对值这样正 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(axis1) | 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( generations20, # 进化代数 population_size2000, # 种群大小 hall_of_fame100, # 名人堂保留的最优个体数 n_components20, # 最终输出的因子数量 function_setfunction_set, metricrank_ic, # 自定义适应度 parsimony_coefficient0.001, # 简约系数控制表达式复杂度 max_samples0.8, # 每代抽样比例 crossover_prob0.7, # 交叉概率 mutation_prob0.1, # 变异概率 p_point_replace0.05, # 变异时替换节点的概率 random_state42, n_jobs-1, # 并行 verbose1 ) 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_samples0.8是 bagging 的思路每代只用 80% 的样本进化剩下的做验证。这在样本量充足时能有效降低过拟合。crossover_prob和mutation_prob之和通常控制在 0.8-0.9剩下的概率是复制。变异概率不宜过高否则退化成随机搜索。3.3 因子有效性检验IC 衰减、换手率和相关性GP 跑出来的因子不能直接用必须过三道检验。第一道是 IC 衰减。计算因子在 T1、T5、T10、T20 的 Rank IC看衰减速度。好的因子应该在 T5 还有一半以上的 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(fFactor {i}: {decay})第二道是换手率。GP 因子容易在相邻截面上剧烈变化导致换手率过高。计算因子值的自相关系数如果 T 日和 T1 日的秩相关低于 0.5说明换手率偏高实盘中交易成本会吃掉大部分收益。第三道是因子间相关性。n_components20输出的 20 个因子如果两两相关都在 0.8 以上那实际上只有 2-3 个独立因子。用pandas.DataFrame.corr()算一下相关矩阵把高相关的因子合并或剔除。4. 避坑与常见问题排查4.1 生成的因子全是 NaN 或 inf现象gp.transform(X)输出的矩阵中大量 NaN 或 inf。原因最常见的是除零和 log 负数。gplearn 的div默认是保护除法但如果你自定义了算子或者用了sqrt、log在负值区域会产生 NaN。另一个原因是输入数据本身有 NaNGP 在进化过程中把 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_size5000、generations50跑了一天还没结束。原因适应度函数计算太慢。如果自定义的时序算子用了 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这种能讲出逻辑的结构。希望帮到你。本文还有配套的精品资源点击获取