基于遗传规划的行业因子挖掘

量化小白

发布于 2023-04-03 20:25:50

1.9K0

发布于 2023-04-03 20:25:50

文章被收录于专栏：量化小白上分记量化小白上分记

序

之前看了worldquant101,一直对遗传规划挖掘因子的套路比较感兴趣，虽然这样挖出来的因子很容易没有什么逻辑，但想尝试一下看看是怎么回事，也懒得自己折腾，就想用现有的模块做一个试试水。看了很多报告和文献。常用的开源模块如下

最后综合各种资料，决定用gplearn来搞，很早之前发了一篇关于gplearn遗传规划的推文:符号回归和遗传规划，大概学了一下gplearn怎么用，有兴趣可以看一看。

鉴于股票数据很大，自己没有想做的非常精细，就直接用29个中信一级行业指数做了，在行业指数上做因子挖掘，难度小很多，最主要的是数据量小，运行速度很快。全文主要代码、报告、数据获取方式见文末。

遗传规划原理和gplearn

原理这部分主要参考文献[1][2]，讲的比较详细，就直接贴图了，细节见文献。

关于gplearn的说明，细节可以看之前的文章和文献，这里给一张参数说明表

这张表只总结了主要参数，还有几个参数没有提到

比如feature_name是对输入的变量进行命名，如果不指定，最终输出结果的变量名为X0、X1、X2...，max_samples指定做out-of-bag test的样本比例，1表示不做。

parsimony_coefficient用来惩罚过于复杂的公式，这里的惩罚办法也比较朴实：

惩罚后的适应度（fitness) = 惩罚前的适应度（raw fitness) - parsimony_coefficient*公式深度（depth)

parsimony_coefficient的设定方法官方文档建议做CV，也没有深究了。

遗传规划下的行业量价因子挖掘

本文使用中信一级行业指数进行行业因子挖掘，基于gplearn，需要完成的内容包括：

运算符(function set)定义

主要参考下表

自定义运算符部分代码如下，不完全同上文

def _protected_division(x1, x2):
    """Closure of division (x1/x2) for zero denominator."""
    with np.errstate(divide='ignore', invalid='ignore'):
        return np.where(np.abs(x2) > 1e-10 ,np.divide(x1, x2), 1.)

def _protected_sqrt(x1):
    """Closure of square root for negative arguments."""
    return np.sqrt(np.abs(x1))

def _protected_log(x1):
    """Closure of log for zero arguments."""
    with np.errstate(divide='ignore', invalid='ignore'):
        return np.where(np.abs(x1) > 1e-10, np.log(np.abs(x1)), 0.)

def _protected_inverse(x1):
    """Closure of log for zero arguments."""
    with np.errstate(divide='ignore', invalid='ignore'):
        return np.where(np.abs(x1) > 1e-10, 1. / x1, 0.)

def _sigmoid(x1):
    """Special case of logistic function to transform to probabilities."""
    with np.errstate(over='ignore', under='ignore'):
        return 1 / (1 + np.exp(-x1))

def gp_add(x,y):
    
    return x + y

def gp_sub(x,y):
    
    return x - y

def gp_mul(x,y):
    
    return x * y

def gp_div(x,y):
    return _protected_division(x, y)

def gp_sqrt(data):
    
    return _protected_sqrt(data)

def gp_log(data):
    
    return _protected_log(data)

def gp_neg(data):
    
    return np.negative(data)

def gp_inv(data):
    
    return _protected_inverse(data)

def gp_abs(data):
    
    return np.abs(data)

def gp_sin(data):
    
    return np.sin(data)

def gp_cos(data):
    
    return np.cos(data)

def gp_tan(data):
    
    return np.tan(data)

def gp_sig(data):
    
    return _sigmoid(data)



# make_function函数群

gp_add = make_function(function=gp_add, name='gp_add', arity=2)
gp_sub = make_function(function=gp_sub, name='gp_sub', arity=2)
gp_mul = make_function(function=gp_mul, name='gp_mul', arity=2)
gp_div = make_function(function=gp_div, name='gp_div', arity=2)
 
gp_sqrt = make_function(function=gp_sqrt, name='gp_sqrt', arity=1)
gp_log = make_function(function=gp_log, name='gp_log', arity=1)
gp_neg = make_function(function=gp_neg, name='gp_neg', arity=1)
gp_inv = make_function(function=gp_inv, name='gp_inv', arity=1)
gp_abs = make_function(function=gp_abs, name='gp_abs', arity=1)
gp_sin = make_function(function=gp_sin, name='gp_sin', arity=1)
gp_cos = make_function(function=gp_cos, name='gp_cos', arity=1)
gp_tan = make_function(function=gp_tan, name='gp_tan', arity=1)
gp_sig = make_function(function=gp_sig, name='gp_sig', arity=1)


def _rolling_rank(data):
    value = rankdata(data)[-1]
    
    return value


def _rolling_prod(data):
    
    return np.prod(data)

自定义函数在make_functions时会有简单的测试，如果通不过会报错，下图为make_function源代码中测试部分的代码，定义函数时要考虑这一点。

原始因子指定

使用的原始因子为指数的量价数据，即开盘价、收盘价、最高价、最低价、成交量、成交额，日频数据。

适应度定义

尝试了两种适应度定义，IC均值的绝对值、ICIR绝对值。计算上，预测未来5天的收益率，用每一期预测结果计算IC，把所有的IC拼在一起算IC均值和ICIR。

结果分析

样本区间为20051231-20190630，取前2000个交易日作为样本内，即200512-201304，之后的交易日作为样本外。每次取不同的随机数种子，就可以生成多组不同的结果。

参数定义上，generations设定较小，主要是考虑到设置太大，生成的因子会更难解释。其他参数定义比较随意。

est_gp = SymbolicTransformer(population_size= 1000,hall_of_fame = 40,tournament_size = 50,n_components = 20,
                             feature_names = fname,init_method = 'grow',
                       generations= 3, stopping_criteria= 2,function_set = funsets,init_depth = (1,4),
                       const_range = (-1,1),
                       p_crossover=0.7, p_subtree_mutation=0.01,
                       p_hoist_mutation=0.05, p_point_mutation=0.01,
                       p_point_replace = 0.1,
                       max_samples = 1, verbose = 1,
                       parsimony_coefficient=0.01, random_state = 5,
                       metric= icir,
                       n_jobs=3)# 构建一个遗传进化的类

设定好参数后进行优化，过程如下

上图表示每一代群体的平均结果和最有个体的结果，fitness是适应度，实际上是ICIR绝对值，可以看出，随着代数更新，fitness增大，说明有一定作用。

每次fit保留最优的20个因子定义，只给出其中一组的结果如下

其中raw_fitness是ICIR绝对值，fitness是惩罚后的，expression为因子表达式，depth和length分别为公式深度和长度。有个问题是很多公式都重复了，这个gplearn没有办法避免，手动去重吧。

按照第一个因子的定义计算因子后，算因子的累计IC曲线如下：

可以看出，样本内（2014年4月以前），因子IC比较稳定，2016年之前也比较稳定，但是2016年以后，IC非常不稳定，说明过拟合了或者后来因子失效了。尝试了多个种子后发现这个现象是普遍存在的，但也会有少数因子在样本外仍然有一定作用，所以需要大量的实验来寻找好的因子,或者想别的办法避免过拟合。另一方面也确实不容易通过定义找因子的逻辑，这个没什么办法避免。

可以优化的点

个人只是出于兴趣做一个实验，没有打算深究，优化上，可以从以下几个点考虑：