financialkit.pde_methods#

financialkit.pde_methods hub · History · Examples · API reference

Finite-difference methods for the Black-Scholes equation: the theta-method family (Schwartz’s explicit scheme with its stability check, Laasonen’s implicit scheme and Crank-Nicolson with Rannacher’s implicit start) on uniform or Tavella-Randall spot grids, with tridiagonal systems solved by scipy.linalg.solve_banded() or the Thomas algorithm; knock-out barriers on a barrier-aligned grid; American options as linear complementarity problems solved by projected SOR, the Brennan-Schwartz algorithm or the penalty method; von Neumann amplification factors, the explicit log-spot scheme as a trinomial tree, the transformation to the heat equation and Richardson extrapolation; and the Peaceman-Rachford ADI scheme for two-asset options.

financialkit.pde_methods: finite-difference methods for option pricing.

The Black-Scholes equation on a spot grid by the theta-method: Schwartz’s explicit scheme with its stability check, Laasonen’s implicit scheme and Crank-Nicolson with Rannacher’s implicit start, the tridiagonal systems solved by scipy.linalg.solve_banded() or a hand-written Thomas algorithm, uniform and Tavella-Randall grids, and knock-out barriers on a barrier-aligned grid. American options by Cryer’s projected SOR, the Brennan-Schwartz direct solve and Forsyth-Vetzal’s penalty iteration. The theory around the schemes: von Neumann amplification factors, the explicit log-spot scheme as a trinomial tree, the transformation to the heat equation and Richardson extrapolation. Two assets by the Peaceman-Rachford ADI scheme.

class financialkit.pde_methods.ADIResult(price, S1, S2, V, n_time, method='peaceman_rachford')[source]#

Bases: object

Output of a two-factor alternating-direction implicit solve.

Parameters:
S1: ndarray#

Grid of the first asset.

Type:

ndarray, shape (n1 + 1,)

S2: ndarray#

Grid of the second asset.

Type:

ndarray, shape (n2 + 1,)

V: ndarray#

Values today, V[i, j] at (S1[i], S2[j]).

Type:

ndarray, shape (n1 + 1, n2 + 1)

method: str = 'peaceman_rachford'#

Scheme used.

n_time: int#

Number of full time steps.

price: float#

Option value at the spots (S1, S2).

class financialkit.pde_methods.CrankNicolson(n_space=200, n_time=200, rannacher_steps=0, **kwargs)[source]#

Bases: ThetaScheme

The Crank-Nicolson scheme (theta = 1/2), optionally with Rannacher’s start.

Second-order accurate in time and unconditionally stable, but its amplification factor tends to -1 for the stiffest modes, so a payoff kink leaves undamped oscillations in delta and gamma when the time step is large; rannacher_steps=2 removes them.

Parameters:

  • n_space (int) – Spot intervals and time steps.

  • n_time (int) – Spot intervals and time steps.

  • rannacher_steps (int) – Initial steps replaced by two implicit half steps each.

  • **kwargs – Any other ThetaScheme option except theta.

Examples

>>> res = CrankNicolson(400, 100, rannacher_steps=2).price(100.0, 100.0, 1.0, 0.05, 0.2)
>>> round(res.price, 2)
10.45
class financialkit.pde_methods.ExplicitScheme(n_space=100, n_time=1000, grid='uniform', concentration=0.1, s_max=None)[source]#

Bases: ThetaScheme

Schwartz’s explicit finite-difference scheme (theta = 0).

Each step is a matrix-vector product, and American exercise is a simple max with the payoff. Stable only while \(\Delta t \lesssim 1/(\sigma^2 (M-1)^2)\) on a uniform grid; the result carries the check in FDResult.stability.

Parameters:

Examples

Hull, 10th ed., Example 15.6 (S = 42, K = 40, r = 10%, sigma = 20%, six months):

>>> res = ExplicitScheme(320, 2500, s_max=160.0).price(42.0, 40.0, 0.5, 0.10, 0.20)
>>> round(res.price, 2), res.stability.stable
(4.76, True)
class financialkit.pde_methods.FDResult(price, delta, gamma, S, V, method, times=<factory>, surface=<factory>, exercise_boundary=<factory>, lcp_iterations=0, stability=None)[source]#

Bases: object

Output of a one-factor finite-difference valuation.

Parameters:
S: ndarray#

Spot grid, boundaries included.

Type:

ndarray, shape (M + 1,)

V: ndarray#

Option values on the grid today (t = 0).

Type:

ndarray, shape (M + 1,)

delta: float#

\(\partial V / \partial S\) at the spot, from the same spline.

exercise_boundary: ndarray#

For an American put, the highest grid spot at each time where the option is exercised (for a call, the lowest); nan where no node is exercised. Empty for European options.

Type:

ndarray, shape (N + 1,)

gamma: float#

\(\partial^2 V / \partial S^2\) at the spot, from the same spline.

lcp_iterations: int = 0#

Total iterations of the linear complementarity solver (American only).

method: str#

Scheme used, e.g. "explicit", "implicit", "crank_nicolson".

price: float#

Option value at the spot, interpolated from the grid by a cubic spline.

stability: StabilityResult | None = None#

Stability check of the time step, for explicit schemes.

surface: ndarray#

surface[k] holds the values at times[k].

Type:

ndarray, shape (N + 1, M + 1)

times: ndarray#

Calendar times t from 0 to T.

Type:

ndarray, shape (N + 1,)

class financialkit.pde_methods.FiniteDifferencePricer(n_space=200, n_time=200)[source]#

Bases: ABC

A finite-difference solver for European and American vanilla options.

Parameters:

  • n_space (int) – Number of spot intervals \(M\).

  • n_time (int) – Number of time steps \(N\).

method: str = ''#
abstractmethod price(S, K, T, r, sigma, q=0.0, option_type='call', american=False)[source]#

Value a vanilla option on the grid.

Parameters:

  • S (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • K (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • T (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • r (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • sigma (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • q (float) – Continuous dividend yield.

  • option_type (str)

  • american (bool) – Allow early exercise at every time step.

Return type:

FDResult

Returns:

FDResult

class financialkit.pde_methods.HeatEquationResult(price, delta, gamma, S, V, method, times=<factory>, surface=<factory>, exercise_boundary=<factory>, lcp_iterations=0, stability=None, x=<factory>, tau=<factory>, u=<factory>, alpha=0.0, beta=0.0)[source]#

Bases: FDResult

An FDResult that also keeps the solve in heat-equation variables.

Parameters:
alpha: float = 0.0#

Exponent \(\alpha\) in \(V = K e^{\alpha x + \beta\tau} u\).

beta: float = 0.0#

Exponent \(\beta\) in \(V = K e^{\alpha x + \beta\tau} u\).

tau: ndarray#

Scaled times to maturity \(\tau = \sigma^2 (T - t) / 2\), ascending.

Type:

ndarray, shape (N + 1,)

u: ndarray#

Heat-equation solution, u[n] at tau[n].

Type:

ndarray, shape (N + 1, M + 1)

x: ndarray#

Log-moneyness grid \(x = \ln(S/K)\).

Type:

ndarray, shape (M + 1,)

class financialkit.pde_methods.ImplicitScheme(n_space=200, n_time=200, **kwargs)[source]#

Bases: ThetaScheme

Laasonen’s fully implicit scheme (theta = 1).

Unconditionally stable and free of oscillations, but only first-order accurate in time.

Parameters:

  • n_space (int) – Spot intervals and time steps.

  • n_time (int) – Spot intervals and time steps.

  • **kwargs – Any other ThetaScheme option except theta.

Examples

>>> res = ImplicitScheme(400, 800).price(100.0, 100.0, 1.0, 0.05, 0.2, option_type="put")
>>> round(res.price, 2)  # Black-Scholes: 5.57
5.57
class financialkit.pde_methods.LCPSolution(x, iterations, method)[source]#

Bases: object

Solution of one linear complementarity problem \(Ax \ge b,\ x \ge g,\ (Ax - b)^\top (x - g) = 0\).

Parameters:
iterations: int#

Iterations used (1 for a direct method).

method: str#

"psor", "brennan_schwartz" or "penalty".

x: ndarray#

The solution.

Type:

ndarray, shape (m,)

class financialkit.pde_methods.PeacemanRachfordADI(n1=100, n2=100, n_time=50, width=5.0)[source]#

Bases: object

Two-asset European option by the Peaceman-Rachford ADI scheme.

Parameters:

  • n1 (int) – Number of log-spot intervals for each asset (rounded up to even, so the spots are nodes).

  • n2 (int) – Number of log-spot intervals for each asset (rounded up to even, so the spots are nodes).

  • n_time (int) – Number of full time steps.

  • width (float) – Half-width of each grid in standard deviations \(\sigma_i\sqrt{T}\).

Examples

Margrabe’s exchange option, worth Black-Scholes with spot \(S_1\), strike \(S_2\), rate \(q_2\) and volatility \(\sqrt{\sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2} = 0.2646\):

>>> adi = PeacemanRachfordADI(100, 100, 50)
>>> res = adi.price(100.0, 100.0, 1.0, 0.05, 0.2, 0.3, 0.5, lambda s1, s2: np.maximum(s1 - s2, 0.0))
>>> round(res.price, 1)  # Margrabe: 10.5
10.5
price(S1, S2, T, r, sigma1, sigma2, rho, payoff, q1=0.0, q2=0.0)[source]#

Value a European payoff payoff(S1, S2) at maturity T.

Parameters:

  • S1 (float) – Spots.

  • S2 (float) – Spots.

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

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

  • sigma1 (float) – Volatilities.

  • sigma2 (float) – Volatilities.

  • rho (float) – Correlation of the two Brownian motions.

  • payoff (Callable[[ndarray, ndarray], ndarray]) – Vectorised payoff(S1, S2) on meshgrid arrays.

  • q1 (float) – Dividend yields.

  • q2 (float) – Dividend yields.

Return type:

ADIResult

Returns:

ADIResult

class financialkit.pde_methods.RichardsonResult(values, extrapolated, order, ratio, observed_order=nan)[source]#

Bases: object

A sequence of approximations on refined grids and its Richardson extrapolation.

Parameters:
extrapolated: float#

Extrapolated value from the two finest approximations.

observed_order: float = nan#

Order estimated from the three finest values (nan with fewer than three).

order: float#

Error order \(p\) used in the extrapolation.

ratio: float#

Refinement ratio between consecutive grids.

values: ndarray#

Approximations, coarsest first.

Type:

ndarray, shape (k,)

class financialkit.pde_methods.StabilityResult(number, limit, stable, scheme='', quantity='')[source]#

Bases: object

Outcome of a stability check for one time step.

Parameters:
limit: float#

Largest stable value of number (inf if the scheme is unconditionally stable).

number: float#

a time step, or a dimensionless diffusion number.

Type:

The quantity checked

quantity: str = ''#

"dt" or "diffusion".

scheme: str = ''#

Scheme the limit applies to.

stable: bool#

Whether number <= limit.

class financialkit.pde_methods.ThetaScheme(n_space=200, n_time=200, theta=0.5, rannacher_steps=0, grid='uniform', concentration=0.1, s_max=None, american_method='psor', solver='banded', omega=1.2)[source]#

Bases: FiniteDifferencePricer

The theta-method for European, American and knock-out vanilla options.

Parameters:

  • n_space (int) – Number of spot intervals \(M\).

  • n_time (int) – Number of time steps \(N\).

  • theta (float) – Implicitness in [0, 1].

  • rannacher_steps (int) – Initial steps replaced by two implicit half steps each.

  • grid (str) – Spot grid; "tavella_randall" clusters the nodes at the strike.

  • concentration (float) – Width of the Tavella-Randall cluster as a fraction of the strike.

  • s_max (float | None) – Upper edge of the grid (default default_s_max()).

  • american_method (str) – Linear complementarity solver for American implicit steps.

  • solver (str) – Tridiagonal solver for European implicit steps.

  • omega (float) – PSOR over-relaxation parameter.

Examples

>>> res = ThetaScheme(400, 200, theta=0.5).price(100.0, 100.0, 1.0, 0.05, 0.2)
>>> round(res.price, 2)  # Black-Scholes: 10.45
10.45
price(S, K, T, r, sigma, q=0.0, option_type='call', american=False)[source]#

Value a vanilla option on a grid from 0 to s_max.

Parameters:

  • S (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • K (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • T (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • r (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • sigma (float) – Spot, strike, maturity (years), risk-free rate and volatility.

  • q (float) – Continuous dividend yield.

  • option_type (str)

  • american (bool) – Allow early exercise at every time step.

Return type:

FDResult

Returns:

FDResult

price_barrier(S, K, H, T, r, sigma, q=0.0, option_type='call', barrier_type='down-and-out')[source]#

Value a continuously monitored knock-out option with a node on the barrier.

The grid ends exactly at the barrier H, where the value is set to zero (no rebate), so the barrier condition is imposed without interpolation error.

Parameters:

  • S (float) – Spot, strike, barrier, maturity, rate and volatility.

  • K (float) – Spot, strike, barrier, maturity, rate and volatility.

  • H (float) – Spot, strike, barrier, maturity, rate and volatility.

  • T (float) – Spot, strike, barrier, maturity, rate and volatility.

  • r (float) – Spot, strike, barrier, maturity, rate and volatility.

  • sigma (float) – Spot, strike, barrier, maturity, rate and volatility.

  • q (float) – Continuous dividend yield.

  • option_type (str)

  • barrier_type (str)

Return type:

FDResult

Returns:

FDResult

spot_grid(s_min, s_max, K)[source]#

The spot grid this scheme uses on [s_min, s_max] for strike K.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n_space + 1,)

class financialkit.pde_methods.TrinomialProbabilities(p_down, p_mid, p_up, discount)[source]#

Bases: object

Weights of the explicit log-spot scheme read as a trinomial tree.

Parameters:
discount: float#

One-step discount factor \(1 / (1 + r \Delta t)\).

property is_probability: bool#

Whether all three weights are non-negative.

Type:

bool

p_down: float#

Weight on the node one log-step down.

p_mid: float#

Weight on the same node.

p_up: float#

Weight on the node one log-step up.

property total: float#

Sum of the three weights (1 by construction).

Type:

float

class financialkit.pde_methods.VonNeumannResult(xi, g, theta)[source]#

Bases: object

Amplification factor \(g(\xi)\) of a scheme over the Fourier modes.

Parameters:
g: ndarray#

Complex amplification factor of each mode per step.

Type:

ndarray, shape (n,)

property max_modulus: float#

\(\max_\xi |g(\xi)|\).

Type:

float

property stable: bool#

Whether no mode grows, \(\max |g| \le 1\) (to rounding).

Type:

bool

theta: float#

Implicitness of the theta-scheme (0 explicit, 1/2 Crank-Nicolson, 1 implicit).

xi: ndarray#

Phase angles \(\xi = k \Delta x\) in \([0, \pi]\).

Type:

ndarray, shape (n,)

financialkit.pde_methods.black_scholes_operator(S, sigma, r, q=0.0)[source]#

Tridiagonal coefficients of the Black-Scholes operator at the interior nodes.

Parameters:

  • S (ArrayLike) – Strictly increasing spot grid, boundaries included.

  • sigma (float) – Volatility, risk-free rate and dividend yield.

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

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

Return type:

tuple[ndarray, ndarray, ndarray]

Returns:

lower, diag, upper (ndarray, shape (M - 1,)) – (LV)_i = lower[i] V_{i} + diag[i] V_{i+1} + upper[i] V_{i+2} for interior node i + 1; lower[0] and upper[-1] multiply the boundary values.

Examples

On a uniform grid the coefficients are Wilmott’s \(\tfrac12(\sigma^2 i^2 \mp (r - q) i)\) and \(-(\sigma^2 i^2 + r)\):

>>> lo, di, up = black_scholes_operator([0.0, 1.0, 2.0], 0.2, 0.05)
>>> [round(float(c[0]), 6) for c in (lo, di, up)]
[-0.005, -0.09, 0.045]
financialkit.pde_methods.brennan_schwartz(lower, diag, upper, rhs, obstacle, exercise='low')[source]#

Brennan and Schwartz’s direct solve of a tridiagonal LCP.

Parameters:

  • lower (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • diag (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • upper (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • rhs (ArrayLike)

  • obstacle (ArrayLike)

  • exercise (str) – End of the grid where the exercise region lies: "low" for a put, "high" for a call.

Return type:

LCPSolution

Returns:

LCPSolution

Examples

>>> sol = brennan_schwartz([0, -1, -1], [2, 2, 2], [-1, -1, 0], [1, 0, 1], [1.5, 0, 0])
>>> sol.x.round(6).tolist()
[1.5, 1.333333, 1.166667]
financialkit.pde_methods.default_s_max(S, K, T, sigma)[source]#

Upper edge of the spot grid: \(\max(S, K) \max(4, e^{5\sigma\sqrt{T}})\).

Four times the larger of spot and strike keeps the strike on a node of any uniform grid with a multiple of four intervals; the exponential term reaches five standard deviations for long or volatile options.

Parameters:
Return type:

float

Returns:

float

Examples

>>> default_s_max(100.0, 100.0, 1.0, 0.2)
400.0
financialkit.pde_methods.diffusion_number(sigma, dt, dx)[source]#

The diffusion number \(\sigma^2 \Delta t / \Delta x^2\) of an explicit log-spot scheme.

The explicit scheme in \(x = \ln S\) (equivalently, the trinomial tree it defines) is stable, with a non-negative middle probability, when this number is at most 1.

Parameters:
Return type:

StabilityResult

Returns:

StabilityResult

Examples

>>> diffusion_number(0.2, 0.01, 0.02).number
1.0
financialkit.pde_methods.discounted_intrinsic(S, K, tau, r, q=0.0, option_type='call')[source]#

Value of a vanilla option far from the strike: \(\max(\pm(S e^{-q\tau} - K e^{-r\tau}), 0)\).

This is the Black-Scholes price in the limits \(S \to 0\) and \(S \to \infty\), and serves as the Dirichlet boundary condition.

Parameters:

  • S (ArrayLike)

  • K (float) – Strike, time to maturity, rate and dividend yield.

  • tau (float) – Strike, time to maturity, rate and dividend yield.

  • r (float) – Strike, time to maturity, rate and dividend yield.

  • q (float) – Strike, time to maturity, rate and dividend yield.

  • option_type (str)

Return type:

ndarray

Returns:

ndarray

financialkit.pde_methods.explicit_stability(S, dt, sigma, r, q=0.0)[source]#

Check an explicit step against the largest stable step on a spot grid.

The explicit step multiplies the interior values by \(I + \Delta t A\), with \(A\) the discrete Black-Scholes operator, so errors stay bounded only while \(|1 + \Delta t \lambda| \le 1\) for every eigenvalue \(\lambda\) of \(A\). The limit returned is the largest such step, \(\min_\lambda -2\operatorname{Re}\lambda / |\lambda|^2\); on a uniform grid from 0 it is close to \(1/(\sigma^2 (M - 1)^2)\), the CFL-type bound that ties the time step to the square of the spot step.

Parameters:
Return type:

StabilityResult

Returns:

StabilityResult

Examples

>>> res = explicit_stability(np.linspace(0.0, 400.0, 101), 0.002, 0.2, 0.05)
>>> round(res.limit, 5), res.stable
(0.00279, True)
financialkit.pde_methods.explicit_trinomial_probabilities(sigma, r, dt, dx, q=0.0)[source]#

Weights of the explicit log-spot scheme.

Parameters:

  • sigma (float) – Volatility and rate.

  • r (float) – Volatility and rate.

  • dt (float) – Time step and log-spot step.

  • dx (float) – Time step and log-spot step.

  • q (float) – Dividend yield.

Return type:

TrinomialProbabilities

Returns:

TrinomialProbabilities

Examples

With \(\Delta x = \sigma\sqrt{3\Delta t}\) the middle weight is 2/3:

>>> p = explicit_trinomial_probabilities(0.2, 0.05, 0.01, 0.2 * np.sqrt(0.03))
>>> round(p.p_mid, 12), round(p.total, 12), p.is_probability
(0.666666666667, 1.0, True)
financialkit.pde_methods.heat_equation_exponents(r, sigma, q=0.0)[source]#

The exponents \((\alpha, \beta)\) that turn Black-Scholes into \(u_\tau = u_{xx}\).

Parameters:
Return type:

tuple[float, float]

Returns:

alpha, beta (float)

Examples

With \(r = \sigma^2/2\) the drift term vanishes, \(k = 1\):

>>> [round(v, 12) for v in heat_equation_exponents(0.02, 0.2)]
[0.0, -1.0]
financialkit.pde_methods.heat_equation_price(S, K, T, r, sigma, q=0.0, option_type='call', n_space=400, n_time=100, theta=0.5, rannacher_steps=2, width=6.0)[source]#

Price a European option by solving the heat equation and transforming back.

Parameters:

  • S (float) – Spot, strike, maturity, rate and volatility.

  • K (float) – Spot, strike, maturity, rate and volatility.

  • T (float) – Spot, strike, maturity, rate and volatility.

  • r (float) – Spot, strike, maturity, rate and volatility.

  • sigma (float) – Spot, strike, maturity, rate and volatility.

  • q (float) – Continuous dividend yield.

  • option_type (str)

  • n_space (int) – Number of intervals of the uniform \(x\) grid (rounded up to even, so the spot is a node).

  • n_time (int) – Number of steps in \(\tau\).

  • theta (float) – Implicitness of the heat-equation scheme.

  • rannacher_steps (int) – Implicit start steps (see theta_march()).

  • width (float) – Half-width of the grid in standard deviations \(\sigma\sqrt{T}\), added to the distance between spot and strike.

Return type:

HeatEquationResult

Returns:

HeatEquationResult

Examples

>>> res = heat_equation_price(100.0, 100.0, 1.0, 0.05, 0.2)
>>> round(res.price, 2)  # Black-Scholes: 10.45
10.45
financialkit.pde_methods.interpolate_at(S_grid, V, S)[source]#

Value, delta and gamma at S from a cubic spline through the grid values.

Parameters:
Return type:

tuple[float, float, float]

Returns:

price, delta, gamma (float)

financialkit.pde_methods.log_explicit_price(S, K, T, r, sigma, q=0.0, option_type='call', american=False, n_time=200, dx=None)[source]#

Value a vanilla option by the explicit log-spot scheme, a trinomial tree on a fixed grid.

The grid has \(2N + 1\) nodes centred on the spot, so the value at the spot never feels the boundaries: it equals the root of the trinomial tree with the same weights.

Parameters:

  • S (float) – Spot, strike, maturity, rate and volatility.

  • K (float) – Spot, strike, maturity, rate and volatility.

  • T (float) – Spot, strike, maturity, rate and volatility.

  • r (float) – Spot, strike, maturity, rate and volatility.

  • sigma (float) – Spot, strike, maturity, rate and volatility.

  • q (float) – Continuous dividend yield.

  • option_type (str)

  • american (bool) – Exercise early whenever the payoff exceeds the continuation value.

  • n_time (int) – Number of time steps \(N\).

  • dx (float | None) – Log-spot step (default \(\sigma\sqrt{3\Delta t}\)).

Return type:

FDResult

Returns:

FDResult – stability holds the diffusion number \(\sigma^2\Delta t/\Delta x^2\).

Examples

>>> res = log_explicit_price(100.0, 100.0, 1.0, 0.05, 0.2, n_time=500)
>>> round(res.price, 2)  # Black-Scholes: 10.45
10.45
financialkit.pde_methods.penalty_solve(lower, diag, upper, rhs, obstacle, x0=None, penalty=100000000.0, tol=1e-10, max_iter=50)[source]#

Forsyth and Vetzal’s penalty iteration for a tridiagonal LCP.

Iterates \((A + P^k) x^{k+1} = b + P^k g\), where \(P^k\) is diagonal with \(\rho\) on the nodes where \(x^k < g\) and 0 elsewhere, until the set of penalised nodes stops changing. The solution satisfies the constraint to within \(O(1/\rho)\).

Parameters:

  • lower (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • diag (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • upper (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • rhs (ArrayLike)

  • obstacle (ArrayLike)

  • x0 (ArrayLike | None) – Starting point (default: the unconstrained solution).

  • penalty (float) – Penalty factor \(\rho\).

  • tol (float) – Relative change below which the iteration stops.

  • max_iter (int)

Return type:

LCPSolution

Returns:

LCPSolution

Examples

>>> sol = penalty_solve([0, -1, -1], [2, 2, 2], [-1, -1, 0], [1, 0, 1], [1.5, 0, 0])
>>> sol.x.round(6).tolist()
[1.5, 1.333333, 1.166667]
financialkit.pde_methods.psor(lower, diag, upper, rhs, obstacle, x0=None, omega=1.2, tol=1e-08, max_iter=1000)[source]#

Cryer’s projected SOR for a tridiagonal LCP.

Parameters:

  • lower (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • diag (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • upper (ArrayLike) – Tridiagonal matrix (lower[0] and upper[-1] ignored).

  • rhs (ArrayLike)

  • obstacle (ArrayLike) – Lower bound \(g\) on the solution.

  • x0 (ArrayLike | None) – Starting point (default: the obstacle).

  • omega (float) – Over-relaxation parameter in (0, 2).

  • tol (float) – Stop when no component moves by more than tol in a sweep.

  • max_iter (int)

Return type:

LCPSolution

Returns:

LCPSolution

Examples

With a low obstacle the constraint is inactive and PSOR solves A x = b:

>>> sol = psor([0, -1, -1], [2, 2, 2], [-1, -1, 0], [1, 0, 1], [-5, -5, -5])
>>> sol.x.round(6).tolist()
[1.0, 1.0, 1.0]
financialkit.pde_methods.richardson_extrapolation(values, ratio=2.0, order=None)[source]#

Extrapolate approximations on successively refined grids to zero step size.

Parameters:

  • values (ArrayLike) – Approximations on grids refined by ratio each time, coarsest first (at least two).

  • ratio (float) – Refinement ratio between consecutive grids.

  • order (float | None) – Error order \(p\). If omitted, it is estimated from the three finest values when there are three or more, and taken as 2 otherwise.

Return type:

RichardsonResult

Returns:

RichardsonResult

Examples

Values with error \(h^2\) extrapolate exactly:

>>> res = richardson_extrapolation([1.0 + 0.04, 1.0 + 0.01, 1.0 + 0.0025])
>>> round(res.extrapolated, 12), round(res.observed_order, 12)
(1.0, 2.0)
financialkit.pde_methods.solve_tridiagonal(lower, diag, upper, rhs, solver='banded')[source]#

Solve a tridiagonal system with scipy.linalg.solve_banded() or the Thomas kernel.

Parameters:

  • lower (ArrayLike) – lower[0] and upper[-1] are ignored.

  • diag (ArrayLike) – lower[0] and upper[-1] are ignored.

  • upper (ArrayLike) – lower[0] and upper[-1] are ignored.

  • rhs (ArrayLike) – Right-hand side(s); several columns need solver="banded".

  • solver (str)

Return type:

ndarray

Returns:

ndarray, shape like rhs

financialkit.pde_methods.tavella_randall_grid(s_min, s_max, n, center, alpha)[source]#

Spot grid clustered around center by Tavella and Randall’s sinh map.

Parameters:

  • s_min (float) – Grid edges.

  • s_max (float) – Grid edges.

  • n (int) – Number of intervals.

  • center (float) – Point to cluster the nodes around, usually the strike.

  • alpha (float) – Width of the cluster; uniform as alpha -> inf.

Return type:

ndarray

Returns:

ndarray, shape (n + 1,)

Examples

>>> g = tavella_randall_grid(0.0, 400.0, 40, center=100.0, alpha=10.0)
>>> float(g[0]), float(g[-1])
(0.0, 400.0)
>>> bool(np.diff(g).min() < 2.0 < 10.0 < np.diff(g).max())
True
financialkit.pde_methods.theta_march(V0, lower, diag, upper, T, n_time, boundary, theta=0.5, rannacher_steps=0, obstacle=None, lcp=None, solver='banded')[source]#

March grid values from the payoff (tau = 0) to tau = T.

Parameters:

  • V0 (ndarray) – Payoff on the grid, boundaries included.

  • lower (ndarray) – Spatial operator at the interior nodes, as returned by black_scholes_operator().

  • diag (ndarray) – Spatial operator at the interior nodes, as returned by black_scholes_operator().

  • upper (ndarray) – Spatial operator at the interior nodes, as returned by black_scholes_operator().

  • T (float) – Time to maturity.

  • n_time (int) – Number of full time steps.

  • boundary (Callable[[float], tuple[float, float]]) – boundary(tau) returns the Dirichlet values at both ends.

  • theta (float) – Implicitness in [0, 1].

  • rannacher_steps (int) – Number of initial steps replaced by two implicit half steps each.

  • obstacle (ndarray | None) – Exercise value; given, the option is American.

  • lcp (Callable[[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray], LCPSolution] | None) – Linear complementarity solver for implicit American steps.

  • solver (str) – Tridiagonal solver for European implicit steps.

Return type:

tuple[ndarray, int]

Returns:

  • surface (ndarray, shape (n_time + 1, M + 1)) – surface[n] holds the values at tau = n * T / n_time.

  • lcp_iterations (int) – Total iterations of the complementarity solver.

financialkit.pde_methods.thomas_solve(lower, diag, upper, rhs)[source]#

Solve the tridiagonal system A x = rhs by the Thomas algorithm.

Parameters:

  • lower (ArrayLike) – Sub-diagonal, lower[i] multiplying x[i - 1]; lower[0] is ignored.

  • diag (ArrayLike) – Main diagonal.

  • upper (ArrayLike) – Super-diagonal, upper[i] multiplying x[i + 1]; upper[-1] is ignored.

  • rhs (ArrayLike)

Return type:

ndarray

Returns:

ndarray, shape (m,)

Examples

The second-difference matrix tridiag(-1, 2, -1):

>>> x = thomas_solve([0, -1, -1], [2, 2, 2], [-1, -1, 0], [1, 0, 1])
>>> x.round(12).tolist()
[1.0, 1.0, 1.0]
financialkit.pde_methods.tridiagonal_matrix(lower, diag, upper)[source]#

Assemble the sparse matrix of a tridiagonal system in the kernels’ full-length convention.

Parameters:

  • lower (ArrayLike) – lower[0] and upper[-1] are ignored.

  • diag (ArrayLike) – lower[0] and upper[-1] are ignored.

  • upper (ArrayLike) – lower[0] and upper[-1] are ignored.

Return type:

csr_matrix

Returns:

scipy.sparse.csr_matrix, shape (m, m)

financialkit.pde_methods.uniform_grid(s_min, s_max, n)[source]#

Evenly spaced spot grid with n intervals.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n + 1,)

financialkit.pde_methods.von_neumann_amplification(sigma, r, dt, dx, theta=0.5, q=0.0, n_modes=201)[source]#

Amplification factor of the log-spot theta-scheme over the Fourier modes.

Parameters:

  • sigma (float) – Volatility and rate.

  • r (float) – Volatility and rate.

  • dt (float) – Time step and log-spot step.

  • dx (float) – Time step and log-spot step.

  • theta (float) – Implicitness: 0 explicit, 1/2 Crank-Nicolson, 1 implicit.

  • q (float) – Dividend yield.

  • n_modes (int) – Number of phase angles in \([0, \pi]\).

Return type:

VonNeumannResult

Returns:

VonNeumannResult

Examples

Above a diffusion number of 1 the explicit scheme amplifies the highest mode; Crank-Nicolson does not:

>>> von_neumann_amplification(0.2, 0.0, 0.012, 0.02, theta=0.0).stable
False
>>> von_neumann_amplification(0.2, 0.0, 0.012, 0.02, theta=0.5).stable
True
On this page
Show Source