financialkit.stochastic_volatility#

financialkit.stochastic_volatility hub · History · Examples · API reference

Local volatility: Cox’s constant elasticity of variance with Schroder’s noncentral chi-square prices, Dupire’s local volatility from a price surface, and its simulation through financialkit.integrators.euler_maruyama(). Stochastic volatility: the Hull-White series and mixing Monte Carlo, and the Heston model with its two-probability formula. Jumps: Merton’s jump diffusion (series formula), variance gamma (Madan-Seneta, Madan-Carr-Chang), Kou’s double-exponential jumps and Bates. Fourier pricing for every characteristic-function model: Carr-Madan’s FFT, Lewis’s single integral and the COS method. Smiles: Hagan’s SABR formula and calibration, Gatheral’s SVI with Durrleman’s condition, and Lee’s moment formula. Simulation: full-truncation, partial-truncation and Higham-Mao Euler schemes and Andersen’s quadratic-exponential scheme for Heston.

financialkit.stochastic_volatility: beyond constant volatility.

Local volatility: Cox’s constant elasticity of variance with Schroder’s noncentral chi-square prices, and Dupire’s local volatility recovered from a price surface and simulated through financialkit.integrators.euler_maruyama(). Stochastic volatility: Hull and White’s mixing formula, and the Heston model with its two-probability formula. Jumps: Merton’s jump diffusion (series formula), the Madan-Seneta and Madan-Carr-Chang variance gamma process, Kou’s double-exponential jumps and Bates’s Heston-with-jumps. Fourier pricing for every characteristic-function model: Carr-Madan’s FFT, Lewis’s single integral and Fang-Oosterlee’s COS method. Smiles: Hagan’s SABR formula and calibration, Gatheral’s SVI with Durrleman’s arbitrage check, and Lee’s moment formula for the wings. Simulation: truncated Euler schemes (Lord-Koekkoek-van Dijk) and Andersen’s quadratic-exponential scheme for Heston.

class financialkit.stochastic_volatility.BatesModel(v0, kappa, theta, xi, rho, lam, mu_j, delta)[source]#

Bases: CharacteristicFunctionModel

Bates’s (1996) model: Heston’s stochastic variance plus Merton’s lognormal jumps.

Parameters:

Examples

With no jumps it is Heston:

>>> b = BatesModel(0.04, 2.0, 0.04, 0.3, -0.7, lam=0.0, mu_j=-0.1, delta=0.1)
>>> h = HestonModel(0.04, 2.0, 0.04, 0.3, -0.7)
>>> bool(abs(b.price(100.0, 95.0, 1.0, 0.03) - h.price(100.0, 95.0, 1.0, 0.03)) < 1e-12)
True
char_func(u, T)[source]#

Characteristic function \(\phi(u) = E[e^{i u X_T}]\).

Parameters:

  • u (ArrayLike) – Arguments; complex values are allowed wherever the moment generating function exists, which the Fourier pricers need.

  • T (float) – Maturity in years.

Return type:

ndarray

Returns:

ndarray of complex – With \(\phi(0) = 1\) and \(\phi(-i) = E[e^{X_T}] = 1\).

name: str = 'Bates'#
class financialkit.stochastic_volatility.BlackScholesModel(sigma)[source]#

Bases: CharacteristicFunctionModel

Black-Scholes as a characteristic-function model, the reference for the Fourier pricers.

\(X_T = -\tfrac12\sigma^2 T + \sigma W_T\), so \(\phi(u) = \exp(-\tfrac12\sigma^2 T (i u + u^2))\).

Parameters:

sigma (float) – Volatility.

Examples

>>> from financialkit.options import black_scholes_price
>>> m = BlackScholesModel(0.2)
>>> round(m.price(100.0, 100.0, 1.0, 0.05), 6) == round(black_scholes_price(100.0, 100.0, 1.0, 0.05, 0.2), 6)
True
char_func(u, T)[source]#

Characteristic function \(\phi(u) = E[e^{i u X_T}]\).

Parameters:

  • u (ArrayLike) – Arguments; complex values are allowed wherever the moment generating function exists, which the Fourier pricers need.

  • T (float) – Maturity in years.

Return type:

ndarray

Returns:

ndarray of complex – With \(\phi(0) = 1\) and \(\phi(-i) = E[e^{X_T}] = 1\).

name: str = 'Black-Scholes'#
class financialkit.stochastic_volatility.CharacteristicFunctionModel[source]#

Bases: ABC

A model for \(X_T = \ln(S_T / F_T)\) known through its characteristic function.

Subclasses implement char_func(); pricing, implied volatilities and cumulants follow from it.

abstractmethod char_func(u, T)[source]#

Characteristic function \(\phi(u) = E[e^{i u X_T}]\).

Parameters:

  • u (ArrayLike) – Arguments; complex values are allowed wherever the moment generating function exists, which the Fourier pricers need.

  • T (float) – Maturity in years.

Return type:

ndarray

Returns:

ndarray of complex – With \(\phi(0) = 1\) and \(\phi(-i) = E[e^{X_T}] = 1\).

cumulants(T, h=0.02)[source]#

First, second and fourth cumulants of \(X_T\).

By central differences of the cumulant generating function \(\kappa(s) = \ln \phi(-i s)\) at \(s = 0\). The COS method uses them to size its truncation interval.

Parameters:

  • T (float) – Maturity in years.

  • h (float) – Step in \(s\); the moment generating function must exist on \([-2h, 2h]\).

Return type:

tuple[float, float, float]

Returns:

c1, c2, c4 (float)

implied_vol(S, K, T, r, q=0.0, method='cos')[source]#

Black-Scholes implied volatilities of the model’s prices.

Out-of-the-money options are inverted (puts below the forward, calls above), where the time value is the whole price.

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

name: str = ''#
price(S, K, T, r, q=0.0, option_type='call', method='cos')[source]#

European option prices under the model.

Parameters:

  • S (float) – Spot.

  • K (ArrayLike) – Strikes.

  • T (float) – Maturity, risk-free rate and dividend yield.

  • r (float) – Maturity, risk-free rate and dividend yield.

  • q (float) – Maturity, risk-free rate and dividend yield.

  • option_type (str)

  • method (str) – cos_price(), lewis_price(), heston_price() or carr_madan_price().

Return type:

ndarray | float

Returns:

float or ndarray

class financialkit.stochastic_volatility.FFTResult(strikes, calls, alpha, eta)[source]#

Bases: object

Call prices on the Carr-Madan log-strike grid.

Parameters:
alpha: float#

Damping exponent of the modified call price \(e^{\alpha k} C(k)\).

calls: ndarray#

European call prices at those strikes.

Type:

ndarray, shape (N,)

eta: float#

Spacing of the integration grid; the log-strike spacing is \(2\pi / (N \eta)\).

strikes: ndarray#

Strikes \(e^{k_j}\) on the FFT grid.

Type:

ndarray, shape (N,)

class financialkit.stochastic_volatility.HestonModel(v0, kappa, theta, xi, rho)[source]#

Bases: CharacteristicFunctionModel

Heston’s (1993) square-root stochastic volatility model.

Parameters:

  • v0 (float) – Initial variance.

  • kappa (float) – Speed of mean reversion of the variance.

  • theta (float) – Long-run variance.

  • xi (float) – Volatility of variance.

  • rho (float) – Correlation between the stock and variance Brownian motions.

Examples

>>> m = HestonModel(v0=0.04, kappa=2.0, theta=0.04, xi=0.3, rho=-0.7)
>>> m.feller_satisfied
True
>>> complex(m.char_func(-1j, 1.0)).real  # E[S_T / F_T] = 1
1.0
char_func(u, T)[source]#

Characteristic function of \(\ln(S_T/F_T)\), in the “little trap” form of Albrecher et al. (2007).

That form keeps the complex logarithm on its principal branch for all maturities, unlike Heston’s original expression.

Return type:

ndarray

Parameters:

  • u (ArrayLike)

  • T (float)

property feller_satisfied: bool#

Whether \(2\kappa\theta \ge \xi^2\), so the variance never reaches 0.

Type:

bool

name: str = 'Heston'#
variance_mean(t)[source]#

Expected variance \(E[v_t] = \theta + (v_0 - \theta) e^{-\kappa t}\).

Return type:

ndarray | float

Parameters:

t (ArrayLike)

class financialkit.stochastic_volatility.HestonPaths(times, spot, variance, scheme)[source]#

Bases: object

Simulated spot and variance paths under the Heston model.

Parameters:
scheme: str#

"full_truncation", "partial_truncation", "higham_mao" or "qe".

Type:

The discretisation

spot: ndarray#

Spot paths.

Type:

ndarray, shape (n_paths, n_steps + 1)

times: ndarray#

The time grid.

Type:

ndarray, shape (n_steps + 1,)

variance: ndarray#

Variance paths (negative values can appear under the truncation schemes, which only use their positive part).

Type:

ndarray, shape (n_paths, n_steps + 1)

class financialkit.stochastic_volatility.KouModel(sigma, lam, p, eta1, eta2)[source]#

Bases: CharacteristicFunctionModel

Kou’s (2002) double-exponential jump diffusion.

Log jump sizes are \(\mathrm{Exp}(\eta_1)\) upward with probability \(p\) and \(-\mathrm{Exp}(\eta_2)\) downward otherwise.

Parameters:

  • sigma (float) – Diffusive volatility.

  • lam (float) – Jump intensity.

  • p (float) – Probability that a jump is upward.

  • eta1 (float) – Rate of upward jumps; must exceed 1 so that \(E[J] < \infty\).

  • eta2 (float) – Rate of downward jumps.

Examples

>>> m = KouModel(0.16, 1.0, 0.4, 10.0, 5.0)
>>> m.critical_moments()
(9.0, 5.0)
char_func(u, T)[source]#

Characteristic function \(\phi(u) = E[e^{i u X_T}]\).

Parameters:

  • u (ArrayLike) – Arguments; complex values are allowed wherever the moment generating function exists, which the Fourier pricers need.

  • T (float) – Maturity in years.

Return type:

ndarray

Returns:

ndarray of complex – With \(\phi(0) = 1\) and \(\phi(-i) = E[e^{X_T}] = 1\).

critical_moments()[source]#

Lee’s critical moments \((\eta_1 - 1, \eta_2)\): the exponential tails cap the moments.

Return type:

tuple[float, float]

Returns:

p_right, q_left (float)

property mean_jump: float#

\(\zeta = E[e^Y] - 1\).

Type:

float

name: str = 'Kou'#
sample_terminal(T, n, rng=None)[source]#

Draws of \(X_T\); the up and down jump totals are gamma given their counts.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

class financialkit.stochastic_volatility.LocalVolSurface(strikes, maturities, vol)[source]#

Bases: object

A local volatility \(\sigma(K, T)\) on a strike-maturity grid.

Parameters:
maturities: ndarray#

Increasing maturities.

Type:

ndarray, shape (n_T,)

strikes: ndarray#

Increasing strikes.

Type:

ndarray, shape (n_K,)

vol: ndarray#

Local volatility at each node.

Type:

ndarray, shape (n_T, n_K)

class financialkit.stochastic_volatility.MertonJumpModel(sigma, lam, mu_j, delta)[source]#

Bases: CharacteristicFunctionModel

Merton’s (1976) jump diffusion: GBM plus compound Poisson lognormal jumps.

\(dS/S_- = (r - q - \lambda\bar k)\,dt + \sigma\,dW + (J - 1)\,dN\), with \(\ln J \sim N(\mu_J, \delta^2)\) and \(\bar k = E[J] - 1 = e^{\mu_J + \delta^2/2} - 1\).

Parameters:

  • sigma (float) – Diffusive volatility.

  • lam (float) – Jump intensity (jumps per year).

  • mu_j (float) – Mean log jump size.

  • delta (float) – Standard deviation of the log jump size.

Examples

>>> m = MertonJumpModel(0.2, 1.0, -0.1, 0.15)
>>> round(m.mean_jump, 6)
-0.084926
char_func(u, T)[source]#

Characteristic function \(\phi(u) = E[e^{i u X_T}]\).

Parameters:

  • u (ArrayLike) – Arguments; complex values are allowed wherever the moment generating function exists, which the Fourier pricers need.

  • T (float) – Maturity in years.

Return type:

ndarray

Returns:

ndarray of complex – With \(\phi(0) = 1\) and \(\phi(-i) = E[e^{X_T}] = 1\).

property mean_jump: float#

\(\bar k = e^{\mu_J + \delta^2/2} - 1\), the mean relative jump.

Type:

float

name: str = 'Merton'#
sample_terminal(T, n, rng=None)[source]#

Draws of \(X_T = \ln(S_T/F_T)\): Gaussian diffusion plus a Poisson number of Gaussian jumps.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

class financialkit.stochastic_volatility.SABRFit(alpha, beta, rho, nu, F, T, rmse)[source]#

Bases: object

SABR parameters fitted to a smile at one expiry.

Parameters:
F: float#

Forward.

T: float#

Expiry in years.

alpha: float#

Initial volatility level \(\alpha\).

beta: float#

CEV exponent \(\beta\) (held fixed in the fit).

implied_vol(K)[source]#

Hagan’s implied volatility at strikes K with the fitted parameters.

Return type:

ndarray | float

Parameters:

K (ArrayLike)

nu: float#

Volatility of volatility \(\nu\).

rho: float#

Correlation between the forward and the volatility.

rmse: float#

Root-mean-square implied volatility error of the fit.

class financialkit.stochastic_volatility.SVIFit(a, b, rho, m, sigma, T, rmse)[source]#

Bases: object

Raw SVI parameters fitted to total implied variance at one expiry.

Parameters:
T: float#

Expiry in years.

a: float#

Vertical level of the total variance.

b: float#

Slope of the wings.

implied_vol(k)[source]#

Implied volatility \(\sqrt{w(k) / T}\) at log-moneyness k.

Return type:

ndarray | float

Parameters:

k (ArrayLike)

m: float#

Horizontal shift.

rho: float#

Asymmetry, in \((-1, 1)\).

rmse: float#

Root-mean-square total variance error of the fit.

sigma: float#

Smoothness at the vertex.

total_variance(k)[source]#

Total implied variance \(w(k)\) at log-moneyness k.

Return type:

ndarray | float

Parameters:

k (ArrayLike)

class financialkit.stochastic_volatility.VarianceGammaModel(sigma, nu, theta=0.0)[source]#

Bases: CharacteristicFunctionModel

Variance gamma: Brownian motion with drift run on a gamma clock (Madan-Seneta 1990, Madan-Carr-Chang 1998).

\(X_T = \omega T + \theta G_T + \sigma W_{G_T}\), with \(G_T \sim \Gamma(T/\nu, \nu)\) and \(\omega = \ln(1 - \theta\nu - \tfrac12\sigma^2\nu)/\nu\) so that \(E[e^{X_T}] = 1\). theta = 0 is the symmetric Madan-Seneta process; theta < 0 skews it to the left.

Parameters:

  • sigma (float) – Volatility of the Brownian motion in business time.

  • nu (float) – Variance rate of the gamma clock (excess kurtosis).

  • theta (float) – Drift in business time (skewness).

Examples

>>> m = VarianceGammaModel(sigma=0.12, nu=0.2, theta=-0.14)
>>> round(m.price(100.0, 90.0, 1.0, 0.1), 6)  # Fang-Oosterlee (2008), Table 7
19.099355
char_func(u, T)[source]#

Characteristic function \(\phi(u) = E[e^{i u X_T}]\).

Parameters:

  • u (ArrayLike) – Arguments; complex values are allowed wherever the moment generating function exists, which the Fourier pricers need.

  • T (float) – Maturity in years.

Return type:

ndarray

Returns:

ndarray of complex – With \(\phi(0) = 1\) and \(\phi(-i) = E[e^{X_T}] = 1\).

critical_moments()[source]#

Lee’s critical moments: the largest \(p\) with \(E[S_T^{1+p}] < \infty\) and \(q\) with \(E[S_T^{-q}] < \infty\).

The moment generating function \(E[e^{sX}]\) is finite exactly while \(1 - \theta\nu s - \tfrac12\sigma^2\nu s^2 > 0\), for every maturity.

Return type:

tuple[float, float]

Returns:

p_right, q_left (float)

name: str = 'variance gamma'#
property omega: float#

Martingale correction \(\omega\) per unit time.

Type:

float

pdf(x, T)[source]#

Closed-form density of \(X_T\) (Madan, Carr and Chang 1998, eq. 23).

\[f(z) = \frac{2 e^{\theta z/\sigma^2}}{\nu^{T/\nu}\sqrt{2\pi}\,\sigma\,\Gamma(T/\nu)} \left(\frac{z^2}{2\sigma^2/\nu + \theta^2}\right)^{\frac{T}{2\nu} - \frac14} K_{\frac{T}{\nu} - \frac12}\!\left(\frac{\sqrt{z^2 (2\sigma^2/\nu + \theta^2)}}{\sigma^2}\right),\]

with \(z = x - \omega T\) and \(K\) the modified Bessel function of the second kind.

Parameters:
Return type:

ndarray

Returns:

ndarray

sample_terminal(T, n, rng=None)[source]#

Draws of \(X_T\) by subordination: a gamma time, then a normal given it.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

financialkit.stochastic_volatility.carr_madan_fft(model, S, T, r, q=0.0, alpha=1.5, N=4096, eta=0.25)[source]#

Call prices on a grid of log strikes by the Carr-Madan (1999) FFT.

The damped call \(c(k) = e^{\alpha k} C(k)\) has Fourier transform

\[\psi(v) = \frac{e^{-rT} \phi_{\ln S_T}(v - (\alpha + 1) i)}{\alpha^2 + \alpha - v^2 + i (2\alpha + 1) v},\]

and the inverse transform at \(N\) log strikes spaced \(\lambda = 2\pi/(N\eta)\) around \(\ln F\) is one FFT of \(\psi\) sampled at \(v_j = j\eta\) with Simpson weights.

Parameters:
Return type:

FFTResult

Returns:

FFTResult

Examples

>>> from financialkit.options import black_scholes_price
>>> res = carr_madan_fft(BlackScholesModel(0.2), 100.0, 1.0, 0.05)
>>> i = abs(res.strikes - 100.0).argmin()
>>> bool(abs(res.calls[i] - black_scholes_price(100.0, res.strikes[i], 1.0, 0.05, 0.2)) < 1e-6)
True
financialkit.stochastic_volatility.carr_madan_price(model, S, K, T, r, q=0.0, option_type='call', alpha=1.5, N=4096, eta=0.25)[source]#

European option prices at arbitrary strikes, by cubic interpolation of carr_madan_fft().

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

financialkit.stochastic_volatility.cev_local_vol(S, sigma, beta)[source]#

CEV local volatility \(\sigma S^{\beta - 1}\) (proportional, as in Black-Scholes).

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

Examples

>>> cev_local_vol(100.0, 2.0, 0.5)
0.2
financialkit.stochastic_volatility.cev_price(S, K, T, r, sigma, beta, q=0.0, option_type='call')[source]#

European option prices under Cox’s CEV model (Schroder 1989).

For \(dS = (r - q) S\,dt + \sigma S^{\beta}\,dW\) with \(0 < \beta < 1\), writing \(Q(z; k, \lambda)\) for the noncentral chi-square CDF,

\[C = S e^{-qT}\left[1 - Q(a; b + 2, c)\right] - K e^{-rT} Q(c; b, a),\]

with \(a = (K e^{-(r-q)T})^{2(1-\beta)} / ((1-\beta)^2 v)\), \(b = 1/(1 - \beta)\), \(c = S^{2(1-\beta)} / ((1-\beta)^2 v)\) and \(v = \sigma^2 [e^{2(r-q)(\beta-1)T} - 1] / (2(r-q)(\beta-1))\) (\(\sigma^2 T\) when \(r = q\)). For \(\beta > 1\) the roles swap: \(C = S e^{-qT}[1 - Q(c; -b, a)] - K e^{-rT} Q(a; 2 - b, c)\). Puts follow by parity.

Parameters:

  • S (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • K (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • T (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • r (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • sigma (float) – Volatility scale; the local volatility is \(\sigma S^{\beta - 1}\).

  • beta (float) – Elasticity exponent, positive and not 1 (that is Black-Scholes).

  • q (ArrayLike) – Dividend yield.

  • option_type (str)

Return type:

ndarray | float

Returns:

float or ndarray

Examples

Hull, 10th ed., Section 27.1, uses the same formula; with \(\beta = 0.5\) and \(\sigma = 2\) the at-the-money local volatility is 20%, and the price is close to Black-Scholes at 20%:

>>> round(cev_price(100.0, 100.0, 1.0, 0.05, 2.0, 0.5), 4)
10.4539
financialkit.stochastic_volatility.cos_density(model, x, T, N=None, L=10.0)[source]#

Density of \(X_T = \ln(S_T / F_T)\) by the cosine expansion.

\[f(x) \approx \frac{2}{b - a} \sum_{k=0}^{N-1}{}' \operatorname{Re}\!\left[\phi(u_k) e^{-i u_k a}\right] \cos(u_k (x - a))\]
Parameters:
Return type:

ndarray

Returns:

ndarray

Examples

>>> from scipy.stats import norm
>>> f = cos_density(BlackScholesModel(0.2), [0.0], 1.0)
>>> bool(abs(f[0] - norm.pdf(0.0, -0.02, 0.2)) < 1e-10)
True
financialkit.stochastic_volatility.cos_price(model, S, K, T, r, q=0.0, option_type='call', N=None, L=10.0)[source]#

European option prices by the Fang-Oosterlee (2008) COS method.

With \(y = \ln(S_T/K)\) truncated to \([a_y, b_y]\) and \(u_k = k\pi/(b_y - a_y)\),

\[V = e^{-rT} \sum_{k=0}^{N-1}{}' \operatorname{Re}\!\left[\phi(u_k) e^{-i u_k a}\right] V_k,\]

where \(a = a_y - \ln(F/K)\) is the interval’s left end for \(X_T\), the prime halves the first term, and \(V_k\) are the cosine coefficients of the payoff, in closed form. Puts are priced by the series (their payoff is bounded) and calls by parity, as Fang and Oosterlee recommend.

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

Examples

Fang and Oosterlee (2008), Section 5.2: Heston with S = K = 100, T = 1, r = 0 has the reference call price 5.785155450.

>>> from financialkit.stochastic_volatility import HestonModel
>>> m = HestonModel(v0=0.0175, kappa=1.5768, theta=0.0398, xi=0.5751, rho=-0.5711)
>>> round(cos_price(m, 100.0, 100.0, 1.0, 0.0), 6)
5.785155
financialkit.stochastic_volatility.cos_terms(a, b, per_unit=40.0, minimum=256)[source]#

Default number of COS terms: a power of two of at least per_unit per unit of interval width.

A fixed count is not enough: the fourth cumulant of a long-dated Heston return with a large volatility of variance widens the truncation interval tenfold, and the series needs proportionally more terms to resolve the density on it.

Parameters:

  • a (float) – Truncation interval.

  • b (float) – Truncation interval.

  • per_unit (float) – Terms per unit of \(b - a\).

  • minimum (int)

Return type:

int

Returns:

int

Examples

>>> cos_terms(-3.0, 3.0), cos_terms(-37.0, 37.0)
(256, 4096)
financialkit.stochastic_volatility.dupire_local_vol(strikes, maturities, call_prices, r, q=0.0)[source]#

Dupire’s (1994) local volatility from a grid of call prices.

\[\sigma^2(K, T) = \frac{\partial_T C + (r - q) K\,\partial_K C + q C}{\tfrac12 K^2\,\partial_{KK} C}\]

Derivatives are second-order finite differences (numpy.gradient(), non-uniform grids allowed). Nodes where the numerator or the butterfly \(\partial_{KK}C\) is not positive (an arbitrage, or numerical noise in the wings) get nan.

Parameters:

  • strikes (ArrayLike) – Increasing strikes.

  • maturities (ArrayLike) – Increasing maturities (at least two).

  • call_prices (ArrayLike) – European call prices on the grid.

  • r (float) – Risk-free rate and dividend yield.

  • q (float) – Risk-free rate and dividend yield.

Return type:

LocalVolSurface

Returns:

LocalVolSurface

Examples

Black-Scholes prices give back the flat volatility:

>>> import numpy as np
>>> from financialkit.options import black_scholes_price
>>> K, T = np.linspace(60, 140, 81), np.linspace(0.5, 1.5, 21)
>>> C = black_scholes_price(100.0, K[None, :], T[:, None], 0.03, 0.25)
>>> lv = dupire_local_vol(K, T, C, 0.03)
>>> bool(np.nanmax(np.abs(lv.vol[5:-5, 20:-20] - 0.25)) < 2e-3)
True
financialkit.stochastic_volatility.durrleman_g(k, a, b, rho, m, sigma)[source]#

Durrleman’s condition for SVI: the implied density is \(g(k)\,n(d_2)/\sqrt{w}\).

\[g(k) = \left(1 - \frac{k w'}{2 w}\right)^2 - \frac{w'^2}{4}\left(\frac1w + \frac14\right) + \frac{w''}{2},\]

so a slice is free of butterfly arbitrage exactly when \(g \ge 0\) (Gatheral and Jacquier 2014).

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

Examples

>>> bool(durrleman_g(0.0, 0.04, 0.1, -0.3, 0.0, 0.2) > 0)
True
financialkit.stochastic_volatility.heston_mc_price(model, S, K, T, r, q=0.0, option_type='call', n_steps=100, n_paths=10000, scheme='qe', rng=None)[source]#

European option price under Heston by Monte Carlo on heston_paths().

The standard error measures only the statistical error; the discretisation bias of the scheme comes on top of it.

Parameters:
Return type:

MonteCarloResult

Returns:

MonteCarloResult

financialkit.stochastic_volatility.heston_paths(model, S, T, r, q=0.0, n_steps=100, n_paths=10000, scheme='qe', rng=None, psi_c=1.5)[source]#

Simulate Heston spot and variance paths.

Parameters:

  • model (HestonModel)

  • S (float) – Spot, horizon, risk-free rate and dividend yield.

  • T (float) – Spot, horizon, risk-free rate and dividend yield.

  • r (float) – Spot, horizon, risk-free rate and dividend yield.

  • q (float) – Spot, horizon, risk-free rate and dividend yield.

  • n_steps (int)

  • n_paths (int)

  • scheme (str) – Andersen’s QE scheme, or an Euler scheme on \((\ln S, v)\) with \(v\) in the coefficients replaced by \(v^+\) in both drift and diffusion (full truncation, Lord et al. 2010), in the diffusion only (partial truncation, Deelstra and Delbaen 1998), or by \(|v|\) in the diffusion (Higham and Mao 2005).

  • rng (Generator | int | None)

  • psi_c (float) – QE switching threshold between the quadratic and exponential branches.

Return type:

HestonPaths

Returns:

HestonPaths

Examples

>>> m = HestonModel(v0=0.04, kappa=1.5, theta=0.04, xi=0.5, rho=-0.7)
>>> res = heston_paths(m, 100.0, 1.0, 0.03, n_steps=20, n_paths=20_000, rng=0)
>>> bool(abs(res.spot[:, -1].mean() / (100 * np.exp(0.03)) - 1) < 5e-3)
True
financialkit.stochastic_volatility.heston_price(model, S, K, T, r, q=0.0, option_type='call')[source]#

European prices by Heston’s (1993) two-probability formula.

\[C = S e^{-qT} P_1 - K e^{-rT} P_2, \qquad P_j = \frac12 + \frac1\pi \int_0^\infty \operatorname{Re}\!\left[\frac{e^{i u \ln(F/K)} f_j(u)}{i u}\right] du,\]

with \(f_2 = \phi\) (the exercise probability under the risk-neutral measure) and \(f_1(u) = \phi(u - i)\) (under the stock measure). It works for any characteristic-function model, not only Heston’s.

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

Examples

Andersen (2008), Case I (Feller condition badly violated): the at-the-money call over ten years is worth 13.0847.

>>> m = HestonModel(v0=0.04, kappa=0.5, theta=0.04, xi=1.0, rho=-0.9)
>>> round(heston_price(m, 100.0, 100.0, 10.0, 0.0), 4)
13.0847
financialkit.stochastic_volatility.hull_white_mixing_mc(S, K, T, r, sigma0, xi, q=0.0, option_type='call', n_steps=100, n_paths=20000, rng=None)[source]#

Hull-White price by mixing: average Black-Scholes prices over simulated average variances.

Variance paths \(dV = \xi V\,dW\) come from financialkit.integrators.gbm_paths(); each path’s time-average variance (trapezoidal rule) sets a Black-Scholes volatility, and the price is the mean of those Black-Scholes prices. This conditional Monte Carlo has far less variance than simulating the stock.

Parameters:
Return type:

MonteCarloResult

Returns:

MonteCarloResult

Examples

>>> res = hull_white_mixing_mc(100.0, 100.0, 1.0, 0.0, 0.2, 0.3, n_paths=4000, rng=0)
>>> bool(abs(res.estimate - hull_white_price(100.0, 100.0, 1.0, 0.0, 0.2, 0.3)) < 4 * res.std_error)
True
financialkit.stochastic_volatility.hull_white_price(S, K, T, r, sigma0, xi, q=0.0, option_type='call')[source]#

Hull-White (1987) series for an option under uncorrelated lognormal variance.

With \(dV = \xi V\,dW\) independent of the stock and \(\bar V = \frac1T\int_0^T V_t\,dt\), the price is \(E[C_{BS}(\bar V)]\); expanding about \(V_0 = \sigma_0^2\),

\[f \approx C_{BS}(V_0) + \tfrac12 C''(V_0)\,\operatorname{Var}\bar V + \tfrac16 C'''(V_0)\,E[(\bar V - V_0)^3],\]

with the derivatives in the variance in closed form. The truncated series is accurate while \(\xi^2 T\) is small (below about 0.1) and drifts away from the exact mixing price beyond that; compare hull_white_mixing_mc().

Parameters:

  • S (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • K (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • T (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • r (ArrayLike) – Spot, strike, maturity and risk-free rate.

  • sigma0 (float) – Initial volatility, \(\sqrt{V_0}\).

  • xi (float) – Volatility of the variance.

  • q (ArrayLike) – Dividend yield.

  • option_type (str)

Return type:

ndarray | float

Returns:

float or ndarray

Examples

Without volatility of variance it is Black-Scholes:

>>> hull_white_price(100.0, 100.0, 1.0, 0.05, 0.2, 0.0) == black_scholes_price(100.0, 100.0, 1.0, 0.05, 0.2)
True
financialkit.stochastic_volatility.lee_wing_slope(p)[source]#

Lee’s (2004) moment formula: the limiting wing slope of \(w(k)/|k|\).

\[\beta = 2 - 4\left(\sqrt{p^2 + p} - p\right) \in [0, 2],\]

with \(p\) the critical moment: \(\tilde p\) for the right wing (\(E[S_T^{1+p}] < \infty\)), \(\tilde q\) for the left (\(E[S_T^{-q}] < \infty\)). Infinite moments of all orders give a slope of 0; no moments beyond the first give the maximum, 2.

Parameters:

p (ArrayLike) – Critical moment, non-negative (inf allowed).

Return type:

ndarray | float

Returns:

float or ndarray

Examples

>>> lee_wing_slope([0.0, np.inf])
array([2., 0.])
financialkit.stochastic_volatility.lewis_price(model, S, K, T, r, q=0.0, option_type='call')[source]#

European option prices by Lewis’s (2001) single-integral formula.

\[C = S e^{-qT} - \frac{\sqrt{F K}\, e^{-rT}}{\pi} \int_0^\infty \operatorname{Re}\!\left[e^{i u \ln(F/K)} \phi\!\left(u - \tfrac{i}{2}\right)\right] \frac{du}{u^2 + \tfrac14}\]

The integrand decays like \(u^{-2}\) even when \(\phi\) does not decay, so the formula is robust for pure-jump models. Puts by parity.

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

Examples

>>> from financialkit.options import black_scholes_price
>>> m = BlackScholesModel(0.25)
>>> abs(lewis_price(m, 100.0, 110.0, 0.5, 0.03) - black_scholes_price(100.0, 110.0, 0.5, 0.03, 0.25)) < 1e-8
True
financialkit.stochastic_volatility.local_vol_paths(surface, S, T, r, q=0.0, n_steps=100, n_paths=10000, rng=None, dW=None)[source]#

Simulate \(dS = (r - q) S\,dt + \sigma(S, t) S\,dW\) on a local-volatility grid.

Euler-Maruyama on \(\ln S\), through financialkit.integrators.euler_maruyama(); the volatility is interpolated bilinearly and held flat outside the grid. Nodes that are nan in surface are filled from their nearest finite neighbour in strike first.

Parameters:
Return type:

tuple[ndarray, ndarray]

Returns:

  • times (ndarray, shape (n_steps + 1,))

  • spot (ndarray, shape (n_paths, n_steps + 1))

Examples

>>> import numpy as np
>>> flat = LocalVolSurface(np.array([50.0, 150.0]), np.array([0.0, 1.0]), np.full((2, 2), 0.2))
>>> t, S = local_vol_paths(flat, 100.0, 1.0, 0.05, n_steps=50, n_paths=20_000, rng=0)
>>> bool(abs(S[:, -1].mean() - 100 * np.exp(0.05)) < 0.6)
True
financialkit.stochastic_volatility.merton_jump_price(S, K, T, r, model, q=0.0, option_type='call', n_terms=80)[source]#

Merton’s (1976) series: a Poisson mixture of Black-Scholes prices.

\[C = \sum_{n \ge 0} \frac{e^{-\lambda' T} (\lambda' T)^n}{n!}\, C_{BS}(S, K, T, r_n, \sigma_n, q),\]

with \(\lambda' = \lambda(1 + \bar k)\), \(\sigma_n^2 = \sigma^2 + n\delta^2/T\) and \(r_n = r - \lambda\bar k + n \ln(1 + \bar k)/T\).

Parameters:
Return type:

ndarray | float

Returns:

float or ndarray

Examples

Without jumps the series is Black-Scholes:

>>> m = MertonJumpModel(0.2, 0.0, -0.1, 0.1)
>>> merton_jump_price(100.0, 100.0, 1.0, 0.05, m) == black_scholes_price(100.0, 100.0, 1.0, 0.05, 0.2)
True
financialkit.stochastic_volatility.sabr_calibrate(F, strikes, T, vols, beta=0.5)[source]#

Fit \((\alpha, \rho, \nu)\) to a smile with \(\beta\) fixed, by scipy.optimize.least_squares().

Parameters:

  • F (float) – Forward.

  • strikes (ArrayLike) – Strikes and market implied volatilities.

  • vols (ArrayLike) – Strikes and market implied volatilities.

  • T (float) – Expiry.

  • beta (float) – CEV exponent, conventionally chosen rather than fitted.

Return type:

SABRFit

Returns:

SABRFit

Examples

>>> import numpy as np
>>> K = np.linspace(70, 130, 13)
>>> fit = sabr_calibrate(100.0, K, 1.0, sabr_implied_vol(100.0, K, 1.0, 2.0, 0.5, -0.3, 0.4))
>>> round(fit.alpha, 4), round(fit.rho, 4), round(fit.nu, 4)
(2.0, -0.3, 0.4)
financialkit.stochastic_volatility.sabr_implied_vol(F, K, T, alpha, beta, rho, nu)[source]#

Hagan et al.’s (2002) Black implied volatility of the SABR model.

\[\sigma_B = \frac{\alpha}{(FK)^{(1-\beta)/2}\left[1 + \frac{(1-\beta)^2}{24}\ell^2 + \frac{(1-\beta)^4}{1920}\ell^4\right]} \frac{z}{x(z)}\left[1 + \left(\frac{(1-\beta)^2 \alpha^2}{24 (FK)^{1-\beta}} + \frac{\rho\beta\nu\alpha}{4 (FK)^{(1-\beta)/2}} + \frac{2 - 3\rho^2}{24}\nu^2\right) T\right],\]

with \(\ell = \ln(F/K)\), \(z = (\nu/\alpha)(FK)^{(1-\beta)/2}\ell\) and \(x(z) = \ln[(\sqrt{1 - 2\rho z + z^2} + z - \rho)/(1 - \rho)]\) (\(z/x(z) = 1\) at the money).

Parameters:

  • F (ArrayLike) – Forward, strike and expiry.

  • K (ArrayLike) – Forward, strike and expiry.

  • T (ArrayLike) – Forward, strike and expiry.

  • alpha (float) – SABR parameters: initial volatility, CEV exponent in [0, 1], correlation and volatility of volatility.

  • beta (float) – SABR parameters: initial volatility, CEV exponent in [0, 1], correlation and volatility of volatility.

  • rho (float) – SABR parameters: initial volatility, CEV exponent in [0, 1], correlation and volatility of volatility.

  • nu (float) – SABR parameters: initial volatility, CEV exponent in [0, 1], correlation and volatility of volatility.

Return type:

ndarray | float

Returns:

float or ndarray

Examples

With \(\beta = 1\) and no volatility of volatility it is lognormal with volatility \(\alpha\):

>>> sabr_implied_vol(100.0, [80.0, 100.0, 120.0], 1.0, 0.2, 1.0, 0.0, 0.0)
array([0.2, 0.2, 0.2])
financialkit.stochastic_volatility.svi_fit(k, w, T)[source]#

Fit raw SVI to total implied variances by scipy.optimize.least_squares().

Several starting points are tried and the best fit kept. Bounds keep \(b \ge 0\), \(|\rho| < 1\) and \(\sigma > 0\).

Parameters:

  • k (ArrayLike) – Log-moneyness and total implied variance \(\sigma_{BS}^2 T\).

  • w (ArrayLike) – Log-moneyness and total implied variance \(\sigma_{BS}^2 T\).

  • T (float) – Expiry, stored for SVIFit.implied_vol().

Return type:

SVIFit

Returns:

SVIFit

Examples

>>> import numpy as np
>>> k = np.linspace(-0.5, 0.5, 21)
>>> fit = svi_fit(k, svi_total_variance(k, 0.01, 0.12, -0.6, 0.05, 0.15), 1.0)
>>> bool(fit.rmse < 1e-8)
True
financialkit.stochastic_volatility.svi_total_variance(k, a, b, rho, m, sigma)[source]#

Gatheral’s raw SVI total implied variance \(w(k)\).

Parameters:

  • k (ArrayLike) – Log-moneyness \(\ln(K/F)\).

  • a (float) – Raw SVI parameters.

  • b (float) – Raw SVI parameters.

  • rho (float) – Raw SVI parameters.

  • m (float) – Raw SVI parameters.

  • sigma (float) – Raw SVI parameters.

Return type:

ndarray | float

Returns:

float or ndarray

Examples

>>> round(svi_total_variance(0.0, 0.02, 0.1, -0.5, 0.0, 0.2), 12)
0.04
financialkit.stochastic_volatility.truncation_range(model, T, L=10.0)[source]#

Fang-Oosterlee truncation interval \(c_1 \pm L\sqrt{c_2 + \sqrt{c_4}}\) for \(X_T\).

Parameters:
Return type:

tuple[float, float]

Returns:

a, b (float)

On this page
Show Source