追踪聪明钱 -- 基于TuShare数据的计算

优矿上有一篇文章 "追踪聪明钱 - A股市场交易的微观结构初探", 它根据方正金工的研报『跟踪聪明钱：从分钟行情数据到选股因子』， 基于聪明因子对A股数据进行了自己的回测。具体参加如下链接。

追踪聪明钱 - A股市场交易的微观结构初探

【方正金工专题】跟踪聪明钱：从分钟行情数据到选股因子

优矿上面的实现使用了一些优矿自定义的没有开源的模块和函数，而且优矿上面的很多数据是要收费的，所以这里我用自己的方法来实现文章 "追踪聪明钱 - A股市场交易的微观结构初探" 里面实现的东西。 在此也感谢两篇文章的作者让我从中有所学习。

聪明钱在交易过程中，往往呈现“单笔订单数量更大，订单报价更为激进”的特征，所以可以使用以下指标S来衡量每一分钟交易的“聪明程度”：

St=|Rt|Vt−−√

其中，Rt 为第t分钟的涨跌幅，Vt 为第t分钟的成交量。指标 St 的值越大，则表示该分钟的交易越“聪明”。借助指标S，我们可以通过以下方法筛选聪明钱的交易：

• 对于特定股票、特定时段的所有分钟行情数据，将其按照指标 S 从大到小进行排序；
• 在按 S 排好序的分钟成交量数据中，将成交量累积占比前 20% 视为聪明钱的交易；

下面实例展示了如何从某一天的最后30分钟时间的分钟线计算出所谓20%的聪明钱的交易：

# 拿取一只股票当日的分钟线数据，用以实例展示聪明钱的筛选过程  
  
import numpy as np  
import pandas as pd  
import tushare as ts  
import datetime  
import time  
import tushare as ts  
import os  
import matplotlib  
import matplotlib.pyplot as plt  
  
#用来在matplotlib画图中显示中文，否则中文会被显示成将会是乱码  
#定义自定义字体，文件名是系统中文字体  
myfont = matplotlib.font_manager.FontProperties(fname='C:\\Windows\\Fonts\\msyh.ttc') #自己系统中的中文字体，这里我选的微软雅黑  
#解决负号'-'显示为方块的问题   
matplotlib.rcParams['axes.unicode_minus']=False  

#global variables
data_dir = 'D:\\python_study\\stock_hist_data\\'
strategies_dir = 'D:\\python_study\\strageties\\'  
  
#从文件读取一只股票某一天的1分钟线  
symbol = '603160'  
date=datetime.date(2017,11,6)  
data_dir = 'D:\\python_study\\stock_hist_data\\'  
file_dir =data_dir + symbol+'\\'+str(date.year)+'\\'+str(date.month)  
min_file=file_dir+'\\'+symbol+'_'+str(date)+'_1min_data.h5'  #1分钟线数据文件已经在本地生成，生成方法参加之前的文章  
hdf5_file=pd.HDFStore(min_file, 'r')  
min_df=hdf5_file['data']   
hdf5_file.close()  
  
#在一分钟线基础上计算每分钟的Smart因子  
bar_data = min_df.reset_index()  
bar_data['ticker']=symbol  
bar_data = bar_data[['ticker','time','close','volume']]  
bar_data['return'] = bar_data['close'].pct_change()# 计算每分钟的收益率  
bar_data=bar_data[bar_data['volume']>0]  #删除volume为零的数据，为下一步计算做准备  
bar_data['smartS'] = np.abs(bar_data['return'])/np.sqrt(bar_data['volume'])*100000  # 计算每分钟成交量的聪明度 S 因子值  
bar_data['volume'] = bar_data['volume']/10000.0      # 成交量除以一万，方便画图  
bar_data = bar_data.tail(30)    #选取最后30个数据来做展示  
bar_data['barNo'] = np.arange(1, len(bar_data)+1)      # 给样本分钟线一个序号  
bar_data.index = np.arange(len(bar_data))   
  
#画图1：某一时间段内成交量和 S 因子的分钟线展示  
fig = plt.figure(figsize=(14,10))  
#拥有多个子图时，你会经常看到不同轴域的标签叠在一起, tight_layout可以自动解决这个问题  
#tight_layout会自动调整子图参数，使之填充整个图像区域。  
fig.set_tight_layout(True)                      
ax1 = fig.add_subplot(211)  #添加2行1列subplot中的第1个子图  
ax1.bar(bar_data.index, bar_data.volume, align='center', width=0.4) #子图中添加柱状图,参数1为bar的x坐标轴，参数2为柱的高度,width为柱的宽度
ax1.set_xlabel(u"样本分钟线序号",fontproperties=myfont, fontsize=16)  
ax1.set_ylabel(u"成交量(万)",fontproperties=myfont, fontsize=16) 
ax1.set_xlim(left=-1, right=len(bar_data)) #（ax1上设置x轴上下限） 
ax1.set_xticks(bar_data.index.values) #x轴刻度
ax1.set_xticklabels(bar_data.barNo.values) #x轴刻度标签
ax1.grid()  #显示网格  
ax1.set_title(u"蓝色柱子(左轴)为成交量， 红色圆点(右轴)为 S 因子",fontproperties=myfont, fontsize=16)  


ax2 = ax1.twinx() #Create a new Axes instance with an invisible x-axis and an independent y-axis positioned opposite to the original one  
ax2.plot(bar_data.index, bar_data.smartS, 'o', color='r') #plot x and y using blue circle markers  
ax2.set_ylabel(u"S 因子值 (乘100000)",fontproperties=myfont, fontsize=16)   
ax2.set_ylim(bottom=-0.5*max(bar_data.smartS)) #（ax2上设置y轴上下限）  
 

#画图2： 聪明钱的选择过程   
bar_data = bar_data.sort('smartS', ascending=False) #按smartS的值从大到小排序  
bar_data.index = np.arange(len(bar_data))   
bar_data['accumVolPct'] = bar_data['volume'].cumsum()*1.0/bar_data['volume'].sum() #累积量占总量百分比  
  
ax3 = fig.add_subplot(212)  
bar_data_1 = bar_data[bar_data.accumVolPct<0.2] #smartS大的前20%成交  
bar_data_2 = bar_data[bar_data.accumVolPct>0.2] #剩余的成交  
ax3.bar(bar_data_1.index, bar_data_1.volume, color='r', align='center', width=0.4) #bar_data1和bar_data2显示到同一个subplot中，一个红色一个默认色  
ax3.bar(bar_data_2.index, bar_data_2.volume, align='center', width=0.4)   
ax3.set_xlabel(u"样本分钟线序号", fontproperties=myfont, fontsize=16)  
ax3.set_ylabel(u"成交量(万)", fontproperties=myfont, fontsize=16) 
ax3.set_xlim(left=-1, right=len(bar_data))  
ax3.set_xticks(bar_data.index.values) #x轴刻度
ax3.set_xticklabels(bar_data.barNo.values) #x轴刻度标签
ax3.grid()  
ax3.set_title(u"蓝色柱子(左轴)为成交量， 绿色曲线(右轴)为累积成交量占比", fontproperties=myfont, fontsize=16) 

ax4 = ax3.twinx()  
ax4.plot(bar_data.index, bar_data.accumVolPct, '-o', color='g') #ax2显示, '-o'中的'-'表示有连线
ax4.set_ylabel(u"成交量累积占比", fontproperties=myfont, fontsize=16)  

结果图如下：

如上所示，首先对于30条样本分钟线计算S因子(上图)；其次以S因子由大到小的顺序重新对这些分钟线排序，并按此顺序计算成交量累积占比(下图)，截取S因子最大的前20%成交量所包含的分钟线(下图中的红色柱子)作为聪明钱。

---

如上划分找到聪明钱之后，我们就可以通过这些聪明钱的交易数据来构造聪明钱的情绪因子Q：

Q=VWAPsmartVWAPall

其中，VWAPsmart是聪明钱的成交量加权平均价，VWAPall是所有交易的成交量加权平均价。不难看出，因子Q 实际上反映了在该时间段中聪明钱参与交易的相对价位。之所以将其称为聪明钱的情绪因子，是因为：

• Q越大，表明聪明钱的交易越倾向于出现在价格较高处，这是逢高出货的表现，反映了聪明钱的悲观态度；
• Q越小，则表明聪明钱的交易多出现在价格较低处，这是逢低吸筹的表现，是乐观的情绪。

#函数
#获取给定symbol某一天的历史分钟线
#date格式为datetime.date(2017,11,6)
#eg. getOneMinuteForDate('603160', datetime.date(2017,11,6))
def getOneMinuteForDate(symbol, date):
    global data_dir
    dir=data_dir + symbol + '\\'+ str(date.year) + '\\' + str(date.month)
    minfile=dir+'\\'+symbol+'_'+str(date)[0:10]+'_1min_data.h5'
    data = pd.DataFrame()
    if os.path.exists(minfile):
        hdf5_file=pd.HDFStore(minfile, 'r')
        data=hdf5_file['data'] 
        print "Successfully read 1min file: "+ minfile
        hdf5_file.close()
    return data

#函数
#获取给定symbol和一段时间的历史分钟线数据
def getOneMinuteForDateRange(symbol, dates):
    # 拼凑拿取某只股票的历史分钟线数据
    data = pd.DataFrame()
    for date in dates:
        tmp_data = getOneMinuteForDate(symbol, date)
        data = data.append(tmp_data)
    return data
    
#函数
#返回一段时间内的交易所日历
#date格式为datetime.date(2017,11,6)
def get_exchange_calendar(begin_date, end_date):
    file_name = 'D:\\python_study\\stock_hist_data\\exchange_calendar.h5'
    hdf5_file=pd.HDFStore(file_name, 'r')  
    df=hdf5_file['data']   
    hdf5_file.close()
    df=df.set_index('calendarDate')
    df=df[begin_date:end_date].reset_index() 
    #df['calendarDate'] = pd.to_datetime(df['calendarDate'])
    #df['calendarDate'] = df['calendarDate'].astype(object)
    #df['calendarDate'] = df['calendarDate'].to_pydatetime()
    return df

#函数
#返回沪深300所有股票代码
def get_all_stock_id():
    stock_info=ts.get_hs300s()
    return stock_info['code'].values
    
    
#函数
#计算一只股票在一段时间内的Q因子
#@symbol:股票代码,eg.603160
#@dates: 需要纳入计算的工作日列表
#@q_dates: 需要纳入计算的工作日中的每月最后一个交易日列表
#@window: 计算月度Q因子时从每月最后一个交易日回推天数
#eg.
'''
    cal_dates=get_exchange_calendar(datetime.date(2017,8,1), datetime.date(2017,11,1))
    cal_dates = cal_dates[cal_dates['isOpen']==1]
    all_dates = cal_dates['calendarDate'].values.tolist()                         # 工作日列表
    q_dates = cal_dates[cal_dates['isMonthEnd']==1]['calendarDate'].values.tolist() # 每月最后一个工作日列表
    #getMktSmartMoneyQFactor('603160', all_dates, q_dates, 10)
    getMktSmartMoneyQFactor('600074', all_dates, q_dates, 10)
'''
def getMktSmartMoneyQFactor(symbol, dates, q_dates, window):
    # 计算Q因子
    data = pd.DataFrame(index=q_dates,columns=['Q'])
    
    bar_data = getOneMinuteForDateRange(symbol, dates)    # 拿取该股票dates的所有分钟线数据
    if not len(bar_data)>0:
        return data
    bar_data['accumAdjFactor'] = 1 #目前没有前复权因子，故这里全部假设为1， 以后这里需要改进读正确的取复权因子
    bar_data['adjClosePrice'] = bar_data['close'] * bar_data['accumAdjFactor']
    bar_data['adjVolume'] = bar_data['volume'] / bar_data['accumAdjFactor']
    bar_data['return'] = bar_data['adjClosePrice'].pct_change()*100
    bar_data['S'] = abs(bar_data['return'])/np.sqrt(bar_data['adjVolume'])   # 计算S指标
    bar_data['tradeDate']=bar_data.index.date   #bar_data的index为time类型eg. 2017-11-06 15:00:00 
    
    for i in range(len(q_dates)):                     # 对每个q_dates计算Q因子算是那个月的Q因子
        end = q_dates[i]                          # 此次计算Q因子的日期
        i_begin = dates.index(end) - window + 1    
        begin = dates[i_begin]                      # 此次计算日期往前退windows天的日期
        tmp_data = bar_data[(bar_data['tradeDate']>=begin) & (bar_data['tradeDate']<=end)]
        
        #计算总的成交平均价
        # 排除停牌影响,四分之一以上时间停牌时,不计算Q因子
        if len(tmp_data) == 0 or len(tmp_data[tmp_data['volume']==0]) > window*60:  
            continue
        try:
            vwap_all = tmp_data['amount'].sum()/tmp_data['adjVolume'].sum()      # 总的成交平均价
        except:
            continue
        tmp_data = tmp_data[(tmp_data['S']>0) & (tmp_data['S']<np.inf)] #np.inf为无限大正数
        if len(tmp_data) == 0:   # 排除停牌影响，或者排除有时候tmp_data长度为0的情况 
            continue
        
        #计算聪明钱的成交平均价
        tmp_data.sort_values('S',inplace=True,ascending=False)
        tmp_data['accumVol'] = tmp_data['adjVolume'].cumsum()
        smart_vol = tmp_data['accumVol'].values[-1] * 0.20      # 聪明钱成交量占比为20%
        tmp_data = tmp_data[tmp_data['accumVol'] <= smart_vol]
        try:
            vwap_smart = tmp_data['amount'].sum()/tmp_data['adjVolume'].sum()      # 聪明钱的成交平均价
        except:
            continue
            
        #计算聪明钱的成交情绪因子
        data['Q'][end] = vwap_smart / vwap_all    # 聪明钱的成交情绪因子Q
   
    return data

       
    
# 对所有股票池中的股票计算一段日期内股票历史Q因子，并保存到file中去以供下一步使用
# eg.
'''
symbols=['603160', '603858', '603993', '601997']
begin_date=datetime.date(2017,8,1)
end_date=datetime.date(2017,11,1)
window=10
file_name=strategies_dir + 'smart_money\\q_factor.csv'
getMonthlyQFactorAll(symbols, begin_date, end_date, window, file_name)
'''
def getMonthlyQFactorAll(symbols, begin_date, end_date, window, file_name):
    # 拿取上海证券交易所日历
    cal_dates = get_exchange_calendar(begin_date, end_date) #eg. calc_dates=get_exchange_calendar(datetime.date(2017,8,1), datetime.date(2017,11,1))
    cal_dates = cal_dates[cal_dates['isOpen']==1]
    all_dates = cal_dates['calendarDate'].values.tolist()                          # 工作日列表
    q_dates = cal_dates[cal_dates['isMonthEnd']==1]['calendarDate'].values.tolist() # 每月最后一个工作日列表
    
    print 'Q factor will be calculated for ' + str(len(symbols)) + ' stocks:'
    count = 0
    secs_time = 0
    start_time = time.time()
    
    q_factor = pd.DataFrame()   # 保存计算出来的聪明钱Q因子
    q_factor.to_csv(file_name)
    for stk in symbols:        # 对每一只股票分别计算历史月度Q因子                        
        # 计算聪明钱Q因子时，只需计算每月最后一个交易日的Q因子即可
        tmp_q = getMktSmartMoneyQFactor(stk, all_dates, q_dates, window).reset_index()
        tmp_q.columns = ['tradeDate', stk]
        #把日期格式转换为str格式，这样存储和读取csv之间避免处理date类型
        tmp_q['tradeDate']=tmp_q['tradeDate'].astype('string')
        
        q_factor = pd.read_csv(file_name)
        if q_factor.empty:
            q_factor = tmp_q
        else:
            q_factor = q_factor[q_factor.columns[1:]]
            q_factor = q_factor.merge(tmp_q, on='tradeDate', how='outer')
            
        q_factor = q_factor.sort_values('tradeDate')
        q_factor.to_csv(file_name)
            
        # 打印进度部分
        if count > 0 and count % 10 == 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)
        count += 1
    return q_factor

def geReturnsAll(universe, begin, end, window, file_name):
    # 计算各股票历史区间前瞻回报率
    
    # 拿取上海证券交易所日历
    cal_dates = DataAPI.TradeCalGet(exchangeCD=u"XSHG", beginDate=begin, endDate=end).sort('calendarDate')
    cal_dates = cal_dates[cal_dates['isOpen']==1]
    all_dates = cal_dates['calendarDate'].values.tolist()                          # 工作日列表
    
    print str(window) + ' days forward returns will be calculated for ' + str(len(universe)) + ' stocks:'
    count = 0
    secs_time = 0
    start_time = time.time()
    
    ret_data = pd.DataFrame()   # 保存计算出来的收益率数据
    ret_data.to_csv(file_name)
    for stk in universe:        # 对每一只股票分别计算历史window天前望收益率     
        tmp_ret_data = DataAPI.MktEqudAdjGet(secID=stk, beginDate=begin, endDate=end, 
                                             field='tradeDate,preClosePrice,closePrice')    # 拿取所以工作日的分钟线数据
        # 计算前向收益率
        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', stk]
        
        ret_data = pd.read_csv(file_name)
        if ret_data.empty:
            ret_data = tmp_ret_data
        else:
            ret_data = ret_data[ret_data.columns[1:]]
            ret_data = ret_data.merge(tmp_ret_data, on='tradeDate', how='outer')
            
        ret_data = ret_data.sort('tradeDate')
        ret_data.to_csv(file_name)
            
        # 打印进度部分
        count += 1
        if count > 0 and count % 50 == 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)
    return ret_data
    
 
begin_date=datetime.date(2017,8,1)   # 开始日期，涉及分钟线，计算速度缓慢，这里只针对过去四个月计算作为示例
end_date=datetime.date(2017,11,1)   # 结束日期

symbols=get_all_stock_id()
window_before = 10                # 向前检查smart money的天数
window_return = 20                # 向后检查涨跌幅的天数
smart_money_Q_file_name = strategies_dir + 'smart_money\\q_factor.csv'

# ----------- 计算Q因子部分 ----------------
start_time = time.time()
q_hist_data = getMonthlyQFactorAll(symbols, begin_date,end_date, window_before, smart_money_Q_file_name)
finish_time = time.time()
print ''
print str(finish_time-start_time) + ' seconds elapsed in total.'

# ----------- 计算股票前瞻收益率部分 ----------------
print '======================='
start_time = time.time()
forward_returns_data = geReturnsAll(symbols, begin=begin_date, end=end_date, window=window_return, file_name='SmartMoneyQ_sample_Returns.csv')
finish_time = time.time()
print ''
print str(finish_time-start_time) + ' seconds elapsed in total.'