最近研究了下最优风险资产组合这个题目。本小白在金融领域是个纯粹的初学者，开始的时候，有点不知所措。

后来在网上找了下计算最优风险资产组合、有效边界、资本市场线的资料以及程序实现。发现资料虽然不少，但是系统的讲解还是不多。也有系统讲解的书，但是感觉并不是很通俗，需要预先打好基础。

所以，经过了一段时间的学习，归纳了参考书的内容，觉得很有必要做下总结笔记，温故而知新。

本文为学习笔记，肯定会有错误或者主观臆想之处，请多多包涵。。。
尽管如此，还是切入正题吧。

大概的思路是：首先获取10支股票数据，然后根据Markowitz均值-方差投资组合理论，分析10支股票的收益率和波动率，生成大量随机风险组合，从中可以发现夏普率较高的随机组合。然后通过这些随机组合，计算有效边界。最后，引入无风险资产，根据有效边界，画出资本市场线。资本市场线与有效边界相切的点即为最优风险资产组合所在的点。该点的组合权重分配、组合收益率和波动率即为我们希望获得的结果。因为波动率代表风险，所以最优风险资产组合是研究收益率和波动率的。

Python是非常适合进行数据分析的语言。本文使用Anaconda3(64 bit)的Python发行版。

该版本集合了大量常用的Python包，比如我们之后会用到的numpy、pandas、matplotlib、scipy.optimize和scipy.interpolate。除此之外，为了获得股票数据，我们需要额外安装tushare包。在Anaconda Prompt (命令行工具)里面，键入pip install tushare，就可以装好此包。Anaconda3(64 bit)里面集成了Spyder编辑器，可以用它运行代码并查看结果。

作为刚接触Python和Spyder不久的人，我觉得这套工具确实赞。左边编辑代码，运行。右边的IPython Console里面就可以随时交互式查看各变量结果。

获取股票数据

为了获得股票数据，本来可以使用pandas_datareader从Yahoo下载数据。
可是不知为何，Yahoo突然不提供股票数据下载服务。于是我只好使用国内的tushare。

加载库：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import tushare as ts
import scipy.optimize as sco
import scipy.interpolate as sci
from scipy import stats

然后选10支股票。随意选了几个好股票：

             '000651',       ##格力电器   
             '600519',       ##贵州茅台  
             '601318',       ##中国平安               
             '000858',       ##五粮液  
             '600887',       ##伊利股份  
             '000333',       ##美的集团   
             '601166',       ##兴业银行  
             '600036',       ##招商银行  
             '601328',       ##交通银行  
             '600104'        ##上汽集团  

tushare的get_hist_data函数获得的是不复权数据。我们需要复权数据可以用get_h_data函数。
get_h_data('000651', start='2014-05-28', end='2017-05-26')可以获得股票000651/格力电器的前复权数据。返回的数据是dataframe格式。dataframe可以看成是二维数据表。提取其中的’close’列，即可获得收盘价。

df = pd.DataFrame()
df = ts.get_h_data('000651', start='2014-05-28', end='2017-05-26') #前复权
s000651 = df['close']
s000651.name = '000651'

依葫芦画瓢，可以获得所有10支股票的收盘价。注意每个dataframe的一列是一个Series。将很多Series组合起来，就是一个dataframe。
将所有10支股票的收盘价格组合起来，形成我们待研究的数据：

data = pd.DataFrame({'000651':s000651, '600519':s600519, '601318':s601318, '000858':s000858, '600887':s600887, '000333':s000333,
                     '601166':s601166, '600036':s600036, '601328':s601328, '600104':s600104})

如果某支股票某天有缺失值，则去掉这个日期的所有股票信息，这样让所有股票的日期保持一致。

data = data.dropna()

可以查看每只股票日期下标是否一致等信息。

print(data.info())

可以将此收盘价表保存成Excel格式文件，则每次只需读此本地文件，而不必每次依靠tushare下载数据。

data.to_excel('stock_data2.xlsx')

前5行大致如下：

再利用pandas读取10支股票的excel数据。

data = pd.read_excel('stock_data2.xlsx')

data.index = data['date'].tolist() #将date列拷贝一份，设为index列
data.pop('date') #移除date列，因为date列已经成为index列了，以后只用看index

将个股价格与其初始值比较并且用100标准化，得出个股在相同初始条件的情况下的走势情况。

(data / data.ix[0] * 100).plot(figsize=(10, 8), grid=True) 

初步统计分析

现在需要计算个股的收益率。金融计算收益率的时候大部分用对数收益率 (Log Return) 而不是用算数收益率。用当天的收盘价与前一天相比较。

#计算对数收益率。金融计算收益率的时候大部分用对数收益率 (Log Return) 而不是用算数收益率
log_returns = np.log(data / data.shift(1))

前5行大致如下：

画出每只股票收益率的直方图，了解一下分布情况。

log_returns.hist(bins=50, figsize=(12, 9))

直方图

可以看到每支股票的分布形状很类似正态分布，但是正态性不强。

投资组合优化

尽管如此，Markowitz均值-方差投资组合理论需要假设正态分布收益率。而投资组合的风险取决于投资各组合中资产收益率的相关性。这样，年化收益率和协方差矩阵就是我们需要计算的。

#使用对数收益率为收益率
rets = log_returns
#计算年化收益率
year_ret = rets.mean() * 252
#计算协方差矩阵
year_volatility = rets.cov() * 252

年化收益率：

000333    0.642380
000651    0.512573
000858    0.601155
600036    0.497143
600104    0.470550
600519    0.703341
600887    0.324943
601166    0.328918
601318    0.457654
601328    0.329218

协方差矩阵：

协方差矩阵

根据理论，风险需要分散，每个股票都会有一定比例的投资权重。
一个资产组合的收益率（均值），为组合中个股收益率（均值）的权重之和。
比如对于三个资产（股票）的组合均值收益率情况：

而三个资产（股票）的组合的方差（波动情况）则为：

对于有多个资产的资产组合的均值和方差可以描述为：

我们使用随机数来随机模拟权重分配。即随机模拟权重分配给10个股票。

#我们一共有10支股票
number_of_assets = 10
#生成10个随机数
weights = np.random.random(number_of_assets)
#将10个随机数归一化，每一份就是权重，权重之和为1
weights /= np.sum(weights)

假设有5000组投资组合，每一组投资组合由一组随机权重组成，则可以这样产生5000组投资组合：

  portfolio_returns = []
  portfolio_volatilities = []
  for p in range (5000):
      weights = np.random.random(number_of_assets)
      weights /= np.sum(weights)
      portfolio_returns.append(np.sum(rets.mean() * weights) * 252)
      portfolio_volatilities.append(np.sqrt(np.dot(weights.T, 
                        np.dot(rets.cov() * 252, weights))))
    

得到5000个组合收益率和波动率

portfolio_returns = np.array(portfolio_returns)
portfolio_volatilities = np.array(portfolio_volatilities)

将5000个组合收益率和波动率用散点图描述，并且计算每个夏普率
哪个点的夏普率越高，证明哪个组合的权重分配越优，虽然还没有找到最优的。

作图：

plt.figure(figsize=(9, 5)) #作图大小
plt.scatter(portfolio_volatilities, portfolio_returns, c=portfolio_returns / portfolio_volatilities, marker='o') #画散点图
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')

5000个组合收益率、波动率和夏普率散点图

每个点对应某个投资组合，该点有其对应的收益率和波动率（标准差），其颜色为对应的夏普率。可见，越往左上方，夏普率越高。

我们现在再描述一下某个投资组合。这个投资组合是这样的一个函数，即输入权重分配，输出该组合的收益率、波动率和夏普率。函数定义如下：

def statistics(weights):        
    #根据权重，计算资产组合收益率/波动率/夏普率。
    #输入参数
    #==========
    #weights : array-like 权重数组
    #权重为股票组合中不同股票的权重    
    #返回值
    #=======
    #pret : float
    #      投资组合收益率
    #pvol : float
    #      投资组合波动率
    #pret / pvol : float
    #    夏普率，为组合收益率除以波动率，此处不涉及无风险收益率资产
    #

    weights = np.array(weights)
    pret = np.sum(rets.mean() * weights) * 252
    pvol = np.sqrt(np.dot(weights.T, np.dot(rets.cov() * 252, weights)))
    return np.array([pret, pvol, pret / pvol])

假设某投资组合的那个点的夏普率最高，即可认为该组合是夏普率最优的投资组合。夏普率最优的那个点不一定是刚才5000个点里面的其中一个，而是需要通过优化算法，找到一个恰当的权重分配输入，导致输出的夏普率最大。

我们之前引入了优化算法包import scipy.optimize as sco，使用其中的最小化优化算法sco.minimize。

则我们的最大化夏普率问题可以转变为最小化负的夏普率问题。定义输入权重分配，返回负夏普率的函数为：

def min_func_sharpe(weights):
    return -statistics(weights)[2]

优化算法即为：

opts = sco.minimize(min_func_sharpe, number_of_assets * [1. / number_of_assets,], method='SLSQP',  bounds=bnds, constraints=cons)

number_of_assets * [1. / number_of_assets,]为权重初始值[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]，'SLSQP'为Sequential Least Squares Programming方法。

边界条件为：
bnds = tuple((0, 1) for x in range(number_of_assets)) 即每个权重需要在0到1之间。
约束条件为：
cons = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1})即权重之和为1。

结果opts['x']即为最大夏普率的投资组合的权重分配。精确到小数点后三位opts['x'].round(3)，结果为：

[ 0.     0.199  0.     0.152  0.     0.553  0.     0.     0.097  0.   ]

statistics(opts['x']).round(3)可以获得最大夏普率的投资组合的收益率、波动率和夏普率。结果为：

[ 0.61   0.369  1.656]

我们的结论为：假定最大夏普率的风险投资组合是最优风险风险投资组合。由我们10支股票组成的例子可以知道。19.9%的权重买000651格力电器，15.2%的权重买600036招商银行，53.3%的权重买600519贵州茅台，9.7%的权重买601318中国平安。这样的组合，根据以往数据，可以分析得出组合年化收益率为61%，波动率为36.9%，夏普率为1.656。

依葫芦画瓢，如果我们想知道最小方差的投资组合，则可以定义这样的函数：

def min_func_variance(weights):
return statistics(weights)[1] ** 2

求解：

optv = sco.minimize(min_func_variance, number_of_assets * [1. / number_of_assets,], method='SLSQP', bounds=bnds, constraints=cons)

有效边界

现在我们可以研究有效边界。有效边界实际上就是，有效边界上的每一个点即为给定收益率情况下拥有最小波动率的投资组合的点。所以计算有效边界上的点，可以描述成，已知该点的收益率，求权重组合，使得该点波动率最小。

输入权重，输出波动率的函数为：

def min_func_port(weights):
    return statistics(weights)[1]  

给定收益率为，从0.35到0.65之间30等份的值。制造一个线性空间给我们的收益率数组，预设好我们的目标收益率:

target_returns = np.linspace(0.35, 0.65, 30)

目标波动率是我们需要求解的

target_volatilities = []

针对每个target_returns中的收益率，使得目标波动率最小。
此处的约束多了一个，即目标收益率已经给定，也就是target_returns空间中的值。
minimize函数返回了对象res，res中挑选fun成员，即为最优函数输出值statistics(weights)[1]，也就是最小波动率。

for tret in target_returns:
    cons = ({'type': 'eq', 'fun': lambda x:  statistics(x)[0] - tret},
            {'type': 'eq', 'fun': lambda x:  np.sum(x) - 1})
    res = sco.minimize(min_func_port, number_of_assets * [1. / number_of_assets,], method='SLSQP',
                       bounds=bnds, constraints=cons)
    target_volatilities.append(res['fun'])

作图：

#画散点图
plt.figure(figsize=(9, 5))
#圆点为随机资产组合
plt.scatter(portfolio_volatilities, portfolio_returns,
            c=portfolio_returns / portfolio_volatilities, marker='o')
#叉叉为有效边界            
plt.scatter(target_volatilities, target_returns,
            c=target_returns / target_volatilities, marker='x')
#红星为夏普率最大值的资产组合            
plt.plot(statistics(opts['x'])[1], statistics(opts['x'])[0],
         'r*', markersize=15.0)
#黄星为最小方差的资产组合            
plt.plot(statistics(optv['x'])[1], statistics(optv['x'])[0],
         'y*', markersize=15.0)
            # minimum variance portfolio
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')

30个叉叉点为有效边界。红星为夏普率最大值的资产组合。黄星为最小波动率的资产组合。

资本市场线

资本市场线穿过无风险资产收益率的点和有效边界相切。对于有效边界而言，上半部分才是有意义的。或者说，上半部分才是有效边界，下半部分是无效边界。这需要找到最小波动率的值，这个点是一个拐点。

如下，找到target_volatilities最小值的下标

ind = np.argmin(target_volatilities)

通过下标，可以提取边界上半部分的点的集合

upper_half_volatilities = target_volatilities[ind:]
upper_half_returns = target_returns[ind:]

此时，我们对有效边界上半部分的点的集合进行插值

tck = sci.splrep(upper_half_volatilities, upper_half_returns)
#tck参数用于构造有效边界函数f(x)
def f(x):
    #有效边界函数 (样条函数逼近).
    return sci.splev(x, tck, der=0)

#同时也构造有效边界函数f(x)的一阶导数函数df(x)
def df(x):
    #有效边界函数f(x)的一阶导数函数 
    return sci.splev(x, tck, der=1)

假设无风险资产收益率为2%，即risk_free_return=0.02
需要满足三个条件:

条件1： t(0) = a, a = risk_free_return, risk_free_return=0.02 (2%)
条件2： t(x) = a + b*x
条件3： dt(x) = b

输入参数p, p=(a,b,x)集合， 即p[0]=a, p[1]=b, p[2]=x
eq1 , eq2 , eq3 为满足上述条件的等式
最后目的使
0 = risk_free_return - p[0]
0 = risk_free_return + p[1] * p[2] - f(p[2])
0 = p[1] - df(p[2])

def equations(p, risk_free_return=0.02):
    eq1 = risk_free_return - p[0]
    eq2 = risk_free_return + p[1] * p[2] - f(p[2])
    eq3 = p[1] - df(p[2])
    return eq1, eq2, eq3

输入a,b,x初始值，返回a,b,x的求解值

opt = sco.fsolve(equations, [0.01, 0.50, 0.15])

结果为

[ 0.02 ,  1.60220549,  0.37100371]

a求得为0.02，b为1.60。
X波动率为0.37，基本上约等于夏普率最大值的资产组合处的波动率值（红星处）

作图：

plt.figure(figsize=(9, 5))
#圆点为随机资产组合
plt.scatter(portfolio_volatilities, portfolio_returns,
            c=(portfolio_returns - 0.02) / portfolio_volatilities, marker='o')
#绿色线为有效边界
plt.plot(upper_half_volatilities, upper_half_returns, 'g', lw=4.0)

#设定资本市场线CML的x范围从0到0.6           
cml_x = np.linspace(0.0, 0.6)
#带入公式a+b*x求得y,作图
plt.plot(cml_x, opt[0] + opt[1] * cml_x, lw=1.5)
#标出资本市场线与有效边界的切点，红星处            
plt.plot(opt[2], f(opt[2]), 'r*', markersize=15.0) 
plt.grid(True)
plt.axhline(0, color='k', ls='--', lw=2.0)
plt.axvline(0, color='k', ls='--', lw=2.0)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharpe ratio')

image.png

绿色线条插值函数为有效边界。红星处为最优风险组合。直线为资本市场线。红星处的点和之前的结论一样，所以红星处的风险投资权重分配也一样。

置信区间

最后研究一下置信空间。这部分我不确定做得对不对，还需请教。。。

以最优风险资产权重分配作为投资组合的股票权重分配，从2014-06-09开始到3年期末。投资组合每天的收益率为actual_opt_portfolio_returns：

opt_weights = opts['x'].round(3)  # opt_weights是最优风险组合权重
actual_opt_portfolio_returns = np.sum(rets * opt_weights, axis=1)  #axis=1 是按行求和  一行中的每一列相加

针对投资组合过往3年中的每一天，都观察当天之前的7日年化收益率。则可以计算出每天考量的7日年化收益率（单利情况）
seven_day_year_returns：

seven_day_year_returns = 252 * (actual_opt_portfolio_returns.shift(7) + actual_opt_portfolio_returns.shift(6) + actual_opt_portfolio_returns.shift(5) + actual_opt_portfolio_returns.shift(4) + actual_opt_portfolio_returns.shift(3) + actual_opt_portfolio_returns.shift(2) + actual_opt_portfolio_returns.shift(1)) / 7

去掉开头的空值。

seven_day_year_returns = seven_day_year_returns.dropna()   

最后计算置信空间confidence_range。

confidence_range = stats.t.interval(0.99, len(seven_day_year_returns)-1, np.mean(seven_day_year_returns), stats.sem(seven_day_year_returns))

结论是置信度为0.99时，7日年化收益率均值置信区间为 （0.34，0.85）