Markowitz投资组合优化的Python应用
我们都知道，马科维茨（Harry M. Markowitz）因其在1952年提出投资组合选择理论（Portfolio Selection）获得了诺贝尔经济学奖。他将投资组合的价格变化作为随机变量、以其均值衡量收益、以其方差衡量风险，把投资组合中各种证券所占比例作为变量，从而求出收益一定，方差（风险）最小的投资组合。而有效边界就是在收益-风险约束条件下能够以最小的风险取得最大的收益的各种证券的集合。
关于投资组合有效边界的数学证明相对复杂，所以这里我们利用Python进行可视化的分析和验证。
目录：
一、程序包准备
二、选择股票代码并进行数据的获取、处理、查看和可视化
三、计算资产组合的预期收益、均值和协方差
四、用蒙特卡洛模拟产生大量随机组合
五、投资组合优化1——Sharpe最大
六、投资组合优化2——方差最小
七、投资组合的有效前沿
程序包准备
我们使用开源财经数据接口包Tushare来获取相关的股票数据。如果没有该数据包，则通过pip installtushare的代码进行下载安装。（该程序包的安装和详细应用可以看往期推送内容）。
代码：
import tushare as ts #开源财经数据接口包import pandas as pdimport numpy as npimport statsmodels.api as sm #统计计算import scipy.stats as scs #科学计算import matplotlib.pyplot as plt #绘图工具
选择股票代码并进行数据的获取、处理、查看和可视化
我们这里选择沪深300（hs300），中国平安（601318），众泰汽车（000980），京东方A（000725）四只股票。从Tushare收集四只股票从2016-01-01至2019-11-01时期内的收盘价，并按照时间顺序导入Dataframe对象中。
代码：
data=pd.DataFrame()data0=ts.get_hist_data('hs300','2016-01-01','2019-11-01')data0=data0['close'] #取收盘价data0=data0[::-1] #按日期由小到大data['hs300']=data0data1=ts.get_hist_data('601318','2016-01-01','2019-11-01')data1=data1['close'] #取收盘价data1=data1[::-1] #按日期由小到大data['601318']=data1data2=ts.get_hist_data('000980','2016-01-01','2019-11-01')data2=data2['close'] #取收盘价data2=data2[::-1] #按日期由小到大data['000980']=data2data3=ts.get_hist_data('000725','2016-01-01','2019-11-01')data3=data3['close'] #取收盘价data3=data3[::-1] #按日期由小到大data['000725']=data3
然后对数据进行清理、查看和可视化。
代码：
data = data. dropna () #清理空白数据data. head () #显示前五行的数据(data / data. ix[0]*100) . plot (figsize = (8,4)) #figsize是用来设置图像大小的
显示前五行数据的结果为：
数据可视化后的结果：
计算资产组合的预期收益、均值和协方差
计算投资组合的协方差是构建资产组合过程中的核心部分。运用pandas内置方法可以生成协方差矩阵。我们先来看看四只不同证券的均值和协方差。
代码：
returns = np.log (data/ data. shift(1))returns. mean()*252#计算均值，这里乘以一个系数可以让后面可视化的图形更显著、更容易观察returns. cov() #计算协方差returns. corr() #计算相关系数
关于均值、协方差和相关系数，得到的结果如下：
从上可见，各证券之间的相关系数其实有一些偏大，如果要做出更精准的投资组合，建议使用相关系数小的各个证券，此处仅为参考。
然后，我们给不同资产随机分配初始权重。假设不允许建立空头头寸，所有的权重系数要在0~1之间。
代码：
noa=4weights=np. random. random(noa)weights/=np. sum(weights)weights
得到的结果为 array ([0.24034142, 0.2274179 , 0.2840282 , 0.24821248])
因为是随机分配权重，所以每一次操作结果都不一样。
从而我们计算预期组合的收益、组合方差和组合标准差。
代码：
np. sum(returns. mean()*weights)np. dot(weights.T,np.dot(returns. cov(), weights))np. sqrt(np.dot(weights.T,np.dot(returns. cov(), weights)))
得到的结果为：
用蒙特卡洛模拟产生大量随机组合
给定了一个投资组合，如何找到风险和收益平衡的位置呢？下面将通过一次蒙特卡洛模拟产生大量随机的权重，并记录随机组合的预期收益和方差。
代码：
port_returns=[]port_variance=[]for p in range(4000):weights=np.random.random(noa)weights/=np.sum(weights)port_returns.append(np.sum(returns.mean()*252*weights))port_variance.append(np.sqrt(np.dot(weights.T,np.dot(returns.cov()*252,weights))))
这里我们假设无风险利率为1.5%，从而可以画出大量随机组合下预期收益和方差的散点图。
代码：
port_returns=np.array(port_returns)port_variance=np.array(port_variance)risk_free=0.015plt.figure(figsize=(10,5)) #设置图表大小plt.scatter(port_variance,port_returns,c=(port_returns-risk_free)/port_variance,marker='o')plt.grid(True)plt.xlabel('expected volatility')plt.ylabel('expected return')plt.colorbar(label='Sharpe ratio')
得到的结果如下：
投资组合优化1——Sharpe最大
在众多可选择的投资组合中，我们应该尽可能地进行优化和筛选。最优化投资组合其实就是一个约束最优化的问题。
首先，我们先来了解一下夏普比率，夏普比率是一个可以同时对收益与风险加以综合考虑的三大经典指标之一，我们在进行选择投资组合时，作为理性投资人，应该“在固定所能承受的风险下，追求最大的报酬；或在固定的预期报酬下，追求最低的风险”。即选择SharpRatio最大的投资组合。
此处，我们建立statistics函数来记录投资组合的收益、方差和夏普比率。
代码：
def statistics(weights):weights=np.array(weights)port_returns=np.sum(returns.mean()*weights)*252port_variance=np.sqrt(np.dot(weights.T,np.dot(returns.cov()*252,weights)))returnnp.array([port_returns,port_variance,port_returns/port_variance])import scipy.optimize as scodef min_sharpe(weights):return-statistics(weights)[2] #最小化夏普比率的负值cons=({'type':'eq','fun':lambda x:np.sum(x)-1})bnds=tuple((0,1) for x in range(noa)) #权重在0-1之间opts=sco.minimize(min_sharpe,noa*[1./noa,],method='SLSQP',bounds=bnds,constraints=cons) #这里简单地使用平均分布optsopts['x'].round(3) #可以得到最优组合权重向量statistics(opts['x']).round(3) #可以得到sharpe最大的组合3个统计数据
结果：
投资组合优化2——方差最小
下面我们将通过方差最小来选出最优投资组合。
代码：
def min_variance(weights):returnstatistics(weights)[1]optv=sco.minimize(min_variance,noa*[1./noa,],method='SLSQP',bounds=bnds,constraints=cons)optvopts['x'].round(3) #方差最小的最优组合权重向量及组合的统计数据statistics(opts['x']).round(3) #得到预期收益率、波动率和夏普指数
得到结果如下
投资组合的有效前沿
有效边界由既定的目标收益率下方差最小的投资组合构成。在最优化时采用两个约束：（1）给定目标收益率（2）投资组合权重和为1
代码：
def min_variance(weights):returnstatistics(weights)[1]target_returns=np.linspace(0.0,0.5,50)target_variance=[]for tar in target_returns:cons=({'type':'eq','fun':lambdax:statistics(x)[0]-tar},{'type':'eq','fun':lambda x:np.sum(x)-1})res=sco.minimize(min_variance,noa*[1./noa,],method='SLSQP',bounds=bnds,constraints=cons)target_variance.append(res['fun'])target_variance=np.array(target_variance) #注意：这句代码在循环之外，切勿放在循环内！否则会造成多次的循环和报错。
下面是最优化结果的展示。
plt.figure(figsize=(12,6))plt.scatter(port_variance,port_returns,c=port_returns/port_variance,marker='o')#圆圈o：蒙特卡洛模拟随机产生的组合plt.scatter(target_variance,target_returns,c=target_returns/target_variance,marker='x')#叉号x：有效边界plt.plot(statistics(opts['x'])[1],statistics(opts['x'])[0],'r*',markersize=15.0)#红星：标记最高SHARPE组合plt.plot(statistics(optv['x'])[1],statistics(optv['x'])[0],'y*',markersize=15.0)#黄星：标记最小方差组合plt.grid(True)plt.xlabel('expected volatility')plt.ylabel('expected return')plt.colorbar(label='Sharpe ratio')
最终得到有效前沿图像（如果使用更多的相关性不大的证券作为投资组合，拟合效果会更好）：
作者：李嘉楠
编辑校对：郭通