Pricing ( wraquant.price ) ¶

Example

>>> black_scholes(100, 100, 1.0, 0.05, 0.2)
10.450...
>>> black_scholes(100, 100, 1.0, 0.05, 0.2, 'put')
5.573...

Notes

The model assumes constant volatility and no dividends. For dividend-paying stocks, reduce S by the present value of expected dividends or use the Merton (1973) continuous-dividend extension (replace S with S * exp(-q * T)).

See also

binomial_tree: Lattice method supporting American exercise. monte_carlo_option: Simulation-based pricing. implied_volatility: Invert this formula to recover vol from price.

References

Black, F. & Scholes, M. (1973). The Pricing of Options and Corporate Liabilities. Journal of Political Economy 81(3).

binomial_tree(S, K, T, r, sigma, n_steps, option_type='call', style='european')[source]¶

Price an option using the Cox-Ross-Rubinstein (CRR) binomial tree.

The binomial tree discretises the asset price evolution into up/down moves and computes the option value by backward induction. Unlike Black-Scholes, the tree supports American-style early exercise. Use this when you need American option prices or when the analytical formula is unavailable.

At each step the asset moves up by u = exp(sigma * sqrt(dt)) or down by d = 1/u with risk-neutral probability p = (exp(r*dt) - d) / (u - d). The tree converges to the Black-Scholes price as n_steps increases.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized, continuously compounded).
sigma (float) – Volatility of the underlying (annualized).
n_steps (int) – Number of time steps in the tree. 200–500 steps give good accuracy for most applications.
option_type (str | OptionType, default: 'call') – 'call' or 'put'.
style (str | OptionStyle, default: 'european') – 'european' or 'american'.

Returns:

The binomial tree option price. For European options: this converges to the Black-Scholes price; for American options the price is at least as large as the European price due to the early exercise premium.

Return type:

Example

>>> binomial_tree(100, 100, 1.0, 0.05, 0.2, 200)
10.4...
>>> binomial_tree(100, 100, 1.0, 0.05, 0.2, 200, 'put', 'american')
6.08...

Notes

Convergence is O(1/n_steps). For faster convergence consider Richardson extrapolation (average the n and 2n step prices).

See also

black_scholes: Analytical European option pricing. monte_carlo_option: Simulation-based pricing.

References

Cox, J.C., Ross, S.A. & Rubinstein, M. (1979). Option Pricing: A Simplified Approach. Journal of Financial Economics 7(3).

monte_carlo_option(S, K, T, r, sigma, n_sims, n_steps, option_type='call', seed=None)[source]¶

Price a European option using Monte Carlo simulation under GBM.

Simulates n_sims asset price paths under risk-neutral geometric Brownian motion, computes the terminal payoff on each path, and returns the discounted average. Use Monte Carlo when the payoff is path-dependent, the underlying follows a non-standard process, or when you need a flexible framework that can be extended to exotic options.

The standard error of the estimate decreases as \(1/\sqrt{n\_sims}\), so 100 000 paths typically give 2–3 significant digits of accuracy.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized, continuously compounded).
sigma (float) – Volatility of the underlying (annualized).
n_sims (int) – Number of simulation paths. More paths reduce variance but increase computation time.
n_steps (int) – Number of time steps per path. For European vanilla options a single step suffices, but path-dependent payoffs require finer discretisation.
option_type (str | OptionType, default: 'call') – 'call' or 'put'.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

The Monte Carlo estimate of the option price. This: is an unbiased estimator that converges to the Black-Scholes price for European vanilla options under GBM.

Return type:

Example

>>> monte_carlo_option(100, 100, 1.0, 0.05, 0.2, 100000, 252, seed=42)
10.4...

Notes

For variance reduction consider antithetic variates or control variates (using the Black-Scholes price as a control).

See also

black_scholes: Closed-form European pricing (faster, exact). binomial_tree: Lattice method for American options.

delta(S, K, T, r, sigma, option_type='call')[source]¶

Compute the Black-Scholes delta (sensitivity to underlying price).

Delta measures how much the option price changes for a one-unit change in the underlying price. It is the first derivative of the option price with respect to S and also represents the hedge ratio (number of shares to hold per option to be delta-neutral).

\[\Delta_{\text{call}} = \Phi(d_1), \quad \Delta_{\text{put}} = \Phi(d_1) - 1\]

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Delta value. Call delta is in [0, 1]; put delta: is in [-1, 0]. At-the-money forward options have delta near 0.5 (call) or -0.5 (put).

Return type:

Example

>>> delta(100, 100, 1.0, 0.05, 0.2, 'call')
0.636...
>>> delta(100, 100, 1.0, 0.05, 0.2, 'put')
-0.363...

See also

gamma: Rate of change of delta. all_greeks: Compute all Greeks in a single call.

gamma(S, K, T, r, sigma)[source]¶

Compute the Black-Scholes gamma (rate of change of delta).

Gamma is the second derivative of option price with respect to S. It measures the convexity of the option position and is identical for calls and puts. High gamma means the hedge ratio changes rapidly, requiring frequent rebalancing.

\[\Gamma = \frac{\phi(d_1)}{S\,\sigma\,\sqrt{T}}\]

Gamma is highest for at-the-money options near expiry.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).

Returns:

Gamma value (always non-negative). Expressed in: price units per unit squared move in the underlying.

Return type:

Example

>>> gamma(100, 100, 1.0, 0.05, 0.2)
0.018...

See also

delta: First-order sensitivity to the underlying. all_greeks: Compute all Greeks in a single call.

theta(S, K, T, r, sigma, option_type='call')[source]¶

Compute the Black-Scholes theta (time decay per year).

Theta measures the rate at which the option loses value as time passes, holding all else constant. It is typically negative for long option positions (time works against the holder). Divide by 252 to get the approximate daily time decay.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Theta value per year. Typically negative for long: options. For daily theta, divide by 252 (or 365 for calendar days).

Return type:

Example

>>> theta(100, 100, 1.0, 0.05, 0.2, 'call')
-6.41...

See also

vega: Sensitivity to volatility. all_greeks: Compute all Greeks in a single call.

vega(S, K, T, r, sigma)[source]¶

Compute the Black-Scholes vega (sensitivity to volatility).

Vega measures how much the option price changes for a one-unit (100 percentage point) change in implied volatility. It is identical for calls and puts and is always non-negative. To get the price change for a 1 percentage point vol move, multiply the result by 0.01.

\[\mathcal{V} = S\,\phi(d_1)\,\sqrt{T}\]

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).

Returns:

Vega per unit change in sigma (always non-negative).: Multiply by 0.01 for the price change per 1% vol move.

Return type:

Example

>>> vega(100, 100, 1.0, 0.05, 0.2)
37.52...

See also

implied_volatility: Invert the BS formula using vega as the: Newton-Raphson derivative.

all_greeks: Compute all Greeks in a single call.

rho(S, K, T, r, sigma, option_type='call')[source]¶

Compute the Black-Scholes rho (sensitivity to interest rate).

Rho measures how much the option price changes for a one-unit (100 percentage point) change in the risk-free rate. Multiply by 0.01 for the price change per 1% rate move. Rho is typically the smallest of the Greeks for short-dated equity options but becomes significant for long-dated or fixed-income options.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Rho per unit change in r. Positive for calls,: negative for puts.

Return type:

Example

>>> rho(100, 100, 1.0, 0.05, 0.2, 'call')
53.23...

See also

all_greeks: Compute all Greeks in a single call.

all_greeks(S, K, T, r, sigma, option_type='call')[source]¶

Compute all Black-Scholes Greeks in a single efficient call.

Computes d1 and d2 once and derives all five Greeks, avoiding the redundant computation that would occur when calling each Greek function individually. Use this when you need a full risk profile for a position.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Dictionary with keys 'delta',: 'gamma', 'theta', 'vega', 'rho'. See the individual Greek functions for interpretation of each value.

Return type:

Example

>>> greeks = all_greeks(100, 100, 1.0, 0.05, 0.2, 'call')
>>> sorted(greeks.keys())
['delta', 'gamma', 'rho', 'theta', 'vega']

See also

delta, gamma, theta, vega, rho: Individual Greek computations.

implied_volatility(market_price, S, K, T, r, option_type='call', tol=1e-08, max_iter=200, initial_guess=0.3)[source]¶

Compute implied volatility using Newton-Raphson iteration.

Solves for \(\sigma\) such that \(\text{BS}(S, K, T, r, \sigma) = \text{market\_price}\) by iterating \(\sigma_{n+1} = \sigma_n - (C(\sigma_n) - C_{\text{mkt}}) / \mathcal{V}(\sigma_n)\), where \(\mathcal{V}\) is the Black-Scholes vega.

Use implied volatility to translate observed option prices into a single number that summarises the market’s view of future uncertainty. Comparing implied vols across strikes and expiries reveals the volatility smile/skew.

When vega is too small (deep in/out of the money), the solver automatically falls back to a robust bisection method.

Parameters:

market_price (float) – Observed market price of the option.
S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.
tol (float, default: 1e-08) – Convergence tolerance on the price difference.
max_iter (int, default: 200) – Maximum number of Newton iterations.
initial_guess (float, default: 0.3) – Starting volatility estimate. 0.3 (30%) is a reasonable default for most equity options.

Returns:

Implied volatility (annualized). For example, 0.20: corresponds to 20% annualized volatility.

Return type:

Raises:

PricingError – If the solver does not converge within max_iter iterations.

Example

>>> price = black_scholes(100, 100, 1.0, 0.05, 0.2)
>>> implied_volatility(price, 100, 100, 1.0, 0.05)
0.2000...

See also

vol_smile: Compute implied vols across multiple strikes. vol_surface: Build a full implied volatility surface.

vol_smile(strikes, market_prices, S, T, r, option_type='call')[source]¶

Compute a volatility smile from market prices at multiple strikes.

For a fixed expiry, the implied volatility typically varies across strikes, forming a “smile” or “skew” pattern. This function computes the implied vol at each strike by inverting the Black-Scholes formula, producing the raw data for smile analysis, model calibration, or visualization.

Parameters:

strikes (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of strike prices.
market_prices (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Corresponding market option prices (same length as strikes).
S (float) – Current underlying price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Dictionary with keys:

'strikes' – 1D array of strike prices.
'implied_vols' – 1D array of implied volatilities corresponding to each strike.

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Raises:

ValueError – If strikes and market_prices differ in length.
PricingError – If any individual implied vol fails to converge.

Example

>>> strikes = [90, 95, 100, 105, 110]
>>> prices = [12.5, 9.0, 6.0, 3.8, 2.1]
>>> result = vol_smile(strikes, prices, 100, 0.5, 0.05)
>>> len(result['implied_vols'])
5

See also

implied_volatility: Single-strike implied vol solver. vol_surface: Extend the smile across multiple expiries.

vol_surface(strikes, expiries, market_prices, S, r, option_type='call')[source]¶

Construct a volatility surface from market prices across strikes and expiries.

The volatility surface is the central object in options trading: it maps every (strike, expiry) pair to an implied volatility. Use the surface for interpolation to price off-market options, for calibrating stochastic volatility models (Heston, SABR), or for detecting relative-value opportunities.

Parameters:

strikes (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of strike prices (columns of the surface).
expiries (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of times to expiry in years (rows of the surface).
market_prices (Sequence[Sequence[float]] | ndarray[tuple[Any, ...], dtype[floating]]) – 2D array of option prices, shape (len(expiries), len(strikes)).
S (float) – Current underlying price.
r (float) – Risk-free interest rate (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Dictionary with keys:

'strikes' – 1D array of strikes.
'expiries' – 1D array of expiries.
'implied_vols' – 2D array of implied vols, shape (len(expiries), len(strikes)).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Raises:

ValueError – If the shape of market_prices does not match (len(expiries), len(strikes)).

Example

>>> strikes = [95, 100, 105]
>>> expiries = [0.25, 0.5]
>>> prices = [[6.0, 3.5, 1.5], [8.0, 5.5, 3.5]]
>>> surface = vol_surface(strikes, expiries, prices, 100, 0.05)
>>> surface['implied_vols'].shape
(2, 3)

See also

vol_smile: Single-expiry volatility smile. implied_volatility: Single-point implied vol solver.

bond_price(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the price of a fixed-rate bond as the PV of its cash flows.

Discounts each coupon payment and the face-value redemption at the yield to maturity. When ytm == coupon_rate, the bond prices at par. Use this to mark positions to market given a yield, or to compute the theoretical price for comparison against the quoted market price.

\[P = \sum_{t=1}^{n} \frac{C}{(1+y)^t} + \frac{F}{(1+y)^n}\]

where \(C = F \cdot c / f\) is the periodic coupon, \(y = \text{ytm} / f\), \(n\) is the number of periods, and \(f\) is the payment frequency.

Parameters:

face_value (float) – Par/face value of the bond (e.g., 1000).
coupon_rate (float) – Annual coupon rate (e.g., 0.05 for 5%).
ytm (float) – Yield to maturity (annualized, e.g., 0.05 for 5%).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Present value (dirty price) of the bond. When: ytm < coupon_rate the bond trades at a premium (price > face); when ytm > coupon_rate it trades at a discount.

Return type:

Example

>>> bond_price(1000, 0.05, 0.05, 10, 2)
1000.0
>>> bond_price(1000, 0.05, 0.06, 10, 2)
925.6...

See also

bond_yield: Solve for the yield given price. duration: Macaulay duration for interest-rate sensitivity.

bond_yield(price, face_value, coupon_rate, periods, freq=2, tol=1e-10, max_iter=200)[source]¶

Compute the yield to maturity (YTM) of a bond via Newton’s method.

YTM is the single discount rate that equates the present value of a bond’s cash flows to its market price. It is the internal rate of return assuming all coupons are reinvested at the same rate. Use YTM to compare bonds with different coupons and maturities on a common basis.

Parameters:

price (float) – Market price of the bond.
face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).
tol (float, default: 1e-10) – Convergence tolerance on the price residual.
max_iter (int, default: 200) – Maximum Newton-Raphson iterations.

Returns:

Annualized yield to maturity. For a par bond: (price == face_value) the YTM equals the coupon rate.

Return type:

Raises:

ValueError – If the solver does not converge within max_iter iterations.

Example

>>> bond_yield(950, 1000, 0.05, 10, 2)
0.059...

See also

bond_price: Price a bond given a yield. zero_rate: Extract a zero rate from a zero-coupon bond.

duration(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the Macaulay duration of a fixed-rate bond.

Macaulay duration is the weighted-average time to receive the bond’s cash flows, where the weights are the present-value fractions. It measures the bond’s sensitivity to parallel shifts in the yield curve: a duration of 5 years means a 1% yield increase causes approximately a 5% price decrease.

\[D = \frac{1}{P} \sum_{t=1}^{n} \frac{t}{f} \cdot \frac{CF_t}{(1 + y)^t}\]

Parameters:

face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
ytm (float) – Yield to maturity (annualized).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Macaulay duration in years. Zero-coupon bonds have: duration equal to their maturity; coupon bonds always have duration less than maturity.

Return type:

Example

>>> duration(1000, 0.05, 0.05, 10, 2)
4.08...

See also

modified_duration: Duration divided by (1 + y/freq). convexity: Second-order interest-rate sensitivity.

modified_duration(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the modified duration of a fixed-rate bond.

Modified duration is the more practically useful sensitivity measure: the percentage price change for a one-unit change in yield.

\[D_{\text{mod}} = \frac{D_{\text{mac}}}{1 + y/f}\]

A modified duration of 4 means a 100 bp (1%) yield increase causes approximately a 4% price decline.

Parameters:

face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
ytm (float) – Yield to maturity (annualized).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Modified duration in years.

Return type:

Example

>>> modified_duration(1000, 0.05, 0.05, 10, 2)
3.98...

See also

duration: Macaulay duration. convexity: Second-order correction for large yield changes.

convexity(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the convexity of a fixed-rate bond.

Convexity is the second-order sensitivity of bond price to yield changes. Duration alone gives a linear approximation; convexity adds the curvature correction, improving accuracy for large yield moves.

The price change for a yield shift \(\Delta y\) is approximately:

\[\frac{\Delta P}{P} \approx -D_{\text{mod}} \cdot \Delta y + \tfrac{1}{2}\,C \cdot (\Delta y)^2\]

where \(C\) is convexity. Positive convexity means the bond gains more from a yield decrease than it loses from an equal yield increase.

Parameters:

face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
ytm (float) – Yield to maturity (annualized).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Convexity measure (in years squared, scaled by: freq^2).

Return type:

Example

>>> convexity(1000, 0.05, 0.05, 10, 2)
19.4...

See also

duration: First-order interest-rate sensitivity. modified_duration: Percentage price change per unit yield change.

zero_rate(price, face_value, periods)[source]¶

Compute the zero-coupon (spot) rate from a zero-coupon bond price.

Solves \(P = F \cdot e^{-r \cdot T}\) for r, giving the continuously compounded rate that equates the bond price to the present value of the face value. Zero rates are the building blocks of the term structure; use them to bootstrap a yield curve or to discount future cash flows.

Parameters:

price (float) – Current market price of the zero-coupon bond.
face_value (float) – Par/face value of the bond.
periods (float) – Time to maturity in years.

Returns:

Annualized continuously compounded zero rate.

Return type:

Raises:

ValueError – If price or periods is not positive.

Example

>>> zero_rate(950, 1000, 1.0)
0.051...

See also

bootstrap_zero_curve: Bootstrap a full zero curve from par rates. discount_factor: Convert a zero rate to a discount factor.

bootstrap_zero_curve(maturities, par_rates, freq=2)[source]¶

Bootstrap zero (spot) rates from par coupon rates.

Bootstrapping is the standard technique for constructing a zero-coupon yield curve from observed par bond yields. Starting from the shortest maturity, each successive zero rate is extracted by stripping previously bootstrapped discount factors from the coupon payments.

The resulting zero rates can be used to discount arbitrary cash flows, compute forward rates, or price swaps and bonds.

Parameters:

maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of maturities in years, evenly spaced at 1/freq (e.g., 0.5, 1.0, 1.5 for semiannual).
par_rates (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of par coupon rates (annualized) for each maturity.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Array of continuously compounded zero rates: corresponding to each maturity.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Raises:

ValueError – If maturities and par_rates differ in length.

Example

>>> mats = [0.5, 1.0, 1.5, 2.0]
>>> pars = [0.04, 0.042, 0.044, 0.046]
>>> zeros = bootstrap_zero_curve(mats, pars, freq=2)
>>> len(zeros)
4

See also

forward_rate: Derive forward rates from zero rates. discount_factor: Convert a zero rate to a discount factor.

interpolate_curve(maturities, rates, target_maturities, method='cubic')[source]¶

Interpolate a yield curve at target maturities.

Yield curves are observed at discrete maturities but often needed at arbitrary points. This function supports three interpolation strategies, each with different smoothness and no-arbitrage properties.

Parameters:

maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Known maturities (sorted ascending).
rates (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Known rates at those maturities.
target_maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Maturities at which to interpolate.
method (str, default: 'cubic') –
Interpolation method:
- 'linear' – piecewise linear (fast, may have kinks).
- 'cubic' – natural cubic spline (smooth, default).
- 'flat_forward' – constant forward rate between knots (no-arbitrage, standard in fixed income).

Returns:

Interpolated rates at target maturities.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Raises:

ValueError – If method is not recognised.

Example

>>> mats = [0.5, 1.0, 2.0, 5.0]
>>> rates = [0.03, 0.035, 0.04, 0.045]
>>> interpolate_curve(mats, rates, [0.75, 1.5, 3.0])
array([...])

See also

bootstrap_zero_curve: Build the curve from par rates. forward_rate: Extract a forward rate from the curve.

forward_rate(zero_rates, maturities, t1, t2)[source]¶

Compute the forward rate between two future dates.

The forward rate is the rate that can be locked in today for borrowing or lending between two future dates. Under no-arbitrage:

\[f(t_1, t_2) = \frac{r_2\,t_2 - r_1\,t_1}{t_2 - t_1}\]

Uses cubic spline interpolation to obtain zero rates at t1 and t2 from the provided curve.

Parameters:

zero_rates (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of continuously compounded zero rates.
maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Corresponding maturities.
t1 (float) – Start of the forward period (years).
t2 (float) – End of the forward period (years).

Returns:

Continuously compounded forward rate between t1: and t2.

Return type:

Raises:

ValueError – If t2 <= t1.

Example

>>> forward_rate([0.04, 0.045], [1.0, 2.0], 1.0, 2.0)
0.05

See also

bootstrap_zero_curve: Build the zero curve from par rates. discount_factor: Convert a rate to a discount factor.

discount_factor(rate, maturity)[source]¶

Compute the discount factor from a continuously compounded rate.

The discount factor is the present value of one unit of currency received at a future date:

\[DF = e^{-r \cdot T}\]

Discount factors are the fundamental building blocks for pricing any fixed-income instrument.

Parameters:

rate (float) – Continuously compounded interest rate.
maturity (float) – Time to maturity in years.

Returns:

Discount factor in (0, 1] for positive rates.

Return type:

Example

>>> discount_factor(0.05, 1.0)
0.9512...
>>> discount_factor(0.0, 5.0)
1.0

See also

zero_rate: Extract a rate from a zero-coupon bond price. bootstrap_zero_curve: Build a full discount factor curve.

vg_european_fft(spot, strike, rf_rate, T, sigma, nu, theta, n_fft=4096)[source]¶

Price a European call under the Variance Gamma model via FFT.

Combines the VG characteristic function with the Carr-Madan FFT method to price a European call option. The VG model captures skewness and excess kurtosis in the return distribution through three parameters.

Use this when the Black-Scholes assumption of Gaussian returns is too restrictive and you observe heavier tails or asymmetry in the market-implied distribution.

Parameters:

spot (float) – Current spot price.
strike (float) – Strike price.
rf_rate (float) – Risk-free rate (annualised).
T (float) – Time to maturity in years.
sigma (float) – VG volatility parameter (controls overall dispersion).
nu (float) – Variance rate of the Gamma subordinator (controls kurtosis; as nu -> 0 the model converges to Black-Scholes).
theta (float) – Drift of the VG process (controls skewness; negative = left-skewed returns).
n_fft (int, default: 4096) – FFT grid size (default 4096).

Returns:

European call price.

Return type:

Example

>>> price = vg_european_fft(100, 100, 0.05, 1.0, 0.2, 0.5, -0.1)
>>> price > 0
True

See also

nig_european_fft: NIG model pricing via FFT. fft_option_price: Generic FFT pricer for any characteristic fn. vg_characteristic: Reusable VG characteristic function.

nig_european_fft(spot, strike, rf_rate, T, alpha, beta, mu, delta, n_fft=4096)[source]¶

Price a European call under the Normal Inverse Gaussian model via FFT.

Combines the NIG characteristic function with the Carr-Madan FFT method. NIG features semi-heavy tails (heavier than Gaussian, lighter than Cauchy) and is widely used in credit, commodity, and FX markets.

Parameters:

spot (float) – Current spot price.
strike (float) – Strike price.
rf_rate (float) – Risk-free rate (annualised).
T (float) – Time to maturity in years.
alpha (float) – Tail heaviness parameter (alpha > 0). Larger values produce lighter tails.
beta (float) – Skewness parameter (-alpha < beta < alpha). Negative beta produces left-skewed returns.
mu (float) – Location parameter.
delta (float) – Scale parameter (delta > 0).
n_fft (int, default: 4096) – FFT grid size (default 4096).

Returns:

European call price.

Return type:

Example

>>> price = nig_european_fft(100, 100, 0.05, 1.0, 15.0, -3.0, 0.0, 0.5)
>>> price > 0
True

See also

vg_european_fft: Variance Gamma model pricing via FFT. fft_option_price: Generic FFT pricer for any characteristic fn. nig_characteristic: Reusable NIG characteristic function.

fft_option_price(char_fn, spot, strike, rf_rate, T, n_fft=4096, alpha_damp=1.5)[source]¶

Price European call option(s) via the Carr-Madan FFT method.

The Carr-Madan (1999) approach uses the Fast Fourier Transform to efficiently evaluate option prices across a grid of log-strikes in a single FFT pass. The method works with any model for which the characteristic function of \(\ln(S_T)\) is known, making it the universal workhorse for Fourier-based option pricing.

The core idea: multiply the modified characteristic function by Simpson’s rule weights, apply the FFT, and interpolate to the desired strike(s).

Parameters:

char_fn (Callable[[ndarray], ndarray]) – Characteristic function of \(\ln(S_T)\) under the risk-neutral measure. Signature: char_fn(u) -> complex ndarray.
spot (float) – Current spot price.
strike (float | ndarray) – Strike price(s). Accepts a scalar or an array for vectorised pricing.
rf_rate (float) – Risk-free rate (annualised, continuously compounded).
T (float) – Time to maturity (years).
n_fft (int, default: 4096) – Number of FFT points (default 4096). Must be a power of 2 for optimal FFT performance.
alpha_damp (float, default: 1.5) – Carr-Madan dampening parameter (default 1.5). Controls the integrand’s decay; values in [1, 2] work well for most models.

Returns:

European call price(s). Returns a float when: strike is scalar, otherwise an ndarray.

Return type:

float | ndarray

Example

>>> import numpy as np
>>> sigma, r, T = 0.2, 0.05, 1.0
>>> log_S = np.log(100.0)
>>> drift = (r - 0.5 * sigma**2) * T
>>> char_fn = lambda u: np.exp(1j*u*(log_S + drift) - 0.5*sigma**2*T*u**2)
>>> price = fft_option_price(char_fn, 100.0, 100.0, r, T)
>>> 9.0 < price < 12.0
True

See also

cos_method: Alternative Fourier method (faster for single strikes). characteristic_function_price: Unified interface to FFT and COS.

References

Carr, P. & Madan, D.B. (1999). Option Valuation Using the Fast Fourier Transform. Journal of Computational Finance 2(4).

cos_method(char_fn, spot, strike, rf_rate, T, n_terms=256, L=10.0)[source]¶

Price a European call option using the COS (Fourier-cosine) method.

The COS method of Fang & Oosterlee (2008) expands the risk-neutral density in a cosine series and integrates analytically against the payoff coefficients. It is faster than the FFT method for pricing a single option because it avoids the full FFT pass and works directly with the characteristic function at specific frequencies.

The price is computed as:

\[C = e^{-rT} \sum_{k=0}^{N-1} ' \text{Re}\bigl[\phi(u_k)\, e^{-iu_k a}\bigr]\,V_k\]

where \(u_k = k\pi/(b-a)\) and \(V_k\) are the cosine coefficients of the call payoff.

Parameters:

char_fn (Callable[[ndarray], ndarray]) – Characteristic function of \(\ln(S_T)\). Signature: char_fn(u) -> complex ndarray.
spot (float) – Current spot price.
strike (float) – Strike price.
rf_rate (float) – Risk-free rate (annualised).
T (float) – Time to maturity in years.
n_terms (int, default: 256) – Number of cosine expansion terms (default 256). 64–256 terms suffice for most models.
L (float, default: 10.0) – Truncation range controlling the integration domain [log(spot) - L, log(spot) + L] (default 10).

Returns:

European call price.

Return type:

See also

fft_option_price: FFT-based pricing (better for strike grids). characteristic_function_price: Unified interface.

References

Fang, F. & Oosterlee, C.W. (2008). A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions. SIAM Journal on Scientific Computing 31(2).

characteristic_function_price(char_fn, spot, strike, rf, T, method='fft', n_terms=4096)[source]¶

Price a European call option given a characteristic function.

This is the universal interface for Fourier-based option pricing. Any model whose characteristic function of log(S_T) is known in closed form can be priced through this single function, unifying Black-Scholes, Heston, Variance Gamma, NIG, CGMY, and others under one consistent API.

The connection to the pricing PDE is via the Feynman-Kac theorem: the characteristic function encodes the full risk-neutral distribution of the log-price, and the Fourier inversion recovers the density (or directly the option price via the Carr-Madan or COS method).

Parameters:

char_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]) – Characteristic function of log(S_T) under the risk-neutral measure. Signature: char_fn(u) -> complex ndarray where u is a real-valued array of frequencies.
spot (float) – Current spot price of the underlying.
strike (float) – Strike price.
rf (float) – Risk-free interest rate (annualised, continuously compounded).
T (float) – Time to maturity in years.
method (str, default: 'fft') – Pricing method – "fft" (Carr-Madan, default) or "cos" (Fang-Oosterlee cosine expansion).
n_terms (int, default: 4096) – Number of terms/grid points. For FFT this is the FFT size (default 4096); for COS this is the number of cosine terms.

Return type:

Returns:

European call option price.

Example

>>> import numpy as np
>>> # Black-Scholes via characteristic function
>>> S, K, r, T, sigma = 100.0, 100.0, 0.05, 1.0, 0.2
>>> log_S = np.log(S)
>>> drift = (r - 0.5 * sigma**2) * T
>>> char_fn = lambda u: np.exp(1j*u*(log_S + drift) - 0.5*sigma**2*T*u**2)
>>> price = characteristic_function_price(char_fn, S, K, r, T)
>>> 9.0 < price < 12.0
True

heston_characteristic(v0, kappa, theta, sigma_v, rho, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the Heston model.

The Heston (1993) stochastic volatility model specifies:

\[\begin{split}dS_t &= r\,S_t\,dt + \sqrt{v_t}\,S_t\,dW_t^1 \\ dv_t &= \kappa(\theta - v_t)\,dt + \sigma_v\sqrt{v_t}\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

The characteristic function of \(\ln(S_T)\) has a known closed-form expression involving complex exponentials of the model parameters.

When to use Heston vs Black-Scholes:

Use Black-Scholes when the implied vol surface is flat (no smile/skew) or as a quick approximation.
Use Heston when you observe a volatility smile/skew in the market and need to match observed option prices across strikes. The negative \(\rho\) parameter captures the leverage effect (negative correlation between returns and vol).
Heston is the standard model for equity index options and is widely used for calibration to the vol surface.

Parameters:

v0 (float) – Initial variance \(v_0\).
kappa (float) – Mean reversion speed of variance.
theta (float) – Long-run variance level.
sigma_v (float) – Volatility of variance (vol of vol).
rho (float) – Correlation between price and variance Brownian motions.
rf (float) – Risk-free rate.
T (float) – Time to maturity in years.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) that returns the characteristic function evaluated at frequencies u.

Example

>>> char_fn = heston_characteristic(
...     v0=0.04, kappa=2.0, theta=0.04, sigma_v=0.3,
...     rho=-0.7, rf=0.05, T=1.0, spot=100.0,
... )
>>> import numpy as np
>>> val = char_fn(np.array([0.0]))
>>> np.abs(val[0])  # phi(0) = 1
1.0

References

Heston, S.L. (1993). A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial Studies 6(2), 327-343.

vg_characteristic(sigma, nu, theta_vg, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the Variance Gamma model.

The Variance Gamma (VG) process models log-returns as a Brownian motion evaluated at a random (Gamma-distributed) time. This introduces both skewness and excess kurtosis into the return distribution while remaining analytically tractable.

\[\phi(u) = \exp\bigl(iu(r + \omega)T\bigr) \cdot \left(1 - iu\,\theta\,\nu + \tfrac{1}{2}\sigma^2\nu\,u^2\right)^{-T/\nu}\]

where \(\omega = \frac{1}{\nu}\ln(1 - \theta\nu - \sigma^2\nu/2)\) is the convexity correction.

Parameters:

sigma (float) – Volatility parameter of the VG process.
nu (float) – Variance rate of the Gamma subordinator. Controls kurtosis; as nu -> 0 the model converges to Black-Scholes.
theta_vg (float) – Drift parameter of the VG process. Controls skewness; negative values produce left-skewed returns (typical for equity markets).
rf (float) – Risk-free rate.
T (float) – Time to maturity.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) returning the characteristic function.

Example

>>> char_fn = vg_characteristic(0.2, 0.5, -0.1, 0.05, 1.0, spot=100)
>>> import numpy as np
>>> np.abs(char_fn(np.array([0.0]))[0])
1.0

nig_characteristic(alpha, beta, mu, delta, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the NIG model.

The Normal Inverse Gaussian (NIG) distribution arises as a normal variance-mean mixture with inverse Gaussian mixing distribution. It features semi-heavy tails (heavier than Gaussian, lighter than Cauchy) and is popular for modelling credit spreads, commodity returns, and FX markets.

\[\phi(u) = \exp\bigl(iu(r + \omega + \mu)T + \delta T\bigl(\gamma_0 - \sqrt{\alpha^2 - (\beta + iu)^2}\bigr)\bigr)\]

where \(\gamma_0 = \sqrt{\alpha^2 - \beta^2}\) and \(\omega = \delta(\sqrt{\alpha^2 - (\beta+1)^2} - \gamma_0) - \mu\) ensures the martingale condition.

Parameters:

alpha (float) – Tail heaviness / steepness parameter (alpha > 0). Larger alpha = lighter tails.
beta (float) – Asymmetry parameter (-alpha < beta < alpha). Negative beta = left skew.
mu (float) – Location parameter.
delta (float) – Scale parameter (delta > 0).
rf (float) – Risk-free rate.
T (float) – Time to maturity.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) returning the characteristic function.

Example

>>> char_fn = nig_characteristic(15.0, -3.0, 0.0, 0.5, 0.05, 1.0,
...                              spot=100)
>>> import numpy as np
>>> np.abs(char_fn(np.array([0.0]))[0])
1.0

cgmy_characteristic(C, G, M, Y_param, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the CGMY model.

The CGMY model (Carr, Geman, Madan & Yor, 2002) generalises the Variance Gamma process with an additional parameter Y controlling the fine structure of jumps:

Y < 0: finite activity (finitely many jumps in any interval)
0 <= Y < 1: infinite activity, finite variation
1 <= Y < 2: infinite activity, infinite variation

\[\phi(u) = \exp\bigl(iu(r + \omega)T + C\,T\,\Gamma(-Y)\bigl[ (M - iu)^Y - M^Y + (G + iu)^Y - G^Y \bigr]\bigr)\]

where \(\omega = -C\,\Gamma(-Y)[(M-1)^Y - M^Y + (G+1)^Y - G^Y]\) is the convexity correction.

Parameters:

C (float) – Overall activity level (C > 0).
G (float) – Rate of exponential decay of the positive jump density (G > 0). Larger G = fewer large positive jumps.
M (float) – Rate of exponential decay of the negative jump density (M > 0). Larger M = fewer large negative jumps.
Y_param (float) – Fine structure parameter (Y < 2, Y != 0, Y != 1). See classification above.
rf (float) – Risk-free rate.
T (float) – Time to maturity.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) returning the characteristic function.

Example

>>> char_fn = cgmy_characteristic(1.0, 5.0, 10.0, 0.5,
...                               0.05, 1.0, spot=100)
>>> import numpy as np
>>> np.abs(char_fn(np.array([0.0]))[0])
1.0

References

Carr, P., Geman, H., Madan, D.B. & Yor, M. (2002). The Fine Structure of Asset Returns: An Empirical Investigation. Journal of Business 75(2), 305-332.

fbsde_european(spot, payoff_fn, drift_fn, vol_fn, rf, T, n_steps=100, n_paths=10000, seed=None)[source]¶

Solve a forward-backward SDE for European option pricing.

The FBSDE system couples the asset dynamics (forward SDE) with the pricing equation (backward SDE):

\[\begin{split}\text{Forward:} \quad dX_t &= \mu(X_t)\,dt + \sigma(X_t)\,dW_t \\ \text{Backward:} \quad dY_t &= -f(t,X_t,Y_t,Z_t)\,dt + Z_t\,dW_t \\ \text{Terminal:} \quad Y_T &= g(X_T)\end{split}\]

Under risk-neutral pricing the driver is f = -r * Y (discounting). The forward SDE is discretised with Euler-Maruyama and the backward component is solved via least-squares regression at each time step (Longstaff-Schwartz style), projecting the continuation value onto polynomial basis functions of the forward process.

The Z process recovered from regression is proportional to sigma(X) * Delta, giving the hedge ratio directly from the BSDE without any finite-difference bumping.

When drift_fn = lambda x: r * x and vol_fn = lambda x: sigma * x with constant sigma, this reduces to the Black-Scholes model and the price converges to the analytical formula.

Parameters:

spot (float) – Initial value of the forward process X_0.
payoff_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Terminal condition g(X_T). Callable mapping an array of terminal values (shape (n_paths,)) to payoffs.
drift_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Drift coefficient mu(x) of the forward SDE. Callable mapping an array of current values to drift values.
vol_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Diffusion coefficient sigma(x) of the forward SDE. Callable mapping an array of current values to volatility values.
rf (float) – Risk-free rate (used as the BSDE driver: f = -r * Y).
T (float) – Time to maturity in years.
n_steps (int, default: 100) – Number of Euler-Maruyama time steps.
n_paths (int, default: 10000) – Number of Monte Carlo paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

price – estimated option price (mean of Y_0).
delta – estimated hedge ratio at t=0 (mean of Z_0 / vol(X_0)).
paths – forward process paths, shape (n_steps + 1, n_paths).
price_process – Y process, shape (n_steps + 1, n_paths).

Return type:

Example

>>> import numpy as np
>>> # Black-Scholes European call
>>> S0, K, r, sigma, T = 100.0, 100.0, 0.05, 0.2, 1.0
>>> payoff = lambda x: np.maximum(x - K, 0.0)
>>> drift = lambda x: r * x
>>> vol = lambda x: sigma * x
>>> result = fbsde_european(S0, payoff, drift, vol, r, T,
...                         n_steps=100, n_paths=50000, seed=42)
>>> 8.0 < result['price'] < 13.0
True

deep_bsde(dim, payoff_fn, drift_fn, vol_fn, rf, T, n_steps=50, n_paths=4096, n_epochs=200, lr=0.001, seed=None)[source]¶

Deep BSDE solver for high-dimensional derivative pricing.

Implements the algorithm of E, Han & Jentzen (2017) for solving BSDEs using deep neural networks. This is the state-of-the-art method for pricing derivatives on baskets of >3 assets, where PDE grid methods suffer from the curse of dimensionality.

The key idea: parameterise the initial value Y_0 and the Z process at each time step with neural networks (or simpler function approximators). Then minimise the terminal loss:

\[\mathcal{L} = \mathbb{E}\bigl[\bigl| Y_T^{\theta} - g(X_T) \bigr|^2\bigr]\]

where Y_T^theta is obtained by rolling forward the discretised BSDE with the learned Z networks.

When torch is available, uses a proper neural network with Adam optimiser. Otherwise falls back to a simplified linear approximation with scipy.optimize.minimize.

Use for problems with >3 dimensions where PDE methods fail due to the curse of dimensionality.

Parameters:

dim (int) – Number of underlying assets (spatial dimension).
payoff_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Terminal condition g(X_T). Callable mapping an array of shape (n_paths, dim) to payoffs of shape (n_paths,).
drift_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Drift mu(X). Callable mapping (n_paths, dim) to (n_paths, dim).
vol_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Diffusion sigma(X). Callable mapping (n_paths, dim) to (n_paths, dim). Assumed diagonal diffusion.
rf (float) – Risk-free rate.
T (float) – Time to maturity in years.
n_steps (int, default: 50) – Number of time steps.
n_paths (int, default: 4096) – Number of Monte Carlo paths per epoch.
n_epochs (int, default: 200) – Number of training epochs.
lr (float, default: 0.001) – Learning rate.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

price – estimated option price (Y_0).
delta – estimated hedge ratios at t=0, shape (dim,).
loss_history – training loss at each epoch.

Return type:

Example

>>> import numpy as np
>>> # 1D Black-Scholes call (use fbsde_european for 1D; this is for demo)
>>> payoff = lambda x: np.maximum(np.mean(x, axis=1) - 100, 0.0)
>>> drift = lambda x: 0.05 * x
>>> vol = lambda x: 0.2 * x
>>> result = deep_bsde(1, payoff, drift, vol, 0.05, 1.0,
...                    n_steps=20, n_paths=512, n_epochs=50, seed=42)
>>> result['price'] > 0
True

References

E, W., Han, J., & Jentzen, A. (2017). Deep Learning-Based Numerical Methods for High-Dimensional Parabolic Partial Differential Equations and Backward Stochastic Differential Equations. Communications in Mathematics and Statistics, 5(4).

reflected_bsde(spot, payoff_fn, drift_fn, vol_fn, rf, T, n_steps=100, n_paths=10000, seed=None)[source]¶

Solve a reflected BSDE (RBSDE) for American option pricing.

American options require the price process Y_t to stay above the intrinsic value (the obstacle) at all times. This is captured by adding a non-decreasing “reflection” process K_t to the BSDE:

\[Y_t = g(X_T) + \int_t^T f(s, X_s, Y_s, Z_s)\,ds - \int_t^T Z_s\,dW_s + K_T - K_t\]

with the constraint \(Y_t \geq h(t, X_t)\) (obstacle) and \(K\) increases only when \(Y_t = h(t, X_t)\).

This implementation uses the penalisation method: at each backward step, the continuation value is lifted to at least the exercise value, effectively penalising paths that would violate the obstacle. This converges to the true reflected BSDE solution as the time grid is refined.

The resulting price is always >= the European option price, with the difference representing the early exercise premium.

Parameters:

spot (float) – Initial value of the forward process X_0.
payoff_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Payoff function g(x) used both as terminal condition and exercise (obstacle) value at each time step.
drift_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Drift coefficient mu(x) of the forward SDE.
vol_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Diffusion coefficient sigma(x) of the forward SDE.
rf (float) – Risk-free rate.
T (float) – Time to maturity in years.
n_steps (int, default: 100) – Number of time steps.
n_paths (int, default: 10000) – Number of Monte Carlo paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

price – estimated American option price.
exercise_boundary – array of exercise boundary estimates at each time step, shape (n_steps + 1,).
optimal_stopping_time – mean optimal stopping time.

Return type:

Example

>>> import numpy as np
>>> S0, K, r, sigma, T = 100.0, 100.0, 0.05, 0.2, 1.0
>>> payoff = lambda x: np.maximum(K - x, 0.0)  # put
>>> drift = lambda x: r * x
>>> vol = lambda x: sigma * x
>>> result = reflected_bsde(S0, payoff, drift, vol, r, T,
...                         n_steps=100, n_paths=20000, seed=42)
>>> result['price'] > 0
True

geometric_brownian_motion(S0, mu, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate paths of geometric Brownian motion (GBM).

GBM is the standard model for equity price dynamics and underlies the Black-Scholes framework. The SDE is:

\[dS_t = \mu\,S_t\,dt + \sigma\,S_t\,dW_t\]

The solution is log-normal: \(S_T = S_0 \exp\bigl((\mu - \sigma^2/2)T + \sigma W_T\bigr)\).

Use GBM for quick scenario generation, Monte Carlo option pricing, or as a baseline against more complex processes (Heston, jump diffusion).

Parameters:

S0 (float) – Initial price.
mu (float) – Drift rate (annualized). Under the risk-neutral measure, set mu = r (the risk-free rate).
sigma (float) – Volatility (annualized).
T (float) – Time horizon in years.
n_steps (int) – Number of time steps. 252 for daily resolution.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated price paths. Column 0 is the initial price S0.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = geometric_brownian_motion(100, 0.05, 0.2, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

heston: Stochastic volatility extension of GBM. jump_diffusion: GBM with Poisson-driven jumps.

heston(S0, v0, mu, kappa, theta, sigma_v, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Heston (1993) stochastic volatility model.

The Heston model captures the leverage effect (negative correlation between returns and volatility) and generates realistic volatility smiles. The coupled SDEs are:

\[\begin{split}dS_t &= \mu\,S_t\,dt + \sqrt{v_t}\,S_t\,dW_t^1 \\ dv_t &= \kappa(\theta - v_t)\,dt + \sigma_v\sqrt{v_t}\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

Uses the full truncation scheme (replace negative variance with zero before computing diffusion) to ensure non-negative variance in the discretisation.

Use Heston when you need to capture the volatility smile or when constant-volatility GBM is insufficient.

Parameters:

S0 (float) – Initial price.
v0 (float) – Initial variance (e.g., 0.04 for 20% vol).
mu (float) – Drift rate (annualized).
kappa (float) – Mean reversion speed of variance.
theta (float) – Long-run variance level.
sigma_v (float) – Volatility of variance (vol of vol).
rho (float) – Correlation between price and variance Brownian motions. Typically negative for equities (-0.7 to -0.3).
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

(price_paths, vol_paths), each of: shape (n_paths, n_steps + 1).

Return type:

tuple[ndarray[tuple[Any, ...], dtype[double]], ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> prices, vols = heston(100, 0.04, 0.05, 2.0, 0.04, 0.3, -0.7, 1.0, 252, 1000, seed=42)
>>> prices.shape
(1000, 253)

Notes

The Feller condition \(2\kappa\theta \geq \sigma_v^2\) ensures variance stays strictly positive in continuous time. Even when violated, the full truncation scheme keeps the discretisation non-negative.

See also

geometric_brownian_motion: Constant-vol baseline. heston_characteristic: Analytical characteristic function for

Fourier-based Heston option pricing.

References

Heston, S.L. (1993). A Closed-Form Solution for Options with Stochastic Volatility. Review of Financial Studies 6(2).

jump_diffusion(S0, mu, sigma, lam, jump_mean, jump_std, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Merton (1976) jump-diffusion model.

Extends GBM with random jumps driven by a compound Poisson process. This captures sudden large moves (crashes, earnings surprises) that GBM cannot reproduce:

\[\frac{dS}{S} = (\mu - \lambda k)\,dt + \sigma\,dW + J\,dN\]

where \(N\) is a Poisson process with intensity \(\lambda\), \(\ln(1+J) \sim N(\mu_J, \sigma_J^2)\), and \(k = E[e^J - 1]\) is the drift compensation.

Use jump diffusion when you need to model fat tails and sudden discontinuities in asset prices, or to capture the implied volatility smile steepening at short maturities.

Parameters:

S0 (float) – Initial price.
mu (float) – Drift rate (annualized, before jump compensation).
sigma (float) – Diffusion volatility (annualized).
lam (float) – Jump intensity (expected number of jumps per year).
jump_mean (float) – Mean of log-jump size. Negative values model downward crashes.
jump_std (float) – Standard deviation of log-jump size.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated price paths.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = jump_diffusion(100, 0.05, 0.2, 1.0, -0.1, 0.15, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

geometric_brownian_motion: Diffusion-only baseline. heston: Stochastic volatility without jumps.

References

Merton, R.C. (1976). Option Pricing When Underlying Stock Returns Are Discontinuous. Journal of Financial Economics 3.

ornstein_uhlenbeck(x0, theta, mu, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Ornstein-Uhlenbeck (OU) mean-reverting process.

The OU process is the canonical continuous-time mean-reverting model, widely used for interest rates, pairs-trading spreads, and volatility dynamics:

\[dX_t = \theta\,(\mu - X_t)\,dt + \sigma\,dW_t\]

This implementation uses the exact discretisation (not Euler) so it is accurate for any step size.

Use OU for modelling quantities that revert to a long-run level: interest rate spreads, log-volatility, or cointegrated pairs.

Parameters:

x0 (float) – Initial value.
theta (float) – Mean reversion speed. Higher values mean faster reversion. The half-life is \(\ln(2) / \theta\).
mu (float) – Long-run mean level.
sigma (float) – Volatility (diffusion coefficient).
T (float) – Time horizon.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated paths.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = ornstein_uhlenbeck(0.05, 5.0, 0.03, 0.01, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

cir_process: Mean-reverting process with non-negative constraint. simulate_vasicek: OU-based interest rate model with bond pricing.

cir_process(x0, kappa, theta, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Cox-Ingersoll-Ross (CIR) process.

The CIR process is the standard mean-reverting, non-negative model for interest rates and variance:

\[dX_t = \kappa(\theta - X_t)\,dt + \sigma\sqrt{X_t}\,dW_t\]

The \(\sqrt{X}\) diffusion term ensures the volatility decreases as the process approaches zero, preventing negative values when the Feller condition \(2\kappa\theta \geq \sigma^2\) holds.

This is a plain-array simulation. For analytics (bond prices, Feller diagnostics), use simulate_cir() instead.

Parameters:

x0 (float) – Initial value (must be positive).
kappa (float) – Mean reversion speed.
theta (float) – Long-run mean level.
sigma (float) – Volatility coefficient.
T (float) – Time horizon.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated paths. Values are kept non-negative via reflection.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = cir_process(0.04, 2.0, 0.04, 0.1, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

simulate_cir: CIR with Feller diagnostics and bond pricing. ornstein_uhlenbeck: Gaussian mean-reverting process (can go

negative).

simulate_sabr(f0, sigma0, alpha, beta, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the SABR stochastic volatility model.

The SABR model (Hagan et al., 2002) is the industry standard for modelling interest rate and FX volatility smiles:

\[\begin{split}dF_t &= \sigma_t\,F_t^{\beta}\,dW_t^1 \\ d\sigma_t &= \alpha\,\sigma_t\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

Key parameters:

\(\beta \in [0, 1]\) controls the backbone: 0 = normal, 1 = lognormal, in between = CEV.
\(\alpha\) is the vol-of-vol.
\(\rho\) controls the skew direction.

When to use SABR: Primarily for interest rate derivatives (swaptions, caps/floors) and FX options where the smile shape is well-characterised by the SABR formula.

Calibration guidance:

sigma0: Calibrate from the ATM implied volatility. For swaptions, this is the ATM normal or lognormal vol.
alpha: Controls the curvature of the smile. Higher alpha = more curvature. Typically 0.2 – 0.6.
beta: Often fixed by convention (beta=0.5 for rates, beta=1.0 for FX) and the other parameters calibrated.
rho: Controls the skew. rho < 0 gives a “normal” skew (downside puts more expensive). rho > 0 reverses it. Typically -0.5 to 0.0 for rates.

When to use Heston vs SABR:

SABR for rates/FX: the SABR formula gives a fast, analytical approximation for implied volatility that is widely used on trading desks.
Heston for equities: the Heston model better captures the mean-reverting nature of equity volatility.
Use simulation (this function) when you need full path dynamics, e.g., for path-dependent payoffs or for evaluating the approximation quality.

Parameters:

f0 (float) – Initial forward rate / price.
sigma0 (float) – Initial stochastic volatility level.
alpha (float) – Volatility of volatility (vol-of-vol).
beta (float) – CEV exponent (0 = normal, 1 = lognormal).
rho (float) – Correlation between forward and vol Brownian motions.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

forwards – simulated forward paths, shape (n_steps + 1, n_paths).
vols – simulated stochastic vol paths, shape (n_steps + 1, n_paths).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> # Simulate a 5% swaption with typical SABR params
>>> result = simulate_sabr(0.05, 0.3, 0.4, 0.5, -0.3, 1.0, 252, 1000, seed=42)
>>> result['forwards'].shape
(253, 1000)

References

Hagan, P.S., Kumar, D., Lesniewski, A.S. & Woodward, D.E. (2002). Managing Smile Risk. Wilmott Magazine, 84-108.

simulate_rough_bergomi(spot, xi, eta, H, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the rough Bergomi stochastic volatility model.

The rough Bergomi model (Bayer, Friz & Gatheral, 2016) is the state-of-the-art for capturing the rough nature of volatility. Empirical studies show that log-volatility behaves like a fractional Brownian motion (fBM) with Hurst parameter H ~ 0.1, much rougher than classical models (H = 0.5 for standard BM).

The model specifies:

\[\begin{split}dS_t &= \sqrt{v_t}\,S_t\,dW_t^1 \\ v_t &= \xi\,\mathcal{E}\bigl(\eta\,\tilde{W}_t^H\bigr)\end{split}\]

where \(\tilde{W}^H\) is a (Riemann-Liouville) fractional Brownian motion with Hurst parameter \(H < 0.5\), and \(\mathcal{E}\) is the Wick exponential \(\exp(x - \frac{1}{2}\text{Var}(x))\).

The fBM is simulated via the Cholesky method on the exact covariance matrix, which is accurate but O(n_steps^2) in memory.

When to use rough Bergomi: When you need to capture the term structure of at-the-money skew, which decays like \(T^{H-1/2}\) as observed in markets. Classical models (Heston, SABR) cannot reproduce this power-law behaviour.

Heston vs SABR vs rough Bergomi – when to use which:

Heston: Best for equity options when you need a tractable model with mean-reverting vol. Has analytical characteristic function for fast pricing. Cannot fit the short-maturity skew well.
SABR: Best for rates/FX desks. Has a fast closed-form implied vol approximation. Industry standard for swaptions.
Rough Bergomi: State-of-the-art for capturing the empirical power-law decay of ATM skew across maturities. Most realistic for equity index options but expensive to simulate (O(n_steps^2) for the Cholesky approach).

Calibration guidance:

H ~ 0.1: The empirical consensus for equity index volatility (Gatheral, Jaisson & Rosenbaum, 2018).
xi: Set to the current level of implied variance (e.g., 0.04 for 20% vol).
eta ~ 1.5-2.5: Controls the vol-of-vol. Calibrate to the overall level of the smile.
rho ~ -0.7 to -0.9: Leverage effect.

Parameters:

spot (float) – Initial spot price.
xi (float) – Forward variance level (flat forward variance curve). E.g., 0.04 for 20% annualized vol.
eta (float) – Volatility of volatility. Typical range 1.0 – 3.0.
H (float) – Hurst parameter of fractional Brownian motion. Must be in (0, 0.5) for the “rough” regime. Empirical estimates cluster around 0.1.
rho (float) – Correlation between spot and vol Brownian motions. Typically -0.7 to -0.9 for equity indices.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps. O(n_steps^2) memory for the Cholesky covariance matrix, so keep moderate (100-500).
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

prices – simulated asset price paths, shape (n_steps + 1, n_paths).
variances – simulated instantaneous variance paths, shape (n_steps + 1, n_paths).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> # Simulate SPX-like dynamics with rough vol
>>> result = simulate_rough_bergomi(
...     spot=4000, xi=0.04, eta=1.9, H=0.1, rho=-0.7,
...     T=0.5, n_steps=100, n_paths=1000, seed=42)
>>> result['prices'].shape
(101, 1000)

References

Bayer, C., Friz, P. & Gatheral, J. (2016). Pricing Under Rough Volatility. Quantitative Finance 16(6), 887-904.

Gatheral, J., Jaisson, T. & Rosenbaum, M. (2018). Volatility is Rough. Quantitative Finance 18(6), 933-949.

simulate_3_2_model(spot, v0, kappa, theta, epsilon, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the 3/2 stochastic volatility model.

The 3/2 model features vol-of-vol that increases with the variance level, unlike Heston where vol-of-vol is constant. This makes it better suited for markets where volatility becomes extremely volatile during high-vol regimes (e.g., VIX dynamics):

\[\begin{split}dS_t &= \sqrt{v_t}\,S_t\,dW_t^1 \\ dv_t &= \kappa\,v_t(\theta - v_t)\,dt + \varepsilon\,v_t^{3/2}\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

The key feature is the \(v^{3/2}\) diffusion term: when variance is high, the variance process itself becomes much more volatile, capturing the empirical observation that vol-of-vol spikes during market crises.

When to use the 3/2 model: For VIX derivatives, or when the Heston model under-estimates the vol-of-vol in high-vol regimes. Also useful when calibrating to options on variance swaps.

Parameters:

spot (float) – Initial asset price.
v0 (float) – Initial variance.
kappa (float) – Mean reversion speed (in the V dynamics).
theta (float) – Long-run variance target.
epsilon (float) – Vol-of-vol parameter (coefficient on V^{3/2}).
rho (float) – Correlation between asset and variance Brownian motions.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

prices – simulated asset price paths, shape (n_steps + 1, n_paths).
variances – simulated variance paths, shape (n_steps + 1, n_paths).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> result = simulate_3_2_model(100.0, 0.04, 2.0, 0.04, 0.5, -0.7,
...                             1.0, 252, 1000, seed=42)
>>> result['prices'].shape
(253, 1000)
>>> result['variances'].shape
(253, 1000)

simulate_cir(r0, kappa, theta, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Cox-Ingersoll-Ross short-rate model with analytics.

The CIR model is the standard mean-reverting, non-negative process for interest rates:

\[dr_t = \kappa(\theta - r_t)\,dt + \sigma\sqrt{r_t}\,dW_t\]

The Feller condition \(2\kappa\theta \geq \sigma^2\) ensures that the process stays strictly positive. When violated, zero is accessible but reflecting.

Unlike cir_process() (which returns a plain array), this function also computes: - The Feller condition diagnostic - Analytical zero-coupon bond prices P(0, T) - The mean reversion check

When to use CIR: For modelling interest rates, credit intensities, or stochastic variance (it is the variance process in Heston). The guaranteed non-negativity is crucial for these applications.

Parameters:

r0 (float) – Initial short rate (must be positive).
kappa (float) – Mean reversion speed.
theta (float) – Long-run mean level.
sigma (float) – Volatility coefficient.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

paths – simulated rate paths, shape (n_steps + 1, n_paths).
params – dictionary of model parameters including: - kappa, theta, sigma – model parameters. - feller_satisfied – whether 2*kappa*theta >= sigma^2. - feller_ratio – 2*kappa*theta / sigma^2 (>= 1 means satisfied).

Return type:

Example

>>> result = simulate_cir(0.05, 0.5, 0.04, 0.1, 10.0, 252, 1000, seed=42)
>>> result['paths'].shape
(253, 1000)
>>> result['params']['feller_satisfied']
True

simulate_vasicek(r0, kappa, theta, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Vasicek interest rate model with bond pricing.

The Vasicek (1977) model is the simplest mean-reverting Gaussian model for the short rate:

\[dr_t = \kappa(\theta - r_t)\,dt + \sigma\,dW_t\]

Being Gaussian, the short rate can go negative – this was historically seen as a drawback but is now relevant in the era of negative rates.

The model admits closed-form solutions for zero-coupon bond prices:

\[P(0, T) = A(T)\,\exp(-B(T)\,r_0)\]

where:

\[\begin{split}B(T) &= \frac{1 - e^{-\kappa T}}{\kappa} \\ A(T) &= \exp\left[ \left(\theta - \frac{\sigma^2}{2\kappa^2}\right)(B(T) - T) - \frac{\sigma^2}{4\kappa}B(T)^2 \right]\end{split}\]

When to use Vasicek: Quick-and-dirty interest rate modelling, or when negative rates are possible/desired. For non-negative rates, use CIR instead.

Parameters:

r0 (float) – Initial short rate.
kappa (float) – Mean reversion speed.
theta (float) – Long-run mean level.
sigma (float) – Volatility.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

paths – simulated rate paths, shape (n_steps + 1, n_paths).
bond_prices – analytical zero-coupon bond prices P(0, t) at each time step, shape (n_steps + 1,).
yield_curve – continuously compounded yields y(0, t) = -ln(P(0,t))/t at each time step (excluding t=0), shape (n_steps,).

Return type:

Example

>>> result = simulate_vasicek(0.05, 0.5, 0.04, 0.01, 10.0, 100, 1000, seed=42)
>>> result['paths'].shape
(101, 1000)
>>> len(result['bond_prices'])
101
>>> len(result['yield_curve'])
100

References

Vasicek, O. (1977). An Equilibrium Characterisation of the Term Structure. Journal of Financial Economics 5(2), 177-188.

quantlib_bond(face, coupon, maturity, yield_curve, settlement=None, frequency=2)[source]¶

Price a fixed-rate bond using QuantLib.

Parameters:

face (float) – Face (par) value of the bond.
coupon (float) – Annual coupon rate (e.g. 0.05 for 5 %).
maturity (date) – Bond maturity date.
yield_curve (list[tuple[date, float]]) – Pairs of (date, zero_rate) defining the yield curve.
settlement (date | None, default: None) – Settlement date. Defaults to today when None.
frequency (int, default: 2) – Coupon payment frequency per year (2 = semi-annual).

Returns:

Dictionary containing:

clean_price – clean price of the bond.
dirty_price – dirty (full) price of the bond.
ytm – yield to maturity.
duration – Macaulay duration.
convexity – convexity.

Return type:

quantlib_option(spot, strike, vol, rf, maturity, option_type='call')[source]¶

Price a European option using QuantLib’s analytic engine.

Parameters:

spot (float) – Current spot price of the underlying.
strike (float) – Strike price.
vol (float) – Annualised volatility (e.g. 0.20 for 20 %).
rf (float) – Risk-free interest rate (annual, continuously compounded).
maturity (float) – Time to maturity in years.
option_type (str, default: 'call') – 'call' or 'put'.

Returns:

Dictionary containing:

price – option price.
delta – option delta.
gamma – option gamma.
theta – option theta (per year).
vega – option vega (per 1 % vol change).
rho – option rho.

Return type:

quantlib_yield_curve(dates, rates)[source]¶

Construct a yield curve using QuantLib.

Parameters:

dates (list[date]) – Curve node dates (the first date is the reference date).
rates (list[float]) – Zero rates corresponding to each date.

Returns:

Dictionary containing:

reference_date – the curve reference date.
max_date – the latest date on the curve.
discount_factors – list of discount factors at each node.
forward_rates – list of instantaneous forward rates at each node.
curve – the raw ql.ZeroCurve object.

Return type:

financepy_option(spot, strike, vol, rf, maturity)[source]¶

Price a European equity option using FinancePy.

Parameters:

spot (float) – Current spot price.
strike (float) – Strike price.
vol (float) – Annualised volatility.
rf (float) – Risk-free rate (annual).
maturity (float) – Time to maturity in years.

Returns:

Dictionary containing:

call_price – European call price.
put_price – European put price.
call_delta – call delta.
put_delta – put delta.

Return type:

rateslib_swap(notional, fixed_rate, float_spread, maturity)[source]¶

Price a plain vanilla interest rate swap using rateslib.

Parameters:

notional (float) – Notional principal amount.
fixed_rate (float) – Fixed leg coupon rate (e.g. 0.03 for 3 %).
float_spread (float) – Spread over the floating index in basis points.
maturity (str) – Tenor string (e.g. '5Y', '10Y').

Returns:

Dictionary containing:

npv – net present value of the swap.
fixed_rate – the fixed rate used.
float_spread – the floating spread used (bps).
maturity – the tenor string.
notional – the notional amount.

Return type:

vollib_implied_vol(price, spot, strike, rf, maturity, option_type='call')[source]¶

Compute implied volatility using py-vollib.

Parameters:

price (float) – Observed option market price.
spot (float) – Current spot price of the underlying.
strike (float) – Strike price.
rf (float) – Risk-free rate (annual, continuously compounded).
maturity (float) – Time to maturity in years.
option_type (str, default: 'call') – 'call' or 'put'.

Returns:

Dictionary containing:

implied_vol – implied volatility.
option_type – the option type used.
price – the input option price.

Return type:

sdeint_solve(drift_fn, diffusion_fn, y0, tspan, dt=0.01, method='euler')[source]¶

Solve a stochastic differential equation using sdeint.

For SDEs of the form:

dY = f(Y, t) dt + g(Y, t) dW

where f is the drift and g is the diffusion coefficient.

sdeint provides Ito-SDE solvers with strong convergence guarantees. Use Euler-Maruyama for quick simulations and Milstein when you need higher-order accuracy (order-1.0 strong convergence vs order-0.5 for Euler).

Parameters:

drift_fn (callable) – f(Y, t) -> array. Drift coefficient. Must accept the current state Y and time t and return an array of the same shape as Y.
diffusion_fn (callable) – g(Y, t) -> array. Diffusion coefficient. For scalar SDEs, return a 1-D array of shape (1,). For vector SDEs, return a 2-D array of shape (d, m) where d is the state dimension and m is the number of Wiener processes.
y0 (float | ndarray) – Initial condition. Scalar or 1-D array.
tspan (tuple[float, float]) – (t0, t1) time interval.
dt (float, default: 0.01) – Time step size.
method (str, default: 'euler') –
Integration method:
- 'euler' – Euler-Maruyama (strong order 0.5).
- 'sri2' – Roessler SRI2 Stochastic Runge-Kutta (strong order 1.0). Suitable for arbitrary noise.

Returns:

Dictionary containing:

times – 1-D array of time points.
paths – solution array of shape (n_steps, d) where d is the state dimension.
final_values – state at the terminal time.

Return type:

Example

>>> import numpy as np
>>> from wraquant.price.integrations import sdeint_solve
>>> # Geometric Brownian Motion: dS = mu*S*dt + sigma*S*dW
>>> result = sdeint_solve(
...     drift_fn=lambda y, t: 0.05 * y,
...     diffusion_fn=lambda y, t: 0.2 * y,
...     y0=100.0,
...     tspan=(0.0, 1.0),
...     dt=0.01,
... )
>>> result["paths"].shape[0]  # number of time steps
101

Notes

Reference: Kloeden & Platen (1992). Numerical Solution of Stochastic Differential Equations. Springer.

See also

wraquant.price.stochastic: Built-in process simulators (GBM, Heston, SABR, etc.) that don’t require sdeint.

Options Pricing¶

Options pricing models: Black-Scholes, binomial tree, and Monte Carlo.

Provides pure numpy/scipy implementations of standard option pricing models for European and American options.

black_scholes(S, K, T, r, sigma, option_type='call')[source]¶

Price a European option using the Black-Scholes closed-form formula.

For a call:

\[C = S\,\Phi(d_1) - K\,e^{-rT}\,\Phi(d_2)\]

where \(d_1 = \frac{\ln(S/K) + (r + \sigma^2/2)T}{\sigma\sqrt{T}}\) and \(d_2 = d_1 - \sigma\sqrt{T}\).

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years. At T <= 0 the intrinsic value is returned.
r (float) – Risk-free interest rate (annualized, continuously compounded).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

The Black-Scholes option price. Always non-negative.: For deep out-of-the-money options the price approaches zero; for deep in-the-money options it approaches the intrinsic value discounted at the risk-free rate.

Return type:

Example

>>> black_scholes(100, 100, 1.0, 0.05, 0.2)
10.450...
>>> black_scholes(100, 100, 1.0, 0.05, 0.2, 'put')
5.573...

Notes

See also

binomial_tree: Lattice method supporting American exercise. monte_carlo_option: Simulation-based pricing. implied_volatility: Invert this formula to recover vol from price.

References

Black, F. & Scholes, M. (1973). The Pricing of Options and Corporate Liabilities. Journal of Political Economy 81(3).

binomial_tree(S, K, T, r, sigma, n_steps, option_type='call', style='european')[source]¶

Price an option using the Cox-Ross-Rubinstein (CRR) binomial tree.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized, continuously compounded).
sigma (float) – Volatility of the underlying (annualized).
n_steps (int) – Number of time steps in the tree. 200–500 steps give good accuracy for most applications.
option_type (str | OptionType, default: 'call') – 'call' or 'put'.
style (str | OptionStyle, default: 'european') – 'european' or 'american'.

Returns:

The binomial tree option price. For European options: this converges to the Black-Scholes price; for American options the price is at least as large as the European price due to the early exercise premium.

Return type:

Example

>>> binomial_tree(100, 100, 1.0, 0.05, 0.2, 200)
10.4...
>>> binomial_tree(100, 100, 1.0, 0.05, 0.2, 200, 'put', 'american')
6.08...

Notes

Convergence is O(1/n_steps). For faster convergence consider Richardson extrapolation (average the n and 2n step prices).

See also

black_scholes: Analytical European option pricing. monte_carlo_option: Simulation-based pricing.

References

Cox, J.C., Ross, S.A. & Rubinstein, M. (1979). Option Pricing: A Simplified Approach. Journal of Financial Economics 7(3).

monte_carlo_option(S, K, T, r, sigma, n_sims, n_steps, option_type='call', seed=None)[source]¶

Price a European option using Monte Carlo simulation under GBM.

The standard error of the estimate decreases as \(1/\sqrt{n\_sims}\), so 100 000 paths typically give 2–3 significant digits of accuracy.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized, continuously compounded).
sigma (float) – Volatility of the underlying (annualized).
n_sims (int) – Number of simulation paths. More paths reduce variance but increase computation time.
n_steps (int) – Number of time steps per path. For European vanilla options a single step suffices, but path-dependent payoffs require finer discretisation.
option_type (str | OptionType, default: 'call') – 'call' or 'put'.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

The Monte Carlo estimate of the option price. This: is an unbiased estimator that converges to the Black-Scholes price for European vanilla options under GBM.

Return type:

Example

>>> monte_carlo_option(100, 100, 1.0, 0.05, 0.2, 100000, 252, seed=42)
10.4...

Notes

For variance reduction consider antithetic variates or control variates (using the Black-Scholes price as a control).

See also

black_scholes: Closed-form European pricing (faster, exact). binomial_tree: Lattice method for American options.

Greeks¶

Option Greeks via analytical Black-Scholes formulas.

Provides closed-form solutions for delta, gamma, theta, vega, and rho for European options under the Black-Scholes model.

delta(S, K, T, r, sigma, option_type='call')[source]¶

Compute the Black-Scholes delta (sensitivity to underlying price).

\[\Delta_{\text{call}} = \Phi(d_1), \quad \Delta_{\text{put}} = \Phi(d_1) - 1\]

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Delta value. Call delta is in [0, 1]; put delta: is in [-1, 0]. At-the-money forward options have delta near 0.5 (call) or -0.5 (put).

Return type:

Example

>>> delta(100, 100, 1.0, 0.05, 0.2, 'call')
0.636...
>>> delta(100, 100, 1.0, 0.05, 0.2, 'put')
-0.363...

See also

gamma: Rate of change of delta. all_greeks: Compute all Greeks in a single call.

gamma(S, K, T, r, sigma)[source]¶

Compute the Black-Scholes gamma (rate of change of delta).

\[\Gamma = \frac{\phi(d_1)}{S\,\sigma\,\sqrt{T}}\]

Gamma is highest for at-the-money options near expiry.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).

Returns:

Gamma value (always non-negative). Expressed in: price units per unit squared move in the underlying.

Return type:

Example

>>> gamma(100, 100, 1.0, 0.05, 0.2)
0.018...

See also

delta: First-order sensitivity to the underlying. all_greeks: Compute all Greeks in a single call.

theta(S, K, T, r, sigma, option_type='call')[source]¶

Compute the Black-Scholes theta (time decay per year).

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Theta value per year. Typically negative for long: options. For daily theta, divide by 252 (or 365 for calendar days).

Return type:

Example

>>> theta(100, 100, 1.0, 0.05, 0.2, 'call')
-6.41...

See also

vega: Sensitivity to volatility. all_greeks: Compute all Greeks in a single call.

vega(S, K, T, r, sigma)[source]¶

Compute the Black-Scholes vega (sensitivity to volatility).

\[\mathcal{V} = S\,\phi(d_1)\,\sqrt{T}\]

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).

Returns:

Vega per unit change in sigma (always non-negative).: Multiply by 0.01 for the price change per 1% vol move.

Return type:

Example

>>> vega(100, 100, 1.0, 0.05, 0.2)
37.52...

See also

implied_volatility: Invert the BS formula using vega as the: Newton-Raphson derivative.

all_greeks: Compute all Greeks in a single call.

rho(S, K, T, r, sigma, option_type='call')[source]¶

Compute the Black-Scholes rho (sensitivity to interest rate).

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Rho per unit change in r. Positive for calls,: negative for puts.

Return type:

Example

>>> rho(100, 100, 1.0, 0.05, 0.2, 'call')
53.23...

See also

all_greeks: Compute all Greeks in a single call.

all_greeks(S, K, T, r, sigma, option_type='call')[source]¶

Compute all Black-Scholes Greeks in a single efficient call.

Parameters:

S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
sigma (float) – Volatility of the underlying (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Dictionary with keys 'delta',: 'gamma', 'theta', 'vega', 'rho'. See the individual Greek functions for interpretation of each value.

Return type:

Example

>>> greeks = all_greeks(100, 100, 1.0, 0.05, 0.2, 'call')
>>> sorted(greeks.keys())
['delta', 'gamma', 'rho', 'theta', 'vega']

See also

delta, gamma, theta, vega, rho: Individual Greek computations.

Volatility Surfaces¶

Implied volatility and volatility surface construction.

Provides Newton-Raphson implied volatility solver and tools for building volatility smiles and surfaces from market data.

implied_volatility(market_price, S, K, T, r, option_type='call', tol=1e-08, max_iter=200, initial_guess=0.3)[source]¶

Compute implied volatility using Newton-Raphson iteration.

When vega is too small (deep in/out of the money), the solver automatically falls back to a robust bisection method.

Parameters:

market_price (float) – Observed market price of the option.
S (float) – Current underlying price.
K (float) – Strike price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.
tol (float, default: 1e-08) – Convergence tolerance on the price difference.
max_iter (int, default: 200) – Maximum number of Newton iterations.
initial_guess (float, default: 0.3) – Starting volatility estimate. 0.3 (30%) is a reasonable default for most equity options.

Returns:

Implied volatility (annualized). For example, 0.20: corresponds to 20% annualized volatility.

Return type:

Raises:

PricingError – If the solver does not converge within max_iter iterations.

Example

>>> price = black_scholes(100, 100, 1.0, 0.05, 0.2)
>>> implied_volatility(price, 100, 100, 1.0, 0.05)
0.2000...

See also

vol_smile: Compute implied vols across multiple strikes. vol_surface: Build a full implied volatility surface.

vol_smile(strikes, market_prices, S, T, r, option_type='call')[source]¶

Compute a volatility smile from market prices at multiple strikes.

Parameters:

strikes (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of strike prices.
market_prices (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Corresponding market option prices (same length as strikes).
S (float) – Current underlying price.
T (float) – Time to expiration in years.
r (float) – Risk-free interest rate (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Dictionary with keys:

'strikes' – 1D array of strike prices.
'implied_vols' – 1D array of implied volatilities corresponding to each strike.

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Raises:

ValueError – If strikes and market_prices differ in length.
PricingError – If any individual implied vol fails to converge.

Example

>>> strikes = [90, 95, 100, 105, 110]
>>> prices = [12.5, 9.0, 6.0, 3.8, 2.1]
>>> result = vol_smile(strikes, prices, 100, 0.5, 0.05)
>>> len(result['implied_vols'])
5

See also

implied_volatility: Single-strike implied vol solver. vol_surface: Extend the smile across multiple expiries.

vol_surface(strikes, expiries, market_prices, S, r, option_type='call')[source]¶

Construct a volatility surface from market prices across strikes and expiries.

Parameters:

strikes (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of strike prices (columns of the surface).
expiries (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of times to expiry in years (rows of the surface).
market_prices (Sequence[Sequence[float]] | ndarray[tuple[Any, ...], dtype[floating]]) – 2D array of option prices, shape (len(expiries), len(strikes)).
S (float) – Current underlying price.
r (float) – Risk-free interest rate (annualized).
option_type (str | OptionType, default: 'call') – 'call' or 'put'.

Returns:

Dictionary with keys:

'strikes' – 1D array of strikes.
'expiries' – 1D array of expiries.
'implied_vols' – 2D array of implied vols, shape (len(expiries), len(strikes)).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Raises:

ValueError – If the shape of market_prices does not match (len(expiries), len(strikes)).

Example

>>> strikes = [95, 100, 105]
>>> expiries = [0.25, 0.5]
>>> prices = [[6.0, 3.5, 1.5], [8.0, 5.5, 3.5]]
>>> surface = vol_surface(strikes, expiries, prices, 100, 0.05)
>>> surface['implied_vols'].shape
(2, 3)

See also

vol_smile: Single-expiry volatility smile. implied_volatility: Single-point implied vol solver.

Fixed Income¶

Fixed income pricing: bonds, yields, duration, and convexity.

Provides present-value based bond pricing, yield-to-maturity solvers, and duration/convexity risk measures.

bond_price(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the price of a fixed-rate bond as the PV of its cash flows.

\[P = \sum_{t=1}^{n} \frac{C}{(1+y)^t} + \frac{F}{(1+y)^n}\]

where \(C = F \cdot c / f\) is the periodic coupon, \(y = \text{ytm} / f\), \(n\) is the number of periods, and \(f\) is the payment frequency.

Parameters:

face_value (float) – Par/face value of the bond (e.g., 1000).
coupon_rate (float) – Annual coupon rate (e.g., 0.05 for 5%).
ytm (float) – Yield to maturity (annualized, e.g., 0.05 for 5%).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Present value (dirty price) of the bond. When: ytm < coupon_rate the bond trades at a premium (price > face); when ytm > coupon_rate it trades at a discount.

Return type:

Example

>>> bond_price(1000, 0.05, 0.05, 10, 2)
1000.0
>>> bond_price(1000, 0.05, 0.06, 10, 2)
925.6...

See also

bond_yield: Solve for the yield given price. duration: Macaulay duration for interest-rate sensitivity.

bond_yield(price, face_value, coupon_rate, periods, freq=2, tol=1e-10, max_iter=200)[source]¶

Compute the yield to maturity (YTM) of a bond via Newton’s method.

Parameters:

price (float) – Market price of the bond.
face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).
tol (float, default: 1e-10) – Convergence tolerance on the price residual.
max_iter (int, default: 200) – Maximum Newton-Raphson iterations.

Returns:

Annualized yield to maturity. For a par bond: (price == face_value) the YTM equals the coupon rate.

Return type:

Raises:

ValueError – If the solver does not converge within max_iter iterations.

Example

>>> bond_yield(950, 1000, 0.05, 10, 2)
0.059...

See also

bond_price: Price a bond given a yield. zero_rate: Extract a zero rate from a zero-coupon bond.

duration(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the Macaulay duration of a fixed-rate bond.

\[D = \frac{1}{P} \sum_{t=1}^{n} \frac{t}{f} \cdot \frac{CF_t}{(1 + y)^t}\]

Parameters:

face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
ytm (float) – Yield to maturity (annualized).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Macaulay duration in years. Zero-coupon bonds have: duration equal to their maturity; coupon bonds always have duration less than maturity.

Return type:

Example

>>> duration(1000, 0.05, 0.05, 10, 2)
4.08...

See also

modified_duration: Duration divided by (1 + y/freq). convexity: Second-order interest-rate sensitivity.

modified_duration(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the modified duration of a fixed-rate bond.

Modified duration is the more practically useful sensitivity measure: the percentage price change for a one-unit change in yield.

\[D_{\text{mod}} = \frac{D_{\text{mac}}}{1 + y/f}\]

A modified duration of 4 means a 100 bp (1%) yield increase causes approximately a 4% price decline.

Parameters:

face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
ytm (float) – Yield to maturity (annualized).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Modified duration in years.

Return type:

Example

>>> modified_duration(1000, 0.05, 0.05, 10, 2)
3.98...

See also

duration: Macaulay duration. convexity: Second-order correction for large yield changes.

convexity(face_value, coupon_rate, ytm, periods, freq=2)[source]¶

Compute the convexity of a fixed-rate bond.

The price change for a yield shift \(\Delta y\) is approximately:

\[\frac{\Delta P}{P} \approx -D_{\text{mod}} \cdot \Delta y + \tfrac{1}{2}\,C \cdot (\Delta y)^2\]

where \(C\) is convexity. Positive convexity means the bond gains more from a yield decrease than it loses from an equal yield increase.

Parameters:

face_value (float) – Par/face value of the bond.
coupon_rate (float) – Annual coupon rate.
ytm (float) – Yield to maturity (annualized).
periods (int) – Total number of coupon periods remaining.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Convexity measure (in years squared, scaled by: freq^2).

Return type:

Example

>>> convexity(1000, 0.05, 0.05, 10, 2)
19.4...

See also

duration: First-order interest-rate sensitivity. modified_duration: Percentage price change per unit yield change.

zero_rate(price, face_value, periods)[source]¶

Compute the zero-coupon (spot) rate from a zero-coupon bond price.

Parameters:

price (float) – Current market price of the zero-coupon bond.
face_value (float) – Par/face value of the bond.
periods (float) – Time to maturity in years.

Returns:

Annualized continuously compounded zero rate.

Return type:

Raises:

ValueError – If price or periods is not positive.

Example

>>> zero_rate(950, 1000, 1.0)
0.051...

See also

bootstrap_zero_curve: Bootstrap a full zero curve from par rates. discount_factor: Convert a zero rate to a discount factor.

Yield Curves¶

Yield curve construction and interpolation.

Provides bootstrapping, interpolation, forward rate calculation, and discount factor utilities for yield curve analysis.

bootstrap_zero_curve(maturities, par_rates, freq=2)[source]¶

Bootstrap zero (spot) rates from par coupon rates.

The resulting zero rates can be used to discount arbitrary cash flows, compute forward rates, or price swaps and bonds.

Parameters:

maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of maturities in years, evenly spaced at 1/freq (e.g., 0.5, 1.0, 1.5 for semiannual).
par_rates (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of par coupon rates (annualized) for each maturity.
freq (int, default: 2) – Coupon payments per year (default 2 for semiannual).

Returns:

Array of continuously compounded zero rates: corresponding to each maturity.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Raises:

ValueError – If maturities and par_rates differ in length.

Example

>>> mats = [0.5, 1.0, 1.5, 2.0]
>>> pars = [0.04, 0.042, 0.044, 0.046]
>>> zeros = bootstrap_zero_curve(mats, pars, freq=2)
>>> len(zeros)
4

See also

forward_rate: Derive forward rates from zero rates. discount_factor: Convert a zero rate to a discount factor.

interpolate_curve(maturities, rates, target_maturities, method='cubic')[source]¶

Interpolate a yield curve at target maturities.

Parameters:

maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Known maturities (sorted ascending).
rates (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Known rates at those maturities.
target_maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Maturities at which to interpolate.
method (str, default: 'cubic') –
Interpolation method:
- 'linear' – piecewise linear (fast, may have kinks).
- 'cubic' – natural cubic spline (smooth, default).
- 'flat_forward' – constant forward rate between knots (no-arbitrage, standard in fixed income).

Returns:

Interpolated rates at target maturities.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Raises:

ValueError – If method is not recognised.

Example

>>> mats = [0.5, 1.0, 2.0, 5.0]
>>> rates = [0.03, 0.035, 0.04, 0.045]
>>> interpolate_curve(mats, rates, [0.75, 1.5, 3.0])
array([...])

See also

bootstrap_zero_curve: Build the curve from par rates. forward_rate: Extract a forward rate from the curve.

forward_rate(zero_rates, maturities, t1, t2)[source]¶

Compute the forward rate between two future dates.

The forward rate is the rate that can be locked in today for borrowing or lending between two future dates. Under no-arbitrage:

\[f(t_1, t_2) = \frac{r_2\,t_2 - r_1\,t_1}{t_2 - t_1}\]

Uses cubic spline interpolation to obtain zero rates at t1 and t2 from the provided curve.

Parameters:

zero_rates (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Array of continuously compounded zero rates.
maturities (Sequence[float] | ndarray[tuple[Any, ...], dtype[floating]]) – Corresponding maturities.
t1 (float) – Start of the forward period (years).
t2 (float) – End of the forward period (years).

Returns:

Continuously compounded forward rate between t1: and t2.

Return type:

Raises:

ValueError – If t2 <= t1.

Example

>>> forward_rate([0.04, 0.045], [1.0, 2.0], 1.0, 2.0)
0.05

See also

bootstrap_zero_curve: Build the zero curve from par rates. discount_factor: Convert a rate to a discount factor.

discount_factor(rate, maturity)[source]¶

Compute the discount factor from a continuously compounded rate.

The discount factor is the present value of one unit of currency received at a future date:

\[DF = e^{-r \cdot T}\]

Discount factors are the fundamental building blocks for pricing any fixed-income instrument.

Parameters:

rate (float) – Continuously compounded interest rate.
maturity (float) – Time to maturity in years.

Returns:

Discount factor in (0, 1] for positive rates.

Return type:

Example

>>> discount_factor(0.05, 1.0)
0.9512...
>>> discount_factor(0.0, 5.0)
1.0

See also

zero_rate: Extract a rate from a zero-coupon bond price. bootstrap_zero_curve: Build a full discount factor curve.

Stochastic Processes¶

Stochastic process simulation.

Pure numpy implementations of common stochastic processes used in quantitative finance: GBM, Heston, Merton jump-diffusion, Ornstein-Uhlenbeck, Cox-Ingersoll-Ross, SABR, rough Bergomi, 3/2 model, and Vasicek.

geometric_brownian_motion(S0, mu, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate paths of geometric Brownian motion (GBM).

GBM is the standard model for equity price dynamics and underlies the Black-Scholes framework. The SDE is:

\[dS_t = \mu\,S_t\,dt + \sigma\,S_t\,dW_t\]

The solution is log-normal: \(S_T = S_0 \exp\bigl((\mu - \sigma^2/2)T + \sigma W_T\bigr)\).

Use GBM for quick scenario generation, Monte Carlo option pricing, or as a baseline against more complex processes (Heston, jump diffusion).

Parameters:

S0 (float) – Initial price.
mu (float) – Drift rate (annualized). Under the risk-neutral measure, set mu = r (the risk-free rate).
sigma (float) – Volatility (annualized).
T (float) – Time horizon in years.
n_steps (int) – Number of time steps. 252 for daily resolution.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated price paths. Column 0 is the initial price S0.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = geometric_brownian_motion(100, 0.05, 0.2, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

heston: Stochastic volatility extension of GBM. jump_diffusion: GBM with Poisson-driven jumps.

heston(S0, v0, mu, kappa, theta, sigma_v, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Heston (1993) stochastic volatility model.

The Heston model captures the leverage effect (negative correlation between returns and volatility) and generates realistic volatility smiles. The coupled SDEs are:

\[\begin{split}dS_t &= \mu\,S_t\,dt + \sqrt{v_t}\,S_t\,dW_t^1 \\ dv_t &= \kappa(\theta - v_t)\,dt + \sigma_v\sqrt{v_t}\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

Uses the full truncation scheme (replace negative variance with zero before computing diffusion) to ensure non-negative variance in the discretisation.

Use Heston when you need to capture the volatility smile or when constant-volatility GBM is insufficient.

Parameters:

S0 (float) – Initial price.
v0 (float) – Initial variance (e.g., 0.04 for 20% vol).
mu (float) – Drift rate (annualized).
kappa (float) – Mean reversion speed of variance.
theta (float) – Long-run variance level.
sigma_v (float) – Volatility of variance (vol of vol).
rho (float) – Correlation between price and variance Brownian motions. Typically negative for equities (-0.7 to -0.3).
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

(price_paths, vol_paths), each of: shape (n_paths, n_steps + 1).

Return type:

tuple[ndarray[tuple[Any, ...], dtype[double]], ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> prices, vols = heston(100, 0.04, 0.05, 2.0, 0.04, 0.3, -0.7, 1.0, 252, 1000, seed=42)
>>> prices.shape
(1000, 253)

Notes

See also

geometric_brownian_motion: Constant-vol baseline. heston_characteristic: Analytical characteristic function for

Fourier-based Heston option pricing.

References

Heston, S.L. (1993). A Closed-Form Solution for Options with Stochastic Volatility. Review of Financial Studies 6(2).

jump_diffusion(S0, mu, sigma, lam, jump_mean, jump_std, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Merton (1976) jump-diffusion model.

Extends GBM with random jumps driven by a compound Poisson process. This captures sudden large moves (crashes, earnings surprises) that GBM cannot reproduce:

\[\frac{dS}{S} = (\mu - \lambda k)\,dt + \sigma\,dW + J\,dN\]

where \(N\) is a Poisson process with intensity \(\lambda\), \(\ln(1+J) \sim N(\mu_J, \sigma_J^2)\), and \(k = E[e^J - 1]\) is the drift compensation.

Use jump diffusion when you need to model fat tails and sudden discontinuities in asset prices, or to capture the implied volatility smile steepening at short maturities.

Parameters:

S0 (float) – Initial price.
mu (float) – Drift rate (annualized, before jump compensation).
sigma (float) – Diffusion volatility (annualized).
lam (float) – Jump intensity (expected number of jumps per year).
jump_mean (float) – Mean of log-jump size. Negative values model downward crashes.
jump_std (float) – Standard deviation of log-jump size.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated price paths.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = jump_diffusion(100, 0.05, 0.2, 1.0, -0.1, 0.15, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

geometric_brownian_motion: Diffusion-only baseline. heston: Stochastic volatility without jumps.

References

Merton, R.C. (1976). Option Pricing When Underlying Stock Returns Are Discontinuous. Journal of Financial Economics 3.

ornstein_uhlenbeck(x0, theta, mu, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Ornstein-Uhlenbeck (OU) mean-reverting process.

The OU process is the canonical continuous-time mean-reverting model, widely used for interest rates, pairs-trading spreads, and volatility dynamics:

\[dX_t = \theta\,(\mu - X_t)\,dt + \sigma\,dW_t\]

This implementation uses the exact discretisation (not Euler) so it is accurate for any step size.

Use OU for modelling quantities that revert to a long-run level: interest rate spreads, log-volatility, or cointegrated pairs.

Parameters:

x0 (float) – Initial value.
theta (float) – Mean reversion speed. Higher values mean faster reversion. The half-life is \(\ln(2) / \theta\).
mu (float) – Long-run mean level.
sigma (float) – Volatility (diffusion coefficient).
T (float) – Time horizon.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated paths.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = ornstein_uhlenbeck(0.05, 5.0, 0.03, 0.01, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

cir_process: Mean-reverting process with non-negative constraint. simulate_vasicek: OU-based interest rate model with bond pricing.

cir_process(x0, kappa, theta, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Cox-Ingersoll-Ross (CIR) process.

The CIR process is the standard mean-reverting, non-negative model for interest rates and variance:

\[dX_t = \kappa(\theta - X_t)\,dt + \sigma\sqrt{X_t}\,dW_t\]

The \(\sqrt{X}\) diffusion term ensures the volatility decreases as the process approaches zero, preventing negative values when the Feller condition \(2\kappa\theta \geq \sigma^2\) holds.

This is a plain-array simulation. For analytics (bond prices, Feller diagnostics), use simulate_cir() instead.

Parameters:

x0 (float) – Initial value (must be positive).
kappa (float) – Mean reversion speed.
theta (float) – Long-run mean level.
sigma (float) – Volatility coefficient.
T (float) – Time horizon.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

Array of shape (n_paths, n_steps + 1) with: simulated paths. Values are kept non-negative via reflection.

Return type:

ndarray[tuple[Any, ...], dtype[double]]

Example

>>> paths = cir_process(0.04, 2.0, 0.04, 0.1, 1.0, 252, 1000, seed=42)
>>> paths.shape
(1000, 253)

See also

simulate_cir: CIR with Feller diagnostics and bond pricing. ornstein_uhlenbeck: Gaussian mean-reverting process (can go

negative).

simulate_sabr(f0, sigma0, alpha, beta, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the SABR stochastic volatility model.

The SABR model (Hagan et al., 2002) is the industry standard for modelling interest rate and FX volatility smiles:

\[\begin{split}dF_t &= \sigma_t\,F_t^{\beta}\,dW_t^1 \\ d\sigma_t &= \alpha\,\sigma_t\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

Key parameters:

\(\beta \in [0, 1]\) controls the backbone: 0 = normal, 1 = lognormal, in between = CEV.
\(\alpha\) is the vol-of-vol.
\(\rho\) controls the skew direction.

When to use SABR: Primarily for interest rate derivatives (swaptions, caps/floors) and FX options where the smile shape is well-characterised by the SABR formula.

Calibration guidance:

sigma0: Calibrate from the ATM implied volatility. For swaptions, this is the ATM normal or lognormal vol.
alpha: Controls the curvature of the smile. Higher alpha = more curvature. Typically 0.2 – 0.6.
beta: Often fixed by convention (beta=0.5 for rates, beta=1.0 for FX) and the other parameters calibrated.
rho: Controls the skew. rho < 0 gives a “normal” skew (downside puts more expensive). rho > 0 reverses it. Typically -0.5 to 0.0 for rates.

When to use Heston vs SABR:

SABR for rates/FX: the SABR formula gives a fast, analytical approximation for implied volatility that is widely used on trading desks.
Heston for equities: the Heston model better captures the mean-reverting nature of equity volatility.
Use simulation (this function) when you need full path dynamics, e.g., for path-dependent payoffs or for evaluating the approximation quality.

Parameters:

f0 (float) – Initial forward rate / price.
sigma0 (float) – Initial stochastic volatility level.
alpha (float) – Volatility of volatility (vol-of-vol).
beta (float) – CEV exponent (0 = normal, 1 = lognormal).
rho (float) – Correlation between forward and vol Brownian motions.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

forwards – simulated forward paths, shape (n_steps + 1, n_paths).
vols – simulated stochastic vol paths, shape (n_steps + 1, n_paths).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> # Simulate a 5% swaption with typical SABR params
>>> result = simulate_sabr(0.05, 0.3, 0.4, 0.5, -0.3, 1.0, 252, 1000, seed=42)
>>> result['forwards'].shape
(253, 1000)

References

Hagan, P.S., Kumar, D., Lesniewski, A.S. & Woodward, D.E. (2002). Managing Smile Risk. Wilmott Magazine, 84-108.

simulate_rough_bergomi(spot, xi, eta, H, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the rough Bergomi stochastic volatility model.

The model specifies:

\[\begin{split}dS_t &= \sqrt{v_t}\,S_t\,dW_t^1 \\ v_t &= \xi\,\mathcal{E}\bigl(\eta\,\tilde{W}_t^H\bigr)\end{split}\]

where \(\tilde{W}^H\) is a (Riemann-Liouville) fractional Brownian motion with Hurst parameter \(H < 0.5\), and \(\mathcal{E}\) is the Wick exponential \(\exp(x - \frac{1}{2}\text{Var}(x))\).

The fBM is simulated via the Cholesky method on the exact covariance matrix, which is accurate but O(n_steps^2) in memory.

Heston vs SABR vs rough Bergomi – when to use which:

Heston: Best for equity options when you need a tractable model with mean-reverting vol. Has analytical characteristic function for fast pricing. Cannot fit the short-maturity skew well.
SABR: Best for rates/FX desks. Has a fast closed-form implied vol approximation. Industry standard for swaptions.
Rough Bergomi: State-of-the-art for capturing the empirical power-law decay of ATM skew across maturities. Most realistic for equity index options but expensive to simulate (O(n_steps^2) for the Cholesky approach).

Calibration guidance:

H ~ 0.1: The empirical consensus for equity index volatility (Gatheral, Jaisson & Rosenbaum, 2018).
xi: Set to the current level of implied variance (e.g., 0.04 for 20% vol).
eta ~ 1.5-2.5: Controls the vol-of-vol. Calibrate to the overall level of the smile.
rho ~ -0.7 to -0.9: Leverage effect.

Parameters:

spot (float) – Initial spot price.
xi (float) – Forward variance level (flat forward variance curve). E.g., 0.04 for 20% annualized vol.
eta (float) – Volatility of volatility. Typical range 1.0 – 3.0.
H (float) – Hurst parameter of fractional Brownian motion. Must be in (0, 0.5) for the “rough” regime. Empirical estimates cluster around 0.1.
rho (float) – Correlation between spot and vol Brownian motions. Typically -0.7 to -0.9 for equity indices.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps. O(n_steps^2) memory for the Cholesky covariance matrix, so keep moderate (100-500).
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

prices – simulated asset price paths, shape (n_steps + 1, n_paths).
variances – simulated instantaneous variance paths, shape (n_steps + 1, n_paths).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> # Simulate SPX-like dynamics with rough vol
>>> result = simulate_rough_bergomi(
...     spot=4000, xi=0.04, eta=1.9, H=0.1, rho=-0.7,
...     T=0.5, n_steps=100, n_paths=1000, seed=42)
>>> result['prices'].shape
(101, 1000)

References

Bayer, C., Friz, P. & Gatheral, J. (2016). Pricing Under Rough Volatility. Quantitative Finance 16(6), 887-904.

Gatheral, J., Jaisson, T. & Rosenbaum, M. (2018). Volatility is Rough. Quantitative Finance 18(6), 933-949.

simulate_3_2_model(spot, v0, kappa, theta, epsilon, rho, T, n_steps, n_paths, seed=None)[source]¶

Simulate the 3/2 stochastic volatility model.

\[\begin{split}dS_t &= \sqrt{v_t}\,S_t\,dW_t^1 \\ dv_t &= \kappa\,v_t(\theta - v_t)\,dt + \varepsilon\,v_t^{3/2}\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

When to use the 3/2 model: For VIX derivatives, or when the Heston model under-estimates the vol-of-vol in high-vol regimes. Also useful when calibrating to options on variance swaps.

Parameters:

spot (float) – Initial asset price.
v0 (float) – Initial variance.
kappa (float) – Mean reversion speed (in the V dynamics).
theta (float) – Long-run variance target.
epsilon (float) – Vol-of-vol parameter (coefficient on V^{3/2}).
rho (float) – Correlation between asset and variance Brownian motions.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

prices – simulated asset price paths, shape (n_steps + 1, n_paths).
variances – simulated variance paths, shape (n_steps + 1, n_paths).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[double]]]

Example

>>> result = simulate_3_2_model(100.0, 0.04, 2.0, 0.04, 0.5, -0.7,
...                             1.0, 252, 1000, seed=42)
>>> result['prices'].shape
(253, 1000)
>>> result['variances'].shape
(253, 1000)

simulate_cir(r0, kappa, theta, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Cox-Ingersoll-Ross short-rate model with analytics.

The CIR model is the standard mean-reverting, non-negative process for interest rates:

\[dr_t = \kappa(\theta - r_t)\,dt + \sigma\sqrt{r_t}\,dW_t\]

The Feller condition \(2\kappa\theta \geq \sigma^2\) ensures that the process stays strictly positive. When violated, zero is accessible but reflecting.

Unlike cir_process() (which returns a plain array), this function also computes: - The Feller condition diagnostic - Analytical zero-coupon bond prices P(0, T) - The mean reversion check

When to use CIR: For modelling interest rates, credit intensities, or stochastic variance (it is the variance process in Heston). The guaranteed non-negativity is crucial for these applications.

Parameters:

r0 (float) – Initial short rate (must be positive).
kappa (float) – Mean reversion speed.
theta (float) – Long-run mean level.
sigma (float) – Volatility coefficient.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

paths – simulated rate paths, shape (n_steps + 1, n_paths).
params – dictionary of model parameters including: - kappa, theta, sigma – model parameters. - feller_satisfied – whether 2*kappa*theta >= sigma^2. - feller_ratio – 2*kappa*theta / sigma^2 (>= 1 means satisfied).

Return type:

Example

>>> result = simulate_cir(0.05, 0.5, 0.04, 0.1, 10.0, 252, 1000, seed=42)
>>> result['paths'].shape
(253, 1000)
>>> result['params']['feller_satisfied']
True

simulate_vasicek(r0, kappa, theta, sigma, T, n_steps, n_paths, seed=None)[source]¶

Simulate the Vasicek interest rate model with bond pricing.

The Vasicek (1977) model is the simplest mean-reverting Gaussian model for the short rate:

\[dr_t = \kappa(\theta - r_t)\,dt + \sigma\,dW_t\]

Being Gaussian, the short rate can go negative – this was historically seen as a drawback but is now relevant in the era of negative rates.

The model admits closed-form solutions for zero-coupon bond prices:

\[P(0, T) = A(T)\,\exp(-B(T)\,r_0)\]

where:

\[\begin{split}B(T) &= \frac{1 - e^{-\kappa T}}{\kappa} \\ A(T) &= \exp\left[ \left(\theta - \frac{\sigma^2}{2\kappa^2}\right)(B(T) - T) - \frac{\sigma^2}{4\kappa}B(T)^2 \right]\end{split}\]

When to use Vasicek: Quick-and-dirty interest rate modelling, or when negative rates are possible/desired. For non-negative rates, use CIR instead.

Parameters:

r0 (float) – Initial short rate.
kappa (float) – Mean reversion speed.
theta (float) – Long-run mean level.
sigma (float) – Volatility.
T (float) – Time horizon in years.
n_steps (int) – Number of time steps.
n_paths (int) – Number of simulation paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

paths – simulated rate paths, shape (n_steps + 1, n_paths).
bond_prices – analytical zero-coupon bond prices P(0, t) at each time step, shape (n_steps + 1,).
yield_curve – continuously compounded yields y(0, t) = -ln(P(0,t))/t at each time step (excluding t=0), shape (n_steps,).

Return type:

Example

>>> result = simulate_vasicek(0.05, 0.5, 0.04, 0.01, 10.0, 100, 1000, seed=42)
>>> result['paths'].shape
(101, 1000)
>>> len(result['bond_prices'])
101
>>> len(result['yield_curve'])
100

References

Vasicek, O. (1977). An Equilibrium Characterisation of the Term Structure. Journal of Financial Economics 5(2), 177-188.

Levy Process Pricing¶

Option pricing under Lévy process models.

Implements FFT-based (Carr-Madan) and COS-method pricing for European options under Variance Gamma, NIG, and arbitrary characteristic function models. All implementations are pure numpy/scipy.

vg_european_fft(spot, strike, rf_rate, T, sigma, nu, theta, n_fft=4096)[source]¶

Price a European call under the Variance Gamma model via FFT.

Use this when the Black-Scholes assumption of Gaussian returns is too restrictive and you observe heavier tails or asymmetry in the market-implied distribution.

Parameters:

spot (float) – Current spot price.
strike (float) – Strike price.
rf_rate (float) – Risk-free rate (annualised).
T (float) – Time to maturity in years.
sigma (float) – VG volatility parameter (controls overall dispersion).
nu (float) – Variance rate of the Gamma subordinator (controls kurtosis; as nu -> 0 the model converges to Black-Scholes).
theta (float) – Drift of the VG process (controls skewness; negative = left-skewed returns).
n_fft (int, default: 4096) – FFT grid size (default 4096).

Returns:

European call price.

Return type:

Example

>>> price = vg_european_fft(100, 100, 0.05, 1.0, 0.2, 0.5, -0.1)
>>> price > 0
True

See also

nig_european_fft: NIG model pricing via FFT. fft_option_price: Generic FFT pricer for any characteristic fn. vg_characteristic: Reusable VG characteristic function.

nig_european_fft(spot, strike, rf_rate, T, alpha, beta, mu, delta, n_fft=4096)[source]¶

Price a European call under the Normal Inverse Gaussian model via FFT.

Parameters:

spot (float) – Current spot price.
strike (float) – Strike price.
rf_rate (float) – Risk-free rate (annualised).
T (float) – Time to maturity in years.
alpha (float) – Tail heaviness parameter (alpha > 0). Larger values produce lighter tails.
beta (float) – Skewness parameter (-alpha < beta < alpha). Negative beta produces left-skewed returns.
mu (float) – Location parameter.
delta (float) – Scale parameter (delta > 0).
n_fft (int, default: 4096) – FFT grid size (default 4096).

Returns:

European call price.

Return type:

Example

>>> price = nig_european_fft(100, 100, 0.05, 1.0, 15.0, -3.0, 0.0, 0.5)
>>> price > 0
True

See also

vg_european_fft: Variance Gamma model pricing via FFT. fft_option_price: Generic FFT pricer for any characteristic fn. nig_characteristic: Reusable NIG characteristic function.

fft_option_price(char_fn, spot, strike, rf_rate, T, n_fft=4096, alpha_damp=1.5)[source]¶

Price European call option(s) via the Carr-Madan FFT method.

The core idea: multiply the modified characteristic function by Simpson’s rule weights, apply the FFT, and interpolate to the desired strike(s).

Parameters:

char_fn (Callable[[ndarray], ndarray]) – Characteristic function of \(\ln(S_T)\) under the risk-neutral measure. Signature: char_fn(u) -> complex ndarray.
spot (float) – Current spot price.
strike (float | ndarray) – Strike price(s). Accepts a scalar or an array for vectorised pricing.
rf_rate (float) – Risk-free rate (annualised, continuously compounded).
T (float) – Time to maturity (years).
n_fft (int, default: 4096) – Number of FFT points (default 4096). Must be a power of 2 for optimal FFT performance.
alpha_damp (float, default: 1.5) – Carr-Madan dampening parameter (default 1.5). Controls the integrand’s decay; values in [1, 2] work well for most models.

Returns:

European call price(s). Returns a float when: strike is scalar, otherwise an ndarray.

Return type:

float | ndarray

Example

>>> import numpy as np
>>> sigma, r, T = 0.2, 0.05, 1.0
>>> log_S = np.log(100.0)
>>> drift = (r - 0.5 * sigma**2) * T
>>> char_fn = lambda u: np.exp(1j*u*(log_S + drift) - 0.5*sigma**2*T*u**2)
>>> price = fft_option_price(char_fn, 100.0, 100.0, r, T)
>>> 9.0 < price < 12.0
True

See also

cos_method: Alternative Fourier method (faster for single strikes). characteristic_function_price: Unified interface to FFT and COS.

References

Carr, P. & Madan, D.B. (1999). Option Valuation Using the Fast Fourier Transform. Journal of Computational Finance 2(4).

cos_method(char_fn, spot, strike, rf_rate, T, n_terms=256, L=10.0)[source]¶

Price a European call option using the COS (Fourier-cosine) method.

The price is computed as:

\[C = e^{-rT} \sum_{k=0}^{N-1} ' \text{Re}\bigl[\phi(u_k)\, e^{-iu_k a}\bigr]\,V_k\]

where \(u_k = k\pi/(b-a)\) and \(V_k\) are the cosine coefficients of the call payoff.

Parameters:

char_fn (Callable[[ndarray], ndarray]) – Characteristic function of \(\ln(S_T)\). Signature: char_fn(u) -> complex ndarray.
spot (float) – Current spot price.
strike (float) – Strike price.
rf_rate (float) – Risk-free rate (annualised).
T (float) – Time to maturity in years.
n_terms (int, default: 256) – Number of cosine expansion terms (default 256). 64–256 terms suffice for most models.
L (float, default: 10.0) – Truncation range controlling the integration domain [log(spot) - L, log(spot) + L] (default 10).

Returns:

European call price.

Return type:

See also

fft_option_price: FFT-based pricing (better for strike grids). characteristic_function_price: Unified interface.

References

Fang, F. & Oosterlee, C.W. (2008). A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions. SIAM Journal on Scientific Computing 31(2).

Characteristic Functions¶

Characteristic function methods for option pricing.

Characteristic functions provide a unified framework for pricing European options under a wide variety of models. Given the characteristic function phi(u) of the log-price log(S_T), the call price can be recovered via either:

FFT (Carr-Madan) – Fourier transform of the modified payoff, evaluated efficiently via the Fast Fourier Transform.
COS (Fang-Oosterlee) – Cosine expansion of the risk-neutral density, integrated analytically against the payoff coefficients.

This module provides model-specific characteristic functions for:

Heston stochastic volatility – the workhorse model for equity vol smiles. Captures mean-reverting variance with leverage effect.
Variance Gamma (VG) – a pure-jump Lévy process that nests Black-Scholes as a limiting case. Three parameters control volatility, skewness, and kurtosis independently.
Normal Inverse Gaussian (NIG) – a flexible Lévy process with semi-heavy tails; popular in credit and commodity markets.
CGMY – generalisation of VG with a parameter controlling the fine structure of jumps (finite/infinite activity, finite/infinite variation).

All characteristic functions return callables that can be plugged directly into wraquant.price.levy_pricing.fft_option_price() or wraquant.price.levy_pricing.cos_method().

References

Heston (1993). A Closed-Form Solution for Options with Stochastic Volatility. Review of Financial Studies 6(2), 327-343.
Madan, Carr & Chang (1998). The Variance Gamma Process and Option Pricing. European Finance Review 2, 79-105.
Barndorff-Nielsen (1997). Normal Inverse Gaussian Distributions and Stochastic Volatility Modelling. Scandinavian Journal of Statistics 24(1), 1-13.
Carr, Geman, Madan & Yor (2002). The Fine Structure of Asset Returns: An Empirical Investigation. Journal of Business 75(2).

characteristic_function_price(char_fn, spot, strike, rf, T, method='fft', n_terms=4096)[source]¶

Price a European call option given a characteristic function.

Parameters:

char_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]) – Characteristic function of log(S_T) under the risk-neutral measure. Signature: char_fn(u) -> complex ndarray where u is a real-valued array of frequencies.
spot (float) – Current spot price of the underlying.
strike (float) – Strike price.
rf (float) – Risk-free interest rate (annualised, continuously compounded).
T (float) – Time to maturity in years.
method (str, default: 'fft') – Pricing method – "fft" (Carr-Madan, default) or "cos" (Fang-Oosterlee cosine expansion).
n_terms (int, default: 4096) – Number of terms/grid points. For FFT this is the FFT size (default 4096); for COS this is the number of cosine terms.

Return type:

Returns:

European call option price.

Example

>>> import numpy as np
>>> # Black-Scholes via characteristic function
>>> S, K, r, T, sigma = 100.0, 100.0, 0.05, 1.0, 0.2
>>> log_S = np.log(S)
>>> drift = (r - 0.5 * sigma**2) * T
>>> char_fn = lambda u: np.exp(1j*u*(log_S + drift) - 0.5*sigma**2*T*u**2)
>>> price = characteristic_function_price(char_fn, S, K, r, T)
>>> 9.0 < price < 12.0
True

heston_characteristic(v0, kappa, theta, sigma_v, rho, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the Heston model.

The Heston (1993) stochastic volatility model specifies:

\[\begin{split}dS_t &= r\,S_t\,dt + \sqrt{v_t}\,S_t\,dW_t^1 \\ dv_t &= \kappa(\theta - v_t)\,dt + \sigma_v\sqrt{v_t}\,dW_t^2 \\ \text{corr}(dW^1, dW^2) &= \rho\end{split}\]

The characteristic function of \(\ln(S_T)\) has a known closed-form expression involving complex exponentials of the model parameters.

When to use Heston vs Black-Scholes:

Use Black-Scholes when the implied vol surface is flat (no smile/skew) or as a quick approximation.
Use Heston when you observe a volatility smile/skew in the market and need to match observed option prices across strikes. The negative \(\rho\) parameter captures the leverage effect (negative correlation between returns and vol).
Heston is the standard model for equity index options and is widely used for calibration to the vol surface.

Parameters:

v0 (float) – Initial variance \(v_0\).
kappa (float) – Mean reversion speed of variance.
theta (float) – Long-run variance level.
sigma_v (float) – Volatility of variance (vol of vol).
rho (float) – Correlation between price and variance Brownian motions.
rf (float) – Risk-free rate.
T (float) – Time to maturity in years.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) that returns the characteristic function evaluated at frequencies u.

Example

>>> char_fn = heston_characteristic(
...     v0=0.04, kappa=2.0, theta=0.04, sigma_v=0.3,
...     rho=-0.7, rf=0.05, T=1.0, spot=100.0,
... )
>>> import numpy as np
>>> val = char_fn(np.array([0.0]))
>>> np.abs(val[0])  # phi(0) = 1
1.0

References

Heston, S.L. (1993). A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial Studies 6(2), 327-343.

vg_characteristic(sigma, nu, theta_vg, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the Variance Gamma model.

\[\phi(u) = \exp\bigl(iu(r + \omega)T\bigr) \cdot \left(1 - iu\,\theta\,\nu + \tfrac{1}{2}\sigma^2\nu\,u^2\right)^{-T/\nu}\]

where \(\omega = \frac{1}{\nu}\ln(1 - \theta\nu - \sigma^2\nu/2)\) is the convexity correction.

Parameters:

sigma (float) – Volatility parameter of the VG process.
nu (float) – Variance rate of the Gamma subordinator. Controls kurtosis; as nu -> 0 the model converges to Black-Scholes.
theta_vg (float) – Drift parameter of the VG process. Controls skewness; negative values produce left-skewed returns (typical for equity markets).
rf (float) – Risk-free rate.
T (float) – Time to maturity.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) returning the characteristic function.

Example

>>> char_fn = vg_characteristic(0.2, 0.5, -0.1, 0.05, 1.0, spot=100)
>>> import numpy as np
>>> np.abs(char_fn(np.array([0.0]))[0])
1.0

nig_characteristic(alpha, beta, mu, delta, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the NIG model.

\[\phi(u) = \exp\bigl(iu(r + \omega + \mu)T + \delta T\bigl(\gamma_0 - \sqrt{\alpha^2 - (\beta + iu)^2}\bigr)\bigr)\]

where \(\gamma_0 = \sqrt{\alpha^2 - \beta^2}\) and \(\omega = \delta(\sqrt{\alpha^2 - (\beta+1)^2} - \gamma_0) - \mu\) ensures the martingale condition.

Parameters:

alpha (float) – Tail heaviness / steepness parameter (alpha > 0). Larger alpha = lighter tails.
beta (float) – Asymmetry parameter (-alpha < beta < alpha). Negative beta = left skew.
mu (float) – Location parameter.
delta (float) – Scale parameter (delta > 0).
rf (float) – Risk-free rate.
T (float) – Time to maturity.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) returning the characteristic function.

Example

>>> char_fn = nig_characteristic(15.0, -3.0, 0.0, 0.5, 0.05, 1.0,
...                              spot=100)
>>> import numpy as np
>>> np.abs(char_fn(np.array([0.0]))[0])
1.0

cgmy_characteristic(C, G, M, Y_param, rf, T, spot=1.0)[source]¶

Characteristic function of log(S_T) under the CGMY model.

The CGMY model (Carr, Geman, Madan & Yor, 2002) generalises the Variance Gamma process with an additional parameter Y controlling the fine structure of jumps:

Y < 0: finite activity (finitely many jumps in any interval)
0 <= Y < 1: infinite activity, finite variation
1 <= Y < 2: infinite activity, infinite variation

\[\phi(u) = \exp\bigl(iu(r + \omega)T + C\,T\,\Gamma(-Y)\bigl[ (M - iu)^Y - M^Y + (G + iu)^Y - G^Y \bigr]\bigr)\]

where \(\omega = -C\,\Gamma(-Y)[(M-1)^Y - M^Y + (G+1)^Y - G^Y]\) is the convexity correction.

Parameters:

C (float) – Overall activity level (C > 0).
G (float) – Rate of exponential decay of the positive jump density (G > 0). Larger G = fewer large positive jumps.
M (float) – Rate of exponential decay of the negative jump density (M > 0). Larger M = fewer large negative jumps.
Y_param (float) – Fine structure parameter (Y < 2, Y != 0, Y != 1). See classification above.
rf (float) – Risk-free rate.
T (float) – Time to maturity.
spot (float, default: 1.0) – Current spot price (default 1.0).

Return type:

Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[cdouble]]]

Returns:

Callable char_fn(u) returning the characteristic function.

Example

>>> char_fn = cgmy_characteristic(1.0, 5.0, 10.0, 0.5,
...                               0.05, 1.0, spot=100)
>>> import numpy as np
>>> np.abs(char_fn(np.array([0.0]))[0])
1.0

References

Carr, P., Geman, H., Madan, D.B. & Yor, M. (2002). The Fine Structure of Asset Returns: An Empirical Investigation. Journal of Business 75(2), 305-332.

FBSDE Methods¶

Forward-Backward Stochastic Differential Equation (FBSDE) solvers.

FBSDEs are the modern mathematical framework for derivatives pricing. The forward SDE models the dynamics of the underlying asset, while the backward SDE (BSDE) simultaneously gives the option price Y and the hedge ratio Z at every point in time.

The coupled system is:

Forward: dX_t = mu(t, X_t) dt + sigma(t, X_t) dW_t Backward: dY_t = -f(t, X_t, Y_t, Z_t) dt + Z_t dW_t Terminal: Y_T = g(X_T)

By the Feynman-Kac theorem, the solution (Y_t, Z_t) of the BSDE is related to the PDE solution u(t, x) of the associated parabolic PDE:

Y_t = u(t, X_t), Z_t = sigma(t, X_t) * nabla_x u(t, X_t)

In the Black-Scholes setting the driver f = -r*Y (risk-free discounting) and the terminal condition g is the payoff. The Z process is then sigma * S * Delta, recovering the classical hedge ratio.

References

El Karoui, Peng, Quenez (1997). Backward Stochastic Differential Equations in Finance. Mathematical Finance 7(1), 1-71.
Pardoux & Peng (1990). Adapted solution of a backward stochastic differential equation. Systems & Control Letters 14(1), 55-61.
E, Han, Jentzen (2017). Deep Learning-Based Numerical Methods for High-Dimensional Parabolic PDEs and BSDEs. arXiv:1706.04702.

All implementations are pure numpy/scipy (with optional torch for deep BSDE).

fbsde_european(spot, payoff_fn, drift_fn, vol_fn, rf, T, n_steps=100, n_paths=10000, seed=None)[source]¶

Solve a forward-backward SDE for European option pricing.

The FBSDE system couples the asset dynamics (forward SDE) with the pricing equation (backward SDE):

\[\begin{split}\text{Forward:} \quad dX_t &= \mu(X_t)\,dt + \sigma(X_t)\,dW_t \\ \text{Backward:} \quad dY_t &= -f(t,X_t,Y_t,Z_t)\,dt + Z_t\,dW_t \\ \text{Terminal:} \quad Y_T &= g(X_T)\end{split}\]

The Z process recovered from regression is proportional to sigma(X) * Delta, giving the hedge ratio directly from the BSDE without any finite-difference bumping.

When drift_fn = lambda x: r * x and vol_fn = lambda x: sigma * x with constant sigma, this reduces to the Black-Scholes model and the price converges to the analytical formula.

Parameters:

spot (float) – Initial value of the forward process X_0.
payoff_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Terminal condition g(X_T). Callable mapping an array of terminal values (shape (n_paths,)) to payoffs.
drift_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Drift coefficient mu(x) of the forward SDE. Callable mapping an array of current values to drift values.
vol_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Diffusion coefficient sigma(x) of the forward SDE. Callable mapping an array of current values to volatility values.
rf (float) – Risk-free rate (used as the BSDE driver: f = -r * Y).
T (float) – Time to maturity in years.
n_steps (int, default: 100) – Number of Euler-Maruyama time steps.
n_paths (int, default: 10000) – Number of Monte Carlo paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

price – estimated option price (mean of Y_0).
delta – estimated hedge ratio at t=0 (mean of Z_0 / vol(X_0)).
paths – forward process paths, shape (n_steps + 1, n_paths).
price_process – Y process, shape (n_steps + 1, n_paths).

Return type:

Example

>>> import numpy as np
>>> # Black-Scholes European call
>>> S0, K, r, sigma, T = 100.0, 100.0, 0.05, 0.2, 1.0
>>> payoff = lambda x: np.maximum(x - K, 0.0)
>>> drift = lambda x: r * x
>>> vol = lambda x: sigma * x
>>> result = fbsde_european(S0, payoff, drift, vol, r, T,
...                         n_steps=100, n_paths=50000, seed=42)
>>> 8.0 < result['price'] < 13.0
True

deep_bsde(dim, payoff_fn, drift_fn, vol_fn, rf, T, n_steps=50, n_paths=4096, n_epochs=200, lr=0.001, seed=None)[source]¶

Deep BSDE solver for high-dimensional derivative pricing.

The key idea: parameterise the initial value Y_0 and the Z process at each time step with neural networks (or simpler function approximators). Then minimise the terminal loss:

\[\mathcal{L} = \mathbb{E}\bigl[\bigl| Y_T^{\theta} - g(X_T) \bigr|^2\bigr]\]

where Y_T^theta is obtained by rolling forward the discretised BSDE with the learned Z networks.

When torch is available, uses a proper neural network with Adam optimiser. Otherwise falls back to a simplified linear approximation with scipy.optimize.minimize.

Use for problems with >3 dimensions where PDE methods fail due to the curse of dimensionality.

Parameters:

dim (int) – Number of underlying assets (spatial dimension).
payoff_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Terminal condition g(X_T). Callable mapping an array of shape (n_paths, dim) to payoffs of shape (n_paths,).
drift_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Drift mu(X). Callable mapping (n_paths, dim) to (n_paths, dim).
vol_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Diffusion sigma(X). Callable mapping (n_paths, dim) to (n_paths, dim). Assumed diagonal diffusion.
rf (float) – Risk-free rate.
T (float) – Time to maturity in years.
n_steps (int, default: 50) – Number of time steps.
n_paths (int, default: 4096) – Number of Monte Carlo paths per epoch.
n_epochs (int, default: 200) – Number of training epochs.
lr (float, default: 0.001) – Learning rate.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

price – estimated option price (Y_0).
delta – estimated hedge ratios at t=0, shape (dim,).
loss_history – training loss at each epoch.

Return type:

Example

>>> import numpy as np
>>> # 1D Black-Scholes call (use fbsde_european for 1D; this is for demo)
>>> payoff = lambda x: np.maximum(np.mean(x, axis=1) - 100, 0.0)
>>> drift = lambda x: 0.05 * x
>>> vol = lambda x: 0.2 * x
>>> result = deep_bsde(1, payoff, drift, vol, 0.05, 1.0,
...                    n_steps=20, n_paths=512, n_epochs=50, seed=42)
>>> result['price'] > 0
True

References

reflected_bsde(spot, payoff_fn, drift_fn, vol_fn, rf, T, n_steps=100, n_paths=10000, seed=None)[source]¶

Solve a reflected BSDE (RBSDE) for American option pricing.

American options require the price process Y_t to stay above the intrinsic value (the obstacle) at all times. This is captured by adding a non-decreasing “reflection” process K_t to the BSDE:

\[Y_t = g(X_T) + \int_t^T f(s, X_s, Y_s, Z_s)\,ds - \int_t^T Z_s\,dW_s + K_T - K_t\]

with the constraint \(Y_t \geq h(t, X_t)\) (obstacle) and \(K\) increases only when \(Y_t = h(t, X_t)\).

The resulting price is always >= the European option price, with the difference representing the early exercise premium.

Parameters:

spot (float) – Initial value of the forward process X_0.
payoff_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Payoff function g(x) used both as terminal condition and exercise (obstacle) value at each time step.
drift_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Drift coefficient mu(x) of the forward SDE.
vol_fn (Callable[[ndarray[tuple[Any, ...], dtype[double]]], ndarray[tuple[Any, ...], dtype[double]]]) – Diffusion coefficient sigma(x) of the forward SDE.
rf (float) – Risk-free rate.
T (float) – Time to maturity in years.
n_steps (int, default: 100) – Number of time steps.
n_paths (int, default: 10000) – Number of Monte Carlo paths.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

price – estimated American option price.
exercise_boundary – array of exercise boundary estimates at each time step, shape (n_steps + 1,).
optimal_stopping_time – mean optimal stopping time.

Return type:

Example

>>> import numpy as np
>>> S0, K, r, sigma, T = 100.0, 100.0, 0.05, 0.2, 1.0
>>> payoff = lambda x: np.maximum(K - x, 0.0)  # put
>>> drift = lambda x: r * x
>>> vol = lambda x: sigma * x
>>> result = reflected_bsde(S0, payoff, drift, vol, r, T,
...                         n_steps=100, n_paths=20000, seed=42)
>>> result['price'] > 0
True

Integrations¶

Advanced pricing integrations using optional packages.

Provides wrappers around QuantLib, financePy, rateslib, and py-vollib for bond pricing, option pricing, yield curve construction, and implied volatility computation.

quantlib_bond(face, coupon, maturity, yield_curve, settlement=None, frequency=2)[source]¶

Price a fixed-rate bond using QuantLib.

Parameters:

face (float) – Face (par) value of the bond.
coupon (float) – Annual coupon rate (e.g. 0.05 for 5 %).
maturity (date) – Bond maturity date.
yield_curve (list[tuple[date, float]]) – Pairs of (date, zero_rate) defining the yield curve.
settlement (date | None, default: None) – Settlement date. Defaults to today when None.
frequency (int, default: 2) – Coupon payment frequency per year (2 = semi-annual).

Returns:

Dictionary containing:

clean_price – clean price of the bond.
dirty_price – dirty (full) price of the bond.
ytm – yield to maturity.
duration – Macaulay duration.
convexity – convexity.

Return type:

quantlib_option(spot, strike, vol, rf, maturity, option_type='call')[source]¶

Price a European option using QuantLib’s analytic engine.

Parameters:

spot (float) – Current spot price of the underlying.
strike (float) – Strike price.
vol (float) – Annualised volatility (e.g. 0.20 for 20 %).
rf (float) – Risk-free interest rate (annual, continuously compounded).
maturity (float) – Time to maturity in years.
option_type (str, default: 'call') – 'call' or 'put'.

Returns:

Dictionary containing:

price – option price.
delta – option delta.
gamma – option gamma.
theta – option theta (per year).
vega – option vega (per 1 % vol change).
rho – option rho.

Return type:

quantlib_yield_curve(dates, rates)[source]¶

Construct a yield curve using QuantLib.

Parameters:

dates (list[date]) – Curve node dates (the first date is the reference date).
rates (list[float]) – Zero rates corresponding to each date.

Returns:

Dictionary containing:

reference_date – the curve reference date.
max_date – the latest date on the curve.
discount_factors – list of discount factors at each node.
forward_rates – list of instantaneous forward rates at each node.
curve – the raw ql.ZeroCurve object.

Return type:

financepy_option(spot, strike, vol, rf, maturity)[source]¶

Price a European equity option using FinancePy.

Parameters:

spot (float) – Current spot price.
strike (float) – Strike price.
vol (float) – Annualised volatility.
rf (float) – Risk-free rate (annual).
maturity (float) – Time to maturity in years.

Returns:

Dictionary containing:

call_price – European call price.
put_price – European put price.
call_delta – call delta.
put_delta – put delta.

Return type:

rateslib_swap(notional, fixed_rate, float_spread, maturity)[source]¶

Price a plain vanilla interest rate swap using rateslib.

Parameters:

notional (float) – Notional principal amount.
fixed_rate (float) – Fixed leg coupon rate (e.g. 0.03 for 3 %).
float_spread (float) – Spread over the floating index in basis points.
maturity (str) – Tenor string (e.g. '5Y', '10Y').

Returns:

Dictionary containing:

npv – net present value of the swap.
fixed_rate – the fixed rate used.
float_spread – the floating spread used (bps).
maturity – the tenor string.
notional – the notional amount.

Return type:

vollib_implied_vol(price, spot, strike, rf, maturity, option_type='call')[source]¶

Compute implied volatility using py-vollib.

Parameters:

price (float) – Observed option market price.
spot (float) – Current spot price of the underlying.
strike (float) – Strike price.
rf (float) – Risk-free rate (annual, continuously compounded).
maturity (float) – Time to maturity in years.
option_type (str, default: 'call') – 'call' or 'put'.

Returns:

Dictionary containing:

implied_vol – implied volatility.
option_type – the option type used.
price – the input option price.

Return type:

sdeint_solve(drift_fn, diffusion_fn, y0, tspan, dt=0.01, method='euler')[source]¶

Solve a stochastic differential equation using sdeint.

For SDEs of the form:

dY = f(Y, t) dt + g(Y, t) dW

where f is the drift and g is the diffusion coefficient.

Parameters:

drift_fn (callable) – f(Y, t) -> array. Drift coefficient. Must accept the current state Y and time t and return an array of the same shape as Y.
diffusion_fn (callable) – g(Y, t) -> array. Diffusion coefficient. For scalar SDEs, return a 1-D array of shape (1,). For vector SDEs, return a 2-D array of shape (d, m) where d is the state dimension and m is the number of Wiener processes.
y0 (float | ndarray) – Initial condition. Scalar or 1-D array.
tspan (tuple[float, float]) – (t0, t1) time interval.
dt (float, default: 0.01) – Time step size.
method (str, default: 'euler') –
Integration method:
- 'euler' – Euler-Maruyama (strong order 0.5).
- 'sri2' – Roessler SRI2 Stochastic Runge-Kutta (strong order 1.0). Suitable for arbitrary noise.

Returns:

Dictionary containing:

times – 1-D array of time points.
paths – solution array of shape (n_steps, d) where d is the state dimension.
final_values – state at the terminal time.

Return type: