Option prices as probability measurements

economics

Pricing European options from a risk-neutral density, and recovering that density from observed prices.

Published

2026-08-28

Modified

2026-08-28

Why option prices hold a distribution

A European call gives its holder the right to buy a stock at a fixed price on a fixed date. Call the fixed price the strike \(K\), and the stock price on that date \(S_T\). The holder exercises only when it pays. So the payoff is \(\max(S_T-K,\,0)\).

That payoff depends on one thing: where \(S_T\) lands. So the price of the call today is an average of payoffs over the places \(S_T\) might land. Each place contributes its payoff, weighted by what the market pays for a payoff delivered there.

Those weights are a density. Section 4 derives the pricing integral, but the shape of it is already clear:

\[ C(K)=e^{-rT}\int_K^{\infty}(s-K)\,q(s)\,ds . \]

The market does not tell you \(q\). It tells you \(C(K)\), at the strikes it happens to quote.

One price is a number. A chain is a shape.

A single quote is one number, and one number cannot pin down a curve. But nobody quotes a single strike. An expiry comes with a row of them, stacked across the range of plausible outcomes.

Each strike weighs the same unknown density differently. A low strike pays off across almost the whole range. A high strike pays off only at the far right end. The gap between two neighbouring strikes prices the band between them. Quote enough strikes and the prices box in the shape of \(q\).

This is the reason to read a chain as a distribution rather than a price list. The distribution is already in there. It is written in an awkward coordinate system, and the work is getting it out.

One volatility number throws the shape away

Traders usually compress a quote into its implied volatility: the volatility you must feed to Black-Scholes to get the quoted price back. It is a good unit for comparing quotes. It also assumes the answer. Black-Scholes assumes \(S_T\) is lognormal, which fixes the shape and leaves only its width free.

You can see the assumption fail. Feed the model each strike in turn and the implied volatility comes out different at each one. That pattern is the volatility smile, and it is the market saying the true shape is not lognormal. The smile is the distortion. The density is the thing being distorted.

Three features a single number hides:

Skew. Traders pay up for protection against a fall, so weight piles into the left tail. The distribution leans.
Fat tails. The market prices a large move as more likely than a lognormal allows.
Two peaks. Before a court ruling, a merger vote or a drug trial, the market holds two stories at once. A single volatility blurs them into one wide hump. A density shows both.

What you can do once you have \(q\)

Price a payoff nobody quotes. Any European payoff at that expiry is an integral against \(q\). Section 10 does this.
Read a probability off the market. The slope of the call curve is a discounted tail probability, so you can ask what the market charges for the stock finishing below 80, and get a number.
Watch the shape move. Recover \(q\) each day and you can see a tail thicken, or a second peak grow, before an event.

Two warnings

\(q\) is a risk-neutral density, not a forecast. It mixes what people expect with what they will pay to avoid pain. Crash insurance sells for more than the crash is likely, so the left tail is heavier than any honest forecast. Read \(q\) as a price, not a prediction.

And the recovery is unstable. Getting \(q\) from \(C\) takes two derivatives of the price curve, and differentiating twice amplifies noise. Quotes are discrete, wide and stale. Difference them naively and you get negative probabilities. Sections 5 to 9 are about handling that.

This notebook develops both directions of the relationship between terminal-price distributions and European option prices:

Forward problem: choose a risk-neutral probability density \(q(s)\), then compute option prices.
Inverse problem: observe option prices at a finite set of strikes, then infer a density consistent with those prices and no-arbitrage constraints.

The progression is:

probability mass versus probability density;
why a point has probability zero in a continuous model;
expectation as an integral against a density;
numerical and Monte Carlo option pricing;
the Breeden–Litzenberger identity;
the instability of taking two derivatives of noisy quotes;
interpolation and smoothing;
constrained, regularized density recovery;
a non-trivial skewed and multimodal option chain.

All monetary quantities are in generic price units. The methods apply equally to SEK-denominated underlyings.

Core references

Black, F. and Scholes, M. (1973), The Pricing of Options and Corporate Liabilities.
Breeden, D. T. and Litzenberger, R. H. (1978), Prices of State-Contingent Claims Implicit in Option Prices.
Bahra, B. (Bank of England), Implied Risk-Neutral Probability Density Functions from Option Prices.
Malz, A. M. (Federal Reserve Bank of New York), A Simple and Reliable Way to Compute Option-Based Risk-Neutral Distributions.

This notebook is educational rather than trading advice.


import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.integrate import cumulative_trapezoid
from scipy.interpolate import UnivariateSpline
from scipy.optimize import lsq_linear
from scipy.stats import lognorm, norm

warnings.filterwarnings("ignore")

random_generator = np.random.default_rng(20260725)
data_path = Path("synthetic_option_chain.csv")
if not data_path.exists():
    data_path = Path("/mnt/data/synthetic_option_chain.csv")

1. Probability mass and probability density

For a discrete random variable, probabilities can be assigned directly to points:

\[ P(X=x_i)=p_i, \qquad \sum_i p_i=1. \]

For a continuously distributed variable, the primitive object is instead a probability measure over sets. When that measure has a density \(f(x)\) relative to length, probabilities are obtained from areas:

\[ P(a<X<b)=\int_a^b f(x)\,dx. \]

The units clarify the name:

\[ [f(x)] = \frac{\text{probability}}{\text{unit of }x}. \]

Thus \(f(x)\) is a local concentration or density of probability, just as mass density is mass per unit volume.


# Discrete example: probability mass can sit directly on individual outcomes.
die_outcomes = np.arange(1, 7)
die_probabilities = np.full(6, 1.0 / 6.0)

plt.figure(figsize=(8, 4.5))
plt.bar(die_outcomes, die_probabilities)
plt.xlabel("Outcome")
plt.ylabel("Probability mass")
plt.title("Discrete distribution: each point can carry positive probability")
plt.ylim(0.0, 0.22)
plt.show()


# Continuous example: probability is area under a density.
x_grid = np.linspace(-4.0, 4.0, 2001)
normal_density = norm.pdf(x_grid)

interval_lower = -0.8
interval_upper = 0.8
inside_interval = (x_grid >= interval_lower) & (x_grid <= interval_upper)

plt.figure(figsize=(8, 4.5))
plt.plot(x_grid, normal_density)
plt.fill_between(
    x_grid[inside_interval],
    normal_density[inside_interval],
    alpha=0.35,
)
plt.axvline(0.0, linestyle="--")
plt.xlabel("x")
plt.ylabel("Probability per unit x")
plt.title("Continuous distribution: probability is an area")
plt.show()

interval_probability = norm.cdf(interval_upper) - norm.cdf(interval_lower)
print(f"P({interval_lower} < X < {interval_upper}) = {interval_probability:.6f}")
print(f"P(X = 0) = 0 exactly for a continuous normal distribution")

P(-0.8 < X < 0.8) = 0.576289
P(X = 0) = 0 exactly for a continuous normal distribution

Why is the probability of one point exactly zero?

For a continuous density \(f\),

\[ P(X=x_0) = P(X\in\{x_0\}) = \int_{\{x_0\}} f(x)\,dx = 0, \]

because a singleton set has zero length.

This does not mean that \(x_0\) is impossible. In probability theory:

an impossible event has probability zero;
but an event with probability zero need not be impossible.

A continuously distributed draw still produces one exact value. Before the draw, however, every pre-specified point has zero probability.

A useful limiting intuition is a shrinking bin:

\[ P(x_0-\tfrac{h}{2}<X<x_0+\tfrac{h}{2}) \approx f(x_0)h. \]

For every \(h>0\), this probability is positive. As \(h\to0\), it converges to exactly zero. Standard real analysis does not retain a positive “infinitesimal” probability at the limit.


center = 0.0
bin_widths = np.logspace(0, -6, 13)
bin_probabilities = np.array(
    [
        norm.cdf(center + width / 2.0)
        - norm.cdf(center - width / 2.0)
        for width in bin_widths
    ]
)
local_approximations = norm.pdf(center) * bin_widths

shrinking_bins = pd.DataFrame(
    {
        "bin_width": bin_widths,
        "exact_probability": bin_probabilities,
        "density_times_width": local_approximations,
    }
)
display(shrinking_bins)

plt.figure(figsize=(8, 4.5))
plt.loglog(bin_widths, bin_probabilities, marker="o", label="Exact bin probability")
plt.loglog(bin_widths, local_approximations, linestyle="--", label=r"$f(0)h$")
plt.xlabel("Bin width h")
plt.ylabel("Probability")
plt.title("A shrinking interval has probability tending to zero")
plt.legend()
plt.show()

	bin_width	exact_probability	density_times_width
0	1.000000	3.829249e-01	3.989423e-01
1	0.316228	1.256329e-01	1.261566e-01
2	0.100000	3.987761e-02	3.989423e-02
3	0.031623	1.261514e-02	1.261566e-02
4	0.010000	3.989406e-03	3.989423e-03
5	0.003162	1.261566e-03	1.261566e-03
6	0.001000	3.989423e-04	3.989423e-04
7	0.000316	1.261566e-04	1.261566e-04
8	0.000100	3.989423e-05	3.989423e-05
9	0.000032	1.261566e-05	1.261566e-05
10	0.000010	3.989423e-06	3.989423e-06
11	0.000003	1.261566e-06	1.261566e-06
12	0.000001	3.989423e-07	3.989423e-07

2. Terminal stock-price densities

Let \(S_T\) denote the underlying price at maturity \(T\). A terminal density \(q(s)\) satisfies

\[ q(s)\ge 0, \qquad \int_0^\infty q(s)\,ds=1. \]

Under the risk-neutral measure \(Q\), no-arbitrage also implies the discounted mean condition

\[ e^{-rT}\mathbb{E}_Q[S_T] = S_0 e^{-\delta T}, \]

or equivalently,

\[ \mathbb{E}_Q[S_T] = F_{0,T}=S_0e^{(r-\delta)T}. \]

Here:

\(S_0\): spot price today;
\(r\): continuously compounded risk-free rate;
\(\delta\): continuous dividend yield;
\(F_{0,T}\): no-arbitrage forward price;
\(q(s)\): risk-neutral density, not necessarily the real-world forecasting density.


spot = 100.0
rate = 0.03
dividend_yield = 0.00
maturity = 0.75
volatility = 0.24

forward = spot * np.exp((rate - dividend_yield) * maturity)
sigma_t = volatility * np.sqrt(maturity)
log_mean = np.log(forward) - 0.5 * sigma_t**2

terminal_grid = np.linspace(1.0, 220.0, 3000)
risk_neutral_density = lognorm.pdf(
    terminal_grid,
    s=sigma_t,
    scale=np.exp(log_mean),
)

density_area = np.trapz(risk_neutral_density, terminal_grid)
density_mean = np.trapz(
    terminal_grid * risk_neutral_density,
    terminal_grid,
)

plt.figure(figsize=(8, 4.5))
plt.plot(terminal_grid, risk_neutral_density)
plt.axvline(spot, linestyle="--", label="Spot")
plt.axvline(forward, linestyle=":", label="Forward")
plt.xlabel(r"Terminal price $S_T$")
plt.ylabel("Probability density")
plt.title("Lognormal risk-neutral terminal density")
plt.legend()
plt.show()

print(f"Integrated density: {density_area:.8f}")
print(f"Mean terminal price: {density_mean:.4f}")
print(f"No-arbitrage forward: {forward:.4f}")

Integrated density: 0.99992441
Mean terminal price: 102.2580
No-arbitrage forward: 102.2755

3. Forward direction: density \(\rightarrow\) option prices

A European call with strike \(K\) pays

\[ (S_T-K)^+ = \max(S_T-K,0) \]

at maturity. Its no-arbitrage price is the discounted risk-neutral expectation

\[ C(K,T) = e^{-rT}\mathbb{E}_Q[(S_T-K)^+] = e^{-rT}\int_K^\infty (s-K)q(s)\,ds. \]

This is a linear map from the density \(q\) to the call-price curve \(C\). Every strike uses a different “measurement function” \((s-K)^+\).


def price_calls_from_density(
    strike_values,
    state_values,
    density_values,
    risk_free_rate,
    time_to_maturity,
):
    strike_values = np.asarray(strike_values, dtype=float)
    prices = np.empty_like(strike_values)
    discount_factor = np.exp(-risk_free_rate * time_to_maturity)

    for index, strike in enumerate(strike_values):
        payoff = np.maximum(state_values - strike, 0.0)
        prices[index] = discount_factor * np.trapz(
            payoff * density_values,
            state_values,
        )

    return prices


strike_grid = np.linspace(50.0, 160.0, 111)
call_prices_by_integration = price_calls_from_density(
    strike_grid,
    terminal_grid,
    risk_neutral_density,
    rate,
    maturity,
)

plt.figure(figsize=(8, 4.5))
plt.plot(strike_grid, call_prices_by_integration)
plt.xlabel("Strike K")
plt.ylabel("Call price")
plt.title("Forward map: a density produces a decreasing, convex call curve")
plt.show()


def black_scholes_call(
    spot_price,
    strike_price,
    risk_free_rate,
    dividend_rate,
    volatility_value,
    time_to_maturity,
):
    strike_price = np.asarray(strike_price, dtype=float)
    sqrt_time = np.sqrt(time_to_maturity)
    d_1 = (
        np.log(spot_price / strike_price)
        + (
            risk_free_rate
            - dividend_rate
            + 0.5 * volatility_value**2
        )
        * time_to_maturity
    ) / (volatility_value * sqrt_time)
    d_2 = d_1 - volatility_value * sqrt_time

    return (
        spot_price
        * np.exp(-dividend_rate * time_to_maturity)
        * norm.cdf(d_1)
        - strike_price
        * np.exp(-risk_free_rate * time_to_maturity)
        * norm.cdf(d_2)
    )


call_prices_closed_form = black_scholes_call(
    spot,
    strike_grid,
    rate,
    dividend_yield,
    volatility,
    maturity,
)

comparison = pd.DataFrame(
    {
        "strike": strike_grid,
        "density_integration": call_prices_by_integration,
        "black_scholes": call_prices_closed_form,
        "absolute_error": np.abs(
            call_prices_by_integration - call_prices_closed_form
        ),
    }
)
display(comparison.iloc[::20])
print(
    "Maximum absolute discrepancy:",
    comparison["absolute_error"].max(),
)

	strike	density_integration	black_scholes	absolute_error
0	50.0	51.100097	51.113497	0.013400
20	70.0	31.775386	31.787309	0.011923
40	90.0	15.190877	15.201320	0.010443
60	110.0	5.323790	5.332753	0.008963
80	130.0	1.430788	1.438274	0.007486
100	150.0	0.315265	0.321274	0.006009

Maximum absolute discrepancy: 0.01339992560969705

Monte Carlo is another forward method

If \(S_T^{(1)},\ldots,S_T^{(N)}\) are sampled from the terminal density, then

\[ C(K,T) \approx e^{-rT}\frac{1}{N}\sum_{i=1}^N(S_T^{(i)}-K)^+. \]

Numerical integration uses the density explicitly. Monte Carlo uses samples drawn from it. Both compute the same expectation.


number_of_paths = 300_000
normal_draws = random_generator.standard_normal(number_of_paths)
terminal_samples = np.exp(
    log_mean + sigma_t * normal_draws
)

selected_strikes = np.array([80.0, 100.0, 120.0, 145.0])
monte_carlo_prices = []
monte_carlo_standard_errors = []

for strike in selected_strikes:
    discounted_payoffs = np.exp(-rate * maturity) * np.maximum(
        terminal_samples - strike,
        0.0,
    )
    monte_carlo_prices.append(discounted_payoffs.mean())
    monte_carlo_standard_errors.append(
        discounted_payoffs.std(ddof=1) / np.sqrt(number_of_paths)
    )

monte_carlo_table = pd.DataFrame(
    {
        "strike": selected_strikes,
        "monte_carlo": monte_carlo_prices,
        "standard_error": monte_carlo_standard_errors,
        "black_scholes": black_scholes_call(
            spot,
            selected_strikes,
            rate,
            dividend_yield,
            volatility,
            maturity,
        ),
    }
)
display(monte_carlo_table)

	strike	monte_carlo	standard_error	black_scholes
0	80.0	22.824731	0.035584	22.846306
1	100.0	9.327656	0.026354	9.345057
2	120.0	2.843485	0.015183	2.848645
3	145.0	0.468854	0.006079	0.473850

4. Why derivatives of call prices reveal probabilities

Write

\[ C(K)=e^{-rT}\int_K^\infty (s-K)q(s)\,ds. \]

Differentiate with respect to strike. Under standard smoothness conditions,

\[ \frac{\partial C}{\partial K} = -e^{-rT}\int_K^\infty q(s)\,ds = -e^{-rT}Q(S_T>K). \]

Therefore the call-price slope is a discounted tail probability.

Differentiate again:

\[ \boxed{ q(K)=e^{rT}\frac{\partial^2 C}{\partial K^2} } \]

This is the Breeden–Litzenberger identity. Because \(q(K)\ge0\), an arbitrage-free call-price curve must be convex in strike.

A narrow butterfly,

\[ C(K-h)-2C(K)+C(K+h), \]

is the finite-difference version of this second derivative. Its terminal payoff is concentrated near \(K\), so its price measures local probability mass.


butterfly_center = 100.0
butterfly_widths = [20.0, 10.0, 5.0, 2.0]
payoff_grid = np.linspace(60.0, 140.0, 1000)

plt.figure(figsize=(8, 4.5))
for width in butterfly_widths:
    butterfly_payoff = (
        np.maximum(payoff_grid - (butterfly_center - width), 0.0)
        - 2.0 * np.maximum(payoff_grid - butterfly_center, 0.0)
        + np.maximum(payoff_grid - (butterfly_center + width), 0.0)
    )
    normalized_payoff = butterfly_payoff / width**2
    plt.plot(
        payoff_grid,
        normalized_payoff,
        label=f"h = {width:g}",
    )

plt.xlabel(r"Terminal price $S_T$")
plt.ylabel(r"Butterfly payoff divided by $h^2$")
plt.title("Narrow butterflies localize around one strike")
plt.legend()
plt.show()


# Recover the known lognormal density from a very fine, noise-free call curve.
fine_strikes = np.linspace(25.0, 210.0, 2001)
fine_call_prices = black_scholes_call(
    spot,
    fine_strikes,
    rate,
    dividend_yield,
    volatility,
    maturity,
)
first_derivative = np.gradient(fine_call_prices, fine_strikes)
second_derivative = np.gradient(first_derivative, fine_strikes)
recovered_density_fine = np.exp(rate * maturity) * second_derivative

true_density_at_strikes = lognorm.pdf(
    fine_strikes,
    s=sigma_t,
    scale=np.exp(log_mean),
)

plt.figure(figsize=(8, 4.5))
plt.plot(fine_strikes, true_density_at_strikes, label="True density")
plt.plot(
    fine_strikes,
    recovered_density_fine,
    linestyle="--",
    label="Recovered from call curvature",
)
plt.xlim(35.0, 190.0)
plt.xlabel("Terminal price / strike")
plt.ylabel("Density")
plt.title("Breeden–Litzenberger works on a dense, noise-free call curve")
plt.legend()
plt.show()

5. Why the inverse problem is ill-conditioned

Real option data differ from the ideal curve:

strikes are discrete and often irregular;
quotes contain bid–ask spreads, stale prices and rounding;
liquid strikes cover only a finite range;
numerical differentiation amplifies high-frequency noise;
small convexity violations can create negative “densities.”

A second finite difference divides by \(h^2\):

\[ C''(K) \approx \frac{C(K-h)-2C(K)+C(K+h)}{h^2}. \]

Consequently, quote errors of size \(\varepsilon\) can become curvature errors of order \(\varepsilon/h^2\). Finer strike spacing does not automatically solve the problem unless quote noise also falls sufficiently quickly.


# Show explicitly how noise is amplified by a second derivative.
noise_levels = [0.0, 0.005, 0.02]
demonstration_strikes = np.linspace(60.0, 150.0, 181)
clean_demonstration_calls = black_scholes_call(
    spot,
    demonstration_strikes,
    rate,
    dividend_yield,
    volatility,
    maturity,
)

plt.figure(figsize=(8, 4.5))
for noise_level in noise_levels:
    noisy_calls = clean_demonstration_calls + random_generator.normal(
        0.0,
        noise_level,
        size=clean_demonstration_calls.size,
    )
    recovered = np.exp(rate * maturity) * np.gradient(
        np.gradient(noisy_calls, demonstration_strikes),
        demonstration_strikes,
    )
    plt.plot(
        demonstration_strikes,
        recovered,
        label=f"Price noise SD = {noise_level:g}",
    )

plt.xlabel("Strike")
plt.ylabel("Recovered density")
plt.title("Taking two derivatives strongly amplifies quote noise")
plt.ylim(-0.12, 0.16)
plt.legend()
plt.show()

6. Data-driven synthetic option chain

The supplied CSV contains a reproducible synthetic chain generated from a deliberately non-lognormal risk-neutral distribution:

a dominant ordinary regime;
a lower-price crash regime;
a small upside-jump regime;
irregular strikes;
heteroscedastic quote noise and bid–ask spreads.

Because this is synthetic, the true density is known and can be used to evaluate the inverse methods objectively. In real data, it is not known.

Data file: synthetic_option_chain.csv


option_chain = pd.read_csv(data_path)
display(option_chain)

spot = float(option_chain["spot"].iloc[0])
rate = float(option_chain["rate"].iloc[0])
maturity = float(option_chain["maturity_years"].iloc[0])
observed_strikes = option_chain["strike"].to_numpy()
observed_calls = option_chain["mid"].to_numpy()
quote_half_spreads = (
    option_chain["ask"].to_numpy()
    - option_chain["bid"].to_numpy()
) / 2.0

plt.figure(figsize=(8, 4.5))
plt.errorbar(
    observed_strikes,
    observed_calls,
    yerr=quote_half_spreads,
    fmt="o",
    capsize=3,
)
plt.xlabel("Strike")
plt.ylabel("Call price")
plt.title("Observed synthetic calls with bid–ask ranges")
plt.show()

	strike	bid	ask	mid	clean_model_price	spot	rate	maturity_years
0	55.0	46.030757	46.469083	46.249920	46.311194	100.0	0.03	0.75
1	60.0	41.320031	41.724696	41.522363	41.530257	100.0	0.03	0.75
2	65.0	36.542127	36.912767	36.727447	36.824670	100.0	0.03	0.75
3	70.0	31.992686	32.328916	32.160801	32.219914	100.0	0.03	0.75
4	75.0	27.617876	27.919305	27.768591	27.750637	100.0	0.03	0.75
5	80.0	23.328868	23.595135	23.462001	23.469073	100.0	0.03	0.75
6	85.0	19.354430	19.585269	19.469850	19.447752	100.0	0.03	0.75
7	90.0	15.648000	15.843310	15.745655	15.770872	100.0	0.03	0.75
8	95.0	12.429605	12.589513	12.509559	12.516198	100.0	0.03	0.75
9	100.0	9.650259	9.775145	9.712702	9.735781	100.0	0.03	0.75
10	102.5	8.464105	8.596697	8.530401	8.529513	100.0	0.03	0.75
11	105.0	7.380640	7.521121	7.450881	7.444148	100.0	0.03	0.75
12	107.5	6.440566	6.589140	6.514853	6.475572	100.0	0.03	0.75
13	110.0	5.546526	5.703411	5.624969	5.617821	100.0	0.03	0.75
14	115.0	4.145643	4.319853	4.232748	4.204255	100.0	0.03	0.75
15	120.0	3.044722	3.237217	3.140969	3.135153	100.0	0.03	0.75
16	125.0	2.229479	2.441185	2.335332	2.339179	100.0	0.03	0.75
17	130.0	1.653085	1.884844	1.768965	1.751115	100.0	0.03	0.75
18	140.0	0.874703	1.148628	1.011666	0.993756	100.0	0.03	0.75
19	150.0	0.376619	0.694696	0.535657	0.567329	100.0	0.03	0.75
20	165.0	0.000000	0.354722	0.161413	0.234320	100.0	0.03	0.75
21	180.0	0.000000	0.230539	0.002057	0.084173	100.0	0.03	0.75


# Recreate the hidden true density solely for validation plots.
validation_state_grid = np.linspace(20.0, 220.0, 4001)

def lognormal_density_from_forward(
    state_values,
    component_forward,
    volatility_value,
    time_to_maturity,
):
    sigma_time = volatility_value * np.sqrt(time_to_maturity)
    component_log_mean = (
        np.log(component_forward) - 0.5 * sigma_time**2
    )
    return lognorm.pdf(
        state_values,
        s=sigma_time,
        scale=np.exp(component_log_mean),
    )


mixture_weights = np.array([0.78, 0.17, 0.05])
original_component_forwards = np.array([104.5, 78.0, 151.0])
component_volatilities = np.array([0.19, 0.28, 0.20])

temporary_density = sum(
    weight
    * lognormal_density_from_forward(
        validation_state_grid,
        component_forward,
        component_volatility,
        maturity,
    )
    for weight, component_forward, component_volatility in zip(
        mixture_weights,
        original_component_forwards,
        component_volatilities,
    )
)
temporary_density /= np.trapz(
    temporary_density,
    validation_state_grid,
)
temporary_mean = np.trapz(
    validation_state_grid * temporary_density,
    validation_state_grid,
)
target_forward = spot * np.exp(rate * maturity)
forward_scale = target_forward / temporary_mean
component_forwards = original_component_forwards * forward_scale

true_mixture_density = sum(
    weight
    * lognormal_density_from_forward(
        validation_state_grid,
        component_forward,
        component_volatility,
        maturity,
    )
    for weight, component_forward, component_volatility in zip(
        mixture_weights,
        component_forwards,
        component_volatilities,
    )
)
true_mixture_density /= np.trapz(
    true_mixture_density,
    validation_state_grid,
)

plt.figure(figsize=(8, 4.5))
plt.plot(validation_state_grid, true_mixture_density)
plt.xlabel(r"Terminal price $S_T$")
plt.ylabel("Density")
plt.title("Hidden validation density: skewed with secondary regimes")
plt.show()

7. Inverse method A: smooth the call curve, then differentiate

One direct implementation of Breeden–Litzenberger is:

fit a smooth function \(\widehat C(K)\) to observed call prices;
compute \(\widehat C''(K)\);
set \(\widehat q(K)=e^{rT}\widehat C''(K)\).

The smoothing parameter controls the bias–variance trade-off:

too little smoothing follows quote noise and creates oscillatory or negative density;
too much smoothing erases genuine skew or multimodality.

An unconstrained spline does not automatically enforce monotonicity, convexity, non-negativity, normalization or the correct forward mean.


evaluation_strikes = np.linspace(
    observed_strikes.min(),
    observed_strikes.max(),
    700,
)

smoothing_settings = {
    "nearly interpolating": 0.0005,
    "moderate smoothing": 0.10,
    "heavy smoothing": 1.00,
}

spline_density_results = {}

plt.figure(figsize=(8, 4.5))
for label, smoothing_scale in smoothing_settings.items():
    spline = UnivariateSpline(
        observed_strikes,
        observed_calls,
        w=1.0 / np.maximum(quote_half_spreads, 0.01),
        s=smoothing_scale * len(observed_strikes),
        k=3,
    )
    density_estimate = (
        np.exp(rate * maturity)
        * spline.derivative(n=2)(evaluation_strikes)
    )
    spline_density_results[label] = density_estimate
    plt.plot(evaluation_strikes, density_estimate, label=label)

plt.plot(
    validation_state_grid,
    true_mixture_density,
    linestyle="--",
    linewidth=2.0,
    label="True density",
)
plt.xlim(observed_strikes.min(), observed_strikes.max())
plt.ylim(-0.03, 0.065)
plt.xlabel("Terminal price / strike")
plt.ylabel("Density estimate")
plt.title("Spline differentiation is sensitive to smoothing")
plt.legend()
plt.show()

8. Inverse method B: solve directly for non-negative probability masses

Discretize terminal prices into grid points \(s_j\) and represent the unknown distribution by probability masses \(p_j\):

\[ p_j\ge0, \qquad \sum_j p_j=1. \]

For observed strikes \(K_i\),

\[ C_i \approx e^{-rT}\sum_j (s_j-K_i)^+p_j. \]

In matrix form,

\[ \mathbf c \approx A\mathbf p. \]

This inverse problem is underdetermined or ill-conditioned because many fine-scale distributions reprice a sparse option chain similarly. Add a smoothness penalty:

\[ \min_{\mathbf p\ge0} \left\| W(A\mathbf p-\mathbf c) \right\|_2^2 + \lambda \left\| D_2\mathbf p \right\|_2^2, \]

subject approximately to

\[ \sum_jp_j=1, \qquad \sum_js_jp_j=F_{0,T}. \]

Here \(D_2\) is a second-difference matrix. The parameter \(\lambda\) is the regularization strength.

This approach uses finite-dimensional linear algebra, but it is a discretization of an inverse problem between function spaces.


def build_second_difference_matrix(number_of_points):
    matrix = np.zeros((number_of_points - 2, number_of_points))
    row_indices = np.arange(number_of_points - 2)
    matrix[row_indices, row_indices] = 1.0
    matrix[row_indices, row_indices + 1] = -2.0
    matrix[row_indices, row_indices + 2] = 1.0
    return matrix


def recover_probability_masses(
    strikes,
    call_prices,
    state_values,
    risk_free_rate,
    time_to_maturity,
    quote_scales,
    regularization_strength,
    constraint_strength=3_000.0,
):
    discount_factor = np.exp(
        -risk_free_rate * time_to_maturity
    )
    payoff_matrix = discount_factor * np.maximum(
        state_values[None, :] - strikes[:, None],
        0.0,
    )

    weights = 1.0 / np.maximum(quote_scales, 0.02)
    weighted_payoff_matrix = weights[:, None] * payoff_matrix
    weighted_prices = weights * call_prices

    second_difference = build_second_difference_matrix(
        len(state_values)
    )

    target_forward = spot * np.exp(
        risk_free_rate * time_to_maturity
    )

    normalization_row = np.ones((1, len(state_values)))
    forward_row = state_values.reshape(1, -1) / target_forward

    augmented_matrix = np.vstack(
        [
            weighted_payoff_matrix,
            np.sqrt(regularization_strength) * second_difference,
            constraint_strength * normalization_row,
            constraint_strength * forward_row,
        ]
    )
    augmented_target = np.concatenate(
        [
            weighted_prices,
            np.zeros(second_difference.shape[0]),
            np.array([constraint_strength]),
            np.array([constraint_strength]),
        ]
    )

    result = lsq_linear(
        augmented_matrix,
        augmented_target,
        bounds=(0.0, np.inf),
        method="trf",
        lsmr_tol="auto",
        max_iter=2_000,
    )

    probability_masses = result.x
    probability_masses /= probability_masses.sum()

    fitted_prices = payoff_matrix @ probability_masses
    return probability_masses, fitted_prices, result


inverse_state_grid = np.linspace(20.0, 220.0, 401)
inverse_grid_width = inverse_state_grid[1] - inverse_state_grid[0]

regularization_strength = 4_000.0
recovered_masses, fitted_calls, optimization_result = (
    recover_probability_masses(
        observed_strikes,
        observed_calls,
        inverse_state_grid,
        rate,
        maturity,
        quote_half_spreads,
        regularization_strength,
    )
)
recovered_density = recovered_masses / inverse_grid_width

print("Optimization success:", optimization_result.success)
print("Recovered total probability:", recovered_masses.sum())
print(
    "Recovered mean terminal price:",
    np.sum(inverse_state_grid * recovered_masses),
)
print("Target forward:", target_forward)

Optimization success: True
Recovered total probability: 1.0
Recovered mean terminal price: 102.27541025190442
Target forward: 102.27550341644461


plt.figure(figsize=(8, 4.5))
plt.plot(
    inverse_state_grid,
    recovered_density,
    label="Regularized inverse estimate",
)
plt.plot(
    validation_state_grid,
    true_mixture_density,
    linestyle="--",
    label="True density",
)
plt.xlabel(r"Terminal price $S_T$")
plt.ylabel("Density")
plt.title("Constrained regularized recovery")
plt.legend()
plt.show()


plt.figure(figsize=(8, 4.5))
plt.errorbar(
    observed_strikes,
    observed_calls,
    yerr=quote_half_spreads,
    fmt="o",
    capsize=3,
    label="Observed quotes",
)
plt.plot(
    observed_strikes,
    fitted_calls,
    marker="x",
    label="Repriced from recovered distribution",
)
plt.xlabel("Strike")
plt.ylabel("Call price")
plt.title("Recovered distribution reprices the observed chain")
plt.legend()
plt.show()

repricing_table = pd.DataFrame(
    {
        "strike": observed_strikes,
        "observed_mid": observed_calls,
        "recovered_price": fitted_calls,
        "residual": fitted_calls - observed_calls,
        "half_spread": quote_half_spreads,
    }
)
display(repricing_table)

	strike	observed_mid	recovered_price	residual	half_spread
0	55.0	46.249920	46.266585	0.016665	0.219163
1	60.0	41.522363	41.498525	-0.023839	0.202333
2	65.0	36.727447	36.738578	0.011131	0.185320
3	70.0	32.160801	32.160320	-0.000481	0.168115
4	75.0	27.768591	27.767086	-0.001505	0.150715
5	80.0	23.462001	23.463768	0.001767	0.133134
6	85.0	19.469850	19.468366	-0.001484	0.115420
7	90.0	15.745655	15.746840	0.001184	0.097655
8	95.0	12.509559	12.508506	-0.001054	0.079954
9	100.0	9.712702	9.714621	0.001919	0.062443
10	102.5	8.530401	8.525796	-0.004605	0.066296
11	105.0	7.450881	7.457765	0.006885	0.070241
12	107.5	6.514853	6.507800	-0.007053	0.074287
13	110.0	5.624969	5.628997	0.004028	0.078442
14	115.0	4.232748	4.231412	-0.001336	0.087105
15	120.0	3.140969	3.141701	0.000732	0.096248
16	125.0	2.335332	2.334964	-0.000368	0.105853
17	130.0	1.768965	1.769086	0.000122	0.115880
18	140.0	1.011666	1.011644	-0.000022	0.136962
19	150.0	0.535657	0.535667	0.000010	0.159039
20	165.0	0.161413	0.161406	-0.000007	0.177361
21	180.0	0.002057	0.002075	0.000018	0.115269

9. Choosing the regularization strength

There is no universally correct \(\lambda\).

A small value prioritizes quote fit and permits rough densities. A large value prioritizes smoothness and may suppress genuine distributional features.

Below, the weighted root-mean-square pricing error is compared with a roughness measure. This is an elementary version of an L-curve diagnostic. In real applications one might use cross-validation, generalized cross-validation, Bayesian priors, information criteria, or economic constraints.


regularization_grid = np.logspace(0, 6, 19)
diagnostics = []

for candidate_strength in regularization_grid:
    candidate_masses, candidate_prices, candidate_result = (
        recover_probability_masses(
            observed_strikes,
            observed_calls,
            inverse_state_grid,
            rate,
            maturity,
            quote_half_spreads,
            candidate_strength,
        )
    )
    weighted_residuals = (
        candidate_prices - observed_calls
    ) / np.maximum(quote_half_spreads, 0.02)
    weighted_rmse = np.sqrt(np.mean(weighted_residuals**2))
    roughness = np.linalg.norm(
        np.diff(candidate_masses, n=2)
    )
    diagnostics.append(
        {
            "regularization_strength": candidate_strength,
            "weighted_rmse": weighted_rmse,
            "roughness": roughness,
        }
    )

diagnostic_table = pd.DataFrame(diagnostics)
display(diagnostic_table)

plt.figure(figsize=(8, 4.5))
plt.loglog(
    diagnostic_table["roughness"],
    diagnostic_table["weighted_rmse"],
    marker="o",
)
for _, row in diagnostic_table.iloc[::4].iterrows():
    plt.annotate(
        f"{row['regularization_strength']:.0e}",
        (row["roughness"], row["weighted_rmse"]),
    )
plt.xlabel("Distribution roughness")
plt.ylabel("Weighted pricing RMSE")
plt.title("L-curve-style regularization diagnostic")
plt.show()

	regularization_strength	weighted_rmse	roughness
0	1.000000	0.028947	0.008106
1	2.154435	0.028947	0.006730
2	4.641589	0.028947	0.006725
3	10.000000	0.028948	0.006716
4	21.544347	0.028951	0.006695
5	46.415888	0.029069	0.006406
6	100.000000	0.029318	0.006012
7	215.443469	0.029723	0.005729
8	464.158883	0.031002	0.005260
9	1000.000000	0.034168	0.004595
10	2154.434690	0.040683	0.003741
11	4641.588834	0.050322	0.002812
12	10000.000000	0.060428	0.002057
13	21544.346900	0.069967	0.001527
14	46415.888336	0.080101	0.001126
15	100000.000000	0.091498	0.000793
16	215443.469003	0.103282	0.000532
17	464158.883361	0.114576	0.000325
18	1000000.000000	0.122802	0.000199

10. Non-trivial use: infer tail probabilities and price unseen payoffs

Once a terminal distribution has been recovered, it can be used to estimate any payoff \(g(S_T)\) that depends only on terminal price:

\[ V_g=e^{-rT}\mathbb{E}_Q[g(S_T)]. \]

We will evaluate:

a cash-or-nothing digital paying 1 if \(S_T<75\);
a corridor digital paying 1 if \(90<S_T<110\);
a capped quadratic payoff \[ \min\left[\left(\frac{S_T-100}{20}\right)^2,4\right]. \]

These payoffs were not included in the calibration. This tests whether the recovered distribution captures economically relevant shape rather than merely reproducing listed calls.

The results remain risk-neutral prices/probabilities. They are not direct real-world event forecasts.


def discounted_expectation_from_masses(
    payoff_values,
    probability_masses,
    risk_free_rate,
    time_to_maturity,
):
    return np.exp(
        -risk_free_rate * time_to_maturity
    ) * np.sum(payoff_values * probability_masses)


true_masses_on_inverse_grid = np.interp(
    inverse_state_grid,
    validation_state_grid,
    true_mixture_density,
)
true_masses_on_inverse_grid = (
    true_masses_on_inverse_grid * inverse_grid_width
)
true_masses_on_inverse_grid /= true_masses_on_inverse_grid.sum()

payoffs = {
    "downside digital, S_T < 75": (
        inverse_state_grid < 75.0
    ).astype(float),
    "corridor digital, 90 < S_T < 110": (
        (inverse_state_grid > 90.0)
        & (inverse_state_grid < 110.0)
    ).astype(float),
    "capped quadratic payoff": np.minimum(
        ((inverse_state_grid - 100.0) / 20.0) ** 2,
        4.0,
    ),
}

rows = []
for payoff_name, payoff_values in payoffs.items():
    recovered_value = discounted_expectation_from_masses(
        payoff_values,
        recovered_masses,
        rate,
        maturity,
    )
    true_value = discounted_expectation_from_masses(
        payoff_values,
        true_masses_on_inverse_grid,
        rate,
        maturity,
    )
    rows.append(
        {
            "payoff": payoff_name,
            "recovered_price": recovered_value,
            "true_price": true_value,
            "absolute_error": abs(recovered_value - true_value),
        }
    )

unseen_payoff_table = pd.DataFrame(rows)
display(unseen_payoff_table)

	payoff	recovered_price	true_price	absolute_error
0	downside digital, S_T < 75	0.101124	0.098722	0.002401
1	corridor digital, 90 < S_T < 110	0.354040	0.365139	0.011099
2	capped quadratic payoff	1.017825	1.017078	0.000747


# Compare the recovered and true cumulative distributions.
recovered_cdf = np.cumsum(recovered_masses)
true_cdf = np.cumsum(true_masses_on_inverse_grid)

plt.figure(figsize=(8, 4.5))
plt.plot(
    inverse_state_grid,
    true_cdf,
    label="True CDF",
)
plt.plot(
    inverse_state_grid,
    recovered_cdf,
    linestyle="--",
    label="Recovered CDF",
)
plt.xlabel(r"Terminal threshold $s$")
plt.ylabel(r"$Q(S_T \leq s)$")
plt.title("CDF comparison emphasizes tail probabilities")
plt.legend()
plt.show()

11. What is and is not identified by option prices?

With calls at every strike and ideal noiseless conditions, the full call curve identifies the risk-neutral terminal distribution.

With finitely many strikes:

only a finite number of linear measurements is observed;
many distributions can fit within bid–ask spreads;
behavior below the lowest and above the highest strike is weakly identified;
regularization and tail assumptions influence the answer;
American exercise, dividends, rates and settlement conventions must be handled correctly;
one maturity identifies only one marginal distribution, not the full stochastic process through time.

This last point is important. Two stochastic models can have the same terminal distribution at one maturity and therefore the same European option prices at that maturity, while having very different path behavior and barrier-option prices.

Summary

The forward and inverse views are:

\[ \boxed{ q(s) \;\xrightarrow{\text{discounted payoff integrals}}\; C(K) } \]

and

\[ \boxed{ C(K_i)+\text{noise} \;\xrightarrow{\text{smoothing or regularized inversion}}\; \widehat q(s). } \]

Probability density is a local rate of probability accumulation. Call options are linear measurements of that density. Curvature of the call-price curve reveals local state-price density, but extracting curvature from sparse noisy quotes is unstable, which is why numerical regularization and no-arbitrage constraints are central.

Exercises

Replace the mixture density with a Student-\(t\)-like distribution on log returns and repeat the forward calculation.
Increase quote noise and inspect how the optimal regularization region changes.
Remove all strikes below 85 and quantify the deterioration in the recovered downside tail.
Add explicit inequality constraints enforcing decreasing and convex call prices before spline differentiation.
Recover a separate distribution for several maturities and examine how difficult it is to enforce calendar consistency.
Compare a density penalty, a log-density penalty and an entropy penalty.
Use put-call parity to combine liquid out-of-the-money puts and calls into one synthetic call curve.