【Python量化】研报复制（一）：《指数高阶矩择时策略》

作者：hzp     Python爱好者社区专栏作者
个人公众号：量化小白上分记

本文为个人对广发研报《指数高阶矩择时策略》的复现，纯代码+结果。

报告认为高阶矩可以刻画资产价格的变化，并且有一定的领先性，可以以此构造指数择时策略，原理见研报（在公众号后台回复“高阶矩”获取研报和代码)

文章为个人对报告的理解，结果并不准确，有问题请指出

python版本：3.6

数据来源：tushare

回测时间段：20050408——20180630

1. 数据及包载入

数据来自tushare

import pandas as pd
import numpy as np
import math
import datetime
import matplotlib.pyplot as plt
import tushare as ts
%matplotlib inline

dateStart = datetime.date(2005,4,8)
dateEnd = datetime.date(2018,6,30)

HS300 = ts.get_k_data('000300', index=True,start = '{}'.format(dateStart),end = '{}'.format(dateEnd))

xticks = np.arange(0,HS300.shape[0],int(HS300.shape[0]/7))

xticklabel = pd.Series(HS300.date[xticks])
plt.figure(figsize=(20,5))
fig = plt.axes()
plt.plot(np.arange(HS300.shape[0]),HS300.close,color = 'red',label = 'close')
fig.set_xticks(xticks)
fig.set_xticklabels(HS300.date[xticks],rotation=60)
plt.show()

def getHighMoment(EarnRate,k,N = 20):
    HighMoment = (EarnRate**k).rolling(window = N).mean()
    for i in range(1,N-1


    
):
        HighMoment[i] = (EarnRate[:i+1]**k).mean()
    return(HighMoment.fillna(0))

def getEMA(moment,alpha,N=120):
    EMA =  pd.DataFrame.ewm(moment,alpha = alpha,adjust = False).mean()
    return(EMA)

报告认为，奇数阶矩对于市场走势具备一定的领先性，即在市场即将出现上升或下降前，有明显数量级变化，我们做出2-7阶矩与指数收盘价的对比图。

其中，日收益率通过每日收盘价计算，高阶矩通过日收益率计算，EMA通过高阶矩计算。

EarnRate = HS300.close.pct_change(1).fillna(0)
Moment_20_2 = getHighMoment(EarnRate,2)
Moment_20_3 = getHighMoment(EarnRate,3)
Moment_20_4 = getHighMoment(EarnRate,4)
Moment_20_5 = getHighMoment(EarnRate,5)
Moment_20_6 = getHighMoment(EarnRate,6)
Moment_20_7 = getHighMoment(EarnRate,7)

五阶矩

HS = plt.subplot(111)
s1 = HS.plot(HS300.close,label = 'HS300')
HS.set_xticks(xticks)
HS.set_xticklabels(HS300.date[xticks],rotation=60)
ax2 = HS.twinx()  
s2 = ax2.plot(Moment_20_5, 'r',label = 'Moment-5')
s = s1+s2
lns = [l.get_label() for l in s]
plt.legend(s,lns)
plt.show()

2至7阶矩，从上往下依次

HS = plt.subplot(611)
HS.plot(HS300.close)
HS.set_xticks(xticks)
HS.set_xticklabels([])
ax2 = HS.twinx()  
ax2.plot(Moment_20_2, 'r')

HS = plt.subplot(612)
HS.plot(HS300.close)
HS.set_xticks(xticks)
HS.set_xticklabels([])
ax2 = HS.twinx()  
ax2.plot(Moment_20_3, 'r')

HS = plt.subplot(613)
HS.plot(HS300.close)
HS.set_xticks(xticks)
HS.set_xticklabels([])
ax2 = HS.twinx()  
ax2.plot(Moment_20_4, 'r')

HS = plt.subplot(614)
HS.plot(HS300.close)
HS.set_xticks(xticks)



    
HS.set_xticklabels([])
ax2 = HS.twinx()  
ax2.plot(Moment_20_5, 'r')

HS = plt.subplot(615)
HS.plot(HS300.close)
HS.set_xticks(xticks)
HS.set_xticklabels([])
ax2 = HS.twinx()  
ax2.plot(Moment_20_6, 'r')


HS = plt.subplot(616)
HS.plot(HS300.close)
HS.set_xticks(xticks)
HS.set_xticklabels(HS300.date[xticks],rotation=60)
ax2 = HS.twinx()  
ax2.plot(Moment_20_7, 'r')
plt.show

<function matplotlib.pyplot.show>

alpha = 0.2,0.4时的EMA

EMA_02 = getEMA(Moment_20_5,0.2)
EMA_04 = getEMA(Moment_20_5,0.4)

HS = plt.subplot(111)
HS.plot(HS300.close)
HS.set_xticks(xticks)
HS.set_xticklabels(HS300.date[xticks],rotation=60


    
)
ax2 = HS.twinx()  
ax2.plot(EMA_02, 'r')
ax3 = ax2.twinx()  
ax3.plot(EMA_04, 'y')
plt.title('Moment-5')
plt.show()

3.回测相关函数定义

切线法寻优函数

输入变量：

• profitrate:日收益率

• Moment：高阶矩

• EMA：EMA4

输出变量：

• num :交易次数

• totalprofit：总收益率

• cumrate：累积收益率

• flag：每天持仓情况 1：多头，-1：空头，0：平仓

def qiexian(profitrate,Moment,EMA):


    

    flag = np.zeros(len(profitrate))
    cumrate = np.ones(len(profitrate))
    for i in [i for i in range(len(profitrate)) if i>1]:
        if EMA[i - 1] > EMA[i-2]:
            flag[i] = 1
        else:
            flag[i] = -1 
    strategy_rate = profitrate * flag + 1
    totalprofit = sum (strategy_rate -1)  

    return totalprofit   

基本高阶矩择时模型中的交易函数

交易逻辑见报告第13页

Sharp：夏普比,年交易天数按250计,MDD：最大回撤率，其他输入输出参数意义同上

def BuySell(profitrate,Moment,EMA):
    flag = np.zeros(len(profitrate))
    cumrate = np.ones(len(profitrate))
    lossflag = 0   # 计算单次择时损失，超过10%平仓
    flagloss = 0  # 单次亏损超过10%时的信号方向

    # 计算最优alpha
    alpha_all = np.arange(0.05,0.5,0.05


    
)

    for i in range(2,len(profitrate)):

        if i%90 ==0 :
            cumrate_all = [qiexian(list(profitrate[i - 90:i]),Moment[i - 90:i].values,
                                   getEMA(Moment,alpha_all[j])[i - 90:i].values) for j in range(len(alpha_all))]            
            # 取使90天内累计收益率最大的alpha作为后续计算的alpha，并计算相应的EMA           
            alpha = alpha_all[np.argmax(cumrate_all)]
            EMA = getEMA(Moment,alpha)    
            EMA = getEMA(Moment,alpha)



        if lossflag -0.1 :
            flag[i] = 0
            flagloss = flag[i-1]
            lossflag = 0
            continue


        if EMA[i - 1] > EMA[i-2] and flagloss != 1:
            flag[i] = 1
            flagloss = 0


        if EMA[i - 1] -2] and flagloss != -1:
            flag[i] = -1
            flagloss = 0

        if flag[i] == flag[i-1] :  
            lossflag = lossflag + profitrate[i]*flag[i]
            lossflag = min(lossflag,0)
        else


    
:
            lossflag = 0


    strategy_rate = profitrate * flag 
    nav = (1+strategy_rate).cumprod()
    cumrate = nav - 1
    totalprofit = nav[len(nav)-1] - 1

    # 交易次数/择时次数
    num = 0
    for i in range(len(flag) - 1):
        if (flag[i+1]!= flag[i]) :
            num+=1
    Sharp = strategy_rate.mean()/strategy_rate.std()* 250**0.5
    MDD = max(1-nav/nav.cummax())
    return(cumrate,num,totalprofit,flag,Sharp,MDD)

拓展后高阶矩交易模型的交易函数

交易逻辑见报告第14页

k：阈值

其他参数含义同上

def BuySell_1(profitrate,Moment,EMA,k):
    flag = np.zeros(len(profitrate))
    cumrate = np.ones(len(profitrate))
    lossflag = 0


    
   # 计算单次择时损失，超过10%平仓
    flagloss = 0  # 单次亏损超过10%时的信号方向

    # 计算最优alpha
    alpha_all = np.arange(0.05,0.5,0.05)
    for i in  range(2,len(profitrate)):

        if i%90 ==0 :
            cumrate_all = [qiexian(list(profitrate[i - 90:i]),Moment[i - 90:i].values,
                                   getEMA(Moment,alpha_all[j])[i - 90:i].values) for j in range(len(alpha_all))]            
            # 取使90天内累计收益率最大的alpha作为后续计算的alpha，并计算相应的EMA           
            alpha = alpha_all[np.argmax(cumrate_all)]
            EMA = getEMA(Moment,alpha)                   
        if lossflag -0.1 :
            flag[i] = 0
            flagloss = flag[i-1]
            lossflag = 0
            continue


        if EMA[i - 1] > EMA[i-2]*(1+k) and flagloss != 1:
            flag[i] = 1
            flagloss = 0


        if EMA[i - 1] -2]*(1-k) and flagloss != -1:
            flag[i] = -1
            flagloss = 0


    


        if EMA[i - 1] > EMA[i-2]*(1-k) and EMA[i - 1] -2]*(1+k) :
            flag[i] = 0
            flagloss = 0

        if flag[i] == flag[i-1] :  
            lossflag = lossflag + profitrate[i]*flag[i]
            lossflag = min(lossflag,0)
        else:
            lossflag = 0


    strategy_rate = profitrate * flag 
    nav = (1+strategy_rate).cumprod()
    cumrate = nav - 1
    totalprofit = nav[len(nav)-1] - 1

    # 交易次数/择时次数
    num = 0
    for i in range(len(flag) - 1):
        if (flag[i+1]!= flag[i]) :
            num+=1
    Sharp = strategy_rate.mean()/strategy_rate.std()* 250**0.5
    MDD = max(1-nav/nav.cummax())
    return(cumrate,num,totalprofit,flag,Sharp,MDD)

4.回测并展示结果

• 回测1：基本的高阶矩择时模型，EMA初始alpha=0.4

• 回测2：阈值k=1%,其他同上

• 回测3：阈值k=2%,其他同上

cumrate,num,totalprofit,flag,Sharp,MDD = BuySell(EarnRate,Moment_20_5,EMA_04)
cumrate_01,num_01,totalprofit_01,flag_02,Sharp_01,MDD_01 = BuySell_1(EarnRate,Moment_20_5,EMA_04,0.01)
cumrate_02,num_02,totalprofit_02,flag_02,Sharp_02,MDD_02 = BuySell_1(EarnRate,Moment_20_5,EMA_04,0.02)

print('回测1交易次数为:%s' % num)
print('回测2交易次数为:%s' % num_01)
print('回测3交易次数为:%s' % num_02)
print('')
print('回测1总收益率为:%s%%' % (totalprofit*100))
print('回测2总收益率为:%s%%' % (totalprofit_01*100))
print('回测3总收益率为:%s%%' % (totalprofit_02*100))
print('')
print('回测1夏普比为:%s' % Sharp)
print('回测2夏普比为:%s' % Sharp_01)
print('回测3夏普比为:%s' % Sharp_02)
print('')
print('回测1最大回撤为:%s%%' % (MDD*100))
print('回测2最大回撤为:%s%%' % (MDD_01*100))
print('回测3最大回撤为:%s%%' % (MDD_02*100))

    
回测1交易次数为:298
回测2交易次数为:353
回测3交易次数为:402

回测1总收益率为:870.434213293%
回测2总收益率为:859.177156583%
回测3总收益率为:571.69369844%

回测1夏普比为:0.8071530425449338
回测2夏普比为:0.8120229809371885
回测3夏普比为:0.7049768658906784

回测1最大回撤为:40.25782265430625%
回测2最大回撤为:34.72891044547697%
回测3最大回撤为:39.44870302175746%

以HS300为基准，回测1累计收益率展示

HS = plt.subplot(111)
HS.plot(HS300.close)
HS.set_xticks(xticks)
HS.set_xticklabels(HS300.date[xticks],rotation=60)
ax2 = HS.twinx()  
ax2.plot(cumrate, 'r')

不同阈值下的累计净值曲线对比

HS = plt.subplot(111)
HS.plot(cumrate,


    
'r',cumrate_01,'g',cumrate_02,'k')
HS.set_xticks(xticks)
HS.set_xticklabels(HS300.date[xticks],rotation=60)
plt.legend(['0%','1%','2%'])
plt.show()

总体来看，策略有明显超额收益的，但与文章结果出入较大，说明代码中可能有一些与报告原意不符的地方，此外，报告只回测到2015年，而从上图可以看出，策略2016年之后出现了较大回撤。

参考文献

广发证券-指数高阶矩择时策略

Python爱好者社区历史文章大合集：

Python爱好者社区历史文章列表（每周append更新一次）

福利：文末扫码立刻关注公众号，“Python爱好者社区”，开始学习Python课程：

关注后在公众号内回复“课程”即可获取：

小编的Python入门免费视频课程！！！

【最新免费微课】小编的Python快速上手matplotlib可视化库！！！

崔老师爬虫实战案例免费学习视频。

陈老师数据分析报告制作免费学习视频。

玩转大数据分析！Spark2.X+Python 精华实战课程免费学习视频。