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代码甚至还不是“那么好”，而是使用了一种相当糟糕的方式来从Standard Normal生成样本！

我正在用Python运行一个超级简单的代码，以通过Monte Carlo为欧式看涨期权定价，而对于10,000个“模拟路径”的融合有多么“糟糕”，我感到惊讶。通常，在C ++甚至VBA中针对此简单问题运行Monte-Carlo时，我会获得更好的收敛性。

我在下面显示代码（该代码摘自教科书“ Python for Finance”，并且在Visual Studio Code中以Python 3.7.7（64位版本）运行）：作为示例，我得到以下结果：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