量化之单因子分析的例子
Posted _feivirus_
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了量化之单因子分析的例子相关的知识，希望对你有一定的参考价值。
优矿上的其他人的帖子，以CLV公式为因子，演示了单因子分析的去极值，标准化，中性化，IC检测，回归法检测，分组法检测的整个过程的例子。
原贴地址:
# coding: utf-8
# # CLV公式 #
# close location value,
# ( (close-day_low) - (day_high - close) ) / (day_high - day_low)
# In[ ]:
#导入包
import time
import pandas as pd
import numpy as np
import DataAPI
import statsmodels.api as sm
from scipy import stats as st
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
font=FontProperties()
# In[ ]:
start = time.time()
# 取1个股票1天的行情数据（用前复权的 ）
DataAPI.MktEqudAdjGet(secID=u"",ticker=u"600276",tradeDate=u"",beginDate=u"20210101",endDate=u"20211101",isOpen="",
field=u"secID,tradeDate,openPrice",pandas="1")
# In[ ]:
#取多个股票的某天的因子（智飞生物300122，招商证券600999，招商蛇口001979，青岛中程300208，大参林603233，恒瑞医药600276）截面数据
DataAPI.MktStockFactorsOneDayGet(tradeDate=u"20150227",
secID=u"",
ticker=u"300122,600999,001979,300208,603233,600276",
field=u"secID,tradeDate,openPrice,DAVOL5",pandas="1")
# In[ ]:
#取一只股票的一段时间的所有因子值
DataAPI.MktStockFactorsDateRangeGet(secID=u"",ticker=u"000001",beginDate=u"20170612",endDate=u"20170616",field=u"",pandas="1")
# In[ ]:
#取两个日期间的所有交易日
date_list = DataAPI.TradeCalGet(exchangeCD=u"XSHG,XSHE",beginDate=u"20211010",endDate=u"20211101",
isOpen=u"1",field=u"",pandas="1")
date_list= date_list[date_list['isOpen']==1]['calendarDate'].values.tolist()
#取股票池
universeD = DynamicUniverse('HS300')
# In[ ]:
#定义CLV函数
def CLV(universe, begin, end, file_name):
# 计算收盘价在k线的位置
count = 0
secs_time = 0
start_time = time.time()
N = 50
factor = pd.DataFrame()
# 每支股票
for stk in universe:
# 得到日交易数据
data = DataAPI.MktEqudAdjGet(secID=stk,
beginDate=begin,
endDate=end,
field='secID,tradeDate,closePrice,highestPrice,lowestPrice,openPrice,preClosePrice')
# data.info()
tmp_factor_data = data.sort('tradeDate')
# 计算CLV因子
tmp_factor_data['CLV'] = ( (tmp_factor_data['closePrice'] - tmp_factor_data['lowestPrice']) - (tmp_factor_data['highestPrice'] - tmp_factor_data['closePrice']))/ (tmp_factor_data['highestPrice'] - tmp_factor_data['lowestPrice'])
# 处理出现一字涨跌停
# 参考 https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.loc.html#pandas.DataFrame.loc 文档
tmp_factor_data.loc[(tmp_factor_data['highestPrice']==tmp_factor_data['lowestPrice']) & (tmp_factor_data['openPrice'] > tmp_factor_data['preClosePrice']),'CLV'] = 1
tmp_factor_data.loc[(tmp_factor_data['highestPrice']==tmp_factor_data['lowestPrice']) & (tmp_factor_data['openPrice'] < tmp_factor_data['preClosePrice']),'CLV'] = -1
tmp_factor_data = tmp_factor_data[['tradeDate','CLV']]
tmp_factor_data.columns = ['tradeDate',stk]
if factor.empty:
factor = tmp_factor_data
else:
factor = factor.merge(tmp_factor_data, on='tradeDate', how='outer')
# 打印进度部分
count+=1
if count>0 and count % N ==0:
finish_time = time.time()
print(count,)
print('\\t'+str((np.round(finish_time-start_time)-secs_time,0)) + ' seconds elapsed')
secs_time = (finish_time - start_time)
factor.to_csv(file_name)
return factor
# In[ ]:
#计算CLV
begin_date = '20211101' # 开始日期
end_date = '20211201' # 结束日期
universe = set_universe('HS300') # 股票池
universe = universe[0:5] # 计算速度缓慢，仅以部分股票池作为sample
# 开始日~结束日，50只股票，计算这个因子，这个因子是，每天对应1个股票1个数，比如 2019-1-2，000001.XSHE股票，对应一个因子值，
# 这个值，不是收益率。
CLV(universe=universe, begin= begin_date, end = end_date, file_name='CLV.csv')
# In[ ]:
#读取CLV计算后的csv文件
CLVfactor = pd.read_csv('CLV.csv')
CLVfactor = CLVfactor.set_index(['tradeDate'])
CLVfactorstack = CLVfactor.stack() # 按照tradeDate，变成一个树状结构
CLVfactorstack.hist(figsize=(12,6), bins=50)
# In[ ]:
#去极值 函数定义
def winsorize_series(se):
q = se.quantile([0.025, 0.975])
if isinstance(q, pd.Series) and len(q) == 2:
se[se < q.iloc[0]] = q.iloc[0]
se[se > q.iloc[1]] = q.iloc[1]
return se
#标准化 函数定义
def standardize_series(se):
se_std = se.std()
se_mean = se.mean()
return (se - se_mean)/se_std
# In[ ]:
#对CLV值去极值，标准化
factor_init = CLVfactorstack.groupby(level='tradeDate').apply(winsorize_series)# 去极值
factor_init = factor_init.groupby(level='tradeDate').apply(standardize_series)# 标准化
factor_init.hist(figsize=(12,6), bins=50)
# # 中性化（可选） #
# 中性化的优点：能够将新因子的已知部分剥离掉，剩下的部分为真正新的因子，这样我们能够真正看明白新因子是否真正是新的，还是仅仅原有因子的线性组合加噪音。
# 中性化的缺点：计算复杂，每一个新因子都要进行一次中性化。作为基准的老因子之间，比如行业和市值因子之间也有相关性。
# 中性化计算步骤：
# - 确定将目标因子进行中性化的因子（市值因子，行业因子等等）
# - 将目标因子做为因变量，将市值因子，行业因子做为自变量
# - 运行回归，将回归的残差项做为中性化后的因子
# In[ ]:
# 计算市值因子
def getMarketValueAll(universe, begin, end, file_name):
# 获取股票历史每日市值
print('MarketValue will be calculated for ' + str(len(universe)) + ' stocks:')
count = 0
secs_time = 0
start_time = time.time()
N = 50
ret_data = pd.DataFrame()
for stk in universe:
data = DataAPI.MktEqudAdjGet(secID=stk, beginDate=begin, endDate=end,field='secID,tradeDate,negMarketValue') # 拿取数据
tmp_ret_data = data.sort('tradeDate')
# 市值部分
tmp_ret_data = tmp_ret_data[['tradeDate', 'negMarketValue']]
tmp_ret_data.columns = ['tradeDate', stk]
if ret_data.empty:
ret_data = tmp_ret_data
else:
ret_data = ret_data.merge(tmp_ret_data, on='tradeDate', how='outer')
# 打印进度部分
count += 1
if count > 0 and count % N == 0:
finish_time = time.time()
print(count,)
print(' ' + str(np.round((finish_time - start_time) - secs_time, 0)) + ' seconds elapsed.')
secs_time = (finish_time - start_time)
ret_data.to_csv(file_name)
return ret_data
# In[ ]:
# ----------- 计算股票流通市值 ----------------
print('=======================')
start_time = time.time()
MVfactor = getMarketValueAll(universe=universe, begin=begin_date, end=end_date, file_name='marketvalue.csv')
finish_time = time.time()
# ----------- 去极值和标准化 ----------------
MVfactor = MVfactor.set_index(['tradeDate'])
MVfactorstack = MVfactor.stack()
MVfactor_init = MVfactorstack.groupby(level='tradeDate').apply(winsorize_series) # 去极值
MVfactor_std = MVfactor_init.groupby(level='tradeDate').apply(standardize_series) # 标准化
MVfactor_std.hist(figsize=(12, 6), bins=50)
# 流通市值因子中性化
# 将个股收益率和因子对齐
comnbinedfactor = pd.concat([CLVfactorstack, MVfactor_std], axis=1, join='inner')
comnbinedfactor = comnbinedfactor.reset_index()
comnbinedfactor.columns = ['tradeDate', 'stkID', 'CLV', 'MVfactor_std']
# 按天进行回归，回归残差作为新因子
CLVneutralizedfactor = pd.DataFrame()
unidate = comnbinedfactor.reset_index().tradeDate.drop_duplicates()
unidate = list(unidate)
#遍历每一天，对当天的所有股票一起做回归
for d in unidate:
# 按照每一天去做回归，过滤出的tempdata，是某一天的所有股票的，'stkID'不同
tempdata = comnbinedfactor.loc[comnbinedfactor['tradeDate'] == d, :]
print(">>>>>>>>>>>>>>>>>>>>>")
print("Calculate Date:",d)
print(tempdata.CLV)
print(tempdata.MVfactor_std)
print("<<<<<<<<<<<<<<<<<<<<")
# 以CLV为因变量，市值因子为自变量，做中性化，取残差部分
model = sm.OLS(np.array(tempdata.CLV), np.array(tempdata.MVfactor_std))
results = model.fit()
CLVneutralizedfactor = CLVneutralizedfactor.append(pd.DataFrame(results.resid))
CLVneutralizedfactor['tradeDate'] = comnbinedfactor.tradeDate
CLVneutralizedfactor['stkID'] = comnbinedfactor.stkID
CLVneutralizedfactor = CLVneutralizedfactor.set_index(['tradeDate', 'stkID']) # 联合主键，日期+股票
print("old columns:",CLVneutralizedfactor.columns)
CLVneutralizedfactor.info()
CLVneutralizedfactor.columns = ['CLVneutralizedfactor'] # 中性因子，是每个股票每天1个
# In[ ]:
#定义 计算未来N日收益率函数
#计算当日，到未来5日后的收益率
def getForwardReturns(universe, begin, end, window, file_name):
"""
每天都计算一下从当日，到当日+5日后的收益率，所以最后5天没有值
"""
# 计算个股历史区间前瞻回报率，未来windows天的回报率
print('计算%d只股票的%d天的回报率' % (len(universe), window))
count = 0
secs_time = 0
start_time = time.time()
N = 50
ret_data = pd.DataFrame()
for stock in universe:
data = DataAPI.MktEqudAdjGet(secID=stock,
beginDate=begin,
endDate=end,
field='secID,tradeDate,closePrice')
tmp_ret_data = data.sort('tradeDate')
# 计算历史窗口的前瞻收益率
tmp_ret_data['forwardReturns'] = tmp_ret_data['closePrice'].shift(-window) / tmp_ret_data['closePrice'] - 1.0
tmp_ret_data = tmp_ret_data[['tradeDate','forwardReturns']]
tmp_ret_data.columns = ['tradeDate',stock]
if ret_data.empty:
ret_data = tmp_ret_data
else:
ret_data = ret_data.merge(tmp_ret_data)
ret_data.to_csv(file_name)
print("保存5天因子收益率到文件：",file_name)
return ret_data
# In[ ]:
window_return = 5
forward_5d_return_data = getForwardReturns(
universe=universe,
begin=begin_date,
end=end_date,
window=window_return,
file_name="ForwardReturns_W5_Rtn.csv"
)
forward_5d_return_data = forward_5d_return_data.set_index(['tradeDate'])
forward_5d_return_data.head(3)
# # IC法 #
# 1.确定持仓周期，我们这里是5天
# 2.计算每天股票对应的未来调仓周期的收益
# 3.计算未来调仓周期收益率和因子收益率之间的相关系数
# In[ ]:
# 计算每天的clv收益率和之后5天的股票收益率的秩的相关系数
ic_data = pd.DataFrame(index=CLVfactor.index,columns=['IC','pValue'])
# 计算相关系数
index = 1
for date in ic_data.index:
tmp_illiq = CLVfactor.ix[date] # 得到日期的clv因子暴露，5只股票的
tmp_ret = forward_5d_return_data.ix[date] # 得到日期的5天后的收益率
corr = pd.DataFrame(tmp_illiq)
ret = pd.DataFrame(tmp_ret)
corr.columns = ['corr']
ret.columns = ['ret']
corr['ret'] = ret['ret']
corr = corr[~np.isnan(corr['corr'])][~np.isnan(corr['ret'])]
if len(corr)<5:
continue
# 参考 https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.spearmanr.html#scipy.stats.spearmanr
ic,p_value = st.spearmanr(corr['corr'],corr['ret']) # 计算秩相关系数 Rank_IC
ic_data['IC'][date] = ic
ic_data['pValue'][date] = p_value
index+=1
print("finish Rank_IC/spearman for date: %r, %d/%d, cnt %d" % (date,index,len(ic_data),corr.size))
print("IC mean：%.4f" % ic_data['IC'].mean())
print("IC median：%.4f" % ic_data['IC'].median())
print("IC %d>0, %d<0" % (len(ic_data[ic_data.IC>0]), len(ic_data[ic_data.IC<0])))
# In[ ]:
# 根据每天的CLV,5日收益率,秩相关系数做图
ic_data = ic_data.dropna()
fig = plt.figure(figsize=(16,6))
ax1 = fig.add_subplot(111)
lns1 = ax1.plot(np.array(ic_data.IC),label='IC')
lns = lns1
labs = [l.get_label() for l in lns]
ax1.legend(lns,labs,bbox_to_anchor=[0.5,0.1],loc='best',mode='',borderaxespad=0.,fontsize=12)
ax1.set_xlabel("date",fontproperties=font,fontsize=16)
ax1.set_ylabel("corr",fontproperties=font,fontsize=16)
ax1.set_title("CLV and five days interest rank corr",fontproperties=font,fontsize=16)
ax1.grid()
# In[ ]:
# 计算IC方法二
forward_5d_return_datastack = forward_5d_return_data.stack()
combineMatrix = pd.concat([CLVfactorstack,forward_5d_return_datastack],axis=1,join='inner')
combineMatrix.columns=['CLV','FiveDayfwdRtn']
DayIC = combineMatrix.groupby(level='tradeDate').corr(method='spearman')
DayIC = DayIC.reset_index().drop(['level_1'],axis=1)
DayIC = pd.DataFrame(DayIC.loc[DayIC.CLV!=1,'CLV'])
DayIC.columns=['IC']
print("IC mean：%.4f" % DayIC.mean())
print("IC median：%.4f" % DayIC.median())
print("IC %d > 0, %d < 0" % (len(DayIC[DayIC.IC>0]),len(DayIC[DayIC.IC<0])))
# # 回归法 #
# 1.首先将因子和未来收益率在界面上对齐（日期、代码）
# 2.将未来的收益率作为因变量，因子作为自变量，回归计算出来的系数作为因子收益率
# 3.计算因子收益率的t值等相关统计量
# In[ ]:
combineMatrix = combineMatrix.reset_index()
combineMatrix.columns = ['tradeDate','stockID','CLV','FiveDayfwdRtn']
print(combineMatrix.info())
print(combineMatrix.head(2))
unidate = comnbinedfactor.reset_index().tradeDate.drop_duplicates()
unidate = list(unidate)
print(unidate[:3])
CLVFactorRtn = pd.DataFrame(columns=['CLVfactorRtn','t_values'],index=unidate)
print(CLVFactorRtn.info())
print(CLVFactorRtn.head(1))
# ## 回归法因子检测 ##
# 1.计算t值绝对值的均值，看t值是不是显著不为0，有效性是>2
# 2.t值绝对值大于2的比例-稳定性(比例大于40%)
# 3.计算因子收益率的时间序列上的t值,是不是显著不为0 -- 风险因子？alpha因子？
# In[ ]:
#回归法检测
print(combineMatrix.head(2))
combineMatrix1 = combineMatrix.reset_index()
print(combineMatrix1.head(2))
combineMatrix1.columns = ['tradeDate','stockID','CLV','FiveDayfwdRtn']
# 按天回归，回归系数作为因子收益率
unidate = comnbinedfactor.reset_index().tradeDate.drop_duplicates()
unidate = list(unidate)
CLVFactorRtn = pd.DataFrame(columns=['CLVfactorRtn','t_values'],index=unidate)
for d in unidate:
tempdata = combineMatrix1.loc[combineMatrix1['tradeDate']==d,:]
tempdata = tempdata.dropna()
if len(tempdata)>0:
model = sm.OLS(np.array(tempdata.FiveDayfwdRtn),
np.array(tempdata.CLV))
results = model.fit()
CLVFactorRtn.loc[d,'CLVfactorRtn'] = results.params[0]
CLVFactorRtn.loc[d, 't_values'] = results.tvalues[0]
# 1.计算t值绝对值的均值，看t值是不是显著不为0，--- 有效性
print("t_value abs value mean:" % (CLVFactorRtn.t_values.abs().mean()))
# 2.t值绝对值序列大于2的比例 --- 稳定性
print("positive IC percent: %.2f" % len(CLVFactorRtn[CLVFactorRtn.t_values.abs()>2]/float(len(CLVFactorRtn))))
# 3.计算因子收益率的时间序列上的t值,是不是显著不为0 -- 风险因子？alpha因子？
print("factor interest mean: %.4f" % CLVFactorRtn.CLVfactorRtn.mean())
print("factor interst std: %.4f" % CLVFactorRtn.CLVfactorRtn.std())
print("factor interst sharp: %.4f" % (CLVFactorRtn.CLVfactorRtn.mean().item()/(CLVFactorRtn.CLVfactorRtn.std().item()+0.0000001)))
# In[ ]:
# 回归法画图
fig = plt.figure(figsize=(16,6))
ax1 = fig.add_subplot(111)
lns1 = ax1.plot(np.array(CLVFactorRtn.CLVfactorRtn.cumsum()),label='IC')
lns = lns1
labs = [l.get_label() for l in lns]
ax1.legend(lns,labs,
bbox_to_anchor=[0.5,0.1],
loc='best',
ncol=3,
mode='',
borderaxespad=0.,
fontsize=12)
ax1.set_xlabel("return intesest",fontproperties=font,fontsize=16)
ax1.set_ylabel("date",fontproperties=font,fontsize=16)
ax1.set_title("CLV return accumulation",fontproperties=font,fontsize=16)
ax1.grid()
"""
这个有正、有负，说明其是风险因子，而不是alpha因子，
"""
# # 分层法 #
# In[ ]:
# 分层法检测
n_quantile = 5
# 统计十分位数
cols_mean = [i+1 for i in range(n_quantile)]
cols = cols_mean
excess_returns_means =pd.DataFrame(index=CLVfactor.index,columns=cols)
# 计算ILLIQ分组的超额收益平均值
for date in excess_returns_means.index:
qt_mean_results = []
# 去ILLIQ中的nan
tmp_CLV = CLVfactor.ix[date].dropna()
tmp_return = forward_5d_return_data.ix[date].dropna()
tmp_return_mean = tmp_return.mean()
pct_quantiles = 1.0/n_quantile
for i in range(n_quantile):
down = tmp_CLV.quantile(pct_quantiles*i)
up = tmp_CLV.quantile(pct_quantiles*(i+1))
i_quantile_index = tmp_CLV[(tmp_CLV<=up) & (tmp_CLV>down)].index
mean_tmp = tmp_return[i_quantile_index].mean()- tmp_return_mean
qt_mean_results.append(mean_tmp)
excess_returns_means.ix[date] = qt_mean_results
excess_returns_means.dropna(inplace=True)
excess_returns_means.tail()
print("一共耗时 ： %.2f 秒" %(time.time() - start))
# In[ ]:
#分层法画图
fig = plt.figure(figsize=(16,6))
ax1 = fig.add_subplot(111)
excess_returns_means_dist = excess_returns_means.mean()
excess_dist_plus = excess_returns_means_dist[excess_returns_means_dist>0]
excess_dist_minus = excess_returns_means_dist[excess_returns_means_dist<0]
lns2 = ax1.bar(excess_dist_plus.index,excess_dist_plus.values, align='center',color='g',width=0.1)
lns3 = ax1.bar(excess_dist_minus.index,excess_dist_minus.values, align='center',color='r',width=0.1)
ax1.set_xlim(left=0.5,right=len(excess_returns_means_dist)+0.5)
ax1.set_ylim(-0.008,0.008)
ax1.set_ylabel("return",fontproperties=font,fontsize=16)
ax1.set_xlabel("5 divisions",fontproperties=font,fontsize=16)
ax1.set_xticks(excess_returns_means_dist.index)
ax1.set_xticklabels([int(x) for x in ax1.get_xticks()],fontproperties=font,fontsize=14)
ax1.set_yticklabels([str(x*100)+"0%" for x in ax1.get_yticks()],fontproperties=font,fontsize=14)
ax1.set_title("ILLIQ factor return interest",fontproperties=font,fontsize=16)
ax1.grid()
print("一共耗时 ： %.2f 秒" %(time.time() - start))
# # 策略回测 #
# In[ ]:
# ---- 回测参数，可编辑 ----
# start = '20190101'
# end = '20211201'
# benchmark = 'ZZ500' # 策略参考基准
# universe = set_universe('HS300')# 股票池
start = '20200101'
end = '20211201'
benchmark = 'ZZ500' # 策略参考基准
universe = set_universe('HS300')
capital_base = 100000 # 投资资金
freq = 'd' # 使用日线进行回测
refresh_rate = 5 # 调仓频率, 表示执行handle_data的时间间隔
CLV_dates = CLVfactor.index.values
# 把回测参数封装到SimulationParameters中，供quick_backtest使用
sim_params = quartz.SimulationParameters(start,end,benchmark,universe,capital_base)
# 获取回测行情数据
idxmap,data =quartz.get_backtest_data(sim_params)
# 运行结果
results_illiq =
# 调整参数(选取股票的ILLIQ因子五分位数，进行快速回测
for quantile_five in range(1,6):
# ---- 策略逻辑部分 ----
commission = Commission(0.0002,0.0002)
def initialize(account):
pass
def handle_data(account): # 单个交易日买入卖出
pre_date = account.previous_date.strftime("%Y-%m-%d")
if pre_date not in CLV_dates: # 只在计算过ILLIQ因子的交易日调仓
return
# 拿取调仓日前一个交易日的CLV因子，并按照相应的无分位选择股票
pre_illiq = CLVfactor.ix[pre_date]
pre_illiq = pre_illiq.dropna()
pre_illiq_min = pre_illiq.quantile((quantile_five-1)*0.2)
pre_illiq_max = pre_illiq.quantile(quantile_five*0.2)
my_univ = pre_illiq[pre_illiq>=pre_illiq_min][pre_illiq<pre_illiq_max].index.values
# 调仓逻辑
univ = [x for x in my_univ if x in account.universe]
# 不在股票池，清仓
for stock in account.valid_secpos:
if stock not in univ:
order_to(stock,0)
# 在目标股票池中，等权买入
for stock in univ:
order_pct_to(stock,1.01/len(univ))
# 把回测逻辑封装到 TradeStrategy中，供quick_backtest调用
strategy = quartz.TradingStrategy(initialize,handle_data)
# 回测部分
bt,acct = quartz.quick_backtest(sim_params,strategy,idxmap,data,refresh_rate,commission)
# 对于回测的结果，可以通过 perf_parse 计算风险指标
perf = quartz.perf_parse(bt,acct)
tmp =
tmp['bt'] = bt
tmp['annualized_return'] = perf['annualized_return']
tmp['volatility'] =perf['volatility']
tmp['max_drawdown'] =perf['max_drawdown']
tmp['alpha'] =perf['alpha']
tmp['beta'] =perf['beta']
tmp['sharp'] =perf['sharp']
tmp['information_ratio'] =perf['information_ratio']
results_illiq[quantile_five] = tmp
以上是关于量化之单因子分析的例子的主要内容，如果未能解决你的问题，请参考以下文章