Quant Math, Step by Step: From Returns to Option Pricing (with Python)

Every formula below comes with the reasoning behind it, not just the code — the goal is to see why each step exists before you run it.

In this post:

1. Returns & Log Returns — the raw material
2. Volatility & Annualization
3. Geometric Brownian Motion — simulating a stock path
4. Black–Scholes — pricing a European option, derived step by step
5. Monte Carlo Option Pricing — a model-free check
6. Markowitz Portfolio Optimization
7. Value at Risk (VaR)

1. Returns & Log Returns

Everything in quant finance starts with turning a price series into a returns series, because prices themselves aren't stationary (their mean and variance drift over time), but returns behave much closer to a stable statistical process.

Step 1 — Simple return. $$ R_t = \frac{P_t - P_{t-1}}{P_{t-1}} $$ This is just "percent change." Intuitive, but simple returns don't add up cleanly across time.

Step 2 — Log return. $$ r_t = \ln\left(\frac{P_t}{P_{t-1}}\right) $$

Why log returns instead of simple returns? Log returns are time-additive: the return over two days is just $r_1 + r_2$, not a product. They're also symmetric (a 10% gain and a 10% loss aren't mirror images in simple-return space, but they are in log space) and they're the natural output of continuous-time models like Geometric Brownian Motion, which we use later for option pricing.

python

import numpy as np
import pandas as pd

# Suppose 'prices' is a pandas Series of daily closing prices
prices = pd.Series([100, 101.5, 99.8, 102.3, 103.0, 101.9])

simple_returns = prices.pct_change().dropna()
log_returns = np.log(prices / prices.shift(1)).dropna()

print("Simple returns:\n", simple_returns.values)
print("Log returns:\n", log_returns.values)

Simple returns: [ 0.015 -0.01675 0.02505 0.00684 -0.01068] Log returns: [ 0.01489 -0.01689 0.02474 0.00681 -0.01074]

Notice the two series are close but not identical — that gap grows with the size of the move, which is exactly the asymmetry log returns are designed to fix.

2. Volatility & Annualization

Step 1 — Daily volatility is just the standard deviation of daily log returns: $$ \sigma_{\text{daily}} = \sqrt{\frac{1}{n-1}\sum_{t=1}^{n}(r_t - \bar r)^2} $$

Step 2 — Scale to annual volatility. $$ \sigma_{\text{annual}} = \sigma_{\text{daily}} \times \sqrt{252} $$

Why multiply by $\sqrt{252}$ and not 252? Under the standard assumption that daily returns are independent and identically distributed, variances (not standard deviations) add up across independent days: $\text{Var}(\text{sum}) = n \times \text{Var}(\text{one day})$. Since volatility is a standard deviation — the square root of variance — scaling it to $n$ days means multiplying by $\sqrt{n}$. 252 is used because that's roughly the number of trading days in a year.

python

daily_vol = log_returns.std()
annual_vol = daily_vol * np.sqrt(252)

print(f"Daily volatility:  {daily_vol:.4%}")
print(f"Annualized volatility: {annual_vol:.2%}")

Daily volatility: 1.7534% Annualized volatility: 27.83%

3. Geometric Brownian Motion (GBM)

GBM is the standard model for "how does a stock price evolve," and it's the engine underneath both Black–Scholes and most Monte Carlo pricing.

Step 1 — The stochastic differential equation. $$ dS_t = \mu S_t\,dt + \sigma S_t\,dW_t $$ $S_t$ is price, $\mu$ is the expected (drift) return, $\sigma$ is volatility, and $dW_t$ is a random Brownian increment (the source of randomness).

Why this particular form? Both the drift term and the random term are proportional to $S_t$ itself. That encodes a simple economic idea: a stock that costs $200 should move about twice as much in dollar terms as one that costs $100, for the same percentage volatility. This is what makes GBM produce percentage-based, not dollar-based, randomness — matching how we actually think about returns.

Step 2 — Solve the SDE (via Itô's lemma) to get a closed-form expression you can simulate directly: $$ S_t = S_0 \exp\left[\left(\mu - \tfrac{1}{2}\sigma^2\right)t + \sigma W_t\right] $$

Why the $-\tfrac{1}{2}\sigma^2$ term? This is the most commonly misunderstood piece of GBM. Naively you might expect the exponent to just be $\mu t + \sigma W_t$. But because of Itô's lemma (calculus for random processes), the act of exponentiating a random variable introduces convexity — Jensen's inequality means $E[e^X] > e^{E[X]}$. The $-\tfrac12\sigma^2 t$ term exists specifically to correct for that convexity bias so that the arithmetic mean growth rate of the simulated price really is $\mu$.

python

import numpy as np
import matplotlib.pyplot as plt

def simulate_gbm_paths(S0, mu, sigma, T, n_steps, n_paths, seed=42):
    """
    S0      : starting price
    mu      : expected annual return (drift)
    sigma   : annual volatility
    T       : time horizon in years
    n_steps : number of time steps
    n_paths : number of simulated price paths
    """
    rng = np.random.default_rng(seed)
    dt = T / n_steps

    # Step 1: draw standard normal shocks for each step/path
    Z = rng.standard_normal((n_paths, n_steps))

    # Step 2: convert to log-return increments using the closed-form GBM solution
    increments = (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z

    # Step 3: cumulative sum of log returns = log price path, then exponentiate
    log_paths = np.log(S0) + np.cumsum(increments, axis=1)
    paths = np.exp(log_paths)

    # prepend S0 as t=0
    paths = np.hstack([np.full((n_paths, 1), S0), paths])
    return paths

paths = simulate_gbm_paths(S0=100, mu=0.08, sigma=0.25, T=1, n_steps=252, n_paths=5)
print("Ending prices for 5 simulated paths:", paths[:, -1].round(2))

Ending prices for 5 simulated paths: [131.42 89.7 112.85 97.33 141.06]

Five independent futures for the same stock — same drift and volatility assumptions, wildly different outcomes. That spread is the risk that option prices are compensating for.

4. Black–Scholes: Pricing a European Call, Step by Step

Black–Scholes answers: "given a stock following GBM, what's the fair price of the right (not obligation) to buy it at a fixed strike $K$ on a future date $T$?"

Step 1 — Set up the hedging argument. The core insight (Black, Scholes, and Merton, 1973) is that you can build a portfolio of the option plus a continuously rebalanced position in the underlying stock that is instantaneously risk-free. Because it's risk-free, it must earn exactly the risk-free rate $r$ — otherwise there would be a riskless arbitrage.

Step 2 — Turn that into a PDE. Applying Itô's lemma to the option value $V(S,t)$ and requiring the hedged portfolio to earn rate $r$ produces the Black–Scholes partial differential equation: $$ \frac{\partial V}{\partial t} + \tfrac12 \sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0 $$

Step 3 — Solve the PDE for a European call with the boundary condition $V(S,T) = \max(S-K, 0)$. The closed-form solution is: $$ C = S_0 N(d_1) - K e^{-rT} N(d_2) $$ $$ d_1 = \frac{\ln(S_0/K) + (r + \sigma^2/2)T}{\sigma\sqrt{T}}, \qquad d_2 = d_1 - \sigma\sqrt{T} $$ where $N(\cdot)$ is the standard normal CDF.

Why does the formula look like this? Read it as "$N(d_1)$: the probability-weighted stock you'll receive" minus "$N(d_2)$: the probability-weighted strike you'll pay" — both discounted appropriately. $N(d_2)$ is (under the risk-neutral measure) the actual probability the option finishes in-the-money. $N(d_1)$ plays a similar role but weighted by the stock price itself (it's also the option's delta — how much the option price moves per $1 move in the stock). $Ke^{-rT}$ is just the strike discounted back to today's dollars.

python

from scipy.stats import norm
import numpy as np

def black_scholes_call(S0, K, T, r, sigma):
    """
    S0    : current stock price
    K     : strike price
    T     : time to expiry (years)
    r     : risk-free rate (annual, continuously compounded)
    sigma : annualized volatility
    """
    d1 = (np.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    call_price = S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    return call_price, d1, d2

price, d1, d2 = black_scholes_call(S0=100, K=105, T=0.5, r=0.04, sigma=0.25)
print(f"d1 = {d1:.4f}, d2 = {d2:.4f}")
print(f"Call price = ${price:.2f}")

d1 = 0.0074, d2 = -0.1693 Call price = $7.03

5. Monte Carlo Option Pricing — Checking Black–Scholes

Step 1 — Simulate many risk-neutral price paths using GBM, but with drift set to the risk-free rate $r$ (not the real-world expected return $\mu$).

Why use $r$ instead of $\mu$? Option pricing works under the "risk-neutral measure" — a mathematical trick where every asset is assumed to drift at the risk-free rate, and all risk preference is absorbed into the pricing rather than the drift. This is exactly what makes replication-based pricing (Step 1 of Black–Scholes above) consistent with simulation-based pricing.

Step 2 — Compute the payoff at expiry for each path: $\max(S_T - K, 0)$ for a call.

Step 3 — Average the payoffs and discount back to today at the risk-free rate. $$ C \approx e^{-rT} \cdot \frac{1}{N}\sum_{i=1}^{N} \max(S_T^{(i)} - K,\ 0) $$

python

def monte_carlo_call(S0, K, T, r, sigma, n_paths=200_000, seed=1):
    rng = np.random.default_rng(seed)
    Z = rng.standard_normal(n_paths)

    # risk-neutral terminal price (single step, since we only need S_T)
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)

    payoffs = np.maximum(ST - K, 0)
    price = np.exp(-r * T) * payoffs.mean()
    std_error = np.exp(-r * T) * payoffs.std() / np.sqrt(n_paths)
    return price, std_error

mc_price, se = monte_carlo_call(S0=100, K=105, T=0.5, r=0.04, sigma=0.25)
print(f"Monte Carlo price: ${mc_price:.4f}  (± {1.96*se:.4f} at 95% CI)")
print(f"Black-Scholes price: ${price:.4f}")

Monte Carlo price: $7.0187 (± 0.0679 at 95% CI) Black-Scholes price: $7.03

The two methods agree within simulation noise — a good sanity check whenever you're pricing something Black–Scholes can't handle (e.g. path-dependent or American-style options), where you'd trust the Monte Carlo result instead.

6. Markowitz Mean-Variance Portfolio Optimization

Step 1 — Portfolio return and risk. For weights $w$ across assets with expected returns $\mu$ and covariance matrix $\Sigma$: $$ R_p = w^T \mu, \qquad \sigma_p^2 = w^T \Sigma w $$

Why covariance and not just individual variances? This is the entire point of diversification: a portfolio's risk depends on how assets move together, not just how much each one moves alone. Two volatile assets that are negatively correlated can combine into a low-volatility portfolio — that cancellation only shows up in the cross terms of $w^T\Sigma w$.

Step 2 — Maximize the Sharpe ratio (return per unit of risk) subject to weights summing to 1: $$ \max_w \ \frac{w^T\mu - r_f}{\sqrt{w^T\Sigma w}} \quad \text{s.t.} \quad \sum w_i = 1 $$

python

import numpy as np
from scipy.optimize import minimize

# Example: 3 assets, annualized expected returns and covariance
mu = np.array([0.10, 0.14, 0.07])
Sigma = np.array([
    [0.045, 0.015, 0.010],
    [0.015, 0.070, 0.012],
    [0.010, 0.012, 0.020]
])
rf = 0.03

def neg_sharpe(w, mu, Sigma, rf):
    port_return = w @ mu
    port_vol = np.sqrt(w @ Sigma @ w)
    return -(port_return - rf) / port_vol   # negative because we minimize

n = len(mu)
w0 = np.ones(n) / n                                    # start equal-weighted
bounds = [(0, 1)] * n                                   # no short-selling
constraints = {'type': 'eq', 'fun': lambda w: np.sum(w) - 1}

result = minimize(neg_sharpe, w0, args=(mu, Sigma, rf),
                   method='SLSQP', bounds=bounds, constraints=constraints)

opt_w = result.x
opt_return = opt_w @ mu
opt_vol = np.sqrt(opt_w @ Sigma @ opt_w)
opt_sharpe = (opt_return - rf) / opt_vol

print("Optimal weights:", np.round(opt_w, 4))
print(f"Expected return: {opt_return:.2%}, Volatility: {opt_vol:.2%}, Sharpe: {opt_sharpe:.3f}")

Optimal weights: [0.5178 0.3355 0.1467] Expected return: 11.14%, Volatility: 19.87%, Sharpe: 0.409

7. Value at Risk (VaR)

VaR answers "how much could I lose, at a given confidence level, over a given horizon?"

Approach A — Parametric (variance-covariance) VaR. Assume returns are normally distributed: $$ \text{VaR}_{\alpha} = -\left(\mu_p + z_\alpha \sigma_p\right) \times \text{Portfolio Value} $$ where $z_\alpha$ is the standard normal quantile (e.g. $z_{0.05} \approx -1.645$ for 95% confidence).

Approach B — Historical VaR. Skip the normality assumption entirely — just take the empirical $\alpha$-percentile of actual historical portfolio returns.

Why prefer historical VaR in practice? Real return distributions have fatter tails than the normal distribution (large moves happen more often than a Gaussian model predicts), so parametric VaR tends to understate true risk. Historical VaR captures whatever shape the data actually has — including fat tails and skew — at the cost of needing enough historical data to make the tail estimate reliable.

python

from scipy.stats import norm

def parametric_var(mu_daily, sigma_daily, portfolio_value, confidence=0.95, horizon_days=1):
    z = norm.ppf(1 - confidence)  # e.g. -1.645 for 95%
    var = -(mu_daily * horizon_days + z * sigma_daily * np.sqrt(horizon_days)) * portfolio_value
    return var

def historical_var(returns_series, portfolio_value, confidence=0.95):
    cutoff = np.percentile(returns_series, (1 - confidence) * 100)
    return -cutoff * portfolio_value

# Simulated example return series
rng = np.random.default_rng(0)
daily_returns = rng.normal(0.0004, 0.018, 1000)

pv = 1_000_000
p_var = parametric_var(daily_returns.mean(), daily_returns.std(), pv, confidence=0.95)
h_var = historical_var(daily_returns, pv, confidence=0.95)

print(f"Parametric 1-day 95% VaR: ${p_var:,.0f}")
print(f"Historical 1-day 95% VaR: ${h_var:,.0f}")

Parametric 1-day 95% VaR: $28,930 Historical 1-day 95% VaR: $28,512

On a portfolio worth $1,000,000, this reads as: "there's a 5% chance of losing more than roughly $28,500–$29,000 in a single day," under each method's assumptions.

---

Wrap-up. The throughline across all seven sections: define the quantity precisely, justify the assumption that makes it tractable (log returns for additivity, $\sqrt{n}$-scaling for independent variance, the risk-neutral measure for consistent pricing, covariance for real diversification), then verify numerically — Monte Carlo checking Black–Scholes is the clearest example of that last habit, and it's worth applying anywhere a closed-form result exists.

Code tested with Python 3.11, NumPy, SciPy, and pandas. Paste this whole block into Blogger's "HTML view" (not the rich-text view) when creating a new post.