API reference

HestonParams(κ, θ, ξ, ρ, v0)

κ mean-reversion speed, θ long-run variance, ξ vol-of-vol, ρ spot-vol correlation, v0 initial variance. Type-parametric (promoting constructor) so ForwardDiff.Dual parameters flow through pricing and calibration.

RoughHestonParams(κ, θ, ξ, ρ, v0, H)

Classical Heston parameters plus the Hurst exponent H ∈ (0, 1/2). H = 1/2 is classical Heston. Type-parametric for ForwardDiff.

SSVIParams(ρ, η, γ)

Global SSVI parameters with power-law φ(θ) = η/(θ^γ (1+θ)^{1−γ}).

SVIParams(a, b, ρ, m, σ)

Raw-SVI slice: w(k) = a + b(ρ(k−m) + √((k−m)² + σ²)) in total variance.

ad_greeks(price, S, r, q, T)

Delta, gamma, theta and rho of any pricer price(S, r, q, T) by ForwardDiff (gamma is a nested derivative). price must be generic in all four arguments. Returns (price, delta, gamma, theta, rho).

atm_total_variance(ks, ivs, T)

ATM total variance σatm²·T, with σatm linearly interpolated at k = 0 from the two quotes bracketing it (nearest quote if k = 0 is not bracketed).

batch_call_prices(ψ, F, disc, Ks, T; iv_hint=0.2, reltail=12.0, osc_pts=8.0,
                  min_nodes=64, max_nodes=1024)

Gil-Pelaez call prices for ALL strikes Ks of one expiry from shared CF evaluations on fixed Gauss-Legendre nodes. ψ is the de-drifted-log-return CF (must accept complex arguments); F the forward, disc = e^{−rT}. Puts follow from exact parity: P = C − disc·(F − K).

bs_delta(S, K, r, q, σ, T; call=true)

∂Price/∂S. Call: e^{−qT}·Φ(d₁). Put: e^{−qT}·(Φ(d₁) − 1) (by parity ∂C/∂S − ∂P/∂S = e^{−qT}).

Degenerate limits return the delta of the intrinsic payoff; the exact-ATM boundary is assigned to the in-the-money side by convention here.

bs_gamma(S, K, r, q, σ, T)

∂²Price/∂S² = e^{−qT}·φ(d₁) / (S·σ·√T). Same for calls and puts (parity is linear in S). Zero in the degenerate limits T ≤ 0 or σ ≤ 0, where the payoff is piecewise linear in S away from the kink.

bs_price(S, K, r, q, σ, T; call=true)

European option price under BSM with continuous carry yield q.

d₁ = (log(S/K) + (r − q + σ²/2)·T) / (σ·√T),   d₂ = d₁ − σ·√T
call = S·e^{−qT}·Φ(d₁) − K·e^{−rT}·Φ(d₂)

Degenerate limits are handled explicitly rather than letting d₁ hit 0/0:

• T ≤ 0 → immediate intrinsic: max(S−K, 0) or max(K−S, 0).
• σ ≤ 0 → deterministic world: discounted intrinsic on the forward F = S·e^{(r−q)T}, i.e. call = e^{−rT}·max(F − K, 0).

bs_rho(S, K, r, q, σ, T; call=true)

∂Price/∂r. Call: K·T·e^{−rT}·Φ(d₂). Put: −K·T·e^{−rT}·Φ(−d₂) (by parity ∂C/∂r − ∂P/∂r = K·T·e^{−rT}). Degenerate limits follow bs_theta.

bs_theta(S, K, r, q, σ, T; call=true)

Time decay ∂Price/∂t = −∂Price/∂T, per year (divide by 365 for per-day).

call: −S·e^{−qT}·φ(d₁)·σ/(2√T) − r·K·e^{−rT}·Φ(d₂) + q·S·e^{−qT}·Φ(d₁)
put:  −S·e^{−qT}·φ(d₁)·σ/(2√T) + r·K·e^{−rT}·Φ(−d₂) − q·S·e^{−qT}·Φ(−d₁)

σ ≤ 0 (T > 0) differentiates the discounted forward intrinsic; T ≤ 0 returns zero. Degenerate cases put the exact-ATM boundary on the in-the-money side, as bs_delta does.

bs_vega(S, K, r, q, σ, T)

∂Price/∂σ = S·e^{−qT}·φ(d₁)·√T. Same for calls and puts: by parity C − P has no σ-dependence, so ∂C/∂σ = ∂P/∂σ.

Vega → 0 both deep ITM and deep OTM (φ(d₁) → 0 as |d₁| → ∞) and as T → 0. That decay is exactly what breaks naive Newton in implied_vol.

business_days(d0, d1; holidays=nothing)

Trading days in the half-open interval (d0, d1]: weekdays that are not holidays. holidays = nothing uses nyse_holidays for every year spanned.

calendar_violations(ws, Ts; ks=range(-1, 1, length=201), tol=1e-10)

Calendar-arbitrage check across slices: ws[i] is a function k ↦ total variance for maturity Ts[i]. Returns the (T_short, T_long, k) triples where total variance DEcreases with maturity; empty means calendar-free.

calibrate_heston(quotes, S; p0=HestonParams(2.0, 0.04, 0.5, -0.5, 0.04),
                 method=:lbfgs, pricer=:adaptive, maxiter=300,
                 price_rtol=1e-7)

Fit classical Heston to a prepared chain with ForwardDiff derivatives. method is :lbfgs (LBFGS on the mean squared residual) or :lm (Levenberg–Marquardt on the residual vector). pricer is :adaptive (per-quote quadgk, price_rtol applies) or :batch (per-expiry fixed-node batch pricer; quotes without an F field get F = S·e^{(r−q)T}). Returns (params, rmse, converged, iterations); rmse is in implied-vol units (multiply by 1e4 for bps).

calibrate_rough_heston(quotes; p0=RoughHestonParams(1.5,0.05,0.35,-0.65,0.011,0.12),
                       method=:lbfgs, maxiter=60, N=96, S=nothing)

Fit rough Heston (6 parameters incl. H) to a prepared chain with ForwardDiff derivatives through the fractional solver. method is :lbfgs or :lm, as in calibrate_heston; S is only needed for quotes without an F field. Same return shape as calibrate_heston.

calibrate_rough_heston_joint(quote_sets; p0=RoughHestonParams(1.5,0.05,0.35,-0.65,0.011,0.12),
                             skew_weight=0.0, maxiter=60, N=96)

Fit rough Heston jointly to several snapshots: (κ, θ, ξ, ρ, H) shared, v0 per snapshot (initialized from each snapshot's shortest-expiry ATM variance). quote_sets is a vector of prepare_chain-style quote vectors (fields T, K, F, r, q, side, mid, iv, vega); a single snapshot is the D = 1 case, where skew_weight > 0 alone adds the skew-slope term.

Returns (params, rmse, skew_rmse, converged, iterations): params holds one RoughHestonParams per snapshot; rmse is the price-residual IV RMSE and skew_rmse the RMS scaled skew mismatch, both in vol units.

feller_ratio(p)

2κθ/ξ². The Feller condition 2κθ ≥ ξ² (ratio ≥ 1) keeps the variance process strictly positive; below 1, v can touch zero. Calibration does not enforce it (see src/calibration.jl) — this is the reporting side. Works for HestonParams and RoughHestonParams alike (same κ, θ, ξ fields).

fit_ssvi(slices; p0=SSVIParams(-0.6, 1.0, 0.4), maxiter=500)

Fit SSVI to slices, a vector of (T, ks, ivs, θ) with θ the slice's ATM total variance (see atm_total_variance). The θ curve is made non-decreasing (running maximum) before fitting, and the parametrization keeps η(1+|ρ|) ≤ 2 and γ ∈ (0, 1/2], so the result is free of static arbitrage by construction. Returns (params, θs, rmse) with rmse the IV RMSE over all quotes.

fit_svi(ks, ivs, T; weights=nothing, maxiter=500)

Fit a raw-SVI slice to implied vols ivs at log-moneyness ks by weighted least squares on implied vol (so the RMSE is comparable to the Heston IV RMSEs). Multi-start LBFGS over a small grid of (m, σ) seeds. Returns (params, rmse); rmse is the (weighted) IV RMSE.

fit_svi_surface(quotes; minquotes=5)

Fit raw SVI per expiry and SSVI globally to prepare_chain-style quotes (fields T, K, F, iv). Returns (slices, ssvi, svi_rmse, ssvi_rmse, butterfly_free, calendar_violations): slices holds (T, params, rmse, n) per expiry, svi_rmse/ssvi_rmse are pooled IV RMSEs, butterfly_free[i] flags each SVI slice, and calendar_violations lists calendar-arbitrage points between consecutive SVI slices.

frac_integral_end(h, r, Δ)

Riemann–Liouville fractional integral I^r h evaluated at the END of the uniform grid carrying h (spacing Δ), by product-trapezoidal quadrature. r = 1 is the trapezoid rule; r → 0 is the identity (returns h[end]).

group_quotes(quotes; S=nothing)

Group prepare_chain-style quotes by expiry for batch pricing. Put mids are converted to call-equivalent mids via exact parity (C = P + disc·(F − K); vega is identical either side). Returns one NamedTuple per expiry with (T, F, r, q, disc, iv_atm, Ks, cmids, vegas).

Quotes without an F field (e.g. hand-built synthetic sets) need the spot S; their forward is then S·e^{(r−q)T}.

heston_batch_loss(x, groups)

Mean of heston_batch_residuals(x, groups)².

heston_batch_residuals(x, groups)

Classical-Heston residuals batch-priced per expiry from group_quotes output — the same fixed-node engine as rough Heston, ~all strikes of an expiry for the price of one CF sweep.

heston_cf(u, T, p::HestonParams)

ψ(u) = E[exp(iu·XT)] for the de-drifted log return XT (see the convention block above). Accepts real or complex u (the pricer evaluates at u − i).

Requires ξ > 0; the ξ → 0 deterministic-variance limit is exercised in tests at ξ = 1e-3 against its closed form rather than at ξ = 0 exactly.

heston_greeks(S, K, r, q, p::HestonParams, T; call=true, rtol=1e-10)

Classical Heston Greeks by AD through heston_price. Returns (price, delta, gamma, theta, rho, vega, sens) where vega = ∂V/∂√v0 and sens = (κ, θ, ξ, ρ, v0) holds ∂V/∂(each model parameter).

heston_loss(x, quotes, S; price_rtol=1e-7)

Mean squared vega-normalized price residual (≈ mean squared IV error) at transformed parameter vector x. quotes need fields (T, K, r, q, side, mid, vega) — prepare_chain output qualifies. Type-generic in x so ForwardDiff Duals flow through.

heston_price(S, K, r, q, p::HestonParams, T; call=true, rtol=1e-9)

Classical Heston price: heston_cf plugged into price_from_cf. The params argument sits where σ sits in bs_price — same call shape, model swapped.

heston_residuals(x, quotes, S; price_rtol=1e-7)

Vega-normalized price residuals (C_model − mid)/vega (≈ IV errors) for classical Heston at transformed parameters x, one per quote, priced by adaptive quadrature. Type-generic in x (ForwardDiff Duals flow through).

implied_vol(price, S, K, r, q, T; call=true, tol=1e-10, maxiter=200)

Invert BSM for σ given an observed price.

Method: Newton's method using the analytic vega, safeguarded by a maintained bisection bracket [lo, hi]. Rationale:

• Price is strictly increasing in σ (vega > 0 for σ, T > 0), so a valid bracket exists and bisection alone would already converge — slowly.
• Newton is quadratically fast near the root but diverges where vega ≈ 0 (deep OTM/ITM, short T): the objective is nearly flat there, so the Newton step Δσ = (model − market)/vega explodes.
• Hybrid: take the Newton step only if it lands strictly inside the current bracket; otherwise bisect. Every iteration shrinks the bracket or Newton-converges, so the method is globally convergent and locally quadratic. (Same philosophy as Brent's method.)

Bounds (with carry yield q): call: max(S·e^{−qT} − K·e^{−rT}, 0) ≤ C < S·e^{−qT} put: max(K·e^{−rT} − S·e^{−qT}, 0) ≤ P < K·e^{−rT} Prices outside these bounds have NO implied vol — that's an input error, not a convergence failure, so we throw a DomainError. (At Stage 4 this becomes the first data-quality filter: stale SPX quotes that violate these bounds land here.) A price at (within tol of) the lower bound corresponds to σ = 0.

levenberg_marquardt(rfun, x0; maxiter=100, λ0=1e-3, gtol=1e-10,
                    xtol=1e-10, ftol=1e-14)

Minimize sum(abs2, rfun(x)) by Levenberg–Marquardt. rfun must be generic in x (the Jacobian is ForwardDiff.jacobian(rfun, x)). Each iteration solves (J'J + λ·diag(J'J)) δ = −J'r; λ shrinks on accepted steps and grows on rejected ones. Converged when the gradient, the step, or the relative decrease of the objective falls below its tolerance.

Returns (minimizer, minimum, iterations, converged) where minimum is the sum of squares.

make_rough_cf(T, p::RoughHestonParams; N=128)

Return the characteristic function u -> ψ(u) with all (α, T, N)-dependent tables precomputed — use this when evaluating many u for one (T, params), e.g. across a quadrature grid.

model_atm_skew(ψ, g)

Linearized model ATM skew ∂σ/∂k for characteristic function ψ on expiry group g (see the file header). AD-safe.

normal_cdf(x)

Φ(x), via the error function: Φ(x) = (1 + erf(x/√2)) / 2. erf has ForwardDiff rules, so Duals flow through.

normal_pdf(x)

φ(x) = exp(−x²/2) / √(2π).

nyse_holidays(year)

Full-day NYSE closures for year under the current rule set: New Year's Day, MLK Day, Presidents' Day, Good Friday, Memorial Day, Juneteenth (from 2022), Independence Day, Labor Day, Thanksgiving, Christmas — with the weekend observance rules (a Saturday New Year's Day is NOT observed on the prior Friday, per NYSE rule). Ad-hoc closures (national days of mourning, weather) are not included; pass them to year_fraction via holidays.

prepare_chain(raw; valuation_date=raw.quote_date, min_days=2, max_years=2.5,
              max_logm=0.40, min_quotes=5, atm_band=0.10,
              daycount=:calendar, holidays=nothing)

Turn a parsed CBOE file into a calibration-ready chain. Returns (quotes, expiries, rejects, spot, valuation_date):

• quotes: one entry per surviving OTM quote — (expiry, T, K, F, r, q, side, bid, ask, mid, iv, vega), with iv from Stage 1's inverter and vega evaluated at that iv (the calibration weight).
• expiries: per-expiry (expiry, T, F, r, q, n) from the parity regression.
• rejects: counted reasons for every discarded row — the data-quality report.

valuation_date matters: weekend/holiday snapshots carry the last session's closing quotes, and for short-dated options the day count is a first-order effect. Pass the last trading date explicitly when the snapshot date isn't it.

daycount = :business measures T in NYSE trading days / 252 (see year_fraction); min_days stays in calendar days. An expiry with zero trading days left is counted under :out_of_window.

prepare_identification_set(quotes; minpts=4)

Group one snapshot's quotes by expiry and precompute each expiry's skew stencil. Returns (groups, skews, stencils): skews[i] is the market ATM skew (nothing where fewer than minpts quotes sit in the ATM window) and stencils[i] the (idx, w) pair turning residuals into a skew mismatch, Σ w·r[idx] (the least-squares slope on k), scaled by σ_atm·√T.

price_from_cf(ψ, S, K, r, q, T; call=true, rtol=1e-9)

European option price by Gil-Pelaez inversion of the characteristic function ψ(u) = E[e^{iuX}], X = log(S_T/F), ψ(−i) = 1. ψ must accept complex arguments (it is evaluated at u − i). Model-independent by construction.

The put comes via parity, P = C − e^{−rT}(F − K), which is exact here because both legs share the same two integrals.

read_cboe(path)

Parse a CBOE delayed-quotes CSV. Returns (spot, quote_date, rows) where each row is (expiry, K, cbid, cask, pbid, pask). Purely structural — no filtering, no judgment; that's prepare_chain's job.

rough_heston_cf(u, T, p::RoughHestonParams; N=128)

One-shot ψ(u). For many u at fixed (T, p) prefer make_rough_cf.

rough_heston_greeks(S, K, r, q, p::RoughHestonParams, T; call=true, N=128)

Rough Heston Greeks by AD through the fractional Riccati solver and the batch pricer. Same return shape as heston_greeks, with H added to sens.

rough_heston_loss(x, groups; N=96)

Mean squared vega-normalized residual for rough Heston at transformed parameters x, batch-priced per expiry group.

rough_heston_mc_price(S, K, r, q, p::RoughHestonParams, T; call=true,
                      nsteps=200, npaths=50_000, rng=nothing)

Monte Carlo price of a European option under rough Heston, with the forward as control variate. Returns (price, stderr).

rough_heston_residuals(x, groups; N=96)

Vega-normalized residuals for rough Heston at transformed parameters x, batch-priced per expiry group.

rough_joint_loss(x, datasets; skew_weight=0.0, N=96, parts=false)

Loss for a joint rough-Heston fit over snapshot datasets (each from prepare_identification_set) at packed parameters x = [log κ, log θ, log ξ, atanh ρ, logit H, log v0₁, …, log v0_D]:

mean_i r_i²  +  skew_weight · mean_e [ (skew_model − skew_mkt)·σ_atm·√T ]²

with ri the vega-normalized price residuals and the skew mismatch taken as the slope of r on k over each expiry's ATM window (see the file header). parts = true returns `(pricemse, skew_mse)` instead of the combined value.

simulate_rough_heston(p::RoughHestonParams, T, nsteps, npaths; rng=nothing)

Simulate npaths paths of the de-drifted log return XT = log(ST/F) under rough Heston with the hybrid Volterra-Euler scheme (see the file header). Returns the vector of X_T. rng is any RNG accepted by randn(rng) (nothing uses the global RNG). Requires H ∈ (0, 1/2].

skew_powerlaw_H(Ts, skews)

Fit |skew| = c·T^{H−1/2} by least squares in log-log space. Returns (H, c). The model-free roughness estimate of scripts/04.

solve_fractional_riccati(c0, c1, c2, α, T, N)

Solve D^α h = c0 + c1 h + c2 h², h(0) = 0 on t_k = kT/N, k = 0..N by the implicit fractional product-trapezoidal method (exact per-step quadratic solve). Returns the trajectory h::Vector (length N+1). Fully generic: coefficients and/or α may be ForwardDiff.Duals.

ssvi_arbitrage_free(p::SSVIParams, θs)

Gatheral–Jacquier sufficient conditions: θ non-decreasing, γ ∈ (0, 1/2], η(1 + |ρ|) ≤ 2.

ssvi_slice(p::SSVIParams, θ)

The raw-SVI slice equal to SSVI at ATM total variance θ, so every SVI tool (svi_iv, svi_butterfly_g, …) applies to SSVI slices too.

ssvi_total_variance(p::SSVIParams, θ, k)

SSVI total variance at ATM total variance θ and log-moneyness k.

svi_butterfly_free(p::SVIParams; ks=range(-2, 2, length=801), tol=0)

true when Durrleman's g(k) ≥ −tol on the grid ks and w > 0 there.

svi_butterfly_g(p::SVIParams, k)

Durrleman's function g(k); the slice is free of butterfly arbitrage iff g ≥ 0 everywhere.

svi_iv(p::SVIParams, k, T)

Implied volatility √(w(k)/T).

svi_total_variance(p::SVIParams, k)

Total implied variance w(k) = σ_iv²·T at log-moneyness k.

year_fraction(d0, d1; daycount=:calendar, holidays=nothing)

Time from d0 to d1 in years. :calendar is days/365; :business is business_days(d0, d1; holidays)/252.