import numpy as np
import pandas as pd
from scipy.optimize import minimize
import datetime as dt
from six import StringIO
from dateutil.parser import parse
import cPickle as pickle
import seaborn as sns
import matplotlib as mpl

body=read_file('全球资产指数行情序列.xlsx')
prices=pd.read_excel(StringIO(body))
# print prices.columns
columns =[u'沪深300',u'中证500',u'恒生指数',u'标普500',u'纳斯达克指数',u'富时100',u'法国CAC40',u'德国DAX',u'日经225',\
        u'上证5年国债(全)',u'上证企债',u'中信沪金商品指数',u'中信沪铜商品指数',u'中信玉米商品指数',\
        u'中信郑棉商品指数', u'中信白糖商品指数',u'南华商品指数',u'布伦特原油']
prices.set_index(u'时间',inplace=True)
prices = prices[columns]['2008-12-31':]
prices.head().append(prices.tail())

mpl.rcParams['font.family']='serif'
mpl.rcParams['axes.unicode_minus']=False # 处理负号
(prices/prices.iloc[0]).plot(figsize=(15,8),grid='on')
pct_daily = prices.pct_change()

fig = plt.figure(figsize= (15,10))
ax = fig.add_subplot(111)
ax = sns.heatmap(pct_daily.corr(),annot=True,annot_kws={'size':9,'weight':'bold'})

#获取指定周期的日期列表 'W、M、Q'
def get_period_date(peroid,start_date, end_date):
    #设定转换周期period_type  转换为周是'W',月'M',季度线'Q',五分钟'5min',12天'12D'
    stock_data = get_price('000001.XSHE',start_date,end_date,'daily',fields=['close'])
    #记录每个周期中最后一个交易日
    stock_data['date']=stock_data.index
    #进行转换，周线的每个变量都等于那一周中最后一个交易日的变量值
    period_stock_data=stock_data.resample(peroid,how='last')
    date=period_stock_data.index
    pydate_array = date.to_pydatetime()
    date_only_array = np.vectorize(lambda s: s.strftime('%Y-%m-%d'))(pydate_array )
    date_only_series = pd.Series(date_only_array)
    start_date = datetime.datetime.strptime(start_date, "%Y-%m-%d")
    start_date=start_date-datetime.timedelta(days=1)
    start_date = start_date.strftime("%Y-%m-%d")
    date_list=date_only_series.values.tolist()
    date_list.insert(0,start_date)
    return date_list

def get_near_tradingday(date):
    return (get_price('000001.XSHG',end_date=date,count=1).index[0])


date_list = (get_period_date('Q','2010-01-01', '2018-10-01'))
# type(date_list[0])
tradingday_list_tot = pct_daily.index.tolist()
tradingday_list = [date for date in pct_daily.index.tolist() if ((date)>=get_near_tradingday(date_list[0]) \
                                                                 and (date)<=get_near_tradingday(date_list[-1]))]
np.array(date_list)
# tradingday_list

def get_smart_weight(pct, method='risk parity', cov_adjusted=False, wts_adjusted=False):
    if cov_adjusted == False:
        #协方差矩阵
        cov_mat = pct.cov()
    else:
        #调整后的半衰协方差矩阵
        cov_mat = pct.iloc[:len(pct)/4].cov()*(1/10.) + pct.iloc[len(pct)/4+1:len(pct)/2].cov()*(2/10.) +\
            pct.iloc[len(pct)/2+1:len(pct)/4*3].cov()*(3/10.) + pct.iloc[len(pct)/4*3+1:].cov()*(4/10.)
    if not isinstance(cov_mat, pd.DataFrame):
        raise ValueError('cov_mat should be pandas DataFrame！')
        
    omega = np.matrix(cov_mat.values)  # 协方差矩阵

    a, b = np.linalg.eig(np.array(cov_mat)) #a为特征值,b为特征向量
    a = matrix(a)
    b = matrix(b)
    # 定义目标函数
    
    def fun1(x):
        tmp = (omega * np.matrix(x).T).A1
        risk = x * tmp/ np.sqrt(np.matrix(x) * omega * np.matrix(x).T).A1[0]
        delta_risk = [sum((i - risk)**2) for i in risk]
        return sum(delta_risk)
    
    def fun2(x):
        tmp = (b**(-1) * omega * np.matrix(x).T).A1
        risk = (b**(-1)*np.matrix(x).T).A1 * tmp/ np.sqrt(np.matrix(x) * omega * np.matrix(x).T).A1[0]
        delta_risk = [sum((i - risk)**2) for i in risk]
        return sum(delta_risk)
    
    # 初始值 + 约束条件 
    x0 = np.ones(omega.shape[0]) / omega.shape[0]  
    bnds = tuple((0,None) for x in x0)
    cons = ({'type':'eq', 'fun': lambda x: sum(x) - 1})
    options={'disp':False, 'maxiter':1000, 'ftol':1e-20}
        
    if method == 'risk parity':
        res = minimize(fun1, x0, bounds=bnds, constraints=cons, method='SLSQP', options=options)
    elif method == 'pc risk parity':
        res = minimize(fun2, x0, bounds=bnds, constraints=cons, method='SLSQP', options=options)
    else:
        raise ValueError('method error！！！')
        
    # 权重调整
    if res['success'] == False:
        # print res['message']
        pass
    wts = pd.Series(index=cov_mat.index, data=res['x'])
    
    if wts_adjusted == True:
        wts[wts < 0.0001]=0.0
        wts = wts / wts.sum()
    elif wts_adjusted == False:
        wts = wts / wts.sum()
    else:
        raise ValueError('wts_adjusted should be True/False！')
        
    risk = pd.Series(wts * (omega * np.matrix(wts).T).A1 / np.sqrt(np.matrix(wts) * omega * np.matrix(wts).T).A1[0],index = cov_mat.index)
    risk[risk<0.0] = 0.0
    return wts,risk

pct_daily_new = pct_daily.reindex(tradingday_list) #回测期间的资产收益率
pct_daily_new.head().append(pct_daily_new.tail())

def get_near_tradingday(date):
    return (get_price('000001.XSHG',end_date=date,count=1).index[0])
exchange_day = map(get_near_tradingday,date_list)
np.array(exchange_day)

weight_rp_df = pd.DataFrame()
risk_rp_df = pd.DataFrame()
weight_pcrp_df = pd.DataFrame()
risk_pcrp_df = pd.DataFrame()
t1=dt.datetime.now()
for date in date_list[:]:
    t2=dt.datetime.now()
    print ('计算到%s,已耗时%s秒'%(date,(t2-t1).seconds))
    date_pre_y = tradingday_list_tot[tradingday_list_tot.index(get_near_tradingday(date))-240]
    weight_and_risk_rp = get_smart_weight(pct_daily.loc[date_pre_y:get_near_tradingday(date)],method='risk parity', cov_adjusted=0, wts_adjusted=1)#.reindex(pct_daily.columns).fillna(0)
    weight_rp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_rp[0]
    risk_rp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_rp[1]
    weight_and_risk_pcrp = get_smart_weight(pct_daily.loc[date_pre_y:get_near_tradingday(date)],method='pc risk parity', cov_adjusted=0, wts_adjusted=1)#.reindex(pct_daily.columns).fillna(0)
    weight_pcrp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_pcrp[0]
    risk_pcrp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_pcrp[1]
    

weight_rp_df= weight_rp_df.T.reindex(tradingday_list).fillna(method='ffill')#风险平价模型的资产权重
weight_rp_df.head()

weight_pcrp_df= weight_pcrp_df.T.reindex(tradingday_list).fillna(method='ffill')#主成分风险平价模型的资产权重
weight_pcrp_df.head()

risk_rp_df= risk_rp_df.T.reindex(tradingday_list).fillna(method='ffill')#风险平价模型的资产风险贡献权重
risk_rp_df.head()

risk_pcrp_df= risk_pcrp_df.T.reindex(tradingday_list).fillna(method='ffill')#主成分风险平价模型的资产风险贡献权重
risk_pcrp_df.head()

pct_daily = pct_daily.reindex(tradingday_list)
pct_daily.head()

net_value=pd.DataFrame()
net_value['rp'] = (weight_rp_df*(pct_daily_new+1)).sum(axis=1).cumprod()
net_value['pcrp'] = (weight_pcrp_df*(pct_daily_new+1)).sum(axis=1).cumprod()

net_value['rp'] = ((weight_rp_df*pct_daily_new).sum(axis=1)+1).cumprod()
net_value['pcrp'] = ((weight_pcrp_df*pct_daily_new).sum(axis=1)+1).cumprod()

net_value[['rp','pcrp']].plot(figsize=(15,6),grid='on',rot=60)

weight_pcrp_df.reindex(exchange_day).plot(figsize=(15,8),rot=60,kind='area',ylim=(0,1))#PCRP模型资产权重图

risk_pcrp_df.reindex(exchange_day).plot(figsize=(15,8),rot=60,kind='area')#PCRP模型资产风险贡献图

weight_rp_df.reindex(exchange_day).plot(figsize=(15,8),rot=60,kind='area',ylim=(0,1))#RP模型资产权重图

risk_rp_df.reindex(exchange_day).plot(figsize=(15,8),rot=60,kind='area')#RP模型资产风险贡献图

pct_daily = prices.pct_change()
weight_wdc_rp_df=pd.DataFrame()
risk_wdc_rp_df=pd.DataFrame()
weight_wdc_pcrp_df=pd.DataFrame()
risk_wdc_pcrp_df=pd.DataFrame()
t1=dt.datetime.now()
for date in date_list[:]:
    t2=dt.datetime.now()
    print ('计算到%s,已耗时%s秒'%(date,(t2-t1).seconds))
    date_pre_y = tradingday_list_tot[tradingday_list_tot.index(get_near_tradingday(date))-240]
    weight_and_risk_wdc_rp = get_smart_weight(pct_daily.loc[date_pre_y:get_near_tradingday(date)],method='risk parity', cov_adjusted=1, wts_adjusted=1)
    #衰减加权的风险平价(wdc_RP)模型资产权重
    weight_wdc_rp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_wdc_rp[0]
    #衰减加权的风险平价(wdc_RP)模型风险贡献比例
    risk_wdc_rp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_wdc_rp[1]
    weight_and_risk_wdc_pcrp = get_smart_weight(pct_daily.loc[date_pre_y:get_near_tradingday(date)],method='pc risk parity', cov_adjusted=1, wts_adjusted=1)
    #衰减加权的主成分风险平价(wdc_PCRP)模型资产权重
    weight_wdc_pcrp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_wdc_pcrp[0]
    #衰减加权的主成分风险平价(wdc_PCRP)模型风险贡献权重
    risk_wdc_pcrp_df.loc[:,get_near_tradingday(date)] = weight_and_risk_wdc_pcrp[1]
    
weight_wdc_rp_df= weight_wdc_rp_df.T.reindex(tradingday_list).fillna(method='ffill')
weight_wdc_pcrp_df= weight_wdc_pcrp_df.T.reindex(tradingday_list).fillna(method='ffill')

#衰减加权的风险平价(wdc_RP)模型净值
net_value['wdc_rp'] = ((weight_wdc_rp_df*pct_daily_new).sum(axis=1)+1).cumprod()
#衰减加权的主成分风险平价(wdc_PCRP)模型净值
net_value['wdc_pcrp'] = ((weight_wdc_pcrp_df*pct_daily_new).sum(axis=1)+1).cumprod()
net_value[['wdc_pcrp','pcrp']].plot(figsize=(15,6),grid='on',rot=60)

#按月末价格进行重采样
prices_M = prices.resample('M',how='last')
#根据价格与均线关系判断择时信号
TF=(prices_M>pd.rolling_mean(prices_M,3)).shift(1)
#将择时信号与交易日进行对齐
TF['near_date']=map(get_near_tradingday,TF.index)
#根据择时信号调整资产权重（若该资产看多:仓位为100%或看空:仓位为0%）
weight_adjusted=(TF.set_index('near_date').reindex(tradingday_list).fillna(method='ffill'))
#无风险收益率
pct_norisk = 1.25e-4
#引入预期走势估计的RP模型净值，等式右边第二项为用无风险资产替换看空仓位部分的收益
net_value['TF_rp']=((weight_adjusted*weight_rp_df*pct_daily_new).sum(axis=1)+\
                    (1-(weight_adjusted*weight_rp_df).sum(axis=1))*pct_norisk+1).cumprod()
#引入预期走势估计的PCRP模型净值，等式右边第二项为用无风险资产替换看空仓位部分的收益
net_value['TF_pcrp']=((weight_adjusted*weight_pcrp_df*pct_daily_new).sum(axis=1)+\
                    (1-(weight_adjusted*weight_pcrp_df).sum(axis=1))*pct_norisk+1).cumprod()

#引入预期走势估计的wdc_RP模型净值，等式右边第二项为用无风险资产替换看空仓位部分的收益
net_value['TF_wdc_rp']=((weight_adjusted*weight_wdc_rp_df*pct_daily_new).sum(axis=1)+\
                    (1-(weight_adjusted*weight_wdc_rp_df).sum(axis=1))*pct_norisk+1).cumprod()
#引入预期走势估计的wdc_PCRP模型净值，等式右边第二项为用无风险资产替换看空仓位部分的收益
net_value['TF_wdc_pcrp']=((weight_adjusted*weight_wdc_pcrp_df*pct_daily_new).sum(axis=1)+\
                    (1-(weight_adjusted*weight_wdc_pcrp_df).sum(axis=1))*pct_norisk+1).cumprod()

net_value[['TF_rp','rp']].plot(figsize=(15,6),grid='on',rot=60)

net_value[['TF_pcrp','pcrp']].plot(figsize=(15,6),grid='on',rot=60)

weight_wdc_rp_df.reindex(exchange_day).plot(figsize=(15,8),rot=60,kind='area',ylim=(0,1))

(weight_adjusted*weight_wdc_rp_df).reindex(exchange_day).plot(figsize=(15,8),rot=60,kind='area',ylim=(0,1))

#使用pickle模块将数据对象保存到文件
save_data={}
save_data['weight_rp'] = weight_rp_df
save_data['weight_pcrp'] = weight_pcrp_df
save_data['weight_wdc_rp'] = weight_wdc_rp_df
save_data['weight_wdc_pcrp'] = weight_wdc_pcrp_df
save_data['risk_rp'] = risk_rp_df
save_data['risk_pcrp'] = risk_pcrp_df
save_data['risk_wdc_rp'] = risk_wdc_rp_df
save_data['risk_wdc_pcrp'] = risk_wdc_pcrp_df
save_data['net_value'] = net_value
content = pickle.dumps(save_data) 
write_file('save_data_v4.pkl', content, append=False)

pkl_file_read = read_file('save_data_v4.pkl')
load_data = pickle.load(StringIO(pkl_file_read))
net_value = load_data['net_value']
weight_rp_df = load_data['weight_rp']
weight_pcrp_df = load_data['weight_pcrp']
weight_wdc_rp_df = load_data['weight_wdc_rp']
weight_wdc_pcrp_df = load_data['weight_wdc_pcrp']

class backtest_result():
    def __init__(self,data):
        self.data = data
        self.total_returns = data.iloc[-1]-1
        self.annualized_returns = data.iloc[-1]**(250./len(data))-1
        self.annualized_volatility = data.pct_change().std()*(250.**0.5)
    def Max_Drawback(self):
        net_value=self.data
        max_value=0
        df_tmp=pd.DataFrame(net_value)
        df_tmp.columns=['value']
        for j in range(0,len(net_value),1):
            max_value=max(max_value,df_tmp.ix[j,'value'])
            df_tmp.ix[j,'drawback']=1-df_tmp.ix[j,'value']/max_value
            drawback=df_tmp['drawback'].max()
        return drawback
    def Sharpe(self):
        net_value=self.data
        bench_pct=0.03
        df_tmp=pd.DataFrame(net_value)
        df_tmp.columns=['value']
        df_tmp['pct']=df_tmp['value'].pct_change()
        annual_pct = df_tmp.ix[-1,'value']**(250./len(df_tmp))-1
        sharpe = (annual_pct-bench_pct)/(df_tmp['pct'].std()*250**0.5)
        return sharpe
    def Calmar(self):
        clamar = self.annualized_returns/self.Max_Drawback()
        return clamar


results=pd.DataFrame()
for c in ['rp','pcrp','wdc_rp','wdc_pcrp','TF_pcrp','wdc_pcrp','TF_wdc_pcrp']:
    nv = net_value[c]
    results.loc['累计收益率',c]=backtest_result(nv).total_returns
    results.loc['年化收益率',c]=backtest_result(nv).annualized_returns
    results.loc['年化波动率',c]=backtest_result(nv).annualized_volatility
    results.loc['最大回撤',c]=backtest_result(nv).Max_Drawback()
    results.loc['夏普比率',c]=backtest_result(nv).Sharpe()
    results.loc['Calmar比率',c]=backtest_result(nv).Calmar()
results