TL;DR: Black-Scholes assumes you can hold the underlying. You can't hold 1,000 barrels of crude oil under your desk. Black-76 fixes this by pricing off the futures price directly. I built a full Monte Carlo pricer in Python with generator-based path simulation, European and Asian payoffs, Greeks via finite differences, and live WTI market data. Here's how it works.
Why Black-Scholes Is Wrong for Commodities
Every quant finance tutorial starts with Black-Scholes. It's elegant, Nobel Prize-winning, and wrong for most commodity derivatives.
Black-Scholes assumes you can buy and hold the underlying asset to hedge your position. That works for equities. Buy 100 shares of Apple and stick them in your brokerage account. But commodities? You're pricing an option on WTI crude oil. You can't store 1,000 barrels in your garage. Storage costs, convenience yields, and delivery logistics make the carry-and-hold assumption break down completely.
Fischer Black figured this out in 1976. His insight: commodity options trade on futures, not spot prices. The futures price already bakes in storage costs, convenience yields, and interest rates. So instead of modelling the spot price with a drift of (r - q), you model the futures price with zero drift, because futures are already forward-looking.
This is the Black-76 model, and it's what professionals actually use for energy, agriculture, and metals derivatives.
The Key Difference: Drift Terms
Here's where most implementations get it wrong. In GBM simulation under Black-76:
# Black-Scholes drift (WRONG for commodities)
drift = (r - 0.5 * sigma**2) * dt
# Black-76 drift (CORRECT for futures)
drift = -0.5 * sigma**2 * dt
The risk-free rate r disappears from the simulation entirely. Why? Because futures prices are already risk-neutral expectations of the future spot price. The r only appears later, when you discount the expected payoff back to today.
This is the single most common bug in commodity option pricers. Use the wrong drift and your prices will systematically overshoot. My test suite catches this explicitly. The mean drift test validates that simulated paths don't drift upward when they shouldn't.
Architecture: Generators, Not Arrays
A naive Monte Carlo pricer generates 100,000 paths, stores them all in a NumPy array, then computes payoffs. That works, but it's wasteful. You never need all paths in memory simultaneously.
I used Python generators instead:
def generate_paths(self, n_paths: int):
"""Yield one simulated price path at a time."""
dt = self.contract.expiry / self.n_steps
drift = -0.5 * self.market.volatility**2 * dt
vol_sqrt_dt = self.market.volatility * np.sqrt(dt)
for _ in range(n_paths):
path = np.empty(self.n_steps + 1)
path[0] = self.market.futures_price
z = np.random.standard_normal(self.n_steps)
for t in range(self.n_steps):
path[t + 1] = path[t] * np.exp(drift + vol_sqrt_dt * z[t])
yield path
Each path is generated, consumed for its payoff, then garbage collected. Memory usage stays constant regardless of how many paths you simulate. The yield keyword turns this into a lazy generator, so the calling code pulls one path at a time through a simple for loop.
The price() method (already written as framework code) orchestrates the whole thing:
def price(self, n_paths: int = 100_000) -> dict:
payoffs = np.array([
self.payoff(path) for path in self.generate_paths(n_paths)
])
discounted = np.exp(-self.market.rate * self.contract.expiry)
price = float(np.mean(payoffs) * discounted)
std_error = float(np.std(payoffs) * discounted / np.sqrt(n_paths))
return {"price": price, "std_error": std_error}
Generate paths. Compute payoffs. Average. Discount. That's Monte Carlo pricing.
European vs Asian Payoffs
European options are simple. Only the final price matters:
def payoff(self, path: np.ndarray) -> float:
terminal = path[-1]
if self.contract.option_type == "call":
return max(terminal - self.contract.strike, 0.0)
return max(self.contract.strike - terminal, 0.0)
But in commodity markets, Asian options dominate. An oil producer hedging their quarterly revenue doesn't care about the price on one specific day. They care about the average price over the period. Asian options pay out based on the arithmetic mean of the path:
def payoff(self, path: np.ndarray) -> float:
average = np.mean(path[1:]) # exclude initial price
if self.contract.option_type == "call":
return max(average - self.contract.strike, 0.0)
return max(self.contract.strike - average, 0.0)
Asian options are cheaper than European options (averaging reduces volatility exposure) and have no closed-form solution, so Monte Carlo is the natural pricing method. This is why the engine architecture matters: the same path generation feeds both payoff types through a clean abstraction.
Greeks via Finite Differences
Greeks measure how the option price changes when inputs move. Delta tells you how much the price moves per $1 change in the underlying. Vega tells you the sensitivity to volatility. In production, these drive hedging decisions.
The analytical approach requires deriving partial derivatives of the pricing formula. The Monte Carlo approach is simpler: bump and reprice.
def compute_greeks(engine, bump_pct=0.01):
base = engine.price()["price"]
# Delta: bump futures price up, reprice, compute slope
# Vega: bump volatility up, reprice, compute slope
# Theta: reduce time to expiry, reprice, compute slope
# ...
Bump each input by a small amount, reprice, and compute the finite difference. The key optimisation is Common Random Numbers (CRN): seed the random number generator identically for both the base and bumped prices. This ensures the difference captures only the input change, not random noise, dramatically reducing variance in the Greek estimates.
Live Market Data
The pricer pulls real WTI crude futures prices and implied volatility from the OVX index via yfinance:
$ python cli.py live --strike 72.50 --expiry 0.5
# Fetches CL=F (WTI futures) and ^OVX (oil volatility index)
# Prices the option with actual market data
No manual data entry. No stale inputs. The CLI fetches current prices, feeds them into the engine, and compares the Monte Carlo result against the Black-76 analytical solution. In testing, the Monte Carlo price consistently lands within 1-2% of the analytical answer with 100,000 paths.
What This Doesn't Do
Honesty about limitations matters more than feature lists:
- Flat volatility only. Real energy markets have a volatility smile/skew. A SABR stochastic volatility model would handle this, but adds significant complexity.
- No American exercise. Early exercise requires Longstaff-Schwartz least-squares Monte Carlo, which is a different beast entirely.
- Single underlying. No spread options, no crack spread modelling.
- Best for ATM options. Deep OTM/ITM options need more paths for convergence. The standard error output tells you when you need more.
These are deliberate scoping decisions, not oversights. The pricer does one thing well: European and Asian options on commodity futures using Black-76.
Try It
The full source code is on GitHub. 29 tests, CI via GitHub Actions, documented CLI.
I also built an interactive Black-76 calculator where you can plug in your own inputs and get instant analytical prices and Greeks, right in the browser.
git clone https://github.com/Z-Mahmood/wti-options-pricer.git
cd wti-options-pricer
pip install -e .
python cli.py price --futures-price 72.50 --strike 72.50 \
--volatility 0.30 --rate 0.05 --expiry 0.5