Monte Carlo 在 Python

Question

编辑以包含 VBA 比较代码

此外，我们知道蒙特卡罗应该收敛的分析值，即 8.021，这使得比较更容易。

Excel VBA 给出 8.067 基于平均 5 次蒙特卡洛模拟 (7.989, 8.187, 8.045, 8.034, 8.075)

Python 基于 5 个 MC（7.913、7.915、8.203、7.739、8.095）和更大的方差得出 7.973！

VBA代码甚至没有"that good"，使用一种相当糟糕的方式从标准正态生成样本！

我正在运行在 Python 中编写一个超级简单的代码，通过 Monte Carlo 为欧式看涨期权定价，我很惊讶 "bad" 与10,000 "simulated paths"。通常，当运行在 C++ 甚至 VBA 中为这个简单的问题使用 Monte-Carlo 时，我会得到更好的收敛。

我显示下面的代码（代码取自教科书"Python for Finance"，我运行在Visual Studio代码下Python 3.7.7，64位版本）：我得到以下结果，例如：运行 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)+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))) + (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

Answer 1

我不认为 Python 或 numpy 内部结构有任何问题，无论您使用什么工具，收敛性绝对应该是相同的。我运行一些具有不同样本大小和不同 sigma 值的模拟。毫不奇怪，事实证明收敛速度在很大程度上受 sigma 值控制，请参见下图。请注意，x 轴是对数刻度！在较大的振荡消失后，在稳定之前会有更多较小的波浪。在 sigma=0.5 时最容易看到。

我绝对不是专家，但我认为最明显的解决方案是增加样本量，正如您提到的那样。很高兴看到来自 C++ 或 VBA 的结果和代码，因为我不知道您对 numpy 和 python 函数有多熟悉。也许有些东西并没有按照你认为的那样去做。

生成剧情的代码（先不谈效率，太烂了）：

import numpy as np
import matplotlib.pyplot as plt

S_0 = 100.  # initial value
K = 105.    # strike
T = 1.0     # time to maturity
r = 0.05    # short rate (constant)

fig = plt.figure()
ax = fig.add_subplot()
plt.xscale('log')

samplesize = np.geomspace(1000, 20000000, 64)
sigmas = np.arange(0, 0.7, 0.1)
for s in sigmas:
    arr = []

    for n in samplesize:

        n = n.astype(int)

        z = np.random.standard_normal(n)

        S_T = S_0 * np.exp((r-0.5*s**2)+np.sqrt(T)*s*z)


        C_T = np.maximum((S_T-K),0) 


        C_0 = np.exp(-r*T)*np.average(C_T) 


        arr.append(C_0)

    ax.scatter(samplesize, arr, label=f'sigma={s:.2f}')

plt.tight_layout()
plt.xlabel('Sample size')
plt.ylabel('Value')
plt.grid()
handles, labels = ax.get_legend_handles_labels()
plt.legend(handles[::-1], labels[::-1], loc='upper left')
plt.show()

加法:

这次您使用 VBA 获得了更接近实际值的结果。但有时你不会。运行domness 的影响在这里太大了。事实上，从低样本数模拟中平均得出 5 个结果是没有意义的。例如，在 Python 中对 50 个不同的模拟进行平均（只有 n=10000，即使你不应该这样做，如果你愿意得到正确的答案）产生 8.025167（± 0.039717，置信水平为 95% ), 非常接近真实解

Monte Carlo 在 Python

Monte Carlo in Python

python

finance

montecarlo