用 Python 中的蒙特卡洛模擬預(yù)測股票收益
蒙特卡洛方法（或蒙特卡洛實(shí)驗(yàn)）是一大類計(jì)算算法，它們依賴于重復(fù)隨機(jī)采樣來獲得數(shù)值結(jié)果?；舅枷胧鞘褂秒S機(jī)性來解決原則上可能是確定性的問題。它們通常用于物理和數(shù)學(xué)問題，并且在難以或不可能使用其他方法時(shí)最有用。Monte Carlo 方法主要用于三個(gè)不同的問題類別：優(yōu)化、數(shù)值積分和從概率分布中生成繪圖。
我們將使用蒙特卡洛模擬來觀察資產(chǎn)價(jià)格隨時(shí)間的潛在變化，假設(shè)它們的每日收益服從正態(tài)分布。這種類型的價(jià)格演變也稱為“隨機(jī)游走”（random walk）。
例如，如果我們想購買一只特定的股票，我們可能想嘗試展望未來并預(yù)測可以以何種概率期望得到何種回報(bào)，或者我們可能有興趣調(diào)查哪些潛在的極端結(jié)果我們可能會(huì)經(jīng)歷以及我們面臨破產(chǎn)風(fēng)險(xiǎn)的程度，或者另一方面，獲得超額回報(bào)的風(fēng)險(xiǎn)。
為了建立模擬，我們需要估計(jì)相關(guān)股票的預(yù)期回報(bào)水平 (mu) 和波動(dòng)率 (vol)。這些數(shù)據(jù)可以從歷史價(jià)格中估計(jì)出來，最簡單的方法是假設(shè)過去的平均回報(bào)和波動(dòng)率水平將持續(xù)到未來。還可以調(diào)整歷史數(shù)據(jù)以考慮投資者觀點(diǎn)或市場制度變化等，但是為了保持簡單并專注于代碼，我們將根據(jù)過去的價(jià)格數(shù)據(jù)設(shè)置簡單的回報(bào)和波動(dòng)率水平。
現(xiàn)在讓我們開始編寫一些代碼并生成我們需要的初始數(shù)據(jù)作為我們蒙特卡羅模擬的輸入。出于說明目的，讓我們看看蘋果公司的股票......
首先我們必須導(dǎo)入必要的模塊，然后開始：
#import necessary packagesimport numpy as npimport mathimport matplotlib.pyplot as pltfrom scipy.stats import normfrom pandas_datareader import data#download Apple price data into DataFrameapple = data.DataReader('AAPL', 'yahoo',start='1/1/2000')#calculate the compound annual growth rate (CAGR) which#will give us our mean return input (mu)days = (apple.index[-1] - apple.index[0]).dayscagr = ((((apple['Adj Close'][-1]) / apple['Adj Close'][1])) ** (365.0/days)) - 1print ('CAGR =',str(round(cagr,4)*100)+"%")mu = cagr#create a series of percentage returns and calculate#the annual volatility of returnsapple['Returns'] = apple['Adj Close'].pct_change()vol = apple['Returns'].std()*sqrt(252)print ("Annual Volatility =",str(round(vol,4)*100)+"%")
結(jié)果如下：
CAGR = 23.09%Annual Volatility = 42.59%
現(xiàn)在我們知道我們的復(fù)合年均增長率是 23.09%，我們的波動(dòng)率輸入 (vol) 是 42.59%——實(shí)際運(yùn)行蒙特卡羅模擬的代碼如下：
#Define VariablesS = apple['Adj Close'][-1] #starting stock price (i.e. last available real stock price)T = 252 #Number of trading daysmu = 0.2309 #Returnvol = 0.4259 #Volatility#create list of daily returns using random normal distributiondaily_returns=np.random.normal((mu/T),vol/math.sqrt(T),T)+1#set starting price and create price series generated by above random daily returnsprice_list = [S]for x in daily_returns:price_list.append(price_list[-1]*x)#Generate Plots - price series and histogram of daily returnsplt.plot(price_list)plt.show()plt.hist(daily_returns-1, 100) #Note that we run the line plot and histogram separately, not simultaneously.plt.show()
此代碼輸出繪圖：
上面的代碼基本上運(yùn)行了一個(gè)交易年度（252 天）內(nèi)潛在價(jià)格序列演變的單一模擬，基于遵循正態(tài)分布的每日收益隨機(jī)的抽取。由第一個(gè)圖表中顯示的單線系列表示。第二個(gè)圖表繪制了一年期間這些隨機(jī)每日收益的直方圖。
現(xiàn)在我們已經(jīng)成功模擬了未來一年的每日價(jià)格數(shù)據(jù)。但實(shí)際上它并沒有讓我們深入了解股票的風(fēng)險(xiǎn)和回報(bào)特征，因?yàn)槲覀冎挥幸粋€(gè)隨機(jī)生成的路徑。實(shí)際價(jià)格完全按照上表所述演變的可能性幾乎為零。
那么你可能會(huì)問這個(gè)模擬有什么意義呢？好吧，真正的洞察力是通過運(yùn)行數(shù)千次、數(shù)萬次甚至數(shù)十萬次模擬獲得的，每次運(yùn)行都會(huì)根據(jù)相同的股票特征（mu 和 vol）產(chǎn)生一系列不同的潛在價(jià)格演變。
我們可以非常簡單地調(diào)整上面的代碼來運(yùn)行多個(gè)模擬。此代碼如下所示。在下面的代碼中，您會(huì)注意到一些事情——首先我刪除了直方圖（我們稍后會(huì)以稍微不同的方式回到這個(gè)問題），并且代碼現(xiàn)在在一個(gè)圖表上繪制多個(gè)價(jià)格系列以顯示信息對(duì)于每個(gè)單獨(dú)的模擬運(yùn)行。
import numpy as npimport mathimport matplotlib.pyplot as pltfrom scipy.stats import norm#Define VariablesS = apple['Adj Close'][-1] #starting stock price (i.e. last available real stock price)T = 252 #Number of trading daysmu = 0.2309 #Returnvol = 0.4259 #Volatility#choose number of runs to simulate - I have chosen 1000for i in range(1000):#create list of daily returns using random normal distributiondaily_returns=np.random.normal(mu/T,vol/math.sqrt(T),T)+1#set starting price and create price series generated by above random daily returnsprice_list = [S]for x in daily_returns:price_list.append(price_list[-1]*x)#plot data from each individual run which we will plot at the endplt.plot(price_list)#show the plot of multiple price series created aboveplt.show()
這為我們提供了以下 1000 個(gè)不同模擬價(jià)格系列的圖:
現(xiàn)在我們可以看到 1000 次不同模擬產(chǎn)生的潛在結(jié)果，考慮到每日收益序列的隨機(jī)性，所有模擬都基于相同的基本輸入。
最終價(jià)格的差價(jià)相當(dāng)大，從大約 45 美元到 500 美元不等！
在當(dāng)前的格式中，由于圖表中充滿了數(shù)據(jù)，因此很難真正清楚地看到正在發(fā)生的事情——所以這就是我們回到之前刪除的直方圖的地方，盡管這次它會(huì)向我們展示 結(jié)束模擬值的分布，而不是單個(gè)模擬的每日收益分布。這次我還模擬了 10,000 次運(yùn)行，以便為我們提供更多數(shù)據(jù)。
同樣，代碼很容易調(diào)整以包含此直方圖。
import numpy as npimport mathimport matplotlib.pyplot as pltfrom scipy.stats import norm#set up empty list to hold our ending values for each simulated price seriesresult = []#Define VariablesS = apple['Adj Close'][-1] #starting stock price (i.e. last available real stock price)T = 252 #Number of trading daysmu = 0.2309 #Returnvol = 0.4259 #Volatility#choose number of runs to simulate - I have chosen 10,000for i in range(10000):#create list of daily returns using random normal distributiondaily_returns=np.random.normal(mu/T,vol/math.sqrt(T),T)+1#set starting price and create price series generated by above random daily returnsprice_list = [S]for x in daily_returns:price_list.append(price_list[-1]*x)#plot data from each individual run which we will plot at the endplt.plot(price_list)#append the ending value of each simulated run to the empty list we created at the beginningresult.append(price_list[-1])#show the plot of multiple price series created aboveplt.show()#create histogram of ending stock values for our mutliple simulationsplt.hist(result,bins=50)plt.show()
輸出如下：
我們現(xiàn)在可以快速計(jì)算分布的平均值以獲得我們的“預(yù)期值”：
#use numpy mean function to calculate the mean of the resultprint(round(np.mean(result),2))
輸出188.41。當(dāng)然，您會(huì)得到略有不同的結(jié)果，因?yàn)檫@些是隨機(jī)每日回報(bào)抽取的模擬。您在每次模擬中包含的路徑或運(yùn)行次數(shù)越多，平均值就越傾向于我們用作“mu”輸入的平均回報(bào)。這是大數(shù)定律的結(jié)果。
我們還可以查看潛在價(jià)格分布的幾個(gè)“分位數(shù)”，以了解非常高或非常低回報(bào)的可能性。
我們可以使用 numpy 的“percentile”函數(shù)來計(jì)算 5% 和 95% 的分位數(shù)：
print("5% quantile =",np.percentile(result,5))print("95% quantile =",np.percentile(result,95))
5% quantile = 85.0268905204829495% quantile = 344.5558966477557
我們現(xiàn)在知道，股價(jià)有 5% 的可能性最終會(huì)低于 85.02 美元，有 5% 的可能性會(huì)高于 344.55 美元。
我們可以開始問自己這樣的問題：“我是否愿意冒 5% 的風(fēng)險(xiǎn)最終獲得價(jià)值低于 63.52 美元的股票，以追逐股價(jià)在 188.41 美元左右的預(yù)期回報(bào)？”
最后要做的就是在直方圖上快速繪制我們剛剛計(jì)算的兩個(gè)分位數(shù)，以給我們一個(gè)直觀的表示。掃描本文最下方二維碼獲取全部完整源碼和Jupyter Notebook 文件打包下載。
plt.hist(result,bins=100)plt.axvline(np.percentile(result,5), color='r', linestyle='dashed', linewidth=2)plt.axvline(np.percentile(result,95), color='r', linestyle='dashed', linewidth=2)plt.show()
↓↓長按掃碼獲取完整源碼↓↓