Python

Black-Scholes Options Pricing in Python

The Black-Scholes model is the foundation of modern options pricing. Published in 1973, it gives a theoretical price for European-style options based on five inputs: stock price, strike price, time to expiry, risk-free rate, and volatility. In this tutorial, we implement Black-Scholes from scratch and calculate the Greeks (Delta, Gamma, Theta, Vega).

Prerequisites

✓Basic Python
✓Understanding of options (calls and puts)
✓Basic statistics (normal distribution)

Step 1.Import dependencies

We only need NumPy for math and scipy.stats for the normal distribution CDF.

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

Step 2.Implement Black-Scholes

d1 and d2 are intermediate values. The call price is the probability-weighted expected payoff discounted to present value.

def black_scholes(S, K, T, r, sigma, option_type="call"):
    """Calculate Black-Scholes option price.
    S: Current stock price
    K: Strike price
    T: Time to expiry (years)
    r: Risk-free rate
    sigma: Volatility
    """
    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, d1, d2

# Example: AAPL at $195, strike $200, 30 days, 5% rate, 25% vol
call_price, d1, d2 = black_scholes(195, 200, 30/365, 0.05, 0.25, "call")
put_price, _, _ = black_scholes(195, 200, 30/365, 0.05, 0.25, "put")

print(f"Call price: ${call_price:.2f}")
print(f"Put price:  ${put_price:.2f}")

Step 3.Calculate the Greeks

Delta measures price sensitivity, Gamma is the rate of change of Delta, Theta is time decay, and Vega is volatility sensitivity.

def greeks(S, K, T, r, sigma):
    _, d1, d2 = black_scholes(S, K, T, r, sigma)

    delta_call = norm.cdf(d1)
    delta_put = delta_call - 1
    gamma = norm.pdf(d1) / (S * sigma * np.sqrt(T))
    theta_call = (-S * norm.pdf(d1) * sigma / (2 * np.sqrt(T))
                  - r * K * np.exp(-r * T) * norm.cdf(d2)) / 365
    vega = S * norm.pdf(d1) * np.sqrt(T) / 100

    return {
        "Delta (Call)": round(delta_call, 4),
        "Delta (Put)": round(delta_put, 4),
        "Gamma": round(gamma, 4),
        "Theta (Call)": round(theta_call, 4),
        "Vega": round(vega, 4)
    }

g = greeks(195, 200, 30/365, 0.05, 0.25)
for k, v in g.items():
    print(f"{k}: {v}")

Step 4.Plot option price vs stock price

The curved line shows the option premium. It always lies above the intrinsic value due to time value.

stock_prices = np.linspace(150, 250, 100)
call_prices = [black_scholes(s, 200, 30/365, 0.05, 0.25, "call")[0] for s in stock_prices]
put_prices = [black_scholes(s, 200, 30/365, 0.05, 0.25, "put")[0] for s in stock_prices]
intrinsic_call = [max(s - 200, 0) for s in stock_prices]

plt.figure(figsize=(12, 6))
plt.plot(stock_prices, call_prices, label="Call Price", color="cyan")
plt.plot(stock_prices, put_prices, label="Put Price", color="orange")
plt.plot(stock_prices, intrinsic_call, "--", label="Intrinsic Value", color="gray", alpha=0.5)
plt.axvline(200, color="white", linestyle=":", alpha=0.3, label="Strike")
plt.xlabel("Stock Price")
plt.ylabel("Option Price")
plt.title("Black-Scholes Option Pricing")
plt.legend()
plt.savefig("black_scholes.png", dpi=150)
plt.show()

Expected Output

Call price: $3.87
Put price:  $8.05
Delta (Call): 0.3842
Delta (Put): -0.6158
Gamma: 0.0285
Theta (Call): -0.1245
Vega: 0.2214

Next Steps

→Build an implied volatility calculator
→Price American options using binomial trees
→Create a volatility surface plotter