Save this block as pop.py and run python3 pop.py. It uses only the standard library. Signed leg quantities determine the piecewise-linear expiry P&L; the code solves each interval between strikes and integrates all profitable pieces with the normal CDF. Adjacent pieces may share a strike endpoint, which does not double-count probability in the continuous case. Edit EXAMPLES and the explicit model inputs to inspect another same-expiry position.
Standalone expiry calculation. The six example outputs and sensitivity rows reproduce the tables above; no network data or hidden defaults for missing contract fields.python
from math import erfc, exp, inf, isfinite, log, sqrt
# Leg: ("call" or "put", signed contract quantity, strike, premium per share).
# All legs have one underlying, one expiry and the same verified multiplier.
def validate(legs, multiplier, fees, slippage):
if not legs or not all(isfinite(x) for x in (multiplier, fees, slippage)):
raise ValueError("Supply legs and finite costs/multiplier")
if multiplier <= 0 or fees < 0 or slippage < 0:
raise ValueError("Multiplier must be positive; costs nonnegative")
for kind, quantity, strike, premium in legs:
if kind not in ("call", "put") or not all(
isfinite(x) for x in (quantity, strike, premium)
):
raise ValueError("Invalid leg")
if quantity == 0 or quantity != int(quantity) or strike <= 0 or premium < 0:
raise ValueError("Use nonzero integer quantity, positive strike, nonnegative premium")
def pnl(price, legs, multiplier, fees, slippage):
result = sum(
q * multiplier * (max(price - k if kind == "call" else k - price, 0) - p)
for kind, q, k, p in legs
) - fees - slippage
if not isfinite(result):
raise ValueError("P&L exceeds numeric range")
return result
def profit_intervals(legs, multiplier=100, fees=0, slippage=0):
validate(legs, multiplier, fees, slippage)
knots = [0.0] + sorted({leg[2] for leg in legs}) + [inf]
intervals = []
# Payoff is affine between strikes; solve a*S + b > 0 on every segment.
for lo, hi in zip(knots, knots[1:]):
a, b = 0.0, -fees - slippage
for kind, q, k, premium in legs:
scale = q * multiplier
b -= scale * premium
if kind == "call" and k <= lo:
a += scale
b -= scale * k
elif kind == "put" and k >= hi:
a -= scale
b += scale * k
if not isfinite(a) or not isfinite(b):
raise ValueError("Payoff coefficients exceed numeric range")
if a > 0:
lo = max(lo, -b / a)
elif a < 0:
hi = min(hi, -b / a)
elif b <= 0:
continue
if lo < hi:
intervals.append((lo, hi))
return intervals
def probability(legs, spot=100, years=30/365, sigma=0.20, drift=0.02,
multiplier=100, fees=0, slippage=0, cash_tolerance=1e-8):
validate(legs, multiplier, fees, slippage)
if not isfinite(cash_tolerance) or cash_tolerance < 0:
raise ValueError("Cash tolerance must be finite and nonnegative")
if not all(isfinite(x) for x in (spot, years, sigma, drift)):
raise ValueError("Model inputs must be finite")
if spot <= 0 or years < 0 or sigma < 0:
raise ValueError("Spot must be positive; time and volatility nonnegative")
if years == 0 or sigma == 0:
try:
terminal = spot if years == 0 else spot * exp(drift * years)
except OverflowError as error:
raise ValueError("Terminal price exceeds numeric range") from error
if not isfinite(terminal) or terminal <= 0:
raise ValueError("Terminal price exceeds numeric range")
return float(pnl(terminal, legs, multiplier, fees, slippage) > cash_tolerance)
center = log(spot) + (drift - sigma*sigma/2) * years
scale = sigma * sqrt(years)
if not isfinite(center) or not isfinite(scale) or scale <= 0:
raise ValueError("Distribution exceeds numeric range")
def z(price):
return -inf if price == 0 else (log(price) - center) / scale
def mass(lo, hi):
a, b = z(lo), z(hi)
# Survival probabilities preserve precision in the right tail.
if a >= 0:
return 0.5 * (erfc(a/sqrt(2)) - erfc(b/sqrt(2)))
return 0.5 * (erfc(-b/sqrt(2)) - erfc(-a/sqrt(2)))
return sum(mass(lo, hi) for lo, hi in profit_intervals(
legs, multiplier, fees, slippage))
EXAMPLES = {
"Long call": ([("call", 1, 100, 3)], 2, 3),
"Long put": ([("put", 1, 100, 3)], 2, 3),
"Short put": ([("put", -1, 100, 3)], 2, 3),
"Bull call spread": ([("call", 1, 100, 5), ("call", -1, 110, 2)], 4, 6),
"Iron condor": ([("put", 1, 90, .5), ("put", -1, 95, 1.5),
("call", -1, 105, 1.5), ("call", 1, 110, .5)], 8, 12),
"Long straddle": ([("call", 1, 100, 3), ("put", 1, 100, 3)], 4, 6),
}
if __name__ == "__main__":
for name, (legs, fees, slippage) in EXAMPLES.items():
result = probability(legs, fees=fees, slippage=slippage)
print(f"{name}: {100*result:.4f}%")
call = EXAMPLES["Long call"][0]
for sigma, drift, fees, slippage in [
(.20, .02, 0, 0), (.20, .02, 2, 3), (.20, .02, 2, 23),
(.15, .02, 2, 3), (.30, .02, 2, 3), (.20, .08, 2, 3),
]:
result = probability(call, sigma=sigma, drift=drift, fees=fees, slippage=slippage)
print(f"sigma={sigma}, drift={drift}, costs={fees+slippage}: {100*result:.4f}%")
