Python Monte-Carlo Convergence

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
• It won’t matter here since maturity=1 but as a side comment you should multiply r-0.5sigma^2 by T or maturity inside the exp.
– Ivan
Apr 23 '20 at 13:38
• Thanks @Ivan, I omitted it because it's 1, but I include it for completeness now, thanks for pointing out. Apr 23 '20 at 14:42

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') • @JanStuller, the convergence has nothing to do with hardware or software. Apr 23 '20 at 16:01

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.