在上一節(jié)課當中主要講了幾個重要的概念麻裁，包括股票的日回報率断国，日回報率的均值，累積回報率鱼鸠，波動(也就是標準差)猛拴，進而以日回報率的均值和波動來計算夏普比率，并說明夏普比率可以用來衡量每份風險背后所具有的收益蚀狰，以此來定量衡量我們投資決策的真實收益大小愉昆。

在這節(jié)課，我們要解決第二個問題麻蹋，給定一定的股票跛溉，比如說我給你四只股票的一系列的時間序列數(shù)據(jù)，以及給你限定10000元的投資，那么如何根據(jù)所給的股票數(shù)據(jù)來劃分每只股票各投資多少占比以求得最大收益芳室，也就是我要計算出股票購買所占資本的權重专肪，這個權重可以實現(xiàn)總體的夏普比率最大化。
總體大綱：

流程

一堪侯、加載數(shù)據(jù)及預處理

1.1 計算累積收益率嚎尤，對數(shù)收益率和算術收益率

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
AAPL=pd.read_csv('AAPL_CLOSE',index_col='Date',parse_dates=True)
AMZN=pd.read_csv('AMZN_CLOSE',index_col='Date',parse_dates=True)
CISCO=pd.read_csv('CISCO_CLOSE',index_col='Date',parse_dates=True)
IBM=pd.read_csv('IBM_CLOSE',index_col='Date',parse_dates=True)
AAPL.head()

輸出結果

接下來我們把這四只股票的收盤價都合在一張dataframe里面

stock=pd.concat([AAPL, AMZN, CISCO, IBM],axis=1)
stock.columns=['AAPL', 'AMZN', 'CISCO', 'IBM']
stock.head()

輸出結果

然后我們計算一下日收益率的平均值

stock.pct_change(1).mean()

輸出結果

我們接下來看一下日收益率的偏相關性所引起的偏方差攀涵。關于這里的偏相關性我是這么理解的画株，因為這是股票市場，并且這幾只股票都是科技股票婶博，它們的市場是趨同的雹拄，因此在這里應該是認為彼此的價格變化是有相互影響的收奔，因此在計算波動的時候，要考慮權重轉變后其他股票對單只股票的影響滓玖。偏方差函數(shù)是cov()

stock.pct_change(1).cov()

在這里以第一行為例子坪哄，AAPL是蘋果股票的代號，它除了有自身波動的影響势篡，即第一列翩肌。也有其他三只股票的影響，也就是其他三列禁悠，這部分的相關性導致的波動即為偏方差念祭。

我們繪制一下股票的累積收益率的變化情況。

normed_ret=stock/stock.iloc[0]
normed_ret.plot(kind='line',figsize=(12,8),grid=True)
plt.legend(loc='best')

我們在這里引入對數(shù)收益率碍侦，原因是對數(shù)收益率相對于我們上面提到的算術收益率來說更加適合模型的計算粱坤，具體原因我覺得知乎這個回答是不錯的，一方面是為了讓我們的數(shù)據(jù)符合模型的平穩(wěn)性假設瓷产，另外一方面應該是方便計算表達站玄，而且可以發(fā)現(xiàn)對數(shù)收益率和算術收益率其實差的不多。關于對數(shù)收益率我還有些不理解的地方濒旦，還得繼續(xù)研究下株旷。我們繪制下對數(shù)收益率和算術收益率的分布吧。

#算術收益率
stock.pct_change(1).plot.hist(bins=100)
plt.legend(loc='best')
#對數(shù)收益率
log_ret=np.log(stock/stock.shift(1))
log_ret.plot.hist(bins=100)
plt.legend(loc='best')

算術收益率的分布

對數(shù)收益率的分布

可以看到差別并不是很大尔邓。

1.2 計算252天的平均收益率的總和以及偏方差矩陣

#252天的四只股票的平均對數(shù)收益率總和
total_mean_ret=log_ret.mean()*252
#252天的四只股票構成的偏方差矩陣
total_covariance=log_ret.cov()*252
print total_mean_ret
print total_covariance

在這里計算資產(chǎn)組合的最優(yōu)化權重晾剖，我們可以用到兩個方法：一個是蒙特卡羅模擬法，另外一個是限制條件下的數(shù)學優(yōu)化方法铃拇。我們分別來看一下钞瀑。

二、蒙特卡羅模擬法

2.1 原理

這個方法聽起來很高大上慷荔，其實方法很接地氣，也很暴力。就是舉出很多很多權重組合（多到幾乎窮舉）显晶，然后計算所有的點的夏普比率贷岸，哪個最高就取哪個點的均值和波動。

2.2 方法實現(xiàn)

2.2.1 計算單個權重的夏普比率

在這里我們先用代碼實現(xiàn)一組投資資產(chǎn)權重的計算夏普比率磷雇，窮舉只用在單次的外圍加上一個遍歷就可以了偿警。

numpy的random函數(shù)用法

#在這里我們用numpy的random.randn函數(shù)實現(xiàn)隨機取四個數(shù)，這個方法是從[0, 1]中取出特定個數(shù)的數(shù)字的方法
weight=np.random.random(4)
print weight

在這里我們看到的確是返回了四個隨機的[0, 1]的值唯笙，但是要注意我們的權重之和是需要等于1的螟蒸，所以要進行歸一化。

normed_weight=weight/np.sum(weight)
print normed_weight

可以看到歸一化以后的權重確實是之和為1了崩掘。然后我們試下用這個權重去計算對應的資產(chǎn)組合的夏普比率七嫌。
按照權重去計算對應的252天的預期回報率和波動（也就是方差）

#預期回報率
exp_ret=np.sum(log_ret.mean()*normed_weight*252)
#波動
volat=np.sqrt(np.dot(normed_weight.T, np.dot(log_ret.cov()*252, normed_weight)))
print exp_ret, volat
#計算夏普比率
SR=exp_ret/volat
print SR

夏普比率

在這里關于波動的計算方式，我是這么理解的后面的一個點乘實際上是計算出每個股票單獨的方差苞慢，然后每只股票的方差又按照權重加和得到最終整體的方差诵原。原因是每只股票波動是跟其他三只股票都是有聯(lián)系的，所以要按照權重先算出每只股票自身的波動挽放。第一步得到的是41的矩陣绍赛，然后前面的weights在轉置以后是14的矩陣，兩者的點乘就是最終經(jīng)過權重轉化的方差辑畦，開方后為標準差吗蚌。

2.2.3 生成多個隨機權重進行模擬

接下來我們就在外面套上一層循環(huán)，生成許多組不同的權重序列纯出，按照上面的方式計算每組權重所得的期望回報和波動蚯妇。

#隨機權重的組數(shù)為15000
num=15000
#初始化總預期收益率和波動以及夏普比率的numpy數(shù)列
exp_rets=np.zeros(num)
exp_vols=np.zeros(num)
exp_SRs=np.zeros(num)
#迭代15000次
for ind in range(num):
    weight=np.random.random(4)
    normed_weight=weight/np.sum(weight)
    #預期回報率
    exp_ret=np.sum(log_ret.mean()*normed_weight*252)
    exp_rets[ind]=exp_ret
    #波動
    volat=np.sqrt(np.dot(normed_weight.T, np.dot(log_ret.cov()*252, normed_weight)))
    exp_vols[ind]=volat
    #計算夏普比率
    SR=exp_ret/volat
    exp_SRs[ind]=SR
#我們找出Sharp Ratio對應最大的Volatility和return那個點，可以用numpy的argmax函數(shù)
max_ind=np.argmax(exp_SRs)
max_SR_vol=exp_vols[max_ind]
max_SR_ret=exp_rets[max_ind]
print 'Sharp Ratio maximum is ', exp_SRs[max_ind]
print '-'*20
print 'Volatility is ', exp_vols[max_ind]
print '-'*20
print 'Return is ', exp_rets[max_ind]

#我們以波動為橫坐標潦刃，總收益率為縱坐標侮措，夏普比率為顏色條，繪制出散點圖的分布乖杠，并且標識出最大夏普比率的那個點
plt.scatter(x=exp_vols, y=exp_rets, c=exp_SRs, alpha=0.5, cmap='coolwarm')
plt.title('Volatility versus Return')
plt.xlabel('Volatility')
plt.ylabel('Return')
plt.colorbar(label='Sharp Ratio')
plt.scatter(max_SR_vol, max_SR_ret, s=50, edgecolors='black')

三個值的輸出結果

可視化模擬結果

在這里我們標識出了那個夏普比率最大的那個點分扎。

三、數(shù)學優(yōu)化方法

上面是通過隨機模擬的方式來找出最優(yōu)的投資配置胧洒，但實際上我們也可以通過邊界條件下的優(yōu)化來找出最優(yōu)的權重畏吓。這里要用到scipy.optimize.minimize方法。

scipy.optimize.minimize方法

這里我們可以看到這個函數(shù)要做的事情就是給定邊界條件

然后計算這個一個或多個變量所決定的標量目標函數(shù)f(x)的最小值卫漫。

從官方文檔可以看到有幾個重要參數(shù)菲饼，我們在這里只傳入幾個參數(shù)：fun最小化的目標函數(shù)，x0是一開始的初始值列赎，method是用于優(yōu)化的方法宏悦，這里選用‘SLSQP’，constraints是約束條件，整體是一個字典饼煞，type有兩種：‘eq’和'ineq'源葫，分別表示后面?zhèn)魅氲募s束條件函數(shù)'fun'等于0或者大于等于0。
在這里要注意兩點砖瞧，一個是我們是要最大化夏普比率的息堂，所以為了用這個minimize方法，我們應該定義目標函數(shù)為夏普比率的負值块促；第二個是約束條件函數(shù)在這里應該是權重之和也就是自變量之和為1荣堰，要稍微調(diào)整成權重之和減1等于0的約束條件形式才滿足minimize方法的要求。

#以權重作為輸入竭翠，計算所得252天的總收益率振坚，總的波動，夏普比率
def get_ret_vol_SR(weight):
    #預期回報率
    exp_ret=np.sum(log_ret.mean()*weight*252)
    #波動
    volat=np.sqrt(np.dot(weight.T, np.dot(log_ret.cov()*252, weight)))
    #計算夏普比率
    SR=exp_ret/volat
    return np.array([exp_ret,  volat, SR])
#定義目標函數(shù)逃片，返回夏普比率的負值
def get_neg_sharp_ratio(weight):
    return get_ret_vol_SR(weight)[2]*-1
#定義限制的邊界函數(shù)屡拨，返回權重之和減1的計算結果，之后會用來跟0作比較褥实，判斷是否相等
def con(weight):
    return np.sum(weight)-1
#初始化四個值為0.25
x0=np.array([0.25,0.25,0.25,0.25])
#設置自變量的取值范圍
bounds=((0,1),(0,1),(0,1),(0,1))
#設置限制條件的字典
constraints={'type':'eq','fun':con }
#進行邊界條件下對目標函數(shù)優(yōu)化的權重計算
from scipy.optimize import minimize
res=minimize(fun=get_neg_sharp_ratio, x0=x0, constraints=constraints, method='SLSQP',bounds=bounds)
print res.x#打印出最優(yōu)化夏普比率后的權重

然后我們接下來根據(jù)計算得到的權重呀狼，來帶入到我們計算三個值的函數(shù)當中求出總的收益率，總的波動损离，以及夏普比率哥艇。

print 'Epxcted return is', get_ret_vol_SR(res.x)[0]
print '-'*20
print 'Epxcted volatility is', get_ret_vol_SR(res.x)[1]
print '-'*20
print 'Epxcted return is', get_ret_vol_SR(res.x)[2]

四、有效邊界

有效邊界可以這么認為僻澎，是一條在給定波動也就是風險下貌踏，所能達到的最大收益的那個點，這種類似的點的組合即為有效邊界窟勃。我們繪制一下這條有效邊界祖乳。

#numpy的linspace方法可以生成一系列梯度漸進的值
lin_y=np.linspace(0,0.3,100)
#然后我們在這里定義的最小化的目標函數(shù)是波動，也就是在確定的return也就是y值以及權重之和為1的限制條件下秉氧，最小化波動值眷昆。
def get_vol(weight):
    return get_ret_vol_SR(weight)[1]
frontier_volatility=[]
#初始化四個值為0.25
x0=np.array([0.25,0.25,0.25,0.25])
#設置自變量的取值范圍
bounds=((0,1),(0,1),(0,1),(0,1))
for y in lin_y:
    constraints2=[{'type':'eq','fun':lambda x:get_ret_vol_SR(x)[1]-y},{'type':'eq','fun':con}]
    res2=minimize(fun=get_vol, x0=x0, constraints=constraints2, method='SLSQP',bounds=bounds)
    frontier_volatility.append(res2.fun)
print frontier_volatility[:10]

輸出結果

我們繪制一下邊界曲線

#我們以波動為橫坐標，總收益率為縱坐標汁咏，夏普比率為顏色條亚斋，繪制出散點圖的分布，并且標識出最大夏普比率的那個點
plt.scatter(x=exp_vols, y=exp_rets, c=exp_SRs, alpha=0.5, cmap='coolwarm')
plt.title('Volatility versus Return')
plt.xlabel('Volatility')
plt.ylabel('Return')
plt.colorbar(label='Sharp Ratio')
plt.scatter(max_SR_vol, max_SR_ret, s=50, edgecolors='black')
plt.plot(frontier_volatility, lin_y, 'g--', linewidth=1.5)

有效邊界

在這里如果從縱向上看攘滩，相同的波動情況下帅刊，有效邊界上面的點能夠達到最大的收益，相同的收益下漂问，有效邊界上面的點能夠達到風險最小赖瞒，因此有效邊界上面的點都是投資配置最優(yōu)化的點女揭。