exotics

A small pricing library for three exotic option types: barrier options, Asian options, and cliquets. Each pricer is a single function in exotics/exotics.py. This document works through the mathematics behind each one: the model, the derivation, and how that derivation maps onto the code.

All three pricers share the same underlying model. The underlying asset follows geometric Brownian motion under the risk neutral measure:

dS_t = (r - q) S_t dt + sigma S_t dW_t

where r is the continuously compounded risk free rate, q is the continuous dividend yield, and sigma is the (constant, flat) volatility. Every closed form and every simulation in this library starts from this dynamic. Nothing here models stochastic volatility, jumps, or a term structure of rates or vol; those are stated limitations, not oversights.

Notation

Throughout:

• S0 is the spot price today.
• K is the strike.
• T is time to maturity in years.
• r, q, sigma as above.
• N(x) is the standard normal CDF, norm.cdf in the code.
• tau denotes time to maturity measured from some point, i.e. tau = T - t, distinct from calendar time t.

Barrier options: price_barrier_option

The pricing PDE

A European contingent claim V(S, t) on this underlying, in the absence of a barrier, satisfies the Black Scholes PDE:

dV/dt + 0.5 * sigma^2 * S^2 * d2V/dS2 + (r - q) * S * dV/dS - r * V = 0

with terminal condition V(S, T) = payoff(S). A barrier feature is a path dependent constraint layered on top: at each monitoring date, if the underlying has crossed the barrier H, the contract is either extinguished (knock out, paying a rebate) or activated (knock in). There is no closed form for this in general once early exercise or discrete monitoring is involved, so the code solves the PDE numerically.

Change of variable

It's more convenient to march forward in time to maturity than backward in calendar time. Define tau = T - t, so tau = 0 is maturity and tau = T is today. Substituting into the PDE flips the sign of the time derivative:

dV/dtau = 0.5 * sigma^2 * S^2 * d2V/dS2 + (r - q) * S * dV/dS - r * V

Now the problem is an initial value problem in tau: start from the terminal payoff at tau = 0 and step forward to tau = T, which corresponds to today. This is exactly what the step loop in _run_grid does: tau_old = stop_times[step-1], tau_new = stop_times[step], and the loop runs until tau reaches T.

Spatial discretization

Lay down a grid of asset prices S_0 = 0 < S_1 < ... < S_n = S_max, evenly spaced with spacing dS, where S_max is S_max_mult * K (three times the strike by default). Index the interior points i = 1, ..., n-1, so S_i = i * dS. Approximate the spatial derivatives at each interior node with central differences:

dV/dS   ~ (V_{i+1} - V_{i-1}) / (2 dS)
d2V/dS2 ~ (V_{i+1} - 2 V_i + V_{i-1}) / dS^2

Substituting these into the PDE and collecting terms in V_{i-1}, V_i, V_{i+1} gives, after simplification (using S_i = i * dS so the dS cancels out of every coefficient):

a_i = 0.25 * dt * (sigma^2 * i^2 - (r - q) * i)
b_i = -0.5 * dt * (sigma^2 * i^2 + r)
c_i = 0.25 * dt * (sigma^2 * i^2 + (r - q) * i)

These are exactly alpha, beta, gamma in _cn_coeffs.

Time discretization: Crank Nicolson

A fully explicit scheme (evaluate the spatial derivatives at the old time level) is only conditionally stable: it blows up unless dt is kept very small relative to dS^2. A fully implicit scheme is unconditionally stable but only first order accurate in time. Crank Nicolson averages the explicit and implicit operators, which is unconditionally stable and second order accurate in both S and tau. Writing the discretized PDE as an operator equation dV/dtau = L V, Crank Nicolson is:

(V^{n+1} - V^n) / dt = 0.5 * L V^{n+1} + 0.5 * L V^n

Substituting the tridiagonal form of L and moving the unknown time level to the left gives two tridiagonal systems, one implicit (left hand side) and one explicit (applied to the known vector):

Implicit (M1):  -a_i V_{i-1}^{n+1} + (1 - b_i) V_i^{n+1} - c_i V_{i+1}^{n+1}
Explicit (M2):   a_i V_{i-1}^n     + (1 + b_i) V_i^n     + c_i V_{i+1}^n

and the system solved at every step is M1 V^{n+1} = M2 V^n + boundary terms. This is precisely the sub, diag, sup (implicit, M1) and M2_sub, M2_diag, M2_sup (explicit, M2) arrays built in _cn_coeffs, and the matrix vector product _tridiag_matvec(M2_sub, M2_diag, M2_sup, V[1:-1]) computes the right hand side.

Boundary conditions

At S = 0, the asset is worthless and stays worthless, so a call is worth 0 and a put is worth its discounted strike K * exp(-r * tau). At S = S_max, far above any reasonable strike, a call behaves like a forward contract (the optionality is negligible) so its value is approximated by S_max * exp(-q * tau) - K * exp(-r * tau), and symmetrically a put is worth 0. This is _boundaries.

Because Crank Nicolson evaluates both time levels, the boundary contribution to the first and last interior equations needs the boundary value at both tau_old and tau_new, which is why the code adds alpha[0] * (V0_old + V0_new) and gamma[-1] * (Vmax_old + Vmax_new) to the right hand side rather than just one of them.

Solving the tridiagonal system: the Thomas algorithm

A tridiagonal linear system can be solved in O(n) time rather than the O(n^3) of general Gaussian elimination, by eliminating the sub diagonal in a forward sweep and then back substituting:

Forward sweep, for i = 1 .. n-1:
    w      = sub_i / diag_{i-1}
    diag_i = diag_i - w * sup_{i-1}
    rhs_i  = rhs_i  - w * rhs_{i-1}

Back substitution:
    x_{n-1} = rhs_{n-1} / diag_{n-1}
    x_i     = (rhs_i - sup_i * x_{i+1}) / diag_i,  for i = n-2 .. 0

This is _thomas_solve, used whenever there is no early exercise decision to make.

American exercise as a linear complementarity problem

Early exercise turns the pricing problem into a linear complementarity problem: at every point on the grid, the value must satisfy V >= payoff(S) and the discretized PDE relation with equality wherever V > payoff(S). Two solvers are provided.

Brennan Schwartz solves this directly, still in O(n), by applying the exercise constraint during the back substitution step rather than iterating:

x_{n-1} = max(rhs_{n-1} / diag_{n-1}, intrinsic_{n-1})
x_i     = max((rhs_i - sup_i * x_{i+1}) / diag_i, intrinsic_i)

The key point is that x_{i+1} used in the recursion for x_i is the already projected value, not the raw linear solve. This is what makes it a valid LCP solver rather than a linear solve followed by a naive clamp, and it is the actual 1977 Brennan Schwartz algorithm, not an approximation of it.

PSOR (projected successive over relaxation) is the iterative alternative, used as a slower cross check. At each grid point it computes a Gauss Seidel update and relaxes it with factor omega, then projects onto the exercise constraint:

y      = (rhs_k - sub_k * x_{k-1} - sup_k * x_{k+1}) / diag_k
x_k    = max(intrinsic_k, x_k + omega * (y - x_k))

repeated until the update is smaller than psor_tol or psor_max_iter is reached. Because this is a nested Python loop over grid points and iterations, it is substantially slower than Brennan Schwartz; it exists as an independent check, not as the default. If it does not converge within psor_max_iter, the code now raises a warning naming how many time steps failed to converge, rather than silently returning an unconverged result.

Barrier monitoring

Discrete (end of day) monitoring is implemented as a projection: after each time step is solved, any grid point on or beyond the barrier is set to the rebate value.

up barrier:   V(S) = rebate  for all S >= H
down barrier: V(S) = rebate  for all S <= H

applied once to the terminal payoff (to catch a barrier breach exactly at maturity) and then again after every subsequent step. Because this projects the whole vector V rather than a single node, it correctly handles the fact that H generally does not sit exactly on a grid point.

For this projection to mean anything, H has to be strictly inside the grid, 0 < H < S_max. If H were at or beyond S_max, the condition S_grid >= H would never be true, the projection would never fire, and an up and out option would silently price as if it had no barrier at all. The code checks this explicitly and raises rather than allowing that silent failure. It also warns (without changing the answer) if S0 already sits on the knocked side of H, since that is either an already dead knock out or an already active knock in, and warns separately if H sits very close to either edge of the grid, since interpolation accuracy degrades there.

Bermudan exercise and the exact time grid

For Bermudan exercise, after solving each time step the code additionally enforces V = max(V, payoff(S)) on any step whose tau coincides with an exercise date T - exercise_date.

Rather than rounding each exercise date to the nearest daily grid point, _build_time_grid merges the uniform daily grid with the exact set of tau values implied by exercise_dates, using numpy.union1d. This produces a nonuniform time grid where every exercise date lands exactly on a step boundary, at the cost of one or two irregularly sized steps near each exercise date. Crank Nicolson does not require a uniform dt, so this is a free improvement in accuracy: the coefficients alpha, beta, gamma are simply recomputed from the local dt at every step rather than once for the whole grid.

In-out parity for knock-in options

A European knock-in and the corresponding knock-out with the same barrier are complementary: exactly one of them ends up being live. With no rebate, this gives the identity

knock_in + knock_out = vanilla

which holds exactly, not approximately, when both legs are priced on the same PDE grid, because the difference of two linear finite difference solutions on identical grids inherits the identity from the continuous problem. The code exploits this directly: knock_in = vanilla - knock_out, computed by calling _run_grid twice.

This identity assumes no early exercise decision is folded into the picture. If the holder can exercise early, the exercise decision itself depends on whether the barrier has been hit, which reintroduces the path dependence that the parity argument cancels out. Pricing that correctly needs a second grid dimension tracking activation state, which is out of scope here, so Bermudan and American knock-in combinations raise NotImplementedError rather than silently returning the (wrong) parity based answer.

Rebates

For a knock-out, the rebate is paid at the moment the barrier is breached, which the projection step already handles: the grid point is simply set to rebate at the monitoring date the breach occurs, and that value then earns the risk free rate for the remainder of the CN evolution automatically, exactly as any other value on the grid would.

For a knock-in, the convention is different: the rebate is paid at maturity only if the barrier is never breached, i.e. only if the option never activates. This is a survival contingent payment, valued as a separate mini pricing problem: a claim worth rebate at maturity everywhere, subject to the same knock-out style barrier (so it is worth 0 wherever the real barrier has already been breached). Its boundary condition on the side away from the barrier is the discounted rebate rebate * exp(-r * tau), since deep in that region survival to maturity is effectively certain. This "survival bond" price is added to the parity based knock-in value:

knock_in_total = (vanilla - knock_out) + survival_bond

which is the _bond_payoff / _bond_bounds branch in the code.

Asian options: price_asian_option

An Asian option's payoff depends on the average of the underlying over a set of observation dates rather than its terminal value alone. With N total observation times t_1 < ... < t_N = T, define the arithmetic average A = (1/N) * sum(S_{t_i}) or the geometric average G = (prod(S_{t_i}))^(1/N). A fixed strike call pays max(A - K, 0); a floating strike call pays max(S_T - A, 0).

Geometric average: exact closed form

The geometric average has a genuinely nice property under GBM: it is itself lognormal, because a product of lognormals is lognormal and a geometric mean is a sum of logs divided by a constant. Under the risk neutral measure, ln(S_{t_i}) = ln(S0) + (r - q - 0.5 sigma^2) t_i + sigma W_{t_i}, so

ln(G) = (1/N) * sum_i ln(S_{t_i}) = ln(S0) + (r - q - 0.5 sigma^2) * mean(t_i) + sigma * (1/N) * sum_i W_{t_i}

The first two terms are deterministic; call their sum m. The last term is a sum of correlated Gaussians. Using Cov(W_s, W_t) = min(s, t), its variance works out to

v = Var(ln G) = sigma^2 / N^2 * sum_i sum_j min(t_i, t_j)

So ln(G) ~ Normal(m, v), exactly, with no approximation. That makes G lognormal, and the expected discounted payoff of a call on a lognormal variable has the same shape as the Black Scholes formula, just with S0 replaced by the lognormal's own expectation E[G] = exp(m + v/2):

d1 = (m - ln(K) + v) / sqrt(v)
d2 = d1 - sqrt(v)
price = exp(-r T) * [E[G] * N(d1) - K * N(d2)]     (call)
price = exp(-r T) * [K * N(-d2) - E[G] * N(-d1)]   (put)

This is _geometric_closed_form. It requires v > 0, which fails only in the degenerate case sigma = 0 or T = 0; the code checks for this and raises rather than dividing by zero.

Arithmetic average: Turnbull-Wakeman moment matching

The arithmetic average of correlated lognormals is not itself lognormal, so there is no exact closed form. The standard practical approximation (Turnbull and Wakeman, and originally Levy) is to match the true first and second moments of the arithmetic average to those of a substitute lognormal variable, then price against that substitute using the same Black Scholes shaped formula.

The first moment is easy, since E[S_{t_i}] = S0 * exp((r - q) t_i) and expectation is linear:

M1 = E[A] = (1/N) * sum_i S0 * exp((r - q) t_i)

The second moment needs the covariance structure of correlated lognormals. For i <= j,

E[S_{t_i} S_{t_j}] = S0^2 * exp((r - q)(t_i + t_j)) * exp(sigma^2 * min(t_i, t_j))

so

M2 = E[A^2] = (1/N^2) * sum_i sum_j E[S_{t_i} S_{t_j}]

Now match these to a lognormal variable X with ln(X) ~ Normal(m, v), which has E[X] = exp(m + v/2) and E[X^2] = exp(2m + 2v). Setting E[X] = M1 and E[X^2] = M2 and solving the two equations for m and v gives

v = ln(M2 / M1^2)
m = ln(M1) - v / 2

and pricing proceeds exactly as in the geometric case, with M1 playing the role of E[G]:

d1 = (m - ln(K) + v) / sqrt(v)
d2 = d1 - sqrt(v)
price = exp(-r T) * [M1 * N(d1) - K * N(d2)]     (call)
price = exp(-r T) * [K * N(-d2) - M1 * N(-d1)]   (put)

This is _moment_matching. It is an approximation, but a very close one in practice, which is why it is the default method for arithmetic, fixed strike Asians. As with the geometric case, v <= 0 is degenerate (again sigma = 0 or T = 0) and raises rather than silently returning nan.

Monte Carlo

The general purpose fallback, used for anything the two closed forms above do not cover, in particular floating strike payoffs. Simulate log returns directly:

S_{t_{i+1}} = S_{t_i} * exp((r - q - 0.5 sigma^2) dt + sigma sqrt(dt) Z_i),   Z_i ~ Normal(0, 1)

vectorized as a cumulative sum of increments in log space, take the arithmetic or geometric average of the simulated path per the averaging_type argument, apply the fixed or floating strike payoff, discount, and average over paths. return_stderr=True additionally returns the Monte Carlo standard error, std(discounted payoff) / sqrt(n_paths), since a point estimate with 500,000 paths is not automatically trustworthy without knowing its own noise.

Seasoning: partially realized averages

All of the above assumes the averaging window starts fresh, at S0, today. In practice an Asian option is usually revalued partway through its life, after some fixings have already occurred. If n_realized fixings have already been observed with a known average realized_avg, and n_obs fixings remain between now and maturity T, the total average over all N = n_realized + n_obs fixings decomposes as a weighted combination of the known past and the unknown future:

A_total = (n_realized * realized_avg + n_obs * A_future) / N
        = w_real * realized_avg + w_future * A_future

with w_real = n_realized / N and w_future = n_obs / N. For a fixed strike call this lets the payoff be rewritten entirely in terms of the future average, which is the only random piece left:

max(A_total - K, 0) = max(w_future * A_future - (K - w_real * realized_avg), 0)
                     = w_future * max(A_future - K_eff, 0),   K_eff = (K - w_real * realized_avg) / w_future

So a seasoned Asian option reduces exactly to an unseasoned one on the remaining n_obs fixings, priced against an adjusted strike K_eff and scaled by w_future. This is applied to both _moment_matching (with M1 and M2 computed only over the remaining fixings) and to Monte Carlo (simulating only the remaining path). If K_eff <= 0, the option is already guaranteed in the money (for a call) or guaranteed worthless (for a put) regardless of what the future fixings do, since the future average is always positive; the code detects this and returns the payoff directly rather than evaluating log(K_eff) on a non-positive number.

For the geometric average the combination is multiplicative rather than additive, but the same idea applies in log space:

ln(G_total) = w_real * ln(realized_avg) + w_future * ln(G_future)

Since ln(G_future) ~ Normal(m_future, v_future) from the unseasoned derivation above, ln(G_total) is normal too, with

m = w_real * ln(realized_avg) + w_future * m_future
v = w_future^2 * v_future

and the same closed form applies with these adjusted m and v.

If no fixings remain (n_obs = 0), there is nothing left to simulate or approximate: the average is exactly realized_avg, so a fixed strike payoff is just max(realized_avg - K, 0) discounted, and a floating strike payoff reduces to a plain vanilla option struck at realized_avg, priced by ordinary Black Scholes.

Cliquets: price_cliquet

A cliquet (ratchet) splits [0, T] into n_periods equal length reset periods. At the end of each period the return over that period is computed, clipped to a local floor and cap, and the clipped returns are combined, either summed or compounded, into a single payout at final maturity. No coupons, no early exercise: everything pays at T.

Decomposing the clip

For a single period with gross return G = S_{tau}/S_0' over that period (where S_0' is the period's starting price) and raw return R = G - 1, the clipped return is

clip(R, floor, cap) = min(max(R, floor), cap)

The identity that makes this tractable in closed form is that a clip can be written as a floor plus two call payoffs:

clip(R, floor, cap) = floor + max(R - floor, 0) - max(R - cap, 0)

Checking the three cases confirms it: if R < floor both max terms are 0 and the result is floor; if floor <= R <= cap the first max term equals R - floor and the second is 0, giving R; if R > cap the first term is R - floor and the second is R - cap, and the difference is cap - floor, giving floor + (cap - floor) = cap. Since R = G - 1, R - floor = G - (1 + floor) and R - cap = G - (1 + cap), so each max term is a call payoff on the gross return G, struck at 1 + floor and 1 + cap respectively.

Expected value of a call on the gross return

Under flat rate and volatility, G = S_tau / S_0' over a period of length tau is lognormal, by the same argument used for the Asian geometric average:

ln(G) ~ Normal(mu, v),   mu = (r - q - 0.5 sigma^2) tau,   v = sigma^2 tau

so the expected discounted-to-the-period-end call payoff on G struck at some level k is, again, a Black Scholes shaped expression with G's own expectation exp(mu + v/2) in place of spot:

d1 = (mu + v - ln(k)) / sqrt(v)
d2 = d1 - sqrt(v)
E[max(G - k, 0)] = exp(mu + v/2) * N(d1) - k * N(d2)

This is _raw_expected_call_payoff(tau, k). It needs v > 0, i.e. sigma > 0 and tau > 0; the degenerate case is checked explicitly rather than left to produce nan.

Assembling one period, then the whole cliquet

Combining the clip identity with the call formula, the expected clipped return for one period is

E[clip(R, floor, cap)] = floor + E[max(G - (1+floor), 0)] - E[max(G - (1+cap), 0)]

which is exactly call_floor - call_cap + floor in _closed_form, using _raw_expected_call_payoff evaluated at k = 1 + floor and k = 1 + cap.

Under additive compounding with flat parameters, every period has the same distribution (each is a fresh forward starting return over an identical length of time), so by linearity of expectation the total expected return is simply n_periods times the single period expectation:

E[total_return] = n_periods * E[clip(R, floor, cap)]

and the price is the notional times the discounted expected return:

price = notional * exp(-r T) * E[total_return]

This is exact given the model's assumptions. It stops being exact under two circumstances, both of which raise NotImplementedError rather than silently returning a wrong number:

Multiplicative compounding introduces a product of clipped returns across periods, and the expectation of a product is not, in general, the product of expectations once the returns are correlated through a shared discount and volatility structure in more general models; even here, where periods are independent, the closed form as written only aggregates the additive case.

A global floor or cap applied to the compounded total return couples the periods together: whether the global cap binds in period 5 depends on what happened in periods 1 through 4, which breaks the clean per-period decomposition entirely.

Monte Carlo

For compounding, or global floors and caps, the code falls back to direct simulation: draw one gross return per period per path, clip locally, combine (sum or compound) across periods, apply the global floor and cap if given, discount, and average. return_stderr=True returns the Monte Carlo standard error alongside the price, as with the Asian pricer.

Where flat volatility actually matters

Every closed form above, in every pricer, leans on sigma and r being constant across the life of the contract. This is a real modeling simplification, especially for the cliquet: each reset period is a forward starting option, and the price of a forward starting option is sensitive to forward volatility, the market's view of volatility over a future interval, not spot volatility today. A single flat sigma says forward vol equals spot vol for every future period, which is not what real volatility surfaces imply. A production desk pricing cliquets would calibrate a stochastic volatility model, such as Heston or SABR, to instruments that actually reveal forward vol and skew. This library prices correctly under its stated assumption; it does not attempt to relax that assumption.