#数据获取
from jqdata import *
import time,datetime
import numpy as np
import pandas as pd
from jqdata import finance

# #获取每周最后一个交易日的申万一级行业指数收盘价
industries_df = get_industries("sw_l1")
deprecated_industries = ['801060','801070','801090','801100','801190','801220']
industries = list(industries_df.index)
start_dates = list(industries_df['start_date'])
trade_days = get_trade_days(start_date='2005-02-03')
ret_dict={}
for code,start_date in zip(industries,start_dates):
    if code in deprecated_industries:
        continue
    #国证A指在2005年2月3日后才有数据
    start_date = trade_days[0]
    q=query(finance.SW1_DAILY_PRICE).filter(finance.SW1_DAILY_PRICE.date>=start_date,finance.SW1_DAILY_PRICE.code==code)
    df=finance.run_query(q)
    if len(df)!=len(trade_days):
        trade_days=trade_days[0:-1]
    df_date=list(df['date'])
    df_date_weekday =np.array([i.weekday() for i in df_date])
    b=df_date_weekday[0:-1]-df_date_weekday[1:]
    c=df['date']
    close=np.array(df[0:-1][b!=-1]['close'])
    ret = close[1:]/close[:-1]-1
    ret_dict[code]=ret
    
ret_dict['date'] = df[0:-1][b!=-1]['date'][1:]

#获取每周最后一个交易日的国证A指
mkt_index = get_price('399317.XSHE',start_date='2005-02-03',end_date=trade_days[-1],fields=['close'])
mkt_close = np.array(mkt_index[0:-1][b!=-1]['close'])
mkt_ret =mkt_close[1:]/mkt_close[0:-1]-1
ret_dict['mkt'] = mkt_ret
ret_df = pd.DataFrame(ret_dict)

ret_df.head(10)

import statsmodels.api as sm
#用N周的数据计算Beta,第N周开始有beta，研报中N为98

#获取行业中文名称
ind_df=get_industries('sw_l1')
ind = ind_df.index
valid_ind = list(set(ind)^set(deprecated_industries))
valid_ind_df = ind_df.loc[valid_ind]
ind_dict = dict(zip(valid_ind_df.index,valid_ind_df['name']))

#计算行业beta
N=98
beta_dict = {}
beta_dict['date'] = ret_dict['date'][N-1:]
for k,ret in ret_dict.items():
    if k!='date'and k!='mkt':
        beta_series = []
        for i in range(len(ret)-N+1):
            x = sm.add_constant(mkt_ret[i:i+N])
            model = sm.OLS(ret[i:i+N],x)
            results = model.fit()
            beta_series.append(results.params[1])
        beta_dict[k] = beta_series

beta_df = pd.DataFrame(beta_dict)        
date = beta_df['date']
beta_df_t = beta_df.T
beta_rank=beta_df_t.iloc[1:].rank().T
plt.figure()
for i in range(1,29):
    print(i)
    #plt.subplot(28,1,i)
    plt.figure(figsize=(10,8))
    code = beta_df.columns[i]
    x = beta_dict['date']
    y1 = beta_dict[code]
    y2 = beta_rank[code]
    fig, ax1 = plt.subplots()
    ax2 = ax1.twinx()
    ax1.plot(x,y1,'g-')
    ax2.plot(x,y2,'b-')
    ax1.set_xlabel(ind_dict[code])
    ax1.set_ylabel("Beta",color='g')

    ax2.set_ylabel("Rank",color='b')

    plt.show()

from scipy.stats import spearmanr
threshold=0.128
   
#threshold为阈值，print_order为是否打印买卖指令
def back_test(threshold,print_order=True):
    corr_list =[]
    #第t周收益与t-1周的beta
    for i in range(len(beta_df)-1):
        corr_list.append(spearmanr(beta_df.iloc[i][1:].astype(float),ret_df.iloc[i+N][:-2].astype(float))[0])
    corr_dict={}
    corr_dict['date']=beta_df['date'][1:]
    corr_dict['corr']=corr_list
    corr_df=pd.DataFrame(corr_dict)
    corr_df['ma_4'] = corr_df['corr'].rolling(window=4,min_periods=1).mean()
    corr_df['signal']=0
    corr_df['signal'][corr_df['ma_4']>threshold]=1
    corr_df['signal'][corr_df['ma_4']<-threshold]=-1

    position_list = [0]
    flag=0
    order_dict = {'date':[],'指令':[]}
    for i in range(1,len(corr_df['signal'])):
        if flag==0:
            if np.array(corr_df['signal'])[i-1]==1 and np.array(corr_df['signal'])[i]==1:
                position_list.append(1)
                flag=1
                order_dict['date'].append(corr_df['date'].iloc[i])
                order_dict['指令'].append('买入')
            else:
                position_list.append(0)
        elif flag==1:
            if np.array(corr_df['signal'])[i-1]==-1 and np.array(corr_df['signal'])[i]==-1:
                position_list.append(0)
                flag=0
                order_dict['date'].append(corr_df['date'].iloc[i])
                order_dict['指令'].append('卖出')
            else:
                position_list.append(1)
    order_df = pd.DataFrame(order_dict)
    if print_order:
        print(order_df)
    corr_df['position']=position_list
    corr_df['ret']=0
    corr_df['ret'].iloc[1:] = np.array(ret_df['mkt'])[N+1:]*np.array(corr_df['position'])[:-1]
    ret = list(corr_df['ret'])
    cum_ret=[]
    for i in range(len(ret)):
        if i==0:
            cum_ret.append(1+ret[i])
        else:
            cum_ret.append((1+ret[i])*cum_ret[i-1])
    corr_df['cum_ret']=cum_ret
    return corr_df

corr_df=back_test(threshold)

#计算各项指标,一年按50周计算
max_cum_ret=[]
cum_ret = np.array(corr_df['cum_ret'])
for i in range(len(cum_ret)):
    max_cum_ret.append(max(cum_ret[:i+1]))
drawdown=1-corr_df['cum_ret']/max_cum_ret

annual_ret = corr_df['cum_ret'].iloc[-1]**(1*50/len(corr_df))-1
vol = corr_df['ret'].std()*np.sqrt(50)
mdd = max(drawdown)
sharpe = annual_ret/vol
summary_df = pd.DataFrame({'年化收益':[annual_ret],'年化波动率':[vol],'最大回撤':[mdd],'夏普比率':[sharpe]})
print(summary_df)

#作图
plt.figure(figsize=(15,6))
mkt_index = get_price('399317.XSHE',start_date='2005-02-03',end_date=trade_days[-1],fields=['close'])
data=corr_df.set_index('date')
combination=data.join(mkt_index)
plt.plot(combination.index,combination['cum_ret'])
plt.plot(combination.index,combination['close']/combination['close'].iloc[0])
plt.show()

thresholds=np.array(range(100,301,1))/1000
sharpe_dict={'阈值':[],'夏普比率':[]}
for threshold in thresholds:
    corr_df=back_test(threshold,print_order=False)
    annual_ret = corr_df['cum_ret'].iloc[-1]**(1*50/len(corr_df))-1
    vol = corr_df['ret'].std()*np.sqrt(50)
    sharpe = annual_ret/vol
    sharpe_dict['阈值'].append(threshold)
    sharpe_dict['夏普比率'].append(sharpe)

sharpe_df = pd.DataFrame(sharpe_dict)
plt.figure(figsize=(20,10))
plt.plot(sharpe_df["阈值"],sharpe_df['夏普比率'])
plt.show()

corr_df=back_test(0.142)

#计算各项指标,一年按50周计算
max_cum_ret=[]
cum_ret = np.array(corr_df['cum_ret'])
for i in range(len(cum_ret)):
    max_cum_ret.append(max(cum_ret[:i+1]))
drawdown=1-corr_df['cum_ret']/max_cum_ret

annual_ret = corr_df['cum_ret'].iloc[-1]**(1*50/len(corr_df))-1
vol = corr_df['ret'].std()*np.sqrt(50)
mdd = max(drawdown)
sharpe = annual_ret/vol
summary_df = pd.DataFrame({'年化收益':[annual_ret],'年化波动率':[vol],'最大回撤':[mdd],'夏普比率':[sharpe]})
print(summary_df)

#作图
plt.figure(figsize=(15,6))
data=corr_df.set_index('date')
combination=data.join(mkt_index)
plt.plot(combination.index,combination['cum_ret'])
plt.plot(combination.index,combination['close']/combination['close'].iloc[0])
plt.show()
corr_df.to_excel('0.142.xlsx')

def combination_back_test(threshold,print_order=True):
#生成均线信号
    corr_list =[]
    #第t周收益与t-1周的beta
    for i in range(len(beta_df)-1):
        corr_list.append(spearmanr(beta_df.iloc[i][1:].astype(float),ret_df.iloc[i+N][:-2].astype(float))[0])
    corr_dict={}
    corr_dict['date']=beta_df['date'][1:]
    corr_dict['corr']=corr_list
    corr_df=pd.DataFrame(corr_dict)
    corr_df['ma_4'] = corr_df['corr'].rolling(window=4,min_periods=1).mean()
    corr_df['signal']=0
    corr_df['signal'][corr_df['ma_4']>threshold]=1
    corr_df['signal'][corr_df['ma_4']<-threshold]=-1
    mkt_index = get_price('399317.XSHE',start_date='2005-02-03',end_date=trade_days[-1],fields=['close'])
    mkt_index['ma_20'] = mkt_index.rolling(window=20).mean()
    data=corr_df.set_index('date',drop=False)
    combination=data.join(mkt_index)
    combination['ma_signal']=(combination['close']>combination['ma_20']).astype(int)
    
    position_list = [0]
    flag=0
    order_dict = {'date':[],'指令':[]}
    for i in range(1,len(corr_df['signal'])):
        if flag==0:
            if np.array(combination['signal'])[i-1]==1 and np.array(combination['signal'])[i]==1:
                if np.array(combination['ma_signal'])[i]==1:
                    position_list.append(1)
                    flag=1
                    order_dict['date'].append(combination['date'].iloc[i])
                    order_dict['指令'].append('买入')
                else:
                    position_list.append(0)
            else:
                position_list.append(0)
        elif flag==1:
            if np.array(combination['signal'])[i-1]==-1 and np.array(combination['signal'])[i]==-1:
                if np.array(combination['ma_signal'])[i]==0:
                    position_list.append(0)
                    flag=0
                    order_dict['date'].append(combination['date'].iloc[i])
                    order_dict['指令'].append('卖出')
                else:
                    position_list.append(1)
            else:
                position_list.append(1)
    if print_order:
        order_df = pd.DataFrame(order_dict)
        print(order_df)
    combination['position']=position_list
    combination['ret']=0
    combination['ret'].iloc[1:] = np.array(ret_df['mkt'])[N+1:]*np.array(combination['position'])[:-1]
    ret = list(combination['ret'])
    cum_ret=[]
    for i in range(len(ret)):
        if i==0:
            cum_ret.append(1+ret[i])
        else:
            cum_ret.append((1+ret[i])*cum_ret[i-1])
    combination['cum_ret']=cum_ret
    return combination

combination = combination_back_test(0.142)
#计算各项指标,一年按50周计算
max_cum_ret=[]
cum_ret = np.array(combination['cum_ret'])
for i in range(len(cum_ret)):
    max_cum_ret.append(max(cum_ret[:i+1]))
drawdown=1-combination['cum_ret']/max_cum_ret

annual_ret = combination['cum_ret'].iloc[-1]**(1*50/len(combination))-1
vol = combination['ret'].std()*np.sqrt(50)
mdd = max(drawdown)
sharpe = annual_ret/vol
summary_df = pd.DataFrame({'年化收益':[annual_ret],'年化波动率':[vol],'最大回撤':[mdd],'夏普比率':[sharpe]})
print(summary_df)

#作图
plt.figure(figsize=(15,6))
#mkt_index = get_price('399317.XSHE',start_date='2005-02-03',end_date=trade_days[-1],fields=['close'])
#data=combination.set_index('date')
#combination=data.join(mkt_index)
plt.plot(combination.index,combination['cum_ret'])
plt.plot(combination.index,combination['close']/combination['close'].iloc[0])
plt.show()

combination

thresholds=np.array(range(100,301,1))/1000
sharpe_dict={'阈值':[],'夏普比率':[]}
for threshold in thresholds:
    corr_df=combination_back_test(threshold,print_order=False)
    annual_ret = corr_df['cum_ret'].iloc[-1]**(1*50/len(corr_df))-1
    vol = corr_df['ret'].std()*np.sqrt(50)
    sharpe = annual_ret/vol
    sharpe_dict['阈值'].append(threshold)
    sharpe_dict['夏普比率'].append(sharpe)

sharpe_df = pd.DataFrame(sharpe_dict)
plt.figure(figsize=(20,10))
plt.plot(sharpe_df["阈值"],sharpe_df['夏普比率'])
plt.show()