[关闭]
@1007477689 2020-07-30T09:42:06.000000Z 字数 8193 阅读 506

研报复制（一）：《指数高阶矩择时策略》

量化

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

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

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

python 版本：3.6
数据来源：tushare
回测时间段：20050408——20180630

1. 数据及包载入

数据来自 tushare

i. 导入必须的 Python 库

  1. import pandas as pd
  2. import numpy as np
  3. import math
  4. import datetime
  5. import matplotlib.pyplot as plt
  6. import tushare as ts
  7. %matplotlib inline

ii. 获取数据

  1. dateStart = datetime.date(2005, 4, 8)
  2. dateEnd = datetime.date(2018, 6, 30)
  3. HS300 = ts.get_k_data(code = '000300', index = True, start = f'{dateStart}', end = f'{dateEnd}')
  4. HS300
  1. >>>        date     open    close   high    low     volume      code
  2.     0   2005-04-08  984.66  1003.45 1003.70 979.53  14762500.0  sh000300
  3.     1   2005-04-11  1003.88 995.42  1008.73 992.77  15936100.0  sh000300
  4.     2   2005-04-12  993.71  978.70  993.71  978.20  10226200.0  sh000300
  5.     3   2005-04-13  987.95  1000.90 1006.50 987.95  16071700.0  sh000300
  6.     4   2005-04-14  1004.64 986.97  1006.42 985.58  12945700.0  sh000300
  7.     5   2005-04-15  982.61  974.08  982.61  971.93  10409000.0  sh000300

iii. 绘制价格走势图

  1. xticks = np.arange(start = 0, stop = HS300.shape[0], step = int(HS300.shape[0]/7))
  2. xticklabel = pd.Series(data = HS300.date[xticks])
  3. plt.figure(figsize = (20,5))
  4. fig = plt.axes()
  5. plt.plot(x = np.arange(HS300.shape[0]), y = HS300.close, color = 'red', label = 'close')
  6. fig.set_xticks(xticks)
  7. fig.set_xticklabels(HS300.date[xticks], rotation = 60)
  8. plt.show()

2. 高阶矩和 概念

i. 基本定义

矩，是统计学中的一个常用的指标，用来反映数据分布的形态特点。矩，也被称为“动差”，它代表总体数据中所有变量值与任意一个给定常熟的差的 次方的算术平均数。

矩有原点矩、中心距等不同类型，在平时的统计中，我们比较常用的是原点矩和中心距。

原点矩是检验变量关于 0 的偏离程度，具体定义如下：

中心矩检验的是变量关于 “期望” 的偏离程度，具体定义如下：

其中，

可见，一阶中心距为 0，二阶中心矩就是方差。

ii. Python 定义

  1. def getHighMoment(EarnRate, k, N = 20):
  2.     ```
  3.     EarnRate: 变量x
  4.     k: 所求矩的阶数k
  5.     ```
  6.     HighMoment = (EarnRate**k).rolling(window = N).mean()
  7.     for i in range(1, N-1):
  8.         HighMoment[i] = (EarnRate[:i+1]**k).mean()
  9.     HighMoment.fillna(value = 0)
  11.     return HighMoment
  1. def getEMA(moment, alpha, N = 120):
  2.     ```
  3.     EMA: Exponential Moving Average
  4.     moment:
  5.     alpha: Specify smoothing factor alpha directly, 0 < alpha ≤ 1.
  6.     ```
  7.     EMA = pd.DataFrame.ewm(moment, alpha = alpha, adjust = False).mean()
  9.     return EMA

3. 高阶矩和 实操

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

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

  1. # EarnRate = HS300.close.pct_change(1).fillna(0)
  2. EarnRate = HS300['close'].pct_change(periods = 1).fillna(value = 0)
  3. EarnRate
  4. # .pct_change() 表示当前第 i 个元素与先前的第 i-1 个元素的相差百分比，
  5. # 当然如果指定 periods=n，表示当前元素与先前 n 个元素的相差百分比。
  1. >>> 0       0.000000
  2.     1      -0.008002
  3.     2      -0.016797
  4.     3       0.022683
  5.     4      -0.013917
  6.     5      -0.013060
  7.     ...
  1. Moment_20_2 = getHighMoment(EarnRate = EarnRate, k = 2)
  2. Moment_20_3 = getHighMoment(EarnRate = EarnRate, k = 3)
  3. Moment_20_4 = getHighMoment(EarnRate = EarnRate, k = 4)
  4. Moment_20_5 = getHighMoment(EarnRate = EarnRate, k = 5)
  5. Moment_20_6 = getHighMoment(EarnRate = EarnRate, k = 6)
  6. Moment_20_7 = getHighMoment(EarnRate = EarnRate, k = 7)

五阶矩

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

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

  1. HS = plt.subplot(611)
  2. HS.plot(HS300.close)
  3. HS.set_xticks(xticks)
  4. HS.set_xticklabels([])
  5. ax2 = HS.twinx()  
  6. ax2.plot(Moment_20_2, 'r')
  1. HS = plt.subplot(612)
  2. HS.plot(HS300.close)
  3. HS.set_xticks(xticks)
  4. HS.set_xticklabels([])
  5. ax2 = HS.twinx()  
  6. ax2.plot(Moment_20_3, 'r')
  1. HS = plt.subplot(613)
  2. HS.plot(HS300.close)
  3. HS.set_xticks(xticks)
  4. HS.set_xticklabels([])
  5. ax2 = HS.twinx()  
  6. ax2.plot(Moment_20_4, 'r')
  1. HS = plt.subplot(614)
  2. HS.plot(HS300.close)
  3. HS.set_xticks(xticks)
  4. HS.set_xticklabels([])
  5. ax2 = HS.twinx()  
  6. ax2.plot(Moment_20_5, 'r')
  1. HS = plt.subplot(615)
  2. HS.plot(HS300.close)
  3. HS.set_xticks(xticks)
  4. HS.set_xticklabels([])
  5. ax2 = HS.twinx()  
  6. ax2.plot(Moment_20_6, 'r')
  1. HS = plt.subplot(616)
  2. HS.plot(HS300.close)
  3. HS.set_xticks(xticks)
  4. HS.set_xticklabels(HS300.date[xticks],rotation=60)
  5. ax2 = HS.twinx()  
  6. ax2.plot(Moment_20_7, 'r')
  1. plt.show
  2. <function matplotlib.pyplot.show>
  1. alpha = 0.2,0.4时的EMA
  2. EMA_02 = getEMA(Moment_20_5,0.2)
  3. EMA_04 = getEMA(Moment_20_5,0.4)
  4. HS = plt.subplot(111)
  5. HS.plot(HS300.close)
  6. HS.set_xticks(xticks)
  7. HS.set_xticklabels(HS300.date[xticks],rotation=60)
  8. ax2 = HS.twinx()  
  9. ax2.plot(EMA_02, 'r')
  10. ax3 = ax2.twinx()  
  11. ax3.plot(EMA_04, 'y')
  12. plt.title('Moment-5')
  13. plt.show()

3.回测相关函数定义

切线法寻优函数

输入变量：

输出变量：

  1. def qiexian(profitrate, Moment, EMA):
  2.     flag = np.zeros(len(profitrate))
  3.     cumrate = np.ones(len(profitrate))
  4.     for i in [i for i in range(len(profitrate)) if i > 1]:
  5.         flag[i] = 1 if EMA[i - 1] > EMA[i-2] else flag[i] = -1 
  6.         # if EMA[i - 1] > EMA[i-2]:
  7.         #    flag[i] = 1
  8.         # else:
  9.         #    flag[i] = -1 
  10.     strategy_rate = profitrate * flag + 1
  11.     totalprofit = sum (strategy_rate - 1)  
  13.     return totalprofit

基本高阶矩择时模型中的交易函数
交易逻辑见报告第13页
Sharp：夏普比,年交易天数按250计,MDD：最大回撤率，其他输入输出参数意义同上

  1. def BuySell(profitrate, Moment, EMA):
  3.     flag = np.zeros(len(profitrate))
  4.     cumrate = np.ones(len(profitrate))
  5.     lossflag = 0   
  6.     # 计算单次择时损失，超过10%平仓
  7.     flagloss = 0  
  8.     # 单次亏损超过10%时的信号方向
  10.     # 计算最优alpha
  11.     alpha_all = np.arange(0.05,0.5,0.05)
  13.     for i in range(2,len(profitrate)):
  15.         if i%90 ==0 :
  16.             cumrate_all = [qiexian(list(profitrate[i - 90:i]),Moment[i - 90:i].values,
  17.                                    getEMA(Moment,alpha_all[j])[i - 90:i].values) for j in range(len(alpha_all))]            
  18.             # 取使90天内累计收益率最大的alpha作为后续计算的alpha，并计算相应的EMA           
  19.             alpha = alpha_all[np.argmax(cumrate_all)]
  20.             EMA = getEMA(Moment, alpha)    
  21.             EMA = getEMA(Moment, alpha)
  23.         if lossflag <  -0.1 :
  24.             flag[i] = 0
  25.             flagloss = flag[i-1]
  26.             lossflag = 0
  27.             continue
  29.         if EMA[i - 1] > EMA[i-2] and flagloss != 1:
  30.             flag[i] = 1
  31.             flagloss = 0
  33.         if EMA[i - 1] < EMA[i-2] and flagloss != -1:
  34.             flag[i] = -1
  35.             flagloss = 0
  37.         if flag[i] == flag[i-1] :  
  38.             lossflag = lossflag + profitrate[i]*flag[i]
  39.             lossflag = min(lossflag,0)
  40.         else:
  41.             lossflag = 0
  43.     strategy_rate = profitrate * flag 
  44.     nav = (1+strategy_rate).cumprod()
  45.     cumrate = nav - 1
  46.     totalprofit = nav[len(nav)-1] - 1
  48.     # 交易次数/择时次数
  49.     num = 0
  50.     for i in range(len(flag) - 1):
  51.         if (flag[i+1]!= flag[i]) :
  52.             num+=1
  53.     Sharp = strategy_rate.mean()/strategy_rate.std()* 250**0.5
  54.     MDD = max(1-nav/nav.cummax())
  56.     return cumrate, num, totalprofit, flag, Sharp, MDD

拓展后高阶矩交易模型的交易函数
交易逻辑见报告第14页
k：阈值
其他参数含义同上

  1. def BuySell_1(profitrate,Moment,EMA,k):
  2.     flag = np.zeros(len(profitrate))
  3.     cumrate = np.ones(len(profitrate))
  4.     lossflag = 0   # 计算单次择时损失，超过10%平仓
  5.     flagloss = 0  # 单次亏损超过10%时的信号方向
  7.     # 计算最优alpha
  8.     alpha_all = np.arange(0.05,0.5,0.05)
  9.     for i in  range(2,len(profitrate)):
  11.         if i%90 ==0 :
  12.             cumrate_all = [qiexian(list(profitrate[i - 90:i]),Moment[i - 90:i].values,
  13.                                    getEMA(Moment,alpha_all[j])[i - 90:i].values) for j in range(len(alpha_all))]            
  14.             # 取使90天内累计收益率最大的alpha作为后续计算的alpha，并计算相应的EMA           
  15.             alpha = alpha_all[np.argmax(cumrate_all)]
  16.             EMA = getEMA(Moment,alpha)                   
  17.         if lossflag <  -0.1 :
  18.             flag[i] = 0
  19.             flagloss = flag[i-1]
  20.             lossflag = 0
  21.             continue
  23.         if EMA[i - 1] > EMA[i-2]*(1+k) and flagloss != 1:
  24.             flag[i] = 1
  25.             flagloss = 0
  27.         if EMA[i - 1] < EMA[i-2]*(1-k) and flagloss != -1:
  28.             flag[i] = -1
  29.             flagloss = 0
  31.         if EMA[i - 1] > EMA[i-2]*(1-k) and EMA[i - 1] < EMA[i-2]*(1+k) :
  32.             flag[i] = 0
  33.             flagloss = 0
  35.         if flag[i] == flag[i-1] :  
  36.             lossflag = lossflag + profitrate[i]*flag[i]
  37.             lossflag = min(lossflag,0)
  38.         else:
  39.             lossflag = 0
  41.     strategy_rate = profitrate * flag 
  42.     nav = (1+strategy_rate).cumprod()
  43.     cumrate = nav - 1
  44.     totalprofit = nav[len(nav)-1] - 1
  46.     # 交易次数/择时次数
  47.     num = 0
  48.     for i in range(len(flag) - 1):
  49.         if (flag[i+1]!= flag[i]) :
  50.             num+=1
  51.     Sharp = strategy_rate.mean()/strategy_rate.std()* 250**0.5
  52.     MDD = max(1-nav/nav.cummax())
  54.     return cumrate, num, totalprofit, flag, Sharp, MDD

4.回测并展示结果

回测1：基本的高阶矩择时模型，EMA 初始 alpha=0.4
回测2：阈值k=1%,其他同上
回测3：阈值k=2%,其他同上

  1. cumrate,num,totalprofit,flag,Sharp,MDD = BuySell(EarnRate,Moment_20_5,EMA_04)
  2. cumrate_01,num_01,totalprofit_01,flag_02,Sharp_01,MDD_01 = BuySell_1(EarnRate,Moment_20_5,EMA_04,0.01)
  3. cumrate_02,num_02,totalprofit_02,flag_02,Sharp_02,MDD_02 = BuySell_1(EarnRate,Moment_20_5,EMA_04,0.02)
  1. print('回测1交易次数为:%s' % num)
  2. print('回测2交易次数为:%s' % num_01)
  3. print('回测3交易次数为:%s' % num_02)
  4. print('')
  5. print('回测1总收益率为:%s%%' % (totalprofit*100))
  6. print('回测2总收益率为:%s%%' % (totalprofit_01*100))
  7. print('回测3总收益率为:%s%%' % (totalprofit_02*100))
  8. print('')
  9. print('回测1夏普比为:%s' % Sharp)
  10. print('回测2夏普比为:%s' % Sharp_01)
  11. print('回测3夏普比为:%s' % Sharp_02)
  12. print('')
  13. print('回测1最大回撤为:%s%%' % (MDD*100))
  14. print('回测2最大回撤为:%s%%' % (MDD_01*100))
  15. 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年之后出现了较大回撤。

6. 参考文献

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

添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注