代码如下
# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import pandas.io.data as web
import matplotlib.pyplot as plt

symbols=['AAPL','MSFT','YHOO','DB','GLD']
noa=len(symbols)
data=pd.DataFrame()
for sym in symbols:
    data[sym]=web.DataReader(sym,data_source='yahoo',end='2014-09-12')['Adj Close']
data.columns=symbols
(data / data.ix[0] *100).plot(figsize=(8,5))
rets=np.log(data/data.shift(1))          #计算收益
print(rets.mean()*252)
print(rets.cov()*252)
init_weights=np.random.random(noa)
init_weights=init_weights/np.sum(weights)        #随机生成的总和为1的权重
print(np.dot(rets.mean().T,init_weights)*252) #组合的收益！
print(np.sqrt(np.dot(weights.T,np.dot(rets.cov()*252,weights)))) #组合的波动率！
#----------------------特别容易忘记收益率，协方差矩阵乘上252！！！----------
#接下来就是蒙特卡洛模拟打点
prets=[]
pvols=[]
for p in range(10000):
    weights=np.random.random(noa)
    weights /= np.sum(weights)
    prets.append(np.sum(rets.mean()*weights*252))
    pvols.append(np.sqrt(np.dot(weights.T,np.dot(rets.cov()*252,weights))))
prets=np.array(prets)          #这里用了np.array主要是因为不用的话就是list
pvols=np.array(pvols)           #list没有array方便计算！！！
plt.figure(figsize=(8,4))
risk_free=0.02


#-------------------------------------------------------------------------

def statistics(weights):
    weights=np.array(weights)
    pret=np.sum(rets.mean()*weights)*252
    pvol=np.dot(weights.T,np.dot(rets.cov()*252,weights))
    return np.array([pret,pvol,pret/pvol])
def min_func_sharp(weights):
    return -statistics(weights)[2]
cons1 = ({'type':'eq', 'fun':lambda x: np.sum(x)-1})
bnds1=tuple((0,1) for x in range(noa))
opts = sco.minimize(min_func_sharp, noa*[1./noa,], method = 'SLSQP', bounds = bnds, constraints = cons1)
print(opts)
print('权重是：',opts['x'].round(3))
print('各种值是：',statistics(opts['x']).round(3))
def min_func_variance(weights):
    return statistics(weights)[1]**2
optv=sco.minimize(min_func_variance,noa*[1./noa,],method='SLSQP',bounds=bnds1,constraints=cons)
print(optv)
print('权重是：',optv['x'].round(3))
print('各种值是：',statistics(optv['x']))
#-----------------------------------------------------------------------
def min_variance(weights):
    return statistics(weights)[1]
target_returns = np.linspace(0.0,0.25,50)
target_variance = []
for tar in target_returns:
    cons2 = ({'type':'eq','fun':lambda weights:statistics(weights)[0]-tar},{'type':'eq','fun':lambda weights:np.sum(weights)-1})
    res = sco.minimize(min_variance, noa*[1./noa,],method = 'SLSQP', bounds = bnds1, constraints = cons2)
    target_variance.append(res['fun'])
target_variance = np.array(target_variance)
plt.scatter(pvols,prets,c=prets-risk_free / pvols, marker='.')  #'c'这个参数是颜色，可以用小数表示
                                                        #这里转换为了夏普比率
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label='Sharp ratio')
plt.scatter(target_variance,target_returns, c = target_returns/target_variance, marker = 'x')

问题是：