时间戳

首页 > 解决方案 > Python 中的蒙特卡洛

python - Python 中的蒙特卡洛

问题描述

编辑以包含 VBA 代码以进行比较

此外，我们知道 Monte-Carlo 应该收敛的解析值 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 中运行一个超级简单的代码，通过蒙特卡洛为欧洲看涨期权定价，我对 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)+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

标签： pythonfinancemontecarlo

解决方案

我不认为 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 获得了更接近实际值的结果。但有时你不会。这里随机性的影响太大了。事实是，从低样本数模拟中仅平均 5 个结果是没有意义的。例如，在 Python 中平均 50 次不同的模拟（只有 n=10000，即使你不应该这样做，如果你愿意得到正确的答案）会产生 8.025167（± 0.039717，置信度为 95%），即非常接近真正的解决方案。

推荐阅读