#!/usr/bin/env python3 """Kelly sizing over the forecast vol range: a self-contained reproduction. Companion script to the Options Lab research page "Kelly sized over the range of volatility". It sizes an option spread with the Kelly criterion when realised volatility to expiry is uncertain, the way the site's strategy suggestions do: 1. The forecast: sizing sigma c = ATM implied vol minus a 1 vol-point premium. 2. The range around it: log RV ~ N(log c - s^2/2, s^2) with s = sqrt(0.064) = 0.253 at every horizon, centred so that c is the range's *mean*; the upper half is widened by k = 1.143 (the "tail haircut") so 7.5 % of the outcomes land above the plain range's 95th percentile instead of 5 %. 3. The smile: a terminal price distribution implied by the chain's IVs (Breeden-Litzenberger on a PCHIP smile, calibrated to the listed calls), scaled so its ATM vol is c. 4. Quadrature: 16 Gauss points on each side of the median. At each point the smile's log-returns are stretched to that vol. 5. Kelly: the log-optimal fraction of the bankroll to stake (the stake is the trade's max loss) on the probability-weighted mixture of all points. The recommendation is half of it, capped at 5 % of the portfolio. Run it with no arguments for a worked example on a synthetic 50ETF-like chain, or pass `--chain my_chain.json` (format: see `load_chain_json`). pip install numpy scipy python kelly_range.py python kelly_range.py --centre 0.1571 # just the range for one sigma Requirements: Python 3.10+, numpy >= 1.26, scipy >= 1.11. For learning only. Not investment advice; past results do not guarantee future ones. MIT License Copyright (c) 2026 Options Lab Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ from __future__ import annotations import argparse import json import math from dataclasses import dataclass from functools import lru_cache from typing import Sequence import numpy as np from numpy.polynomial.hermite_e import hermegauss from scipy.interpolate import PchipInterpolator from scipy.optimize import brentq from scipy.stats import norm # --------------------------------------------------------------------------- # Parameters (as used on the site, fitted on QVIX history 2015-2026) # --------------------------------------------------------------------------- FORECAST_PREMIUM = 0.01 # the forecast: ATM IV minus 1 vol point RANGE_ALPHA = 0.064 # log-variance of realised vol around the forecast LOG_SD = math.sqrt(RANGE_ALPHA) # 0.253, the same at every horizon TAIL_ABOVE_Q95 = 0.075 # share of outcomes above the plain range's 95th pct _Z95 = float(norm.ppf(0.95)) TAIL_K = _Z95 / float(norm.ppf(1.0 - TAIL_ABOVE_Q95)) # 1.143: upper-half widening QUADRATURE_NODES = 31 # 16 per side with the haircut (31 Gauss-Hermite without) VOL_FLOOR = 0.01 MAX_RATIO = 5.0 # nodes above 5x the centre are clipped to it KELLY_CEILING = 1.0 - 1e-6 # never stake the whole bankroll KELLY_CAP = 0.05 # recommended stake = min(half-Kelly, 5 % of the portfolio) MODERATE_FROM = 0.02 # verdict bands on the full fraction STRONG_FROM = 0.10 TRADING_DAYS_PER_YEAR = 252 # IV is quoted on the trading-day clock # Smile -> distribution GRID_POINTS = 1200 GRID_WIDTH_SD = 8.0 MAX_WING_SLOPE = 1.9 CALIBRATION_ITERATIONS = 6 _Z90 = float(norm.ppf(0.9)) # A leg: (option_type "call"/"put", signed quantity (+ long, - short), strike, # entry price per share) Leg = tuple[str, float, float, float] # --------------------------------------------------------------------------- # 1-2. The forecast and the range around it # --------------------------------------------------------------------------- def forecast_sigma(atm_iv: float, premium: float = FORECAST_PREMIUM) -> float: """The sizing sigma: ATM IV minus the long-run premium, floored at 1 %.""" return max(atm_iv - premium, VOL_FLOOR) @lru_cache(maxsize=16) def half_range_rule(n: int) -> tuple[np.ndarray, np.ndarray]: """n-point Gauss rule for the half-normal (weight e^(-z^2/2) on [0, inf)), weights normalised to sum 1. Built by the discretised Stieltjes procedure and the Golub-Welsch eigenvalue method.""" x, wx = np.polynomial.legendre.leggauss(400) z = 6.0 * (x + 1.0) dz = 6.0 * wx * np.exp(-0.5 * z * z) a, b = np.zeros(n), np.zeros(n) p_prev, p = np.zeros_like(z), np.ones_like(z) norm_prev = 1.0 for j in range(n): norm_j = float(dz @ (p * p)) a[j] = float(dz @ (z * p * p)) / norm_j b[j] = norm_j if j == 0 else norm_j / norm_prev p_prev, p = p, (z - a[j]) * p - (b[j] if j > 0 else 0.0) * p_prev norm_prev = norm_j off = np.sqrt(b[1:]) nodes, vecs = np.linalg.eigh(np.diag(a) + np.diag(off, 1) + np.diag(off, -1)) weights = vecs[0] ** 2 return nodes, weights / weights.sum() @dataclass(frozen=True) class VolRange: """Quadrature nodes (sigma_k, w_k) of the forecast range, sum w = 1.""" centre: float sd: float median: float sigmas: np.ndarray weights: np.ndarray tail: bool @property def upside_sd(self) -> float: return self.sd * (TAIL_K if self.tail else 1.0) @property def p10(self) -> float: return self.median * math.exp(-_Z90 * self.sd) @property def p90(self) -> float: return self.median * math.exp(_Z90 * self.upside_sd) def quantile(self, q: float) -> float: """Exact quantile of the (split) lognormal range.""" z = float(norm.ppf(q)) return self.median * math.exp(z * (self.upside_sd if z > 0 else self.sd)) @property def mean(self) -> float: return float(self.weights @ self.sigmas) def forecast_range( centre: float, *, sd: float = LOG_SD, tail: bool = True, nodes: int = QUADRATURE_NODES, ) -> VolRange: """The range Kelly sizes a forecast `centre` over. log RV ~ N(log c - s^2/2, s^2): the -s^2/2 shift makes c the mean. With `tail`, the upper half's log sd is s * TAIL_K (a split normal that keeps the median and the lower half) and each half gets its own Gauss rule, so the kink at the median is integrated exactly.""" if not centre > 0: raise ValueError(f"centre must be positive, got {centre!r}") median = centre * math.exp(-0.5 * sd * sd) if tail: u, wu = half_range_rule((nodes + 1) // 2) z = np.concatenate([-u[::-1], u * TAIL_K]) w = np.concatenate([wu[::-1], wu]) / 2.0 else: z, w = hermegauss(nodes) w = w / w.sum() sigmas = np.clip(median * np.exp(sd * z), VOL_FLOOR, MAX_RATIO * centre) if sd == 0: sigmas = np.full(len(z), float(centre)) return VolRange(float(centre), float(sd), median, sigmas, w, tail) # --------------------------------------------------------------------------- # 3. The smile and the terminal price distribution # --------------------------------------------------------------------------- def black76_call(forward: float, strikes: np.ndarray, sigma: np.ndarray, t: float) -> np.ndarray: """Undiscounted Black-76 call value E[(S_T - K)+].""" sd = np.maximum(sigma * math.sqrt(t), 1e-12) d1 = (np.log(forward / strikes) + 0.5 * sd * sd) / sd return forward * norm.cdf(d1) - strikes * norm.cdf(d1 - sd) def _on_lower_hull(strikes: np.ndarray, calls: np.ndarray, forward: float, far: float) -> np.ndarray: """Strikes whose call price is on the lower convex hull of (0, F), the listed points and (far, 0): no butterfly arbitrage. Others are dropped.""" points = [(0.0, forward), *zip(strikes.tolist(), calls.tolist()), (far, 0.0)] hull: list[tuple[float, float]] = [] for p in points: while len(hull) >= 2: (x1, y1), (x2, y2) = hull[-2], hull[-1] if (y2 - y1) * (p[0] - x1) >= (p[1] - y1) * (x2 - x1): hull.pop() else: break hull.append(p) kept = {x for x, _ in hull} return np.array([k in kept for k in strikes.tolist()], dtype=bool) def _smile_function(k: np.ndarray, iv: np.ndarray, t: float): """sigma(log-moneyness): PCHIP through the points, total variance linear beyond the outermost strikes (slope clamped to [0, 1.9] outward).""" if len(k) == 1: return lambda kg: np.full_like(kg, iv[0]) pchip = PchipInterpolator(k, iv, extrapolate=False) slope = pchip.derivative() w_lo, w_hi = iv[0] ** 2 * t, iv[-1] ** 2 * t s_lo = min(max(2 * iv[0] * float(slope(k[0])) * t, -MAX_WING_SLOPE), 0.0) s_hi = max(min(2 * iv[-1] * float(slope(k[-1])) * t, MAX_WING_SLOPE), 0.0) def sigma(kg: np.ndarray) -> np.ndarray: out = np.empty_like(kg) inner = (kg >= k[0]) & (kg <= k[-1]) out[inner] = pchip(kg[inner]) lo, hi = kg < k[0], kg > k[-1] out[lo] = np.sqrt(np.maximum(w_lo + s_lo * (kg[lo] - k[0]), 1e-12) / t) out[hi] = np.sqrt(np.maximum(w_hi + s_hi * (kg[hi] - k[-1]), 1e-12) / t) return out return sigma def _masses_from_calls(nodes: np.ndarray, calls: np.ndarray) -> np.ndarray: """Breeden-Litzenberger: node mass = change in the call curve's slope.""" slopes = np.diff(calls) / np.diff(nodes) masses = np.empty(len(nodes)) masses[0] = 1.0 + slopes[0] masses[1:-1] = np.diff(slopes) masses[-1] = -slopes[-1] return masses def _calibrate(nodes, masses, forward, strikes, calls): """Smallest relative correction of the masses so they sum to 1, have the forward as mean and reprice every kept strike's call (project and clip).""" a = np.vstack([np.ones_like(nodes), nodes, np.maximum(nodes[None, :] - strikes[:, None], 0.0)]) b = np.concatenate([[1.0, forward], calls]) best, best_residual = masses, float(np.abs(a @ masses - b).max()) p = masses for _ in range(CALIBRATION_ITERATIONS): w = p + 1e-12 gram = (a * w) @ a.T lam = np.linalg.lstsq(gram, b - a @ p, rcond=1e-14)[0] p = np.clip(p + w * (a.T @ lam), 0.0, None) p = p / p.sum() residual = float(np.abs(a @ p - b).max()) if residual < best_residual: best, best_residual = p, residual if residual <= 1e-9 * forward: break return best class Distribution: """A discrete distribution of the price at expiry: nodes and masses.""" def __init__(self, forward: float, t: float, nodes: np.ndarray, masses: np.ndarray, strikes: np.ndarray, ivs: np.ndarray): self.forward, self.t = forward, t self.nodes, self.masses = nodes, masses self.strikes, self.ivs = strikes, ivs @classmethod def from_smile(cls, forward: float, t: float, strikes: Sequence[float], ivs: Sequence[float]) -> "Distribution": """The distribution implied by one IV per strike.""" k_arr = np.asarray(strikes, dtype=np.float64) iv_arr = np.asarray(ivs, dtype=np.float64) ok = np.isfinite(k_arr) & np.isfinite(iv_arr) & (k_arr > 0) & (iv_arr > 0) k_arr, iv_arr = k_arr[ok], iv_arr[ok] order = np.argsort(k_arr) k_arr, iv_arr = k_arr[order], iv_arr[order] k_arr, first = np.unique(k_arr, return_index=True) iv_arr = iv_arr[first] width = GRID_WIDTH_SD * float(iv_arr.max()) * math.sqrt(t) k_all = np.log(k_arr / forward) far = forward * math.exp(max(float(k_all.max()), 0.0) + width) listed = black76_call(forward, k_arr, iv_arr, t) keep = _on_lower_hull(k_arr, listed, forward, far) sigma = _smile_function(k_all[keep], iv_arr[keep], t) lo = min(float(k_all.min()), 0.0) - width hi = max(float(k_all.max()), 0.0) + width grid_k = np.union1d(np.linspace(lo, hi, GRID_POINTS), k_all) nodes = forward * np.exp(grid_k) nodes[np.searchsorted(grid_k, k_all)] = k_arr # strikes exactly on the grid calls = black76_call(forward, nodes, sigma(grid_k), t) masses = np.clip(_masses_from_calls(nodes, calls), 0.0, None) masses = masses / masses.sum() masses = _calibrate(nodes, masses, forward, k_arr[keep], listed[keep]) return cls(forward, t, nodes, masses, k_arr, iv_arr) @classmethod def lognormal(cls, forward: float, t: float, sigma: float) -> "Distribution": """A flat smile.""" return cls.from_smile(forward, t, [forward], [sigma]) def scaled(self, ratio: float) -> "Distribution": """The same smile with every IV multiplied by `ratio` (rebuilt).""" return Distribution.from_smile(self.forward, self.t, self.strikes, self.ivs * ratio) def return_scaled(self, lam: float) -> "Distribution": """Same masses, every log-return log(S/F) stretched by `lam` and re-centred on the forward: a realised-vol scenario that keeps the smile's skew and kurtosis. Arbitrage-free for every lam > 0.""" if lam == 1.0: return self stretched = np.exp(lam * np.log(self.nodes / self.forward)) nodes = self.forward * stretched / float(self.masses @ stretched) return Distribution(self.forward, self.t, nodes, self.masses, self.strikes, self.ivs * lam) def pnl(self, legs: Sequence[Leg]) -> np.ndarray: """Expiry P&L per share at every node.""" out = np.zeros_like(self.nodes) for option_type, qty, strike, price in legs: payoff = np.maximum(self.nodes - strike, 0.0) if option_type == "call" \ else np.maximum(strike - self.nodes, 0.0) out += qty * (payoff - price) return out def model_value(self, option_type: str, strike: float) -> float: payoff = np.maximum(self.nodes - strike, 0.0) if option_type == "call" \ else np.maximum(strike - self.nodes, 0.0) return float(self.masses @ payoff) # --------------------------------------------------------------------------- # 4-5. Kelly over the range # --------------------------------------------------------------------------- def expiry_pnl(legs: Sequence[Leg], price: float) -> float: total = 0.0 for option_type, qty, strike, entry in legs: payoff = max(price - strike, 0.0) if option_type == "call" else max(strike - price, 0.0) total += qty * (payoff - entry) return total def max_loss_per_share(legs: Sequence[Leg]) -> float | None: """The largest loss at expiry (>= 0), or None if the loss is unlimited (net short calls: the P&L keeps falling as the price rises).""" net_calls = sum(qty for option_type, qty, _, _ in legs if option_type == "call") if net_calls < 0: return None strikes = sorted({k for _, _, k, _ in legs}) worst = min(expiry_pnl(legs, s) for s in [0.0, *strikes, 2.0 * strikes[-1]]) return max(-worst, 0.0) def log_optimal_fraction(pnl: np.ndarray, weights: np.ndarray, max_loss: float) -> float: """f* = argmax_f sum_i w_i log(1 + f X_i / L), f in [0, 1). X is the P&L per share at each outcome, w its probability, L the max loss (the stake). 0 when the expected P&L is <= 0.""" if not max_loss > 0: return 0.0 pnl, weights = np.asarray(pnl, float), np.asarray(weights, float) keep = weights > 0 ratio, weights = pnl[keep] / max_loss, weights[keep] if ratio.size == 0: return 0.0 weights = weights / weights.sum() if not float(weights @ ratio) > 0: return 0.0 def slope(f: float) -> float: # derivative of expected log growth return float(weights @ (ratio / (1.0 + f * ratio))) upper = KELLY_CEILING / max(1.0, -float(ratio.min())) if slope(upper) >= 0: return upper return float(brentq(slope, 0.0, upper, xtol=1e-12)) def verdict(fraction: float) -> str: if fraction <= 0: return "negative" if fraction < MODERATE_FROM: return "marginal" if fraction < STRONG_FROM: return "moderate" return "strong" @dataclass class Sizing: full_kelly: float half_kelly_capped: float verdict: str expected_pnl: float # per share, on the distribution the fraction used max_loss: float | None contracts: int = 0 def kelly_sizing( market: Distribution, atm_iv: float, centre: float, legs: Sequence[Leg], *, vol_range: VolRange | None = None, multiplier: int = 10_000, portfolio: float = 1_000_000, ) -> Sizing: """Size `legs`. `market` is the chain's smile at its own ATM IV; it is scaled so its ATM vol is the sizing sigma `centre`. With `vol_range` (centre must match) the fraction is the log-optimal one on the mixture of every quadrature node; without it, on the one point forecast.""" loss = max_loss_per_share(legs) if loss is None: return Sizing(0.0, 0.0, "unlimited_risk", float("nan"), None) point = market.scaled(centre / atm_iv) if vol_range is None: pnl, weights = point.pnl(legs), point.masses else: if abs(vol_range.centre - centre) > 1e-12: raise ValueError("the range must be centred on the sizing sigma") pnls, ws = [], [] for sigma_k, w_k in zip(vol_range.sigmas, vol_range.weights): node = point.return_scaled(float(sigma_k) / centre) pnls.append(node.pnl(legs)) ws.append(w_k * node.masses) pnl, weights = np.concatenate(pnls), np.concatenate(ws) f = log_optimal_fraction(pnl, weights, loss) half = min(f / 2.0, KELLY_CAP) contracts = 0 if half > 0 and loss > 0: raw = half * portfolio / (loss * multiplier) contracts = 1 if 0 < raw < 1 else int(raw) return Sizing(f, half, verdict(f), float(weights @ pnl / weights.sum()), loss, contracts) def kelly_if_vol_fixed(point: Distribution, centre: float, legs: Sequence[Leg], sigma: float) -> float: """Full Kelly if realised vol were known to be `sigma`: the point distribution stretched to it (one node of the mixture).""" loss = max_loss_per_share(legs) if loss is None: return 0.0 node = point.return_scaled(sigma / centre) return log_optimal_fraction(node.pnl(legs), node.masses, loss) # --------------------------------------------------------------------------- # A synthetic chain for the worked example # --------------------------------------------------------------------------- def synthetic_chain(forward: float = 3.05, atm_iv: float = 0.14, trading_days: int = 15, step: float = 0.05, n_side: int = 10, skew: float = -0.30, curvature: float = 1.5) -> tuple[Distribution, np.ndarray, np.ndarray]: """A 50ETF-like chain: IV(k) = ATM + skew*k + curvature*k^2 with k = log(K/F), strikes every `step` around the forward.""" t = trading_days / TRADING_DAYS_PER_YEAR centre = round(forward / step) * step strikes = np.round(centre + step * np.arange(-n_side, n_side + 1), 4) k = np.log(strikes / forward) ivs = np.maximum(atm_iv + skew * k + curvature * k * k, 0.03) return Distribution.from_smile(forward, t, strikes, ivs), strikes, ivs def quote(market: Distribution, option_type: str, strike: float, direction: str, half_spread: float = 0.0002) -> float: """Execution price: the smile's model value, plus half a bid-ask spread when buying, minus it when selling.""" value = market.model_value(option_type, strike) return value + half_spread if direction == "long" else max(value - half_spread, 0.0001) def leg(market: Distribution, option_type: str, direction: str, strike: float) -> Leg: qty = 1.0 if direction == "long" else -1.0 return (option_type, qty, strike, quote(market, option_type, strike, direction)) def load_chain_json(path: str) -> dict: """Your own chain as JSON: {"forward": 3.05, "trading_days": 15, "atm_iv": 0.14, "strikes": [2.8, 2.85, ...], "ivs": [0.17, 0.16, ...], "multiplier": 10000, "legs": [["put", 1, 2.7, 0.0040], ["put", -1, 2.9, 0.0195]]} IVs as decimals, one per strike (the out-of-the-money side's quote); legs are [type, signed quantity, strike, entry price per share].""" with open(path, encoding="utf-8") as f: return json.load(f) # --------------------------------------------------------------------------- # Report # --------------------------------------------------------------------------- def print_range(centre: float) -> None: r = forecast_range(centre) print(f" sigma {centre:.2%}: 80 % range {r.p10:.1%} - {r.p90:.1%}, " f"median {r.median:.2%}, 95th pct {r.quantile(0.95):.1%}, " f"mean of the range {r.mean:.2%} ({r.mean / centre - 1:+.1%} vs sigma)") def print_sizing(name: str, market: Distribution, atm_iv: float, legs: Sequence[Leg], multiplier: int, curve_vols: Sequence[float]) -> None: centre = forecast_sigma(atm_iv) text = " ".join(f"{'+' if q > 0 else '-'}{k:g}{'C' if t == 'call' else 'P'} @ {p:.4f}" for t, q, k, p in legs) print(f"\n{name}: {text}") print(f" ATM IV {atm_iv:.2%}, sizing sigma {centre:.2%}") point = kelly_sizing(market, atm_iv, centre, legs, multiplier=multiplier) mixed = kelly_sizing(market, atm_iv, centre, legs, vol_range=forecast_range(centre), multiplier=multiplier) if mixed.max_loss is None: print(" unlimited risk: Kelly sizes 0") return for label, s in (("one number (point forecast)", point), ("over the vol range", mixed)): print(f" {label:28s} full Kelly {s.full_kelly:.4f} half-Kelly (cap 5 %) " f"{s.half_kelly_capped:.2%} {s.verdict:8s} contracts/1M {s.contracts:3d} " f"EV/contract {s.expected_pnl * multiplier:+.1f}") scaled = market.scaled(centre / atm_iv) curve = [kelly_if_vol_fixed(scaled, centre, legs, v) for v in curve_vols] print(" full Kelly if realised vol were fixed at:") print(" " + " ".join(f"{v:>6.0%}" for v in curve_vols)) print(" " + " ".join(f"{f:6.3f}" for f in curve)) def main() -> int: ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) ap.add_argument("--centre", type=float, help="only print the range for this sizing sigma (decimal)") ap.add_argument("--chain", help="size the legs of your own chain (JSON, see load_chain_json)") args = ap.parse_args() print(f"Range: log sd {LOG_SD:.3f} at every horizon, upper half x {TAIL_K:.3f} " f"({TAIL_ABOVE_Q95:.1%} above the plain 95th percentile)") if args.centre: print_range(args.centre) return 0 for centre in (0.1281, 0.1571, 0.2488): print_range(centre) if args.chain: spec = load_chain_json(args.chain) t = spec["trading_days"] / TRADING_DAYS_PER_YEAR market = Distribution.from_smile(spec["forward"], t, spec["strikes"], spec["ivs"]) legs = [tuple(x) for x in spec["legs"]] print_sizing("Your legs", market, spec["atm_iv"], legs, spec.get("multiplier", 10_000), (0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.25, 0.30)) return 0 atm_iv = 0.14 market, _, _ = synthetic_chain(atm_iv=atm_iv) vols = (0.08, 0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.25, 0.30) print("\nSynthetic 50ETF-like chain: forward 3.05, 15 trading days, ATM IV 14 %, put skew") print_sizing("Bull put spread, near the money", market, atm_iv, [leg(market, "put", "long", 2.85), leg(market, "put", "short", 3.00)], 10_000, vols) print_sizing("Bull put spread, far out of the money", market, atm_iv, [leg(market, "put", "long", 2.70), leg(market, "put", "short", 2.90)], 10_000, vols) print_sizing("Iron condor", market, atm_iv, [leg(market, "put", "long", 2.85), leg(market, "put", "short", 2.95), leg(market, "call", "short", 3.15), leg(market, "call", "long", 3.25)], 10_000, vols) return 0 if __name__ == "__main__": raise SystemExit(main())