import pandas as pd
import numpy as np
import calendar
import matplotlib.pyplot as plt
import matplotlib.dates as mdate
from datetime import datetime
from dateutil.parser import parse
from scipy.optimize import minimize
TOLERANCE = 1e-10


# =============================================================================
# part1 :处理数据 ,通过时间区间[time1,time2]中的数据获取日收益率何周收益率
# =============================================================================
start = '20041231'
end = '20181231'


stock = pd.read_excel('C:\\Users\\Administrator\\Desktop\\Risk parity\\test.xlsx',sheet_name='股票')
bond = pd.read_excel('C:\\Users\\Administrator\\Desktop\\Risk parity\\test.xlsx',sheet_name='债券')
#goods = pd.read_excel('C:\\Users\\Administrator\\Desktop\\Risk parity\\stock&bond.xlsx',sheet_name='南华商品')


Name = [stock , bond] # 所要配置的资产
asset = ['stock' , 'bond'] #记录资产组合中资产的名称,方便后期作图


def _deal_data(Name):
    
    ## 首先实现表格的连接
    if len(Name) == 1 :
        data = Name[0]
    else:
        data = pd.merge(Name[0], Name[1], how='outer',on=['Date'])
    # 两个表格直接连接,多表循环
    i = 2
    while i < len(Name):
        data = pd.merge(data, Name[i], how='outer',on=['Date'])
        i = i + 1
    data = data.sort_values(by = 'Date') 
    data = data.reset_index(drop = True)
    # 返回经时间排序的资产收盘价数据
    
    ## 其次删除空值,选取目标日期
    data.dropna(axis=0 ,how ='any' ,inplace =True)
    data = data.reset_index(drop = True)
    
    data = data[(data['Date'] >= pd.to_datetime(start)) & (data['Date'] <= pd.to_datetime(end))]
    data = data.reset_index(drop = True)
    
    dataAll = pd.date_range(start ,end )
    dataAll = pd.DataFrame(dataAll)
    dataAll.columns = ['Date']
    data = pd.merge(dataAll, data, how='outer',on=['Date'])
    data = data.fillna(method='pad')
    # 返回包含每个日期的收盘价数据,非交易日附之前最近一个交易日数据

    data.dropna(axis=0 ,how ='any' ,inplace =True)
    data = data.reset_index(drop = True)    
    
    data.set_index('Date',inplace=True)
    #直接将指标设为日期以方便后续计算
    
    return data

test = _deal_data(Name)

def _day_return(Name):
    # 计算收盘价的日收益
    data = _deal_data(Name)
    name = data.columns
    for j in range(data.shape[1]):
        data['day_return_'+name[j]] = (data[name[j]]-data[name[j]].shift(1)) / data[name[j]].shift(1)
    # 返回dataFrame,前line列为资产收盘价 ,后line列为标的资产日收益率  
    # 由于是按日历日期计算的,因而由于节假日存在部分零值,为保证波动率的真实性,删除零值
#    data = 
#    data['day_return_中证全债']=4*data['day_return_中证全债']
    line = int(0.5*data.shape[1])
    return data.iloc[1::,line::]
  
def _day_return_leverage(Name,leverage):
    # 计算收盘价的日收益
    # leverage : 表示杠杆的倍数
    data = _deal_data(Name)
    name = data.columns
    for j in range(data.shape[1]):
        data['day_return_'+name[j]] = (data[name[j]]-data[name[j]].shift(1)) / data[name[j]].shift(1)
    # 返回dataFrame,前line列为资产收盘价 ,后line列为标的资产日收益率  
    # 由于是按日历日期计算的,因而由于节假日存在部分零值,为保证波动率的真实性,删除零值    
    data['day_return_中证全债']=leverage*data['day_return_中证全债']-(leverage-1)*0.03/365
    line = int(0.5*data.shape[1])
    
    return data.iloc[1::,line::]

def _period_return(data , period):
    # 计算收盘价的月收益率 
    # 输入: period 对应的是 'W':周 ;'M':月 ;'3M':季度 ;'6M':半年 

    __df= data.resample(period).last()
    dataw = pd.concat([data.iloc[0,0::].T,__df.T],axis=1).T 
    # 返回周末收盘价,同时保留表格第一行数据
    name = dataw.columns    
    for j in range(dataw.shape[1]):
        dataw['return_'+name[j]] = (dataw[name[j]]-dataw[name[j]].shift(1)) / dataw[name[j]].shift(1)
    # 返回dataFrame,前line列为资产收盘价 ,后line列为标的资产日收益率   
    line = int(0.5*dataw.shape[1])
    return dataw.iloc[1::,line::]


def _calculate_cov(data):
    return data.cov().values


return_daily = _day_return(Name) #日收益率
#return_weekly = _weekly_return(Name) # 周收益率

def MaxDrawdown(asset_value):
#     '''最大回撤率'''
    # np.maximum.accumulate():生成一列当日之前历史最高价值的序列。
    # argmax(f(x))是使得 f(x)取得最大值所对应的变量x
    i = np.argmax((np.maximum.accumulate(asset_value) - asset_value) / np.maximum.accumulate(asset_value))  # 结束位置
#    if i == 0:
#        return 0
    j = np.argmax(asset_value[:i])  # 开始位置
    return (asset_value[j] - asset_value[i]) / (asset_value[j])


def asset_statistics_model(Name,asset0):
    #计算资产的各项指标
    
    asset_price_data = _deal_data(Name)
    asset_price_data.columns = asset0
    
    A_return = _period_return(asset_price_data , 'M') #'M'表示当前是一月一调仓
    A_return_yearly = _period_return(asset_price_data , 'Y')
#    A_return_daily = _period_return(asset_price_data , 'D')
    
    total_return = ((A_return + 1).T.values.cumprod(1)).T[-1]-1
    total_return = pd.DataFrame(total_return).T
    total_return.columns = asset0

    annual_return = pow(total_return + 1, 1.0/13) - 1
    
    annual_return_mean = A_return_yearly.mean()
    annual_return_mean = pd.DataFrame(annual_return_mean).T
    annual_return_mean.columns = asset0
    
    annual_return_std = A_return_yearly.std()
    annual_return_std = pd.DataFrame(annual_return_std).T
    annual_return_std.columns = asset0
    
#    annual_Volatility = (A_return_daily.std() / A_return_daily.mean())/ np.sqrt(1/365)
    
    yearly_sharpe = annual_return_mean/annual_return_std
#    yearly_sharpe = pd.DataFrame(yearly_sharpe).T
    yearly_sharpe.columns = asset0    
    
    
    asset_maxDrawdown = []
    for i in range(len(Name)):
        asset_maxDrawdown.append(MaxDrawdown(asset_price_data[asset0[i]]))
    asset_maxDrawdown = pd.DataFrame(asset_maxDrawdown).T
    asset_maxDrawdown.columns = asset0

    asset_statistics = pd.concat([total_return , annual_return , annual_return_mean , annual_return_std , yearly_sharpe , asset_maxDrawdown])
    
    asset_statistics = asset_statistics.T
    asset_statistics.columns = ['total_return' , 'annual_return' , 'annual_return_mean' , 'annual_return_std' , 'yearly_sharpe' , 'asset_maxDrawdown']
    
    return asset_statistics.T


asset_behaviour_statistics = asset_statistics_model(Name , asset)


# =============================================================================
# part2 : Risk parity 二次规划问题求解 
# =============================================================================

#########构建求解二次规划问题的方法


#### step1: 构建二次规划问题并优化
    
def _allocation_risk(weights, covariances):
    # 计算给定资产权重下资产组合的风险
    # 输入: weights:资产权重,行向量(注意使用np.matrix()将list转化为1*n为行向量) ; covrances:资产的协方差矩阵
    portfolio_risk = np.sqrt((weights * covariances * weights.T))[0, 0]    
    return portfolio_risk

def _assets_risk_contribution_to_allocation_risk(weights, covariances):
    # 计算资产组合中单个资产贡献的风险(total risk contribution)
    # 输入:与_allocation_risk()函数一致 ; 输出为 N*1的列向量(N表示资产的数量)
    portfolio_risk = _allocation_risk(weights, covariances)    
    assets_risk_contribution = np.multiply(weights.T,covariances*weights.T)/portfolio_risk
    return assets_risk_contribution

def _risk_budget_objective_error(weights, args):
    # 二次规划问题的目标函数 :weights为规划问题求解的变量 ;args是参数的传递
    # 参数1为资产的协方差阵
    covariances = args[0]
    # 参数2为单个资产对于资产组合的目标贡献率——1/N
    assets_risk_budget = args[1]
    # convert the weights to a matrix
    weights = np.matrix(weights)
    # We calculate the risk of the weights distribution
    portfolio_risk = _allocation_risk(weights,covariances)
    # We calculate the contribution of each asset to the risk of the weights distribution
    assets_risk_contribution = _assets_risk_contribution_to_allocation_risk(weights,covariances)
    # We calculate the desired contribution of each asset to the risk of the weights distribution
    # asset_risk_target 为行向量
    assets_risk_target = np.asmatrix(np.multiply(portfolio_risk, assets_risk_budget))
    # Error between the desired contribution and the calculated contribution of each asset
    error = sum(np.square(assets_risk_contribution - assets_risk_target.T))[0, 0]
    # It returns the calculated error ; the aim is minimize error
#    error = sum(np.square(assets_risk_contribution[0]-assets_risk_contribution[1]))
    return error


def _get_risk_parity_weights(covariances, assets_risk_budget, initial_weights):
    # Restrictions to consider in the optimisation: only long positions whose sum equals 100%
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1.0},{'type': 'ineq', 'fun': lambda x: x})
    # Optimisation process in scipy
    optimize_result = minimize(fun=_risk_budget_objective_error,
                               x0=initial_weights,
                               args=[covariances, assets_risk_budget],
                               method='SLSQP',
                               constraints=constraints,
                               tol=TOLERANCE,
                               options={'disp': False})
    # Recover the weights from the optimised object
    weights = optimize_result.x
    # It returns the optimised weights
    return weights 

#### step2: 确定portfolio中asset_weight

def get_theday_of_last_year(date):
    date_parse = parse(date)
    year = date_parse.year - 1
    month = date_parse.month
    day = date_parse.day
    
    if month==2:
    #此时涉及润2月的情况(普通闰年:能被4整除但不能被100整除的年份为普通闰年。;世纪闰年:能被400整除的为世纪闰年。)
        if (year+1) % 4 == 0:
            if (year+1) % 100 == 0:
                day = day - 2
            else:
                day = day - 1
        
    last_date = datetime(year,month,day)
    return last_date


def get_theday_of_next_year(date):
    date_parse = parse(date)
    year = date_parse.year + 1
    month = date_parse.month
    day = date_parse.day
    
    if month==2:
        #此时涉及润2月的情况(普通闰年:能被4整除但不能被100整除的年份为普通闰年。;世纪闰年:能被400整除的为世纪闰年。)
        if (year+1) % 4 == 0:
            if (year+1) % 100 == 0:
                day = day - 2
            else:
                day = day - 1
                
    next_date = datetime(year,month,day)
    return next_date


def get_nmonth_endday(time,month_number):
    
    """
    获取后n个月最后一天的日期
 输入 time 以字符串表示的日期
 输入 month_number 为后n个月,可正可负
    return: 返回 list 日期
 需要import datatime,parse,calendar
    """
    today = parse(time)
    year = today.year
    month = today.month
    month+=month_number
    while month>12:
        month -= 12
        year+=1
    while month <= 0:
        month += 12
        year -= 1   
    firstDayWeekDay, monthRange = calendar.monthrange(year, month)
    res=datetime(year,month,monthRange).strftime('%Y%m%d')
    return res


def _time_interval_build(number):
    # number:时间区间所包含的月份 该函数通常用于生成月末,季度末,半年末的日历日期
    # 输入:number(int);输出:date_month(DataFrame)
    dateAll = pd.date_range(get_theday_of_next_year(start) , end , freq = 'M')
    dateAll = pd.DataFrame(dateAll)
    date_month = dateAll.iloc[(dateAll.index%number) ==0]        
    date_month.columns = ['Date']
    date_month = date_month.reset_index(drop = True)
    return date_month

# =============================================================================
# def _get_nonzero_value(data):
#     
#     return
# =============================================================================
    
# =============================================================================
# def get_weights(_return,time):
#     
#     # 输入:_return: return on aeest ; time: date of rebalance
#     
#     data = _return.iloc[(_return.index <= time)&(_return.index> get_nmonth_endday(time))]
# 
#     covariances = data.cov().values*250
#     # We calculate the covariance matrixre
#     # The desired contribution of each asset to the portfolio risk: we want all asset to contribute equally
#     assets_risk_budget = [1/data.shape[1]]*data.shape[1]
#     # Initial weights: equally weighted
#     init_weights = [1/data.shape[1]]*data.shape[1]
#     # Optimisation process of weights
#     weights = _get_risk_parity_weights(covariances,assets_risk_budget, init_weights)#两行一列的向量
#     # Convert the weights to a pandas Series
#     weights = pd.Series(weights, index=data.columns, name='weight')
#     
#     return weights
# =============================================================================


def get_weights(_return,time):
    
    # 输入:_return: return on aeest ; time: date of rebalance
    
    data = _return.iloc[(_return.index <= time)&(_return.index> get_theday_of_last_year(time))]
    #data反映了回看周期 ,即以过去多长时间的数据计算协方差阵
    covariances = data.cov().values*250 #*250表示年化
    # We calculate the covariance matrixre
    # The desired contribution of each asset to the portfolio risk: we want all asset to contribute equally
    assets_risk_budget = [1/data.shape[1]]*data.shape[1]
    # Initial weights: equally weighted
    init_weights = [1/data.shape[1]]*data.shape[1]
    # init_weights = [0.2,0.8]
    # Optimisation process of weights
    weights = _get_risk_parity_weights(covariances,assets_risk_budget, init_weights)#两行一列的向量
    # Convert the weights to a pandas Series
    weights = pd.Series(weights, index=data.columns, name='weight')
    
    return weights


# =============================================================================
# #def all_portfolio_weight(number):
# #    # 返回的是投资组合权重的字典序
# #    result={}
# #    date_m = _time_interval_build(number)
# #    
# #    
# #    if number == 6:
# #        _return=return_weekly
# #    else:
# #        _return=return_daily
# #    # 如果调仓频率为半年,则按照周收益率计算协方差阵,若调仓频率为单月或季度,则按照日收益率计算协方差阵
# #    
# #    for i in range(date_m.shape[0]-1):    
# ##    result[date_m.iloc[i+1].values[0]]=pd.DataFrame(get_weights(return_daily , date_m.iloc[i].values[0] , date_m.iloc[i+1].values[0])).iloc[0,0::]
# #        result[date_m.iloc[i+1].values[0]]=pd.DataFrame(get_weights(_return , date_m.iloc[i].values[0] , date_m.iloc[i+1].values[0])).T
# #  
# #    return result
# =============================================================================


def all_portfolio_weight(_return , number):
    # 返回的是投资组合权重的字典序
    result=[]
    result = pd.DataFrame(result)#生成空的dataFrame以方便记录portfolio中asset的权重
    
    date_m = _time_interval_build(number)  # date_m表示的是调仓的频率
    
    result = pd.concat([result ,date_m],axis=1)
    result.set_index('Date',inplace=True) #result的index指标对应须调仓日
        
    
    for i in range(len(_return.columns.values)):
        result[_return.columns.values[i]]=0.0
    # 通过收益率确定对应资产
    
    for i in range(result.shape[0]):    
        for j in range(result.shape[1]):
            result[result.columns.values[j]][date_m.iloc[i].values[0]] = get_weights(_return ,str( date_m.iloc[i].values[0]) )[j]

  
    return result

  
#weight_weekly_semiannually = all_portfolio_weight(return_weekly , 6)
weight_daily_monthly = all_portfolio_weight(return_daily , 1)
weight_daily_3monthly = all_portfolio_weight(return_daily , 3)
# =============================================================================
# part3 : backtest
# 回测期:200601-201908
# =============================================================================

backtest_start = '20051231'
backtest_end = '20181231'

def plt_weight_change(weight):
## 功能:采用堆叠柱状图反映资产权重变化
    #prepare data
    weight.columns=asset
    
    weight_value = [] #记录各资产的权重值,以满足堆叠折线图的绘制过程中,stackplot函数对于纵坐标的要求
    for i in range(len(asset)):
        weight_value.append(weight[asset[i]])
    
    # Draw
#    plt.figure(figsize = (12,6) , dpi = 60 )
    colors = ['powderblue','steelblue','darkturquoise','coral','darkcyan','gold','teal','pink','y','r','snow','b','k','g','orange','c','bisque','brown','lime','aqua']
        
    plt.stackplot(weight.index, weight_value, labels=asset, colors=colors)
    
    plt.legend()
    plt.show()
   
    return

#plt_weight_change(weight_daily_3monthly)


def _rebalance_portfolio_value(Name):
    # 输入:资产组合中各资产的权重和日收益率
    # 计算资产组合的单位净值以及累计收益
    
    _return = _day_return(Name)
    weight = all_portfolio_weight(_return , 1)
    asset_value0 = 1
    asset_daily_value = []
    asset_return = 1
    
    for i in range(weight.shape[0]-1):
        
        asset_value = asset_value0 *weight.iloc[i]
        _return_peiod = _return[(_return.index > weight.index[i]) & (_return.index <= weight.index[i+1])]
        
        for j in range(_return_peiod.shape[0]):
            asset_value = np.multiply( ( _return_peiod.iloc[j]+1) , asset_value)
            asset_daily_value.append(pd.DataFrame(asset_value).T)
        
        asset_return = asset_return * (np.sum(asset_value)/asset_value0)
        asset_value0 = np.sum(asset_value)
    
    asset_daily_value = pd.concat(asset_daily_value)
    asset_daily_value.columns = asset
    
    return [asset_daily_value,asset_return]

#portfolio_value  = _rebalance_portfolio_value (Name)[0]
#test = _rebalance_portfolio_value(Name)


def _rebalance_portfolio_leverage_value(Name,leverage):
    # 输入:资产组合中各资产的权重和日收益率
    # 计算资产组合的单位净值以及累计收益
    
    _return = _day_return_leverage(Name,leverage)
    weight = all_portfolio_weight(_return , 1)
    asset_value0 = 1
    asset_daily_value = []
    asset_return = 1
    
    for i in range(weight.shape[0]-1):
        
        asset_value = asset_value0 *weight.iloc[i]
        _return_peiod = _return[(_return.index > weight.index[i]) & (_return.index <= weight.index[i+1])]
        
        for j in range(_return_peiod.shape[0]):
            asset_value = np.multiply( ( _return_peiod.iloc[j]+1) , asset_value)
            asset_daily_value.append(pd.DataFrame(asset_value).T)
        
        asset_return = asset_return * (np.sum(asset_value)/asset_value0)
        asset_value0 = np.sum(asset_value)
    
    asset_daily_value = pd.concat(asset_daily_value)
    asset_daily_value.columns = asset
    
    return [asset_daily_value,asset_return]

def portfolio_statistics(Name):
    
    portfolio_value = _rebalance_portfolio_leverage_value(Name,3)[0]
    portfolio_value['Asset_sum'] = portfolio_value.apply(lambda x: x.sum(), axis=1)
    portfolio = pd.DataFrame(portfolio_value['Asset_sum'])
    date = pd.DataFrame( portfolio.index.values)
    portfolio = portfolio.reset_index(drop = True)
    portfolio = pd.concat([date,portfolio],axis = 1 )
    portfolio.columns = ['Date','资产组合']
    portfolio_statistics = asset_statistics_model([portfolio],['资产组合'])
 
    return portfolio_statistics


def plt_asset_value(Name):
## 功能:绘制资产以及资产组合的价格变化
    
    # Prepare data
    portfolio_value = _rebalance_portfolio_leverage_value (Name,3)[0]
    portfolio_value['Asset_sum'] = portfolio_value.apply(lambda x: x.sum(), axis=1) # 确定资产组合的每日价值
    
    data = _deal_data(Name).iloc[(_deal_data(Name).index >= portfolio_value.index[0])&(_deal_data(Name).index <= portfolio_value.index[-1])]
    weight = pd.DataFrame(1/data.iloc[0]).T.values
    data = np.multiply(weight,data)    # 归一化处理
    data.columns = asset #data反映的是资产组合中个资产的价值
    
    # Plot 
    
# =============================================================================
#     fig,ax1 = plt.subplots(figsize = (10,8))
#     ax2 = ax1.twinx()    # mirror the ax1   
#     # 由于股票的变动幅度较大,为更好的反映各个资产价格的变化,将股票价格单独设主坐标轴表示
#     ax1.plot(data.index, data[asset[0]] ,'g',label = asset[0])
#     ax2.plot(data.index, data[asset[1]],'b')
#     ax2.plot(portfolio_value.index, portfolio_value['Asset_sum'],'r')
#     plt.show()
# =============================================================================
      
    for i in range(data.shape[1]):
        plt.plot(data.index , data[asset[i]] ,label=asset[i])    
    plt.plot(portfolio_value.index , portfolio_value['Asset_sum'] ,label = 'Portfolio_value')

    # Decorations

    plt.xlabel('Date')
    plt.ylabel('price')
    plt.legend()
    plt.show()   
    return

plt_asset_value(Name)


def risk_contribution(Name):
    # 根据二次规划的解:计算各项资产对于资产组合的实际贡献程度
#    portfolio_value  = _rebalance_portfolio_value (Name) # 资产组合中各资产的净值
    
#    _return = _day_return_leverage(Name,3) #资产组合中各资产的日收益率
    _return = _day_return_leverage(Name,3)
    weight = all_portfolio_weight(_return , 1) #资产组合中各资产的权重
      
    contribution = []
    
    for i in range(weight.shape[0]-1):
         _return_peiod = _return[(_return.index > weight.index[i]) & (_return.index <= weight.index[i+1])]         
         weights = np.matrix(pd.DataFrame(weight.iloc[i]).T)
         covariances = np.matrix(_return_peiod.cov())
         portfolio_risk = _allocation_risk(weights, covariances)
         asset_risk = _assets_risk_contribution_to_allocation_risk(weights, covariances)
         asset_contribution = pd.DataFrame(asset_risk/portfolio_risk).T
#         asset_contribution= pd.concat( weight.index[i] , asset_contribution)
         contribution.append(asset_contribution)
         
    contribution = pd.concat(contribution)
    contribution.columns = asset    
    contribution.index = pd.date_range(weight.index[1] , weight.index[-1] ,freq = 'M')
        
    return contribution 


asset_actual_risk_contribution = risk_contribution(Name)