Skip to content

quanttoolbox.maths

maths.numerical_diff

Python alternatives

Hybrid: numdifftools uses adaptive step sizing and Richardson extrapolation — meaningfully more accurate than this module's fixed-step approach. Worth using where precision matters; keep this module for the magnitude-scaled step convention already wired into econometrics.estimation/whittle.

quanttoolbox.maths.numerical_diff

Numerical gradient, Jacobian, and Hessian, with an adaptive magnitude-scaled step size.

Ported from QuantToolBox/maths/{numerical_gradient,numerical_hessian, numerical_jacobian,sign_operator}.m

Translation notes:

  • numerical_gradient.m handles three cases based on output shape (scalar-output gradient, vector-output-matching-input gradient, and general Jacobian) via a single dispatch; here this is split into numerical_gradient (scalar-valued fun) and numerical_jacobian (vector-valued fun) for clarity, matching how they're actually called elsewhere in this package.
  • The step size dh is scaled per-parameter by max(|x0_i|, 0.01) * sign(x0_i) (falling back to a fixed direction when x0_i == 0), exactly as in the original, so the step is proportional to each parameter's magnitude rather than a fixed absolute value.
  • method="forward" (default) or "central" difference, matching the original's method=1/2.

numerical_gradient(fun, x0, dh=1e-08, method='forward')

Numerical gradient of a scalar-valued function fun at x0 (each parameter perturbed one at a time, holding the others fixed).

Original: maths/numerical_gradient.m (scalar-output case)

Source code in src/quanttoolbox/maths/numerical_diff.py 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def numerical_gradient(
    fun: Callable[[np.ndarray], float], x0: np.ndarray, dh: float = 1e-8, method: str = "forward"
) -> np.ndarray:
    """Numerical gradient of a scalar-valued function fun at x0 (each
    parameter perturbed one at a time, holding the others fixed).

    Original: maths/numerical_gradient.m (scalar-output case)
    """
    x0 = np.asarray(x0, dtype=float).flatten()
    p = x0.shape[0]
    f0 = float(fun(x0))

    step = _step_size(x0, dh)
    x1, x2 = x0 - step, x0 + step
    dx1, dx2 = x0 - x1, x2 - x0

    f1 = np.zeros(p)
    f2 = np.zeros(p)
    for i in range(p):
        xi2 = x0.copy()
        xi2[i] = x2[i]
        f2[i] = float(fun(xi2))
        if method == "central":
            xi1 = x0.copy()
            xi1[i] = x1[i]
            f1[i] = float(fun(xi1))

    if method == "central":
        dx = dx1 + dx2
        return (f2 - f1) / dx
    return (f2 - f0) / dx2

numerical_hessian(fun, x0, dh=6e-05)

Numerical Hessian of a scalar-valued function fun at x0, via second-order finite differences.

Original: maths/numerical_hessian.m

Source code in src/quanttoolbox/maths/numerical_diff.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def numerical_hessian(
    fun: Callable[[np.ndarray], float], x0: np.ndarray, dh: float = 6e-5
) -> np.ndarray:
    """Numerical Hessian of a scalar-valued function fun at x0, via
    second-order finite differences.

    Original: maths/numerical_hessian.m
    """
    x0 = np.asarray(x0, dtype=float).flatten()
    p = x0.shape[0]
    f0 = float(fun(x0))

    step = _step_size(x0, dh)
    x1 = x0 + step
    dx = x1 - x0
    e = np.diag(dx)

    f1 = np.array([fun(x0 + e[:, i]) for i in range(p)])

    f2 = np.zeros((p, p))
    for i in range(p):
        for j in range(i, p):
            f2[i, j] = fun(x0 + e[:, i] + e[:, j])
            if i != j:
                f2[j, i] = f2[i, j]

    return ((f2 - f1[:, None]) - f1[None, :] + f0) / np.outer(dx, dx)

numerical_jacobian(fun, x0, dh=1e-08, method='forward')

Numerical Jacobian of a vector-valued function fun at x0: fun(x0) has shape (n,), x0 has shape (p,), result has shape (n, p).

Original: maths/{numerical_gradient,numerical_jacobian}.m (Jacobian case)

Source code in src/quanttoolbox/maths/numerical_diff.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def numerical_jacobian(
    fun: Callable[[np.ndarray], np.ndarray],
    x0: np.ndarray,
    dh: float = 1e-8,
    method: str = "forward",
) -> np.ndarray:
    """Numerical Jacobian of a vector-valued function fun at x0: fun(x0) has
    shape (n,), x0 has shape (p,), result has shape (n, p).

    Original: maths/{numerical_gradient,numerical_jacobian}.m (Jacobian case)
    """
    x0 = np.asarray(x0, dtype=float).flatten()
    p = x0.shape[0]
    f0 = np.atleast_1d(fun(x0))
    n = f0.shape[0]

    step = _step_size(x0, dh)
    x1, x2 = x0 - step, x0 + step
    dx1, dx2 = x0 - x1, x2 - x0

    f1 = np.zeros((n, p))
    f2 = np.zeros((n, p))
    for i in range(p):
        xi2 = x0.copy()
        xi2[i] = x2[i]
        f2[:, i] = np.atleast_1d(fun(xi2))
        if method == "central":
            xi1 = x0.copy()
            xi1[i] = x1[i]
            f1[:, i] = np.atleast_1d(fun(xi1))

    if method == "central":
        dx = dx1 + dx2
        return (f2 - f1) / dx[None, :]
    return (f2 - f0[:, None]) / dx2[None, :]

sign_operator(x)

Sign function: 1 if x>0, -1 if x<0, 0 if x==0 (equivalent to numpy.sign, provided here for direct call-site compatibility with the original).

Original: maths/sign_operator.m

Source code in src/quanttoolbox/maths/numerical_diff.py 
137
138
139
140
141
142
143
144
145
def sign_operator(x: np.ndarray) -> np.ndarray:
    """Sign function: 1 if x>0, -1 if x<0, 0 if x==0 (equivalent to
    numpy.sign, provided here for direct call-site compatibility with the
    original).

    Original: maths/sign_operator.m
    """
    x = np.asarray(x)
    return (x > 0).astype(float) - (x < 0).astype(float)

Examples

Numerical gradient and Hessian vs. analytical, near a small coordinate — maths/hess2.py 
"""Translated from Examples/maths/hess2.m -- numerical gradient and Hessian
of the same function as grad2.m, at a point with a very small second
coordinate (x2=1e-5), compared against their known analytical forms."""

import numpy as np

from quanttoolbox.maths.numerical_diff import numerical_gradient, numerical_hessian


def fun(x):
    x1, x2 = x[0], x[1]
    return 3 * x1**2 + 6 * x1 + 7 + np.log(x1) + x1 * x2 + x2**2 + np.exp(x2)


def grad_analytical(x):
    x1, x2 = x[0], x[1]
    return np.array([6 * x1 + 6 + 1.0 / x1 + x2, x1 + 2 * x2 + np.exp(x2)])


def hess_analytical(x):
    x1, x2 = x[0], x[1]
    h = np.zeros((2, 2))
    h[0, 0] = 6 - 1.0 / (x1**2)
    h[1, 0] = 1.0
    h[0, 1] = h[1, 0]
    h[1, 1] = 2 + np.exp(x2)
    return h


x0 = np.array([0.5, 0.00001])
g = grad_analytical(x0)
h = hess_analytical(x0)

g1 = numerical_gradient(fun, x0)
print("Gradient, forward difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

h1 = numerical_hessian(fun, x0, dh=6e-5)
print("\nNumerical Hessian:")
print(h1)
print("\nAnalytical Hessian:")
print(h)
Numerical gradient of a two-variable scalar function — maths/grad2.py 
"""Translated from Examples/maths/grad2.m -- numerical vs. analytical
gradient of a scalar function of two variables that mixes polynomial,
log, and exponential terms."""

import numpy as np

from quanttoolbox.maths.numerical_diff import numerical_gradient


def fun(x):
    x1, x2 = x[0], x[1]
    return 3 * x1**2 + 6 * x1 + 7 + np.log(x1) + x1 * x2 + x2**2 + np.exp(x2)


def grad_analytical(x):
    x1, x2 = x[0], x[1]
    return np.array([6 * x1 + 6 + 1.0 / x1 + x2, x1 + 2 * x2 + np.exp(x2)])


x0 = np.array([0.5, 0.001])
g = grad_analytical(x0)

g1 = numerical_gradient(fun, x0)
print("Forward difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

g1 = numerical_gradient(fun, x0, method="central")
print("Central difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))
Numerical gradient of an elementwise function via a sum trick — maths/grad3.py 
"""Translated from Examples/maths/grad3.m -- numerical vs. analytical
gradient of the separable, elementwise function f(x) = x^2 * exp(x^2/3).

The original MATLAB `numerical_gradient` accepts an elementwise-vectorized
function and returns an elementwise gradient. This package's
`numerical_gradient` is scalar-valued only (see `numerical_jacobian` for
the general vector case), so the elementwise function is summed first --
since each term depends on a single x_i, the gradient of the sum w.r.t.
x_i equals the elementwise derivative at x_i (see also grad1.m's
translation in building_blocks.md, which uses the same trick)."""

import numpy as np

from quanttoolbox.maths.numerical_diff import numerical_gradient


def fun_elementwise(x):
    return x**2 * np.exp(x**2 / 3)


def fun_sum(x):
    return np.sum(fun_elementwise(x))


def grad_analytical(x):
    return 2 * x * np.exp(x**2 / 3) + x**2 * (2 * x / 3) * np.exp(x**2 / 3)


x0 = np.array([2.5, 3.0, 3.5])
g = grad_analytical(x0)

g1 = numerical_gradient(fun_sum, x0)
print("Forward difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

g1 = numerical_gradient(fun_sum, x0, method="central")
print("Central difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))
Numerical gradient, same elementwise-sum trick as grad3 — maths/grad4.py 
"""Translated from Examples/maths/grad4.m -- numerical vs. analytical
gradient of the separable, elementwise function f(x) = 2*x * exp(x^2/3).

Same elementwise-via-sum approach as grad3.py's translation (see that
file's docstring for why)."""

import numpy as np

from quanttoolbox.maths.numerical_diff import numerical_gradient


def fun_elementwise(x):
    return 2 * x * np.exp(x**2 / 3)


def fun_sum(x):
    return np.sum(fun_elementwise(x))


def grad_analytical(x):
    return 2 * np.exp(x**2 / 3) + 2 * x * (2 * x / 3) * np.exp(x**2 / 3)


x0 = np.array([2.5, 3.0, 3.5])
g = grad_analytical(x0)

g1 = numerical_gradient(fun_sum, x0)
print("Forward difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

g1 = numerical_gradient(fun_sum, x0, method="central")
print("Central difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))
Numerical gradient, scalar vs. explicit-sum function forms — maths/grad5.py 
"""Translated from Examples/maths/grad5.m -- two equivalent formulations of
the same gradient: fun(x) = 0.5*x'x + exp(x)'*(1/x) (already scalar-valued,
so it needs no summing trick), and fun2(x) = sum(0.5*x_i^2 + exp(x_i)/x_i)
(explicitly written as a sum of separable terms). Both have gradient
g(x) = x + exp(x) .* (1/x - 1/x^2), confirming the two formulations agree
numerically."""

import numpy as np

from quanttoolbox.maths.numerical_diff import numerical_gradient


def fun(x):
    return 0.5 * x @ x + np.exp(x) @ (1.0 / x)


def fun2(x):
    return np.sum(0.5 * x * x + np.exp(x) / x)


def grad_analytical(x):
    return x + np.exp(x) * (1.0 / x - 1.0 / (x * x))


x0 = np.array([1.0, 2.0, 3.0])
g = grad_analytical(x0)

print("fun(x) = 0.5*x'x + exp(x)'*(1/x)")
g1 = numerical_gradient(fun, x0)
print("Forward difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

g1 = numerical_gradient(fun, x0, method="central")
print("Central difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

print("\nfun2(x) = sum(0.5*x^2 + exp(x)/x)  (same gradient, written as a sum)")
g1 = numerical_gradient(fun2, x0)
print("Forward difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

g1 = numerical_gradient(fun2, x0, method="central")
print("Central difference (numerical, analytical, |diff|):")
print(np.column_stack([g1, g, np.abs(g - g1)]))
print("d =", np.max(np.abs(g - g1)))

maths.simulation

Python alternatives

compute_ewma: hybrid — pandas.DataFrame.ewm() is highly optimized but uses a different parameterization (alpha/span/halflife vs. this module's mean-reversion-rate lambda_ + explicit dt); a small translation layer would be needed to switch. GBM simulation and the Riccati/Lyapunov solvers: keep — the latter already route through scipy.linalg.solve_continuous_are/solve_lyapunov (that is the switch), and GBM simulation has no equally-simple standard-library equivalent (QuantLib-Python is a much heavier dependency).

quanttoolbox.maths.simulation

Geometric Brownian motion simulation, EWMA-based mean/vol estimation, volatility targeting, and continuous algebraic Riccati / Lyapunov equation solvers.

Ported from QuantToolBox/maths/{simulate_gbm,simulate_gbm2, simulate_multi_gbm,compute_ewma,momentum_ewma,volatility_target, algebraic_riccati_equation,lyapunov_equation}.m

Translation notes:

  • algebraic_riccati_equation.m/lyapunov_equation.m hand-roll a Schur-decomposition solver and a Kronecker-product linear solve, respectively, for two classic control-theory equations that scipy.linalg already solves directly and robustly (solve_continuous_are, solve_lyapunov) -- used here instead of porting the custom solvers.
  • simulate_multi_gbm.m's body is byte-identical to simulate_gbm.m (single-asset simulation) despite taking a correlation parameter rho that is never used -- this looks like an incomplete/unfinished implementation in the original rather than an intentional simplification. simulate_multi_gbm here instead implements what the name and signature promise: a genuine N-asset correlated GBM simulation via Cholesky decomposition of the correlation matrix (generalizing simulate_gbm2's 2-asset case, which is correctly implemented in the original and is ported as-is).
  • momentum_ewma simulates and analytically decomposes a simple EWMA-momentum trend-following strategy, matching the original's field names conceptually (renamed to snake_case, e.g. V_t -> v_t/result.v_t) via a MomentumEWMAResult dataclass rather than a MATLAB-style long positional return tuple.

algebraic_riccati_equation(a, b, c)

Solve the continuous algebraic Riccati equation A'X + XA - XBX + C = 0 (B assumed symmetric positive semi-definite).

Original: maths/algebraic_riccati_equation.m (reimplemented via scipy.linalg.solve_continuous_are -- see module docstring)

Source code in src/quanttoolbox/maths/simulation.py 
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
def algebraic_riccati_equation(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> np.ndarray:
    """Solve the continuous algebraic Riccati equation A'X + XA - XBX + C = 0
    (B assumed symmetric positive semi-definite).

    Original: maths/algebraic_riccati_equation.m (reimplemented via
    scipy.linalg.solve_continuous_are -- see module docstring)
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    c = np.asarray(c, dtype=float)
    # scipy solves A'X + XA - X @ B_s @ R^-1 @ B_s' @ X + Q = 0. To recover
    # the original's X @ B @ X term (not X @ B @ R^-1 @ B' @ X), factor
    # B = sqrt(B) @ sqrt(B)' (valid since B is symmetric PSD) and pass
    # B_s = sqrt(B), R = I, so B_s @ R^-1 @ B_s' == B exactly.
    b_sqrt = sqrtm(b).real
    return solve_continuous_are(a, b_sqrt, c, np.eye(a.shape[0]))

compute_ewma(prices, lambda_mu, lambda_sigma=None, dt=1.0 / 260)

Exponentially-weighted moving average mean and volatility of a return series (in the mean-reversion-rate parameterization: larger lambda means faster decay toward the new observation).

Original: maths/compute_ewma.m

Source code in src/quanttoolbox/maths/simulation.py 
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
def compute_ewma(
    prices: np.ndarray, lambda_mu: float, lambda_sigma: float | None = None, dt: float = 1.0 / 260
) -> tuple[np.ndarray, np.ndarray]:
    """Exponentially-weighted moving average mean and volatility of a
    return series (in the mean-reversion-rate parameterization: larger
    lambda means faster decay toward the new observation).

    Original: maths/compute_ewma.m
    """
    from quanttoolbox.backtest.returns import price_to_return

    if lambda_sigma is None:
        lambda_sigma = lambda_mu

    prices = np.asarray(prices, dtype=float)
    if prices.ndim == 1:
        prices = prices[:, None]
    n_rows, n_cols = prices.shape

    r_daily = price_to_return(prices, 1)

    mu_t = np.zeros((n_rows, n_cols))
    sigma0 = np.sqrt(1 / dt) * np.nanstd(r_daily, axis=0, ddof=1)
    sigma2_t = np.zeros((n_rows, n_cols))
    sigma2_t[0, :] = sigma0**2

    r_filled = np.where(np.isnan(r_daily), 0.0, r_daily)

    for i in range(1, n_rows):
        mu_t[i, :] = (1 - lambda_mu * dt) * mu_t[i - 1, :] + lambda_mu * r_filled[i, :]
        sigma2_t[i, :] = (1 - lambda_sigma * dt) * sigma2_t[i - 1, :] + lambda_sigma * r_filled[
            i, :
        ] ** 2

    sigma_t = np.sqrt(sigma2_t)
    missing = np.isnan(prices)
    mu_t = np.where(missing, np.nan, mu_t)
    sigma_t = np.where(missing, np.nan, sigma_t)

    return mu_t, sigma_t

lyapunov_equation(a, c)

Solve the Lyapunov equation AX + XA' = C.

Original: maths/lyapunov_equation.m (reimplemented via scipy.linalg.solve_lyapunov -- see module docstring)

Source code in src/quanttoolbox/maths/simulation.py 
382
383
384
385
386
387
388
389
390
def lyapunov_equation(a: np.ndarray, c: np.ndarray) -> np.ndarray:
    """Solve the Lyapunov equation AX + XA' = C.

    Original: maths/lyapunov_equation.m (reimplemented via
    scipy.linalg.solve_lyapunov -- see module docstring)
    """
    a = np.asarray(a, dtype=float)
    c = np.asarray(c, dtype=float)
    return solve_lyapunov(a, c)

momentum_ewma(prices, alpha, lambda_mu, lambda_sigma=None, dt=1.0 / 260, multiplier=1.0)

Simulate and analytically decompose a simple EWMA-momentum trend-following strategy: exposure e_t = alpha * mu_t (position size proportional to the EWMA-estimated drift mu_t), applied with a one-period lag to next-period returns.

Returns both the actually-realized strategy wealth index (v_t, from compounding the lagged-exposure-weighted returns) and a continuous-time-motivated analytical approximation to it, decomposed into a "gamma" component driven by changes in the squared EWMA drift (g_t) and a "theta"/variance-drag component (g_small_t); their product-equivalent sum is v_tilde_t, which approximates v_t.

Original: maths/momentum_ewma.m

Source code in src/quanttoolbox/maths/simulation.py 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
def momentum_ewma(
    prices: np.ndarray,
    alpha: float,
    lambda_mu: float,
    lambda_sigma: float | None = None,
    dt: float = 1.0 / 260,
    multiplier: float = 1.0,
) -> MomentumEWMAResult:
    """Simulate and analytically decompose a simple EWMA-momentum
    trend-following strategy: exposure e_t = alpha * mu_t (position size
    proportional to the EWMA-estimated drift mu_t), applied with a
    one-period lag to next-period returns.

    Returns both the actually-realized strategy wealth index (``v_t``,
    from compounding the lagged-exposure-weighted returns) and a
    continuous-time-motivated analytical approximation to it, decomposed
    into a "gamma" component driven by changes in the squared EWMA drift
    (``g_t``) and a "theta"/variance-drag component (``g_small_t``); their
    product-equivalent sum is ``v_tilde_t``, which approximates ``v_t``.

    Original: maths/momentum_ewma.m
    """
    from quanttoolbox.backtest.returns import price_to_return

    if lambda_sigma is None:
        lambda_sigma = lambda_mu

    prices = np.asarray(prices, dtype=float)
    if prices.ndim == 1:
        prices = prices[:, None]
    n_rows, n_cols = prices.shape

    mu_t, sigma_t = compute_ewma(prices, lambda_mu, lambda_sigma, dt)

    r_st = price_to_return(prices, 1)
    r_st = np.where(np.isnan(r_st), 0.0, r_st)

    e_lag = np.zeros(n_cols)
    r_vt = np.zeros((n_rows, n_cols))
    r_tilde_gt = np.zeros((n_rows, n_cols))
    r_tilde_g_small_t = np.zeros((n_rows, n_cols))
    e_t = np.zeros((n_rows, n_cols))
    mu_lag = mu_t[0, :]

    for k in range(1, n_rows):
        mu_k = mu_t[k, :]
        sigma_k = sigma_t[k, :]
        sr_k = mu_k / sigma_k

        r_vt[k, :] = e_lag * r_st[k, :]
        r_tilde_gt[k, :] = (0.5 * alpha / lambda_mu) * (mu_k**2 - mu_lag**2) * multiplier
        r_tilde_g_small_t[k, :] = (
            alpha
            * sigma_k**2
            * (sr_k**2 * (1 - 0.5 * alpha * sigma_k**2) - 0.5 * lambda_mu)
            * dt
            * multiplier
        )

        e_lag = alpha * mu_k * multiplier
        e_t[k, :] = e_lag
        mu_lag = mu_k

    r_tilde_vt = r_tilde_gt + r_tilde_g_small_t

    v_t = 100 * np.cumprod(1 + r_vt, axis=0)
    g_t = 100 * np.exp(np.cumsum(r_tilde_gt, axis=0))
    g_small_t = 100 * np.exp(np.cumsum(r_tilde_g_small_t, axis=0))
    v_tilde_t = 100 * np.exp(np.cumsum(r_tilde_vt, axis=0))

    return MomentumEWMAResult(
        r_realized=r_st,
        v_t=v_t,
        v_tilde_t=v_tilde_t,
        g_t=g_t,
        g_small_t=g_small_t,
        r_v_t=r_vt,
        r_tilde_v_t=r_tilde_vt,
        r_tilde_g_t=r_tilde_gt,
        r_tilde_g_small_t=r_tilde_g_small_t,
        e_t=e_t,
    )

simulate_gbm(x0, mu, sigma, t, n_paths, rng=None)

Simulate geometric Brownian motion paths (exact scheme, no discretization bias) for a single asset.

Original: maths/simulate_gbm.m

Returns an (n_times, n_paths) array.

Source code in src/quanttoolbox/maths/simulation.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def simulate_gbm(
    x0: float,
    mu: float,
    sigma: float,
    t: np.ndarray,
    n_paths: int,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Simulate geometric Brownian motion paths (exact scheme, no
    discretization bias) for a single asset.

    Original: maths/simulate_gbm.m

    Returns an (n_times, n_paths) array.
    """
    rng = np.random.default_rng() if rng is None else rng
    t = np.asarray(t, dtype=float).flatten()
    n_t = t.shape[0]

    x = np.zeros((n_t, n_paths))
    x_prev = np.full(n_paths, x0)

    for i in range(n_t):
        dt = t[0] if i == 0 else t[i] - t[i - 1]
        k1 = (mu - 0.5 * sigma**2) * dt
        k2 = sigma * np.sqrt(dt)
        u = rng.standard_normal(n_paths)
        x_prev = x_prev * np.exp(k1 + k2 * u)
        x[i, :] = x_prev

    return x

simulate_gbm2(x01, x02, mu1, mu2, sigma1, sigma2, rho, t, n_paths, rng=None)

Simulate two correlated geometric Brownian motion processes (exact scheme).

Original: maths/simulate_gbm2.m

Source code in src/quanttoolbox/maths/simulation.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def simulate_gbm2(
    x01: float,
    x02: float,
    mu1: float,
    mu2: float,
    sigma1: float,
    sigma2: float,
    rho: float,
    t: np.ndarray,
    n_paths: int,
    rng: np.random.Generator | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate two correlated geometric Brownian motion processes (exact
    scheme).

    Original: maths/simulate_gbm2.m
    """
    rng = np.random.default_rng() if rng is None else rng
    t = np.asarray(t, dtype=float).flatten()
    n_t = t.shape[0]

    x1 = np.zeros((n_t, n_paths))
    x2 = np.zeros((n_t, n_paths))
    x1_prev = np.full(n_paths, x01)
    x2_prev = np.full(n_paths, x02)

    rho2 = rho**2

    for i in range(n_t):
        dt = t[0] if i == 0 else t[i] - t[i - 1]
        k1 = (mu1 - 0.5 * sigma1**2) * dt
        k2 = (mu2 - 0.5 * sigma2**2) * dt
        sqrt_dt = np.sqrt(dt)

        u1 = rng.standard_normal(n_paths)
        u2 = rho * u1 + np.sqrt(1 - rho2) * rng.standard_normal(n_paths)

        x1_prev = x1_prev * np.exp(k1 + sigma1 * sqrt_dt * u1)
        x2_prev = x2_prev * np.exp(k2 + sigma2 * sqrt_dt * u2)

        x1[i, :] = x1_prev
        x2[i, :] = x2_prev

    return x1, x2

simulate_multi_gbm(x0, mu, sigma, rho, t, n_paths, rng=None)

Simulate N correlated geometric Brownian motion processes (exact scheme), via Cholesky decomposition of the correlation matrix rho.

Note: implemented as a genuine N-asset correlated simulator (see module docstring for why this differs from the original, whose simulate_multi_gbm.m body is identical to the uncorrelated single-asset simulate_gbm.m despite taking a correlation parameter).

Original: maths/simulate_multi_gbm.m (reimplemented -- see docstring)

Returns an (n_times, n_assets, n_paths) array.

Source code in src/quanttoolbox/maths/simulation.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
def simulate_multi_gbm(
    x0: np.ndarray,
    mu: np.ndarray,
    sigma: np.ndarray,
    rho: np.ndarray,
    t: np.ndarray,
    n_paths: int,
    rng: np.random.Generator | None = None,
) -> np.ndarray:
    """Simulate N correlated geometric Brownian motion processes (exact
    scheme), via Cholesky decomposition of the correlation matrix rho.

    Note: implemented as a genuine N-asset correlated simulator (see
    module docstring for why this differs from the original, whose
    ``simulate_multi_gbm.m`` body is identical to the uncorrelated
    single-asset ``simulate_gbm.m`` despite taking a correlation
    parameter).

    Original: maths/simulate_multi_gbm.m (reimplemented -- see docstring)

    Returns an (n_times, n_assets, n_paths) array.
    """
    rng = np.random.default_rng() if rng is None else rng
    x0 = np.asarray(x0, dtype=float).flatten()
    mu = np.asarray(mu, dtype=float).flatten()
    sigma = np.asarray(sigma, dtype=float).flatten()
    rho = np.asarray(rho, dtype=float)
    t = np.asarray(t, dtype=float).flatten()
    n_assets = x0.shape[0]
    n_t = t.shape[0]

    chol = np.linalg.cholesky(rho)

    x = np.zeros((n_t, n_assets, n_paths))
    x_prev = np.tile(x0[:, None], (1, n_paths))

    for i in range(n_t):
        dt = t[0] if i == 0 else t[i] - t[i - 1]
        k1 = (mu - 0.5 * sigma**2) * dt
        sqrt_dt = np.sqrt(dt)

        z = rng.standard_normal((n_assets, n_paths))
        u = chol @ z

        x_prev = x_prev * np.exp(k1[:, None] + sigma[:, None] * sqrt_dt * u)
        x[i, :, :] = x_prev

    return x

volatility_target(x_t, lambda_, vol_target, min_leverage=0.0, max_leverage=1.0, dt=1.0 / 260, multiplier=1.0)

Apply a volatility-targeting overlay to a price series: scale daily returns by leverage = vol_target / (EWMA volatility), clipped to [min_leverage, max_leverage] and lagged by one period (leverage is set using yesterday's volatility estimate).

Original: maths/volatility_target.m

Source code in src/quanttoolbox/maths/simulation.py 
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def volatility_target(
    x_t: np.ndarray,
    lambda_: float,
    vol_target: float,
    min_leverage: float | None = 0.0,
    max_leverage: float | None = 1.0,
    dt: float = 1.0 / 260,
    multiplier: float = 1.0,
) -> VolatilityTargetResult:
    """Apply a volatility-targeting overlay to a price series: scale daily
    returns by leverage = vol_target / (EWMA volatility), clipped to
    [min_leverage, max_leverage] and lagged by one period (leverage is
    set using yesterday's volatility estimate).

    Original: maths/volatility_target.m
    """
    from quanttoolbox.backtest.returns import price_to_return

    x_t = np.asarray(x_t, dtype=float)
    if x_t.ndim == 1:
        x_t = x_t[:, None]

    _, sigma_t = compute_ewma(x_t, lambda_, lambda_, dt)
    sigma_t = multiplier * sigma_t

    leverage_t = vol_target / sigma_t
    if min_leverage is not None and max_leverage is not None:
        leverage_t = np.clip(leverage_t, min_leverage, max_leverage)
    elif max_leverage is not None:
        leverage_t = np.minimum(leverage_t, max_leverage)
    elif min_leverage is not None:
        leverage_t = np.maximum(leverage_t, min_leverage)

    leverage_lagged = np.full_like(leverage_t, np.nan)
    leverage_lagged[1:] = leverage_t[:-1]
    leverage_lagged[0] = leverage_lagged[1] if leverage_lagged.shape[0] > 1 else np.nan

    r_t = price_to_return(x_t, 1)
    r_filled = np.where(np.isnan(r_t), 0.0, r_t)
    r_levered = leverage_lagged * r_filled

    y_t = 100 * np.cumprod(1 + r_levered, axis=0)
    y_t = np.where(np.isnan(x_t), np.nan, y_t)

    return VolatilityTargetResult(y_t=y_t, sigma_t=sigma_t, leverage_t=leverage_lagged)