Monte Carlo Option Pricing and Risk-Neutral Sampling Errors

# Monte carlo Black scholes option pricer only converges for small sample space


# Monte carlo Black scholes option pricer only converges for small sample space












I am a non-professional interested in quantitative finance, with some knowledge of stochastic calculus and the Black–Scholes model.

I have implemented a Black–Scholes option pricer that computes the expectation of the payoff under the risk-neutral measure, which is constructed using Girsanov’s theorem and approximated computationally using a Monte Carlo method.

The following code prices a European call option. There is a unit test that compares the Monte Carlo price of the European call option with the analytical Black–Scholes formula. For small sample spaces, the Monte Carlo price is close to the analytical solution. However, as the sample space size increases, the Monte Carlo price diverges from the analytical solution. Any advice would be appreciated, including any small comments on code quality or style.

Full code can be found here:https://github.com/p0gster/Black-Scholes-Option-Pricer

The following unit test succeeds for small sample spaces (approximately 10) but fails as the sample space size increases.

```

class TestOptionPricing(unittest.TestCase):
    def test_monte_carlo_vs_black_scholes(self):
        # Parameters
        S0 = 100.0
        K = 100.0
        mu = 0.05
        sigma = 0.2
        r = 0.05
        T = terminal_time  # Use the same terminal_time as in your Monte Carlo code

        # Monte Carlo price using your RandomVariable approach
        payoff_rv = european_call_payoff_rv(S0, mu, sigma, K)
        monte_carlo_price = price_option(payoff_rv, mu, sigma, r)

        # Analytical Black-Scholes price
        d1 = (math.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * math.sqrt(T))
        d2 = d1 - sigma * math.sqrt(T)
        analytical_price = S0 * norm.cdf(d1) - K * math.exp(-r * T) * norm.cdf(d2)

        # Allow some tolerance due to small sample space
        tolerance = 5.0  # increase tolerance if sample_space_size is small
        self.assertAlmostEqual(monte_carlo_price, analytical_price, delta=tolerance)
```

The main body is here:

```

from typing import Callable, List
import unittest
import math
from scipy.stats import norm
import numpy as np

import requests
import pandas as pd

import matplotlib.pyplot as plt

# Sample space is the set of integers from 0 to sample_space_size - 1

# Define the size of the sample space
sample_space_size = 10  # You can change this number
terminal_time=10

def interval(start: int, end: int) -> list:
    """Return a list representing the interval from start to end inclusive."""
    if end < start:
        raise ValueError("End of interval must be greater than or equal to start")
    return list(range(start, end + 1))

sample_space=interval(0,sample_space_size-1)

# Random Variable Abstract Data Type
class RandomVariable:
    def __init__(self, func: Callable[[int], float], sample_space: List[int]):
        """Construct a random variable from a function and a sample space (as a list)."""
        if not callable(func):
            raise TypeError("Function must be callable")
        if not isinstance(sample_space, list):
            raise TypeError("Sample space must be provided as a list")
        self.func = func
        self.sample_space = list(set(sample_space))  # internally store as a unique list

    def evaluate(self, outcome: int) -> float:
        return self.func(outcome)

    def values(self) -> List[float]:
        return [self.func(o) for o in self.sample_space]

# Procedure: normally distributed random variable

def normally_distributed_random_variable(mean: float, variance: float) -> RandomVariable:
    """Return a RandomVariable representing a normally distributed real value.
    The RNG is seeded so the value is reproducible for each outcome.
    Each outcome gets its own sampled value.
    """
    def seed_to_normal_rv_value(seed):
      rng = np.random.default_rng(seed)
      std = variance ** 0.5
      # Precompute values for the discrete sample space
      return rng.normal(mean, std)

    return RandomVariable(lambda sample: seed_to_normal_rv_value(sample),sample_space)

# Probability Measure Abstract Data Type
class ProbabilityMeasure:
    def __init__(self, func: Callable[[int], float]):
        """Construct a probability measure from the sample space to real numbers (probabilities). 
        Note: sum of probabilities must be one, this code does not check this."""
        if not callable(func):
            raise TypeError("Function must be callable")
        self.func = func

    def probability(self, outcome: int) -> float:
        p = self.func(outcome)
        if not (0 <= p <= 1):
            raise ValueError(f"Probability must be between 0 and 1, got {p}")
        return p

    def probabilities(self, space: List[int]) -> List[float]:
        return [self.probability(o) for o in space]

def standard_probability_measure() -> ProbabilityMeasure:
    """
    Return the standard (uniform) probability measure
    on the finite sample space.
    """
    return ProbabilityMeasure(lambda _: 1 / sample_space_size)

def create_risk_neutral_measure(mu: float, sigma: float, r: float) -> ProbabilityMeasure:
    """
    Create a risk-neutral probability measure using Girsanov's theorem.

    Inputs:
    - mu: drift of the underlying
    - sigma: volatility
    - r: risk-free interest rate

    Returns:
    - ProbabilityMeasure instance representing the risk-neutral measure
    """
    global terminal_time, sample_space_size

    # Create Wiener random variable with mean 0 and variance terminal_time
    wiener_rv = normally_distributed_random_variable(0.0, terminal_time)

    # Market price of risk
    theta = (mu - r) / sigma

    # Define the risk-neutral probability function
    def p_star(omega: int) -> float:
        return (1 / sample_space_size) * math.exp(
            -theta * wiener_rv.evaluate(omega)
            - 0.5 * theta**2 * terminal_time
        )

    return ProbabilityMeasure(p_star)

# Expectation function
def expectation(rv: RandomVariable, pm: ProbabilityMeasure) -> float:
    """Return the expectation of the random variable with respect to the probability measure.
    Computed as the sum over the sample space of rv(outcome) * pm(outcome).
    """
    return sum(rv.evaluate(o) * pm.probability(o) for o in rv.sample_space)

def price_option(payoff: RandomVariable,
                 mu: float,
                 sigma: float,
                 r: float) -> float:
    """
    Price an option at time 0 using risk-neutral valuation.

    The option price is the expectation, under the risk-neutral
    probability measure, of the discounted payoff.

    Inputs:
    - payoff: RandomVariable representing the payoff at terminal time T
    - mu: drift of the underlying
    - sigma: volatility
    - r: risk-free interest rate

    Returns:
    - Option price at time 0
    """
    # Construct the risk-neutral probability measure
    pm_star = create_risk_neutral_measure(mu, sigma, r)

    # Discounted expectation under the risk-neutral measure
    discounted_expectation = math.exp(-r * terminal_time) * expectation(payoff, pm_star)

    return discounted_expectation

def european_call_payoff_rv(S0: float, mu: float, sigma: float, K: float) -> RandomVariable:
    """
    Create a RandomVariable representing the payoff of a European call option at terminal time.

    Payoff: (S_T - K)+
    where S_T = S0 * exp((mu - 0.5*sigma^2)*T + sigma*W_T)

    Inputs:
    - S0: initial stock price
    - mu: drift of the underlying
    - sigma: volatility
    - K: strike price

    Returns:
    - RandomVariable representing the call option payoff at T
    """
    # Wiener process at terminal time T
    W_T = normally_distributed_random_variable(0.0, terminal_time)

    # Define the payoff function
    def payoff(omega: int) -> float:
        S_T = S0 * math.exp((mu - 0.5 * sigma**2) * terminal_time + sigma * W_T.evaluate(omega))
        return max(S_T - K, 0.0)

    return RandomVariable(payoff,sample_space)
```