Python 蒙地卡羅收斂
編輯以包含 VBA 程式碼以進行比較
此外,我們知道簡單看漲期權的分析值,即 8.021,蒙地卡羅應該收斂到該值,這使得比較更容易。
Excel VBA 根據平均 5 次蒙特卡羅模擬(7.989、8.187、8.045、8.034、8.075)給出 8.067
Python 基於 5 個 MC(7.913、7.915、8.203、7.739、8.095)給出 7.973 和更大的變異數!
VBA 程式碼使用一種相當糟糕的方式從標準正常生成樣本……
我正在 Python 中執行一個超級簡單的程式碼,通過 Monte Carlo 為歐洲看漲期權定價,我對 10,000 條“模擬路徑”的收斂有多“糟糕”感到驚訝。通常,當在 C++ 甚至 VBA 中針對這個簡單問題執行 Monte-Carlo 時,我會獲得更好的收斂性。
我展示了下面的程式碼(程式碼取自教科書“Python for Finance”,我在 Python 3.7.7,64 位版本下的 Visual Studio Code 中執行):我得到以下結果,例如:Run 1 = 7.913,執行 2 = 7.915,執行 3 = 8.203,執行 4 = 7.739,執行 5 = 8.095,
上述結果相差如此之大,將是不可接受的。如何提高收斂性???(顯然通過執行更多路徑,但正如我所說:對於 10,000 條路徑,結果應該已經收斂得更好):
#MonteCarlo valuation of European Call Option import math import numpy as np #Parameter Values S_0 = 100. # initial value K = 105. # strike T = 1.0 # time to maturity r = 0.05 # short rate (constant) sigma = 0.2 # vol nr_simulations = 10000 #Valuation Algo: # Notice the vectorization below, instead of a loop z = np.random.standard_normal(nr_simulations) # Notice that the S_T below is a VECTOR! S_T = S_0 * np.exp((r-0.5*sigma**2)*T+math.sqrt(T)*sigma*z) #Call option pay-off at maturity (Vector!) C_T = np.maximum((S_T-K),0) # C_0 is a scalar C_0 = math.exp(-r*T)*np.average(C_T) print('Value of the European Call is: ', C_0)
我還包括了 VBA 程式碼,它產生了更好的結果(在我看來):使用下面的 VBA 程式碼,我得到 7.989、8.187、8.045、8.034、8.075。
Option Explicit Sub monteCarlo() ' variable declaration ' stock initial & final values, option pay-off at maturity Dim stockInitial, stockFinal, optionFinal As Double ' r = rate, sigma = volatility, strike = strike price Dim r, sigma, strike As Double 'maturity of the option Dim maturity As Double ' instatiate variables stockInitial = 100# r = 0.05 maturity = 1# sigma = 0.2 strike = 105# ' normal is Standard Normal Dim normal As Double ' randomNr is randomly generated nr via "rnd()" function, between 0 & 1 Dim randomNr As Double ' variable for storing the final result value Dim result As Double Dim i, j As Long, monteCarlo As Long monteCarlo = 10000 For j = 1 To 5 result = 0# For i = 1 To monteCarlo ' get random nr between 0 and 1 randomNr = Rnd() 'max(Rnd(), 0.000000001) ' standard Normal normal = Application.WorksheetFunction.Norm_S_Inv(randomNr) stockFinal = stockInitial * Exp((r - (0.5 * (sigma ^ 2)))*maturity + (sigma * Sqr(maturity) * normal)) optionFinal = max((stockFinal - strike), 0) result = result + optionFinal Next i result = result / monteCarlo result = result * Exp(-r * maturity) Worksheets("sheet1").Cells(j, 1) = result Next j MsgBox "Done" End Sub Function max(ByVal number1 As Double, ByVal number2 As Double) If number1 > number2 Then max = number1 Else max = number2 End If End Function
蒙地卡羅的收斂與您使用的程式語言幾乎沒有關係,並且可以通過分佈和中心極限定理來解釋。您應該能夠用不同的語言實現完全相同的東西。
感受融合的最佳方式是將其形象化。
import matplotlib.pyplot as plt from scipy.stats import norm N = norm.cdf def call(S, K, T, r, vol): d1 = (np.log(S/K) + (r + 0.5*vol**2)*T) / (vol*np.sqrt(T)) d2 = d1 - vol * np.sqrt(T) return S * N(d1) - np.exp(-r * T) * K * N(d2) df = np.exp(-r*T) prices = [df*C_T[:n].mean() for n in range(1,nr_simulations)] plt.plot(prices) plt.axhline(call(S_0, K, T, r, sigma), c='r')