Conditional Monte Carlo for Option Pricing with Random Volatility

# Option price estimates are inconsistent between two separate formulations


# Option price estimates are inconsistent between two separate formulations












I am currently studying slides from Prof. Lech A. Grzelak's lectures on Financial Engineering.

In an analysis on the application of conditional expectations, the student is expected to calculate the value of an option at $ t = t_0$ using a Monte-Carlo (MC) approach on the stock price process v/s Monte-Carlo on a Conditional Expectation-based formulation. The slides suggest that the latter is a more stable estimate while MC on the price process can be more unstable (and computationally expensive). This suggestion is what I have understood about the exercise.

To do this, the professor has given the following SDE: $$ dS(t) = rS(t)dt + JS(t)dW^{\mathbb{Q}}(t)$$ where $J \sim \mathcal{N}(\mu_J, \sigma_j^2)$ acts as a stochastic volatility random variable.

From the reference code provided I was able to understand the following numerical experiment:

#### Numerical Experiment - 1: On the Stability of Estimation

- Write two separate estimation strategies for option pricing: SDE Simulation and Direct Application of the Pay-Off Formula: Write a function to generate sample paths of the price process. Write a function to evaluate pay-off of each path at maturity (i.e., at $t = T$). For a Call Option, this would be $\max(S(T) - K, 0)$. Calculate the average (or expected value of each pay-off) and discount it by $T$. This would then be: $e^{-rT}\mathbb{E}[\max(S(T)-K, 0)]$. This is the price of the option at $t = t_0$ by a pure MC simulation approach. Conditional Expectation over realizations of $J$ Write a function to calculate the Black-Scholes option price (say $f(S_0, K, \sigma, t, T, r)$. Write a function to average the price of the option for each realization of $J$ (i.e., sample from the distribution that defines stochastic volatility). This would look like $\mathbb{E}[f(S_0, K, j, t, T, r)|J = j]$ (I simply use $\sigma = j$). This is the price of the option based on conditional expectation (that requires averaging over the distribution of $J$, instead of simulating the entire stochastic process and discounting the pay off of the function).



Upon running my experiment, I observed the following for the strike price $K = 80$: While the professor's results look like the following:

The code I have written (which is near identical to the professor's code) is as follows:

```
import enum
import scipy.stats as st
import numpy 
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

class OptionType(enum.Enum):
    CALL = 1.00
    PUT = -1.00

def generate_paths(no_of_paths, no_of_steps, S0, T, muJ, sigmaJ, r):

    X = np.zeros([no_of_paths, no_of_steps+1])
    S = np.zeros([no_of_paths, no_of_steps+1])
    time = np.zeros([no_of_steps+1])

    dt = T / float(no_of_steps)
    X[:, 0] = np.log(S0)
    S[:, 0] = S0

    # Simulate the standard normal
    Z = np.random.normal(0, 1, [no_of_paths, no_of_steps])
    # Simulate stochastic volatility
    J = np.random.normal(muJ, sigmaJ, [no_of_paths, no_of_steps])

    for i in range(0, no_of_steps):
        if no_of_paths > 1:
            Z[:, i] = (Z[:, i] - np.mean(Z[:, i])) / np.std(Z[:, i])
        
        X[:, i+1] = X[:, i] + (r - 0.5*J[:, i]**2)*dt + J[:, i]*Z[:, i]*np.sqrt(dt)
        time[i+1] = time[i] + dt
    
    S = np.exp(X)
    paths = {"time": time, "S": S, "X": X, "J": J}
    return paths

def european_option_price_from_MC_paths(CP, S, K, T,r):
    # S is a vector of MC samples at T
    if CP == OptionType.CALL:
        return np.exp(-r*T)*np.mean(np.maximum(S - K, 0))
    if CP == OptionType.PUT:
        return np.exp(-r*T)*np.mean(np.maximum(K - S, 0))

def BS_call_put_option_price(CP, S_0, K, sigma, t, T, r):
    # Ensure K is a numpy array with shape (n, 1) if needed, but allow scalar as well
    K_arr = np.array(K)
    if K_arr.ndim == 0:
        K_arr = K_arr.reshape(1, 1)
    elif K_arr.ndim == 1:
        K_arr = K_arr.reshape(-1, 1)
    else:
        # If K already has correct dimensions, do nothing
        pass
    
    d1 = (np.log(S_0 / K_arr) + (r + 0.5*sigma**2)*(T - t)) / (sigma*np.sqrt(T - t))
    d2 = d1 - sigma*np.sqrt(T - t)

    if CP == OptionType.CALL:
        return st.norm.cdf(d1)*S_0 - st.norm.cdf(d2)*K_arr*np.exp(-r*(T-t))
    if CP == OptionType.PUT:
        return st.norm.cdf(-d2)*K_arr*np.exp(-r*(T-t)) - st.norm.cdf(-d1)*S_0

def call_option_conditional_expectation(no_of_paths, T, S0, K, J, r):

    # sigma at time T:
    J_i = J[:, -1]

    result = np.zeros([no_of_paths])

    for j in range(0, no_of_paths):
        sigma_j = J_i[j]
        result[j] = BS_call_put_option_price(OptionType.CALL, S0, K, sigma_j, 0, T, r)
    
    return np.mean(result)
K = np.array([80])
CP = OptionType.CALL

N_grid = range(0,20000, 100)
N_runs = len(N_grid)

result_MC = np.zeros([N_runs])
result_Cond_Exp = np.zeros([N_runs])

for (i, N) in enumerate(N_grid):
    paths = generate_paths(N, no_of_steps, S_0, T, muJ, sigmaJ, r)
    time_grid = paths["time"]
    S = paths["S"]

    result_MC[i] = european_option_price_from_MC_paths(CP, S, K, T, r)
    
    J = paths["J"]

    result_Cond_Exp[i] = call_option_conditional_expectation(N, T, S_0, K, J, r)
```

#### Questions:

- Why is there about $\sim 6.17$ units of error in MC estimates v/s Conditional Expectation estimates?

- Assuming that the professor's estimates are correct (in the attached picture from his slides), then MC estimates are missing something. What could it be?