-作者：JoeyQ

在上节  中，我们简单介绍了统计套利策略的原理，在本节中，我们将在沪深300各个行业中寻找可以配对的股票，并以此构造投资组合。投资组合更新的频率为每季度，我们主要步骤如下：

1.划分行业：在本研究中，我们主要研究的是银行，券商，医药，钢铁，煤炭，保险，房地产，家电，汽车制造这几个行业，概含了沪深300指数主要的行业

2.收益率相关性：我们选取的相关性阙值为0.7，太高导致价差波动小，太低则均值回归力度较弱。我们筛选出上一季度日收益率相关系数大于0.7的股票对作为备选股票

3.平稳性分析：在满足收益率相关系数大于0.7后，我们计算出价差，同时要求价差在0.05水平下通过单位根检验，这是为了防止上节中价差呈现单边走势的情况，我们要求价差必须围绕0轴波动

4.价差形态过滤：我们希望价差回归均值的速度较快，但由于我们调仓频率为一个季度，时间足够长，因此不需要通过这一步检验

通过以上我们获得所有满足条件的股票对，并按照上一节的股票对交易策略进行交易，即当价差大于1.5时做多以单位A，同时做空$\beta$单位的B，当价差小于1.5时相反，当价差归零时清仓。我们测试了2007年至2017年的收益情况。结果如下：

可以发现我们获得一个很好的收益曲线！夏普值高达6.95！

我们知道金融类股票的波动较小，因此我们剔除了银行，保险和券商股票，再构造投资组合，表现如下：

累计收益率显著提高，回撤控制在较小范围！

改进策略思考¶

我们知道有左侧交易和右侧交易，左侧偏向反转，右侧偏向趋势，我们考虑建立配对交易策略的右侧交易策略，即当价差从上往下穿过1.5阙值线时建仓，同时，当价差趋向0时，回归均值的力度已经转弱，我们考虑价差小于0.5时清仓。具体流程如下：

最后的最后，我们集中比较一下两种策略以及分行业的股票套利策略表现：

以下是近十年两种策略分行业的年度收益情况：

反思：¶

在这两节中，我们构造了统计套利之股票配对交易策略，并进行了沪深300行业的回测。我们假设所有的沪深300均可融资融券，从表现看我们的策略获得一个较为满意的收益，有效对冲了大盘下跌的风险。然而现实中并不是每只股票均可融资融券，因此这里可能存在未被市场参与者未能获取的超额收益。下一步，我们可以尝试针对可以融资融券的股票进行统计套利的回测！

import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom scipy import statsimport statsmodels.api as smimport statsmodels.tsa.stattools as ts

##价差函数def get_jiacha(securityy, securityx,start_date,end_date):pricex=get_price(securityx,start_date,end_date,fields='close',fq='pre')['close']pricey=get_price(securityy,start_date,end_date,fields='close',fq='pre')['close']model = sm.OLS(log(pricey),log(pricex)).fit()slope=model.params[0]st1=log(pricey)-slope*log(pricex)st=(st1-mean(st1))/std(st1)return st,slope

def get_peidui(code_index,stocks,start_date,end_date,corrindex=0.7,stationindex=0.05): ### 获得配对股票 #####输入：行业或者概念的代码，总的股票池，开始时间，结束时间，相关系数阙值，平稳性阙值##输出：所有的配对股票，按照行业分类hangye=[]    for i in arange(len(code_index)):tempt99=[i99 for i99 in get_concept_stocks(code_index[i]) if i99 in stocks]hangye.append(tempt99)store=[]allstore=[]##储存所有配对股票    for i56 in arange(len(hangye)):price=get_price(hangye[i56],start_date,end_date,fields='close')['close']##相关系数阙值for i in arange(len(price.columns)-1):for j in arange(i+1,len(price.columns)):if corrcoef(price[price.columns[i]],price[price.columns[j]])[0][1]>corrindex and \ts.adfuller(get_jiacha(price.columns[i],price.columns[j],start_date,end_date)[0],0)[1]<stationindex: ###平稳性检验以及相关性检验store.append([price.columns[i],price.columns[j]])allstore=allstore+storereturn allstore

def get_cumrate(allstore,start_date,end_date,jiacha_index=1.5):### 计算累计净值 #####输入：allstore所有的配对股票，开始时间，结束时间##输出：每个配对股票的净值bigstore=pd.DataFrame()if len(allstore)<1:bigstore['0']=[1 for i in arange(len(get_price('600001.XSHG',start_date,end_date))-2)]else:for j in arange(len(allstore)):##allstore长度price=get_price(allstore[j],start_date,end_date,fields='close',fq='pre')['close']st,slope=get_jiacha(allstore[j][0],allstore[j][1],start_date,end_date)u=0sec1=[] sec2=[]sec1num=0sec2num=0for i in arange(1,len(st)):if u==0 and st[i]>jiacha_index and st[i]<3:sec1num=-1##做空sec1sec2num=-sec1num*slope ##做多sec2u=ielif (u<>0 and st[i]*st[i-1]<0) or (abs(st[i])>3): ##此时价差变0清仓,或者价差偏离过大达到3清仓sec1num=0sec2num=0u=0elif u==0 and st[i]<-jiacha_index and st[i]>-3:sec1num=1sec2num=-sec1num*slopeu=isec1.append(sec1num)sec2.append(sec2num) #储存的是建仓信号，长度为price-1sec1rate=price[allstore[j][0]].diff(1)[1:]/list(price[allstore[j][0]][:-1])sec2rate=price[allstore[j][1]].diff(1)[1:]/list(price[allstore[j][1]][:-1])sumrate=sec1rate[1:]*sec1[:-1]+sec2rate[1:]*sec2[:-1] #发出信号后一天建仓cumrate=[]tempt=1for i in arange(len(sumrate)):tempt=tempt*(1+sumrate[i])if sumrate[i-1]<>0:tempt=tempt*(1-0.01/255)##融券费用if i >=1:if sumrate[i-1]<>0 and sumrate[i]==0:tempt=tempt*(1-0.001-0.00025) ##平仓手续费   cumrate.append(tempt)#plt.plot(cumrate)#plt.title('配对策略净值变化曲线')bigstore[j]=cumratereturn bigstore

def trade(code_index,stocks,time,plot_index=True,print_index=True):###交易函数 #####输入：行业代码，股票池，时间##输出：净值曲线，总收益，年化收益，夏普比，最大回撤netvalue=[]##储存净值tempt22=1nianhua=[]##储存年末净值for ii11 in arange(len(time)-2):allstore=get_peidui(code_index,stocks,start_date=time[ii11],end_date=time[ii11+1])##获得配对bigstore=get_cumrate(allstore,start_date=time[ii11+1],end_date=time[ii11+2])tempt23=[ii77*tempt22 for ii77 in list(bigstore.mean(axis=1))]netvalue=netvalue+tempt23tempt22=netvalue[-1]nianhua.append(netvalue[-1])###绘制净值图像if plot_index==True:plt.figure(figsize=(20,6))plt.plot(netvalue)timestamp=[time[i*4][:4] for i in arange(len(time)/4)]plt.xticks(np.arange(0,len(netvalue),len(netvalue)/len(timestamp)),timestamp)plt.grid(True)plt.title('配对策略净值')plt.xlabel('年份')plt.ylabel('净值')###计算最大回撤，平均年化收益率，夏普比tt1=[]tt2=[]for i in arange(len(netvalue)):for j in arange(i+1,len(netvalue)):tt1.append((netvalue[i]-netvalue[j])/netvalue[i])tt2.append(max(tt1))huiche=max(tt2)*100 #最大回撤yearrate=[nianhua[i]/nianhua[i-1]-1 for i in arange(1,len(nianhua))]shouyi=[netvalue[ii88+1]/netvalue[ii88]-1 for ii88 in arange(len(netvalue)-1) ]#每日收益sharpe=(mean(shouyi)-0.0285/252)/std(shouyi)*np.sqrt(252)#夏普值#sharpe=(mean(yearrate)-0.04)/std(yearrate) pjnh=((netvalue[-1])**(1./(len(time)/4))-1)*100 #平均年化ljnh=(netvalue[-1]-1)*100 #累计收益##打印if print_index==True:print '-'*30print '%14s %15.5f' % ('累计收益%：',ljnh)print '%14s %15.5f' % ('平均年化%:',pjnh)print '%14s %15.5f' % ('夏普比：' , sharpe)print '%14s %15.5f' % ('最大回撤%:', huiche)##返回一些数值进行对比:累计收益，平均年化，夏普值，最大回撤,年末净值return [ljnh,pjnh,sharpe,huiche,nianhua]

##获取沪深300分行业股票数据##券商:GN780,银行：GN815，房地产GN733，保险GN646，汽车制造GN835，煤炭GN823，钢铁GN727，白酒GN705，家电GN706，化学制药GN750stocks = get_index_stocks('000300.XSHG')hangye=[]code_index=['GN780','GN815','GN733','GN646','GN835','GN823','GN727','GN705','GN706','GN750']for i in arange(len(code_index)):tempt99=[i99 for i99 in get_concept_stocks(code_index[i]) if i99 in stocks]hangye.append(tempt99)

###按照季度生成时间time=[]for i in arange(2007,2018):for j in arange(1,13,3):time.append(str(i)+'-'+str(j)+'-01')

##储存输出指标storeindex=pd.DataFrame()##进行交易storeindex['全部行业']=trade(code_index,stocks,time)

累计收益%：     35998.35195
平均年化%:        70.80432
  夏普比：         6.94854
最大回撤%:         7.31672

##根据经验我们知道银行，保险，券商等金融类股票波动较小，我们把他们剔除后看他们的配对交易收益storeindex['非金融']=trade(['GN733','GN835','GN823','GN727','GN705','GN706','GN750'],stocks,time)

累计收益%：    182156.93730
平均年化%:        97.89090
  夏普比：         5.27588
最大回撤%:        13.94375

##分行业计算相关指标storeindex['券商']=trade(['GN780'],stocks,time,print_index=False,plot_index=False)storeindex['银行']=trade(['GN815'],stocks,time,print_index=False,plot_index=False)storeindex['房地产']=trade(['GN733'],stocks,time,print_index=False,plot_index=False)storeindex['保险']=trade(['GN646'],stocks,time,print_index=False,plot_index=False)storeindex['汽车制造']=trade(['GN835'],stocks,time,print_index=False,plot_index=False)storeindex['煤炭']=trade(['GN823'],stocks,time,print_index=False,plot_index=False)storeindex['钢铁']=trade(['GN727'],stocks,time,print_index=False,plot_index=False)storeindex['白酒']=trade(['GN705'],stocks,time,print_index=False,plot_index=False)storeindex['家电']=trade(['GN706'],stocks,time,print_index=False,plot_index=False)storeindex['化学制药']=trade(['GN750'],stocks,time,print_index=False,plot_index=False)

下面我们尝试改进一下策略，由于价差刚开始大于阙值时可能会会继续增大导致回撤，我们设置新的策略是当价差从大于3变成小于3时买入，同样当价差接近0时持有的价值已经较小，因此我们在价差从大于0.5变成小于0.5后卖出。这其实是把左侧交易变成右侧交易。¶

def get_cumrate2(allstore,start_date,end_date,jiacha_index=1.5):### 计算累计净值 #####输入：allstore所有的配对股票，开始时间，结束时间##输出：每个配对股票的净值bigstore=pd.DataFrame()if len(allstore)<1:bigstore['0']=[1 for i in arange(len(get_price('600001.XSHG',start_date,end_date))-2)]else:for j in arange(len(allstore)):##allstore长度price=get_price(allstore[j],start_date,end_date,fields='close',fq='pre')['close']st,slope=get_jiacha(allstore[j][0],allstore[j][1],start_date,end_date)u=0sec1=[] sec2=[]sec1num=0sec2num=0for i in arange(1,len(st)):if u==0 and st[i-1]>jiacha_index and st[i]<jiacha_index: ##从上往下穿过1.5时开仓sec1num=-1##做空sec1sec2num=-sec1num*slope ##做多sec2u=ielif (u<>0 and abs(st[i-1])>0.5 and abs(st[i])<0.5) or (abs(st[i])>3): ##此时价差变0清仓,或者价差偏离过大达到3清仓sec1num=0sec2num=0u=0elif u==0 and st[i-1]<-jiacha_index and st[i]>-jiacha_index: ##从下往上穿过1.5时开仓sec1num=1sec2num=-sec1num*slopeu=isec1.append(sec1num)sec2.append(sec2num) #储存的是建仓信号，长度为price-1sec1rate=price[allstore[j][0]].diff(1)[1:]/list(price[allstore[j][0]][:-1])sec2rate=price[allstore[j][1]].diff(1)[1:]/list(price[allstore[j][1]][:-1])sumrate=sec1rate[1:]*sec1[:-1]+sec2rate[1:]*sec2[:-1] #发出信号后一天建仓cumrate=[]tempt=1for i in arange(len(sumrate)):tempt=tempt*(1+sumrate[i])if i >=1:if sumrate[i-1]<>0:tempt=tempt*(1-0.01/250)##融券费用if sumrate[i]==0:tempt=tempt*(1-0.001-0.00025)    ##交易费用cumrate.append(tempt)#plt.plot(cumrate)#plt.title('配对策略净值变化曲线')bigstore[j]=cumratereturn bigstore

def trade2(code_index,stocks,time,plot_index=True,print_index=True):###交易函数2 #####输入：行业代码，股票池，时间##输出：净值曲线，总收益，年化收益，夏普比，最大回撤netvalue=[]##储存净值tempt22=1nianhua=[]##储存年末净值for ii11 in arange(len(time)-2):allstore=get_peidui(code_index,stocks,start_date=time[ii11],end_date=time[ii11+1])##获得配对bigstore=get_cumrate2(allstore,start_date=time[ii11+1],end_date=time[ii11+2])tempt23=[ii77*tempt22 for ii77 in list(bigstore.mean(axis=1))]netvalue=netvalue+tempt23tempt22=netvalue[-1]nianhua.append(netvalue[-1])###绘制净值图像if plot_index==True:plt.figure(figsize=(20,6))plt.plot(netvalue)timestamp=[time[i*4][:4] for i in arange(len(time)/4)]plt.xticks(np.arange(0,len(netvalue),len(netvalue)/len(timestamp)),timestamp)plt.grid(True)plt.title('配对策略净值')plt.xlabel('年份')plt.ylabel('净值')###计算最大回撤，平均年化收益率，夏普比tt1=[]tt2=[]for i in arange(len(netvalue)):for j in arange(i+1,len(netvalue)):tt1.append((netvalue[i]-netvalue[j])/netvalue[i])tt2.append(max(tt1))huiche=max(tt2)*100 #最大回撤yearrate=[nianhua[i]/nianhua[i-1]-1 for i in arange(1,len(nianhua))]shouyi=[netvalue[ii88+1]/netvalue[ii88]-1 for ii88 in arange(len(netvalue)-1) ]#每日收益sharpe=(mean(shouyi)-0.0285/252)/std(shouyi)*np.sqrt(250)#夏普值pjnh=((netvalue[-1])**(1./(len(time)/4))-1)*100 #平均年化ljnh=(netvalue[-1]-1)*100 #累计收益##打印if print_index==True:print '-'*30print '%14s %15.5f' % ('累计收益%：',ljnh)print '%14s %15.5f' % ('平均年化%:',pjnh)print '%14s %15.5f' % ('夏普比：' , sharpe)print '%14s %15.5f' % ('最大回撤%:', huiche)##返回一些数值进行对比:累计收益，平均年化，夏普值，最大回撤,年净值return [ljnh,pjnh,sharpe,huiche,nianhua]

##进行新的策略交易storeindex['全部行业2']=trade2(code_index,stocks,time,print_index=False,plot_index=False)storeindex['非金融2']=trade2(['GN733','GN835','GN823','GN727','GN705','GN706','GN750'],stocks,time,print_index=False,plot_index=False)storeindex['券商2']=trade2(['GN780'],stocks,time,print_index=False,plot_index=False)storeindex['银行2']=trade2(['GN815'],stocks,time,print_index=False,plot_index=False)storeindex['房地产2']=trade2(['GN733'],stocks,time,print_index=False,plot_index=False)storeindex['保险2']=trade2(['GN646'],stocks,time,print_index=False,plot_index=False)storeindex['汽车制造2']=trade2(['GN835'],stocks,time,print_index=False,plot_index=False)storeindex['煤炭2']=trade2(['GN823'],stocks,time,print_index=False,plot_index=False)storeindex['钢铁2']=trade2(['GN727'],stocks,time,print_index=False,plot_index=False)storeindex['白酒2']=trade2(['GN705'],stocks,time,print_index=False,plot_index=False)storeindex['家电2']=trade2(['GN706'],stocks,time,print_index=False,plot_index=False)storeindex['化学制药2']=trade2(['GN750'],stocks,time,print_index=False,plot_index=False)

storeindex1=storeindex

##更改行名，同时计算每年的收益storeindex1.index=['累计收益%','平均年化%','夏普值','最大回撤%','年净值']##计算每年收益shinianrate=pd.DataFrame()for ii67 in arange(len(storeindex1.columns)):    uiu=[(storeindex1.iloc[-1][ii67][-1-4*i]/storeindex1.iloc[-1][ii67][-1-4*(i+1)]-1)*100 for i in arange(10)] #2008-2018年化收益shinianrate[storeindex1.columns[ii67]]=uiushinianrate.index=[str(ii864) for ii864 in arange(2017,2007,-1)]

###比较一下改进的策略和未改进策略##生成交错的index【0,12,1,13...】tempt776=[]for i776 in arange(12):tempt776=tempt776+[i776,i776+12]

storeindex1[:-1].iloc[:,tempt776[15:24]]

###近十年各行业标准策略每年收益率情况shinianrate.iloc[:,:11]

###近十年各行业每年收益率情况shinianrate.iloc[:,12:]