financialkit.integrators#

Shared SDE path simulation used by every domain that needs paths: Euler-Maruyama and Milstein schemes, exact geometric Brownian motion and Ornstein-Uhlenbeck sampling, plain, correlated (Cholesky), antithetic and scrambled-Sobol Brownian increments, Brownian-bridge construction, and empirical strong and weak convergence orders.

SDE path simulation shared across financialkit subpackages.

Drift and diffusion callbacks use the f(x, t, params) -> ndarray convention, with x the batch of paths of shape (n_paths, d), so a single compiled callback is reused across parameter values without recompiling (see financialkit.integrators.sde).

Random numbers are drawn in Python from np.random.default_rng(rng) and passed to the kernels as arrays, so results are identical with and without Numba. Every domain simulates paths through this module rather than reimplementing it.

The scheme kernels are compiled per process and deliberately not cached to disk (no cache=True). Each one takes the callbacks themselves as arguments, so a cache entry is keyed on those callbacks’ dispatchers, and Numba cannot re-save an index holding keys whose dispatchers belonged to an earlier process (ReferenceError: underlying object has vanished).

njit() is re-exported here for writing those callbacks: it is numba.njit() when Numba is installed (pip install financialkit[fast]) and a no-op decorator otherwise, with HAS_NUMBA reporting which.

class financialkit.integrators.ConvergenceResult(dts, errors, standard_errors, order)[source]#

Bases: object

Errors of a scheme at several step sizes and the fitted order.

Parameters:

  • dts (NDArray[float64])

  • errors (NDArray[float64])

  • standard_errors (NDArray[float64])

  • order (float)

dts: NDArray[float64]#

Step sizes, coarsest first.

Type:

ndarray of float

errors: NDArray[float64]#

Monte Carlo estimate of the error at each step size.

Type:

ndarray of float

order: float#

Slope of log(errors) against log(dts).

Type:

float

standard_errors: NDArray[float64]#

Standard error of each entry of errors.

Type:

ndarray of float

financialkit.integrators.antithetic_increments(T, n_steps, n_paths, dim=None, rng=None)[source]#

Antithetic Brownian increments (Hammersley-Morton 1956).

The first n_paths // 2 paths are drawn as in brownian_increments(); the second half are their negatives, so every path is paired with its mirror image and the sample mean of dW is zero up to rounding.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_paths, n_steps) or (n_paths, n_steps, dim)

financialkit.integrators.brownian_bridge(z, T)[source]#

Build Brownian paths from standard normals by Brownian-bridge construction.

Column 0 of z sets the terminal value W(T) = sqrt(T) z_0; each later column fills the midpoint of an interval whose ends are already known, from the bridge law

\[W(t_m) \mid W(t_l), W(t_r) \sim N\left(\frac{(t_r - t_m) W(t_l) + (t_m - t_l) W(t_r)}{t_r - t_l},\; \frac{(t_m - t_l)(t_r - t_m)}{t_r - t_l}\right),\]

bisecting breadth-first (Glasserman 2004, Sec. 3.1). The result has the law of Brownian motion for i.i.d. normal z; the point of the ordering is that the first columns carry most of the variance, which is what quasi-Monte Carlo needs.

Parameters:

  • z (ArrayLike) – Standard normal draws (pseudo- or quasi-random).

  • T (float) – Horizon in years.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_paths, n_steps + 1) – W on the uniform grid 0, T/n_steps, ..., T, with W(0) = 0.

Examples

>>> import numpy as np
>>> W = brownian_bridge(np.random.default_rng(0).standard_normal((100_000, 8)), T=2.0)
>>> bool(np.allclose(W.var(axis=0), np.linspace(0.0, 2.0, 9), atol=0.05))
True
financialkit.integrators.brownian_increments(T, n_steps, n_paths, dim=None, rng=None)[source]#

Independent Brownian increments dW ~ N(0, dt).

Parameters:

  • T (float) – Horizon in years.

  • n_steps (int) – Number of time steps, so dt = T / n_steps.

  • n_paths (int) – Number of sample paths.

  • dim (int | None) – Number of independent Brownian motions. None (default) returns a single one without the trailing axis.

  • rng (Generator | int | None) – Generator or seed.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_paths, n_steps) or (n_paths, n_steps, dim)

financialkit.integrators.correlated_increments(corr, T, n_steps, n_paths, rng=None)[source]#

Correlated Brownian increments by Cholesky factorisation.

With corr = L L^T and independent standard normals Z, the increments sqrt(dt) * L Z have covariance corr * dt.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_paths, n_steps, dim)

Raises:

ValueError – If corr is not a valid correlation matrix.

Examples

>>> import numpy as np
>>> dW = correlated_increments([[1.0, -0.7], [-0.7, 1.0]], T=1.0, n_steps=1, n_paths=200_000, rng=0)
>>> round(float(np.corrcoef(dW[:, 0, 0], dW[:, 0, 1])[0, 1]), 2)
-0.7
financialkit.integrators.euler_maruyama(drift, diffusion, x0, T, n_steps, n_paths, params=None, t0=0.0, rng=None, dW=None)[source]#

Simulate SDE paths with the Euler-Maruyama scheme (Maruyama 1955).

\[X_{n+1} = X_n + a(X_n, t_n)\,\Delta t + b(X_n, t_n)\,\Delta W_n\]

Strong order 1/2 and weak order 1 for smooth coefficients (Kloeden-Platen 1992, Thms. 10.2.2 and 14.5.1).

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

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

  • paths (ndarray of float, shape (n_paths, n_steps + 1) or (n_paths, n_steps + 1, d)) – Scalar x0 drops the trailing state axis.

Examples

Geometric Brownian motion dS = mu S dt + sigma S dW:

>>> import numpy as np
>>> from financialkit.integrators import euler_maruyama, njit
>>> @njit
... def drift(x, t, p):
...     return p[0] * x
>>> @njit
... def diffusion(x, t, p):
...     return p[1] * x
>>> times, paths = euler_maruyama(drift, diffusion, 100.0, T=1.0, n_steps=50, n_paths=20_000, params=[0.05, 0.2], rng=0)
>>> paths.shape
(20000, 51)
>>> bool(abs(paths[:, -1].mean() - 100 * np.exp(0.05)) < 0.5)
True
financialkit.integrators.gbm_paths(s0, mu, sigma, T, n_steps, n_paths, rng=None, dW=None)[source]#

Exact paths of geometric Brownian motion dS = mu S dt + sigma S dW.

Uses the closed-form solution (Samuelson 1965; Black-Scholes 1973)

\[S_t = S_0 \exp\left((\mu - \tfrac12 \sigma^2) t + \sigma W_t\right).\]

Under the risk-neutral measure, pass mu = r - q.

Parameters:

  • s0 (float) – Initial price.

  • mu (float) – Drift (continuously compounded, per year).

  • sigma (float) – Volatility (per square-root year).

  • T (float) – Horizon in years.

  • n_steps (int) – Number of uniform steps.

  • n_paths (int) – Number of paths.

  • rng (Generator | int | None) – Generator or seed; ignored if dW is given.

  • dW (ArrayLike | None) – Brownian increments to use, e.g. the same ones fed to euler_maruyama() for a pathwise comparison.

Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

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

  • paths (ndarray of float, shape (n_paths, n_steps + 1))

Examples

>>> import numpy as np
>>> times, S = gbm_paths(100.0, 0.05, 0.2, T=1.0, n_steps=12, n_paths=100_000, rng=0)
>>> bool(abs(S[:, -1].mean() / (100 * np.exp(0.05)) - 1) < 0.005)
True
financialkit.integrators.milstein(drift, diffusion, diffusion_dx, x0, T, n_steps, n_paths, params=None, t0=0.0, rng=None, dW=None)[source]#

Simulate SDE paths with the Milstein scheme (Milstein 1974).

\[X_{n+1} = X_n + a\,\Delta t + b\,\Delta W_n + \tfrac12\, b\, \partial_x b\,(\Delta W_n^2 - \Delta t)\]

Strong order 1 (Kloeden-Platen 1992, Thm. 10.3.5). The correction uses only diffusion_dx, the derivative of each b_i with respect to its own x_i, so it is exact Milstein when each b_i depends on x_i alone (diagonal noise); otherwise the Levy-area cross terms are dropped.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

times, paths – As in euler_maruyama().

financialkit.integrators.njit(*args, **kws)[source]#

Legacy decorator that is equivalent to the preferred API: jit().

See documentation for jit function/decorator for full description.

Return type:

_JITWrapper

financialkit.integrators.ou_paths(x0, kappa, theta, sigma, T, n_steps, n_paths, rng=None)[source]#

Exact paths of the Ornstein-Uhlenbeck process dX = kappa (theta - X) dt + sigma dW.

Samples the Gaussian transition (Uhlenbeck-Ornstein 1930; the Vasicek 1977 short rate)

\[X_{t+\Delta t} = \theta + (X_t - \theta) e^{-\kappa \Delta t} + \sigma \sqrt{\frac{1 - e^{-2\kappa\Delta t}}{2\kappa}}\, Z,\]

which reduces to arithmetic Brownian motion when kappa = 0.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

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

  • paths (ndarray of float, shape (n_paths, n_steps + 1))

financialkit.integrators.sobol_increments(T, n_steps, n_paths, bridge=True, rng=None)[source]#

Quasi-random Brownian increments from a scrambled Sobol sequence.

Each path uses one n_steps-dimensional point of an Owen-scrambled Sobol sequence (scipy.stats.qmc.Sobol), mapped to standard normals by the inverse normal CDF. With bridge=True the normals are assembled by brownian_bridge(), so the first, best distributed Sobol coordinates fix the coarse shape of each path (Caflisch-Moskowitz 1995); otherwise coordinate i drives step i.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_paths, n_steps)

financialkit.integrators.strong_convergence(drift, diffusion, x0, T, exact, n_steps_list=(8, 16, 32, 64, 128), n_paths=2000, params=None, scheme='euler_maruyama', diffusion_dx=None, rng=None)[source]#

Estimate a scheme’s strong order against a known pathwise solution.

Parameters:
Return type:

ConvergenceResult

Returns:

ConvergenceResult – errors are mean absolute (Euclidean, for d > 1) terminal errors.

Examples

>>> import numpy as np
>>> from financialkit.integrators import njit, strong_convergence
>>> @njit
... def drift(x, t, p):
...     return p[0] * x
>>> @njit
... def diffusion(x, t, p):
...     return p[1] * x
>>> exact = lambda W: np.exp((1.0 - 0.5 * 0.5**2) + 0.5 * W)
>>> res = strong_convergence(drift, diffusion, 1.0, 1.0, exact, params=[1.0, 0.5], rng=0)
>>> bool(abs(res.order - 0.5) < 0.15)  # strong order 1/2
True
financialkit.integrators.weak_convergence(drift, diffusion, x0, T, exact_mean, payoff=None, n_steps_list=(4, 8, 16, 32), n_paths=100000, params=None, scheme='euler_maruyama', diffusion_dx=None, rng=None)[source]#

Estimate a scheme’s weak order against a known expectation.

The weak error is usually small next to the Monte Carlo noise, so use coarse steps and many paths, and check that errors stand well clear of standard_errors before trusting order.

Parameters:
Return type:

ConvergenceResult

Returns:

ConvergenceResult – errors are |mean(payoff(X_T)) - exact_mean|.