quanttoolbox.mixtures

mixtures.gaussian_mixture

Python alternatives

Hybrid: sklearn.mixture.GaussianMixture is more numerically robust (covariance regularization, multiple initializations, convergence diagnostics) for the general n-component EM-fitting step (estimate_em_mixture). Keep everything downstream of fitting — VaR/ES, risk contribution, risk budgeting, PDF/skewness under the mixture — since sklearn's GaussianMixture only fits parameters, nothing else.

quanttoolbox.mixtures.gaussian_mixture

Two-component Gaussian mixture models: moments, PDF, simulation, EM estimation, VaR/ES risk measures, and risk budgeting.

Ported from QuantToolBox/mixture/{mixture_moments,mixture_pdf_assets, mixture_pdf_portfolio,mixture_simulate,mixture_skewness, mixture_skewness_portfolio,mixture_univariate_thresholding, mixture_probability_filtering,mixture_compute_var,mixture_compute_es, mixture_compute_rc_var,mixture_compute_rc_es,mixture_compute_rb_var, mixture_compute_rb_es,estimate_em_mixture,logl_em_mixture}.m

Model: a random vector R is pi1-probability drawn from N(mu1, Sigma1) and pi2=1-pi1-probability drawn from N(mu2, Sigma2) -- e.g. a "normal regime" and a "stress regime" for asset returns.

Translation notes:

• mixture_compute_rb_var/mixture_compute_rb_es originally support two solver algorithms: (1) fmincon minimizing the sum of squared deviations between (risk contribution / budget) ratios across assets, and (2) an fminunc log-barrier variant referencing a RB_lagrangian global that's never actually set anywhere in the original codebase (a global exists but nothing ever assigns to it before use -- effectively dead/broken code). Only algorithm (1) is ported here, via scipy.optimize.minimize (SLSQP, budget-constrained).
• MATLAB's global MIXTURE_*/RB_* state-passing to nested objective/constraint functions is replaced by closures capturing the mixture parameters directly.

estimate_em_mixture(data, init, estimate_mixing_weights=True, tol=1e-06, max_iters=1000)

Fit a 2-component Gaussian mixture to data via Expectation-Maximization.

estimate_mixing_weights=True (default): re-estimate pi1/pi2 each iteration (standard EM). False: hold pi1/pi2 fixed at their initial values (e.g. when they're set exogenously, as in the jump-diffusion parameterization).

Original: mixture/{estimate_em_mixture,logl_em_mixture}.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def estimate_em_mixture(
    data: np.ndarray,
    init: MixtureParams,
    estimate_mixing_weights: bool = True,
    tol: float = 1e-6,
    max_iters: int = 1000,
) -> EMMixtureResult:
    """Fit a 2-component Gaussian mixture to data via Expectation-Maximization.

    estimate_mixing_weights=True (default): re-estimate pi1/pi2 each
    iteration (standard EM). False: hold pi1/pi2 fixed at their initial
    values (e.g. when they're set exogenously, as in the jump-diffusion
    parameterization).

    Original: mixture/{estimate_em_mixture,logl_em_mixture}.m
    """
    data_arr = np.asarray(data, dtype=float)
    if data_arr.ndim == 1:
        y = data_arr[:, None]
    else:
        y = data_arr
    y = y[~np.isnan(y).any(axis=1)]
    n_obs = y.shape[0]

    pi1, pi2 = init.pi1, init.pi2
    mu1, mu2 = np.asarray(init.mu1, dtype=float), np.asarray(init.mu2, dtype=float)
    sigma1, sigma2 = np.asarray(init.sigma1, dtype=float), np.asarray(init.sigma2, dtype=float)

    converged = False
    n_iter = 0
    pi1_t = np.full(n_obs, pi1)
    pi2_t = np.full(n_obs, pi2)

    for n_iter in range(1, max_iters + 1):  # noqa: B007 (used after loop)
        old = (pi1, pi2, mu1.copy(), mu2.copy(), sigma1.copy(), sigma2.copy())

        pdf1 = multivariate_normal.pdf(y, mean=mu1, cov=sigma1)
        pdf2 = multivariate_normal.pdf(y, mean=mu2, cov=sigma2)
        sum1, sum2 = pi1 * pdf1, pi2 * pdf2

        pi1_t = sum1 / (sum1 + sum2)
        pi2_t = 1 - pi1_t

        if estimate_mixing_weights:
            pi1, pi2 = float(np.mean(pi1_t)), float(np.mean(pi2_t))

        mu1 = (pi1_t[:, None] * y).sum(axis=0) / pi1_t.sum()
        mu2 = (pi2_t[:, None] * y).sum(axis=0) / pi2_t.sum()

        y1c, y2c = y - mu1, y - mu2
        sigma1 = (pi1_t[:, None, None] * (y1c[:, :, None] * y1c[:, None, :])).sum(
            axis=0
        ) / pi1_t.sum()
        sigma2 = (pi2_t[:, None, None] * (y2c[:, :, None] * y2c[:, None, :])).sum(
            axis=0
        ) / pi2_t.sum()

        diff = np.concatenate(
            [
                [pi1 - old[0], pi2 - old[1]],
                (mu1 - old[2]).flatten(),
                (mu2 - old[3]).flatten(),
                (sigma1 - old[4]).flatten(),
                (sigma2 - old[5]).flatten(),
            ]
        )
        err = float(np.max(np.abs(diff)))
        if err <= tol:
            converged = True
            break

    log_l = float(
        np.sum(
            np.log(
                pi1 * multivariate_normal.pdf(y, mean=mu1, cov=sigma1)
                + pi2 * multivariate_normal.pdf(y, mean=mu2, cov=sigma2)
            )
        )
    )

    fitted = MixtureParams(pi1=pi1, mu1=mu1, sigma1=sigma1, pi2=pi2, mu2=mu2, sigma2=sigma2)
    return EMMixtureResult(
        params=fitted, log_l=log_l, pi1_t=pi1_t, pi2_t=pi2_t, n_iters=n_iter, converged=converged
    )

mixture_compute_es(x, params, alpha)

Expected Shortfall of a portfolio x's return under the mixture distribution, at confidence level alpha.

Original: mixture/mixture_compute_es.m

Returns (ES_mixture, VaR_mixture, ES_Gaussian, VaR_Gaussian).

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_compute_es(
    x: np.ndarray, params: MixtureParams, alpha: float
) -> tuple[float, float, float, float]:
    """Expected Shortfall of a portfolio x's return under the mixture
    distribution, at confidence level alpha.

    Original: mixture/mixture_compute_es.m

    Returns (ES_mixture, VaR_mixture, ES_Gaussian, VaR_Gaussian).
    """
    x = np.asarray(x, dtype=float).flatten()
    p = params

    mu1_x = float(x @ p.mu1)
    sigma1_x = float(np.sqrt(x @ p.sigma1 @ x))
    mu2_x = float(x @ p.mu2)
    sigma2_x = float(np.sqrt(x @ p.sigma2 @ x))

    var_mixture, var_gaussian = mixture_compute_var(x, params, alpha)
    es_gaussian = -mu1_x + sigma1_x * norm.pdf(norm.ppf(alpha)) / (1 - alpha)
    es_mixture = p.pi1 * _psi_function(var_mixture, mu1_x, sigma1_x, alpha) + p.pi2 * _psi_function(
        var_mixture, mu2_x, sigma2_x, alpha
    )

    return es_mixture, var_mixture, es_gaussian, var_gaussian

mixture_compute_rb_es(params, alpha, b=None, x0=None, x_minus=0.0, x_plus=1.0)

ES risk budgeting portfolio under the mixture distribution.

Original: mixture/mixture_compute_rb_es.m (algorithm 1 -- see module docstring for the unported algorithm 2 branch)

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_compute_rb_es(
    params: MixtureParams,
    alpha: float,
    b: np.ndarray | None = None,
    x0: np.ndarray | None = None,
    x_minus: float = 0.0,
    x_plus: float = 1.0,
) -> RiskBudgetingResult:
    """ES risk budgeting portfolio under the mixture distribution.

    Original: mixture/mixture_compute_rb_es.m (algorithm 1 -- see module
    docstring for the unported algorithm 2 branch)
    """
    return _solve_mixture_rb(mixture_compute_rc_es, params, alpha, b, x0, x_minus, x_plus)

mixture_compute_rb_var(params, alpha, b=None, x0=None, x_minus=0.0, x_plus=1.0)

VaR risk budgeting portfolio under the mixture distribution: find weights x (in [x_minus, x_plus], summing to 1) whose VaR risk contributions best match the target budgets b.

Original: mixture/mixture_compute_rb_var.m (algorithm 1 -- see module docstring for the unported algorithm 2 branch)

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_compute_rb_var(
    params: MixtureParams,
    alpha: float,
    b: np.ndarray | None = None,
    x0: np.ndarray | None = None,
    x_minus: float = 0.0,
    x_plus: float = 1.0,
) -> RiskBudgetingResult:
    """VaR risk budgeting portfolio under the mixture distribution: find
    weights x (in [x_minus, x_plus], summing to 1) whose VaR risk
    contributions best match the target budgets b.

    Original: mixture/mixture_compute_rb_var.m (algorithm 1 -- see module
    docstring for the unported algorithm 2 branch)
    """
    return _solve_mixture_rb(mixture_compute_rc_var, params, alpha, b, x0, x_minus, x_plus)

mixture_compute_rc_es(x, params, alpha)

ES risk contribution decomposition under the mixture distribution.

Original: mixture/mixture_compute_rc_es.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_compute_rc_es(
    x: np.ndarray, params: MixtureParams, alpha: float
) -> RiskContributionResult:
    """ES risk contribution decomposition under the mixture distribution.

    Original: mixture/mixture_compute_rc_es.m
    """
    x = np.asarray(x, dtype=float).flatten()
    p = params
    var_result = mixture_compute_rc_var(x, params, alpha)
    var_mixture = var_result.risk

    mu_tilde = p.mu2 - p.mu1
    mu1_x = float(x @ p.mu1)
    sigma1_x = float(np.sqrt(x @ p.sigma1 @ x))
    mu2_x = float(x @ p.mu2)
    sigma2_x = float(np.sqrt(x @ p.sigma2 @ x))

    h1_x = (var_mixture + mu1_x) / sigma1_x
    h2_x = (var_mixture + mu2_x) / sigma2_x

    es = p.pi1 * _psi_function(var_mixture, mu1_x, sigma1_x, alpha) + p.pi2 * _psi_function(
        var_mixture, mu2_x, sigma2_x, alpha
    )

    w1 = p.pi1 * norm.pdf(h1_x) / sigma1_x
    w2 = p.pi2 * norm.pdf(h2_x) / sigma2_x

    delta1_x = (1 + (h1_x / sigma1_x) * var_mixture) * (p.sigma1 @ x) - var_mixture / (w1 + w2) * (
        w1 * (h1_x / sigma1_x) * (p.sigma1 @ x)
        + w2 * ((h2_x / sigma2_x) * (p.sigma2 @ x) - mu_tilde)
    )
    delta2_x = (1 + (h2_x / sigma2_x) * var_mixture) * (p.sigma2 @ x) - var_mixture / (w1 + w2) * (
        w2 * (h2_x / sigma2_x) * (p.sigma2 @ x)
        + w1 * ((h1_x / sigma1_x) * (p.sigma1 @ x) + mu_tilde)
    )

    mr = (
        w1 * delta1_x
        + w2 * delta2_x
        - (p.pi1 * p.mu1 * norm.cdf(-h1_x) + p.pi2 * p.mu2 * norm.cdf(-h2_x))
    )
    mr = mr / (1 - alpha)
    rc = x * mr
    prc = rc / np.sum(rc)

    return RiskContributionResult(
        risk=es, marginal_risk=mr, risk_contribution=rc, pct_risk_contribution=prc
    )

mixture_compute_rc_var(x, params, alpha)

VaR risk contribution decomposition under the mixture distribution.

Original: mixture/mixture_compute_rc_var.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_compute_rc_var(
    x: np.ndarray, params: MixtureParams, alpha: float
) -> RiskContributionResult:
    """VaR risk contribution decomposition under the mixture distribution.

    Original: mixture/mixture_compute_rc_var.m
    """
    x = np.asarray(x, dtype=float).flatten()
    p = params
    var_mixture, _ = mixture_compute_var(x, params, alpha)

    mu1_x = float(x @ p.mu1)
    sigma1_x = float(np.sqrt(x @ p.sigma1 @ x))
    mu2_x = float(x @ p.mu2)
    sigma2_x = float(np.sqrt(x @ p.sigma2 @ x))

    h1_x = (var_mixture + mu1_x) / sigma1_x
    h2_x = (var_mixture + mu2_x) / sigma2_x

    w1 = p.pi1 * norm.pdf(h1_x) / sigma1_x
    w2 = p.pi2 * norm.pdf(h2_x) / sigma2_x

    mr = w1 * ((h1_x / sigma1_x) * (p.sigma1 @ x) - p.mu1) + w2 * (
        (h2_x / sigma2_x) * (p.sigma2 @ x) - p.mu2
    )
    mr = mr / (w1 + w2)
    rc = x * mr
    prc = rc / np.sum(rc)

    return RiskContributionResult(
        risk=var_mixture, marginal_risk=mr, risk_contribution=rc, pct_risk_contribution=prc
    )

mixture_compute_var(x, params, alpha)

Value-at-Risk of a portfolio x's return under the mixture distribution, at confidence level alpha (found via bisection on the mixture CDF).

Original: mixture/mixture_compute_var.m

Returns (VaR_mixture, VaR_Gaussian) -- the latter using only the first (normal-regime) component, for comparison.

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_compute_var(x: np.ndarray, params: MixtureParams, alpha: float) -> tuple[float, float]:
    """Value-at-Risk of a portfolio x's return under the mixture
    distribution, at confidence level alpha (found via bisection on the
    mixture CDF).

    Original: mixture/mixture_compute_var.m

    Returns (VaR_mixture, VaR_Gaussian) -- the latter using only the
    first (normal-regime) component, for comparison.
    """
    x = np.asarray(x, dtype=float).flatten()
    p = params

    mu1_x = float(x @ p.mu1)
    sigma1_x = float(np.sqrt(x @ p.sigma1 @ x))
    mu2_x = float(x @ p.mu2)
    sigma2_x = float(np.sqrt(x @ p.sigma2 @ x))

    var_gaussian = -mu1_x + norm.ppf(alpha) * sigma1_x
    min_var, max_var = var_gaussian / 5, 5 * var_gaussian

    def objective(var: np.ndarray) -> np.ndarray:
        v = float(var)
        h1 = (v + mu1_x) / sigma1_x
        h2 = (v + mu2_x) / sigma2_x
        return np.array(p.pi1 * norm.cdf(h1) + p.pi2 * norm.cdf(h2) - alpha)

    var_mixture = float(bisection(objective, min_var, max_var))
    return var_mixture, var_gaussian

mixture_moments(params)

Mean vector and covariance matrix of the mixture distribution.

Original: mixture/mixture_moments.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_moments(params: MixtureParams) -> tuple[np.ndarray, np.ndarray]:
    """Mean vector and covariance matrix of the mixture distribution.

    Original: mixture/mixture_moments.m
    """
    p = params
    mu_bar = p.pi1 * p.mu1 + p.pi2 * p.mu2
    d_mu = p.mu1 - p.mu2
    sigma_bar = p.pi1 * p.sigma1 + p.pi2 * p.sigma2 + p.pi1 * p.pi2 * np.outer(d_mu, d_mu)
    return mu_bar, sigma_bar

mixture_pdf_assets(y, params)

Marginal PDF of each asset under the mixture distribution, evaluated at y (one column per asset, matching mu1/mu2's dimension).

Original: mixture/mixture_pdf_assets.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_pdf_assets(
    y: np.ndarray, params: MixtureParams
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Marginal PDF of each asset under the mixture distribution, evaluated
    at y (one column per asset, matching mu1/mu2's dimension).

    Original: mixture/mixture_pdf_assets.m
    """
    y = np.atleast_2d(np.asarray(y, dtype=float))
    p = params
    sigma1 = np.sqrt(np.diag(p.sigma1))
    sigma2 = np.sqrt(np.diag(p.sigma2))

    pdf1 = norm.pdf(y, loc=p.mu1, scale=sigma1)
    pdf2 = norm.pdf(y, loc=p.mu2, scale=sigma2)
    pdf = p.pi1 * pdf1 + p.pi2 * pdf2
    return pdf, pdf1, pdf2

mixture_pdf_portfolio(y, x, params)

PDF of a portfolio x's return under the mixture distribution, evaluated at y.

Original: mixture/mixture_pdf_portfolio.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_pdf_portfolio(
    y: np.ndarray, x: np.ndarray, params: MixtureParams
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """PDF of a portfolio x's return under the mixture distribution,
    evaluated at y.

    Original: mixture/mixture_pdf_portfolio.m
    """
    y = np.asarray(y, dtype=float)
    x = np.asarray(x, dtype=float).flatten()
    p = params

    mu1_x = float(x @ p.mu1)
    mu2_x = float(x @ p.mu2)
    sigma1_x = float(np.sqrt(x @ p.sigma1 @ x))
    sigma2_x = float(np.sqrt(x @ p.sigma2 @ x))

    pdf1 = norm.pdf(y, loc=mu1_x, scale=sigma1_x)
    pdf2 = norm.pdf(y, loc=mu2_x, scale=sigma2_x)
    pdf = p.pi1 * pdf1 + p.pi2 * pdf2
    return pdf, pdf1, pdf2

mixture_probability_filtering(r_t, params)

Posterior regime probabilities given an observation r_t (Bayes' rule applied to the mixture likelihood).

Original: mixture/mixture_probability_filtering.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_probability_filtering(
    r_t: np.ndarray, params: MixtureParams
) -> tuple[np.ndarray, np.ndarray]:
    """Posterior regime probabilities given an observation r_t (Bayes'
    rule applied to the mixture likelihood).

    Original: mixture/mixture_probability_filtering.m
    """
    p = params
    pdf1 = multivariate_normal.pdf(r_t, mean=p.mu1, cov=p.sigma1)
    pdf2 = multivariate_normal.pdf(r_t, mean=p.mu2, cov=p.sigma2)

    c1 = p.pi1 * pdf1
    c2 = p.pi2 * pdf2
    pi1_t = c1 / (c1 + c2)
    pi2_t = 1 - pi1_t
    return pi1_t, pi2_t

mixture_simulate(params, n_samples, rng=None)

Simulate n_samples draws from the mixture distribution.

Original: mixture/mixture_simulate.m

Returns (samples, regime) where regime[i]=1 if sample i was drawn from component 1, else 2.

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_simulate(
    params: MixtureParams, n_samples: int, rng: np.random.Generator | None = None
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate n_samples draws from the mixture distribution.

    Original: mixture/mixture_simulate.m

    Returns (samples, regime) where regime[i]=1 if sample i was drawn from
    component 1, else 2.
    """
    rng = np.random.default_rng() if rng is None else rng
    p = params
    n = np.asarray(p.mu1).shape[0]

    is_regime1 = rng.uniform(size=n_samples) <= p.pi1
    u1 = multivariate_normal.rvs(mean=p.mu1, cov=p.sigma1, size=n_samples, random_state=rng)
    u2 = multivariate_normal.rvs(mean=p.mu2, cov=p.sigma2, size=n_samples, random_state=rng)
    if n == 1:
        # scipy's rvs collapses to shape (n_samples,) for a 1-D distribution;
        # reshape to (n_samples, 1) for consistency with the n>1 case
        # (np.atleast_2d would incorrectly prepend a dimension instead).
        u1 = u1.reshape(n_samples, 1)
        u2 = u2.reshape(n_samples, 1)

    u = np.where(is_regime1[:, None], u1, u2)
    regime = np.where(is_regime1, 1, 2)
    return u, regime

mixture_skewness(pi1, mu1, sigma1, pi2, mu2, sigma2)

Mean, standard deviation, and skewness of a scalar 2-component Gaussian mixture.

Original: mixture/mixture_skewness.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_skewness(
    pi1: float, mu1: float, sigma1: float, pi2: float, mu2: float, sigma2: float
) -> tuple[float, float, float]:
    """Mean, standard deviation, and skewness of a scalar 2-component
    Gaussian mixture.

    Original: mixture/mixture_skewness.m
    """
    var1, var2 = sigma1**2, sigma2**2
    d_mu = mu1 - mu2

    mu = pi1 * mu1 + pi2 * mu2
    sigma = np.sqrt(pi1 * var1 + pi2 * var2 + pi1 * pi2 * d_mu**2)
    gamma1 = pi1 * pi2 * ((pi2 - pi1) * d_mu**3 + 3 * d_mu * (var1 - var2))
    gamma1 = gamma1 / sigma**3

    return mu, sigma, gamma1

mixture_skewness_portfolio(x, params)

Mean, standard deviation, and skewness of a portfolio x's return under the mixture distribution.

Original: mixture/mixture_skewness_portfolio.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_skewness_portfolio(x: np.ndarray, params: MixtureParams) -> tuple[float, float, float]:
    """Mean, standard deviation, and skewness of a portfolio x's return
    under the mixture distribution.

    Original: mixture/mixture_skewness_portfolio.m
    """
    x = np.asarray(x, dtype=float).flatten()
    p = params
    mu1_x = float(x @ p.mu1)
    sigma1_x = float(np.sqrt(x @ p.sigma1 @ x))
    mu2_x = float(x @ p.mu2)
    sigma2_x = float(np.sqrt(x @ p.sigma2 @ x))
    return mixture_skewness(p.pi1, mu1_x, sigma1_x, p.pi2, mu2_x, sigma2_x)

mixture_univariate_thresholding(pi1, mu1, sigma1, pi2, mu2, sigma2, pi2_star)

Find the threshold values [y_minus, y_plus] outside of which the posterior probability of being in regime 2 exceeds pi2_star (a univariate regime-classification boundary).

Original: mixture/mixture_univariate_thresholding.m

Source code in src/quanttoolbox/mixtures/gaussian_mixture.py

def mixture_univariate_thresholding(
    pi1: float, mu1: float, sigma1: float, pi2: float, mu2: float, sigma2: float, pi2_star: float
) -> tuple[float, float]:
    """Find the threshold values [y_minus, y_plus] outside of which the
    posterior probability of being in regime 2 exceeds pi2_star (a
    univariate regime-classification boundary).

    Original: mixture/mixture_univariate_thresholding.m
    """
    a = (pi2 * (1 - pi2_star) * sigma1) / (pi1 * pi2_star * sigma2)

    alpha = sigma2**2 - sigma1**2
    beta = mu2 * sigma1**2 - mu1 * sigma2**2
    gamma = mu1**2 * sigma2**2 - mu2**2 * sigma1**2 + 2 * sigma2**2 * sigma1**2 * np.log(a)

    delta = beta**2 - alpha * gamma
    y_minus = -(beta + np.sqrt(delta)) / alpha
    y_plus = (-beta + np.sqrt(delta)) / alpha
    return float(y_minus), float(y_plus)

mixtures.jump_diffusion

Python alternatives

Keep — jump-diffusion-specific risk measures with no general-purpose equivalent.

quanttoolbox.mixtures.jump_diffusion

Jump-diffusion risk measures: thin parameter-transform wrappers around the Gaussian mixture machinery in mixtures.gaussian_mixture, plus lognormal moment/skewness formulas.

Ported from QuantToolBox/mixture/{jump_compute_var,jump_compute_es, jump_compute_rc_var,jump_compute_rc_es,jump_compute_rb_var, jump_compute_rb_es,jump_pdf_assets,jump_pdf_portfolio,jump_simulate, jump_skewness,jump_skewness_portfolio,jump_univariate_thresholding, jump_probability_filtering,lognormal_moments,lognormal_skewness, bivariate_lognormal_skewness}.m

Model: over a short time step dt, returns follow a diffusion (mean mu_bar, covariance Sigma_bar) with a Poisson-arrival jump component (intensity lambda, jump mean mu_tilde, jump covariance Sigma_tilde). This is exactly a 2-component Gaussian mixture with

pi1 = 1 - lambda*dt,  mu1 = mu_bar*dt,           Sigma1 = Sigma_bar*dt
pi2 = lambda*dt,      mu2 = mu_bar*dt + mu_tilde, Sigma2 = Sigma_bar*dt + Sigma_tilde

so every jump_*.m function in the original is just this parameter transform followed by a call into the corresponding mixture_*.m function -- ported here as jump_to_mixture_params plus one-line wrappers around gaussian_mixture's functions, rather than independently re-implementing the same math.

bivariate_lognormal_skewness(mu_x, sigma_x, mu_y, sigma_y, rho)

Skewness of the sum of two correlated lognormal random variables X=exp(N(mu_x,sigma_x^2)), Y=exp(N(mu_y,sigma_y^2)) with correlation rho between the underlying normals.

Original: mixture/bivariate_lognormal_skewness.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def bivariate_lognormal_skewness(
    mu_x: float, sigma_x: float, mu_y: float, sigma_y: float, rho: float
) -> float:
    """Skewness of the sum of two correlated lognormal random variables
    X=exp(N(mu_x,sigma_x^2)), Y=exp(N(mu_y,sigma_y^2)) with correlation
    rho between the underlying normals.

    Original: mixture/bivariate_lognormal_skewness.m
    """
    sigma2_x, sigma2_y = sigma_x**2, sigma_y**2
    mu1_x = np.exp(mu_x + 0.5 * sigma2_x)
    mu1_y = np.exp(mu_y + 0.5 * sigma2_y)

    mu2_x = np.exp(2 * mu_x + sigma2_x) * (np.exp(sigma2_x) - 1)
    mu2_y = np.exp(2 * mu_y + sigma2_y) * (np.exp(sigma2_y) - 1)
    cov_xy = (
        np.exp(mu_x + mu_y + 0.5 * sigma2_x + 0.5 * sigma2_y + rho * sigma_x * sigma_y)
        - mu1_x * mu1_y
    )
    var_xy = mu2_x + mu2_y + 2 * cov_xy

    mu3_x = (np.exp(3 * sigma2_x) - 3 * np.exp(sigma2_x) + 2) * np.exp(3 * mu_x + 1.5 * sigma2_x)
    mu3_y = (np.exp(3 * sigma2_y) - 3 * np.exp(sigma2_y) + 2) * np.exp(3 * mu_y + 1.5 * sigma2_y)

    cov_xxy = (
        np.exp(2 * mu_x + sigma2_x + mu_y + 0.5 * sigma2_y)
        * (np.exp(rho * sigma_x * sigma_y) - 1)
        * (np.exp(sigma2_x + rho * sigma_x * sigma_y) + np.exp(sigma2_x) - 2)
    )
    cov_xyy = (
        np.exp(2 * mu_y + sigma2_y + mu_x + 0.5 * sigma2_x)
        * (np.exp(rho * sigma_x * sigma_y) - 1)
        * (np.exp(sigma2_y + rho * sigma_x * sigma_y) + np.exp(sigma2_y) - 2)
    )

    mu3_xy = mu3_x + mu3_y + 3 * (cov_xxy + cov_xyy)
    return float(mu3_xy / var_xy**1.5)

jump_compute_es(x, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, alpha)

Expected Shortfall under the jump-diffusion model. Original: jump_compute_es.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_compute_es(
    x: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
    alpha: float,
) -> tuple[float, float, float, float]:
    """Expected Shortfall under the jump-diffusion model. Original: jump_compute_es.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_compute_es(x, params, alpha)

jump_compute_rb_es(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, alpha, b=None, x0=None, x_minus=0.0, x_plus=1.0)

ES risk budgeting under the jump-diffusion model. Original: jump_compute_rb_es.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_compute_rb_es(
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
    alpha: float,
    b: np.ndarray | None = None,
    x0: np.ndarray | None = None,
    x_minus: float = 0.0,
    x_plus: float = 1.0,
) -> RiskBudgetingResult:
    """ES risk budgeting under the jump-diffusion model. Original: jump_compute_rb_es.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_compute_rb_es(params, alpha, b=b, x0=x0, x_minus=x_minus, x_plus=x_plus)

jump_compute_rb_var(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, alpha, b=None, x0=None, x_minus=0.0, x_plus=1.0)

VaR risk budgeting under the jump-diffusion model. Original: jump_compute_rb_var.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_compute_rb_var(
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
    alpha: float,
    b: np.ndarray | None = None,
    x0: np.ndarray | None = None,
    x_minus: float = 0.0,
    x_plus: float = 1.0,
) -> RiskBudgetingResult:
    """VaR risk budgeting under the jump-diffusion model. Original: jump_compute_rb_var.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_compute_rb_var(params, alpha, b=b, x0=x0, x_minus=x_minus, x_plus=x_plus)

jump_compute_rc_es(x, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, alpha)

ES risk contribution under the jump-diffusion model. Original: jump_compute_rc_es.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_compute_rc_es(
    x: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
    alpha: float,
) -> RiskContributionResult:
    """ES risk contribution under the jump-diffusion model. Original: jump_compute_rc_es.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_compute_rc_es(x, params, alpha)

jump_compute_rc_var(x, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, alpha)

VaR risk contribution under the jump-diffusion model. Original: jump_compute_rc_var.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_compute_rc_var(
    x: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
    alpha: float,
) -> RiskContributionResult:
    """VaR risk contribution under the jump-diffusion model. Original: jump_compute_rc_var.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_compute_rc_var(x, params, alpha)

jump_compute_var(x, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, alpha)

Value-at-Risk under the jump-diffusion model. Original: jump_compute_var.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_compute_var(
    x: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
    alpha: float,
) -> tuple[float, float]:
    """Value-at-Risk under the jump-diffusion model. Original: jump_compute_var.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_compute_var(x, params, alpha)

jump_pdf_assets(y, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)

Marginal asset PDFs under the jump-diffusion model. Original: jump_pdf_assets.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_pdf_assets(
    y: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Marginal asset PDFs under the jump-diffusion model. Original: jump_pdf_assets.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_pdf_assets(y, params)

jump_pdf_portfolio(y, x, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)

Portfolio return PDF under the jump-diffusion model. Original: jump_pdf_portfolio.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_pdf_portfolio(
    y: np.ndarray,
    x: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Portfolio return PDF under the jump-diffusion model. Original: jump_pdf_portfolio.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_pdf_portfolio(y, x, params)

jump_probability_filtering(r_t, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)

Posterior jump-regime probability given an observation. Original: jump_probability_filtering.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_probability_filtering(
    r_t: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
) -> tuple[np.ndarray, np.ndarray]:
    """Posterior jump-regime probability given an observation. Original: jump_probability_filtering.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_probability_filtering(r_t, params)

jump_simulate(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, n_samples, rng=None)

Simulate returns under the jump-diffusion model. Original: jump_simulate.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_simulate(
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
    n_samples: int,
    rng: np.random.Generator | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate returns under the jump-diffusion model. Original: jump_simulate.m"""
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_simulate(params, n_samples, rng=rng)

jump_skewness(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)

Mean/std/skewness of univariate returns under the jump-diffusion model.

Original: jump_skewness.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_skewness(
    mu_bar: float, sigma_bar: float, mu_tilde: float, sigma_tilde: float, lambda_: float, dt: float
) -> tuple[float, float, float]:
    """Mean/std/skewness of univariate returns under the jump-diffusion model.

    Original: jump_skewness.m
    """
    pi1, mu1, sigma1, pi2, mu2, sigma2 = _jump_univariate_params(
        mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt
    )
    return mixture_skewness(pi1, mu1, sigma1, pi2, mu2, sigma2)

jump_skewness_portfolio(x, mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)

Mean/std/skewness of a portfolio's return under the jump-diffusion model.

Original: jump_skewness_portfolio.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_skewness_portfolio(
    x: np.ndarray,
    mu_bar: np.ndarray,
    sigma_bar: np.ndarray,
    mu_tilde: np.ndarray,
    sigma_tilde: np.ndarray,
    lambda_: float,
    dt: float,
) -> tuple[float, float, float]:
    """Mean/std/skewness of a portfolio's return under the jump-diffusion model.

    Original: jump_skewness_portfolio.m
    """
    params = jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)
    return mixture_skewness_portfolio(x, params)

jump_to_mixture_params(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt)

Convert jump-diffusion parameters (diffusion mean/cov, jump intensity/mean/cov, time step) into the equivalent 2-component Gaussian mixture parameterization.

sigma_bar/sigma_tilde may be given as covariance matrices (for the multivariate case) or as scalar standard deviations (for the univariate skewness/thresholding helpers, which take sigma1/sigma2 directly rather than full covariance matrices).

Original: the parameter-transform preamble shared by every jump_*.m function (see module docstring)

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_to_mixture_params(
    mu_bar: np.ndarray | float,
    sigma_bar: np.ndarray | float,
    mu_tilde: np.ndarray | float,
    sigma_tilde: np.ndarray | float,
    lambda_: float,
    dt: float,
) -> MixtureParams:
    """Convert jump-diffusion parameters (diffusion mean/cov, jump
    intensity/mean/cov, time step) into the equivalent 2-component
    Gaussian mixture parameterization.

    sigma_bar/sigma_tilde may be given as covariance matrices (for the
    multivariate case) or as scalar standard deviations (for the
    univariate skewness/thresholding helpers, which take sigma1/sigma2
    directly rather than full covariance matrices).

    Original: the parameter-transform preamble shared by every
    jump_*.m function (see module docstring)
    """
    mu_bar = np.asarray(mu_bar, dtype=float) if not np.isscalar(mu_bar) else mu_bar
    sigma_bar = np.asarray(sigma_bar, dtype=float) if not np.isscalar(sigma_bar) else sigma_bar
    mu_tilde = np.asarray(mu_tilde, dtype=float) if not np.isscalar(mu_tilde) else mu_tilde
    sigma_tilde = (
        np.asarray(sigma_tilde, dtype=float) if not np.isscalar(sigma_tilde) else sigma_tilde
    )

    pi1 = 1 - lambda_ * dt
    mu1 = mu_bar * dt
    sigma1 = sigma_bar * dt
    pi2 = lambda_ * dt
    mu2 = mu_bar * dt + mu_tilde
    sigma2 = sigma_bar * dt + sigma_tilde

    return MixtureParams(pi1=pi1, mu1=mu1, sigma1=sigma1, pi2=pi2, mu2=mu2, sigma2=sigma2)

jump_univariate_thresholding(mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt, pi2_star)

Univariate jump-regime classification thresholds. Original: jump_univariate_thresholding.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def jump_univariate_thresholding(
    mu_bar: float,
    sigma_bar: float,
    mu_tilde: float,
    sigma_tilde: float,
    lambda_: float,
    dt: float,
    pi2_star: float,
) -> tuple[float, float]:
    """Univariate jump-regime classification thresholds. Original: jump_univariate_thresholding.m"""
    pi1, mu1, sigma1, pi2, mu2, sigma2 = _jump_univariate_params(
        mu_bar, sigma_bar, mu_tilde, sigma_tilde, lambda_, dt
    )
    return mixture_univariate_thresholding(pi1, mu1, sigma1, pi2, mu2, sigma2, pi2_star)

lognormal_moments(mu, sigma)

Mean, std dev, skewness, and excess kurtosis of exp(N(mu, sigma^2)) (a lognormal random variable), via its raw (uncentered) moments.

Original: mixture/lognormal_moments.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def lognormal_moments(
    mu: np.ndarray, sigma: np.ndarray
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Mean, std dev, skewness, and excess kurtosis of exp(N(mu, sigma^2))
    (a lognormal random variable), via its raw (uncentered) moments.

    Original: mixture/lognormal_moments.m
    """
    mu = np.asarray(mu, dtype=float)
    sigma = np.asarray(sigma, dtype=float)

    m = [np.exp(k * mu + (k**2) * sigma**2 / 2) for k in range(1, 5)]

    mu_x = m[0]
    var_x = m[1] - m[0] ** 2
    sigma_x = np.sqrt(var_x)

    gamma1_x = (m[2] - 3 * m[0] * m[1] + 2 * m[0] ** 3) / sigma_x**3
    gamma2_x = (m[3] - 4 * m[0] * m[2] + 6 * m[0] ** 2 * m[1] - 3 * m[0] ** 4) / sigma_x**4

    return mu_x, sigma_x, gamma1_x, gamma2_x

lognormal_skewness(mu, sigma)

Skewness of a lognormal random variable (closed form, in terms of sigma only).

Original: mixture/lognormal_skewness.m

Source code in src/quanttoolbox/mixtures/jump_diffusion.py

def lognormal_skewness(mu: np.ndarray, sigma: np.ndarray) -> np.ndarray:
    """Skewness of a lognormal random variable (closed form, in terms of
    sigma only).

    Original: mixture/lognormal_skewness.m
    """
    sigma2 = np.asarray(sigma, dtype=float) ** 2
    return (np.exp(3 * sigma2) - 3 * np.exp(sigma2) + 2) / (np.exp(sigma2) - 1) ** 1.5