Python中的蒙特卡洛

Question

已编辑以包含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”，并且在Visual Studio Code中以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

经过编辑以包含VBA代码以进行比较另外，我们知道蒙特卡洛应收敛的分析值8.021，这使比较更加容易。 Excel VBA给出8.067 ...

Answer 1

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

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