Монте-Карло в Python

Question

Отредактировано для включения кода VBA для сравнения

Кроме того, мы знаем аналитическое значение, равное 8,021, к которому должен сходиться Монте-Карло, что делает сравнение проще.

Excel VBA выдает 8,067 на основе усреднения 5 симуляций Монте-Карло (7,989, 8,187, 8,045, 8,034, 8,075)

Python дает 7,973 на основе 5 MC (7,913, 7,915, 8,203, 7,739, 8,095) и большей дисперсии!

Код VBA даже не "такой хороший", используя довольно плохой способ получения сэмплов из Standard Normal!

Я использую супер простой код в Python, чтобы оценить European Call Option через Монте-Карло, и я удивлен тем, насколько "плохой" является конвергенция с 10 000 «смоделированных путей». Обычно, когда я запускаю Монте-Карло для этой простой проблемы в C ++ или даже VBA, у меня улучшается сходимость.

Я показываю код ниже (код взят из Учебника "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

Answer 1

Я не думаю, что что-то не так с внутренностями Python или numpy, конвергенция определенно должна быть одинаковой независимо от того, какой инструмент вы используете. Я провел несколько симуляций с разными размерами выборки и различными значениями сигма. Неудивительно, что скорость сходимости сильно зависит от значения сигмы, см. График ниже. Обратите внимание, что ось х находится в логарифмическом масштабе! После того, как большие колебания исчезают, появляются более мелкие волны, прежде чем они стабилизируются. Наиболее легко увидеть при сигме = 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 результатов моделирования с низким числом образцов не имеет смысла. Например, усреднение 50 различных симуляций в Python (только с n = 10000, даже если вы не должны этого делать, если хотите получить правильный ответ), приводит к 8,025167 (± 0,039717 с уровнем достоверности 95%), что очень близко к реальному решению.