Metadata-Version: 2.4
Name: exotics
Version: 1.0.0
Summary: Pricing library for barrier, Asian, and cliquet options under Black-Scholes dynamics
Author-email: jr1concepcion@gmail.com
License-Expression: MIT
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Intended Audience :: Financial and Insurance Industry
Classifier: Topic :: Office/Business :: Financial
Classifier: Operating System :: OS Independent
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Dynamic: license-file

# 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.
