%matplotlib inline
import pandas as pd
import numpy as np
from  tqdm import tqdm
import datetime as dt
from collections import defaultdict
import itertools
import pickle
import warnings
warnings.filterwarnings('ignore')

N=48;n=5
base='000905.XSHG'
START_DATE='2007-1-1'
END_DATE='2019-3-27'

# 生成并保存数据
Codes=get_index_stocks(base,date=END_DATE)

sts=get_extras('is_st', Codes, start_date=START_DATE,\
               end_date=END_DATE,df=True)# 记录ST股票
Codes=sts.T.index[~sts.any(axis=0)].tolist()

def SolveADay(rlist):# 根据一天计算指标
    rlist=np.diff(np.log(rlist))
    RDVar=np.sum(rlist**2)
    RDSkew=np.sqrt(N)*np.sum(rlist**3)/RDVar**1.5
    RDKurt=N*np.sum(rlist**4)/RDVar**2
    return pd.Series([RDVar,RDSkew,RDKurt])

basicFactors=[]
for code in tqdm(Codes):# 按股票代码获取数据
    codePrice=get_price(code, start_date=START_DATE, end_date=END_DATE,skip_paused=True,# 去除停牌
              frequency='5m', fields=['close'])
    codePriceResampled=codePrice['close'].resample('D')
    if (codePriceResampled.count()>0).sum()<242:# 剔除小于一年的数据
        continue
    else:
        codePrice=codePriceResampled.apply(SolveADay).unstack(level=1).dropna()#去除不存在波动的日子
        codePrice.index=pd.MultiIndex.from_product([[code],codePrice.index])
        basicFactors.append(codePrice)

basicFactors=pd.concat(basicFactors)
basicFactors.columns=['RDVar','RDSkew','RDKurt']
basicFactors.index.set_names(['code','date'],inplace=True)
basicFactors.to_csv('zz500c.csv')

basicFactors=pd.read_csv('zz500c.csv',index_col=[0,1])
basicFactors.head()
basicFactors.index.set_levels([
    basicFactors.index.levels[0],
    pd.to_datetime(basicFactors.index.levels[1])
],inplace=True)

gWeek=pd.Grouper(freq='W', level=-1)
# 构建因子指标
RVol=np.sqrt(242/n*basicFactors['RDVar'].groupby(['code',gWeek]).sum())# 计算RVol
RSkew=1/n*basicFactors['RDSkew'].groupby(['code',gWeek]).sum()# 计算RSkew
RKurt=1/n*basicFactors['RDKurt'].groupby(['code',gWeek]).sum()# 计算RKurt
del basicFactors

ax=RVol.hist(bins=50,grid=False)
ax.set_ylabel('Frequency')
ax.set_title('ZZ500 Volatility')

ax=RSkew.hist(bins=50,grid=False)
ax.set_ylabel('Frequency')
ax.set_title('ZZ500 Skewness')

ax=RKurt.hist(bins=50,grid=False)
ax.set_ylabel('Frequency')
ax.set_title('ZZ500 Kurtosis')

def CrossSectionalPercentiles(df):
    df=df.swaplevel().sort_index()
    df=df.groupby('date').apply(lambda x:np.percentile(x,[90,75,50,25,10]))
    df=pd.DataFrame.from_items(zip(df.index, df.values))
    df.index=[90,75,'median',25,10]
    return df.T

CrossSectionalPercentiles(RVol).plot(figsize(8,5))

CrossSectionalPercentiles(RSkew).plot(figsize=(8,5))

CrossSectionalPercentiles(RKurt).plot(figsize(8,5))

# 提取中证500各股周收益率
Codes=list(RVol.index.levels[0])
Price=get_price(Codes, start_date=START_DATE, end_date=END_DATE,
                     skip_paused=False,frequency='1d', fields=['close'])
Price=Price.to_frame()['close']
Price.index.set_names(['date','code'],inplace=True)
Price=np.log(Price.groupby(['code',pd.Grouper(freq='W', level=0)]).last())
Return=np.exp(Price.groupby('code').diff().groupby('code').shift(-1)).swaplevel().sort_index()
del Price

# 提取各股市值
mktCap=[]
for date in Return.index.levels[0]:
    df=get_fundamentals_continuously(
        query(valuation.market_cap).filter(valuation.code.in_(Codes)),
            end_date=date,count=1).to_frame()['market_cap'] 
    df.index.set_levels([date], level=0,inplace=True)
    mktCap.append(df)
mktCap=pd.concat(mktCap)
mktCap.index.set_names(['date','code'],inplace=True)
mktCap.index.set_levels([
    pd.to_datetime(mktCap.index.levels[0]),mktCap.index.levels[1]
],inplace=True)

def GetGroup(df):# 依据指标分组
    df=df.swaplevel().sort_index()
    group=df.groupby('date').apply(lambda x:pd.qcut(x,q=5,precision=10,labels=['Q1','Q2','Q3','Q4','Q5']))
    group.name='group'
    return group
def GroupTestRet(group,mktCap,ret):# 计算分组后各组周收益率
    df=pd.concat([mktCap,group,ret],axis=1)
    df.ffill(inplace=True)
    df['market_cap*return']=df['close']*df['market_cap']
    df=df.groupby(['date','group'])['market_cap*return'].sum()/\
        df.groupby(['date','group'])['market_cap'].sum()
    return df.unstack(level=1)

RVolLevel=GroupTestRet(GetGroup(RVol),mktCap,Return)
RVolLevel.cumprod().plot(figsize=(10,5))

RSkewLevel=GroupTestRet(GetGroup(RSkew),mktCap,Return)
RSkewLevel.cumprod().plot(figsize=(10,5))

RKurtLevel=GroupTestRet(GetGroup(RKurt),mktCap,Return)
RKurtLevel.cumprod().plot(figsize=(10,5))

IC_RSkew=pd.concat([RSkew.swaplevel(),Return],axis=1).\
    groupby('date').corr('pearson').iloc[0::2,-1].reset_index(level=1, drop=True)
IC_RSkew.name='IC'
ax=IC_RSkew.plot(figsize=(16,4),legend=True)
IC_RSkewRolling=IC_RSkew.rolling(12).mean()
IC_RSkewRolling.name='IC average(12 terms)'
IC_RSkewRolling.plot(ax=ax,color='red',legend=True)

aggOperator={'mean':np.mean,
             'std':np.std,
             'min':np.min,
             'max':np.max,
             'negative percentile':lambda x:'%.2f'%(np.sum(x<0)/x.shape[0]*100)+'%'}
IC_RSkewSummary=IC_RSkew.agg(aggOperator)
IC_RSkewSummary=pd.DataFrame(IC_RSkewSummary).T.\
                applymap(lambda x:np.round(x,3) if isinstance(x,float) else x)
IC_RSkewSummary

IC_RSkewSummary=IC_RSkew.resample('Y').agg(aggOperator).unstack(level=1).T
IC_RSkewSummary.set_index(IC_RSkewSummary.index.year,inplace=True)
IC_RSkewSummary=IC_RSkewSummary[list(aggOperator)].\
                applymap(lambda x:np.round(x,3) if isinstance(x,float) else x)
IC_RSkewSummary

def GroupTestRetLongShortCumprod(group,mktCap,ret):# 多空策略分组累计收益率
    df=pd.concat([mktCap,group,ret],axis=1)
    df.ffill(inplace=True)
    df['market_cap*return']=df['close']*df['market_cap']
    df=df.groupby(['date','group'])['market_cap*return'].sum()/\
        df.groupby(['date','group'])['market_cap'].sum()
    df=df.unstack(level=1)
    df.columns=list(df.columns)
    df['Q1-Q5']=1/2*(df['Q1']-df['Q5'])+1
    df=df['Q1-Q5'].cumprod()
    return df

# zz500Ret 周度
zz500Ret=get_price(base,start_date=START_DATE, end_date=END_DATE,fields=['close'])
zz500Ret=zz500Ret.resample('W').last()
zz500Ret=(np.exp(np.log(zz500Ret).diff().shift(-1)))
zz500Ret.ffill(inplace=True)

RSkewLevel.columns=list(RSkewLevel.columns)
RSkewLevel['Q1-Q5_cum']=GroupTestRetLongShortCumprod(GetGroup(RSkew),mktCap,Return)
RSkewLevel['zz500']=zz500Ret

ax=RSkewLevel[['Q1','Q5','zz500']].cumprod().plot(figsize=(10,5))
RSkewLevel['Q1-Q5_cum'].plot(secondary_y=True,ax=ax,label='Q1-Q5',legend=True)

def BacktrackInformation(level,ret):#计算回测所需数据
    # 累计收益率
    accRate=level.resample('Y').last()/\
        level.resample('Y').first()-1
    # 年化收益率
    yldRate=(1+accRate)**(242/5/level.resample('Y').count())-1
    yldRate[-1]=(1+accRate[-1])**(60.5/5/level.resample('Y').count()[-1])-1#2019调整
    # 年化波动率
    sigma=np.sqrt(242/5)*ret.resample('Y').std()
    sigma[-1]*=np.sqrt(60.5/5)/np.sqrt(242/5)#2019调整
    # 最大回撤
    maxDropdown=(-1+np.exp(np.log(level).diff())).resample('Y').max()
    # 信息比率
    informationRatio=yldRate/sigma
    # 整体
    accRateTotal=level.loc[level.last_valid_index()]/level.loc[level.first_valid_index()]-1
    yldRateTotal=(1+accRateTotal)**(242/5/level.count().sum())-1
    sigmaTotal=np.sqrt(242/5)*ret.std()
    maxDropdownTotal=max(maxDropdown)
    informationRatioTotal=yldRateTotal/sigmaTotal
    return accRate,maxDropdown,yldRate,sigma,informationRatio,\
            accRateTotal,maxDropdownTotal,yldRateTotal,sigmaTotal,informationRatioTotal
def BacktrackTable(accRate,maxDropdown,yldRate,sigma,informationRatio,\
                   accRateTotal,maxDropdownTotal,yldRateTotal,sigmaTotal,informationRatioTotal):
    #罗列回测信息
    cols=['accRate','maxDropdown','yldRate','sigma','informationRatio']
    df=pd.DataFrame([accRate,maxDropdown,yldRate,sigma,informationRatio],index=cols).T
    years=list(df.index.year)
    df=df.append(pd.DataFrame(
        [accRateTotal,maxDropdownTotal,yldRateTotal,sigmaTotal,informationRatioTotal],
        index=cols).T)
    df=df.apply(lambda x:list(map(lambda z:"%.2f"%(z*100)+'%',x[:-1]))+["%.3f"%x[-1:]],axis=1)
    df=pd.DataFrame(df.values.tolist(),columns=cols,index=years+['total'])
    return df

ret=1/2*(RSkewLevel['Q1']-RSkewLevel['Q5'])
level=RSkewLevel['Q1-Q5_cum']
BacktrackTable(*BacktrackInformation(level,ret))

ax=RSkewLevel[['Q1','zz500']].cumprod().plot(figsize=(16,5))
RSkewLevel['Q1-zz500']=RSkewLevel['Q1']-RSkewLevel['zz500']
RSkewLevel['Q1-zz500'].plot(ax=ax,secondary_y=True,legend=True)

ret=RSkewLevel['Q1']
RSkewLevel['Q1_cum']=RSkewLevel['Q1'].cumprod()
level=RSkewLevel['Q1_cum']
BacktrackTable(*BacktrackInformation(level,ret))

ret=1/2*(RSkewLevel['Q1']-RSkewLevel['zz500'])
RSkewLevel['Q1-zz500']=(ret+1).cumprod()
level=RSkewLevel['Q1-zz500']
BacktrackTable(*BacktrackInformation(level,ret))

RSkewGroup=GetGroup(RSkew)
Q1_RSkewGroup=RSkewGroup.unstack(level=1)
Q1_RSkewGroup=((Q1_RSkewGroup=='Q1')==True).astype(int)
Q1_mktCap=(mktCap.unstack()*Q1_RSkewGroup)
aggOperator={'mean':np.mean,'max':np.max,'min':np.min,'std':np.std,'acc':np.sum}

Q1_RSkewChangeBuy=(Q1_RSkewGroup.diff()>0).astype(int) #buy in
Q1_RSkewChangeBuy=(Q1_RSkewChangeBuy*Q1_mktCap).sum(axis=1)/Q1_mktCap.sum(axis=1)
Q1_RSkewChangeBuy=Q1_RSkewChangeBuy.shift(-1).iloc[:-1]

Q1_RSkewChangeSell=(Q1_RSkewGroup.diff()<0).astype(int).shift(-1) #sell out
Q1_RSkewChangeSell=(Q1_RSkewChangeSell*Q1_mktCap).sum(axis=1)/Q1_mktCap.sum(axis=1)
Q1_RSkewChangeSell=Q1_RSkewChangeSell.iloc[:-1]

Q1_RSkewChange=pd.DataFrame([Q1_RSkewChangeBuy,Q1_RSkewChangeSell]).min(axis=0)
Q1_RSkewChange=(Q1_RSkewChange.resample('Y').agg(aggOperator)).unstack(level=1).T
Q1_RSkewChange.index=Q1_RSkewChange.index.year

Q1_RSkewChange[list(aggOperator)]