quanttoolbox.linalg

Python alternatives

No standalone public utility module does this. vec/vech/xpnd/commutation/duplication/elimination matrices show up as private internals scattered inside packages like linearmodels, but nothing exposes them cleanly. Keep — genuinely fills a gap.

linalg.special_matrices

quanttoolbox.linalg.special_matrices

Special matrix operations: vec/vech/xpnd, commutation/duplication/elimination.

Ported from QuantToolBox/matrix/{vec,vech,vecr,xpnd,commutation_matrix, duplication_matrix,elimination_matrix,reshapec,reshaper,diagrv,lowmat, upmat,design}.m

Key translation notes (apply throughout this module and beyond):

• MATLAB is column-major; NumPy defaults to row-major (C order). The original vec(X) (column-major flatten) is X.flatten(order="F"). The original vecr(X) (row-major flatten, called "vec by rows" in the MATLAB code) is simply X.flatten(order="C") -- no special casing needed, it's NumPy's default.
• MATLAB is 1-indexed with inclusive slice ends; all loop bounds below are converted to 0-indexed, exclusive-end Python/NumPy equivalents.
• MATLAB "column vectors" (shape (n,1)) are represented here as flat 1-D NumPy arrays (shape (n,)), which is the idiomatic NumPy convention. Downstream ports should not assume a trailing singleton dimension.
• packr (drop NaN rows) is not ported as a standalone utility; where it was only used to strip NaN placeholders from a padded array (as in elimination_matrix), it's inlined via boolean masking.

commutation_matrix(m, n)

Commutation matrix K_(m,n) such that K @ vec(A) == vec(A.T) for an m x n A.

Original: matrix/commutation_matrix.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def commutation_matrix(m: int, n: int) -> np.ndarray:
    """Commutation matrix K_(m,n) such that K @ vec(A) == vec(A.T) for an m x n A.

    Original: matrix/commutation_matrix.m
    """
    p = m * n
    # entry (i, j), 1-indexed: i + j*m  -->  the target column of the 1 in row (i,j)
    idx_matrix = np.arange(1, m + 1)[:, None] + np.arange(0, n * m, m)[None, :]
    v = idx_matrix.flatten(order="C")  # vecr

    k = np.zeros((p, p))
    for i in range(p):
        k[i, v[i] - 1] = 1
    return k

design(v)

Build a 0/1 design (indicator) matrix from a vector of category indices.

v[i] (1-indexed category, values <= 0 produce an all-zero row) selects the column that gets a 1 in row i.

Original: matrix/design.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def design(v: np.ndarray) -> np.ndarray:
    """Build a 0/1 design (indicator) matrix from a vector of category indices.

    v[i] (1-indexed category, values <= 0 produce an all-zero row) selects
    the column that gets a 1 in row i.

    Original: matrix/design.m
    """
    v = np.asarray(v).flatten()
    n = v.shape[0]
    v_rounded = np.round(v).astype(int)
    m = int(v_rounded.max())
    x = np.zeros((n, m))
    for i in range(n):
        if v_rounded[i] > 0:
            x[i, v_rounded[i] - 1] = 1.0
    return x

diagrv(x, v)

Return a copy of x with its diagonal replaced by v.

Original: matrix/diagrv.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def diagrv(x: np.ndarray, v: np.ndarray) -> np.ndarray:
    """Return a copy of x with its diagonal replaced by v.

    Original: matrix/diagrv.m
    """
    y = np.array(x, copy=True, dtype=float)
    np.fill_diagonal(y, v)
    return y

duplication_matrix(m)

Duplication matrix D_m such that D @ vech(A) == vec(A) for symmetric A.

Original: matrix/duplication_matrix.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def duplication_matrix(m: int) -> np.ndarray:
    """Duplication matrix D_m such that D @ vech(A) == vec(A) for symmetric A.

    Original: matrix/duplication_matrix.m
    """
    p = m * (m + 1) // 2
    v = np.arange(1, p + 1)
    a = xpnd(v, method=2)
    vv = vec(a)

    d = np.zeros((m * m, p))
    for i in range(m * m):
        d[i, int(vv[i]) - 1] = 1
    return d

elimination_matrix(m)

Elimination matrix L_m such that L @ vec(A) == vech(A) for square A.

Original: matrix/elimination_matrix.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def elimination_matrix(m: int) -> np.ndarray:
    """Elimination matrix L_m such that L @ vec(A) == vech(A) for square A.

    Original: matrix/elimination_matrix.m
    """
    p = m * (m + 1) // 2
    row = np.arange(1, m + 1)
    col = np.arange(0, m * m, m)
    v = col[:, None] + row[None, :]

    shift_v = np.full((m, m), np.nan)
    for i in range(m):
        length = m - i
        shift_v[i, 0:length] = v[i, i:m]

    flat = shift_v.flatten(order="C")  # vecr
    flat = flat[~np.isnan(flat)]  # packr

    elim = np.zeros((p, m * m))
    for i in range(p):
        elim[i, int(flat[i]) - 1] = 1
    return elim

lowmat(x, v=0.0)

Return the lower-triangular part of x (including diagonal), with the strictly-upper part filled with v.

Original: matrix/lowmat.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def lowmat(x: np.ndarray, v: float = 0.0) -> np.ndarray:
    """Return the lower-triangular part of x (including diagonal), with the
    strictly-upper part filled with v.

    Original: matrix/lowmat.m
    """
    x = np.asarray(x)
    r, c = x.shape
    return np.tril(x) + np.triu(np.full((r, c), v), 1)

reshapec(v, n, m)

Reshape (and recycle/truncate as needed) a vector into an n x m matrix, column-major.

Original: matrix/reshapec.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def reshapec(v: np.ndarray, n: int, m: int) -> np.ndarray:
    """Reshape (and recycle/truncate as needed) a vector into an n x m matrix,
    column-major.

    Original: matrix/reshapec.m
    """
    v = vec(np.asarray(v))
    r = v.shape[0]
    nm = n * m
    if r > nm:
        v = v[:nm]
    elif r < nm:
        nc = int(np.ceil(nm / r))
        v = np.tile(v, nc)[:nm]
    return v.reshape((n, m), order="F")

reshaper(v, n, m)

Reshape (and recycle/truncate as needed) a vector into an n x m matrix, row-major.

Original: matrix/reshaper.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def reshaper(v: np.ndarray, n: int, m: int) -> np.ndarray:
    """Reshape (and recycle/truncate as needed) a vector into an n x m matrix,
    row-major.

    Original: matrix/reshaper.m
    """
    v = vecr(np.asarray(v))
    r = v.shape[0]
    nm = n * m
    if r > nm:
        v = v[:nm]
    elif r < nm:
        nc = int(np.ceil(nm / r))
        v = np.tile(v, nc)[:nm]
    return v.reshape((n, m), order="C")

upmat(x, v=0.0)

Return the upper-triangular part of x (including diagonal), with the strictly-lower part filled with v.

Original: matrix/upmat.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def upmat(x: np.ndarray, v: float = 0.0) -> np.ndarray:
    """Return the upper-triangular part of x (including diagonal), with the
    strictly-lower part filled with v.

    Original: matrix/upmat.m
    """
    x = np.asarray(x)
    r, c = x.shape
    return np.triu(x) + np.tril(np.full((r, c), v), -1)

vec(x)

Column-major ("Fortran order") flatten of a matrix.

Original: matrix/vec.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def vec(x: np.ndarray) -> np.ndarray:
    """Column-major ("Fortran order") flatten of a matrix.

    Original: matrix/vec.m
    """
    return np.asarray(x).flatten(order="F")

vech(a, method=1)

Half-vectorization: stack the lower triangle of a square matrix.

method=1 (default) or 'r': row-wise ordering (matches MATLAB default). method=2 or 'c': column-wise ordering (the more common convention elsewhere -- e.g. this is what most other vech implementations do).

Original: matrix/vech.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def vech(a: np.ndarray, method: int | str = 1) -> np.ndarray:
    """Half-vectorization: stack the lower triangle of a square matrix.

    method=1 (default) or 'r': row-wise ordering (matches MATLAB default).
    method=2 or 'c': column-wise ordering (the more common convention
    elsewhere -- e.g. this is what most other vech implementations do).

    Original: matrix/vech.m
    """
    if isinstance(method, str):
        method = 1 if method == "r" else 2

    a = np.asarray(a)
    r, c = a.shape
    if r != c:
        raise ValueError("vech: matrix not square")

    if method == 1:
        at = a.T
        mask = np.triu(np.ones((r, r), dtype=bool))
        return at.flatten(order="F")[mask.flatten(order="F")]
    else:
        return np.concatenate([a[i:r, i] for i in range(r)])

vecr(x)

Row-major ("C order") flatten of a matrix.

Original: matrix/vecr.m -- MATLAB computed this as vec(x'), which is algebraically identical to a native row-major flatten of x.

Source code in src/quanttoolbox/linalg/special_matrices.py

def vecr(x: np.ndarray) -> np.ndarray:
    """Row-major ("C order") flatten of a matrix.

    Original: matrix/vecr.m -- MATLAB computed this as ``vec(x')``, which
    is algebraically identical to a native row-major flatten of ``x``.
    """
    return np.asarray(x).flatten(order="C")

xpnd(v, method=1)

Inverse of vech: expand a half-vectorized vector back to a symmetric matrix.

method=1 (default) or 'r': row-wise ordering (matches MATLAB default). method=2 or 'c': column-wise ordering.

Original: matrix/xpnd.m

Source code in src/quanttoolbox/linalg/special_matrices.py

def xpnd(v: np.ndarray, method: int | str = 1) -> np.ndarray:
    """Inverse of vech: expand a half-vectorized vector back to a symmetric matrix.

    method=1 (default) or 'r': row-wise ordering (matches MATLAB default).
    method=2 or 'c': column-wise ordering.

    Original: matrix/xpnd.m
    """
    if isinstance(method, str):
        method = 1 if method == "r" else 2

    v = np.asarray(v).flatten()
    p = v.shape[0]
    n_float = (-1 + np.sqrt(1 + 8 * p)) / 2
    n = round(n_float)
    if abs(n_float - n) > 1e-4:
        raise ValueError("xpnd: the vector does not have the right dimension")

    a = np.zeros((n, n))

    if method == 1:
        for i in range(1, n + 1):
            start = i * (i - 1) // 2
            end = i * (i + 1) // 2
            x = v[start:end]
            a[i - 1, 0:i] = x
            a[0:i, i - 1] = x
    else:
        j = 0
        for i in range(n):
            length = n - i
            x = v[j : j + length]
            a[i:n, i] = x
            a[i, i:n] = x
            j += length

    return a

Examples

Black-Litterman view sensitivity across five scenarios — rpb/test_bl2.py

"""Translated from Examples/rpb/test_bl2.m -- Black-Litterman sensitivity
analysis across 6 scenarios (base case + 5 view/uncertainty/tau variants),
each solved as a fixed-risk-aversion MVO portfolio."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.black_litterman import black_litterman_moments, implied_risk_premia
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

x0 = np.array([0.40, 0.30, 0.20, 0.10])
r = 0.03
irp = implied_risk_premia(x0, cov_matrix, sharpe_ratio=0.25)
mu_tilde = r + irp.pi
gamma0 = irp.gamma

scenarios = [
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.07, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.20**2, 0.20**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=0.10,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=0.01,
    ),
]

results = [
    dict(weights=x0, mu=x0 @ mu_tilde, sigma=np.sqrt(x0 @ cov_matrix @ x0), alpha=0.0, te=0.0)
]
for s in scenarios:
    bl = black_litterman_moments(mu_tilde, s["tau"] * cov_matrix, s["P"], s["Q"], s["Omega"])
    mvo = mvo_portfolio(bl.mu_bar, cov_matrix, gamma=gamma0, lb=0.0, ub=1.0)
    alpha = (mvo.weights - x0) @ bl.mu_bar
    te = np.sqrt((mvo.weights - x0) @ cov_matrix @ (mvo.weights - x0))
    results.append(
        dict(weights=mvo.weights, mu=mvo.expected_return, sigma=mvo.volatility, alpha=alpha, te=te)
    )

for i, r_ in enumerate(results):
    print(
        f"scenario {i}: weights={np.round(r_['weights'],4)} mu={round(r_['mu'],5)} te={round(r_['te'],5)}"
    )

Elimination/duplication/commutation matrix identities across sizes — matrix/matrix2.py

"""Translated from Examples/matrix/matrix2.m -- checks seven algebraic
identities relating the elimination, duplication, and commutation
matrices, for M = 1..10.

The original compares matrix products against identity/zero matrices with
exact `==`; since these matrices are built from exact 0/1 selections the
products are exact in floating point for the equality- and sum-based
checks, but `inv(D'*D)` (a genuine matrix inversion) is not bit-exact, so
`np.isclose`/`np.allclose` are used throughout instead of `==` to make the
comparison robust the way an equivalent MATLAB run using a tolerance
would be."""

import numpy as np

from quanttoolbox.linalg.special_matrices import (
    commutation_matrix,
    duplication_matrix,
    elimination_matrix,
)

for M in range(1, 11):
    p = M * (M + 1) // 2

    L = elimination_matrix(M)
    D = duplication_matrix(M)
    K = commutation_matrix(M, M)
    K1 = commutation_matrix(M, 1)
    K2 = commutation_matrix(1, M)

    checks = [
        np.allclose(L @ D, np.eye(p)),
        np.allclose(K @ D, D),
        np.array_equal(K1 == K2, K1 == np.eye(M)),
        np.isclose(np.sum(np.diag(K)), M),
        np.isclose(np.sum(np.diag(D.T @ D)), M**2),
        np.allclose(L @ L.T, np.eye(p)),
        np.isclose(np.sum(np.diag(np.linalg.inv(D.T @ D))), M * (M + 3) / 4),
    ]

    verdict = (
        "The propositions are verified" if all(checks) else "The propositions are NOT verified"
    )
    print(f"M={M}: {verdict}")

Equal risk contribution and box-constrained risk budgeting — rpb/test_box1.py

"""Translated from Examples/rpb/test_box1.m -- ERC and box-constrained
("C-ERC") risk budgeting portfolios at progressively wider bounds around
a starting position."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import erc_portfolio, solve_box_constrained
from quanttoolbox.stats.moments import corr_to_cov

x0 = np.array([0.29, 0.25, 0.23, 0.18, 0.05])
sigma = np.array([0.20, 0.20, 0.25, 0.15, 0.25])
rho = xpnd(
    np.array(
        [1.00, 0.40, 1.00, 0.70, 0.75, 1.00, 0.60, 0.55, 0.90, 1.00, 0.70, 0.60, 0.70, 0.65, 1.00]
    ),
    method=1,
)
cov_matrix = corr_to_cov(sigma, rho)

r1 = erc_portfolio(cov_matrix)
print("ERC weights:", np.round(r1.weights, 4))

for delta in [0.02, 0.07, 0.20]:
    x_minus, x_plus = x0 - delta, x0 + delta
    r = solve_box_constrained(cov_matrix, x_minus=x_minus, x_plus=x_plus, x0=x0)
    print(f"box (delta={delta}) weights:", np.round(r.weights, 4), "converged:", r.converged)

Explicit/implicit constraint round-trip, plus a design() demo — optim/explicit2.py

"""Translated from Examples/optim/explicit2.m -- explicit<->implicit
constraint conversion for a single equality constraint (x[0] = x[1]), plus
a `design` matrix built from a category-index vector."""

import numpy as np

from quanttoolbox.linalg.special_matrices import design
from quanttoolbox.optim.bisection import explicit_to_implicit, implicit_to_explicit

CC = np.array([[1.0, -1.0, 0, 0, 0, 0, 0, 0]])  # constraint: x[0] - x[1] = 0
c = np.array([0.0])

RR, r = explicit_to_implicit(CC, c)
CC2, c2 = implicit_to_explicit(RR, r)

print("R:")
print(RR)
print("\nr:", r)
print("\nC:")
print(CC2)
print("\nc:", c2)

w = np.concatenate([[1.0, 1.0], np.arange(2, 8, dtype=float)])  # seqa(2,1,6)
RR2 = design(w)
r2 = np.zeros(8)

print("\nR (design(w)):")
print(RR2)
print("\nr:", r2)

Larger restricted SUR system, two covariance variants compared — ects/varx4b.py

"""Translated from Examples/ects/varx4b.m -- Judge, Hill, Griffiths,
Lutkepohl & Lee [1988], pages 460-462: restricted SUR estimation of a
3-equation system (log income on log prices + log quantity, 15 stacked
coefficients restricted to 7/9 free parameters) on a 30-observation
embedded dataset, comparing two different restricted-residual Sigma
estimates used as the GLS weighting matrix for the final restricted SUR
fit.

Three steps, exactly as the original: (1) Sigma1 from a 9-free-parameter
restriction (`w`), (2) Sigma2 from a more tightly restricted 7-free-
parameter version of the same `w` (indices 8 and 12 tied to the same
free parameter as index 4), (3) the final restricted SUR fit run once
with each Sigma."""

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate
from quanttoolbox.linalg.special_matrices import design

data = np.array(
    [
        [10.763, 4.474, 6.629, 487.648, 11.632, 13.194, 45.770],
        [13.033, 10.836, 13.774, 364.877, 12.029, 2.181, 13.393],
        [9.244, 5.856, 4.063, 514.037, 8.196, 5.586, 104.819],
        [4.605, 14.010, 3.868, 760.343, 33.908, 5.231, 137.269],
        [13.045, 11.417, 14.922, 421.746, 4.561, 10.930, 15.914],
        [7.706, 8.755, 14.138, 578.214, 17.594, 11.854, 23.667],
        [7.405, 7.317, 4.794, 561.734, 18.842, 17.045, 62.057],
        [7.519, 6.360, 3.768, 301.470, 11.637, 2.682, 52.262],
        [8.764, 4.188, 8.089, 379.636, 7.645, 13.008, 31.916],
        [13.511, 1.996, 2.708, 478.855, 7.881, 19.623, 123.026],
        [4.943, 7.268, 12.901, 433.741, 9.614, 6.534, 26.255],
        [8.360, 5.839, 11.115, 525.702, 9.067, 9.397, 35.540],
        [5.721, 5.160, 11.220, 513.067, 14.070, 13.188, 32.487],
        [7.225, 9.145, 5.810, 408.666, 15.474, 3.340, 45.838],
        [6.617, 5.034, 5.516, 192.061, 3.041, 4.716, 26.867],
        [14.219, 5.926, 3.707, 462.621, 14.096, 17.141, 43.325],
        [6.769, 8.187, 10.125, 312.659, 4.118, 4.695, 24.330],
        [7.769, 7.193, 2.471, 400.848, 10.489, 7.639, 107.017],
        [9.804, 13.315, 8.976, 392.215, 6.231, 9.089, 23.407],
        [11.063, 6.874, 12.883, 377.724, 6.458, 10.346, 18.254],
        [6.535, 15.533, 4.115, 343.552, 8.736, 3.901, 54.895],
        [11.063, 4.477, 4.962, 301.599, 5.158, 4.350, 45.360],
        [4.016, 9.231, 6.294, 294.112, 16.618, 7.371, 25.318],
        [4.759, 5.907, 8.298, 365.032, 11.342, 6.507, 32.852],
        [5.483, 7.077, 9.638, 256.125, 2.903, 3.770, 22.154],
        [7.890, 9.942, 7.122, 184.798, 3.138, 1.360, 20.575],
        [8.460, 7.043, 4.157, 359.084, 15.315, 6.497, 44.205],
        [6.195, 4.142, 10.040, 629.378, 22.240, 10.963, 44.443],
        [6.743, 3.369, 15.459, 306.527, 10.012, 10.140, 13.251],
        [11.977, 4.806, 6.172, 347.488, 3.982, 8.637, 41.845],
    ]
)
data = np.log(data)

p = data[:, [0, 1, 2]]  # Prices
y = data[:, 3]  # Quantities
q = data[:, [4, 5, 6]]  # Income

big_y = q
big_x = np.column_stack([np.ones(30), p, y])

# First, estimate Sigma from the unrestricted least-squares residuals
w1 = np.array([1, 2, 3, 4, 0, 0, 0, 5, 0, 0, 0, 6, 7, 8, 9])
rr1 = design(w1)
r1 = np.zeros(15)
result1 = varx_estimate(big_y, big_x, p=0, restriction=(rr1, r1), method="ls")
sigma1 = result1.sigma

# Or estimate Sigma from the (more tightly) restricted least-squares
# residuals (indices 8 and 12 tied to free parameter 4)
w2 = np.array([1, 2, 3, 4, 0, 0, 0, 4, 0, 0, 0, 4, 5, 6, 7])
rr2 = design(w2)
r2 = np.zeros(15)
result2 = varx_estimate(big_y, big_x, p=0, restriction=(rr2, r2), method="ls")
sigma2 = result2.sigma

# Then, perform the restricted SUR estimation with each Sigma
sur1 = varx_estimate(big_y, big_x, p=0, restriction=(rr2, r2), sigma=sigma1, method="ls")
print("SUR (Sigma1) theta:", np.round(sur1.theta, 5))

sur2 = varx_estimate(big_y, big_x, p=0, restriction=(rr2, r2), sigma=sigma2, method="ls")
print("\nSUR (Sigma2) theta:", np.round(sur2.theta, 5))

Mean-variance frontier: risk-aversion, return, and volatility targets — rpb/test_mvo2.py

"""Translated from Examples/rpb/test_mvo2.m -- Roncalli [2013], "Introduction
to Risk Parity and Budgeting", Example 1 (pages 7-8): the same 4-asset
mean-variance problem evaluated three ways -- the gamma-problem (pick a
risk-aversion, solve directly), the mu-problem (pick a target expected
return, bisect on gamma to hit it), and the sigma-problem (pick a target
volatility, bisect on gamma to hit it). All three route through
`compute_mvo_portfolio.m`'s three branches; here that's
`mvo_frontier`/`mvo_target_portfolio`.

The original passes `lb=0, ub=0` (MATLAB's "use the default -100/100 wide
bounds" sentinel); passed through explicitly here as `lb=-100, ub=100`."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.mean_variance import mvo_frontier, mvo_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

print("1. gamma-problem (page 7)")
gamma_values = np.array([0.00, 0.20, 0.50, 1.00, 2.00, 5.00])
results = mvo_frontier(mu, cov_matrix, gamma_values, lb=-100.0, ub=100.0)
for g, r in zip(gamma_values, results, strict=False):
    print(
        f"  gamma={g:5.2f}  mu={100 * r.expected_return:6.2f}  sigma={100 * r.volatility:6.2f}  "
        f"w={np.round(100 * r.weights, 2)}"
    )

print("\n2. mu-problem (page 8)")
mu_targets = np.array([5.00, 6.00, 7.00, 8.00, 9.00]) / 100
mu_results = mvo_target_portfolio(mu, cov_matrix, mu_targets, problem="mu", lb=-100.0, ub=100.0)
for target, r in zip(mu_targets, mu_results, strict=False):
    print(
        f"  target_mu={100 * target:5.2f}  gamma={r.gamma:6.3f}  mu={100 * r.expected_return:6.2f}  "
        f"sigma={100 * r.volatility:6.2f}  w={np.round(100 * r.weights, 2)}"
    )

print("\n3. sigma-problem (page 8)")
sigma_targets = np.array([15.00, 20.00, 25.00, 30.00, 35.00]) / 100
sigma_results = mvo_target_portfolio(
    mu, cov_matrix, sigma_targets, problem="sigma", lb=-100.0, ub=100.0
)
for target, r in zip(sigma_targets, sigma_results, strict=False):
    print(
        f"  target_sigma={100 * target:5.2f}  gamma={r.gamma:6.3f}  mu={100 * r.expected_return:6.2f}  "
        f"sigma={100 * r.volatility:6.2f}  w={np.round(100 * r.weights, 2)}"
    )

Mean-variance optimization plus ridge/lasso-penalized portfolios — rpb/test_lasso1.py

"""Translated from Examples/rpb/test_lasso1.m -- Roncalli [2013],
"Introduction to Risk Parity and Budgeting", Example 1 (page 53):
compares an unconstrained gamma-problem MVO portfolio against the same
problem with a budget (sum-to-1) constraint, and ridge/lasso-penalized
variants (toward 0 and toward equal-weight) of the unconstrained
problem.

`compute_mvo_portfolio`/`quadprog_ridge`/`quadprog_lasso` map to
`mvo_portfolio`/`solve_qp(..., ridge_penalty=...)`/`solve_qp(...,
lasso_penalty=...)` as established in test_lasso3.py/test_lasso5.py; the
original's `0, 0` sentinel arguments for "no equality constraint" /
"no bounds" map to simply omitting those keyword arguments."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)

# Case gamma-problem
gamma_x = 0.5

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x).weights
x2 = mvo_portfolio(
    mu, cov_matrix, gamma=gamma_x, a_eq=np.ones((1, n)), b_eq=np.array([1.0])
).weights

lambda_ridge = 0.03
s_ridge = lambda_ridge * np.eye(n)
x3 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, np.zeros(n)))
x4 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, x0))

lambda_lasso = 0.03 / 2
s_lasso = lambda_lasso * np.ones(n)
x5 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, np.zeros(n)))
x6 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, x0))

results = 100 * np.column_stack([x1, x2, x3, x4, x5, x6])
print("            x1      x2      x3      x4      x5      x6")
print(np.round(results, 2))

Minimum-variance portfolio under general linear constraints — rpb/test_minvar2.py

"""Translated from Examples/rpb/test_minvar2.m -- minimum-variance
portfolio with general linear equality/inequality constraints and box
bounds."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.mean_variance import minvar_portfolio
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

a_eq = np.array([[1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 0.0, 0.0]])
b_eq = np.array([1.0, 0.0])
c_ineq = np.array([[0.0, 0.0, 0.0, -1.0]])
d_ineq = np.array([-0.90])

r = minvar_portfolio(
    cov_matrix, a_eq=a_eq, b_eq=b_eq, c_ineq=c_ineq, d_ineq=d_ineq, lb=-1.50, ub=2.00
)
print("weights:", np.round(r.weights, 3))
print("volatility:", round(r.volatility, 5))

Mixed ridge+lasso penalties toward two different target vectors — rpb/test_lasso5.py

"""Translated from Examples/rpb/test_lasso5.m -- ridge/lasso/mixed
portfolios with two different ridge/lasso target vectors (equal-weight
and a custom 20/20/30/30 target)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
a_eq, b_eq = np.ones((1, n)), np.array([1.0])
gamma_x = 0.5

y1 = np.full(n, 1 / n)
y2 = np.array([0.20, 0.20, 0.30, 0.30])
S_ridge = np.diag(np.diag(cov_matrix))
lambda_lasso = 0.005

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0).weights
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, y1)
)
x4 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, y1)
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y1),
)
# mixed with DIFFERENT targets for ridge (toward y1) vs lasso (toward y2)
x8 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y2),
)

for name, x in [
    ("MVO", x1),
    ("Ridge->y1", x2),
    ("Lasso->y1", x4),
    ("Mixed(ridge->y1,lasso->y1)", x6),
    ("Mixed(ridge->y1,lasso->y2)", x8),
]:
    print(f"{name}: {np.round(x, 4)}")

Principal component analysis of a 3-asset correlation matrix — stats/pca1.py

"""Translated from Examples/stats/pca1.m -- PCA on a 3-asset correlation
matrix."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.stats.regression.ols import pca

C = xpnd(np.array([1.00, 0.80, 1.00, 0.80, 0.80, 1.00]), method=1)
result = pca(C)
print("eigenvalues:", np.round(result.eigenvalues, 4))
print("quality (variance share):", np.round(result.quality, 4))
print("loadings:\n", np.round(result.loadings, 4))

Restricted seemingly-unrelated-regressions system, two-step GLS — ects/varx4a.py

"""Translated from Examples/ects/varx4a.m -- Judge, Hill, Griffiths,
Lutkepohl & Lee [1988], "Introduction to the Theory and Practice of
Econometrics", pages 444-455: restricted SUR (seemingly-unrelated-
regressions) estimation of a 2-equation system with 10 stacked
coefficients restricted down to 5 free parameters (via `design(w)`), on
a small 20-observation embedded dataset.

Two-step procedure, exactly as the original: first estimate with an
identity residual covariance (`sigma=None`) to obtain a residual
covariance estimate, then re-estimate using that Sigma as the GLS
weighting matrix (`varx_cls` -> `varx_estimate(..., method="ls")` in
both steps; `p=0` means no autoregressive lags, i.e. this is a pure
restricted SUR regression via the VARX machinery, `results.B` here is
just `result.beta` since `p=0` makes `result.phi` empty)."""

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate
from quanttoolbox.linalg.special_matrices import design

data = np.array(
    [
        [40.05292, 1170.6, 97.8, 2.52813, 191.5, 1.8],
        [54.64859, 2015.8, 104.4, 24.91888, 516, 0.8],
        [40.31206, 2803.3, 118, 29.3427, 729, 7.4],
        [84.21099, 2039.7, 156.2, 27.61823, 560.4, 18.1],
        [127.5724, 2256.2, 172.6, 60.35945, 519.9, 23.5],
        [124.8797, 2132.2, 186.6, 50.61588, 628.5, 26.5],
        [96.55514, 1834.1, 220.9, 30.70955, 537.1, 36.2],
        [131.1601, 1588, 287.8, 60.69605, 561.2, 60.8],
        [77.02764, 1749.4, 319.9, 30.00972, 617.2, 84.4],
        [46.96689, 1687.2, 321.3, 42.5075, 626.7, 91.2],
        [100.6597, 2007.7, 319.6, 58.61146, 737.2, 92.4],
        [115.7467, 2208.3, 346, 46.96287, 760.5, 86],
        [114.5826, 1656.7, 456.4, 57.87651, 581.4, 111.1],
        [119.8762, 1604.4, 543.4, 43.22093, 662.3, 130.6],
        [105.5699, 1431.8, 618.3, 22.87143, 583.8, 141.8],
        [148.4266, 1610.5, 647.4, 52.94754, 635.2, 136.7],
        [194.3622, 1819.4, 671.3, 71.2303, 723.8, 129.7],
        [158.2037, 2079.7, 726.1, 61.7255, 864.1, 145.5],
        [163.093, 2371.6, 800.3, 85.13053, 1193.5, 174.8],
        [227.5634, 2759.9, 888.9, 88.27518, 1188.9, 213.5],
    ]
)

y = data[:, [0, 3]]
x = np.column_stack([np.ones(20), data[:, [1, 2, 4, 5]]])

w = np.array([1, 2, 3, 0, 4, 0, 0, 5, 0, 6])
rr = design(w)
r = np.zeros(10)

# First perform a VARX estimation to obtain an estimate of Sigma
result0 = varx_estimate(y, x, p=0, restriction=(rr, r), method="ls")
sigma = result0.sigma

# Then, perform a VARX estimation given the estimated Sigma
result = varx_estimate(y, x, p=0, restriction=(rr, r), sigma=sigma, method="ls")

print("theta:", np.round(result.theta, 5))
print("stderr:", np.round(result.stderr, 5))
print("\nB (= beta, since p=0):")
print(np.round(result.beta, 5))

Restricted VAR via concentrated ML, vs. restricted least squares — ects/varx2c.py

"""Translated from Examples/ects/varx2c.m -- identical setup to
varx2b.py, but estimated via `varx_cml` -> `varx_estimate_cml`
(concentrated/iterated ML) instead of restricted LS."""

import io

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate_cml
from quanttoolbox.linalg.special_matrices import design

_LUTKEPOHL_DATA = """
       601       180       451       415
       602       179       465       421
       603       185       485       434
       604       192       493       448
       611       211       509       459
       612       202       520       458
       613       207       521       479
       614       214       540       487
       621       231       548       497
       622       229       558       510
       623       234       574       516
       624       237       583       525
       631       206       591       529
       632       250       599       538
       633       259       610       546
       634       263       627       555
       641       264       642       574
       642       280       653       574
       643       282       660       586
       644       292       694       602
       651       286       709       617
       652       302       734       639
       653       304       751       653
       654       307       763       668
       661       317       766       679
       662       314       779       686
       663       306       808       697
       664       304       785       688
       671       292       794       704
       672       275       799       699
       673       273       799       709
       674       301       812       715
       681       280       837       724
       682       289       853       746
       683       303       876       758
       684       322       897       779
       691       315       922       798
       692       339       949       816
       693       364       979       837
       694       371       988       858
       701       375      1025       881
       702       432      1063       905
       703       453      1104       934
       704       460      1131       968
       711       475      1137       983
       712       496      1178      1013
       713       494      1211      1034
       714       498      1256      1064
       721       526      1290      1101
       722       519      1314      1102
       723       516      1346      1145
       724       531      1385      1173
       731       573      1416      1216
       732       551      1436      1229
       733       538      1462      1242
       734       532      1493      1267
       741       558      1516      1295
       742       524      1557      1317
       743       525      1613      1355
       744       519      1642      1371
       751       526      1690      1402
       752       510      1759      1452
       753       519      1756      1485
       754       538      1780      1516
       761       549      1807      1549
       762       570      1831      1567
       763       559      1873      1588
       764       584      1897      1631
       771       611      1910      1650
       772       597      1943      1685
       773       603      1976      1722
       774       619      2018      1752
       781       635      2040      1774
       782       658      2070      1807
       783       675      2121      1831
       784       700      2132      1842
"""

data = np.loadtxt(io.StringIO(_LUTKEPOHL_DATA))
data = np.log(data)

investment = data[:, 1]
income = data[:, 2]
consumption = data[:, 3]

n = data.shape[0]
lag_investment = np.full(n, np.nan)
lag_investment[1:] = investment[:-1]

x = np.column_stack([np.ones(n), lag_investment])
y = np.column_stack([income, consumption])

w = np.array([0, 1, 0, 2, 0, 3, 4, 0])
rr = design(w)
r = np.array([1.0, 0, 0, 0, 0, 0, 0, 0])

result = varx_estimate_cml(y, x, p=1, restriction=(rr, r))
b = np.hstack([result.phi, result.beta])

print("               Inc(t-1)   Cons(t-1)    Constant    Inv(t-1)")
print("Inc(t)  ", np.round(b[0], 4))
print("Cons(t) ", np.round(b[1], 4))

Restricted VAR(2) via least squares and concentrated ML — ects/varx1d.py

"""Translated from Examples/ects/varx1d.m -- Lutkepohl [1991], chapter 5:
re-estimates the VAR(2) model from varx1a.py with 14 of its 21
coefficients restricted to exactly 0 (the ones that came out
non-significant in the unrestricted fit), leaving only 7 free
parameters, via both the LS (`varx_cls` -> `varx_estimate(...,
method="ls")`) and concentrated/iterated-ML (`varx_cml` ->
`varx_estimate_cml`) restricted estimators, and compares the two
resulting B = [Phi, beta] matrices.

`RR = design(w)` builds the (21, 7) restriction matrix directly from the
free-parameter index vector `w` (0 = restricted to r=0, k = maps to free
parameter k), exactly matching `quanttoolbox`'s `restriction=(RR, r)`
convention."""

import io

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate, varx_estimate_cml
from quanttoolbox.linalg.special_matrices import design

_LUTKEPOHL_DATA = """
       601       180       451       415
       602       179       465       421
       603       185       485       434
       604       192       493       448
       611       211       509       459
       612       202       520       458
       613       207       521       479
       614       214       540       487
       621       231       548       497
       622       229       558       510
       623       234       574       516
       624       237       583       525
       631       206       591       529
       632       250       599       538
       633       259       610       546
       634       263       627       555
       641       264       642       574
       642       280       653       574
       643       282       660       586
       644       292       694       602
       651       286       709       617
       652       302       734       639
       653       304       751       653
       654       307       763       668
       661       317       766       679
       662       314       779       686
       663       306       808       697
       664       304       785       688
       671       292       794       704
       672       275       799       699
       673       273       799       709
       674       301       812       715
       681       280       837       724
       682       289       853       746
       683       303       876       758
       684       322       897       779
       691       315       922       798
       692       339       949       816
       693       364       979       837
       694       371       988       858
       701       375      1025       881
       702       432      1063       905
       703       453      1104       934
       704       460      1131       968
       711       475      1137       983
       712       496      1178      1013
       713       494      1211      1034
       714       498      1256      1064
       721       526      1290      1101
       722       519      1314      1102
       723       516      1346      1145
       724       531      1385      1173
       731       573      1416      1216
       732       551      1436      1229
       733       538      1462      1242
       734       532      1493      1267
       741       558      1516      1295
       742       524      1557      1317
       743       525      1613      1355
       744       519      1642      1371
       751       526      1690      1402
       752       510      1759      1452
       753       519      1756      1485
       754       538      1780      1516
       761       549      1807      1549
       762       570      1831      1567
       763       559      1873      1588
       764       584      1897      1631
       771       611      1910      1650
       772       597      1943      1685
       773       603      1976      1722
       774       619      2018      1752
       781       635      2040      1774
       782       658      2070      1807
       783       675      2121      1831
       784       700      2132      1842
"""

data = np.loadtxt(io.StringIO(_LUTKEPOHL_DATA))
data = np.log(data)

investment = data[:, 1]
income = data[:, 2]
consumption = data[:, 3]


def _diff1(a: np.ndarray) -> np.ndarray:
    d = np.full_like(a, np.nan)
    d[1:] = a[1:] - a[:-1]
    return d


y = np.column_stack([_diff1(investment), _diff1(income), _diff1(consumption)])
n = y.shape[0]
x = np.ones((n, 1))

# Use the results of the VAR(2) estimation and set non-significant
# estimates equal to 0
w = np.array([1, 0, 0, 0, 0, 2, 0, 3, 4, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 6, 7])
rr = design(w)
r = np.zeros(21)

result1 = varx_estimate(y, x, p=2, restriction=(rr, r), method="ls")
result2 = varx_estimate_cml(y, x, p=2, restriction=(rr, r))

b1 = np.hstack([result1.phi, result1.beta])
b2 = np.hstack([result2.phi, result2.beta])

print("B(LS) =")
print(np.round(b1, 4))

print("\nB(ML) =")
print(np.round(b2, 4))

Restricted-coefficient VAR via restricted least squares — ects/varx2b.py

"""Translated from Examples/ects/varx2b.m -- Lutkepohl [1991], chapter
10: re-estimates the reduced-form dynamic simultaneous-equations model
from varx2a.py with four coefficients restricted: the Inc(t-1)
coefficient in the Inc(t) equation fixed at 1, and the Cons(t-1)/
constant coefficients in the Inc(t) equation and the Inv(t-1)
coefficient in the Cons(t) equation all fixed at 0 -- via `varx_cls` ->
`varx_estimate(..., method="ls")`."""

import io

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate
from quanttoolbox.linalg.special_matrices import design

_LUTKEPOHL_DATA = """
       601       180       451       415
       602       179       465       421
       603       185       485       434
       604       192       493       448
       611       211       509       459
       612       202       520       458
       613       207       521       479
       614       214       540       487
       621       231       548       497
       622       229       558       510
       623       234       574       516
       624       237       583       525
       631       206       591       529
       632       250       599       538
       633       259       610       546
       634       263       627       555
       641       264       642       574
       642       280       653       574
       643       282       660       586
       644       292       694       602
       651       286       709       617
       652       302       734       639
       653       304       751       653
       654       307       763       668
       661       317       766       679
       662       314       779       686
       663       306       808       697
       664       304       785       688
       671       292       794       704
       672       275       799       699
       673       273       799       709
       674       301       812       715
       681       280       837       724
       682       289       853       746
       683       303       876       758
       684       322       897       779
       691       315       922       798
       692       339       949       816
       693       364       979       837
       694       371       988       858
       701       375      1025       881
       702       432      1063       905
       703       453      1104       934
       704       460      1131       968
       711       475      1137       983
       712       496      1178      1013
       713       494      1211      1034
       714       498      1256      1064
       721       526      1290      1101
       722       519      1314      1102
       723       516      1346      1145
       724       531      1385      1173
       731       573      1416      1216
       732       551      1436      1229
       733       538      1462      1242
       734       532      1493      1267
       741       558      1516      1295
       742       524      1557      1317
       743       525      1613      1355
       744       519      1642      1371
       751       526      1690      1402
       752       510      1759      1452
       753       519      1756      1485
       754       538      1780      1516
       761       549      1807      1549
       762       570      1831      1567
       763       559      1873      1588
       764       584      1897      1631
       771       611      1910      1650
       772       597      1943      1685
       773       603      1976      1722
       774       619      2018      1752
       781       635      2040      1774
       782       658      2070      1807
       783       675      2121      1831
       784       700      2132      1842
"""

data = np.loadtxt(io.StringIO(_LUTKEPOHL_DATA))
data = np.log(data)

investment = data[:, 1]
income = data[:, 2]
consumption = data[:, 3]

n = data.shape[0]
lag_investment = np.full(n, np.nan)
lag_investment[1:] = investment[:-1]

x = np.column_stack([np.ones(n), lag_investment])
y = np.column_stack([income, consumption])

# Restrictions: Inc(t-1) coeff in Inc(t) eq = 1; Cons(t-1) coeff and
# constant in Inc(t) eq = 0; Inv(t-1) coeff in Cons(t) eq = 0.
w = np.array([0, 1, 0, 2, 0, 3, 4, 0])
rr = design(w)
r = np.array([1.0, 0, 0, 0, 0, 0, 0, 0])

result = varx_estimate(y, x, p=1, restriction=(rr, r), method="ls")
b = np.hstack([result.phi, result.beta])

print("               Inc(t-1)   Cons(t-1)    Constant    Inv(t-1)")
print("Inc(t)  ", np.round(b[0], 4))
print("Cons(t) ", np.round(b[1], 4))

Ridge, lasso, and mixed-norm penalized portfolios — rpb/test_lasso3.py

"""Translated from Examples/rpb/test_lasso3.m (byte-identical to
test_lasso4.m) -- MVO/ridge/lasso/mixed-penalty portfolios via the
consolidated solve_qp (replacing the original's separate quadprog_ridge/
quadprog_lasso/quadprog_mixed calls -- see optim/quadprog.py docstring)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)
gamma_x = 0.5
a_eq, b_eq = np.ones((1, n)), np.array([1.0])

# gamma-problem (plain MVO)
r1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0)
x1 = r1.weights

# ridge-problem: penalize toward zero, scaled by each asset's own variance
S_ridge = np.diag(np.diag(cov_matrix))
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, x0)
)

# lasso-problem: L1 penalty toward equal weight
lambda_lasso = 0.005
x3 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, x0)
)

# mixed-problem: both ridge and lasso, toward various targets
x4 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x5 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, x0),
)
x7 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, x0),
)

for name, x in [
    ("MVO", x1),
    ("Ridge", x2),
    ("Lasso", x3),
    ("Mixed(0,0)", x4),
    ("Mixed(EW,0)", x5),
    ("Mixed(0,EW)", x6),
    ("Mixed(EW,EW)", x7),
]:
    print(f"{name}: {np.round(x, 4)}")

Risk budgeting toward unequal target budgets — rpb/test_erc3.py

"""Translated from Examples/rpb/test_erc3.m -- Example 17 (page 123),
Roncalli (2013). RB portfolio with unequal target budgets."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import risk_contribution, solve_unconstrained
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.30, 0.10])
rho = xpnd(np.array([1.00, 0.50, 1.00, 0.00, 0.20, 1.00, -0.10, 0.40, 0.70, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
x = np.full(4, 0.25)

rc = risk_contribution(x, cov_matrix)
print("equal-weight risk contribution:", rc.risk, np.round(100 * rc.pct_risk_contribution, 2))

b = np.array([0.20, 0.20, 0.30, 0.30])
r = solve_unconstrained(cov_matrix, b=b, method="ccd")
print("RB weights (target budgets 20/20/30/30):", np.round(r.weights, 4))
print("converged:", r.converged, "n_iters:", r.n_iters)

Risk contribution decomposition and equal-budget risk budgeting — rpb/test_erc2.py

"""Translated from Examples/rpb/test_erc2.m -- Example 7 (page 80),
Roncalli (2013). Risk contribution decomposition + risk budgeting at
various target budgets."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import risk_contribution, solve_unconstrained
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.30, 0.20, 0.15])
rho = xpnd(np.array([1.00, 0.80, 1.00, 0.50, 0.30, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
x = np.array([0.50, 0.20, 0.30])

rc = risk_contribution(x, cov_matrix)
print("fixed-x risk contribution:", rc.risk, np.round(100 * rc.pct_risk_contribution, 2))

r_equal = solve_unconstrained(cov_matrix, b=np.full(3, 1 / 3), method="ccd")
print("equal-budget RB weights:", np.round(r_equal.weights, 4))

r_custom = solve_unconstrained(cov_matrix, b=x, method="ccd")
print("target-x-as-budget RB weights:", np.round(r_custom.weights, 4))

Simultaneous-equations system: OLS, GLS, and LIML compared — ects/varx5a.py

"""Translated from Examples/ects/varx5a.m -- Judge, Hill, Griffiths,
Lutkepohl & Lee [1988], pages 636-663: a 3-equation simultaneous linear
model estimated four ways -- OLS, GLS-2S, GLS-3S, LIML -- on a
20-observation embedded dataset, recovering the structural Gamma/B
matrices and the implied reduced form PI = -B @ inv(Gamma) at each step.

The original manually iterates feasible GLS across steps 1-3 (each
`varx_cls` call re-estimates Sigma, then feeds it into the next call as
the GLS weight) before switching to `varx_cml` (-> `varx_estimate_cml`,
concentrated/iterated ML) for the LIML step, seeded from step 3's Sigma
via `VARX_Tol = 1e-20`. `varx_estimate_cml` here always starts its own
internal iteration from the identity matrix rather than accepting a
seed Sigma -- since concentrated ML converges to the same fixed point
regardless of starting Sigma (only the number of iterations to get there
differs), this is translated as a plain `varx_estimate_cml(..., tol=
1e-20)` call rather than trying to seed it.

`theta(1:9)` reshaped (row-major, truncating/ignoring the rest) into
Gamma, and `theta(10:24)` reshaped into B, are both translated via the
ported `reshaper` (row-major reshape/recycle) exactly as in the
original."""

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate, varx_estimate_cml
from quanttoolbox.linalg.special_matrices import design, reshaper

data = np.array(
    [
        [1, 3.06, 1.34, 8.48, 28, 359.27, 102.96, 578.49],
        [1, 3.19, 1.44, 9.16, 35, 415.76, 114.38, 650.86],
        [1, 3.3, 1.54, 9.9, 37, 435.11, 118.23, 684.87],
        [1, 3.4, 1.71, 11.02, 36, 440.17, 120.45, 680.47],
        [1, 3.48, 1.89, 11.64, 29, 410.66, 116.25, 642.19],
        [1, 3.6, 1.99, 12.73, 47, 530.33, 140.27, 787.41],
        [1, 3.68, 2.22, 13.88, 50, 557.15, 143.84, 818.06],
        [1, 3.72, 2.43, 14.5, 35, 472.8, 128.2, 712.16],
        [1, 3.92, 2.43, 15.47, 33, 471.76, 126.65, 722.23],
        [1, 4.15, 2.31, 16.61, 40, 538.3, 141.05, 811.44],
        [1, 4.35, 2.39, 17.4, 38, 547.76, 143.71, 816.36],
        [1, 4.37, 2.63, 18.83, 37, 539, 142.37, 807.78],
        [1, 4.59, 2.69, 20.62, 56, 677.6, 173.13, 983.53],
        [1, 5.23, 3.35, 23.76, 88, 943.85, 223.21, 1292.99],
        [1, 6.04, 5.81, 26.52, 62, 893.42, 198.64, 1179.64],
        [1, 6.36, 6.38, 27.45, 51, 871, 191.89, 1134.78],
        [1, 7.04, 6.14, 30.28, 29, 793.93, 181.27, 1053.16],
        [1, 7.81, 6.14, 25.4, 22, 850.36, 180.56, 1085.91],
        [1, 8.09, 6.19, 28.84, 38, 967.42, 208.24, 1246.99],
        [1, 9.24, 6.69, 34.36, 41, 1102.61, 235.43, 1401.94],
    ]
)

x0 = data[:, 0:5]
y = data[:, 5:8]
x = np.column_stack([y, x0])

w = np.array([0, 1, 0, 2, 0, 3, 4, 0, 0, 5, 6, 7, 0, 8, 9, 0, 10, 0, 0, 11, 0, 0, 0, 12])
rr = design(w)
r = np.zeros(24)

# OLS estimate, see JHGLL page 660
sigma = None
labels = ("OLS", "GLS-2S", "GLS-3S", "LIML")

for i in range(4):
    if i == 3:
        result = varx_estimate_cml(y, x, p=0, restriction=(rr, r), tol=1e-20)
    else:
        result = varx_estimate(y, x, p=0, restriction=(rr, r), sigma=sigma, method="ls")
    sigma = result.sigma

    theta = result.theta
    g = reshaper(theta, 3, 3) - np.eye(3)
    b = reshaper(theta[9:24], 5, 3)
    rf = -b @ np.linalg.inv(g)

    print("=" * 60, labels[i])
    print("\nGamma:")
    print(np.round(g, 5))
    print("\nB:")
    print(np.round(b, 5))
    print("\nPI (reduced form):")
    print(np.round(rf, 5))
    print()

vech/xpnd (both orderings) and reshaper/reshapec compared — matrix/reshape1.py

"""Translated from Examples/matrix/reshape1.m -- vech/xpnd (both row- and
column-wise orderings) and reshaper/reshapec on a fixed 4x4 matrix.

The original also exercises `shiftr` (a lagged-shift primitive); that
function is not ported (superseded by NumPy roll/slicing idioms -- see
`matrix/shiftr1.m`'s own tracker entry) and is skipped here rather than
reimplemented."""

import numpy as np

from quanttoolbox.linalg.special_matrices import reshapec, reshaper, vech, xpnd

n = 4
col = np.arange(1, n + 1, dtype=float).reshape(-1, 1)  # seqa(1,1,n)
row = (1 + 0.25 * np.arange(n)).reshape(1, -1)  # seqa(1,0.25,n)'
y = col**row
print("y:")
print(y)

v = vech(y)
print("\nvech(y), row-wise (default):", v)
z = xpnd(v)
print("xpnd(v):")
print(z)

v = vech(y, method="c")
print("\nvech(y), column-wise:", v)
z = xpnd(v, method="c")
print("xpnd(v):")
print(z)

print("\nreshaper(y, 3, 2):")
print(reshaper(y, 3, 2))
print("reshaper(y, 5, 1):")
print(reshaper(y, 5, 1))
print("reshaper(y, 1, 5):")
print(reshaper(y, 1, 5))

print("\nreshapec(y, 3, 2):")
print(reshapec(y, 3, 2))
print("reshapec(y, 5, 1):")
print(reshapec(y, 5, 1))
print("reshapec(y, 1, 5):")
print(reshapec(y, 1, 5))