ML-Powered Monte Carlo Option Pricer
Train an XGBoost surrogate model to price options 1000x faster than Monte Carlo.
Introduction
30 Second Summary
When financial firms need to price a stock option, they often run thousands of random simulations to estimate its value. This is accurate but painfully slow. What if a machine learning model could learn the same pricing function and deliver answers instantly?
In this project, you will build an ML model that learns to replicate a Monte Carlo simulation for option pricing. Your model will deliver the same accuracy at a fraction of the speed.
What You'll Build
A complete ML surrogate pipeline where a Monte Carlo option pricing engine generates training data, then an XGBoost regressor learns the pricing function from those samples.
By the end of this project, you'll have:
- A vectorized Monte Carlo pricing engine for European options using Geometric Brownian Motion simulation.
- A trained XGBoost surrogate model that prices options ~1000x faster than the full simulation with less than 5% error.
- Publication-ready benchmark visualizations comparing speed and accuracy side by side.
- Secret Mission: Estimate option Greeks (Delta and Vega) using both your MC engine and ML model.
Are there any prerequisites?
You should be comfortable with basic Python (loops, functions, arrays) and have some familiarity with probability distributions and expected values. No prior finance or machine learning experience is required.
Set Up Your Environment
Your packages are installed and your virtual environment is active. Now you need a workspace where you can build the pipeline incrementally, running one piece at a time.
Jupyter notebooks let you run code in small chunks called cells. This makes them ideal for iterative work like building a simulation engine, generating training data, and visualizing results.
In this step, get ready to:
- Create a Jupyter notebook for the project.
- Verify all required packages import successfully.
Create your Jupyter notebook and verify imports
You need a single notebook file that will hold every cell of the ML pipeline, from Monte Carlo simulation all the way to benchmark charts.
- In VS Code, press Cmd+Shift+P (macOS) or Ctrl+Shift+P (Windows) to open the Command Palette.
- Type Create: New Jupyter Notebook and select it from the list.
- VS Code opens a blank notebook. Save it immediately with Cmd+S (macOS) or Ctrl+S (Windows).
- Name the file mc_ml_pricer.ipynb and save it inside your mc-ml-pricer folder.
- Click Select Kernel in the top right of the notebook.
- Choose the Python environment from your .venv virtual environment. It may appear as Python 3.x ('.venv') in the list.
Why select a kernel?
The kernel is the Python process that runs your notebook cells. By pointing it at your .venv environment, you ensure the notebook can access NumPy, XGBoost, and all the other packages you installed.
- Paste the following code into the first cell of the notebook:
# Cell 1: Imports and Setup
import numpy as np
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
import matplotlib.pyplot as plt
import time
from scipy.stats import norm
print(f"XGBoost version: {xgb.__version__}")
print("All imports successful!")What does this code do?
- numpy provides vectorized math for simulating thousands of stock price paths at once.
- xgboost is the gradient-boosted decision tree library that will learn to predict option prices.
- sklearn supplies train/test splitting and error metrics (RMSE, MAPE) for evaluating your model.
- matplotlib creates the benchmark visualizations (scatter plots, histograms, bar charts).
- scipy.stats.norm provides the cumulative normal distribution function needed for Black-Scholes validation later.
- time measures execution speed so you can compare Monte Carlo vs ML inference.
- Run the cell by pressing Shift+Enter.
- You should see two lines of output below the cell: XGBoost version: 3.x.x (where 3.x.x is your installed version) and All imports successful!.
Seeing a ModuleNotFoundError?
- Make sure you selected the correct kernel. Click Select Kernel in the top right and choose the .venv environment. If you pick a different Python installation, it won't have your packages.
- Confirm your virtual environment is activated in the terminal. You should see (.venv) at the start of your prompt. If not, activate it again and re-select the kernel.
- If scipy is missing specifically, it was installed automatically with scikit-learn. Try running pip install scipy in your terminal, then restart the kernel and re-run the cell.
If you're still stuck, help me fix my ModuleNotFoundError in VS Code Jupyter.
✔️ Awesome, I've got everything!
Great. Double check you've saved your notebook file.
ⓧ I'd like to double check the full code
Here is the complete code for Cell 1. Make sure yours matches exactly:
# Cell 1: Imports and Setup
import numpy as np
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
import matplotlib.pyplot as plt
import time
from scipy.stats import norm
print(f"XGBoost version: {xgb.__version__}")
print("All imports successful!")Your notebook is ready and every library is confirmed working. Next up, you'll build the Monte Carlo pricing engine that simulates thousands of stock price paths to value European options.
Build the Monte Carlo Pricing Engine
Your environment is ready and all libraries are imported. Now you need the core simulation engine that prices options by simulating thousands of possible future stock prices.
A Monte Carlo simulation works by generating random scenarios and averaging the results. For option pricing, you simulate where a stock price might end up at expiration, compute the option payoff for each scenario, then average and discount back to today. The more paths you simulate, the more accurate the estimate.
In this step, get ready to:
- Build a vectorized Monte Carlo pricing function using Geometric Brownian Motion.
- Validate Monte Carlo accuracy against Black-Scholes analytical prices.
Implement the vectorized Monte Carlo pricing function
The Geometric Brownian Motion (GBM) model assumes stock prices follow a log-normal distribution. Instead of looping through paths one by one, you can generate all random samples in a single NumPy call and compute all terminal prices at once.
- Create a new cell in your notebook (below Cell 1) and paste the following code:
# Cell 2: Monte Carlo Pricing Engine
def monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type="call", seed=None):
"""
Price a European option using Monte Carlo simulation with GBM.
Parameters:
S: Current stock price
K: Strike price
T: Time to maturity (years)
r: Risk-free rate
sigma: Volatility
n_paths: Number of simulated paths
option_type: 'call' or 'put'
seed: Random seed for reproducibility
"""
if seed is not None:
rng = np.random.default_rng(seed)
else:
rng = np.random.default_rng()
# Simulate terminal stock prices using GBM
Z = rng.standard_normal(n_paths)
S_T = S * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
# Compute payoffs
if option_type == "call":
payoffs = np.maximum(S_T - K, 0)
else:
payoffs = np.maximum(K - S_T, 0)
# Discount to present value
price = np.exp(-r * T) * np.mean(payoffs)
return price
# Test with a known case
test_price = monte_carlo_price(S=100, K=100, T=1, r=0.05, sigma=0.2, n_paths=100000, seed=42)
print(f"MC Call Price (S=100, K=100, T=1, r=0.05, sigma=0.2): {test_price:.4f}")
print("Expected (Black-Scholes): ~10.45")What does this code do?
The function simulates where a stock price could end up at expiration using the GBM formula. It generates all 100,000 random paths in one vectorized call with rng.standard_normal(n_paths) instead of looping. This is what makes it fast.
The GBM formula S * exp((r - 0.5*sigma^2)*T + sigma*sqrt(T)*Z) models stock prices as a drift term (risk-free rate minus a volatility correction) plus a random shock scaled by volatility.
After computing terminal stock prices, the function calculates payoffs (how much the option would be worth) and discounts them back to today using exp(-r*T). The average discounted payoff is the option price estimate.
- Run the cell by pressing Shift+Enter.
- You should see a call price near 10.45 printed in the output. The exact value will vary slightly due to the random seed.
Seeing a NameError?
Make sure Cell 1 (your imports cell) has already been executed. The function depends on numpy being imported as np. Return to Cell 1, click inside it, and press Shift+Enter to re-run it.
help me fix a NameError in my Monte Carlo pricing notebook
Validate against Black-Scholes analytical prices
Your MC engine looks like it works, but how do you know the prices are correct? The Black-Scholes model provides an exact closed-form solution for European option prices. By comparing your MC output against Black-Scholes across multiple parameter sets, you can confirm your simulation is accurate.
- Create a new cell below Cell 2 and paste the following code:
# Cell 3: Black-Scholes Validation
def black_scholes_price(S, K, T, r, sigma, option_type="call"):
"""
Compute the Black-Scholes analytical price for a European option.
"""
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
else:
price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
return price
# Validate MC against BS
print("\nValidation: MC vs Black-Scholes")
print("-" * 50)
test_cases = [
(100, 100, 1.0, 0.05, 0.2, "call"),
(100, 110, 0.5, 0.03, 0.25, "call"),
(100, 90, 1.5, 0.08, 0.3, "put"),
(120, 100, 0.25, 0.04, 0.15, "call"),
]
for S, K, T, r, sigma, opt_type in test_cases:
bs = black_scholes_price(S, K, T, r, sigma, opt_type)
mc = monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type=opt_type, seed=42)
error = abs(mc - bs) / bs * 100
print(f" {opt_type.upper()} S={S} K={K} T={T}: BS={bs:.4f}, MC={mc:.4f}, Error={error:.2f}%")What does this code do?
The black_scholes_price function computes the exact theoretical price using the Black-Scholes formula. It uses norm.cdf (the cumulative normal distribution function you imported from scipy.stats) to calculate the probability-weighted components d1 and d2.
The validation loop tests 4 different option configurations (varying strike, time to expiry, volatility, and option type). For each case, it compares the MC price against the BS analytical price and prints the percentage error. All errors should be well under 1%.
- Run the cell by pressing Shift+Enter.
- You should see a table of 4 test cases with each error below 1%. This confirms your Monte Carlo engine is producing accurate prices.
Seeing errors above 1%?
Check that your monte_carlo_price function in Cell 2 uses n_paths=100000 as the default. Fewer paths increase variance and can push errors above 1%. Re-run Cell 2 to make sure the latest version of the function is loaded.
help me debug my Monte Carlo vs Black-Scholes validation errors
Why validate against Black-Scholes?
Black-Scholes gives you the mathematically exact answer for European options under certain assumptions. Your MC engine approximates the same answer using random simulation. If both agree (within statistical noise), you know the MC engine is correct.
This matters because in the next steps, you will train an ML model on MC output. Garbage in means garbage out. Validating now ensures your training labels are trustworthy.
✔️ Awesome, I've got everything!
Great. Double check you have saved your notebook.
ⓧ I'd like to double check the full code
Here is the complete notebook code through Cell 3:
# Cell 1: Imports and Setup
import numpy as np
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
import matplotlib.pyplot as plt
import time
from scipy.stats import norm
print(f"XGBoost version: {xgb.__version__}")
print("All imports successful!")
# Cell 2: Monte Carlo Pricing Engine
def monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type="call", seed=None):
if seed is not None:
rng = np.random.default_rng(seed)
else:
rng = np.random.default_rng()
Z = rng.standard_normal(n_paths)
S_T = S * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
if option_type == "call":
payoffs = np.maximum(S_T - K, 0)
else:
payoffs = np.maximum(K - S_T, 0)
price = np.exp(-r * T) * np.mean(payoffs)
return price
test_price = monte_carlo_price(S=100, K=100, T=1, r=0.05, sigma=0.2, n_paths=100000, seed=42)
print(f"MC Call Price (S=100, K=100, T=1, r=0.05, sigma=0.2): {test_price:.4f}")
print("Expected (Black-Scholes): ~10.45")
# Cell 3: Black-Scholes Validation
def black_scholes_price(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
else:
price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
return price
print("\nValidation: MC vs Black-Scholes")
print("-" * 50)
test_cases = [
(100, 100, 1.0, 0.05, 0.2, "call"),
(100, 110, 0.5, 0.03, 0.25, "call"),
(100, 90, 1.5, 0.08, 0.3, "put"),
(120, 100, 0.25, 0.04, 0.15, "call"),
]
for S, K, T, r, sigma, opt_type in test_cases:
bs = black_scholes_price(S, K, T, r, sigma, opt_type)
mc = monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type=opt_type, seed=42)
error = abs(mc - bs) / bs * 100
print(f" {opt_type.upper()} S={S} K={K} T={T}: BS={bs:.4f}, MC={mc:.4f}, Error={error:.2f}%")Your Monte Carlo engine is validated and producing accurate prices. Next up, you will use this engine to generate thousands of training samples that your ML model will learn from.
Generate Training Data
Your Monte Carlo engine is accurate, but pricing a single option requires simulating 50,000 to 100,000 random paths.
To train an ML model that can replicate this pricing function instantly, you need thousands of example input-output pairs. In this step, you will sample random option parameters, compute MC prices for each, and organize the results into a training dataset.
In this step, get ready to:
- Create a training data generator that samples random option parameters and computes MC prices.
- Engineer a moneyness feature to help the model learn faster.
- Split the dataset into training and test sets for model evaluation.
Build the training data generator
An ML model learns by seeing many examples of inputs mapped to correct outputs. Here, the inputs are option parameters (spot price, strike, time, rate, volatility) and the outputs are the MC prices your engine computes.
- Create a new cell below Cell 3 in your notebook and paste the following code:
# Cell 4: Generate Training Data
def generate_training_data(n_samples=5000, n_paths=50000, seed=42):
"""
Generate training data by sampling random option parameters
and computing MC prices for each.
"""
rng = np.random.default_rng(seed)
# Sample random parameters
S_arr = rng.uniform(50, 150, n_samples)
K_arr = rng.uniform(50, 150, n_samples)
T_arr = rng.uniform(0.1, 2.0, n_samples)
r_arr = rng.uniform(0.01, 0.10, n_samples)
sigma_arr = rng.uniform(0.1, 0.5, n_samples)
# Compute MC prices
prices = np.zeros(n_samples)
for i in range(n_samples):
prices[i] = monte_carlo_price(
S=S_arr[i], K=K_arr[i], T=T_arr[i],
r=r_arr[i], sigma=sigma_arr[i],
n_paths=n_paths, option_type="call", seed=seed + i
)
# Build feature matrix with engineered moneyness feature
moneyness = S_arr / K_arr
X = np.column_stack([S_arr, K_arr, T_arr, r_arr, sigma_arr, moneyness])
y = prices
return X, y
print("Generating training data (5000 samples)... This may take 2-3 minutes.")
start = time.time()
X, y = generate_training_data(n_samples=5000, n_paths=50000, seed=42)
elapsed = time.time() - start
print(f"Done! Generated {len(y)} samples in {elapsed:.1f} seconds.")
print(f"Feature matrix shape: {X.shape}")
print(f"Price range: [{y.min():.2f}, {y.max():.2f}]")What does this code do?
The generate_training_data function creates a synthetic dataset by sampling random option parameters from uniform distributions.
- Each parameter (S, K, T, r, sigma) is drawn randomly within realistic ranges.
- For every parameter combination, the function calls monte_carlo_price to compute the "true" option price using 50,000 simulated paths.
- The moneyness ratio (S/K) is appended as a sixth feature to help the model learn the spot-strike relationship.
- The function returns a feature matrix X and a target array y of MC prices.
- Run the cell.
- This will take 2-3 minutes as it computes 5,000 Monte Carlo simulations. You should eventually see output confirming 5000 samples were generated with a feature matrix shape of (5000, 6).
Cell taking longer than 5 minutes?
The data generation runs 5,000 separate MC simulations, so it is genuinely computation-heavy. If it is taking much longer than expected, check that you have not accidentally increased n_paths beyond 50000 or n_samples beyond 5000.
help me debug why my training data generation is running slowly
Understand the moneyness feature
You may have noticed the line moneyness = S_arr / K_arr inside the function. This is a technique called feature engineering.
Why add moneyness as a feature?
Option prices depend heavily on how the spot price (S) relates to the strike price (K). An option with S=100 and K=50 behaves very differently from S=100 and K=150.
The moneyness ratio (S/K) encodes this relationship as a single number. A ratio above 1 means the option is "in the money" (profitable if exercised now). A ratio below 1 means it is "out of the money." This makes the nonlinear relationship between spot and strike easier for the model to learn directly.
Split into training and test sets
To evaluate whether the ML model genuinely learned the pricing function, you need data it has never seen during training. Splitting the dataset into separate training and test portions ensures an honest accuracy measurement.
- Create a new cell below Cell 4 and paste the following code:
# Cell 5: Train/Test Split
feature_names = ["S", "K", "T", "r", "sigma", "moneyness"]
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
print(f"Training set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print(f"Features: {feature_names}")What does this code do?
- feature_names lists the 6 columns in your feature matrix for later reference in plots.
- train_test_split shuffles the data and holds back 20% (1,000 samples) as a test set the model will never train on.
- random_state=42 fixes the shuffle so your results are reproducible every time you run the notebook.
- Run the cell. You should see output confirming 4000 training samples and 1000 test samples.
Seeing a NameError for X or y?
This means Cell 4 did not finish running successfully. Switch back to Cell 4 and re-run it. Wait for the "Done!" message before running Cell 5.
help me fix a NameError in my Jupyter notebook
✔️ Awesome, I've got everything!
Great. Double check you have saved your notebook file.
ⓧ I'd like to double check the full code
Here is the complete notebook up to this point:
# Cell 1: Imports and Setup
import numpy as np
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
import matplotlib.pyplot as plt
import time
from scipy.stats import norm
print(f"XGBoost version: {xgb.__version__}")
print("All imports successful!")
# Cell 2: Monte Carlo Pricing Engine
def monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type="call", seed=None):
if seed is not None:
rng = np.random.default_rng(seed)
else:
rng = np.random.default_rng()
Z = rng.standard_normal(n_paths)
S_T = S * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
if option_type == "call":
payoffs = np.maximum(S_T - K, 0)
else:
payoffs = np.maximum(K - S_T, 0)
price = np.exp(-r * T) * np.mean(payoffs)
return price
test_price = monte_carlo_price(S=100, K=100, T=1, r=0.05, sigma=0.2, n_paths=100000, seed=42)
print(f"MC Call Price (S=100, K=100, T=1, r=0.05, sigma=0.2): {test_price:.4f}")
print("Expected (Black-Scholes): ~10.45")
# Cell 3: Black-Scholes Validation
def black_scholes_price(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
else:
price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
return price
print("\nValidation: MC vs Black-Scholes")
print("-" * 50)
test_cases = [
(100, 100, 1.0, 0.05, 0.2, "call"),
(100, 110, 0.5, 0.03, 0.25, "call"),
(100, 90, 1.5, 0.08, 0.3, "put"),
(120, 100, 0.25, 0.04, 0.15, "call"),
]
for S, K, T, r, sigma, opt_type in test_cases:
bs = black_scholes_price(S, K, T, r, sigma, opt_type)
mc = monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type=opt_type, seed=42)
error = abs(mc - bs) / bs * 100
print(f" {opt_type.upper()} S={S} K={K} T={T}: BS={bs:.4f}, MC={mc:.4f}, Error={error:.2f}%")
# Cell 4: Generate Training Data
def generate_training_data(n_samples=5000, n_paths=50000, seed=42):
rng = np.random.default_rng(seed)
S_arr = rng.uniform(50, 150, n_samples)
K_arr = rng.uniform(50, 150, n_samples)
T_arr = rng.uniform(0.1, 2.0, n_samples)
r_arr = rng.uniform(0.01, 0.10, n_samples)
sigma_arr = rng.uniform(0.1, 0.5, n_samples)
prices = np.zeros(n_samples)
for i in range(n_samples):
prices[i] = monte_carlo_price(
S=S_arr[i], K=K_arr[i], T=T_arr[i],
r=r_arr[i], sigma=sigma_arr[i],
n_paths=n_paths, option_type="call", seed=seed + i
)
moneyness = S_arr / K_arr
X = np.column_stack([S_arr, K_arr, T_arr, r_arr, sigma_arr, moneyness])
y = prices
return X, y
print("Generating training data (5000 samples)... This may take 2-3 minutes.")
start = time.time()
X, y = generate_training_data(n_samples=5000, n_paths=50000, seed=42)
elapsed = time.time() - start
print(f"Done! Generated {len(y)} samples in {elapsed:.1f} seconds.")
print(f"Feature matrix shape: {X.shape}")
print(f"Price range: [{y.min():.2f}, {y.max():.2f}]")
# Cell 5: Train/Test Split
feature_names = ["S", "K", "T", "r", "sigma", "moneyness"]
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
print(f"Training set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print(f"Features: {feature_names}")Your training dataset is ready. Next up, you will feed these 4,000 samples into an XGBoost model and teach it to predict option prices without running a single simulation.
Train the XGBoost Surrogate Model
Your Monte Carlo engine is accurate, but pricing a single option requires simulating 50,000 random paths. Pricing thousands of options would take minutes or hours.
In this step, you will train an XGBoost regression model on your MC-generated training data. The model will learn the pricing function directly from the data, allowing it to predict option prices without running any simulations.
In this step, get ready to:
- Train an XGBoost surrogate model on your Monte Carlo training data.
- Evaluate model accuracy using RMSE and MAPE metrics.
- Visualize feature importance to understand what drives option prices.
Train the XGBoost regressor
XGBoost builds an ensemble of gradient-boosted decision trees. Each tree corrects the errors of the previous ones, making the ensemble progressively more accurate. For tabular data with a small number of features, this approach consistently matches or outperforms neural networks while training much faster.
- In your notebook, create a new cell (Cell 6) and add the model training code by pasting the following:
# Cell 6: Train XGBoost Model
model = xgb.XGBRegressor(
n_estimators=200,
max_depth=6,
learning_rate=0.1,
tree_method="hist",
objective="reg:squarederror",
random_state=42
)
print("Training XGBoost model...")
model.fit(X_train, y_train)
print("Training complete!")What do these hyperparameters mean?
- n_estimators=200 builds 200 decision trees in sequence. Each tree fixes mistakes from the one before it.
- max_depth=6 limits each tree to 6 levels of splits. This prevents overfitting while still capturing nonlinear patterns.
- learning_rate=0.1 controls how much each new tree contributes. Smaller values mean slower but more careful learning.
- tree_method="hist" uses histogram-based splitting for faster training on larger datasets.
- objective="reg:squarederror" tells XGBoost to minimize squared error, the standard loss function for regression.
Do not run this cell yet. You will add evaluation code to the same cell first.
Evaluate model accuracy
Training a model means nothing if you cannot measure how well it performs on data it has never seen. You will evaluate on the held-out test set using two metrics: RMSE (root mean squared error) and MAPE (mean absolute percentage error). Your target is a MAPE below 5%.
- Add this evaluation code directly below the training code in the same Cell 6:
# Evaluate on held-out test set
y_pred = model.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
mape = mean_absolute_percentage_error(y_test, y_pred) * 100
print(f"\nTest Set Performance:")
print(f" RMSE: {rmse:.4f}")
print(f" MAPE: {mape:.2f}%")
if mape < 5:
print(" Target met! MAPE is below 5%.")
else:
print(" Target not met. Consider increasing n_estimators.")What does this code do?
- model.predict(X_test) feeds the 1000 test samples through the trained XGBoost model to get predicted prices.
- RMSE tells you the average dollar error in absolute terms.
- MAPE tells you the average percentage error relative to the true price. A MAPE below 5% means the model is off by less than 5% on average.
- Run Cell 6 by clicking the play button or pressing Shift+Enter.
- You should see "Training complete!" followed by your test set metrics. The MAPE should be below 5%, confirming the model meets the accuracy target.
MAPE above 5%?
If your MAPE is slightly above 5%, try increasing n_estimators to 300 or reducing learning_rate to 0.05. Re-run the cell after making changes.
If you see a NameError for X_train or y_train, make sure you have run all previous cells (Cells 1-5) first.
Help me debug my XGBoost model training
Plot feature importance
XGBoost tracks how much each feature contributed to reducing prediction error across all 200 trees. Plotting this tells you which option parameters matter most for pricing. You would expect volatility (sigma) and moneyness to rank highest, since they have the strongest nonlinear impact on option value.
- Create a new cell (Cell 7) and add the feature importance visualization by pasting this code:
# Cell 7: Feature Importance
importances = model.feature_importances_
sorted_idx = np.argsort(importances)
plt.figure(figsize=(8, 5))
plt.barh(range(len(feature_names)), importances[sorted_idx])
plt.yticks(range(len(feature_names)), [feature_names[i] for i in sorted_idx])
plt.xlabel("Feature Importance (Gain)")
plt.title("XGBoost Feature Importance for Option Pricing")
plt.tight_layout()
plt.savefig("feature_importance.png", dpi=150)
plt.show()What does this code do?
- model.feature_importances_ returns an array of importance scores (one per feature) based on how much each feature improved predictions across all trees.
- np.argsort(importances) sorts the features from least to most important, so the horizontal bar chart reads naturally from bottom to top.
- plt.savefig("feature_importance.png", dpi=150) saves a publication-quality image to your project folder.
- Run Cell 7.
- You should see a horizontal bar chart with six bars. The features sigma and moneyness should appear near the top as the most important predictors.
Why do sigma and moneyness dominate?
Volatility (sigma) directly controls how widely the simulated stock price can move. Higher volatility means a wider distribution of terminal prices, which dramatically changes option value.
Moneyness (S/K) captures whether the option is in-the-money, at-the-money, or out-of-the-money. This single ratio encodes the most critical pricing signal more efficiently than raw S and K values alone.
Not seeing a chart?
If the plot does not appear inline, make sure you are running in a Jupyter notebook environment (not a plain Python script). Check that import matplotlib.pyplot as plt ran successfully in Cell 1.
If you see a NameError for model, go back and run Cell 6 first. Cell 7 depends on the trained model from the previous cell.
Help me troubleshoot my feature importance plot
✔️ Awesome, I've got everything!
Great. Double check you have saved your notebook with both Cell 6 and Cell 7 executed.
ⓧ I'd like to double check the full code
Here is the complete code for both cells added in this step:
# Cell 6: Train XGBoost Model
model = xgb.XGBRegressor(
n_estimators=200,
max_depth=6,
learning_rate=0.1,
tree_method="hist",
objective="reg:squarederror",
random_state=42
)
print("Training XGBoost model...")
model.fit(X_train, y_train)
print("Training complete!")
# Evaluate
y_pred = model.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
mape = mean_absolute_percentage_error(y_test, y_pred) * 100
print(f"\nTest Set Performance:")
print(f" RMSE: {rmse:.4f}")
print(f" MAPE: {mape:.2f}%")
if mape < 5:
print(" Target met! MAPE is below 5%.")
else:
print(" Target not met. Consider increasing n_estimators.")# Cell 7: Feature Importance
importances = model.feature_importances_
sorted_idx = np.argsort(importances)
plt.figure(figsize=(8, 5))
plt.barh(range(len(feature_names)), importances[sorted_idx])
plt.yticks(range(len(feature_names)), [feature_names[i] for i in sorted_idx])
plt.xlabel("Feature Importance (Gain)")
plt.title("XGBoost Feature Importance for Option Pricing")
plt.tight_layout()
plt.savefig("feature_importance.png", dpi=150)
plt.show()Your surrogate model is trained and verified. Next up, you will put it to the ultimate test: a head-to-head speed and accuracy benchmark against the full Monte Carlo engine.
Benchmark Speed and Accuracy
Your XGBoost model is trained and achieving less than 5% MAPE on the test set. The feature importance plot confirms sigma and moneyness drive predictions the most. But accuracy alone does not justify replacing a simulation with an ML model.
The whole point of a surrogate model is speed. In this step, you will measure exactly how much faster the XGBoost model is compared to the Monte Carlo engine, then visualize accuracy and error distributions so you can prove the tradeoff is worth it.
In this step, get ready to:
- Run a head-to-head speed benchmark timing 1000 MC prices vs 1000 ML predictions.
- Create an accuracy scatter plot showing how closely ML predictions match MC prices.
- Create an error distribution histogram proving most predictions fall within a tight error band.
Time MC pricing vs ML inference
To prove the speedup is real, you need to time both approaches on the same workload. You will generate 1000 random option parameter sets, price all of them with the Monte Carlo engine (50,000 paths each), then predict all of them with the trained XGBoost model.
- Create a new cell in your notebook (Cell 8) and add the speed benchmark code:
# Cell 8: Speed Benchmark
n_benchmark = 1000
benchmark_params = np.random.default_rng(99).uniform(
low=[50, 50, 0.1, 0.01, 0.1],
high=[150, 150, 2.0, 0.10, 0.5],
size=(n_benchmark, 5)
)
benchmark_moneyness = benchmark_params[:, 0] / benchmark_params[:, 1]
benchmark_X = np.column_stack([benchmark_params, benchmark_moneyness])
# Time MC pricing
print(f"Benchmarking {n_benchmark} option prices...")
start_mc = time.time()
mc_prices = np.array([
monte_carlo_price(
S=benchmark_params[i, 0], K=benchmark_params[i, 1],
T=benchmark_params[i, 2], r=benchmark_params[i, 3],
sigma=benchmark_params[i, 4], n_paths=50000
)
for i in range(n_benchmark)
])
time_mc = time.time() - start_mc
# Time ML pricing
start_ml = time.time()
ml_prices = model.predict(benchmark_X)
time_ml = time.time() - start_ml
speedup = time_mc / time_ml
print(f"\nMonte Carlo: {time_mc:.2f} seconds for {n_benchmark} prices")
print(f"XGBoost ML: {time_ml:.4f} seconds for {n_benchmark} prices")
print(f"Speedup: {speedup:.0f}x faster")What does this code do?
The benchmark generates 1000 random option parameter sets using np.random.default_rng(99) for reproducibility. Each set has S, K, T, r, and sigma drawn from the same ranges as the training data.
The MC block loops through all 1000 options, running 50,000 simulated paths per price. The ML block calls model.predict() once on the entire batch. The speedup ratio divides MC time by ML time to show how many times faster the surrogate model is.
- Run the cell. This will take 1-2 minutes while the MC engine prices all 1000 options.
- Check the printed output. You should see the MC engine taking tens of seconds while XGBoost finishes in milliseconds, with a speedup of at least 100x.
Speedup lower than 100x?
If your speedup is below 100x, your machine may be running background tasks that slow the timing. Close other applications and rerun the cell. The MC loop is CPU-bound, so heavy background processes affect it directly.
help me understand why my speed benchmark shows a lower-than-expected speedup ratio
Plot the accuracy scatter
A speedup means nothing if the predictions are wrong. The accuracy scatter plot places MC prices on the x-axis and ML predictions on the y-axis. If the model is accurate, every point sits on the diagonal line where predicted equals actual.
- Create a new cell (Cell 9) and add the accuracy scatter plot code:
# Cell 9: Accuracy Scatter Plot
plt.figure(figsize=(8, 6))
plt.scatter(mc_prices, ml_prices, alpha=0.4, s=10, color="steelblue")
min_val = min(mc_prices.min(), ml_prices.min())
max_val = max(mc_prices.max(), ml_prices.max())
plt.plot([min_val, max_val], [min_val, max_val], "r--", linewidth=2, label="Perfect prediction")
plt.xlabel("Monte Carlo Price ($)")
plt.ylabel("XGBoost Predicted Price ($)")
plt.title("ML Surrogate vs Monte Carlo: Option Pricing Accuracy")
plt.legend()
plt.tight_layout()
plt.savefig("accuracy_scatter.png", dpi=150)
plt.show()What does this code do?
The scatter plot uses the 1000 benchmark prices from the previous cell. Each dot represents one option. The red dashed line is the "perfect prediction" diagonal where predicted price equals MC price exactly.
Tight clustering around this diagonal means the ML model learned the pricing function well. Dots far from the line would indicate systematic over- or under-pricing for certain parameter regions.
- Run the cell.
- You should see a scatter plot where nearly all points sit tightly along the red dashed diagonal. The chart is saved as accuracy_scatter.png in your project folder.
Points scattered far from the diagonal?
If you see significant scatter away from the diagonal, your model may not have trained correctly. Return to Cell 6 and verify the MAPE is below 5%. If MAPE is high, increase n_estimators to 300 and retrain.
help me debug why my accuracy scatter plot shows poor ML predictions
Visualize the error distribution
The scatter plot shows overall fit, but you also want to know the shape of your errors. A histogram of percentage errors reveals whether the model has a bias (shifted mean) or fat tails (extreme outliers).
- Create a new cell (Cell 10) and add the error distribution code:
# Cell 10: Error Distribution
pct_errors = (ml_prices - mc_prices) / mc_prices * 100
# Filter out extreme outliers from near-zero prices
mask = mc_prices > 0.5
pct_errors_filtered = pct_errors[mask]
plt.figure(figsize=(8, 5))
plt.hist(pct_errors_filtered, bins=50, edgecolor="black", alpha=0.7, color="steelblue")
plt.axvline(x=0, color="red", linestyle="--", linewidth=2)
plt.xlabel("Percentage Error (%)")
plt.ylabel("Frequency")
plt.title(f"Error Distribution (Mean={pct_errors_filtered.mean():.2f}%, Std={pct_errors_filtered.std():.2f}%)")
plt.tight_layout()
plt.savefig("error_distribution.png", dpi=150)
plt.show()
print(f"\nError Statistics (excluding near-zero prices):")
print(f" Mean Error: {pct_errors_filtered.mean():.2f}%")
print(f" Std Error: {pct_errors_filtered.std():.2f}%")
print(f" Max Error: {np.abs(pct_errors_filtered).max():.2f}%")Why filter out near-zero prices?
Options that are deep out-of-the-money have MC prices close to zero (e.g., $0.02). A tiny absolute error like $0.01 becomes a 50% relative error when the true price is near zero. Filtering with mc_prices > 0.5 removes these misleading outliers so the histogram shows meaningful error patterns.
- Run the cell.
- You should see a bell-shaped histogram centered near 0% with most errors falling within +/- 5%. The standard deviation should be below 3%. The chart is saved as error_distribution.png in your project folder.
Histogram not centered at zero?
A shifted histogram indicates systematic bias. Check that you ran the training data generation (Cell 4) with the same seed as the benchmark. Also confirm you did not accidentally retrain the model between generating training data and running this cell.
help me understand why my error distribution is biased away from zero
✔️ Awesome, I've got everything!
Great. Double check you have saved all three cells and can see accuracy_scatter.png and error_distribution.png in your project folder.
ⓧ I'd like to double check the full code
Here is the complete notebook code up to this point (Cells 1-10):
# Cell 1: Imports and Setup
import numpy as np
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
import matplotlib.pyplot as plt
import time
from scipy.stats import norm
print(f"XGBoost version: {xgb.__version__}")
print("All imports successful!")
# Cell 2: Monte Carlo Pricing Engine
def monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type="call", seed=None):
if seed is not None:
rng = np.random.default_rng(seed)
else:
rng = np.random.default_rng()
Z = rng.standard_normal(n_paths)
S_T = S * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
if option_type == "call":
payoffs = np.maximum(S_T - K, 0)
else:
payoffs = np.maximum(K - S_T, 0)
price = np.exp(-r * T) * np.mean(payoffs)
return price
test_price = monte_carlo_price(S=100, K=100, T=1, r=0.05, sigma=0.2, n_paths=100000, seed=42)
print(f"MC Call Price (S=100, K=100, T=1, r=0.05, sigma=0.2): {test_price:.4f}")
print("Expected (Black-Scholes): ~10.45")
# Cell 3: Black-Scholes Validation
def black_scholes_price(S, K, T, r, sigma, option_type="call"):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == "call":
price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
else:
price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
return price
print("\nValidation: MC vs Black-Scholes")
print("-" * 50)
test_cases = [
(100, 100, 1.0, 0.05, 0.2, "call"),
(100, 110, 0.5, 0.03, 0.25, "call"),
(100, 90, 1.5, 0.08, 0.3, "put"),
(120, 100, 0.25, 0.04, 0.15, "call"),
]
for S, K, T, r, sigma, opt_type in test_cases:
bs = black_scholes_price(S, K, T, r, sigma, opt_type)
mc = monte_carlo_price(S, K, T, r, sigma, n_paths=100000, option_type=opt_type, seed=42)
error = abs(mc - bs) / bs * 100
print(f" {opt_type.upper()} S={S} K={K} T={T}: BS={bs:.4f}, MC={mc:.4f}, Error={error:.2f}%")
# Cell 4: Generate Training Data
def generate_training_data(n_samples=5000, n_paths=50000, seed=42):
rng = np.random.default_rng(seed)
S_arr = rng.uniform(50, 150, n_samples)
K_arr = rng.uniform(50, 150, n_samples)
T_arr = rng.uniform(0.1, 2.0, n_samples)
r_arr = rng.uniform(0.01, 0.10, n_samples)
sigma_arr = rng.uniform(0.1, 0.5, n_samples)
prices = np.zeros(n_samples)
for i in range(n_samples):
prices[i] = monte_carlo_price(
S=S_arr[i], K=K_arr[i], T=T_arr[i],
r=r_arr[i], sigma=sigma_arr[i],
n_paths=n_paths, option_type="call", seed=seed + i
)
moneyness = S_arr / K_arr
X = np.column_stack([S_arr, K_arr, T_arr, r_arr, sigma_arr, moneyness])
y = prices
return X, y
print("Generating training data (5000 samples)... This may take 2-3 minutes.")
start = time.time()
X, y = generate_training_data(n_samples=5000, n_paths=50000, seed=42)
elapsed = time.time() - start
print(f"Done! Generated {len(y)} samples in {elapsed:.1f} seconds.")
print(f"Feature matrix shape: {X.shape}")
print(f"Price range: [{y.min():.2f}, {y.max():.2f}]")
# Cell 5: Train/Test Split
feature_names = ["S", "K", "T", "r", "sigma", "moneyness"]
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)
print(f"Training set: {X_train.shape[0]} samples")
print(f"Test set: {X_test.shape[0]} samples")
print(f"Features: {feature_names}")
# Cell 6: Train XGBoost Model
model = xgb.XGBRegressor(
n_estimators=200,
max_depth=6,
learning_rate=0.1,
tree_method="hist",
objective="reg:squarederror",
random_state=42
)
print("Training XGBoost model...")
model.fit(X_train, y_train)
print("Training complete!")
y_pred = model.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
mape = mean_absolute_percentage_error(y_test, y_pred) * 100
print(f"\nTest Set Performance:")
print(f" RMSE: {rmse:.4f}")
print(f" MAPE: {mape:.2f}%")
if mape < 5:
print(" Target met! MAPE is below 5%.")
else:
print(" Target not met. Consider increasing n_estimators.")
# Cell 7: Feature Importance
importances = model.feature_importances_
sorted_idx = np.argsort(importances)
plt.figure(figsize=(8, 5))
plt.barh(range(len(feature_names)), importances[sorted_idx])
plt.yticks(range(len(feature_names)), [feature_names[i] for i in sorted_idx])
plt.xlabel("Feature Importance (Gain)")
plt.title("XGBoost Feature Importance for Option Pricing")
plt.tight_layout()
plt.savefig("feature_importance.png", dpi=150)
plt.show()
# Cell 8: Speed Benchmark
n_benchmark = 1000
benchmark_params = np.random.default_rng(99).uniform(
low=[50, 50, 0.1, 0.01, 0.1],
high=[150, 150, 2.0, 0.10, 0.5],
size=(n_benchmark, 5)
)
benchmark_moneyness = benchmark_params[:, 0] / benchmark_params[:, 1]
benchmark_X = np.column_stack([benchmark_params, benchmark_moneyness])
print(f"Benchmarking {n_benchmark} option prices...")
start_mc = time.time()
mc_prices = np.array([
monte_carlo_price(
S=benchmark_params[i, 0], K=benchmark_params[i, 1],
T=benchmark_params[i, 2], r=benchmark_params[i, 3],
sigma=benchmark_params[i, 4], n_paths=50000
)
for i in range(n_benchmark)
])
time_mc = time.time() - start_mc
start_ml = time.time()
ml_prices = model.predict(benchmark_X)
time_ml = time.time() - start_ml
speedup = time_mc / time_ml
print(f"\nMonte Carlo: {time_mc:.2f} seconds for {n_benchmark} prices")
print(f"XGBoost ML: {time_ml:.4f} seconds for {n_benchmark} prices")
print(f"Speedup: {speedup:.0f}x faster")
# Cell 9: Accuracy Scatter Plot
plt.figure(figsize=(8, 6))
plt.scatter(mc_prices, ml_prices, alpha=0.4, s=10, color="steelblue")
min_val = min(mc_prices.min(), ml_prices.min())
max_val = max(mc_prices.max(), ml_prices.max())
plt.plot([min_val, max_val], [min_val, max_val], "r--", linewidth=2, label="Perfect prediction")
plt.xlabel("Monte Carlo Price ($)")
plt.ylabel("XGBoost Predicted Price ($)")
plt.title("ML Surrogate vs Monte Carlo: Option Pricing Accuracy")
plt.legend()
plt.tight_layout()
plt.savefig("accuracy_scatter.png", dpi=150)
plt.show()
# Cell 10: Error Distribution
pct_errors = (ml_prices - mc_prices) / mc_prices * 100
mask = mc_prices > 0.5
pct_errors_filtered = pct_errors[mask]
plt.figure(figsize=(8, 5))
plt.hist(pct_errors_filtered, bins=50, edgecolor="black", alpha=0.7, color="steelblue")
plt.axvline(x=0, color="red", linestyle="--", linewidth=2)
plt.xlabel("Percentage Error (%)")
plt.ylabel("Frequency")
plt.title(f"Error Distribution (Mean={pct_errors_filtered.mean():.2f}%, Std={pct_errors_filtered.std():.2f}%)")
plt.tight_layout()
plt.savefig("error_distribution.png", dpi=150)
plt.show()
print(f"\nError Statistics (excluding near-zero prices):")
print(f" Mean Error: {pct_errors_filtered.mean():.2f}%")
print(f" Std Error: {pct_errors_filtered.std():.2f}%")
print(f" Max Error: {np.abs(pct_errors_filtered).max():.2f}%")You now have hard evidence that your ML surrogate model delivers the same pricing accuracy as Monte Carlo simulation at a fraction of the computational cost. Next up, you will take this even further in the Secret Mission by computing option Greeks using your trained model.
Secret mission
Estimate Greeks with Your ML Model
Your ML model can predict prices instantly. But can it also predict how sensitive those prices are to changes in the underlying parameters? In this secret mission, you'll compute option Greeks (Delta and Vega) using both your MC engine and your ML surrogate, proving the model captures the pricing surface well enough to compute derivatives.
Clean Up Your Resources
Clean Up Your Resources
This project runs entirely locally, so there are no ongoing costs or cloud resources to worry about. Decide whether to keep your project files for future experimentation, or delete them to free up disk space.
Resources you used:
- mc-ml-pricer project folder (contains your notebook, virtual environment, and generated chart images).
Keep everything running
No action needed. Choose this if you plan to keep experimenting with your Monte Carlo engine or XGBoost model.
- Your mc-ml-pricer folder stays on your Desktop with the notebook, virtual environment, and all generated charts.
- The virtual environment (.venv) holds all your installed packages (NumPy, XGBoost, scikit-learn, matplotlib, Jupyter).
- To resume work later, open the folder in VS Code and activate the virtual environment again.
Pause - I'll come back to this later
Close VS Code to free up memory. Your project files remain on disk and are ready to pick back up any time.
- Close VS Code.
- When you return, open the mc-ml-pricer folder in VS Code and reactivate your virtual environment with .venv\Scripts\activate (Windows) or source .venv/bin/activate (macOS/Linux).
Delete - I don't want to use this again
Remove all project files, the virtual environment, and generated charts permanently.
- Close VS Code.
- Delete the project folder by running this command in your terminal.
macOS
cd ~/Desktop
rm -rf mc-ml-pricerWindows
cd ~\Desktop
rmdir /s /q mc-ml-pricer- Confirm the folder is gone by running ls (macOS) or dir (Windows) on your Desktop. You should no longer see mc-ml-pricer in the listing.
That's a Wrap!
That's a Wrap!
Nice work! You just built a complete ML surrogate pipeline that replaces a computationally expensive Monte Carlo simulation with near-instant predictions.
This "simulation-to-ML" pattern is the same technique used in production systems at quant firms, autonomous vehicle companies, and physics simulations.
You've learned how to:
- Build a vectorized Monte Carlo pricing engine for European options using Geometric Brownian Motion, validated against Black-Scholes analytical prices.
- Train an XGBoost surrogate model on synthetic training data, achieving less than 5% MAPE on a held-out test set.
- Benchmark ML inference speed at 100-1000x faster than full Monte Carlo simulation, with publication-ready visualizations proving accuracy.
- Secret Mission: Estimated option Greeks (Delta and Vega) using both MC finite differences and your ML model, proving the surrogate captures the pricing surface well enough to compute derivatives.
Ready to quiz yourself?