Python Monte-Carlo Convergence

Question

Edited to include VBA code for comparison

Also, we know the analytical value of the simple Call option, which is 8.021, towards which the Monte-Carlo should converge, which makes the comparison easier.

Excel VBA gives 8.067 based on averaging 5 Monte-Carlo simulations (7.989, 8.187, 8.045, 8.034, 8.075)

Python gives 7.973 based on 5 MCs (7.913, 7.915, 8.203, 7.739, 8.095) and a larger Variance!

The VBA code is using a rather bad way to produce samples from Standard Normal...

I am running a super simple code in Python to price European Call Option via Monte Carlo, and I am surprised at how "bad" the convergence is with 10,000 "simulated paths". Usually, when running a Monte-Carlo for this simple problem in C++ or even VBA, I get better convergence.

I show the code below (the code is taken from Textbook "Python for Finance" and I run in in Visual Studio Code under Python 3.7.7, 64-bit version): I get the following results, as an example: Run 1 = 7.913, Run 2 = 7.915, Run 3 = 8.203, Run 4 = 7.739, Run 5 = 8.095,

Results such as the above, that differ by so much, would be unacceptable. How can the convergence be improved??? (Obviously by running more paths, but as I said: for 10,000 paths, the result should already have converged much better):

#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)*T+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)

I also include VBA code, which produces slightly better results (in my opinion): with the VBA code below, I get 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)))*maturity + (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

The convergence of your monte carlo has little to do with the programming language you are using and is explained by the distribution and the central limit theorem. You should be able to implement the exact same thing in different languages.

The best way to get a feel for the convergence would be to visualize it.

import matplotlib.pyplot as plt
from scipy.stats import norm
N = norm.cdf

def call(S, K, T, r, vol):
    d1 = (np.log(S/K) + (r + 0.5*vol**2)*T) / (vol*np.sqrt(T))
    d2 = d1 - vol * np.sqrt(T)
    return S * N(d1) - np.exp(-r * T) * K * N(d2)

df = np.exp(-r*T)
prices = [df*C_T[:n].mean() for n in range(1,nr_simulations)]

plt.plot(prices)
plt.axhline(call(S_0, K, T, r, sigma), c='r')

Answer 2

In order to test the convergence you can add the computation of confidence interval at a given probability. We add in the code the computation of the variance and the interval.

$V\_T = \frac{1}{M-1} \times \sum (e^{-rT}C\_T - C\_0)^2$

$I_M = [C_0-q\times \sqrt{\frac{V\_T}{M}}; C_0+q\times \sqrt{\frac{V\_T}{M}}]$

With $M$ the number of simulation and a confidence level of 95% for $q = 2$ and 99% for $q = 3$

V_T = 1/(nr_simulations-1)*sum((np.exp(-r*T)*C_T-C_0)**2)
q=3
I_m = C_0 - q*math.sqrt(V_T/nr_simulations)
I_M = C_0 + q*math.sqrt(V_T/nr_simulations)
print("Price between " + str(I_m) + " and " + str(I_M) + " with probability 99%.")

Result for $q = 3$ $I_M = [7.59; 8.37]$

So the values

Run 1 = 7.913, Run 2 = 7.915, Run 3 = 8.203, Run 4 = 7.739, Run 5 = 8.095,

are coherent with a confidence interval at 99% for 10000 simulations.