如何用Python找到投資時的最佳組合比例
關注公眾號，回覆“20190327”，即可得到本教程程式碼
現代投資組合理論（Modern Portfolio Theory，MPT）告訴我們投資者應該分散投資來實現最小化風險最大化投資回報。大鄧剛開始學習這方面知識，用了將近一天的時候才搞懂MPT理論的推導，順便複習了部分高中數學知識，這樣會讓我們更加有新信心的去使用自己編寫的程式碼。現在我們從實戰開始接觸理論。
一、資產組合理論導論
1.1 風險厭惡（Risk aversion）
在投資組合理論中，我們常常使用方差來刻畫資產的風險。這裡舉個例子，方便大家理解。
假設你現在將要進行投資，有兩種策略：
- 資產A給你帶來的收益是200元，或者0元。每種情況發生的概率也是各50%
- 資產B給你帶來的收益是400元，或者-200元（虧損）。每種情況發生的概率也是各50%
資產的數學期望收益
兩種投資策略的波動性-標準差
如果你是風險偏好者，你可能會選擇B，因為B的潛在的最大收益最大。但是MPT理論認為大多數的投資者是風險厭惡者，不喜歡玩心跳，所以更傾向於選擇A。
1.2資產組合池(portfolio)
假設我們將採用分散投資，每種資產的比例為 ``w1、w2、w3...wn``，我們知道所有投資比例之和為1，即 ``w1+w2+w3+...+wn=1`` 。假設R0、R1、R2...Rn分別代表每種資產的收益，則資產組合投資收益為
我們預期的投資組合收益為
1.3 相關性
在計算資產組合風險前，我們需要先回憶一下高中數學中的方差和相關性的計算方法。方差和相關性主要是刻畫任意兩個變數之間的線性關係。
X和Y的方差計算公式
1.4 風險
現在我們將要計算資產組合的風險，這裡使用 ``資產組合收益的方差`` 來刻畫投資風險。
本來大鄧直接看到推導完的結果，但是忽略了中間過程，心裡怎麼也不相信。所以花了很多時間用來找推導過程的教程和視訊，終於找到一份比較詳細的，現在貼給大家看。我們將該公式中的加總用累加zigma符號表示，一步步的進行推導
現在我們先看看簡單的例子，如何將累加和的平方展開？
我們再將上面的公式抽象一下
因此，資產組合收益的方差可以表示為
最後我們將上面的符號全部用向量（或矩陣）來表示，這樣後期我們使用python寫程式碼時候可以直接使用numpy表達。
二、實戰
2.1 載入資料
這裡我參考JoinQuant中的一篇文章，其選取的股票的收益比較符合正太分佈，我就直接拿來用吧。
'000651',##格力電器 '600519',##貴州茅臺 '601318',##中國平安 '000858',##五糧液 '600887',##伊利股份 '000333',##美的集團 '601166',##興業銀行 '601328',##交通銀行 '600104'##上汽集團
我們使用tushare來獲取股票資料，tushare的get_hist_data(stock_code, start, end)函式獲取stock_code從start到end期間內的所有交易資料，返回的是dataframe型別。這裡我們都是用close列的資料，即只用收盤價資料。
import tushare as ts
def fetch_stock_data(stock_code, stock_name, start, end):
df = ts.get_hist_data(stock_code, start=start, end=end) #前復權
df = df.close
df.name = stock_name
return df
geli = fetch_stock_data('000651', '格力電器', '2016-05-28', '2018-03-26')
geli.head()
Run
date 2018-03-2648.13 2018-03-2348.88 2018-03-2249.90 2018-03-2151.70 2018-03-2052.20 Name: 格力電器, dtype: float64
現在我們將所有股票都的收盤資料都裝進pandas.DataFrame中，每一列代表一隻股票，每一行代表一天的收盤價
import pandas as pd
data = pd.DataFrame({'格力電器':fetch_stock_data('000651', '格力電器', '2016-05-28', '2018-03-26'),
'貴州茅臺':fetch_stock_data('600519', '貴州茅臺', '2016-05-28', '2018-03-26'),
'中國平安':fetch_stock_data('601318', '中國平安', '2016-05-28', '2018-03-26'),
'五糧液':fetch_stock_data('000858', '五糧液', '2016-05-28', '2018-03-26'),
'伊利股份':fetch_stock_data('600887', '伊利股份', '2016-05-28', '2018-03-26'),
'美的集團':fetch_stock_data('000333', '美的集團', '2016-05-28', '2018-03-26'),
'興業銀行':fetch_stock_data('601166', '興業銀行', '2016-05-28', '2018-03-26'),
'交通銀行':fetch_stock_data('601328', '交通銀行', '2016-05-28', '2018-03-26'),
'上汽集團':fetch_stock_data('600104', '上汽集團', '2016-05-28', '2018-03-26')})
data = data.dropna()
data.head()
data.to_excel('stock_data.xlsx')
我將整理好的資料儲存到了 ``stock_data.xlsx`` 中，方便大家即使無法使用tushare也有資料可供分析。
import pandas as pd
data = pd.read_excel('stock_data.xlsx')
data.head()
將每個股票價格與最初始（2016年10月24日）的價格作比較，並據此得到之後的股價走勢圖
#將date列從newdata中踢出
date = data.pop('date')
#data.iloc[0, :] 選取第一行的資料
newdata = (data/data.iloc[0, :])*100
昨天我們分享的視覺化使用的pyplotz庫，該庫支援所有基於matplotlib的視覺化庫。今天我們用plotly這個動態視覺化庫來作圖，正好複習下前面分享的plotly-express庫，注意plotly-express是基於plotly
from plotly.offline import init_notebook_mode, iplot, plot import plotly.graph_objs as go init_notebook_mode() stocks = ['格力電器', '貴州茅臺', '中國平安', '五糧液', '伊利股份', '美的集團', '興業銀行', '交通銀行', '上汽集團'] def trace(df, date, stock): return go.Scatter(x = date, #橫座標日期 y = df[stock], name=stock)#縱座標為股價與（2016年10月24日）的比值 data = [trace(newdata,date,stock) for stock in stocks] iplot(data)
2.2 計算不同股票的均值、協方差
每年有252個交易日，用每日收益率乘以252得到年華收益率。現在需要計算每隻股票的收益率，在金融領域中我們一般使用對數收益率。這裡體現了pandas的強大，df.pct_change()直接就能得到股票收益率
import numpy as np log_returns = np.log(newdata.pct_change()+1) log_returns = log_returns.dropna() log_returns.mean()*252
Run
格力電器0.582497 貴州茅臺0.648666 中國平安0.524031 五糧液0.561282 伊利股份0.352049 美的集團0.560136 興業銀行0.024539 交通銀行0.068539 上汽集團0.289760 dtype: float64
2.3 進行正態檢驗
馬科維茨的投資組合理論需要滿足收益率符合正態分佈，scipy.stats庫為我們提供了正態性測試函式
scipy.stats.normaltest 測試樣本是否與正態分佈不同，返回p值。
import scipy.stats as scs
def normality_test(array):
print('Norm test p-value %14.3f' % scs.normaltest(array)[1])
for stock in stocks:
print('\nResults for {}'.format(stock))
print('-'*32)
log_data = np.array(log_returns[stock])
normality_test(log_data)
Run
Results for 格力電器 -------------------------------- Norm test p-value0.000 Results for 貴州茅臺 -------------------------------- Norm test p-value0.000 Results for 中國平安 -------------------------------- Norm test p-value0.000 Results for 五糧液 -------------------------------- Norm test p-value0.051 Results for 伊利股份 -------------------------------- Norm test p-value0.000 Results for 美的集團 -------------------------------- Norm test p-value0.453 Results for 興業銀行 -------------------------------- Norm test p-value0.000 Results for 交通銀行 -------------------------------- Norm test p-value0.000 Results for 上汽集團 -------------------------------- Norm test p-value0.000
從上面的檢驗中，伊利股份和興業銀行股票收益率不符合正態分佈。
2.4 投資組合預期收益率、波動率
我們先隨機生成一維投資組合權重向量（長度為9，與股票數量相等），因為中國股市的不允許賣空，所以投資組合權重向量中的數值必須在0到1之間。
weights = np.random.random(9) weights /= np.sum(weights) weights
array([0.13403144, 0.11703939, 0.14125659, 0.02530677, 0.10496042, 0.16106127, 0.05155371, 0.10361131, 0.16117911])
投資組合預期收益率等於每隻股票的權重與其對應股票的年化收益率的乘積。
np.dot(weights, log_returns.mean())*252
0.4035947984701051
投資組合波動率（方差）
np.dot(weights, np.dot(log_returns.cov()*252, weights))
0.035798322938178584
投資組合收益的年化風險（標準差）
np.sqrt(np.dot(weights, np.dot(log_returns.cov()*252, weights)))
0.18920444745877033
2.5 隨機生成大量的投資組合權重
生成1000種隨機的投資組合，即權重weights的尺寸為（1000\*9）。
import matplotlib.pyplot as plt
%matplotlib inline
port_returns = []
port_variance = []
for p in range(1000):
weights = np.random.random(9)
weights /=np.sum(weights)
port_returns.append(np.sum(log_returns.mean()*252*weights))
port_variance.append(np.sqrt(np.dot(weights.T, np.dot(log_returns.cov()*252, weights))))
port_returns = np.array(port_returns)
port_variance = np.array(port_variance)
#無風險利率設定為3%
risk_free = 0.03
plt.figure(figsize=(8, 6))
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')
2.5.1 投資組合優化1—夏普率最大
建立statss函式來記錄重要的投資組合統計資料（收益，方差和夏普比）。scipy.optimize可以提供給我們最小優化演算法，而最大化夏普率可以轉化為最小化負的夏普率。
import scipy.optimize as sco
def stats(weights):
weights = np.array(weights)
port_returns = np.sum(log_returns.mean()*weights)*252
port_variance = np.sqrt(np.dot(weights.T, np.dot(log_returns.cov()*252,weights)))
return np.array([port_returns, port_variance, port_returns/port_variance])
#最小化夏普指數的負值
def min_sharpe(weights):
return -stats(weights)[2]
#給定初始權重
x0 = 9*[1./9]
#權重（某股票持倉比例）限制在0和1之間。
bnds = tuple((0,1) for x in range(9))
#權重（股票持倉比例）的總和為1。
cons = ({'type':'eq', 'fun':lambda x: np.sum(x)-1})
#優化函式呼叫中忽略的唯一輸入是起始引數列表(對權重的初始猜測)。我們簡單的使用平均分佈。
opts = sco.minimize(min_sharpe,
x0,
method = 'SLSQP',
bounds = bnds,
constraints = cons)
opts
Run
fun: -2.7277947674404794 jac: array([ 1.08440101e-01, -4.65214252e-05,3.44902277e-04,7.20474541e-01, 5.15412062e-01, -2.46465206e-05,3.42730999e-01, -1.48534775e-04, -5.43147326e-04]) message: 'Optimization terminated successfully.' nfev: 111 nit: 10 njev: 10 status: 0 success: True x: array([1.07588996e-16, 4.93897426e-01, 2.37878143e-01, 6.47455750e-17, 5.74725095e-17, 1.14655596e-01, 8.60056411e-17, 6.88721816e-02, 8.46966532e-02])
最優投資組合權重向量，小數點保留3位
opts['x'].round(3)
array([0., 0.494, 0.238, 0., 0., 0.115, 0., 0.069, 0.085])
sharpe最大的組合3個統計資料分別為：
stats(opts['x']).round(3) array([0.539, 0.197, 2.728])
2.5.2 投資組合優化2——方差最小
接下來，我們通過方差最小來選出最優投資組合。
#但是我們定義一個函式對 方差進行最小化 def min_variance(weights): return stats(weights)[1] optv = sco.minimize(min_variance, x0, method = 'SLSQP', bounds = bnds, constraints = cons) optv
Run
fun: 0.1260429626044057 jac: array([0.15127511, 0.12627605, 0.1532943 , 0.14743594, 0.12581278, 0.1258321 , 0.12612321, 0.1258938 , 0.12583541]) message: 'Optimization terminated successfully.' nfev: 88 nit: 8 njev: 8 status: 0 success: True x: array([0.00000000e+00, 2.13881938e-01, 0.00000000e+00, 8.07576429e-18, 3.06402571e-02, 1.40860179e-02, 3.22191709e-01, 3.65127109e-01, 5.40729695e-02])
方差最小的最優投資組合權重向量
optv['x'].round(3)
array([0., 0.214, 0., 0., 0.031, 0.014, 0.322, 0.365, 0.054])
得到的投資組合預期收益率、波動率和夏普指數
stats(optv['x']).round(3)
array([0.206, 0.126, 1.634])
2.5.3 組合的有效邊界
有效邊界是由一系列既定的目標收益率下方差最小的投資組合點組成的。在最優化時採用兩個約束，1.給定目標收益率，2.投資組合權重和為1。
def min_variance(weights):
return statistics(weights)[1]
#在不同目標收益率水平（target_returns）迴圈時，最小化的一個約束條件會變化。
target_returns = np.linspace(0.0,0.5,50)
target_variance = []
for tar in target_returns:
#給定限制條件：給定收益率、投資組合權重之和為1
cons = ({'type':'eq','fun':lambda x:stats(x)[0]-tar},{'type':'eq','fun':lambda x:np.sum(x)-1})
res = sco.minimize(min_variance, x0, method = 'SLSQP', bounds = bnds, constraints = cons)
target_variance.append(res['fun'])
target_variance = np.array(target_variance)
下面是最優化結果的展示。
叉號：構成的曲線是有效前沿（目標收益率下最優的投資組合）
紅星：sharpe最大的投資組合
黃星：方差最小的投資組合
plt.figure(figsize = (8,4))
#圓點：隨機生成的投資組合散佈的點
plt.scatter(port_variance, port_returns, c = port_returns/port_variance,marker = 'o')
#叉號：投資組合有效邊界
plt.scatter(target_variance,target_returns, c = target_returns/target_variance, marker = 'x')
#紅星：標記夏普率最大的組合點
plt.plot(stats(opts['x'])[1], stats(opts['x'])[0], 'r*', markersize = 15.0)
#黃星：標記方差最小投資組合點
plt.plot(stats(optv['x'])[1], stats(optv['x'])[0], 'y*', markersize = 15.0)
plt.grid(True)
plt.xlabel('expected volatility')
plt.ylabel('expected return')
plt.colorbar(label = 'Sharpe ratio')
Run
從黃色五角星到紅色五角星是投資最有效的組合，這一系列的點所組成的邊界就叫做 投資有效邊界 。 這條邊界的特點是同樣的風險的情況下獲得的收益最大，同樣的收益水平風險是最小的。從這條邊界也印證了風險與收益成正比，要想更高的收益率就請承擔更大的風險，但如果落在投資有效邊界上，價效比最高。