from jqdata import *
from jqlib.technical_analysis import *
from jqfactor import *
from scipy import stats
from statsmodels import regression
import datetime
import pandas as pd
import numpy as np
import statsmodels.api as sm
import warnings
warnings.filterwarnings('ignore')


#获取每月最后一个交易日（其余频次是获取第一个交易日，暂未改动）
def get_tradeday_list(start,end,frequency=None,count=None):
    if count != None:
        df = get_price('000001.XSHG',end_date=end,count=count)
    else:
        df = get_price('000001.XSHG',start_date=start,end_date=end)
    if frequency == None or frequency =='day':
        return df.index
    else:
        df['year-month'] = [str(i)[0:7] for i in df.index]
        if frequency == 'month':
            return df.drop_duplicates('year-month', keep = 'last').index
        elif frequency == 'quarter':
            df['month'] = [str(i)[5:7] for i in df.index]
            df = df[(df['month']=='01') | (df['month']=='04') | (df['month']=='07') | (df['month']=='10') ]
            return df.drop_duplicates('year-month').index
        elif frequency =='halfyear':
            df['month'] = [str(i)[5:7] for i in df.index]
            df = df[(df['month']=='01') | (df['month']=='07')]
            return df.drop_duplicates('year-month').index 

#获取股票池
def get_stock(stockPool, begin_date):
    if stockPool == 'HS300':
        stockList = get_index_stocks('000300.XSHG', begin_date)
    elif stockPool == 'ZZ500':
        stockList = get_index_stocks('399905.XSHE', begin_date)
    elif stockPool == 'ZZ800':
        stockList = get_index_stocks('399906.XSHE', begin_date)   
    elif stockPool == 'A':
        stockList = list(get_all_securities(['stock'], date = begin_date).index)
    
    #剔除ST股
    st_data = get_extras('is_st', stockList, count = 1,end_date = begin_date)
    stockList = [stock for stock in stockList if not st_data[stock][0]]
    #剔除*st股票
    stockList = [stock for stock in stockList if '*' not in get_security_info(stock).display_name]
    #剔除上市不足6个月的新股
    stockList = delete_new(stockList, begin_date, n = 183)
    #剔除停牌
    suspended_info_df = get_price(stockList, end_date = begin_date, count = 1, frequency = 'daily', fields = 'paused')['paused']
    stockList = [stock for stock in stockList if suspended_info_df[stock][0] == 0]
    
    return stockList

#剔除新股
def delete_new(stocks, beginDate, n = 365):
    stockList = []
    for stock in stocks:
        start_date = get_security_info(stock).start_date
        if start_date < (beginDate - datetime.timedelta(days = n)):
            stockList.append(stock)
    return stockList

def get_factor_data(stocks, start_date, end_date):
    if factor == 'IVol':
        factor_data = FF3(stocks, start_date, end_date)
    elif factor == 'IVR':
        factor_data = FF3(stocks, start_date, end_date)
    elif factor == 'PriceDelay':
        factor_data = PD(stocks, start_date, end_date)
    elif factor == 'adjTurnover':
        factor_data = adjTurnover(stocks, end_date)
        
    return factor_data

#求Fama-French三因子模型特质波动率
def FF3(stocks, start_date, end_date, rf = 0.04):
    LoS=len(stocks)
    #查询三因子的语句
    q = query(
        valuation.code,
        valuation.circulating_market_cap,
        (balance.total_owner_equities/valuation.circulating_market_cap/100000000.0).label("BP"),
        #indicator.roe,
        #balance.total_assets.label("Inv")
    ).filter(
        valuation.code.in_(stocks)
    )
    
    df = get_fundamentals(q, start_date)
    df.index = df['code']
    del df['code']
    #中性化
    #df = neutralize(df, how = ['sw_l1', 'market_cap'], date = end_date, axis = 0)
    n = int(LoS/3)

    #选出特征股票组合
    S=df.sort_values('circulating_market_cap').index.tolist()[: n]
    B=df.sort_values('circulating_market_cap').index.tolist()[LoS - n:]
    L=df.sort_values('BP').index.tolist()[: n]
    H=df.sort_values('BP').index.tolist()[LoS - n:]
    
    df5 = df['circulating_market_cap']
    
    # 获得样本期间的股票价格并计算日收益率
    df2 = get_price(stocks, start_date = start_date, end_date = end_date, fields=['close'])['close']
    df4 = df2.pct_change()
    df4.dropna(how ='all', inplace = True)
    df4.fillna(0, inplace = True)
    #求因子的值，按流通市值加权
    SMB = list(np.dot(df4[S], df5.loc[S] / df5.loc[S].sum()) - np.dot(df4[B], df5.loc[B] / df5.loc[B].sum()))
    HML = list(np.dot(df4[H], df5.loc[H] / df5.loc[H].sum()) - np.dot(df4[L], df5.loc[L] / df5.loc[L].sum()))
    
    #用股票池，流通市值为权重作为市场基准
    df6 = df5.loc[df4.columns]
    df6.fillna(df5.mean(), inplace = True)
    RM = list(np.dot(df4, df6 / df6.sum()) - rf/252)
    if len(SMB) > len(RM):
        SMB.drop(SMB.index[0], inplace = True)
        HML.drop(HML.index[0], inplace = True)

    #将因子们计算好并且放好
    X = pd.DataFrame({"RM":RM, "SMB":SMB, "HML":HML})
    # 对样本数据进行线性回归并计算残差
    t_scores=[0.0] * LoS
    for i in range(LoS):
        t_stock = stocks[i]
        t_r = linreg(X, df4[t_stock] - rf/252)
        t_scores[i] = t_r
    
    #这个scores就是残差 
    scores = pd.DataFrame({'score': t_scores}, index = stocks)
    
    return scores

#9
# 辅助线性回归的函数
# 输入:X:回归自变量 Y:回归因变量 完美支持list,array,DataFrame等三种数据类型
# 输出:参数估计结果-list类型
def linreg(X,Y):
    X=sm.add_constant(array(X, dtype=float))
    Y=array(Y)
    results = sm.OLS(Y, X).fit()
    
    if factor == 'IVol':
        return results.resid.std() * sqrt(252)
    elif factor == 'IVR':
        return (1 - results.rsquared)
    elif factor == 'PriceDelay':
        return results.rsquared
    elif factor == 'adjTurnover':
        return results.resid

#计算价格时滞
def PD(stocks, minus1_date, end_date, rf = 0.04):
    
    X1 = MKT(stocks, minus1_date, end_date, n = 0)
    X2 = MKT(stocks, minus1_date, end_date, n = 3)
    quote = get_price(stocks, start_date = minus1_date, end_date = end_date, fields=['close'])['close']
    ret = quote.pct_change()
    ret.dropna(how ='all', inplace = True)
    ret.fillna(0, inplace = True)
    ret = ret - rf/252
    #OLS计算hetero
    PriceDelay = pd.DataFrame(index = ['PriceDelay'])
    for stock in ret.columns:
        Y1 = ret[stock][:len(X1)]
        Y2 = ret[stock][:len(X2)]
        r0 = linreg(X1 - rf/252, Y1)
        r3 = linreg(X2 - rf/252, Y2)
        PriceDelay[stock] = 1 - r0/r3
    
    PriceDelay = PriceDelay.T

    #返回特质波动率vol
    return PriceDelay

#计算滞后n期的市场因子
def MKT(stocks, minus1_date, end_date, n):
    if n == 0:
        period = [minus1_date, end_date]
    else:
        period = list(get_tradeday_list(start = None, end = minus1_date, frequency = 'month', count = n * 24).date)
        period.append(end_date)
    
    X = pd.DataFrame(columns = period[1:], index = range(24))
    for i in range(len(period) - 1):
        #设置统计范围
        quote = get_price(stocks, start_date = period[i], end_date = period[i + 1], fields=['close'])['close']
        ret = quote.pct_change()
        ret.dropna(how ='all', inplace = True)
        ret.fillna(0, inplace = True)
        #构造市场基准收益：流通市值加权
        q = query(valuation.circulating_market_cap, valuation.code).filter(valuation.code.in_(stocks))
        df = get_fundamentals(q, period[i])
        df.index = df['code']
        del df['code']
        #df = (df * 1e8).apply(np.log)
        df = df/df.sum()
        
        if len(ret.columns) > len(df):
            ret = ret[df.index.tolist()]
        elif len(ret.columns) < len(df):
            df = df.loc[ret.columns.tolist()]
        else:
            pass
        
        final = pd.DataFrame(np.dot(ret, df), columns = ['final'])
        X[period[i + 1]] = final['final']
        
    X.dropna(how = 'any', inplace = True)
    
    return X

#市值调整换手
def adjTurnover(stocks, end_date):
    VOL20 = get_factor_values(stocks, ['VOL20'], end_date = end_date, count = 1)['VOL20'].T
    VOL20.replace(0, np.nan, inplace = True)
    VOL20 = VOL20.apply(np.log)
    VOL20.columns = ['VOL20']
    q = query(valuation.code, valuation.circulating_market_cap).filter(valuation.code.in_(stocks))
    df = get_fundamentals(q, end_date)
    df.index = df['code']
    del df['code']
    VOL20['circulating_market_cap'] = (df['circulating_market_cap'] * 1e8).apply(np.log)
    VOL20.dropna(how = 'any', inplace = True)
    

    X = VOL20['circulating_market_cap']
    Y = VOL20['VOL20']
    l = linreg(X,Y)

    adjTurnover = pd.DataFrame(l, index = VOL20.index, columns = ['adjTurnover'])
    
    return adjTurnover

#因子值数据处理，最后一个参数决定是否做市值行业中性化，默认做
def data_preprocessing(factor_data, stockList, industry_code, date, ind_neutral = True):
    #去极值，参数来自光大研报
    factor_data = winsorize_med(factor_data, scale = 3, inf2nan = False, axis = 0)
    #缺失值处理
    nan_ratio = 1 - factor_data.count() / len(factor_data)
    if nan_ratio[0] > 0.2:
        print(str(date) + '数据缺失率为：' + str(nan_ratio))
    factor_data = replace_nan_indu(factor_data, stockList, industry_code, date)
    #标准化处理
    factor_data = standardlize(factor_data, axis = 0)
    #中性化处理
    if ind_neutral == True:
        factor_data = neutralize(factor_data, how = ['sw_l1', 'market_cap'], date = date, axis = 0)
    
    return factor_data

#缺失值处理
def replace_nan_indu(factor_data, stockList, industry_code, date):
    #把nan用行业中位数代替，依然会有nan，此时用所有股票中位数代替
    i_Constituent_Stocks = {}
    data_temp = pd.DataFrame(index = industry_code, columns = factor_data.columns)
    for i in industry_code:
        temp = get_industry_stocks(i, date)
        i_Constituent_Stocks[i] = list(set(temp).intersection(set(stockList)))
        data_temp.loc[i] = nanmedian(array(factor_data.loc[i_Constituent_Stocks[i], :]))
    for factor in data_temp.columns:
        #行业缺失值用所有行业中位数代替
        null_industry = list(data_temp.loc[pd.isnull(data_temp[factor]), factor].keys())
        for i in null_industry:
            data_temp.loc[i, factor] = nanmedian(array(data_temp[factor]))
        #查询空值所在位置并填充
        null_stock = list(factor_data.loc[pd.isnull(factor_data[factor]), factor].keys())
        for i in null_stock:
            industry = get_industry_name(i_Constituent_Stocks, i)
            if industry:
                factor_data.loc[i, factor] = data_temp.loc[industry[0], factor] 
            else:
                factor_data.loc[i, factor] = nanmedian(array(factor_data[factor]))
    
    return factor_data

#取股票对应行业
def get_industry_name(i_Constituent_Stocks, value):
    return [k for k, v in i_Constituent_Stocks.items() if value in v]

#获取所需数据并存储
def get_df(period):
    #申万更新以前
    industry_old_code=['801010','801020','801030','801040','801050','801080','801110','801120','801130','801140','801150',\
                    '801160','801170','801180','801200','801210','801230']
    #申万更新以后
    industry_new_code = ['801010','801020','801030','801040','801050','801080','801110','801120','801130','801140','801150',\
                    '801160','801170','801180','801200','801210','801230','801710','801720','801730','801740','801750',\
                   '801760','801770','801780','801790','801880','801890']
    
    #时段内全部因子值及当期收益率
    df = pd.DataFrame(index = get_stock(universe, period[-1]))

    for i in range(1, len(period) - 1):
        if period[i] < datetime.datetime.date(datetime.datetime.strptime('2014-02-21', '%Y-%m-%d')):
            industry_code = industry_old_code
        else:
            industry_code = industry_new_code
        all_stock = get_stock(universe, period[i])
        #获取数据
        data = get_factor_data(all_stock, period[i - 1], period[i])

        #处理数据
        p_data = data_preprocessing(data, all_stock, industry_code, period[i], ind_neutral = False)
        p_data.columns = [str(period[i]) + factor]
        df = pd.concat([df, p_data], axis = 1)
        price = get_price(p_data.index.tolist(), start_date = period[i], end_date = period[i + 1], fields = ['close'])['close']
        
        #无法相对于中证全指的超额收益，因为2005年还没有
        ret = pd.DataFrame((price.iloc[-1, :] / price.iloc[0, :]) - 1, columns = [str(period[i]) + '收益率'])
        ret.fillna(nanmedian(ret), inplace = True)
        df = pd.concat([df, ret], axis = 1)

    df.dropna(how = 'all', inplace = True)
    
    
    return df

#计算IC
def get_IC(df):
    
    #IC值序列
    IC = []
    S = []

    for i in range(1, len(period) - 1):
        raw = pd.DataFrame(index = df.index)
        raw['0'] = df[str(period[i]) + factor]
        raw['1'] = df[str(period[i]) + '收益率']
        raw.dropna(inplace = True)
        ic, p = stats.spearmanr(raw, axis = 0, nan_policy = 'omit')
        IC.append(ic)
        #显著性检验
        if p < 0.05:
            if ic > 0:
                S.append('pos')
            elif ic < 0:
                S.append('neg')
            else:
                S.append('non')
        else:
            S.append('non')
            
        
    return IC, S

#汇总结果
def draw_conclusion(period):    
    df = get_df(period)
    #f, t_f = get_f_t(df)
    IC, S = get_IC(df)
    conclusion = pd.DataFrame(index = [factor])
    conclusion['IC均值'] = nanmean(IC)
    conclusion['IC正显著比例'] = S.count('pos') / len(S)
    conclusion['IC负显著比例'] = S.count('neg') / len(S)
    conclusion['abs(IC)>0.02比例'] = len([i for i in IC if abs(i) > 0.02]) / len(IC)
    conclusion['IR'] = nanmean(IC) / std(IC) 
    print(conclusion.T)
    IC = pd.DataFrame(IC, index = period[1: len(period) - 1], columns = ['IC'])
    IC['index'] = range(len(IC))
    IC['S'] = S
    IC['pos'] = IC.loc[IC['S'] == 'pos', 'IC']
    IC['neg'] = IC.loc[IC['S'] == 'neg', 'IC']
    IC['non'] = IC.loc[IC['S'] == 'non', 'IC']
    fig1 = plt.figure(2, figsize=(30,5))
    plt.scatter(IC['index'], IC['pos'], c = 'r')
    plt.scatter(IC['index'], IC['neg'], c = 'b')
    plt.scatter(IC['index'], IC['non'], c = 'grey')
    plt.grid()
    plt.tick_params(labelsize = 20)
    plt.rcParams['axes.unicode_minus'] = False 
    #画多空组合收益率图
    r = get_Return_per_day(period, df)
    return_line(r, period)
    
    return None

#多空一十组合
def get_Return_per_day(period, df):

    #IC = pd.DataFrame()
    r = []
    for i in range(1, len(period) - 1):
        
        abin_sorted = df[str(period[i]) + factor].dropna().sort_values()
        #分组，取第一组和第十组
        n = int(len(abin_sorted)/10)
        long = abin_sorted.index.tolist()[0:n]
        short = abin_sorted.index.tolist()[-n:]

        #多头
        long_close = df[str(period[i]) + '收益率'].loc[long]
        #空头
        short_close = df[str(period[i]) + '收益率'].loc[short]

        #计算收益率
        r.append(long_close.mean() - short_close.mean())

    return r

#画多空组合收益率图
def return_line(r, period):
    r = np.array(r) + 1
    p = []
    m = 1
    p.append(m)
    for i in range(len(r)):
        m *= r[i]
        p.append(m)
        
    d = period[0: len(p)]
    p = pd.Series(p, index = d)
    fig1 = plt.figure(1, figsize=(30,5))
    plt.plot(p.index, p)
    plt.grid()
    plt.tick_params(labelsize = 20)
    plt.rcParams['axes.unicode_minus'] = False 
    
    return None

#控制获取数据时段（因为数据库只到2005年1月1日，且IC计算要前一个月的数据，因此从1月底开始计算）
start_date = '2005-01-31'
end_date = '2019-02-27'
period = get_tradeday_list(start = start_date, end = end_date, frequency = 'month').date


#所需股票池，可选'HS300'，'ZZ500'，'ZZ800'，'A'
global universe
universe = 'A'

#所需因子IVol
global factor
factor = 'IVol'
draw_conclusion(period)

#所需因子IVR
global factor
factor = 'IVR'
draw_conclusion(period)

#所需因子adjTurnover
global factor
factor = 'adjTurnover'
draw_conclusion(period)

#控制获取数据时段（因为数据库只到2005年1月1日，且IC计算要前一个月的数据，因此从1月底开始计算）
start_date = '2005-04-30'
end_date = '2019-02-27'
period = get_tradeday_list(start = start_date, end = end_date, frequency = 'month').date


#所需股票池，可选'HS300'，'ZZ500'，'ZZ800'，'A'
global universe
universe = 'A'

分组回测用自己的代码，数据处理使回测与研究模块保持一致

#所需因子PriceDelay
global factor
factor = 'PriceDelay'
draw_conclusion(period)

#回测用代码
code = '''
# 导入函数库
from jqdata import *
from jqfactor import *
import datetime as dt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels import regression


# 初始化函数，设定基准等等
def initialize(context):
    # 设定中证全指作为基准
    set_benchmark({'000002.XSHG': 0.5, '399107.XSHE': 0.5})
    # 开启动态复权模式(真实价格)
    set_option('use_real_price', True)
    log.set_level('order', 'error')
    
    #第几组
    g.group = 1
    g.method = 'adjTurnover'

    ## 运行函数（reference_security为运行时间的参考标的；传入的标的只做种类区分，因此传入'000300.XSHG'或'510300.XSHG'是一样的）
      # 开盘前运行
    run_monthly(before_market_open, monthday = -1, time='before_open', reference_security = '399107.XSHE')
      # 开盘时运行
    run_monthly(market_open, monthday = -1, time='open', reference_security = '399107.XSHE')


## 开盘前运行函数
def before_market_open(context):
    #设置排序方式
    if g.method == 'BP':
        g.asc = False
    else:
        g.asc = True
    #设置滑点、手续费
    set_slip_fee(context)
    #取全A作为股票池
    all_stocks = list(get_all_securities(['stock'], date = context.current_dt).index)
    feasible_stocks = set_feasible_stocks(context, all_stocks)

    if g.method == 'IVol' or g.method == 'IVR':
        factor = FF3(feasible_stocks, context.current_dt)
    elif g.method == 'PriceDelay':
        factor = PD(feasible_stocks, context.current_dt)
    elif g.method == 'adjTurnover':
        factor = adjTurnover(feasible_stocks, context.current_dt)
          
    
    #排序
    factor = factor.sort(g.method, ascending = g.asc)

    n = int(len(factor)/10)
    #分组取样
    if g.group == 10:
        g.tobuy_list = factor.index[(g.group - 1) * n :]
    else:
        g.tobuy_list = factor.index[(g.group - 1) * n : g.group * n]
    
    
#1
#设置可行股票池，剔除st、停牌股票，输入日期
def set_feasible_stocks(context, stockList):
    #剔除ST股
    st_data = get_extras('is_st', stockList, count = 1, end_date = context.current_dt)
    stockList = [stock for stock in stockList if not st_data[stock][0]]
    #剔除*st股票
    stockList = [stock for stock in stockList if '*' not in get_security_info(stock).display_name]
    #剔除上市不足6月的新股
    stockList = delete_new(stockList, context.current_dt, n = 183)
    #剔除停牌
    suspended_info_df = get_price(stockList, end_date = context.current_dt, count = 1, frequency = 'daily', fields = 'paused')['paused']
    stockList = [stock for stock in stockList if suspended_info_df[stock][0] == 0]
    
    return stockList   

#剔除新股
def delete_new(stocks, beginDate, n = 365):
    stockList = []
    for stock in stocks:
        start_date = get_security_info(stock).start_date
        if start_date < dt.datetime.date(beginDate - dt.timedelta(days = n)):
            stockList.append(stock)
    return stockList
    

# 根据不同的时间段设置滑点与手续费
def set_slip_fee(context):
    # 将滑点设置为0
    set_slippage(FixedSlippage(0)) 
    # 根据不同的时间段设置手续费
    dt=context.current_dt
    
    if dt>datetime.datetime(2013,1, 1):
        set_commission(PerTrade(buy_cost=0.0003, 
                                sell_cost=0.0013, 
                                min_cost=5)) 
        
    elif dt>datetime.datetime(2011,1, 1):
        set_commission(PerTrade(buy_cost=0.001, 
                                sell_cost=0.002, 
                                min_cost=5))
            
    elif dt>datetime.datetime(2009,1, 1):
        set_commission(PerTrade(buy_cost=0.002, 
                                sell_cost=0.003, 
                                min_cost=5))
                
    else:
        set_commission(PerTrade(buy_cost=0.003, 
                                sell_cost=0.004, 
                                min_cost=5))


## 开盘时运行函数
def market_open(context):
    #调仓，先卖出股票
    for stock in context.portfolio.long_positions:
        if stock not in g.tobuy_list:
            order_target_value(stock, 0)
    #再买入新股票
    total_value = context.portfolio.total_value # 获取总资产
    for i in range(len(g.tobuy_list)):
        value = total_value / len(g.tobuy_list) # 确定每个标的的权重

        order_target_value(g.tobuy_list[i], value) # 调整标的至目标权重
    
    #查看本期持仓股数
    print(len(context.portfolio.long_positions))
    

#求Fama-French三因子模型特质波动率
def FF3(stocks, end_date, rf = 0.04):
    LoS=len(stocks)
    #查询三因子/五因子的语句
    q = query(
        valuation.code,
        valuation.circulating_market_cap,
        (balance.total_owner_equities/valuation.circulating_market_cap/100000000.0).label("BP"),
        #indicator.roe,
        #balance.total_assets.label("Inv")
    ).filter(
        valuation.code.in_(stocks)
    )
    
    start_date = list(get_tradeday_list(start = None, end = end_date, frequency = 'month', count = 24).date)[0]
    df = get_fundamentals(q, start_date)
    df.index = df['code']
    del df['code']
    #中性化
    #df = neutralize(df, how = ['sw_l1', 'market_cap'], date = start_date, axis = 0)

    #选出特征股票组合
    S=df.sort('circulating_market_cap').index.tolist()[:LoS/3]
    B=df.sort('circulating_market_cap').index.tolist()[LoS-LoS/3:]
    L=df.sort('BP').index.tolist()[:LoS/3]
    H=df.sort('BP').index.tolist()[LoS-LoS/3:]
    
    df5 = df['circulating_market_cap']
    
    # 获得样本期间的股票价格并计算日收益率
    df2 = get_price(stocks, start_date = start_date, end_date = end_date, fields=['close'])['close']
    df4 = df2.pct_change()
    df4.dropna(how ='all', inplace = True)
    df4.fillna(0, inplace = True)
    #求因子的值，按流通市值加权
    SMB = list(np.dot(df4[S], df5.loc[S] / df5.loc[S].sum()) - np.dot(df4[B], df5.loc[B] / df5.loc[B].sum()))
    HML = list(np.dot(df4[H], df5.loc[H] / df5.loc[H].sum()) - np.dot(df4[L], df5.loc[L] / df5.loc[L].sum()))
    
    #用股票池，流通市值为权重作为市场基准
    df6 = df5.loc[df4.columns]
    df6.fillna(df5.mean(), inplace = True)
    RM = list(np.dot(df4, df6 / df6.sum()) - rf/252)
    if len(SMB) > len(RM):
        SMB.drop(SMB.index[0], inplace = True)
        HML.drop(HML.index[0], inplace = True)

    #将因子们计算好并且放好
    X=pd.DataFrame({"RM":RM, "SMB":SMB, "HML":HML})
    # 对样本数据进行线性回归并计算ai
    t_scores=[0.0] * LoS
    for i in range(LoS):
        t_stock = stocks[i]
        t_r = linreg(X, df4[t_stock] - rf/252)
        t_scores[i] = t_r
    
    scores = pd.DataFrame({g.method: t_scores}, index = stocks)
    
    return scores


#9
# 辅助线性回归的函数
# 输入:X:回归自变量 Y:回归因变量 完美支持list,array,DataFrame等三种数据类型
# 输出:参数估计结果-list类型
def linreg(X,Y):
    X=sm.add_constant(array(X))
    Y=array(Y)
    results = sm.OLS(Y, X).fit()
    
    if g.method == 'IVol':
        return results.resid.std() * sqrt(252)
    elif g.method == 'IVR':
        return (1 - results.rsquared)
    elif g.method == 'PriceDelay':
        return results.rsquared
    elif g.method == 'adjTurnover':
        return results.resid

    
#计算价格时滞
def PD(stocks, end_date, rf = 0.04):
    
    minus1_date = list(get_tradeday_list(start = None, end = end_date, frequency = 'month', count = 24).date)[0]
    X1 = MKT(stocks, minus1_date, end_date, n = 0)
    X2 = MKT(stocks, minus1_date, end_date, n = 3)
    quote = get_price(stocks, start_date = minus1_date, end_date = end_date, fields=['close'])['close']
    ret = quote.pct_change()
    ret.dropna(how ='all', inplace = True)
    ret.fillna(0, inplace = True)
    ret = ret - rf/252
    #OLS计算hetero
    PriceDelay = pd.DataFrame(index = [g.method])
    for stock in ret.columns:
        Y1 = ret[stock][:len(X1)]
        Y2 = ret[stock][:len(X2)]
        r0 = linreg(X1 - rf/252, Y1)
        r3 = linreg(X2 - rf/252, Y2)
        PriceDelay[stock] = 1 - r0/r3
    
    PriceDelay = PriceDelay.T

    #返回价格时滞
    return PriceDelay

#计算滞后n期的市场因子
def MKT(stocks, minus1_date, end_date, n):
    if n == 0:
        period = [minus1_date, end_date]
    else:
        period = list(get_tradeday_list(start = None, end = minus1_date, frequency = 'month', count = n * 24).date)
        period.append(end_date)
    
    X = pd.DataFrame(columns = period[1:], index = range(24))
    for i in range(len(period) - 1):
        #设置统计范围
        quote = get_price(stocks, start_date = period[i], end_date = period[i + 1], fields=['close'])['close']
        ret = quote.pct_change()
        ret.dropna(how ='all', inplace = True)
        ret.fillna(0, inplace = True)
        #构造市场基准收益：流通市值加权
        q = query(valuation.circulating_market_cap, valuation.code).filter(valuation.code.in_(stocks))
        df = get_fundamentals(q, period[i])
        df.index = df['code']
        del df['code']
        #df = (df * 1e8).apply(np.log)
        df = df/df.sum()
        
        final = pd.DataFrame(np.dot(ret, df), columns = ['final'])
        X[period[i + 1]] = final['final']
        
    X.dropna(how = 'any', inplace = True)
    
    return X
    
def get_tradeday_list(start,end,frequency=None,count=None):
    if count != None:
        df = get_price('000001.XSHG',end_date=end,count=count)
    else:
        df = get_price('000001.XSHG',start_date=start,end_date=end)
    if frequency == None or frequency =='day':
        return df.index
    else:
        df['year-month'] = [str(i)[0:7] for i in df.index]
        if frequency == 'month':
            return df.drop_duplicates('year-month', take_last = True).index
        elif frequency == 'quarter':
            df['month'] = [str(i)[5:7] for i in df.index]
            df = df[(df['month']=='01') | (df['month']=='04') | (df['month']=='07') | (df['month']=='10') ]
            return df.drop_duplicates('year-month').index
        elif frequency =='halfyear':
            df['month'] = [str(i)[5:7] for i in df.index]
            df = df[(df['month']=='01') | (df['month']=='07')]
            return df.drop_duplicates('year-month').index 
            
            
#流通市值调整换手
def adjTurnover(stocks, end_date):
    start_date = list(get_tradeday_list(start = None, end = end_date, frequency = 'month', count = 24).date)[0]
    VOL20 = get_factor_values(stocks, ['VOL20'], end_date = end_date, count = 1)['VOL20'].T
    VOL20.replace(0, np.nan, inplace = True)
    VOL20 = VOL20.apply(np.log)
    VOL20.columns = ['VOL20']
    q = query(valuation.code, valuation.circulating_market_cap).filter(valuation.code.in_(stocks))
    df = get_fundamentals(q, end_date)
    df.index = df['code']
    del df['code']
    VOL20['circulating_market_cap'] = (df['circulating_market_cap'] * 1e8).apply(np.log)
    VOL20.dropna(how = 'any', inplace = True)

    X = VOL20['circulating_market_cap']
    Y = VOL20['VOL20']
    l = linreg(X,Y)

    adjTurnover = pd.DataFrame(l, index = VOL20.index, columns = ['adjTurnover'])
    
    return adjTurnover
    
## 收盘后运行函数
def after_market_close(context):
    pass
    
'''

下个模块里method默认为'IVR'，需要其他因子的时候更改这个变量至'IVol', 'PriceDelay', 'adjTurnover'即可。按顺序一次一个，重复运行下面4个模块

#研究调用回测
created_bt_ids = []
method = 'IVR'
#取全部组别进行回测
for i in range(1, 11):
    algorithm_id = "e1f0bf5d8c9d75c38497c9b00a01eb9d" #用自己的策略ID
    extra_vars = {'group': i, 'method': method}

    params = {
                "algorithm_id": algorithm_id,
                "start_date": "2005-05-31",
                "end_date": "2019-02-27",
                "frequency": "day",
                "initial_cash": "30000000",
                "initial_positions": None,
                "extras": extra_vars, 
                "name" : method + ' group:' + str(i) 
                }
    created_bt_ids.append(create_backtest(code = code, **params))

#等上个模块的回测全部运行结束再运行这个模块

#先获取基准收益和交易日信息
gt = get_backtest(backtest_id = created_bt_ids[0])
res = gt.get_results()
b_return = []
t = []
for r in res:
    b_return.append(r['benchmark_returns'])
    t.append(r['time'])

#建立df存储数据   
data = pd.DataFrame(index = t)
data['b_return'] = b_return

#填入不同参数下的收益数据
for i in range(len(created_bt_ids)):
    gt = get_backtest(backtest_id = created_bt_ids[i])
    res = gt.get_results()
    name = gt.get_params()['name'] 
    s_return = []
    for r in res:
        s_return.append(r['returns'])
    
    data[name] = s_return

#年化超额收益图
dat = (data.iloc[-1, :] + 1)
((dat / dat[0]) ** (1/14) - 1).drop('b_return').plot.bar(figsize = (7.5, 5))

#夏普比率图
sharpe = []
max_drawdown = []
information_ratio = []
name = range(1, 11)
for i in range(len(created_bt_ids)):
    gt = get_backtest(backtest_id = created_bt_ids[i])
    sharpe.append(gt.get_risk()['sharpe'])
    max_drawdown.append(gt.get_risk()['max_drawdown']) 
    information_ratio.append(gt.get_risk()['information'])

sharpe = pd.Series(sharpe, index = name)
max_drawdown = pd.Series(max_drawdown, index = name)
information_ratio = pd.Series(information_ratio, index = name)
sharpe.plot.bar(figsize = (7.5, 5))

因子收益率回归解释

#获取全部期因子收益率
def get_RegExp_df(period):
    
    IVol = []
    IVR= [] 
    PriceDelay = []
    adjT = [] 
    MKT = [] 
    SMB = [] 
    HML = [] 
    MOM = []
    for i in range(1, len(period) - 1):
        stocks = get_stock(universe, period[i - 1])
        LoS = len(stocks)
        #查询因子的语句
        q = query(
            valuation.code,
            valuation.circulating_market_cap,
            (balance.total_owner_equities/valuation.circulating_market_cap/100000000.0).label("BP"),
        ).filter(
            valuation.code.in_(stocks)
        )

        df = get_fundamentals(q, period[i - 1])
        df.index = df['code']
        del df['code']
        #中性化
        #df = neutralize(df, how = ['sw_l1', 'market_cap'], date = end, axis = 0)
        n = int(LoS/3)

        #选出特征股票组合
        S=df.sort_values('circulating_market_cap').index.tolist()[:n]
        B=df.sort_values('circulating_market_cap').index.tolist()[LoS - n:]
        L=df.sort_values('BP').index.tolist()[:n]
        H=df.sort_values('BP').index.tolist()[LoS - n:]

        #df5 = (df['circulating_market_cap'] * 1e8).apply(np.log)
        df5 = df['circulating_market_cap']
        # 获得样本期间的股票价格并计算月收益率
        df2 = get_price(stocks, start_date = period[i - 1], end_date = period[i], fields=['close'])['close']
        df4 = df2.iloc[-1, :] / df2.iloc[0, :] - 1
        df4.fillna(0, inplace = True)

        #动量（反转）因子
        prior_date = list(get_tradeday_list(start = None, end = period[i - 1], frequency = 'month', count = 24).date)[0]
        df22 = get_price(stocks, start_date = prior_date, end_date = period[i - 1], fields=['close'])['close']
        df42 = pd.DataFrame(df22.iloc[-1, :] / df22.iloc[0, :] - 1, columns = ['ret'])
        df42.fillna(0, inplace = True)
        LO = df42.sort_values('ret').index.tolist()[: n]
        W = df42.sort_values('ret').index.tolist()[LoS - n:]

        #求因子的值，按流通市值加权
        SMB.append(np.dot(df4[S], df5.loc[S] / df5.loc[S].sum()) - np.dot(df4[B], df5.loc[B] / df5.loc[B].sum()))
        HML.append(np.dot(df4[H], df5.loc[H] / df5.loc[H].sum()) - np.dot(df4[L], df5.loc[L] / df5.loc[L].sum()))

        #价格数据股票比财务数据多
        df6 = df5.loc[df4.index]
        df6.fillna(df5.mean(), inplace = True)
        MOM.append(np.dot(df4[LO], df6.loc[LO] / df6.loc[LO].sum()) - np.dot(df4[W], df6.loc[W] / df6.loc[W].sum()))
        #用股票池，流通市值为权重作为市场基准
        MKT.append(np.dot(df4, df6 / df6.sum()))
        
        global factor
        factor = 'IVol'
        factor_data = FF3(stocks, period[i - 1], period[i])
        factor_data.columns = [factor]
        first = factor_data.sort_values(factor).index.tolist()[:n]
        last = factor_data.sort_values(factor).index.tolist()[LoS - n:]
        IVol.append(np.dot(df4[first], df6.loc[first] / df6.loc[first].sum()) - np.dot(df4[last], df6.loc[last] / df6.loc[last].sum()))
        factor = 'IVR'
        factor_data = FF3(stocks, period[i - 1], period[i])
        factor_data.columns = [factor]
        first = factor_data.sort_values(factor).index.tolist()[:n]
        last = factor_data.sort_values(factor).index.tolist()[LoS - n:]
        IVR.append(np.dot(df4[first], df6.loc[first] / df6.loc[first].sum()) - np.dot(df4[last], df6.loc[last] / df6.loc[last].sum()))
        factor = 'PriceDelay'
        factor_data = PD(stocks, period[i - 1], period[i])
        factor_data.columns = [factor]
        first = factor_data.sort_values(factor).index.tolist()[:n]
        last = factor_data.sort_values(factor).index.tolist()[LoS - n:]
        PriceDelay.append(np.dot(df4[first], df6.loc[first] / df6.loc[first].sum()) - np.dot(df4[last], df6.loc[last] / df6.loc[last].sum()))
        factor = 'adjTurnover'
        factor_data = adjTurnover(stocks, period[i])
        factor_data.columns = [factor]
        first = factor_data.sort_values(factor).index.tolist()[:n]
        last = factor_data.sort_values(factor).index.tolist()[LoS - n:]
        adjT.append(np.dot(df4[first], df6.loc[first] / df6.loc[first].sum()) - np.dot(df4[last], df6.loc[last] / df6.loc[last].sum()))
    
    
    #将因子们计算好并且放好
    f_data = pd.DataFrame({"MKT": MKT, "SMB": SMB, "HML": HML, "MOM": MOM,
                                'IVol': IVol, 'IVR': IVR, 'PriceDelay': PriceDelay, 'adjTurnover': adjT})
    f_data.dropna(how = 'any', inplace = True)
    
    return f_data

#控制获取数据时段（因为数据库只到2005年1月1日，且IC计算要前一个月的数据，因此从1月底开始计算）
start_date = '2005-04-30'
end_date = '2019-02-27'
period = get_tradeday_list(start = start_date, end = end_date, frequency = 'month').date

RegExp_df = get_RegExp_df(period)

#线性回归解释，Yf是被解释变量（只能有一个），Xf是用来解释的变量列表（一个或多个，不能包含Yf）
#可选变量列表['IVol', 'IVR', 'PriceDelay', 'adjTurnover', 'MKT', 'SMB', 'HML', 'MOM']
def RegExp(df, Yf, Xf):
    
    X = sm.add_constant(array(df[Xf], dtype=float))
    Y = array(df[Yf])
    results = sm.OLS(Y, X).fit()
    print(results.summary())
    
    return None

RegExp(RegExp_df, Yf = 'IVol', Xf = ['IVR', 'adjTurnover'])

RegExp(RegExp_df, Yf = 'IVR', Xf = ['IVol', 'PriceDelay', 'adjTurnover'])

RegExp(RegExp_df, Yf = 'IVR', Xf = ['MKT', 'SMB', 'HML', 'MOM', 'IVol', 'PriceDelay', 'adjTurnover'])

RegExp(RegExp_df, Yf = 'PriceDelay', Xf = ['IVR'])

RegExp(RegExp_df, Yf = 'adjTurnover', Xf = ['IVol', 'IVR', 'PriceDelay'])

RegExp(RegExp_df, Yf = 'adjTurnover', Xf = ['MKT', 'SMB', 'HML', 'MOM', 'IVol', 'IVR', 'PriceDelay'])

量化交易吧 / 数理科学 帖子：3354223 新帖：56

复现东方证券研报--投机、交易行为与股票收益

专门套利发表于：5 月 10 日 03：02回复(1)

结论概述

被投机个股的交易行为特征

交易行为类指标

1、特质波动率

2、特异度

3、价格时滞

4、流通市值调整换手率

5、小结

因子收益率回归检查替代关系

风险提示

全部回复

0/140

粉丝:913

帖子数:0

粉丝:733

帖子数:0

粉丝:675

帖子数:391

量化课程

热门标签

删除回复

确认要删除这篇文章么？

举报用户

信息提示

该文章已删除

设置置顶

完成设置【置顶】！

设置置顶

已取消设置【置顶】！

设置精华

完成设置【精华】！

设置精华

已取消设置【精华】！

审核信息

该文章已审核通过

审核信息

您已设置该文章审核不通过

举报成功

您已举报成功

用户登录

移动帖子

创建私信

屏蔽提示

确认要屏蔽该用户么？

屏蔽回复

您已对该用户实现屏蔽

信息回复

已发送成功

量化交易吧 / 数理科学帖子：3354223 新帖：56