Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Common Probability Distributions in Finance › Discrete Distributions: Binomial and Poisson

Michael BrenndoerferOctober 21, 202524 min read

Quantitative FinanceAll chapters

Part 1 of 15Mathematical Foundations

Chapter 3Common Probability Distributions in Finance › Discrete Distributions: Binomial and Poisson

Reading section 4 of 8 · Discrete Distributions: Binomial and Poisson · 24 minSection 4 of 8 · 24 min

Section 4 of 8 · 24 min read

Common Probability Distributions in Finance › Discrete Distributions: Binomial and PoissonOne section · Sections in this chapter4 / 8 ↓

Sections in this chapter8

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Not all financial phenomena are continuous. Events like defaults, dividend announcements, or large jumps in prices are discrete occurrences. The binomial and Poisson distributions model these situations and give a framework for counting how many times something happens rather than measuring continuous quantities. These distributions are particularly important in credit risk, where we count defaults in a portfolio, and in derivatives pricing, where discrete lattice models approximate continuous price processes.

The Binomial DistributionLink Copied

The binomial distribution models the number of successes in a fixed number of independent trials, each with the same probability of success. The classic example is flipping a coin multiple times and counting heads, but the same mathematics applies to questions such as: How many stocks in a portfolio will outperform? How many corporate bonds will default? How many days this month will have positive returns?

Binomial Distribution

The probability distribution of the number of successes $k$ in $n$ independent Bernoulli trials, where each trial has success probability $p$ . Used in discrete option pricing models and credit risk modeling.

If $X \sim \text{Binomial}(n, p)$ , the probability mass function is:

P(X = k) = \binom{n}{k} p^k (1-p)^{n-k}, \quad k = 0, 1, \ldots, n

This formula answers the question: if we have $n$ independent trials, each with success probability $p$ , what is the probability of getting exactly $k$ successes? The answer decomposes into two multiplicative components: the number of ways the successes can be arranged, and the probability of any particular arrangement.

Each component of this formula serves a specific purpose:

$k$ : the number of successes we're computing the probability for
$n$ : total number of independent trials
$p$ : probability of success on each trial
$\binom{n}{k} = \frac{n!}{k!(n-k)!}$ : the binomial coefficient, representing the number of ways to choose $k$ successes among $n$ trials
$p^k$ : success-probability factor for $k$ specified success positions
$(1-p)^{n-k}$ : failure-probability factor for $n-k$ specified failure positions

The formula multiplies the number of arrangements by the probability of each specific arrangement. Intuitively, if you flip a coin $n$ times with probability $p$ of heads, the PMF tells you the probability of getting exactly $k$ heads: you need $k$ successes (probability $p^k$ ), $n-k$ failures (probability $(1-p)^{n-k}$ ), and the binomial coefficient counts how many orderings produce that outcome. The key insight is that because the trials are independent, the probability of any specific sequence of outcomes is simply the product of the individual probabilities, and the binomial coefficient tells us how many such sequences lead to exactly $k$ successes.

The mean and variance have elegant formulas that reveal the distribution's structure:

E[X] = np, \quad \text{Var}(X) = np(1-p)

These formulas encode intuitive relationships:

$E[X]$ : the expected number of successes
$\text{Var}(X)$ : the variance in the number of successes
$n$ : total number of trials
$p$ : probability of success on each trial
$(1-p)$ : probability of failure; multiplying by $p$ gives each Bernoulli trial's variance

The expected number of successes, $np$ , is simply the number of trials times the success probability, exactly what intuition suggests. The variance formula $np(1-p)$ is more subtle: it shows that variance is maximized when $p = 0.5$ (complete uncertainty about each trial's outcome) and decreases toward zero as $p$ approaches 0 or 1 (near certainty). This makes sense because if $p$ is very close to 1, almost all trials succeed, leaving little room for variability in the total count.

Out[12]:

Visualization

Bar chart of binomial PMFs for n=20 at p=0.2, 0.5, and 0.8. — Binomial probability mass functions for fixed n=20 and varying success probability p, showing how the mean and skewness shift with p.

Line plot of binomial PMFs for p=0.3 at n=10, 30, and 50. — Binomial probability mass functions for fixed p=0.3 and varying n. On the count scale, the mean and absolute spread increase with n; the success proportion X/n becomes more concentrated around p.

Binomial Option Pricing ModelLink Copied

The Cox-Ross-Rubinstein (CRR) binomial model prices options by assuming that at each time step, the stock price either moves up by a factor $u$ or down by a factor $d$ . This turns a continuous-time model into a tree of possible outcomes that can be analyzed one step at a time. Under shared assumptions for a non-dividend-paying stock and a vanilla European option, its price converges to the Black-Scholes price as the number of time steps increases.

After $n$ steps, the stock price has experienced $k$ up moves and $n-k$ down moves, giving:

S_T = S_0 u^k d^{n-k}

This formula shows how the final price depends entirely on the number of up moves, not their specific timing. The components are:

$S_T$ : stock price at maturity after $n$ steps
$S_0$ : initial stock price
$u$ : up factor (price multiplier for an up move)
$d$ : down factor (price multiplier for a down move)
$k$ : number of up moves
$n-k$ : number of down moves

Whether the stock rises first and then falls, or falls first and then rises, the final price is the same as long as the total number of up and down moves matches. This recombination supports efficient pricing for vanilla European calls and puts, whose payoff depends only on the final price. European exercise alone does not rule out path dependence. Claims such as Asian options require the tree state to be augmented with path information; the simple terminal-price state is not sufficient.

The number of up moves follows a binomial distribution under the risk-neutral measure, with probability $q$ of an up move:

q = \frac{e^{r\Delta t} - d}{u - d}

This risk-neutral probability is derived from the pricing model, not estimated as the physical probability of an up move. It need not equal that physical probability, although the two can coincide. It makes the modeled expected stock growth equal to risk-free growth. The components are:

$q$ : risk-neutral probability of an up move
$r$ : continuously compounded annual risk-free rate
$\Delta t$ : length of each time step in years
$u$ : up factor (stock price multiplier for an up move)
$d$ : down factor (stock price multiplier for a down move)
$e^{r\Delta t}$ : the growth factor for a risk-free investment over one time step

This formula is derived from the no-arbitrage condition: the expected stock return under the risk-neutral measure must equal the risk-free rate. Setting $q \cdot u + (1-q) \cdot d = e^{r\Delta t}$ and solving for $q$ yields this result. In this complete one-period binomial model, a replicable option is priced from the risk-neutral probability rather than the physical up-move probability. This conclusion depends on no arbitrage, market completeness, and the ability to form the replicating portfolio; it is not a universal statement about every derivative market.

Assuming $d<u$ , this $q$ is a strict probability and the one-period stock-and-bond model is arbitrage-free only when

d < e^{r\Delta t} < u,

which is equivalent to $0<q<1$ . Inputs outside this range do not define valid CRR risk-neutral probabilities.

Out[13]:

Visualization

Recombining binomial tree diagram showing stock prices at each node. — A 4-step binomial tree showing possible stock price paths. The tree recombines because an up move followed by a down move gives the same price as a down move followed by an up move. Node labels show the stock price at each node.

In[14]:

Code

import numpy as np
from scipy import stats


def binomial_option_price(S0, K, T, r, sigma, n_steps, option_type="call"):
    """
    Price a vanilla European option on a non-dividend-paying stock using CRR.

    Parameters:
        S0: Initial stock price
        K: Strike price
        T: Time to expiration (years)
        r: Continuously compounded annual risk-free rate
        sigma: Volatility (annualized)
        n_steps: Number of time steps
        option_type: 'call' or 'put'

    Returns:
        Option price

    This instructional implementation uses double precision. Finite parameters
    and the no-arbitrage condition do not guarantee representable terminal
    prices or accurate tiny residuals. The scaled branches avoid selected
    intermediate failures; they are not an arbitrary-range accuracy guarantee.
    """
    if option_type not in {"call", "put"}:
        raise ValueError("option_type must be 'call' or 'put'")
    if not all(np.isfinite(value) for value in (S0, K, T, r, sigma)):
        raise ValueError("Pricing parameters must be finite")
    if S0 <= 0 or K <= 0 or T <= 0 or sigma <= 0:
        raise ValueError("Require S0, K, T and sigma to be positive")
    if (
        isinstance(n_steps, (bool, np.bool_))
        or not isinstance(n_steps, (int, np.integer))
        or n_steps <= 0
    ):
        raise ValueError("n_steps must be a positive integer")

    dt = T / n_steps

    # Keep the ordinary path, with time scaling if T/n rounds to zero.
    log_up = (
        sigma * np.sqrt(dt)
        if dt > 0
        else _scaled_product(sigma, np.sqrt(T), np.sqrt(n_steps))
    )
    rate_step = r * dt if dt > 0 else _scaled_product(r, T, n_steps)
    if not np.isfinite(log_up) or log_up <= 0 or not np.isfinite(rate_step):
        raise FloatingPointError(
            "CRR step coefficients are outside double-precision range"
        )
    if not -log_up < rate_step < log_up:
        raise ValueError("CRR no-arbitrage condition requires d < exp(r*dt) < u")
    if log_up > 350 or min(rate_step + log_up, log_up - rate_step) < 1e-8 * log_up:
        # Algebraically scale by exp(log_up); retain the small probability
        # independently instead of subtracting it from a rounded one.
        denominator = -np.expm1(-2 * log_up)
        log_q = (
            rate_step
            - log_up
            + np.log(-np.expm1(-rate_step - log_up))
            - np.log(denominator)
        )
        q = np.exp(log_q)
        down_probability = -np.expm1(rate_step - log_up) / denominator
        log_down_probability = np.log(-np.expm1(rate_step - log_up)) - np.log(
            denominator
        )
    else:
        q = (np.expm1(rate_step) - np.expm1(-log_up)) / (2 * np.sinh(log_up))
        down_probability = 1 - q
        log_q = np.log(q)
        log_down_probability = np.log1p(-q)
    use_down_probability = (
        q == 1
        or (log_up > 350 or min(rate_step + log_up, log_up - rate_step) < 1e-8 * log_up)
        and down_probability < q
    )

    # Calculate all possible final stock prices
    k = np.arange(n_steps + 1)
    log_moves = (2 * k - n_steps) * log_up
    with np.errstate(over="ignore", under="ignore", invalid="ignore"):
        final_prices = S0 * np.exp(log_moves)
        if not np.all(np.isfinite(final_prices) & (final_prices > 0)):
            final_prices = np.exp(np.log(S0) + log_moves)
    if not np.all(np.isfinite(final_prices) & (final_prices > 0)):
        raise FloatingPointError("CRR terminal prices exceed floating-point range")

    # Calculate payoffs
    if option_type == "call":
        payoffs = np.maximum(final_prices - K, 0)
    else:
        payoffs = np.maximum(K - final_prices, 0)

    # Near the strike, subtracting rounded price levels can erase the payoff.
    relative_spot_gap = (S0 - K) / K if abs(S0 - K) < 0.5 * K else np.inf
    log_spot_strike = (
        np.log1p(relative_spot_gap)
        if abs(relative_spot_gap) < 0.5
        else np.log(S0) - np.log(K)
    )
    log_moneyness = log_spot_strike + log_moves
    near_strike = np.abs(log_moneyness) <= 1e-6
    if option_type == "call":
        payoffs[near_strike] = K * np.maximum(np.expm1(log_moneyness[near_strike]), 0)
    else:
        payoffs[near_strike] = K * np.maximum(-np.expm1(log_moneyness[near_strike]), 0)

    # Risk-neutral probabilities
    # A stable PMF avoids overflowing binomial coefficients in large trees.
    if use_down_probability:
        probs = stats.binom.pmf(n_steps - k, n_steps, down_probability)
    else:
        probs = stats.binom.pmf(k, n_steps, q)

    # Expected payoff under risk-neutral measure, discounted
    with np.errstate(over="ignore", under="ignore", invalid="ignore"):
        option_price = np.exp(-r * T) * np.sum(probs * payoffs)
    positive_payoff = payoffs > 0
    if not np.isfinite(option_price) or np.any((probs == 0) & positive_payoff):
        # The same terminal expectation, with probability/payoff/discount
        # combined before exponentiation if an intermediate loses range.
        from scipy.special import gammaln, logsumexp

        if not np.any(positive_payoff):
            return 0.0
        log_probs = (
            gammaln(n_steps + 1)
            - gammaln(k + 1)
            - gammaln(n_steps - k + 1)
            + k * log_q
            + (n_steps - k) * log_down_probability
        )
        log_expectation = logsumexp(
            log_probs[positive_payoff] + np.log(payoffs[positive_payoff]) - r * T
        )
        with np.errstate(over="ignore", under="ignore"):
            option_price = np.exp(log_expectation)
    if not np.isfinite(option_price):
        raise FloatingPointError("Option price exceeds floating-point representability")

    return option_price


# Price a call option
S0 = 100
K = 100
T = 0.5  # 6 months
r = 0.05
sigma = 0.20

Out[15]:

Console

Binomial Model Call Option Prices
-----------------------------------
Steps:   10  |  Price: $6.7494
Steps:   50  |  Price: $6.8605
Steps:  100  |  Price: $6.8746
Steps:  200  |  Price: $6.8817
Steps:  500  |  Price: $6.8859
-----------------------------------
Black-Scholes:  |  Price: $6.8887

As the number of steps increases, the CRR binomial price converges to the Black-Scholes price under their shared assumptions. The error can oscillate because the strike sits differently relative to successive lattices. First-order convergence, on the scale of $1/n$ , has been established for the European-call setting; this is not a uniform rate guarantee for arbitrary payoffs and parameter choices. The convergence reflects a mathematical connection between the binomial random walk and Brownian motion, formalized in Donsker's theorem.

Out[16]:

Visualization

Line plot showing binomial price converging to Black-Scholes price with increasing steps. — Convergence of the binomial option pricing model to the Black-Scholes benchmark under their shared model assumptions. The oscillation around the continuous-model limit decreases as n grows.

The Poisson DistributionLink Copied

The Poisson distribution models the number of events occurring in a fixed interval when events happen at a constant average rate and independently of each other. Unlike the binomial distribution, which assumes a fixed number of trials, the Poisson distribution allows for any number of events, making it useful for modeling rare occurrences over a continuous interval.

Poisson Distribution

A discrete probability distribution for an event count with expected interval count $\lambda\geq0$ . For a homogeneous Poisson process with intensity $a$ events per unit time, the count over an interval of length $T$ has parameter $\lambda=aT$ . It provides a baseline for rare independent defaults, jump counts, and arrivals.

If $X \sim \text{Poisson}(\lambda)$ , the probability mass function is:

P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \ldots

This formula has a clear structure when we trace its derivation from first principles. Imagine dividing the interval into many tiny subintervals, each with a small probability of containing an event. As the number of subintervals goes to infinity while keeping the expected total number of events fixed at $\lambda$ , we obtain the Poisson distribution.

The components of this formula are:

$k$ : the number of events we're computing the probability for
$\lambda$ : expected number of events in the interval, equal to intensity times interval length
$e^{-\lambda}$ : the probability of zero events, which is a normalization factor
$\lambda^k$ : the intensity factor for $k$ events; together with $e^{-\lambda}$ , it makes $P(X=k)$ largest at $\lambda=k$ for $k>0$
$k!$ : accounts for the indistinguishability of event orderings

The Poisson distribution arises as the limit of a binomial distribution when $n\to\infty$ and $p\to0$ while $np=\lambda$ remains constant. More generally, a Poisson approximation requires independent rare-event opportunities whose largest individual probability vanishes while the sum of their probabilities approaches $\lambda$ . A finite expected count alone is not enough: dependence or a non-negligible individual probability can produce a different count distribution. Credit portfolios with small independent default probabilities and jump models with independent small-interval arrivals can satisfy the approximation. Dependence-driven clustering requires a richer model; occasional groups of arrivals can also occur by chance under a Poisson process.

Both the mean and variance equal $\lambda$ :

E[X] = \text{Var}(X) = \lambda

This property deserves emphasis:

$E[X]$ : the expected number of events
$\text{Var}(X)$ : the variance in the number of events
$\lambda$ : expected interval count, not a per-unit-time rate unless the interval is one time unit

The equality of mean and variance is a distinctive property of the Poisson distribution. When analyzing count data, if the observed variance significantly exceeds the mean (overdispersion), this suggests the Poisson model may be inadequate. Dependence is one possible cause, but rate heterogeneity, omitted covariates, and excess zeros can also produce overdispersion. For instance, defaults may cluster because borrowers share economic conditions even when individual default mechanisms differ.

Out[17]:

Visualization

Line plot of Poisson PMFs for λ equal to 1, 3, 7, and 15. — Poisson probability mass functions for different expected interval counts λ, showing how the distribution becomes more symmetric as λ increases.

Bar chart comparing Poisson(5) probabilities to binomial approximations with n equal to 10, 50, and 200. — Poisson distribution as the rare-event limit of the binomial distribution, comparing Poisson(λ=5) to Binomial(n, λ/n) for increasing n.

Applications of the Poisson DistributionLink Copied

The Poisson distribution is particularly useful for modeling:

Jump processes in asset prices: The Merton jump-diffusion model adds jumps whose arrivals follow a Poisson process and whose sizes follow a separate distribution to geometric Brownian motion. The stock price follows GBM between jumps and moves discontinuously at jump times. This combination represents continuous diffusion-like fluctuations together with occasional discontinuous price moves.
Credit defaults in a portfolio: A large collection of independent, rare default opportunities can admit a Poisson count approximation under the limiting conditions above. Shared economic conditions can make defaults dependent, so a Poisson approximation is not justified by a large portfolio alone.
Order arrivals: A homogeneous Poisson process is a baseline arrival model with independent exponential waiting times. Chance clusters can occur, but the model does not encode self-excitation or time-varying intraday intensity. Those features require richer arrival models.

In[18]:

Code

import numpy as np


def simulate_jump_diffusion(
    S0, mu, sigma, lambda_jump, mu_jump, sigma_jump, T, n_steps, n_paths, rng=None
):
    """
    Simulate Merton jump-diffusion price paths.

    Parameters:
        S0: Initial stock price
        mu: Drift of the diffusion component
        sigma: Volatility of the diffusion component
        lambda_jump: Average number of jumps per year (Poisson intensity)
        mu_jump: Mean of log-jump size
        sigma_jump: Std dev of log-jump size
        T: Time horizon in years
        n_steps: Number of time steps
        n_paths: Number of paths to simulate
        rng: Optional local NumPy RandomState; defaults to a fresh seed-42 stream

    Returns:
        Array of simulated price paths

    Finite inputs do not guarantee representable coefficients, Poisson counts
    or prices. Numerical range failures are not evidence of model invalidity.
    The returned time-zero column directly preserves the floating-point
    representation of S0, not an arbitrary integer's exact precision.
    """
    if not all(
        np.isfinite(value)
        for value in (S0, mu, sigma, lambda_jump, mu_jump, sigma_jump, T)
    ):
        raise ValueError("Simulation parameters must be finite")
    if S0 <= 0 or min(sigma, lambda_jump, sigma_jump, T) < 0:
        raise ValueError("Require S0>0 and nonnegative scales, intensity and horizon")
    if any(
        isinstance(value, (bool, np.bool_))
        or not isinstance(value, (int, np.integer))
        or value <= 0
        for value in (n_steps, n_paths)
    ):
        raise ValueError("n_steps and n_paths must be positive integers")
    dt = T / n_steps
    rng = np.random.RandomState(42) if rng is None else rng

    # Diffusion component
    Z_diffusion = rng.standard_normal((n_paths, n_steps))
    if T == 0:
        diffusion = np.zeros_like(Z_diffusion)
    else:
        ordinary_step = dt > 0 and sigma <= np.sqrt(np.finfo(float).max)
        if ordinary_step:
            with np.errstate(over="ignore", invalid="ignore"):
                diffusion = (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(
                    dt
                ) * Z_diffusion
            ordinary_step = np.all(np.isfinite(diffusion))
        if not ordinary_step:
            volatility_step = _scaled_product(sigma, np.sqrt(T), np.sqrt(n_steps))
            drift_step = _scaled_product(mu, T, n_steps)
            variance_step = _scaled_product(volatility_step, volatility_step, 2)
            diffusion = drift_step - variance_step + volatility_step * Z_diffusion

    # Jump component
    # Number of jumps in each interval (Poisson distributed)
    interval_intensity = (
        lambda_jump * dt
        if dt > 0 or T == 0
        else _scaled_product(lambda_jump, T, n_steps)
    )
    n_jumps = rng.poisson(interval_intensity, (n_paths, n_steps))

    # Log jump sizes (normal)
    jump_sizes = np.zeros((n_paths, n_steps))
    for i in range(n_paths):
        for j in range(n_steps):
            if n_jumps[i, j] > 0:
                jumps = rng.normal(mu_jump, sigma_jump, n_jumps[i, j])
                jump_sizes[i, j] = np.sum(jumps)

    # Combine diffusion and jumps
    log_returns = diffusion + jump_sizes

    # Calculate prices
    log_prices = np.zeros((n_paths, n_steps + 1))
    log_prices[:, 0] = np.log(S0)
    log_prices[:, 1:] = np.log(S0) + np.cumsum(log_returns, axis=1)

    if T == 0:
        return np.full((n_paths, n_steps + 1), S0, dtype=float)
    with np.errstate(over="raise", under="ignore", invalid="raise"):
        prices = np.exp(log_prices)
    prices[:, 0] = S0
    if not np.all(np.isfinite(prices) & (prices > 0)):
        raise FloatingPointError("Price paths exceed floating-point representability")
    return prices


# Compare GBM vs Jump-Diffusion
S0 = 100
mu = 0.10
sigma = 0.20
T = 1.0
n_steps = 252
n_paths = 10000

# Pure GBM paths
gbm_paths = simulate_gbm_paths(S0, mu, sigma, T, n_steps, n_paths)

# Jump-diffusion paths
jd_rng = np.random.RandomState(42)
jd_paths = simulate_jump_diffusion(
    S0,
    mu,
    sigma,
    lambda_jump=2,  # 2 jumps per year on average
    mu_jump=-0.05,  # Jumps tend to be negative
    sigma_jump=0.10,
    T=T,
    n_steps=n_steps,
    n_paths=n_paths,
    rng=jd_rng,
)

gbm_returns = np.diff(np.log(gbm_paths), axis=1).flatten()
jd_returns = np.diff(np.log(jd_paths), axis=1).flatten()

Out[19]:

Console

Return Distribution Comparison
==================================================
Statistic                     GBM  Jump-Diffusion
--------------------------------------------------
Mean                     0.000312       -0.000083
Std Dev                  0.012609        0.016096
Skewness                  -0.0006         -3.1474
Excess Kurtosis           -0.0002         55.5761
Min                       -0.0605         -0.4391
Max                        0.0611          0.3572

With the negatively centered jump-size distribution used here, the simulated jump-diffusion returns have a more negative sample skewness and a larger fourth standardized moment than the simulated GBM returns. Changing the jump-size distribution can change the sign of skewness. Likewise, higher kurtosis alone does not determine a unique allocation of probability among the center, shoulders, and tails; the histogram and quantiles must be examined directly.

Out[20]:

Visualization

Overlaid full-range histograms on a logarithmic density axis, including simulated jump-diffusion daily log-return extrema near -43.9% and 35.7%. — Full-range histograms of simulated daily log returns from geometric Brownian motion and Merton jump-diffusion on a logarithmic density axis. All simulated observations are shown: the GBM sample ranges from about -6.1% to 6.1%, while the jump-diffusion sample ranges from about -43.9% to 35.7%.

Q-Q scatter plot of simulated GBM and jump-diffusion returns versus theoretical normal quantiles. — Q-Q plot against the normal distribution for the simulated GBM and jump-diffusion returns, showing their different departures from normality.