Python - Problem of random numbers in MC simulation

Question

I am interested in estimating the price of a European Call Option using the Montecarlo simulation, to get a good approximation of the analytical Black Scholes formula, so a very simple task. Performing the simulation, I have noticed one thing: if I change the seed for the random numbers, I get different results. Is there any procedure that allows me to get the same price, regardless of the chosen seed? (maybe not exactly the same price, but some strategy to reduce the effect of chosen seed).

As we can see from this pictures, here we have the difference in absolute value between the analytical price and estimated one, with different seed. Here we can see that the best seed is the number 100.

X-axis = Value of spot prices in interval [1,1.2]

Y-axis = difference in absolute value between estimated option and analytical one

number of simulation = 200000

K = 1

r = q = 0

sigma = 0.5

Here I have attached the code used to generate the MC simulation for different kind of seed.

from math import log, sqrt, pi, exp
from scipy.stats import norm
from datetime import datetime, date
import numpy as np
import pandas as pd
from pandas import DataFrame

def Option_MC(s0, K, T, j, sigma, r, q, seed):
    Rand = np.random.RandomState()
    Rand.seed(seed)
    S = np.ndarray(shape=(2, j), dtype=np.double)
    S[0] = s0
    Brownian_Motion = Rand.normal(-.5 * sigma * sigma * T, sigma * sqrt(T), (1, j))
    S[1] = S[0] * np.exp(Brownian_Motion)
    P = np.maximum(S[1] - K, 0)
    mean = np.average(P)
    return mean

def d1(S,K,T,r,sigma):
    return (log(S/K)+(r+sigma**2/2.)*T)/(sigma*sqrt(T))
def d2(S,K,T,r,sigma):
    return d1(S,K,T,r,sigma)-sigma*sqrt(T)
def bs_call(S,K,T,r,sigma):
    return S*norm.cdf(d1(S,K,T,r,sigma))-K*exp(-r*T)*norm.cdf(d2(S,K,T,r,sigma))

analytical = []
for i in arrays:
    analytical.append(bs_call(i,1.,1,0,0.5))

arrays = np.linspace(1,1.2,100)
seeds = [1,10,100,1000]
price_tot = []
for seed in seeds:
    price = []
    for i in arrays:
        price.append(Option_MC(i,1.,1,200000,0.5,0,0,seed))
    price_tot.append(price)

for p,seed in zip(price_tot,seeds):
    plt.plot(arrays,abs(np.array(p)-np.array(analytical)),label=f'Seed = {seed}')
plt.legend()
plt.show()

for p,seed in zip(price_tot,seeds):
    plt.plot(arrays,p,label=f'Seed = {seed}')
plt.legend()
plt.show()

Answer 1

The quickest and easiest (for you, not the computer) is to increase the number of simulations. When I use 1 million sims the error is smaller:

You can find more sophisticated techniques to reduce variance of your estimate in Glasserman's "Monte Carlo methods in financial engineering". For option pricing you can look at Control Variates and Antithetic Variates, they are relatively easy to implement.

Answer 2

I believe that Python's Rand uses some kind of Linear Congruential Generator to recover "random" numbers.

Think about the random generators as algorithms (deterministic ones) that generate a sequence of numbers in [0, 1]. For example: $$0.0234, \ 0.382, \ 0.375432, \ 0.392, \ 0.785, \ ...$$ Those numbers are then "transformed" (see for example inversion of CDF or AR) to Normal "random" drawings.

Think about the seed as the "starting point" of this [0, 1] sequence. Since the algorithm generating the sequence is deterministic, giving the same seed will result always in the same sequence. Changing the seed, you change the starting point and so the whole sequence. If you change the whole sequence your result will always be different.

What you want to achieve is the Monte Carlo convergence. The greater the number of simulations (not always) the lower the standard error of the samples random variables.

Tip: try to avoid using low seeds. At least 1000+ should be used.