quanttoolbox.econometrics

econometrics.estimation (OLS/GMM/ML, Wald test)

Python alternatives

Hybrid: statsmodels.sandbox.regression.gmm.GMM, statsmodels.base.model.GenericLikelihoodModel, or the linearmodels package (more modern GMM/IV support) are more mature for standard use cases. Keep this module for the explicit theta = RR @ gamma + r linear-restriction interface, which none of the alternatives expose as directly.

quanttoolbox.econometrics.estimation

Constrained OLS, GMM, and Maximum Likelihood estimation with linear parameter restrictions.

Ported from QuantToolBox/ects/{ols_estimation,ols_constrained_estimation, gmm_estimation,gmm_constrained_estimation,ml_estimation, ml_constrained_estimation}.m

Translation notes:

• All three estimators share the same "linear restriction" pattern: the full parameter vector theta = RR @ gamma + r for a free parameter gamma, where RR/r default to the identity/zero (i.e. no restriction). This is ported as-is via an optional restriction=(RR, r) argument.
• The unconstrained convenience wrappers (ols_estimation.m, gmm_estimation.m, ml_estimation.m -- each just a nargin-dispatch stub calling the "constrained" version with default arguments) are not ported as separate functions; call the functions below with restriction=None for the same effect.
• GMM/ML optimization: MATLAB's fminunc/fmincon (with a trust-region-vs-quasi-Newton branch depending on whether analytical gradients/Hessians were supplied) is replaced throughout by scipy.optimize.minimize (BFGS by default, Newton-CG if an analytical Hessian is supplied) -- both are standard local optimizers for smooth objectives and the choice of exact algorithm doesn't change what's being estimated.
• GMM's iteratively-updated efficient weighting matrix uses a Bartlett-kernel (Newey-West-style) HAC covariance of the moments, ported directly (_bartlett_covariance).
• ML's three covariance estimator options (Hessian-based, OPG, and heteroskedasticity-consistent "sandwich") are all ported, selected via cov="hessian"|"opg"|"hc".
• MATLAB's global Print_Results/GMM_*/ML_* blocks are replaced by quanttoolbox.config.EstimationConfig.

gmm_estimation(moments_fn, sv, restriction=None, weights=1.0, jacobian_fn=None, weight_matrix=None, n_lags=0, config=None)

Iterated (two-step-and-beyond) efficient GMM estimation: theta = RR @ gamma + r, minimizing g(theta)' @ inv(W) @ g(theta) where g is the sample average of moments_fn(theta), re-estimating the efficient weighting matrix W (Bartlett/HAC, or a fixed weight_matrix if given) each iteration until convergence.

Original: ects/{gmm_estimation,gmm_constrained_estimation}.m

Source code in src/quanttoolbox/econometrics/estimation.py

def gmm_estimation(
    moments_fn: Callable[[np.ndarray], np.ndarray],
    sv: np.ndarray,
    restriction: tuple[np.ndarray, np.ndarray] | None = None,
    weights: np.ndarray | float = 1.0,
    jacobian_fn: Callable[[np.ndarray], np.ndarray] | None = None,
    weight_matrix: np.ndarray | None = None,
    n_lags: int = 0,
    config: EstimationConfig | None = None,
) -> GMMEstimationResult:
    """Iterated (two-step-and-beyond) efficient GMM estimation:
    theta = RR @ gamma + r, minimizing g(theta)' @ inv(W) @ g(theta) where
    g is the sample average of moments_fn(theta), re-estimating the
    efficient weighting matrix W (Bartlett/HAC, or a fixed weight_matrix if
    given) each iteration until convergence.

    Original: ects/{gmm_estimation,gmm_constrained_estimation}.m
    """
    if config is None:
        config = EstimationConfig()

    sv = np.asarray(sv, dtype=float).flatten()
    n_params = sv.shape[0]
    rr, r = _default_restriction(n_params) if restriction is None else restriction
    rr = np.asarray(rr, dtype=float)
    r = np.asarray(r, dtype=float).flatten()
    n_free = rr.shape[1]

    theta0 = rr @ sv + r
    h0 = np.atleast_2d(moments_fn(theta0))
    if h0.shape[0] == 1 and h0.shape[1] != 1:
        h0 = h0.reshape(-1, h0.shape[1]) if h0.ndim == 2 else h0
    m = h0.shape[1]
    n_obs = h0.shape[0]
    valid0 = ~np.isnan(h0).any(axis=1)
    n_valid = int(np.sum(valid0))
    df = n_valid - n_free

    if m < n_free:
        nan_p = np.full(n_free, np.nan)
        return GMMEstimationResult(
            theta=nan_p,
            stderr=nan_p,
            vcv=np.full((n_free, n_free), np.nan),
            q_min=np.nan,
            jacobian=np.full((m, n_free), np.nan),
            j_test=np.nan,
            j_test_pvalue=np.nan,
            df=df,
            n_obs=n_obs,
            n_obs_valid=n_valid,
            converged=False,
            n_iters=0,
        )

    inv_w = np.eye(m) if weight_matrix is None else np.asarray(weight_matrix, dtype=float)
    w = (
        np.broadcast_to(weights, (n_obs,))
        if np.isscalar(weights)
        else np.asarray(weights, dtype=float)
    )

    def _avg_moments(theta: np.ndarray) -> np.ndarray:
        h = np.atleast_2d(moments_fn(theta))
        h = w[:, None] * h
        valid = ~np.isnan(h).any(axis=1)
        return h[valid].mean(axis=0)

    def _objective(gamma: np.ndarray) -> float:
        theta = rr @ gamma + r
        g = _avg_moments(theta)
        return float(g @ inv_w @ g)

    gamma = sv.copy()
    converged = False
    n_iter = 0
    for n_iter in range(1, config.max_iters + 1):  # noqa: B007 (used after loop)
        result = minimize(_objective, gamma, method="BFGS")
        gamma_new = result.x
        if np.max(np.abs(gamma_new - gamma)) < config.tol:
            gamma = gamma_new
            converged = True
            break
        gamma = gamma_new

        if weight_matrix is None:
            theta = rr @ gamma + r
            h = np.atleast_2d(moments_fn(theta))
            h = w[:, None] * h
            valid = ~np.isnan(h).any(axis=1)
            inv_w = np.linalg.pinv(_bartlett_covariance(h[valid], n_lags))

    theta = rr @ gamma + r
    g_final = _avg_moments(theta)
    q_min = float(g_final @ inv_w @ g_final)

    if jacobian_fn is not None:
        d = jacobian_fn(theta)
    else:
        d = _numerical_jacobian(_avg_moments, theta)
    d = d @ rr

    try:
        vcv_gamma = np.linalg.inv(d.T @ inv_w @ d)
    except np.linalg.LinAlgError:
        vcv_gamma = np.full((n_free, n_free), np.nan)

    vcv = rr @ vcv_gamma @ rr.T / n_valid
    stderr = np.sqrt(np.diag(vcv))

    if m > n_free:
        j_test = n_valid * q_min
        j_test_pvalue = float(chi2.sf(j_test, m - n_free))
    else:
        j_test, j_test_pvalue = 0.0, 0.0

    return GMMEstimationResult(
        theta=theta,
        stderr=stderr,
        vcv=vcv,
        q_min=q_min,
        jacobian=d,
        j_test=j_test,
        j_test_pvalue=j_test_pvalue,
        df=df,
        n_obs=n_obs,
        n_obs_valid=n_valid,
        converged=converged,
        n_iters=n_iter,
    )

ml_estimation(logpdf_fn, sv, restriction=None, weights=1.0, jacobian_fn=None, hessian_fn=None, cov='hessian', config=None)

Maximum likelihood estimation: theta = RR @ gamma + r, maximizing the sum of logpdf_fn(theta) (a per-observation log-density) over the free parameter gamma.

cov="hessian" (default): asymptotic covariance from the (numerical or analytical) Hessian of the log-likelihood. cov="opg": outer-product- of-gradients estimator. cov="hc": heteroskedasticity-consistent "sandwich" estimator (Hessian and OPG combined).

Original: ects/{ml_estimation,ml_constrained_estimation}.m

Source code in src/quanttoolbox/econometrics/estimation.py

def ml_estimation(
    logpdf_fn: Callable[[np.ndarray], np.ndarray],
    sv: np.ndarray,
    restriction: tuple[np.ndarray, np.ndarray] | None = None,
    weights: np.ndarray | float = 1.0,
    jacobian_fn: Callable[[np.ndarray], np.ndarray] | None = None,
    hessian_fn: Callable[[np.ndarray], np.ndarray] | None = None,
    cov: str = "hessian",
    config: EstimationConfig | None = None,
) -> MLEstimationResult:
    """Maximum likelihood estimation: theta = RR @ gamma + r, maximizing the
    sum of logpdf_fn(theta) (a per-observation log-density) over the free
    parameter gamma.

    cov="hessian" (default): asymptotic covariance from the (numerical or
    analytical) Hessian of the log-likelihood. cov="opg": outer-product-
    of-gradients estimator. cov="hc": heteroskedasticity-consistent
    "sandwich" estimator (Hessian and OPG combined).

    Original: ects/{ml_estimation,ml_constrained_estimation}.m
    """
    if config is None:
        config = EstimationConfig()

    sv = np.asarray(sv, dtype=float).flatten()
    n_params = sv.shape[0]
    rr, r = _default_restriction(n_params) if restriction is None else restriction
    rr = np.asarray(rr, dtype=float)
    r = np.asarray(r, dtype=float).flatten()
    n_free = rr.shape[1]

    logl0 = np.atleast_1d(logpdf_fn(rr @ sv + r))
    n_obs = logl0.shape[0]
    w = (
        np.broadcast_to(weights, (n_obs,))
        if np.isscalar(weights)
        else np.asarray(weights, dtype=float)
    )

    def _neg_sum_logl(gamma: np.ndarray) -> float:
        theta = rr @ gamma + r
        logl = np.atleast_1d(logpdf_fn(theta))
        wl = w * logl
        return float(-np.sum(wl[~np.isnan(wl)]))

    jac_for_opt = None
    if jacobian_fn is not None:

        def jac_for_opt(gamma: np.ndarray) -> np.ndarray:
            theta = rr @ gamma + r
            j = np.atleast_2d(jacobian_fn(theta)) @ rr
            j = w[:, None] * j
            return -np.sum(j[~np.isnan(j).any(axis=1)], axis=0)

    result = minimize(_neg_sum_logl, sv, method="BFGS", jac=jac_for_opt)
    gamma = result.x
    theta = rr @ gamma + r

    logl = np.atleast_1d(logpdf_fn(theta))
    valid = ~np.isnan(logl)
    n_valid = int(np.sum(valid))
    df = n_valid - n_free
    sum_logl = float(np.sum(logl[valid]))

    def _neg_sum_logl_full(t: np.ndarray) -> float:
        ll = np.atleast_1d(logpdf_fn(t))
        wl = w * ll
        return float(-np.sum(wl[~np.isnan(wl)]))

    if cov in ("opg", "hc"):
        g = jacobian_fn(theta) if jacobian_fn is not None else _numerical_jacobian(logpdf_fn, theta)
        g = w[:, None] * g
        g_valid = g[valid] @ rr
        gg_valid = g_valid.T @ g_valid

    if cov == "opg":
        try:
            vcv_free = np.linalg.inv(gg_valid)
        except np.linalg.LinAlgError:
            vcv_free = np.full((n_free, n_free), np.nan)
        cov_type = "opg"
    else:
        h = (
            hessian_fn(theta)
            if hessian_fn is not None
            else -_numerical_hessian(_neg_sum_logl_full, theta)
        )
        h_valid = rr.T @ h @ rr
        try:
            inv_h_valid = np.linalg.inv(h_valid)
        except np.linalg.LinAlgError:
            inv_h_valid = np.full((n_free, n_free), np.nan)

        if cov == "hc":
            vcv_free = inv_h_valid @ gg_valid @ inv_h_valid
            cov_type = "hc"
        else:
            vcv_free = -inv_h_valid
            cov_type = "hessian"

    vcv = rr @ vcv_free @ rr.T
    stderr = np.sqrt(np.diag(vcv))

    residuals_full = np.full(n_obs, np.nan)
    residuals_full[valid] = logl[valid]

    return MLEstimationResult(
        theta=theta,
        stderr=stderr,
        vcv=vcv,
        log_l=residuals_full,
        sum_log_l=sum_logl,
        df=df,
        n_obs=n_obs,
        n_obs_valid=n_valid,
        cov_type=cov_type,
        converged=result.success,
    )

ols_estimation(y, x, restriction=None, weights=None)

Weighted, linearly-restricted OLS: beta = RR @ gamma + r, estimating the free parameter gamma by weighted least squares.

Original: ects/{ols_estimation,ols_constrained_estimation}.m

Source code in src/quanttoolbox/econometrics/estimation.py

def ols_estimation(
    y: np.ndarray,
    x: np.ndarray,
    restriction: tuple[np.ndarray, np.ndarray] | None = None,
    weights: np.ndarray | None = None,
) -> OLSEstimationResult:
    """Weighted, linearly-restricted OLS: beta = RR @ gamma + r, estimating
    the free parameter gamma by weighted least squares.

    Original: ects/{ols_estimation,ols_constrained_estimation}.m
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    n_obs = y.shape[0]
    n_params = x.shape[1]

    rr, r = _default_restriction(n_params) if restriction is None else restriction
    rr = np.asarray(rr, dtype=float)
    r = np.asarray(r, dtype=float).flatten()
    n_free = rr.shape[1]

    w = np.ones(n_obs) if weights is None else np.asarray(weights, dtype=float).flatten()

    valid = ~np.isnan(y) & ~np.isnan(x).any(axis=1)
    valid_idx = np.where(valid)[0]
    y_v, x_v, w_v = y[valid], x[valid], w[valid]
    n_valid = y_v.shape[0]

    wx = w_v[:, None] * x_v
    xwx = x_v.T @ wx
    xx = rr.T @ xwx @ rr
    rxw = (wx @ rr).T
    xy = rxw @ (y_v - x_v @ r)

    inv_xx = np.linalg.inv(xx)
    gamma = inv_xx @ xy
    beta = rr @ gamma + r

    u = y_v - x_v @ beta
    rss = np.sum(w_v * u**2)

    df_y = n_valid - 1
    df_x = n_free - 1
    df_u = df_y - df_x

    sigma2 = rss / df_u
    sigma = np.sqrt(sigma2)
    vcv = sigma2 * (rr @ inv_xx @ rr.T)
    stderr = np.sqrt(np.diag(vcv))

    t_stat = beta / np.where(stderr == 0, np.nan, stderr)
    p_value = 2 * t.sf(np.abs(t_stat), df_u)

    tss = np.sum(w_v * y_v**2)
    r_squared = 1 - rss / tss
    r_squared_adj = 1 - (rss / df_u) / (tss / df_y)

    y_bar = np.sum(w_v * y_v) / np.sum(w_v)
    yc = y_v - y_bar
    tss_c = np.sum(w_v * yc**2)
    r_squared_c = 1 - rss / tss_c
    r_squared_c_adj = 1 - (rss / df_u) / (tss_c / df_y)

    with np.errstate(divide="ignore", invalid="ignore"):
        f_stat = (r_squared_c / df_x) / ((1 - r_squared_c) / df_u)
    f_pvalue = f.sf(f_stat, df_x, df_u)

    residuals = np.full(n_obs, np.nan)
    residuals[valid_idx] = u

    return OLSEstimationResult(
        beta=beta,
        stderr=stderr,
        vcv=vcv,
        residuals=residuals,
        t_stat=t_stat,
        p_value=p_value,
        r_squared=r_squared,
        r_squared_adj=r_squared_adj,
        r_squared_centered=r_squared_c,
        r_squared_centered_adj=r_squared_c_adj,
        sigma=sigma,
        sigma2=sigma2,
        f_stat=f_stat,
        f_pvalue=f_pvalue,
        df_residual=df_u,
        n_obs=n_obs,
        n_obs_valid=n_valid,
    )

wald_test(constraint_fn, theta, vcv, n_obs)

Wald test of the (possibly nonlinear) hypothesis constraint_fn(theta) == 0.

Original: ects/wald_test.m

Source code in src/quanttoolbox/econometrics/estimation.py

def wald_test(
    constraint_fn: Callable[[np.ndarray], np.ndarray],
    theta: np.ndarray,
    vcv: np.ndarray,
    n_obs: int,
) -> WaldTestResult:
    """Wald test of the (possibly nonlinear) hypothesis constraint_fn(theta) == 0.

    Original: ects/wald_test.m
    """
    theta = np.asarray(theta, dtype=float).flatten()
    vcv = np.asarray(vcv, dtype=float)
    c0 = np.atleast_1d(constraint_fn(theta))
    g = c0.shape[0]

    jac = _numerical_jacobian(constraint_fn, theta)

    chi2_stat = float(c0 @ np.linalg.inv(jac @ vcv @ jac.T) @ c0)
    chi2_pvalue = float(chi2.sf(chi2_stat, g))

    f_stat = chi2_stat / g
    f_pvalue = float(f.sf(f_stat, g, n_obs))

    return WaldTestResult(
        chi2_stat=chi2_stat, chi2_pvalue=chi2_pvalue, f_stat=f_stat, f_pvalue=f_pvalue
    )

Examples

Beta-distribution MLE for LGD modeling — ects/ml3.py

"""Translated from Examples/ects/ml3.m -- Beta-distribution MLE (LGD --
Loss Given Default -- modeling), comparing the Hessian/OPG/HC covariance
estimators, on 13 observed LGD values."""

import numpy as np
from scipy.stats import beta as beta_dist

from quanttoolbox.econometrics.estimation import ml_estimation

lgd = np.array([0.68, 0.90, 0.22, 0.45, 0.17, 0.25, 0.89, 0.65, 0.75, 0.56, 0.87, 0.92, 0.46])

sv = np.array([1.0, 1.0])


def logpdf(theta):
    theta = np.sqrt(theta**2)  # force positivity, matches the original
    a, b = theta[0], theta[1]
    return np.log(beta_dist.pdf(lgd, a, b))


r_hess = ml_estimation(logpdf, sv, cov="hessian")
r_opg = ml_estimation(logpdf, sv, cov="opg")
r_hc = ml_estimation(logpdf, sv, cov="hc")

print("theta (a, b):", np.round(r_hess.theta, 4))
print("\nML covariance matrix (Hessian):")
print(np.round(r_hess.vcv, 5))
print("\nML covariance matrix (OPG):")
print(np.round(r_opg.vcv, 5))
print("\nML covariance matrix (heteroskedasticity-consistent):")
print(np.round(r_hc.vcv, 5))

Local Linear Trend model fit by time-domain ML, vs. Whittle — ects/kalman2c.py

"""Translated from Examples/ects/kalman2c.m -- Harvey [1990], pages
89-90: time-domain (exact) maximum-likelihood estimation of the Local
Linear Trend model's (sigma_epsilon, sigma_eta, sigma_zeta) on the same
61-observation Gnp.asc series as kalman2a.py/kalman2b.py, then runs the
Kalman filter with the estimated parameters (numeric core only, plot
dropped).

As in kalman1b.py/kalman1c.py, the original's `theta = sqrt(theta.^2)`
inside `LLT_ml` is translated as `np.abs(theta)`."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

y = np.array(
    [
        116.8,
        120.0,
        123.2,
        130.2,
        131.4,
        125.6,
        124.5,
        134.3,
        135.2,
        151.8,
        146.4,
        139.0,
        127.8,
        147.0,
        165.9,
        165.5,
        179.4,
        190.0,
        189.8,
        190.9,
        203.6,
        183.5,
        169.3,
        144.2,
        141.5,
        154.3,
        169.5,
        193.0,
        203.2,
        192.9,
        209.4,
        227.2,
        263.7,
        297.8,
        337.1,
        361.3,
        355.2,
        312.6,
        309.9,
        323.7,
        324.1,
        255.3,
        383.4,
        395.1,
        412.8,
        406.0,
        438.0,
        446.1,
        452.5,
        447.3,
        475.9,
        487.7,
        497.2,
        529.8,
        551.0,
        581.1,
        617.8,
        658.1,
        675.2,
        706.6,
        724.7,
    ],
    dtype=float,
)
nobs = y.shape[0]

a0 = np.array([y[0], 0.0])
p0 = np.zeros((2, 2))


def _ssm(sigma_epsilon: float, sigma_eta: float, sigma_zeta: float) -> StateSpaceModel:
    return StateSpaceModel(
        z=np.array([[1.0, 0.0]]),
        d=np.array([0.0]),
        h=np.array([[sigma_epsilon**2]]),
        t=np.array([[1.0, 1.0], [0.0, 1.0]]),
        c=np.array([0.0, 0.0]),
        r=np.eye(2),
        q=np.diag([sigma_eta**2, sigma_zeta**2]),
    )


def llt_ml(theta: np.ndarray) -> np.ndarray:
    sigma_epsilon, sigma_eta, sigma_zeta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta, sigma_zeta)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    return result.log_l


sv = 3.0 * np.ones(3)
ml_result = ml_estimation(llt_ml, sv)
theta = ml_result.theta

sigma_epsilon, sigma_eta, sigma_zeta = np.abs(theta)
print("theta (sigma_epsilon, sigma_eta, sigma_zeta):", np.round(theta, 4))
print("log-likelihood:", round(ml_result.sum_log_l, 4))

ssm = _ssm(sigma_epsilon, sigma_eta, sigma_zeta)
result = kalman_filter(ssm, y[:, None], a0, p0)

t = np.arange(1909, 1909 + nobs)
print("\nt, y, level a(t|t-1), slope a(t|t-1) -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, result.a_pred])[:10], 3))
print(np.round(np.column_stack([t, y, result.a_pred])[-10:], 3))

ML recovery of all free state-space model parameters — ects/kalman3d.py

"""Translated from Examples/ects/kalman3d.m -- maximum-likelihood
recovery of all 11 free parameters (Z diagonal, d[0], H diagonal, T,
c, Q) of the general state-space model used in kalman3a.py/kalman3b.py/
kalman3c.py, starting the optimizer from the model's own true parameter
values (as the original does) and comparing the recovered estimates
against them.

`theta([4 5 11]) = sqrt(theta([4 5 11]).^2)` in the original enforces
H[0,0], H[1,1], Q > 0 (indices 4, 5, 11 in MATLAB's 1-based numbering);
translated as `np.abs(theta[[3, 4, 10]])` here. `a0`/`P0` are
recomputed from the model's own steady state at every trial `theta`,
exactly as in `ssm_ml_fun`."""

import io

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter, steady_state

_KALMAN3_DATA = """
12.004191367990 5.082325511092
15.665734306703 6.508713750061
13.009699927942 3.153318564843
15.335177481404 3.972550347263
12.508230936087 5.832925230371
14.085284431457 3.966144631579
13.880704275605 4.341577981101
17.299643540570 5.261307380309
12.247804130775 3.316255057591
13.432744976803 2.048672540611
11.622880573894 7.676694410981
11.788856472009 3.495128971097
12.917862956802 6.138211990221
14.317625056476 5.482944588328
16.847846044264 4.682035523793
11.967832572522 4.767366961079
14.974022160932 4.887977730878
15.016692818794 5.719026158573
14.424733030162 4.591756017103
16.042748794722 5.727057061524
14.885645813926 5.669898953266
14.755022585609 6.118419596523
13.536221951032 3.602206553287
11.596395495070 3.633464391525
11.604635094167 4.411888939754
12.390307519857 6.690385683215
15.305430497460 2.781199059932
15.402296198009 2.968476593783
16.168052982092 7.873695626866
15.981037574492 4.520455068664
19.740491145225 4.774787134276
17.774423955212 8.041963181636
16.375436108110 8.144223054978
19.684112696330 7.149207939792
18.032141720118 5.323073136597
11.342074743746 3.524220139858
15.349548300273 3.689291210812
18.349772674601 4.527378947789
15.530812985344 2.950467678878
17.225938859816 5.287363621403
16.579428949350 5.842125841512
16.351459570867 5.529566663511
16.770518442192 5.821908296493
14.106124772987 3.111896017411
15.726291282545 4.762045600136
17.699211885524 6.004068704830
17.634176507780 8.031641953772
14.932325325953 4.317652273000
15.296105656851 3.679922559759
13.240165786848 4.406086727283
17.378133569049 6.483437034615
17.546097415114 7.379225581477
16.276632663535 4.752764527387
14.051049656425 4.949933424766
14.994843414610 4.064828477724
11.049136259484 5.860639265902
18.042263203996 5.894551158525
17.083048154770 3.911844584559
12.821851106294 4.515592694496
14.945384650234 2.067169324453
14.714474843610 5.152085003822
14.110763908126 3.007624240753
13.982036948425 4.638848725578
11.939397400605 2.060874318059
14.067684713478 5.489049813033
14.631150348061 4.833850145498
16.059681645779 3.547100987793
16.607934346919 5.691142902747
16.445288786249 4.048342876977
15.500879280182 5.105351576079
13.597040366017 4.629353517066
15.379042268757 6.350990982664
15.308553895014 3.938070863008
19.443989876593 4.343043602184
15.106704622680 4.941879404882
14.905124865490 5.317340054023
16.437928609781 4.966803648347
13.449475918200 4.371207295415
14.742788936007 6.125527718659
15.335008061543 5.647248674467
15.737287573398 3.226021249963
14.236120737521 4.285589175699
14.522833146604 3.212970566847
12.261176341209 5.803765210148
14.510640575162 5.453329322308
16.696274014570 4.779056934387
17.445925211911 4.377813743756
12.684520208790 5.293449220792
15.365482189411 5.775810076543
15.144318098993 4.350539527521
15.530949743147 5.853153253807
13.539603397937 3.598410097115
14.161037667699 5.142596818775
14.355274468497 4.362791955901
11.799576339582 4.597257887517
13.564221566673 4.706566007960
17.016076203763 5.190421267339
16.480283255317 3.844473935740
12.416545487942 3.250372697960
13.752720516863 4.947745693575
"""

y = np.loadtxt(io.StringIO(_KALMAN3_DATA))

ssm0 = StateSpaceModel(
    z=np.eye(2),
    d=np.array([10.0, 0.0]),
    h=np.array([[2.0, 0.0], [0.0, 1.0]]),
    t=np.array([[0.5, 0.3], [0.0, 0.2]]),
    c=np.array([1.0, 4.0]),
    r=np.array([[1.0], [1.0]]),
    q=np.array([[1.0]]),
)
a_bar0, p_bar0 = steady_state(ssm0)
result0 = kalman_filter(ssm0, y, a_bar0, p_bar0)
print("Value of the log-likelihood function:", round(float(np.sum(result0.log_l)), 4))


def ssm_ml_fun(theta: np.ndarray) -> np.ndarray:
    theta = theta.copy()
    theta[[3, 4, 10]] = np.abs(theta[[3, 4, 10]])

    ssm = StateSpaceModel(
        z=np.array([[theta[0], 0.0], [0.0, theta[1]]]),
        d=np.array([theta[2], 0.0]),
        h=np.array([[theta[3], 0.0], [0.0, theta[4]]]),
        t=np.array([[theta[5], theta[6]], [0.0, theta[7]]]),
        c=np.array([theta[8], theta[9]]),
        r=np.array([[1.0], [1.0]]),
        q=np.array([[theta[10]]]),
    )
    a_bar, p_bar = steady_state(ssm)
    result = kalman_filter(ssm, y, a_bar, p_bar)
    return result.log_l


sv = np.array([1.0, 1.0, 10.0, 2.0, 1.0, 0.5, 0.3, 0.2, 1.0, 4.0, 1.0])

ml_result = ml_estimation(ssm_ml_fun, sv)
theta_hat = ml_result.theta

print("\n True           Estimated")
print(" values         values")
print(np.round(np.column_stack([sv, theta_hat]), 5))

Numerical vs. analytical Jacobian/Hessian across covariance estimators — ects/ml4.py

"""Translated from Examples/ects/ml4.m -- Gaussian-MLE with numerical vs.
analytical Jacobian/Hessian, comparing all three covariance estimators
under both, on a 10-observation dataset. Note theta[4] is parametrized as
sigma^2 directly here (not sigma), matching the original's own warning
comment."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation

data = np.array(
    [
        [1.1, 1.0, 2.4, 3.6, 0.3],
        [20.4, 1.0, 1.1, 3.8, 5.9],
        [17.1, 1.0, 5.1, 6.3, 6.1],
        [30.9, 1.0, 2.7, 2.4, 9.5],
        [22.2, 1.0, 3.3, 3.0, 7.4],
        [9.1, 1.0, 1.0, 5.4, 4.9],
        [39.2, 1.0, 9.6, 2.8, 8.1],
        [3.1, 1.0, 2.9, 4.4, 1.0],
        [7.2, 1.0, 4.2, 5.6, 1.7],
        [27.6, 1.0, 8.1, 1.7, 5.4],
    ]
)
y = data[:, 0]
x = data[:, 1:5]


def logpdf(theta):
    beta = theta[:4]
    sigma2 = theta[4]
    u = y - x @ beta
    return -0.5 * np.log(2 * np.pi) - 0.5 * np.log(sigma2) - 0.5 * (u * u) / sigma2


def jacobian(theta):
    beta = theta[:4]
    sigma2 = theta[4]
    sigma4 = sigma2**2
    u = y - x @ beta
    g_beta = (x * u[:, None]) / sigma2
    g_sigma2 = -1 / (2 * sigma2) + (u * u) / (2 * sigma4)
    return np.column_stack([g_beta, g_sigma2])


def hessian(theta):
    beta = theta[:4]
    sigma2 = theta[4]
    sigma4 = sigma2**2
    sigma6 = sigma2**3
    u = y - x @ beta
    n = y.shape[0]

    h11 = -(x.T @ x) / sigma2
    h12 = -(x.T @ u) / sigma4
    h22 = n / (2 * sigma4) - (u @ u) / sigma6
    h = np.zeros((5, 5))
    h[:4, :4] = h11
    h[:4, 4] = h12
    h[4, :4] = h12
    h[4, 4] = h22
    return h


sv = np.ones(5)

results = {}
for cov in ("hessian", "opg", "hc"):
    results[f"num_{cov}"] = ml_estimation(logpdf, sv, cov=cov)
    results[f"ana_{cov}"] = ml_estimation(
        logpdf, sv, cov=cov, jacobian_fn=jacobian, hessian_fn=hessian
    )

for label, key in [
    ("Numerical Hessian", "num_hessian"),
    ("Numerical OPG", "num_opg"),
    ("Numerical heteroskedasticity-consistent", "num_hc"),
    ("Analytical Hessian", "ana_hessian"),
    ("Analytical OPG", "ana_opg"),
    ("Analytical heteroskedasticity-consistent", "ana_hc"),
]:
    print(f"\nML covariance matrix ({label}):")
    print(np.round(results[key].vcv, 5))

OLS vs. Gaussian-MLE covariance: Hessian, OPG, and HC estimators — ects/ml2.py

"""Translated from Examples/ects/ml2.m -- OLS covariance matrix vs.
Gaussian-MLE covariance matrix under all three `ml_estimation` cov
estimators (Hessian, OPG, heteroskedasticity-consistent sandwich), on a
10-observation dataset."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation, ols_estimation

data = np.array(
    [
        [1.5, 1.0, 2.4, 3.6, 0.3],
        [20.4, 1.0, 1.1, 3.8, 5.9],
        [17.1, 1.0, 5.1, 6.3, 6.1],
        [30.9, 1.0, 2.7, 2.4, 9.5],
        [22.2, 1.0, 3.3, 3.0, 7.4],
        [9.1, 1.0, 1.0, 5.4, 4.9],
        [39.2, 1.0, 9.6, 2.8, 8.1],
        [3.1, 1.0, 2.9, 4.4, 1.0],
        [7.2, 1.0, 4.2, 5.6, 1.7],
        [27.6, 1.0, 8.1, 1.7, 5.4],
    ]
)
y = data[:, 0]
x = data[:, 1:5]

r1 = ols_estimation(y, x)
sigma1 = r1.sigma
theta1 = np.concatenate([r1.beta, [sigma1]])


def logpdf(theta):
    beta = theta[:4]
    sigma2 = theta[4] ** 2
    u = y - x @ beta
    return -0.5 * np.log(2 * np.pi) - 0.5 * np.log(sigma2) - 0.5 * (u * u) / sigma2


r_hess = ml_estimation(logpdf, theta1, cov="hessian")
r_opg = ml_estimation(logpdf, theta1, cov="opg")
r_hc = ml_estimation(logpdf, theta1, cov="hc")

print("OLS covariance matrix:")
print(np.round(r1.vcv, 5))
print("\nML covariance matrix (Hessian):")
print(np.round(r_hess.vcv, 5))
print("\nML covariance matrix (OPG):")
print(np.round(r_opg.vcv, 5))
print("\nML covariance matrix (heteroskedasticity-consistent):")
print(np.round(r_hc.vcv, 5))

OLS vs. Gaussian-MLE covariance: Hessian, OPG, and HC standard errors — stats/cov1.py

"""Translated from Examples/stats/cov1.m -- OLS in closed form vs. maximum
likelihood (Gaussian log-likelihood, same model), comparing standard
errors from the Hessian, OPG, and "sandwich" (HC) covariance estimators.

The original's `ml_robust_vcv` helper computing 5 partially-redundant
covariance variants (2 Hessian-based, 1 pure-OPG, 1 sandwich, plus a
duplicate) isn't ported as a standalone function; this package's
`ml_estimation(..., cov=...)` covers the same three conceptual estimators
directly (`cov="hessian"`, `"opg"`, `"hc"`), used here instead.

The original draws x/u from MATLAB's unseeded `rand`/`randn`; a fixed seed
(`np.random.default_rng(0)`) is substituted here. `ml_ols.m` (the
per-observation log-density) is inlined below rather than kept as a
separate file, matching its role as a helper, not a standalone example."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation

rng = np.random.default_rng(0)
n = 100
beta_true = np.array([1.0, 2.0])
sigma_true = 0.20
x = 10 * (rng.random((n, 2)) - 0.5)
u = sigma_true * rng.standard_normal(n)
y = x @ beta_true + u

beta_hat = np.linalg.inv(x.T @ x) @ (x.T @ y)
u_hat = y - x @ beta_hat
sigma_hat = np.std(u_hat, ddof=1)
cov_beta = (sigma_hat**2) * np.linalg.inv(x.T @ x)
stderr_ols = np.sqrt(np.diag(cov_beta))

print("OLS: beta_true, beta_hat, stderr")
print(np.round(np.column_stack([beta_true, beta_hat, stderr_ols]), 4))


def logpdf(theta, y, x):
    k = x.shape[1]
    beta = theta[:k]
    sigma2 = theta[k] ** 2
    u = y - x @ beta
    return -0.5 * np.log(2 * np.pi) - 0.5 * np.log(sigma2) - 0.5 * (u**2) / sigma2


sv = np.concatenate([beta_true, [sigma_true]])

r_hess = ml_estimation(lambda th: logpdf(th, y, x), sv, cov="hessian")
r_opg = ml_estimation(lambda th: logpdf(th, y, x), sv, cov="opg")
r_hc = ml_estimation(lambda th: logpdf(th, y, x), sv, cov="hc")

print("\nMLE (Gaussian log-likelihood): sv, theta, stderr(hessian), stderr(opg), stderr(hc)")
print(
    np.round(
        np.column_stack([sv, r_hess.theta, r_hess.stderr, r_opg.stderr, r_hc.stderr]),
        4,
    )
)
print("\nconverged:", r_hess.converged, r_opg.converged, r_hc.converged)

OLS, LAD, quantile, and SVM regression compared — svm/svm6.py

"""Translated from Examples/svm/svm6.m -- comparison of OLS, LAD, and
SVM-regression (LS-SVM and epsilon-SVM, at two different C values)
estimates on the same 10-obs dataset used in ols1.m/robust1.m."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression
from quanttoolbox.svm.svm import svm_regression_primal

data = np.array(
    [
        [1.5, 1.0, 2.4, 3.6, 0.3],
        [20.4, 1.0, 1.1, 3.8, 5.9],
        [17.1, 1.0, 5.1, 6.3, 6.1],
        [30.9, 1.0, 2.7, 2.4, 9.5],
        [22.2, 1.0, 3.3, 3.0, 7.4],
        [9.1, 1.0, 1.0, 5.4, 4.9],
        [39.2, 1.0, 9.6, 2.8, 8.1],
        [3.1, 1.0, 2.9, 4.4, 1.0],
        [7.2, 1.0, 4.2, 5.6, 1.7],
        [27.6, 1.0, 8.1, 1.7, 5.4],
    ]
)
y = data[:, 0]
x_full = data[:, 1:5]

beta_ols = ols_estimation(y, x_full).beta
beta_lad = lad_regression(y, x_full).beta
beta_lad2, _, _ = quantile_regression(y, x_full, tau=0.5)

x = x_full[:, 1:4]  # drop the intercept and 4th column for SVM (3 predictors)

r_ls1 = svm_regression_primal(y, x, c=1)
r_eps1 = svm_regression_primal(y, x, c=1, epsilon=1)
r_ls2 = svm_regression_primal(y, x, c=100000)
r_eps2 = svm_regression_primal(y, x, c=100000, epsilon=0)

print("OLS:", np.round(beta_ols, 3))
print("LAD:", np.round(beta_lad, 3))
print("Quantile(0.5):", np.round(beta_lad2, 3))
print("SVM LS (C=1):", round(r_ls1.beta0, 3), np.round(r_ls1.beta, 3))
print("SVM eps (C=1, eps=1):", round(r_eps1.beta0, 3), np.round(r_eps1.beta, 3))
print("SVM LS (C=1e5):", round(r_ls2.beta0, 3), np.round(r_ls2.beta, 3))
print("SVM eps (C=1e5, eps=0):", round(r_eps2.beta0, 3), np.round(r_eps2.beta, 3))

OLS, median, LAD, and Huber regression compared — ects/robust2.py

"""Translated from Examples/ects/robust2.m -- OLS vs. median (quantile,
alpha=0.5) regression, LAD regression, and Huber regression, on a
simulated 3-predictor dataset.

The original draws x/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import huber_regression, lad_regression

rng = np.random.default_rng(0)
n, k = 100, 3
beta_true = np.arange(1, k + 1, dtype=float)
sigma = 0.20

x = rng.random((n, k))
y = x @ beta_true + sigma * rng.standard_normal(n)

r_ols = ols_estimation(y, x)
u_ols = r_ols.residuals
print("OLS beta:", np.round(r_ols.beta, 4))

alpha = 0.50
beta_med, _, _ = quantile_regression(y, x, alpha)
print("\nMedian regression beta:", np.round(beta_med, 4))

r_lad = lad_regression(y, x)
print("\nLAD regression beta:", np.round(r_lad.beta, 4), "converged:", r_lad.converged)

c = np.quantile(np.abs(u_ols), 0.90)
r_huber = huber_regression(y, x, c=c)
print(
    f"\nHuber regression (c={c:.4f}) beta:",
    np.round(r_huber.beta, 4),
    "converged:",
    r_huber.converged,
)

OLS, MLE, and GMM under an optional linear restriction — ects/gmm1.py

"""Translated from Examples/ects/gmm1.m -- linear regression estimated
three ways (OLS, Gaussian MLE, GMM with moment conditions matching OLS's
normal equations plus a second-moment condition), each in an unconstrained
and a beta[3]=1 restricted form, on a 10-observation dataset with one
missing y value.

The original's `theta7` (a *second*, independently-coded way of enforcing
the same beta[3]=1 restriction, by hard-wiring it inside the moment
function instead of using the RR/r restriction machinery) is not
re-translated separately -- `theta6` below already demonstrates the
restriction via `gmm_estimation`'s own `restriction=` parameter, the same
mechanism `ols_estimation`/`ml_estimation` use for `theta2`/`theta4`."""

import numpy as np

from quanttoolbox.econometrics.estimation import gmm_estimation, ml_estimation, ols_estimation

data = np.array(
    [
        [np.nan, 1.0, 2.4, 3.6, 0.3],
        [20.4, 1.0, 1.1, 3.8, 5.9],
        [17.1, 1.0, 5.1, 6.3, 6.1],
        [30.9, 1.0, 2.7, 2.4, 9.5],
        [22.2, 1.0, 3.3, 3.0, 7.4],
        [9.1, 1.0, 1.0, 5.4, 4.9],
        [39.2, 1.0, 9.6, 2.8, 8.1],
        [3.1, 1.0, 2.9, 4.4, 1.0],
        [7.2, 1.0, 4.2, 5.6, 1.7],
        [27.6, 1.0, 8.1, 1.7, 5.4],
    ]
)
y = data[:, 0]
x = data[:, 1:5]

# Unconstrained OLS, and OLS with beta[3] restricted to 1
r1 = ols_estimation(y, x)
sigma1 = r1.sigma
rr_ols = np.eye(4, 3)
r_ols = np.array([0.0, 0.0, 0.0, 1.0])
r2 = ols_estimation(y, x, restriction=(rr_ols, r_ols))
sigma2 = r2.sigma

theta1 = np.concatenate([r1.beta, [sigma1]])
theta2 = np.concatenate([r2.beta, [sigma2]])
sv = theta1 + 1.0


def logpdf(theta):
    beta = theta[:4]
    sigma2_ = theta[4] ** 2
    u = y - x @ beta
    return -0.5 * np.log(2 * np.pi) - 0.5 * np.log(sigma2_) - 0.5 * (u * u) / sigma2_


def moments(theta):
    beta = theta[:4]
    sigma2_ = theta[4] ** 2
    u = y - x @ beta
    h = np.zeros((u.shape[0], 5))
    h[:, 0] = u
    h[:, 1] = u * u - sigma2_
    h[:, 2:5] = u[:, None] * x[:, 1:4]
    return h


rr_ml = np.eye(5, 4)
rr_ml[3, 3] = 0
rr_ml[4, 3] = 1
r_ml = np.array([0.0, 0.0, 0.0, 1.0, 0.0])

r3 = ml_estimation(logpdf, sv)
r4 = ml_estimation(logpdf, sv[:4], restriction=(rr_ml, r_ml))
r5 = gmm_estimation(moments, sv)
r6 = gmm_estimation(moments, sv[:4], restriction=(rr_ml, r_ml))

print("Unconstrained: OLS, MLE, GMM")
print(np.round(np.column_stack([theta1, r3.theta, r5.theta]), 5))
print("\nbeta[3]=1 restricted: OLS, MLE, GMM")
print(np.round(np.column_stack([theta2, r4.theta, r6.theta]), 5))

OLS/median/LAD plus quantile regression via IRLS and exact LP — ects/robust3.py

"""Translated from Examples/ects/robust3.m -- OLS vs. median regression,
LAD regression, and quantile regression at alpha=0.90, both via the IRLS
M-estimator (`quantile_m_regression`) and the exact LP formulation
(`quantile_regression`), on a simulated 3-predictor dataset.

The original draws x/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression, quantile_m_regression

rng = np.random.default_rng(0)
n, k = 100, 3
beta_true = np.arange(1, k + 1, dtype=float)
sigma = 0.20

x = rng.random((n, k))
y = x @ beta_true + sigma * rng.standard_normal(n)

r_ols = ols_estimation(y, x)
u_ols = r_ols.residuals
print("OLS beta:", np.round(r_ols.beta, 4))

alpha = 0.50
beta_med, _, _ = quantile_regression(y, x, alpha)
print("\nMedian regression beta:", np.round(beta_med, 4))

r_lad = lad_regression(y, x)
print("\nLAD regression beta:", np.round(r_lad.beta, 4), "converged:", r_lad.converged)

alpha = 0.90
r_qm = quantile_m_regression(y, x, alpha)
print(f"\nQuantile (alpha={alpha}) regression via IRLS M-estimator, beta:", np.round(r_qm.beta, 4))

beta_q, _, _ = quantile_regression(y, x, alpha)
print(f"Quantile (alpha={alpha}) regression via exact LP, beta:        ", np.round(beta_q, 4))

OLS/SVM-LS and quantile/SVM-epsilon regression on synthetic data — svm/svm8.py

"""Translated from Examples/svm/svm8.m -- OLS, LAD, quantile, and SVM
regression (primal, at two very different C values, LS- and
epsilon-insensitive) compared on a larger (n=1000) simulated dataset.

The original draws x/beta/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here for reproducibility
-- this is the fixed-seed substitution the tracker notes svm8.m as still
needing."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression
from quanttoolbox.svm.svm import svm_regression_primal

rng = np.random.default_rng(0)
n, k = 1000, 4
x = rng.random((n, k))
beta_true = 5 * rng.random(k)
u = 0.20 * rng.standard_normal(n)
beta0_true = -3.0
y = beta0_true + x @ beta_true + u

x_design = np.column_stack([np.ones(n), x])

beta_ols = ols_estimation(y, x_design).beta
beta_lad = lad_regression(y, x_design).beta
beta_lad2, _, _ = quantile_regression(y, x_design, tau=0.5)

# SVM regression (primal), C=1
r_ls1 = svm_regression_primal(y, x, c=1)  # LS-SVM
beta_svm_ls = np.concatenate([[r_ls1.beta0], r_ls1.beta])
r_eps1 = svm_regression_primal(y, x, c=1, epsilon=1)  # epsilon-SVM
beta_svm_epsilon = np.concatenate([[r_eps1.beta0], r_eps1.beta])

# SVM regression (primal), C=1000 (effectively unregularized)
r_ls2 = svm_regression_primal(y, x, c=1000)  # LS-SVM
beta_svm_ls2 = np.concatenate([[r_ls2.beta0], r_ls2.beta])
r_eps2 = svm_regression_primal(y, x, c=1000, epsilon=0)  # epsilon-SVM, epsilon=0
beta_svm_epsilon2 = np.concatenate([[r_eps2.beta0], r_eps2.beta])

results = np.column_stack(
    [beta_ols, beta_lad, beta_lad2, beta_svm_ls, beta_svm_epsilon, beta_svm_ls2, beta_svm_epsilon2]
)
print("Comparison of OLS, LAD, quantile(0.5), and SVM estimates")
print(
    "columns: OLS, LAD, Quantile(0.5), SVM-LS(C=1), SVM-eps(C=1,eps=1), SVM-LS(C=1000), SVM-eps(C=1000,eps=0)"
)
print(np.round(results, 3))

print("\nOLS vs. SVM-LS at large C (expect near-identical -- unregularized limit):")
print(np.round(np.column_stack([beta_ols, beta_svm_ls2]), 3))

print(
    "\nQuantile(0.5) vs. SVM-eps at large C, epsilon=0 (expect near-identical -- both approach LAD):"
)
print(np.round(np.column_stack([beta_lad2, beta_svm_epsilon2]), 3))

Time-domain ML of Kalman noise variances, two starting states — ects/kalman1b.py

"""Translated from Examples/ects/kalman1b.m -- Harvey [1990],
"Forecasting, Structural Time Series and the Kalman Filter", pages
89-90: time-domain (exact) maximum-likelihood estimation of the
local-level model's (sigma_epsilon, sigma_eta) on the same 71-
observation Purse.asc series as panel1.py/whittle1.py, comparing two
choices of the Kalman filter's initial state a0 used *during
estimation* (0, vs. y[0]) -- then re-running the filter with each
estimated theta (using a0 = y[0], P0 = 0 for both, matching the
original's separate post-estimation filter block) to compare the
resulting one-step-ahead predictions.

The original's `theta = sqrt(theta.^2)` inside `LL_ml` (a roundabout way
of enforcing sigma_epsilon, sigma_eta > 0 while keeping the objective
smooth for the optimizer) is translated literally as `np.abs(theta)`.
`ml_estimation`'s `logpdf_fn` here returns the Kalman filter's
per-observation `log_l` array (the natural per-observation log-density
for this model), rather than a single pre-summed scalar (which
`ml_estimation` would also accept, but only as a degenerate 1-observation
case)."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

y = np.array(
    [
        10,
        15,
        10,
        10,
        12,
        10,
        7,
        17,
        10,
        14,
        8,
        17,
        14,
        18,
        3,
        9,
        11,
        10,
        6,
        12,
        14,
        10,
        25,
        29,
        33,
        33,
        12,
        19,
        16,
        19,
        19,
        12,
        34,
        15,
        36,
        29,
        26,
        21,
        17,
        19,
        13,
        20,
        24,
        12,
        6,
        14,
        6,
        12,
        9,
        11,
        17,
        12,
        8,
        14,
        14,
        12,
        5,
        8,
        10,
        3,
        16,
        8,
        8,
        7,
        12,
        6,
        10,
        8,
        10,
        5,
        7,
    ],
    dtype=float,
)
nobs = y.shape[0]


def _ssm(sigma_epsilon: float, sigma_eta: float) -> StateSpaceModel:
    return StateSpaceModel(
        z=np.array([[1.0]]),
        d=np.array([0.0]),
        h=np.array([[sigma_epsilon**2]]),
        t=np.array([[1.0]]),
        c=np.array([0.0]),
        r=np.array([[1.0]]),
        q=np.array([[sigma_eta**2]]),
    )


def ll_ml(theta: np.ndarray, a0: float) -> np.ndarray:
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], np.array([a0]), np.array([[0.0]]))
    return result.log_l


sv = np.array([3.0, 1.0])

theta1 = ml_estimation(lambda theta: ll_ml(theta, 0.0), sv).theta
theta2 = ml_estimation(lambda theta: ll_ml(theta, y[0]), sv).theta

print("theta1 (sigma_epsilon, sigma_eta), a0=0 during estimation:", np.round(theta1, 4))
print("theta2 (sigma_epsilon, sigma_eta), a0=y[0] during estimation:", np.round(theta2, 4))

a0 = np.array([y[0]])
p0 = np.array([[0.0]])

y_cond = np.zeros((nobs, 2))
for i, theta in enumerate((theta1, theta2)):
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    y_cond[:, i] = result.y_pred[:, 0]

t = np.arange(nobs)
print("\nt, y, y(t|t-1) [theta1], y(t|t-1) [theta2] -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, y_cond])[:10], 3))
print(np.round(np.column_stack([t, y, y_cond])[-10:], 3))

Time-domain vs. frequency-domain (Whittle) Kalman ML — ects/kalman1c.py

"""Translated from Examples/ects/kalman1c.m -- Harvey [1990],
"Forecasting, Structural Time Series and the Kalman Filter", pages
89-90: compares time-domain (exact) maximum likelihood against
frequency-domain (Whittle) maximum likelihood estimation of the
local-level model's (sigma_epsilon, sigma_eta), on the same Purse.asc
series as kalman1b.py/panel1.py/whittle1.py.

As in the original, the time-domain likelihood is evaluated with a0 = 0,
P0 = 0 during estimation, but the final filter run used to compare the
two estimates' one-step-ahead predictions uses a0 = y[0], P0 = 0 (the
original's own, slightly inconsistent, choice -- preserved here rather
than "fixed")."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter
from quanttoolbox.econometrics.whittle import whittle_local_level

y = np.array(
    [
        10,
        15,
        10,
        10,
        12,
        10,
        7,
        17,
        10,
        14,
        8,
        17,
        14,
        18,
        3,
        9,
        11,
        10,
        6,
        12,
        14,
        10,
        25,
        29,
        33,
        33,
        12,
        19,
        16,
        19,
        19,
        12,
        34,
        15,
        36,
        29,
        26,
        21,
        17,
        19,
        13,
        20,
        24,
        12,
        6,
        14,
        6,
        12,
        9,
        11,
        17,
        12,
        8,
        14,
        14,
        12,
        5,
        8,
        10,
        3,
        16,
        8,
        8,
        7,
        12,
        6,
        10,
        8,
        10,
        5,
        7,
    ],
    dtype=float,
)
nobs = y.shape[0]


def _ssm(sigma_epsilon: float, sigma_eta: float) -> StateSpaceModel:
    return StateSpaceModel(
        z=np.array([[1.0]]),
        d=np.array([0.0]),
        h=np.array([[sigma_epsilon**2]]),
        t=np.array([[1.0]]),
        c=np.array([0.0]),
        r=np.array([[1.0]]),
        q=np.array([[sigma_eta**2]]),
    )


def ll_ml(theta: np.ndarray) -> np.ndarray:
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], np.array([0.0]), np.array([[0.0]]))
    return result.log_l


sv = np.array([3.0, 1.0])

theta1 = ml_estimation(ll_ml, sv).theta
theta2 = whittle_local_level(y, sv).theta

print("theta1 (sigma_epsilon, sigma_eta), time-domain ML:", np.round(theta1, 4))
print("theta2 (sigma_epsilon, sigma_eta), frequency-domain (Whittle) ML:", np.round(theta2, 4))

a0 = np.array([y[0]])
p0 = np.array([[0.0]])

y_cond = np.zeros((nobs, 2))
for i, theta in enumerate((theta1, theta2)):
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    y_cond[:, i] = result.y_pred[:, 0]

t = np.arange(nobs)
print("\nt, y, y(t|t-1) [time-domain ML], y(t|t-1) [Whittle ML] -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, y_cond])[:10], 3))
print(np.round(np.column_stack([t, y, y_cond])[-10:], 3))

Time-varying-coefficient model with variances estimated by ML — ects/kalman4b.py

"""Translated from Examples/ects/kalman4b.m -- same simulated
time-varying-coefficient regression as kalman4a.py, but instead of using
the true (sigma_epsilon, sigma_beta1, sigma_beta2), estimates them by
maximum likelihood (starting from sv = [1, 1, 1]) and re-runs the
time-varying Kalman filter with the estimated variances, comparing the
recovered beta_t path against the simulated true path.

As in kalman1b.py/kalman2c.py/kalman3d.py, `theta = sqrt(theta.^2)` is
translated as `np.abs(theta)`."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

rng = np.random.default_rng(0)
n_t = 200

sigma1 = 0.5
sigma2 = 0.25
sigma = 1.0

beta1 = np.cumsum(sigma1 * rng.standard_normal(n_t))
beta2 = np.cumsum(sigma2 * rng.standard_normal(n_t))
beta = np.column_stack([beta1, beta2])
x = rng.random((n_t, 2))

y = np.sum(x * beta, axis=1) + sigma * rng.standard_normal(n_t)

a0 = np.zeros(2)
p0 = np.zeros((2, 2))

z = x[None, :, :].transpose(0, 2, 1)  # (1, 2, nT)
d = np.zeros((1, n_t))
t_mat = np.repeat(np.eye(2)[:, :, None], n_t, axis=2)
c = np.zeros((2, n_t))
r = np.repeat(np.eye(2)[:, :, None], n_t, axis=2)


def _ssm(sigma_epsilon: float, sigma_beta1: float, sigma_beta2: float) -> StateSpaceModel:
    h = np.full((1, 1, n_t), sigma_epsilon**2)
    q = np.repeat(np.diag([sigma_beta1**2, sigma_beta2**2])[:, :, None], n_t, axis=2)
    return StateSpaceModel(z=z, d=d, h=h, t=t_mat, c=c, r=r, q=q)


def ssm_logl(theta: np.ndarray) -> np.ndarray:
    sigma_epsilon, sigma_beta1, sigma_beta2 = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_beta1, sigma_beta2)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    return result.log_l


sv = np.ones(3)
ml_result = ml_estimation(ssm_logl, sv)
theta_hat = np.abs(ml_result.theta)

sigma_epsilon_hat, sigma1_hat, sigma2_hat = theta_hat
print(
    "Estimated (sigma_epsilon, sigma_beta1, sigma_beta2):",
    np.round(theta_hat, 4),
    "-- true:",
    [sigma, sigma1, sigma2],
)

ssm_hat = _ssm(sigma_epsilon_hat, sigma1_hat, sigma2_hat)
result = kalman_filter(ssm_hat, y[:, None], a0, p0)
at = result.a_filt

t = np.arange(1, n_t + 1)
print("\nt, true beta1, filtered beta1, true beta2, filtered beta2 -- first/last 10:")
print(np.round(np.column_stack([t, beta[:, 0], at[:, 0], beta[:, 1], at[:, 1]])[:10], 4))
print(np.round(np.column_stack([t, beta[:, 0], at[:, 0], beta[:, 1], at[:, 1]])[-10:], 4))

VAR(0) equivalence to plain OLS — ects/varx3.py

"""Translated from Examples/ects/varx3.m -- a trivial 5-observation
example comparing plain OLS against `varx_estimate(..., p=0)` on the
same regressor matrix. With no autoregressive lags (p=0), VARX
estimation collapses to OLS, so `theta`/`beta` from the second call
should match `beta` from the first."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.econometrics.var import varx_estimate

y = np.array([2, 3, 1, 7, 5], dtype=float)

x = np.array(
    [
        [1, 3, 2],
        [2, 3, 1],
        [7, 1, 7],
        [5, 3, 1],
        [3, 5, 5],
    ],
    dtype=float,
)

ols_result = ols_estimation(y, np.column_stack([np.ones(5), x]))
print("OLS beta:", np.round(ols_result.beta, 4))

varx_result = varx_estimate(y[:, None], np.column_stack([np.ones(5), x]), p=0, method="ls")
print("\nVARX(p=0) theta:", np.round(varx_result.theta, 4))
print("VARX(p=0) beta:", np.round(varx_result.beta, 4))

Wald test for no Granger-causality between VAR variables — ects/varx1b.py

"""Translated from Examples/ects/varx1b.m -- Lutkepohl [1991], section
3.6: Wald test for no Granger-causality from income/consumption to
investment in the VAR(2) model estimated in varx1a.py (same
log-differenced Lutkepohl.asc data), i.e. H0: theta[4] = theta[7] =
theta[13] = theta[16] = 0 (1-indexed positions in the 21-element stacked
coefficient vector).

`RR = [design([4;7;13;16]) zeros(4,5)]` is only used to build the
constraint function `C(theta) = RR @ theta` -- it is not a restriction
passed to `varx_estimate` here, so it's translated as a direct 4-element
index/select into theta rather than via `design`."""

import io

import numpy as np

from quanttoolbox.econometrics.estimation import wald_test
from quanttoolbox.econometrics.var import varx_estimate

_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))

result = varx_estimate(y, x, p=2, method="ls")
theta = result.theta
vcv = result.vcv

# Test for no Granger-causality from income/consumption to investment:
# theta[4] = theta[7] = theta[13] = theta[16] = 0 (1-indexed)
idx = np.array([4, 7, 13, 16]) - 1


def wald_function(th: np.ndarray) -> np.ndarray:
    return th[idx]


wald_result = wald_test(wald_function, theta, vcv, n_obs=data.shape[0])

print("Wald test -- no Granger-causality from income/consumption to investment:")
print("chi2 stat:", round(wald_result.chi2_stat, 4), " p-value:", round(wald_result.chi2_pvalue, 4))
print("F stat:   ", round(wald_result.f_stat, 4), " p-value:", round(wald_result.f_pvalue, 4))

Wald test for no instantaneous causality in a VAR — ects/varx1c.py

"""Translated from Examples/ects/varx1c.m -- Lutkepohl [1991], section
3.6: Wald test for no instantaneous causality from income/consumption to
investment, on the ML estimate of the VAR(2) model from varx1a.py/
varx1b.py's data. `varx_ls`/`varx_ml` map to `varx_estimate(...,
method="ls"|"ml")`; with method="ml" the returned `theta` is the
21-element stacked coefficient vector followed by the 6-element vech of
the Cholesky factor of Sigma (`varx_estimate`'s documented ML
convention), so theta[23]/theta[24] (1-indexed) fall within that
appended Cholesky block -- exactly as in the original."""

import io

import numpy as np

from quanttoolbox.econometrics.estimation import wald_test
from quanttoolbox.econometrics.var import varx_estimate

_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))

result = varx_estimate(y, x, p=2, method="ml")
theta = result.theta
vcv = result.vcv

# Test for no instantaneous causality from income/consumption to investment:
# theta[23] = theta[24] = 0 (1-indexed)
idx = np.array([23, 24]) - 1


def wald_function(th: np.ndarray) -> np.ndarray:
    return th[idx]


wald_result = wald_test(wald_function, theta, vcv, n_obs=data.shape[0])

print("Wald test -- no instantaneous causality from income/consumption to investment:")
print("chi2 stat:", round(wald_result.chi2_stat, 4), " p-value:", round(wald_result.chi2_pvalue, 4))
print("F stat:   ", round(wald_result.f_stat, 4), " p-value:", round(wald_result.f_pvalue, 4))

econometrics.var

Python alternatives

Switch to statsmodels.tsa.api.VAR for the unrestricted case — comprehensive lag-order selection, impulse response functions, forecast error variance decomposition, forecasting. Keep varx_estimate's linear-restriction support (a_eq/b_eq) — statsmodels' VAR doesn't support arbitrary parameter restrictions.

quanttoolbox.econometrics.var

VAR/VARX estimation (with exogenous regressors and linear parameter restrictions) and lag-order selection.

Ported from QuantToolBox/ects/{varx_cls,varx_cml,varx_ls,varx_ml, varx_order,var_constrained_estimation_onestep, varx_constrained_estimation_onestep}.m

Translation notes:

• var_constrained_estimation_onestep.m and varx_constrained_estimation_onestep.m are near-identical (the latter just returns a few extra result fields); only the superset (varx_estimate here) is ported, and var_* convenience wrappers become calls to varx_estimate with no exogenous regressors.
• varx_ls/varx_ml (thin dispatchers to the one-step estimator with a fixed identity starting covariance) and varx_cls (one-step LS with a user-supplied restriction) are not ported as separate functions -- call varx_estimate(..., method="ls"|"ml") directly.
• varx_cml (concentrated/iterated ML: re-estimate Sigma from the residuals and re-run one-step ML until Sigma stabilizes) is ported as varx_estimate_cml.
• The estimator itself uses the standard vec/GLS formulation: vec(theta) = (R'(Z Z' ⊗ Sigma¹)R)¹ R'(Z ⊗ Sigma¹) vec(Y), built on the vec/vech/duplication_matrix/elimination_matrix/ commutation_matrix helpers from quanttoolbox.linalg.special_matrices ported earlier.

varx_estimate(y, x=None, p=1, restriction=None, sigma=None, method='ls', compute_cov=True)

One-step VARX(p) estimation: y_t = Phi_1 y_{t-1} + ... + Phi_p y_{t-p} + beta @ x_t + u_t, with optional linear restrictions on the stacked coefficient vector theta = RR @ gamma + r, via GLS using the given (or identity) residual covariance Sigma.

method="ls" (default): residual covariance normalized by (n_usable - K*p - L) degrees of freedom. method="ml": normalized by n_usable (and the returned theta additionally includes the vech of the Cholesky factor of Sigma, matching the original's convention).

Original: ects/{var_constrained_estimation_onestep, varx_constrained_estimation_onestep}.m

Source code in src/quanttoolbox/econometrics/var.py

def varx_estimate(
    y: np.ndarray,
    x: np.ndarray | None = None,
    p: int = 1,
    restriction: tuple[np.ndarray, np.ndarray] | None = None,
    sigma: np.ndarray | None = None,
    method: str = "ls",
    compute_cov: bool = True,
) -> VARXResult:
    """One-step VARX(p) estimation: y_t = Phi_1 y_{t-1} + ... + Phi_p y_{t-p}
    + beta @ x_t + u_t, with optional linear restrictions on the stacked
    coefficient vector theta = RR @ gamma + r, via GLS using the given
    (or identity) residual covariance Sigma.

    method="ls" (default): residual covariance normalized by
    (n_usable - K*p - L) degrees of freedom. method="ml": normalized by
    n_usable (and the returned theta additionally includes the vech of
    the Cholesky factor of Sigma, matching the original's convention).

    Original: ects/{var_constrained_estimation_onestep,
    varx_constrained_estimation_onestep}.m
    """
    y = np.asarray(y, dtype=float)
    n_obs, k = y.shape

    if x is None:
        x_arr = np.zeros((n_obs, 1))
        n_exog = 0
    else:
        x_arr = np.asarray(x, dtype=float)
        if x_arr.ndim == 1:
            x_arr = x_arr[:, None]
        n_exog = x_arr.shape[1]

    n_params = k * (k * p + n_exog)
    rr, r = (np.eye(n_params), np.zeros(n_params)) if restriction is None else restriction
    rr = np.asarray(rr, dtype=float)
    r = np.asarray(r, dtype=float).flatten()

    if sigma is None:
        sigma_mat = np.eye(k)
    else:
        sigma_mat = np.asarray(sigma, dtype=float)

    valid = ~np.isnan(y).any(axis=1) & ~np.isnan(x_arr).any(axis=1)
    valid_idx = np.where(valid)[0]
    y_v = y[valid]
    x_v = x_arr[valid]
    n_valid = y_v.shape[0]
    n_usable = n_valid - p

    z = np.zeros((p * k, n_valid))
    for i in range(1, p + 1):
        w = np.full((n_valid, k), np.nan)
        w[i:] = y_v[:-i]
        z[(i - 1) * k : i * k, :] = w.T

    if n_exog > 0:
        z = np.vstack([z, x_v.T])

    y_stacked = y_v[p:].T
    z = z[:, p:]
    vec_y = y_stacked.flatten(order="F")

    inv_sigma = np.linalg.pinv(sigma_mat)
    y_c = np.kron(z.T, np.eye(k)) @ r
    y_c = vec_y - y_c

    w1 = z @ z.T
    w2 = np.kron(w1, inv_sigma)
    w3 = np.kron(z, inv_sigma)
    w4 = rr.T @ w2 @ rr
    w5 = np.linalg.pinv(w4)
    theta_c = w5 @ rr.T @ w3 @ y_c
    theta = rr @ theta_c + r

    bb = reshapec(theta, k, k * p + n_exog)
    phi = bb[:, : k * p]
    beta = bb[:, k * p :] if n_exog > 0 else None

    u = y_stacked - bb @ z

    residuals = np.full((n_obs, k), np.nan)
    residuals[valid_idx[p:]] = u.T

    sigma_hat = u @ u.T
    if method == "ml":
        sigma_hat = sigma_hat / n_usable
    else:
        sigma_hat = sigma_hat / (n_usable - k * p - n_exog)

    log_l = float(
        -0.5 * n_usable * np.log(np.linalg.det(sigma_hat))
        - 0.5 * n_usable * k * (np.log(2 * np.pi) + 1)
    )

    n_params_free = rr.shape[1]
    n_params_total = theta.shape[0]

    if method == "ml":
        chol = np.linalg.cholesky(sigma_hat)
        p_star = vech(chol, method=2)
        theta = np.concatenate([theta, p_star])
        n_params_total = theta.shape[0]

    if not compute_cov:
        return VARXResult(
            theta=theta,
            stderr=None,
            vcv=None,
            log_l=log_l,
            phi=phi,
            beta=beta,
            residuals=residuals,
            sigma=sigma_hat,
            n_obs=n_obs,
            n_obs_valid=n_valid,
            n_obs_usable=n_usable,
            n_params=n_params_total,
            n_params_free=n_params_free,
            df=n_usable - n_params_free,
            n_lags=p,
            k=k,
            l_exog=n_exog,
        )

    inv_sigma = np.linalg.pinv(sigma_hat)
    w2 = np.kron(w1, inv_sigma)
    w4 = rr.T @ w2 @ rr
    w5 = np.linalg.pinv(w4)
    vcv1 = rr @ w5 @ rr.T
    n_param1 = vcv1.shape[0]

    if method == "ml":
        d_k = duplication_matrix(k)
        l_k = elimination_matrix(k)
        k_kk = commutation_matrix(k, k)

        chol = np.linalg.cholesky(sigma_hat)
        h_mat = l_k @ (np.eye(k**2) + k_kk) @ np.kron(chol, np.eye(k)) @ l_k.T
        h_mat = np.linalg.pinv(h_mat)

        d_star = np.linalg.pinv(d_k.T @ d_k) @ d_k.T
        vcv2 = 2 * d_star @ np.kron(sigma_hat, sigma_hat) @ d_star.T
        vcv2 = vcv2 / n_usable
        vcv2 = h_mat @ vcv2 @ h_mat.T
        n_param2 = vcv2.shape[0]

        vcv = np.zeros((n_param1 + n_param2, n_param1 + n_param2))
        vcv[:n_param1, :n_param1] = vcv1
        vcv[n_param1:, n_param1:] = vcv2
    else:
        vcv = vcv1

    stderr = np.sqrt(np.diag(vcv))

    return VARXResult(
        theta=theta,
        stderr=stderr,
        vcv=vcv,
        log_l=log_l,
        phi=phi,
        beta=beta,
        residuals=residuals,
        sigma=sigma_hat,
        n_obs=n_obs,
        n_obs_valid=n_valid,
        n_obs_usable=n_usable,
        n_params=n_params_total,
        n_params_free=n_params_free,
        df=n_usable - n_params_free,
        n_lags=p,
        k=k,
        l_exog=n_exog,
    )

varx_estimate_cml(y, x=None, p=1, restriction=None, tol=1e-08, max_iters=100)

Concentrated (iterated) maximum likelihood VARX estimation: alternate one-step ML estimation and re-estimating Sigma from the residuals until det(Sigma) stabilizes.

Original: ects/varx_cml.m

Source code in src/quanttoolbox/econometrics/var.py

def varx_estimate_cml(
    y: np.ndarray,
    x: np.ndarray | None = None,
    p: int = 1,
    restriction: tuple[np.ndarray, np.ndarray] | None = None,
    tol: float = 1e-8,
    max_iters: int = 100,
) -> VARXResult:
    """Concentrated (iterated) maximum likelihood VARX estimation: alternate
    one-step ML estimation and re-estimating Sigma from the residuals
    until det(Sigma) stabilizes.

    Original: ects/varx_cml.m
    """
    y_arr = np.asarray(y, dtype=float)
    k = y_arr.shape[1]
    sigma = np.eye(k)

    for _ in range(max_iters):
        result = varx_estimate(y, x, p, restriction, sigma, method="ml", compute_cov=False)
        new_sigma = result.sigma
        if abs(np.linalg.det(sigma) - np.linalg.det(new_sigma)) < tol:
            sigma = new_sigma
            break
        sigma = new_sigma

    return varx_estimate(y, x, p, restriction, sigma, method="ml", compute_cov=True)

varx_order(y, x, p_max)

Select the VARX lag order by evaluating BIC, AICa (alpha=3), AICc, SIC, FPE, AIC, and HQ information criteria across p=0..p_max (or the explicit lag values in p_max, if given as an array).

Original: ects/varx_order.m

Source code in src/quanttoolbox/econometrics/var.py

def varx_order(y: np.ndarray, x: np.ndarray | None, p_max: int | np.ndarray) -> VARXOrderResult:
    """Select the VARX lag order by evaluating BIC, AICa (alpha=3), AICc,
    SIC, FPE, AIC, and HQ information criteria across p=0..p_max (or the
    explicit lag values in p_max, if given as an array).

    Original: ects/varx_order.m
    """
    if np.isscalar(p_max):
        p_values = np.arange(0, int(p_max) + 1)
    else:
        p_values = np.atleast_1d(p_max)
    n_lags = p_values.shape[0]
    criteria = np.zeros((n_lags, 7))

    for i, p in enumerate(p_values):
        result = varx_estimate(y, x, int(p), method="ml", compute_cov=False)
        t = result.n_obs_usable
        k = result.k
        log_det_sigma = np.log(np.linalg.det(result.sigma))

        criteria[i, 0] = log_det_sigma + p * k**2 * np.log(t) / t  # BIC
        criteria[i, 1] = log_det_sigma + 3 * p * k**2 / t  # AICa
        criteria[i, 2] = log_det_sigma + (1 + p * k**2 / t) / (1 - (p * k**2 - 2) / t)  # AICc
        criteria[i, 3] = log_det_sigma + k * np.log(1 + 2 * (p * k**2 + 1) / t)  # SIC
        criteria[i, 4] = log_det_sigma + k * np.log((t + p * k + 1) / (t - p * k - 1))  # FPE
        criteria[i, 5] = log_det_sigma + 2 * p * k**2 / t  # AIC
        criteria[i, 6] = log_det_sigma + 2 * p * k**2 * np.log(np.log(t)) / t  # HQ

    optimal_idx = np.argmin(criteria, axis=0)
    optimal_p = p_values[optimal_idx]

    return VARXOrderResult(p_values=p_values, criteria=criteria, optimal_p=optimal_p)

Examples

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))

Reduced form of a dynamic simultaneous-equations system — ects/varx2a.py

"""Translated from Examples/ects/varx2a.m -- Lutkepohl [1991], chapter
10: estimation of the reduced form (10.2.4) of a system of dynamic
simultaneous equations relating West German Income and Consumption to
their own lag and lagged Investment, on the same 76-row Lutkepohl.asc
data as varx1a.py (log-levels here, not log-differences).

`varx_ls(y, x, 1)` with `x = [ones(76,1) lag1(Investment)]` and p=1
means the regressors are the endogenous lag y_{t-1} = [Income(t-1),
Consumption(t-1)] (added automatically by `varx_estimate` for p=1) plus
the exogenous [constant, Investment(t-1)] -- exactly the 4-column layout
the original prints under the header "Inc(t-1)  Cons(t-1)  Constant
Inv(t-1)". `results.B` (not a field on `VARXResult`) is reconstructed as
`np.hstack([phi, beta])`, matching that same column order."""

import io

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate

_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])

result = varx_estimate(y, x, p=1, 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))

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))

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()

VAR lag-order selection — ects/varx1e.py

"""Translated from Examples/ects/varx1e.m -- Lutkepohl [1991], chapter 5:
lag-order selection (BIC/AICa/AICc/SIC/FPE/AIC/HQ) for a VAR model on the
same log-differenced West German investment/income/consumption data as
varx1a.py, over p = 1..5.

`varx_order(y, 1, seqa(1, 1, 5))` in the original again passes the scalar
`1` for the exogenous regressor `x`, i.e. "just a constant"; `seqa(1,
1, 5)` is the sequence 1, 2, 3, 4, 5 (`np.arange(1, 6)`)."""

import io

import numpy as np

from quanttoolbox.econometrics.var import varx_order

_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


d_investment = _diff1(investment)
d_income = _diff1(income)
d_consumption = _diff1(consumption)
y = np.column_stack([d_investment, d_income, d_consumption])

n = y.shape[0]
x = np.ones((n, 1))

result = varx_order(y, x, np.arange(1, 6))

print("p values:", result.p_values)
print("\ncriteria (columns: BIC, AICa, AICc, SIC, FPE, AIC, HQ):")
print(np.round(result.criteria, 4))
print("\noptimal p per criterion:", result.optimal_p)

VAR(0) equivalence to plain OLS — ects/varx3.py

"""Translated from Examples/ects/varx3.m -- a trivial 5-observation
example comparing plain OLS against `varx_estimate(..., p=0)` on the
same regressor matrix. With no autoregressive lags (p=0), VARX
estimation collapses to OLS, so `theta`/`beta` from the second call
should match `beta` from the first."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.econometrics.var import varx_estimate

y = np.array([2, 3, 1, 7, 5], dtype=float)

x = np.array(
    [
        [1, 3, 2],
        [2, 3, 1],
        [7, 1, 7],
        [5, 3, 1],
        [3, 5, 5],
    ],
    dtype=float,
)

ols_result = ols_estimation(y, np.column_stack([np.ones(5), x]))
print("OLS beta:", np.round(ols_result.beta, 4))

varx_result = varx_estimate(y[:, None], np.column_stack([np.ones(5), x]), p=0, method="ls")
print("\nVARX(p=0) theta:", np.round(varx_result.theta, 4))
print("VARX(p=0) beta:", np.round(varx_result.beta, 4))

VAR(2) on log-differenced macroeconomic data — ects/varx1a.py

"""Translated from Examples/ects/varx1a.m -- Lutkepohl [1991],
"Introduction to Multiple Time Series Analysis", section 3.2.3:
estimation of a VAR(2) process on log-differenced West German
investment/income/consumption data (1960Q1-1978Q4, the first 76 of the
92 rows in Lutkepohl.asc, embedded directly below per the project's
convention).

Note (preserved from the original): the constant is the FIRST column of
B in Lutkepohl's book, but ends up the LAST column of B here, since
`varx_estimate` places the autoregressive (Phi) columns before the
exogenous (beta) columns in its stacked coefficient layout -- the same
convention the original MATLAB program already followed.

`varx_ls(y, 1, 2)` in the original passes the scalar `1` as the
exogenous regressor `x`, which is the classic GAUSS/MATLAB-toolbox
shorthand for "just a constant" (equivalent to `ones(n, 1)`); that is
translated literally as `x = np.ones((n, 1))` below. `lag1` (first row
undefined/NaN, replaced by `varx_estimate`'s own NaN-row-dropping logic)
is translated as a plain first-difference with a leading NaN."""

import io

import numpy as np

from quanttoolbox.econometrics.var import varx_estimate

_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


d_investment = _diff1(investment)
d_income = _diff1(income)
d_consumption = _diff1(consumption)
y = np.column_stack([d_investment, d_income, d_consumption])

n = y.shape[0]
x = np.ones((n, 1))

result = varx_estimate(y, x, p=2, method="ls")

print("-" * 60, "Sigma", "-" * 1)
print(np.round(1e5 * result.sigma, 4))

b = np.hstack([result.phi, result.beta])
print("\n" + "-" * 60, "B", "-" * 1)
print(np.round(b, 4))

print("\n" + "-" * 60, "Phi", "-" * 1)
print(np.round(result.phi, 4))

print("\n" + "-" * 60, "beta", "-" * 1)
print(np.round(result.beta, 4))

Wald test for no Granger-causality between VAR variables — ects/varx1b.py

"""Translated from Examples/ects/varx1b.m -- Lutkepohl [1991], section
3.6: Wald test for no Granger-causality from income/consumption to
investment in the VAR(2) model estimated in varx1a.py (same
log-differenced Lutkepohl.asc data), i.e. H0: theta[4] = theta[7] =
theta[13] = theta[16] = 0 (1-indexed positions in the 21-element stacked
coefficient vector).

`RR = [design([4;7;13;16]) zeros(4,5)]` is only used to build the
constraint function `C(theta) = RR @ theta` -- it is not a restriction
passed to `varx_estimate` here, so it's translated as a direct 4-element
index/select into theta rather than via `design`."""

import io

import numpy as np

from quanttoolbox.econometrics.estimation import wald_test
from quanttoolbox.econometrics.var import varx_estimate

_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))

result = varx_estimate(y, x, p=2, method="ls")
theta = result.theta
vcv = result.vcv

# Test for no Granger-causality from income/consumption to investment:
# theta[4] = theta[7] = theta[13] = theta[16] = 0 (1-indexed)
idx = np.array([4, 7, 13, 16]) - 1


def wald_function(th: np.ndarray) -> np.ndarray:
    return th[idx]


wald_result = wald_test(wald_function, theta, vcv, n_obs=data.shape[0])

print("Wald test -- no Granger-causality from income/consumption to investment:")
print("chi2 stat:", round(wald_result.chi2_stat, 4), " p-value:", round(wald_result.chi2_pvalue, 4))
print("F stat:   ", round(wald_result.f_stat, 4), " p-value:", round(wald_result.f_pvalue, 4))

Wald test for no instantaneous causality in a VAR — ects/varx1c.py

"""Translated from Examples/ects/varx1c.m -- Lutkepohl [1991], section
3.6: Wald test for no instantaneous causality from income/consumption to
investment, on the ML estimate of the VAR(2) model from varx1a.py/
varx1b.py's data. `varx_ls`/`varx_ml` map to `varx_estimate(...,
method="ls"|"ml")`; with method="ml" the returned `theta` is the
21-element stacked coefficient vector followed by the 6-element vech of
the Cholesky factor of Sigma (`varx_estimate`'s documented ML
convention), so theta[23]/theta[24] (1-indexed) fall within that
appended Cholesky block -- exactly as in the original."""

import io

import numpy as np

from quanttoolbox.econometrics.estimation import wald_test
from quanttoolbox.econometrics.var import varx_estimate

_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))

result = varx_estimate(y, x, p=2, method="ml")
theta = result.theta
vcv = result.vcv

# Test for no instantaneous causality from income/consumption to investment:
# theta[23] = theta[24] = 0 (1-indexed)
idx = np.array([23, 24]) - 1


def wald_function(th: np.ndarray) -> np.ndarray:
    return th[idx]


wald_result = wald_test(wald_function, theta, vcv, n_obs=data.shape[0])

print("Wald test -- no instantaneous causality from income/consumption to investment:")
print("chi2 stat:", round(wald_result.chi2_stat, 4), " p-value:", round(wald_result.chi2_pvalue, 4))
print("F stat:   ", round(wald_result.f_stat, 4), " p-value:", round(wald_result.f_pvalue, 4))

econometrics.kalman

Python alternatives

Switch to statsmodels.tsa.statespace.MLEModel/kalman_filter for anything beyond simple filtering — smoothing, built-in MLE fitting, diffuse initialization, a compiled Cython backend. Keep this module for simple, transparent filtering or minimal-dependency use.

quanttoolbox.econometrics.kalman

Linear-Gaussian state-space models and the Kalman filter.

Ported from QuantToolBox/ects/{state_space_model,ssm_set,ssm_steady_state, Kalman_filtering}.m

Model convention (matching the original):

measurement:  y_t = Z_t @ a_t + d_t + eps_t,   eps_t ~ N(0, H_t)
transition:   a_t = T_t @ a_{t-1} + c_t + R_t @ eta_t,   eta_t ~ N(0, Q_t)

Translation notes:

• state_space_model/ssm_set (construct and validate a state-space model, optionally time-varying via 3-D arrays) are consolidated into the StateSpaceModel dataclass, whose __post_init__ performs the same dimension-consistency checks the original returned as a retcode flag.
• Kalman_filtering.m's main loop is ported directly; MATLAB's try/catch around the innovation-covariance inverse (falling back to Moore-Penrose pseudo-inverse if singular) is replicated exactly with numpy.linalg.inv/numpy.linalg.pinv.

StateSpaceModel(z, d, h, t, c, r, q) dataclass

Linear-Gaussian state-space model. Pass time-invariant matrices (2-D arrays) for a constant model, or time-varying matrices (3-D arrays, last axis = time) for a time-varying one.

Original: ects/state_space_model.m (+ ects/ssm_set.m for allocating empty time-varying matrix blocks -- just use np.zeros((n, m, nobs)) etc. directly in Python instead of a dedicated allocator function)

kalman_filter(ssm, y, a0, p0)

Run the Kalman filter for the given state-space model, observations y, and initial state mean/covariance (a0, P0).

Original: ects/Kalman_filtering.m

Source code in src/quanttoolbox/econometrics/kalman.py

def kalman_filter(
    ssm: StateSpaceModel, y: np.ndarray, a0: np.ndarray, p0: np.ndarray
) -> KalmanFilterResult:
    """Run the Kalman filter for the given state-space model, observations
    y, and initial state mean/covariance (a0, P0).

    Original: ects/Kalman_filtering.m
    """
    y = np.asarray(y, dtype=float)
    n, m = ssm.n, ssm.m
    nobs = ssm.nobs if ssm.nobs is not None else y.shape[0]

    if y.shape != (nobs, n):
        raise ValueError("kalman_filter: y has the wrong dimensions")

    if not ssm.time_varying:
        zt, dt, ht, tt, ct, rt, qt = ssm.z, ssm.d, ssm.h, ssm.t, ssm.c, ssm.r, ssm.q
        rqr = rt @ qt @ rt.T

    at, pt = np.asarray(a0, dtype=float).copy(), np.asarray(p0, dtype=float).copy()

    y_pred = np.zeros((nobs, n))
    v = np.zeros((nobs, n))
    f = np.zeros((n, n, nobs))
    a_pred = np.zeros((nobs, m))
    p_pred = np.zeros((m, m, nobs))
    a_filt = np.zeros((nobs, m))
    p_filt = np.zeros((m, m, nobs))
    log_l = np.zeros(nobs)

    for i in range(nobs):
        yt = y[i]

        if ssm.time_varying:
            zt, dt, ht = ssm.z[:, :, i], ssm.d[:, i], ssm.h[:, :, i]
            tt, ct, rt, qt = ssm.t[:, :, i], ssm.c[:, i], ssm.r[:, :, i], ssm.q[:, :, i]
            rqr = rt @ qt @ rt.T

        # prediction
        at1 = tt @ at + ct
        pt1 = tt @ pt @ tt.T + rqr

        # innovation
        yt1 = zt @ at1 + dt
        vt = yt - yt1

        # updating
        ft = zt @ pt1 @ zt.T + ht
        try:
            inv_ft = np.linalg.inv(ft)
        except np.linalg.LinAlgError:
            inv_ft = np.zeros((n, n)) if np.allclose(ft, 0) else np.linalg.pinv(ft)

        a_mat = pt1 @ zt.T
        b_mat = a_mat @ inv_ft

        at = at1 + b_mat @ vt
        pt = pt1 - b_mat @ a_mat.T

        det_ft = np.linalg.det(ft)
        if det_ft <= 0:
            det_ft = 1e-10

        y_pred[i] = yt1
        v[i] = vt
        f[:, :, i] = ft
        a_pred[i] = at1
        p_pred[:, :, i] = pt1
        a_filt[i] = at
        p_filt[:, :, i] = pt

        log_l[i] = -(n / 2) * np.log(2 * np.pi) - 0.5 * np.log(det_ft) - 0.5 * vt @ inv_ft @ vt

    return KalmanFilterResult(
        y_pred=y_pred,
        v=v,
        f=f,
        a_pred=a_pred,
        p_pred=p_pred,
        a_filt=a_filt,
        p_filt=p_filt,
        log_l=log_l,
    )

steady_state(ssm)

Steady-state (unconditional) mean and covariance of the state, for a time-invariant, stable state-space model: a_bar = (I-T)^-1 c, vec(P_bar) = (I - T ⊗ T)^-1 vec(R Q R').

Original: ects/ssm_steady_state.m

Source code in src/quanttoolbox/econometrics/kalman.py

def steady_state(ssm: StateSpaceModel) -> tuple[np.ndarray, np.ndarray]:
    """Steady-state (unconditional) mean and covariance of the state, for a
    time-invariant, stable state-space model: a_bar = (I-T)^-1 c,
    vec(P_bar) = (I - T ⊗ T)^-1 vec(R Q R').

    Original: ects/ssm_steady_state.m
    """
    if ssm.time_varying:
        raise ValueError("steady_state: model is time-varying")

    m = ssm.m
    w0 = np.eye(m) - ssm.t
    if np.allclose(w0, 0):
        raise ValueError("steady_state: the model is not stable (T == I)")
    try:
        w0_inv = np.linalg.inv(w0)
    except np.linalg.LinAlgError as exc:
        raise ValueError("steady_state: the model is not stable (I - T singular)") from exc

    a_bar = w0_inv @ ssm.c

    if np.allclose(ssm.q, 0):
        return a_bar, np.zeros((m, m))

    w1 = np.kron(ssm.t, ssm.t)
    w2 = np.eye(m * m) - w1
    try:
        w2_inv = np.linalg.inv(w2)
    except np.linalg.LinAlgError as exc:
        raise ValueError("steady_state: the model is not stable (I - T⊗T singular)") from exc

    w3 = w2_inv @ vec(ssm.r @ ssm.q @ ssm.r.T)
    p_bar = reshapec(w3, m, m)
    return a_bar, p_bar

Examples

Filtered state with a 95% confidence band — ects/kalman3c.py

"""Translated from Examples/ects/kalman3c.m -- same model and data as
kalman3a.py; extracts the filtered first state component alpha(t) =
a(t)[0] and its 95% confidence band alpha(t) +/- 1.96*sqrt(P(t)[0,0])
(numeric core only, the original's plot is dropped)."""

import io

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter, steady_state

_KALMAN3_DATA = """
12.004191367990 5.082325511092
15.665734306703 6.508713750061
13.009699927942 3.153318564843
15.335177481404 3.972550347263
12.508230936087 5.832925230371
14.085284431457 3.966144631579
13.880704275605 4.341577981101
17.299643540570 5.261307380309
12.247804130775 3.316255057591
13.432744976803 2.048672540611
11.622880573894 7.676694410981
11.788856472009 3.495128971097
12.917862956802 6.138211990221
14.317625056476 5.482944588328
16.847846044264 4.682035523793
11.967832572522 4.767366961079
14.974022160932 4.887977730878
15.016692818794 5.719026158573
14.424733030162 4.591756017103
16.042748794722 5.727057061524
14.885645813926 5.669898953266
14.755022585609 6.118419596523
13.536221951032 3.602206553287
11.596395495070 3.633464391525
11.604635094167 4.411888939754
12.390307519857 6.690385683215
15.305430497460 2.781199059932
15.402296198009 2.968476593783
16.168052982092 7.873695626866
15.981037574492 4.520455068664
19.740491145225 4.774787134276
17.774423955212 8.041963181636
16.375436108110 8.144223054978
19.684112696330 7.149207939792
18.032141720118 5.323073136597
11.342074743746 3.524220139858
15.349548300273 3.689291210812
18.349772674601 4.527378947789
15.530812985344 2.950467678878
17.225938859816 5.287363621403
16.579428949350 5.842125841512
16.351459570867 5.529566663511
16.770518442192 5.821908296493
14.106124772987 3.111896017411
15.726291282545 4.762045600136
17.699211885524 6.004068704830
17.634176507780 8.031641953772
14.932325325953 4.317652273000
15.296105656851 3.679922559759
13.240165786848 4.406086727283
17.378133569049 6.483437034615
17.546097415114 7.379225581477
16.276632663535 4.752764527387
14.051049656425 4.949933424766
14.994843414610 4.064828477724
11.049136259484 5.860639265902
18.042263203996 5.894551158525
17.083048154770 3.911844584559
12.821851106294 4.515592694496
14.945384650234 2.067169324453
14.714474843610 5.152085003822
14.110763908126 3.007624240753
13.982036948425 4.638848725578
11.939397400605 2.060874318059
14.067684713478 5.489049813033
14.631150348061 4.833850145498
16.059681645779 3.547100987793
16.607934346919 5.691142902747
16.445288786249 4.048342876977
15.500879280182 5.105351576079
13.597040366017 4.629353517066
15.379042268757 6.350990982664
15.308553895014 3.938070863008
19.443989876593 4.343043602184
15.106704622680 4.941879404882
14.905124865490 5.317340054023
16.437928609781 4.966803648347
13.449475918200 4.371207295415
14.742788936007 6.125527718659
15.335008061543 5.647248674467
15.737287573398 3.226021249963
14.236120737521 4.285589175699
14.522833146604 3.212970566847
12.261176341209 5.803765210148
14.510640575162 5.453329322308
16.696274014570 4.779056934387
17.445925211911 4.377813743756
12.684520208790 5.293449220792
15.365482189411 5.775810076543
15.144318098993 4.350539527521
15.530949743147 5.853153253807
13.539603397937 3.598410097115
14.161037667699 5.142596818775
14.355274468497 4.362791955901
11.799576339582 4.597257887517
13.564221566673 4.706566007960
17.016076203763 5.190421267339
16.480283255317 3.844473935740
12.416545487942 3.250372697960
13.752720516863 4.947745693575
"""

y = np.loadtxt(io.StringIO(_KALMAN3_DATA))
nobs = y.shape[0]

ssm = StateSpaceModel(
    z=np.eye(2),
    d=np.array([10.0, 0.0]),
    h=np.array([[2.0, 0.0], [0.0, 1.0]]),
    t=np.array([[0.5, 0.3], [0.0, 0.2]]),
    c=np.array([1.0, 4.0]),
    r=np.array([[1.0], [1.0]]),
    q=np.array([[1.0]]),
)

a_bar, p_bar = steady_state(ssm)
result = kalman_filter(ssm, y, a_bar, p_bar)

alpha = result.a_filt[:, 0]
cov = result.p_filt[0, 0, :]
lower = alpha - 1.96 * np.sqrt(cov)
upper = alpha + 1.96 * np.sqrt(cov)

t = np.arange(1, nobs + 1)
print("t, alpha(t), lower 95% bound, upper 95% bound -- first/last 10 observations:")
print(np.round(np.column_stack([t, alpha, lower, upper])[:10], 4))
print(np.round(np.column_stack([t, alpha, lower, upper])[-10:], 4))

General state-space filter from a steady-state start — ects/kalman3a.py

"""Translated from Examples/ects/kalman3a.m -- runs the Kalman filter for
a general fixed, time-invariant 2-state / 2-observation state-space
model on the 100-observation Kalman3.asc series (embedded directly
below), starting from the steady-state mean/covariance of the state,
and prints the filtered state a(t) and its covariance P(t) for the
first 10 observations."""

import io

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter, steady_state

_KALMAN3_DATA = """
12.004191367990 5.082325511092
15.665734306703 6.508713750061
13.009699927942 3.153318564843
15.335177481404 3.972550347263
12.508230936087 5.832925230371
14.085284431457 3.966144631579
13.880704275605 4.341577981101
17.299643540570 5.261307380309
12.247804130775 3.316255057591
13.432744976803 2.048672540611
11.622880573894 7.676694410981
11.788856472009 3.495128971097
12.917862956802 6.138211990221
14.317625056476 5.482944588328
16.847846044264 4.682035523793
11.967832572522 4.767366961079
14.974022160932 4.887977730878
15.016692818794 5.719026158573
14.424733030162 4.591756017103
16.042748794722 5.727057061524
14.885645813926 5.669898953266
14.755022585609 6.118419596523
13.536221951032 3.602206553287
11.596395495070 3.633464391525
11.604635094167 4.411888939754
12.390307519857 6.690385683215
15.305430497460 2.781199059932
15.402296198009 2.968476593783
16.168052982092 7.873695626866
15.981037574492 4.520455068664
19.740491145225 4.774787134276
17.774423955212 8.041963181636
16.375436108110 8.144223054978
19.684112696330 7.149207939792
18.032141720118 5.323073136597
11.342074743746 3.524220139858
15.349548300273 3.689291210812
18.349772674601 4.527378947789
15.530812985344 2.950467678878
17.225938859816 5.287363621403
16.579428949350 5.842125841512
16.351459570867 5.529566663511
16.770518442192 5.821908296493
14.106124772987 3.111896017411
15.726291282545 4.762045600136
17.699211885524 6.004068704830
17.634176507780 8.031641953772
14.932325325953 4.317652273000
15.296105656851 3.679922559759
13.240165786848 4.406086727283
17.378133569049 6.483437034615
17.546097415114 7.379225581477
16.276632663535 4.752764527387
14.051049656425 4.949933424766
14.994843414610 4.064828477724
11.049136259484 5.860639265902
18.042263203996 5.894551158525
17.083048154770 3.911844584559
12.821851106294 4.515592694496
14.945384650234 2.067169324453
14.714474843610 5.152085003822
14.110763908126 3.007624240753
13.982036948425 4.638848725578
11.939397400605 2.060874318059
14.067684713478 5.489049813033
14.631150348061 4.833850145498
16.059681645779 3.547100987793
16.607934346919 5.691142902747
16.445288786249 4.048342876977
15.500879280182 5.105351576079
13.597040366017 4.629353517066
15.379042268757 6.350990982664
15.308553895014 3.938070863008
19.443989876593 4.343043602184
15.106704622680 4.941879404882
14.905124865490 5.317340054023
16.437928609781 4.966803648347
13.449475918200 4.371207295415
14.742788936007 6.125527718659
15.335008061543 5.647248674467
15.737287573398 3.226021249963
14.236120737521 4.285589175699
14.522833146604 3.212970566847
12.261176341209 5.803765210148
14.510640575162 5.453329322308
16.696274014570 4.779056934387
17.445925211911 4.377813743756
12.684520208790 5.293449220792
15.365482189411 5.775810076543
15.144318098993 4.350539527521
15.530949743147 5.853153253807
13.539603397937 3.598410097115
14.161037667699 5.142596818775
14.355274468497 4.362791955901
11.799576339582 4.597257887517
13.564221566673 4.706566007960
17.016076203763 5.190421267339
16.480283255317 3.844473935740
12.416545487942 3.250372697960
13.752720516863 4.947745693575
"""

y = np.loadtxt(io.StringIO(_KALMAN3_DATA))
nobs = y.shape[0]

ssm = StateSpaceModel(
    z=np.eye(2),
    d=np.array([10.0, 0.0]),
    h=np.array([[2.0, 0.0], [0.0, 1.0]]),
    t=np.array([[0.5, 0.3], [0.0, 0.2]]),
    c=np.array([1.0, 4.0]),
    r=np.array([[1.0], [1.0]]),
    q=np.array([[1.0]]),
)

a_bar, p_bar = steady_state(ssm)
result = kalman_filter(ssm, y, a_bar, p_bar)

at = result.a_filt
pt = result.p_filt

for t in range(10):
    print("-" * 60, f"t = {t + 1}")
    print("a(t) =")
    print(np.round(at[t], 7))
    print("P(t) =")
    print(np.round(pt[:, :, t], 7))
    print()

Local Linear Trend Kalman filter, fixed variance parameters — ects/kalman2a.py

"""Translated from Examples/ects/kalman2a.m -- Harvey [1990],
"Forecasting, Structural Time Series and the Kalman Filter", pages
89-90: runs the Kalman filter for a Local Linear Trend (LLT) model
(state = [level, slope]) on 61 annual US GNP observations (1909-1969,
Gnp.asc, embedded directly below), with fixed
(sigma_epsilon, sigma_eta1, sigma_eta2) = (1, 0.5, 0.5) -- numeric core
only, the original's dual-axes plot (level + slope) is dropped.

`results.a_cond` (the filtered-ahead state a(t|t-1), used for the level
sub-plot) maps to `result.a_pred`; `results.a` (the filtered state
a(t|t), computed but only referenced for completeness in the original)
maps to `result.a_filt`."""

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

y = np.array(
    [
        116.8,
        120.0,
        123.2,
        130.2,
        131.4,
        125.6,
        124.5,
        134.3,
        135.2,
        151.8,
        146.4,
        139.0,
        127.8,
        147.0,
        165.9,
        165.5,
        179.4,
        190.0,
        189.8,
        190.9,
        203.6,
        183.5,
        169.3,
        144.2,
        141.5,
        154.3,
        169.5,
        193.0,
        203.2,
        192.9,
        209.4,
        227.2,
        263.7,
        297.8,
        337.1,
        361.3,
        355.2,
        312.6,
        309.9,
        323.7,
        324.1,
        255.3,
        383.4,
        395.1,
        412.8,
        406.0,
        438.0,
        446.1,
        452.5,
        447.3,
        475.9,
        487.7,
        497.2,
        529.8,
        551.0,
        581.1,
        617.8,
        658.1,
        675.2,
        706.6,
        724.7,
    ],
    dtype=float,
)
nobs = y.shape[0]

sigma_epsilon = 1.0
sigma_eta1 = 0.5
sigma_eta2 = 0.5

ssm = StateSpaceModel(
    z=np.array([[1.0, 0.0]]),
    d=np.array([0.0]),
    h=np.array([[sigma_epsilon**2]]),
    t=np.array([[1.0, 1.0], [0.0, 1.0]]),
    c=np.array([0.0, 0.0]),
    r=np.eye(2),
    q=np.diag([sigma_eta1**2, sigma_eta2**2]),
)

a0 = np.array([y[0], 0.0])
p0 = np.zeros((2, 2))
result = kalman_filter(ssm, y[:, None], a0, p0)

t = np.arange(1909, 1909 + nobs)
print("t, y, level a(t|t-1), slope a(t|t-1) -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, result.a_pred])[:10], 3))
print(np.round(np.column_stack([t, y, result.a_pred])[-10:], 3))
print("\nsum log-likelihood:", round(float(np.sum(result.log_l)), 4))

Local Linear Trend model fit by time-domain ML, vs. Whittle — ects/kalman2c.py

"""Translated from Examples/ects/kalman2c.m -- Harvey [1990], pages
89-90: time-domain (exact) maximum-likelihood estimation of the Local
Linear Trend model's (sigma_epsilon, sigma_eta, sigma_zeta) on the same
61-observation Gnp.asc series as kalman2a.py/kalman2b.py, then runs the
Kalman filter with the estimated parameters (numeric core only, plot
dropped).

As in kalman1b.py/kalman1c.py, the original's `theta = sqrt(theta.^2)`
inside `LLT_ml` is translated as `np.abs(theta)`."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

y = np.array(
    [
        116.8,
        120.0,
        123.2,
        130.2,
        131.4,
        125.6,
        124.5,
        134.3,
        135.2,
        151.8,
        146.4,
        139.0,
        127.8,
        147.0,
        165.9,
        165.5,
        179.4,
        190.0,
        189.8,
        190.9,
        203.6,
        183.5,
        169.3,
        144.2,
        141.5,
        154.3,
        169.5,
        193.0,
        203.2,
        192.9,
        209.4,
        227.2,
        263.7,
        297.8,
        337.1,
        361.3,
        355.2,
        312.6,
        309.9,
        323.7,
        324.1,
        255.3,
        383.4,
        395.1,
        412.8,
        406.0,
        438.0,
        446.1,
        452.5,
        447.3,
        475.9,
        487.7,
        497.2,
        529.8,
        551.0,
        581.1,
        617.8,
        658.1,
        675.2,
        706.6,
        724.7,
    ],
    dtype=float,
)
nobs = y.shape[0]

a0 = np.array([y[0], 0.0])
p0 = np.zeros((2, 2))


def _ssm(sigma_epsilon: float, sigma_eta: float, sigma_zeta: float) -> StateSpaceModel:
    return StateSpaceModel(
        z=np.array([[1.0, 0.0]]),
        d=np.array([0.0]),
        h=np.array([[sigma_epsilon**2]]),
        t=np.array([[1.0, 1.0], [0.0, 1.0]]),
        c=np.array([0.0, 0.0]),
        r=np.eye(2),
        q=np.diag([sigma_eta**2, sigma_zeta**2]),
    )


def llt_ml(theta: np.ndarray) -> np.ndarray:
    sigma_epsilon, sigma_eta, sigma_zeta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta, sigma_zeta)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    return result.log_l


sv = 3.0 * np.ones(3)
ml_result = ml_estimation(llt_ml, sv)
theta = ml_result.theta

sigma_epsilon, sigma_eta, sigma_zeta = np.abs(theta)
print("theta (sigma_epsilon, sigma_eta, sigma_zeta):", np.round(theta, 4))
print("log-likelihood:", round(ml_result.sum_log_l, 4))

ssm = _ssm(sigma_epsilon, sigma_eta, sigma_zeta)
result = kalman_filter(ssm, y[:, None], a0, p0)

t = np.arange(1909, 1909 + nobs)
print("\nt, y, level a(t|t-1), slope a(t|t-1) -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, result.a_pred])[:10], 3))
print(np.round(np.column_stack([t, y, result.a_pred])[-10:], 3))

Local Linear Trend model fit by Whittle (frequency-domain) ML — ects/kalman2b.py

"""Translated from Examples/ects/kalman2b.m -- Harvey [1990], pages
89-90: frequency-domain (Whittle) maximum-likelihood estimation of the
Local Linear Trend model's (sigma_epsilon, sigma_eta, sigma_zeta) on the
same 61-observation Gnp.asc series as kalman2a.py, then runs the Kalman
filter with the estimated parameters (numeric core only, plot dropped).

The original calls `whittle_local_linear_trend` twice -- once with
`WHITTLE_algorithm = 1` (BFGS) and once with `WHITTLE_algorithm = 3`
("scoring method") -- keeping only the second result. As established in
`whittle1.py`/`whittle.py`'s module docstring, the port only has a single
optimizer path (BFGS), so there is just one call here."""

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter
from quanttoolbox.econometrics.whittle import whittle_local_linear_trend

y = np.array(
    [
        116.8,
        120.0,
        123.2,
        130.2,
        131.4,
        125.6,
        124.5,
        134.3,
        135.2,
        151.8,
        146.4,
        139.0,
        127.8,
        147.0,
        165.9,
        165.5,
        179.4,
        190.0,
        189.8,
        190.9,
        203.6,
        183.5,
        169.3,
        144.2,
        141.5,
        154.3,
        169.5,
        193.0,
        203.2,
        192.9,
        209.4,
        227.2,
        263.7,
        297.8,
        337.1,
        361.3,
        355.2,
        312.6,
        309.9,
        323.7,
        324.1,
        255.3,
        383.4,
        395.1,
        412.8,
        406.0,
        438.0,
        446.1,
        452.5,
        447.3,
        475.9,
        487.7,
        497.2,
        529.8,
        551.0,
        581.1,
        617.8,
        658.1,
        675.2,
        706.6,
        724.7,
    ],
    dtype=float,
)
nobs = y.shape[0]

# BFGS method (very sensitive to the starting values -- kept as-is)
sv = 15.0 * np.ones(3)
result_est = whittle_local_linear_trend(y, sv)
theta = result_est.theta

sigma_epsilon, sigma_eta, sigma_zeta = theta
print("theta (sigma_epsilon, sigma_eta, sigma_zeta):", np.round(theta, 4))

ssm = StateSpaceModel(
    z=np.array([[1.0, 0.0]]),
    d=np.array([0.0]),
    h=np.array([[sigma_epsilon**2]]),
    t=np.array([[1.0, 1.0], [0.0, 1.0]]),
    c=np.array([0.0, 0.0]),
    r=np.eye(2),
    q=np.diag([sigma_eta**2, sigma_zeta**2]),
)

a0 = np.array([y[0], 0.0])
p0 = np.zeros((2, 2))
result = kalman_filter(ssm, y[:, None], a0, p0)

t = np.arange(1909, 1909 + nobs)
print("\nt, y, level a(t|t-1), slope a(t|t-1) -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, result.a_pred])[:10], 3))
print(np.round(np.column_stack([t, y, result.a_pred])[-10:], 3))

Local-level Kalman filter (Harvey 1990) — ects/panel1.py

"""Translated from Examples/ects/panel1.m -- despite its filename and the
tracker's earlier note ("panel-data example; no ported panel module
yet"), this is actually a local-level (random-walk-plus-noise) Kalman
filter example (Harvey 1990, pp.89-90) applied to a 71-observation
"Purse" series -- it doesn't touch any panel-data functionality. Since
`quanttoolbox.econometrics.kalman` *is* ported, this is translated
properly rather than left as a gap; the tracker entry is corrected
accordingly.

The 71 observations from Purse.asc are embedded directly below, matching
the project's convention of self-contained example scripts (numeric core
only; the original's plot of y_t vs. the one-step-ahead prediction is
dropped)."""

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

y = np.array(
    [
        10,
        15,
        10,
        10,
        12,
        10,
        7,
        17,
        10,
        14,
        8,
        17,
        14,
        18,
        3,
        9,
        11,
        10,
        6,
        12,
        14,
        10,
        25,
        29,
        33,
        33,
        12,
        19,
        16,
        19,
        19,
        12,
        34,
        15,
        36,
        29,
        26,
        21,
        17,
        19,
        13,
        20,
        24,
        12,
        6,
        14,
        6,
        12,
        9,
        11,
        17,
        12,
        8,
        14,
        14,
        12,
        5,
        8,
        10,
        3,
        16,
        8,
        8,
        7,
        12,
        6,
        10,
        8,
        10,
        5,
        7,
    ],
    dtype=float,
)
nobs = y.shape[0]

sigma_epsilon = 5.0
sigma_eta = 2.15

ssm = StateSpaceModel(
    z=np.array([[1.0]]),
    d=np.array([0.0]),
    h=np.array([[sigma_epsilon**2]]),
    t=np.array([[1.0]]),
    c=np.array([0.0]),
    r=np.array([[1.0]]),
    q=np.array([[sigma_eta**2]]),
)

a0 = np.array([y[0]])
p0 = np.array([[0.0]])

result = kalman_filter(ssm, y[:, None], a0, p0)
y_cond = result.y_pred[:, 0]

t = np.arange(nobs)
print("t, y, y(t|t-1) -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, y_cond])[:10], 3))
print(np.round(np.column_stack([t, y, y_cond])[-10:], 3))
print("\nsum log-likelihood:", round(float(np.sum(result.log_l)), 4))

ML recovery of all free state-space model parameters — ects/kalman3d.py

"""Translated from Examples/ects/kalman3d.m -- maximum-likelihood
recovery of all 11 free parameters (Z diagonal, d[0], H diagonal, T,
c, Q) of the general state-space model used in kalman3a.py/kalman3b.py/
kalman3c.py, starting the optimizer from the model's own true parameter
values (as the original does) and comparing the recovered estimates
against them.

`theta([4 5 11]) = sqrt(theta([4 5 11]).^2)` in the original enforces
H[0,0], H[1,1], Q > 0 (indices 4, 5, 11 in MATLAB's 1-based numbering);
translated as `np.abs(theta[[3, 4, 10]])` here. `a0`/`P0` are
recomputed from the model's own steady state at every trial `theta`,
exactly as in `ssm_ml_fun`."""

import io

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter, steady_state

_KALMAN3_DATA = """
12.004191367990 5.082325511092
15.665734306703 6.508713750061
13.009699927942 3.153318564843
15.335177481404 3.972550347263
12.508230936087 5.832925230371
14.085284431457 3.966144631579
13.880704275605 4.341577981101
17.299643540570 5.261307380309
12.247804130775 3.316255057591
13.432744976803 2.048672540611
11.622880573894 7.676694410981
11.788856472009 3.495128971097
12.917862956802 6.138211990221
14.317625056476 5.482944588328
16.847846044264 4.682035523793
11.967832572522 4.767366961079
14.974022160932 4.887977730878
15.016692818794 5.719026158573
14.424733030162 4.591756017103
16.042748794722 5.727057061524
14.885645813926 5.669898953266
14.755022585609 6.118419596523
13.536221951032 3.602206553287
11.596395495070 3.633464391525
11.604635094167 4.411888939754
12.390307519857 6.690385683215
15.305430497460 2.781199059932
15.402296198009 2.968476593783
16.168052982092 7.873695626866
15.981037574492 4.520455068664
19.740491145225 4.774787134276
17.774423955212 8.041963181636
16.375436108110 8.144223054978
19.684112696330 7.149207939792
18.032141720118 5.323073136597
11.342074743746 3.524220139858
15.349548300273 3.689291210812
18.349772674601 4.527378947789
15.530812985344 2.950467678878
17.225938859816 5.287363621403
16.579428949350 5.842125841512
16.351459570867 5.529566663511
16.770518442192 5.821908296493
14.106124772987 3.111896017411
15.726291282545 4.762045600136
17.699211885524 6.004068704830
17.634176507780 8.031641953772
14.932325325953 4.317652273000
15.296105656851 3.679922559759
13.240165786848 4.406086727283
17.378133569049 6.483437034615
17.546097415114 7.379225581477
16.276632663535 4.752764527387
14.051049656425 4.949933424766
14.994843414610 4.064828477724
11.049136259484 5.860639265902
18.042263203996 5.894551158525
17.083048154770 3.911844584559
12.821851106294 4.515592694496
14.945384650234 2.067169324453
14.714474843610 5.152085003822
14.110763908126 3.007624240753
13.982036948425 4.638848725578
11.939397400605 2.060874318059
14.067684713478 5.489049813033
14.631150348061 4.833850145498
16.059681645779 3.547100987793
16.607934346919 5.691142902747
16.445288786249 4.048342876977
15.500879280182 5.105351576079
13.597040366017 4.629353517066
15.379042268757 6.350990982664
15.308553895014 3.938070863008
19.443989876593 4.343043602184
15.106704622680 4.941879404882
14.905124865490 5.317340054023
16.437928609781 4.966803648347
13.449475918200 4.371207295415
14.742788936007 6.125527718659
15.335008061543 5.647248674467
15.737287573398 3.226021249963
14.236120737521 4.285589175699
14.522833146604 3.212970566847
12.261176341209 5.803765210148
14.510640575162 5.453329322308
16.696274014570 4.779056934387
17.445925211911 4.377813743756
12.684520208790 5.293449220792
15.365482189411 5.775810076543
15.144318098993 4.350539527521
15.530949743147 5.853153253807
13.539603397937 3.598410097115
14.161037667699 5.142596818775
14.355274468497 4.362791955901
11.799576339582 4.597257887517
13.564221566673 4.706566007960
17.016076203763 5.190421267339
16.480283255317 3.844473935740
12.416545487942 3.250372697960
13.752720516863 4.947745693575
"""

y = np.loadtxt(io.StringIO(_KALMAN3_DATA))

ssm0 = StateSpaceModel(
    z=np.eye(2),
    d=np.array([10.0, 0.0]),
    h=np.array([[2.0, 0.0], [0.0, 1.0]]),
    t=np.array([[0.5, 0.3], [0.0, 0.2]]),
    c=np.array([1.0, 4.0]),
    r=np.array([[1.0], [1.0]]),
    q=np.array([[1.0]]),
)
a_bar0, p_bar0 = steady_state(ssm0)
result0 = kalman_filter(ssm0, y, a_bar0, p_bar0)
print("Value of the log-likelihood function:", round(float(np.sum(result0.log_l)), 4))


def ssm_ml_fun(theta: np.ndarray) -> np.ndarray:
    theta = theta.copy()
    theta[[3, 4, 10]] = np.abs(theta[[3, 4, 10]])

    ssm = StateSpaceModel(
        z=np.array([[theta[0], 0.0], [0.0, theta[1]]]),
        d=np.array([theta[2], 0.0]),
        h=np.array([[theta[3], 0.0], [0.0, theta[4]]]),
        t=np.array([[theta[5], theta[6]], [0.0, theta[7]]]),
        c=np.array([theta[8], theta[9]]),
        r=np.array([[1.0], [1.0]]),
        q=np.array([[theta[10]]]),
    )
    a_bar, p_bar = steady_state(ssm)
    result = kalman_filter(ssm, y, a_bar, p_bar)
    return result.log_l


sv = np.array([1.0, 1.0, 10.0, 2.0, 1.0, 0.5, 0.3, 0.2, 1.0, 4.0, 1.0])

ml_result = ml_estimation(ssm_ml_fun, sv)
theta_hat = ml_result.theta

print("\n True           Estimated")
print(" values         values")
print(np.round(np.column_stack([sv, theta_hat]), 5))

Same state-space filter via the time-varying code path — ects/kalman3b.py

"""Translated from Examples/ects/kalman3b.m -- same model and data as
kalman3a.py, but broadcasts the (constant) system matrices out to
explicit per-period 3-D arrays and runs the Kalman filter through the
time-varying code path, to confirm it reproduces the time-invariant
path's results exactly (the whole point of the original's exercise)."""

import io

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter, steady_state

_KALMAN3_DATA = """
12.004191367990 5.082325511092
15.665734306703 6.508713750061
13.009699927942 3.153318564843
15.335177481404 3.972550347263
12.508230936087 5.832925230371
14.085284431457 3.966144631579
13.880704275605 4.341577981101
17.299643540570 5.261307380309
12.247804130775 3.316255057591
13.432744976803 2.048672540611
11.622880573894 7.676694410981
11.788856472009 3.495128971097
12.917862956802 6.138211990221
14.317625056476 5.482944588328
16.847846044264 4.682035523793
11.967832572522 4.767366961079
14.974022160932 4.887977730878
15.016692818794 5.719026158573
14.424733030162 4.591756017103
16.042748794722 5.727057061524
14.885645813926 5.669898953266
14.755022585609 6.118419596523
13.536221951032 3.602206553287
11.596395495070 3.633464391525
11.604635094167 4.411888939754
12.390307519857 6.690385683215
15.305430497460 2.781199059932
15.402296198009 2.968476593783
16.168052982092 7.873695626866
15.981037574492 4.520455068664
19.740491145225 4.774787134276
17.774423955212 8.041963181636
16.375436108110 8.144223054978
19.684112696330 7.149207939792
18.032141720118 5.323073136597
11.342074743746 3.524220139858
15.349548300273 3.689291210812
18.349772674601 4.527378947789
15.530812985344 2.950467678878
17.225938859816 5.287363621403
16.579428949350 5.842125841512
16.351459570867 5.529566663511
16.770518442192 5.821908296493
14.106124772987 3.111896017411
15.726291282545 4.762045600136
17.699211885524 6.004068704830
17.634176507780 8.031641953772
14.932325325953 4.317652273000
15.296105656851 3.679922559759
13.240165786848 4.406086727283
17.378133569049 6.483437034615
17.546097415114 7.379225581477
16.276632663535 4.752764527387
14.051049656425 4.949933424766
14.994843414610 4.064828477724
11.049136259484 5.860639265902
18.042263203996 5.894551158525
17.083048154770 3.911844584559
12.821851106294 4.515592694496
14.945384650234 2.067169324453
14.714474843610 5.152085003822
14.110763908126 3.007624240753
13.982036948425 4.638848725578
11.939397400605 2.060874318059
14.067684713478 5.489049813033
14.631150348061 4.833850145498
16.059681645779 3.547100987793
16.607934346919 5.691142902747
16.445288786249 4.048342876977
15.500879280182 5.105351576079
13.597040366017 4.629353517066
15.379042268757 6.350990982664
15.308553895014 3.938070863008
19.443989876593 4.343043602184
15.106704622680 4.941879404882
14.905124865490 5.317340054023
16.437928609781 4.966803648347
13.449475918200 4.371207295415
14.742788936007 6.125527718659
15.335008061543 5.647248674467
15.737287573398 3.226021249963
14.236120737521 4.285589175699
14.522833146604 3.212970566847
12.261176341209 5.803765210148
14.510640575162 5.453329322308
16.696274014570 4.779056934387
17.445925211911 4.377813743756
12.684520208790 5.293449220792
15.365482189411 5.775810076543
15.144318098993 4.350539527521
15.530949743147 5.853153253807
13.539603397937 3.598410097115
14.161037667699 5.142596818775
14.355274468497 4.362791955901
11.799576339582 4.597257887517
13.564221566673 4.706566007960
17.016076203763 5.190421267339
16.480283255317 3.844473935740
12.416545487942 3.250372697960
13.752720516863 4.947745693575
"""

y = np.loadtxt(io.StringIO(_KALMAN3_DATA))
nobs = y.shape[0]

z_t = np.eye(2)
d_t = np.array([10.0, 0.0])
h_t = np.array([[2.0, 0.0], [0.0, 1.0]])
t_t = np.array([[0.5, 0.3], [0.0, 0.2]])
c_t = np.array([1.0, 4.0])
r_t = np.array([[1.0], [1.0]])
q_t = np.array([[1.0]])

ssm_invariant = StateSpaceModel(z=z_t, d=d_t, h=h_t, t=t_t, c=c_t, r=r_t, q=q_t)

ssm_variant = StateSpaceModel(
    z=np.repeat(z_t[:, :, None], nobs, axis=2),
    d=np.repeat(d_t[:, None], nobs, axis=1),
    h=np.repeat(h_t[:, :, None], nobs, axis=2),
    t=np.repeat(t_t[:, :, None], nobs, axis=2),
    c=np.repeat(c_t[:, None], nobs, axis=1),
    r=np.repeat(r_t[:, :, None], nobs, axis=2),
    q=np.repeat(q_t[:, :, None], nobs, axis=2),
)

a_bar, p_bar = steady_state(ssm_invariant)
result = kalman_filter(ssm_variant, y, a_bar, p_bar)

at = result.a_filt
pt = result.p_filt

for t in range(10):
    print("-" * 60, f"t = {t + 1}")
    print("a(t) =")
    print(np.round(at[t], 7))
    print("P(t) =")
    print(np.round(pt[:, :, t], 7))
    print()

Time-domain ML of Kalman noise variances, two starting states — ects/kalman1b.py

"""Translated from Examples/ects/kalman1b.m -- Harvey [1990],
"Forecasting, Structural Time Series and the Kalman Filter", pages
89-90: time-domain (exact) maximum-likelihood estimation of the
local-level model's (sigma_epsilon, sigma_eta) on the same 71-
observation Purse.asc series as panel1.py/whittle1.py, comparing two
choices of the Kalman filter's initial state a0 used *during
estimation* (0, vs. y[0]) -- then re-running the filter with each
estimated theta (using a0 = y[0], P0 = 0 for both, matching the
original's separate post-estimation filter block) to compare the
resulting one-step-ahead predictions.

The original's `theta = sqrt(theta.^2)` inside `LL_ml` (a roundabout way
of enforcing sigma_epsilon, sigma_eta > 0 while keeping the objective
smooth for the optimizer) is translated literally as `np.abs(theta)`.
`ml_estimation`'s `logpdf_fn` here returns the Kalman filter's
per-observation `log_l` array (the natural per-observation log-density
for this model), rather than a single pre-summed scalar (which
`ml_estimation` would also accept, but only as a degenerate 1-observation
case)."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

y = np.array(
    [
        10,
        15,
        10,
        10,
        12,
        10,
        7,
        17,
        10,
        14,
        8,
        17,
        14,
        18,
        3,
        9,
        11,
        10,
        6,
        12,
        14,
        10,
        25,
        29,
        33,
        33,
        12,
        19,
        16,
        19,
        19,
        12,
        34,
        15,
        36,
        29,
        26,
        21,
        17,
        19,
        13,
        20,
        24,
        12,
        6,
        14,
        6,
        12,
        9,
        11,
        17,
        12,
        8,
        14,
        14,
        12,
        5,
        8,
        10,
        3,
        16,
        8,
        8,
        7,
        12,
        6,
        10,
        8,
        10,
        5,
        7,
    ],
    dtype=float,
)
nobs = y.shape[0]


def _ssm(sigma_epsilon: float, sigma_eta: float) -> StateSpaceModel:
    return StateSpaceModel(
        z=np.array([[1.0]]),
        d=np.array([0.0]),
        h=np.array([[sigma_epsilon**2]]),
        t=np.array([[1.0]]),
        c=np.array([0.0]),
        r=np.array([[1.0]]),
        q=np.array([[sigma_eta**2]]),
    )


def ll_ml(theta: np.ndarray, a0: float) -> np.ndarray:
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], np.array([a0]), np.array([[0.0]]))
    return result.log_l


sv = np.array([3.0, 1.0])

theta1 = ml_estimation(lambda theta: ll_ml(theta, 0.0), sv).theta
theta2 = ml_estimation(lambda theta: ll_ml(theta, y[0]), sv).theta

print("theta1 (sigma_epsilon, sigma_eta), a0=0 during estimation:", np.round(theta1, 4))
print("theta2 (sigma_epsilon, sigma_eta), a0=y[0] during estimation:", np.round(theta2, 4))

a0 = np.array([y[0]])
p0 = np.array([[0.0]])

y_cond = np.zeros((nobs, 2))
for i, theta in enumerate((theta1, theta2)):
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    y_cond[:, i] = result.y_pred[:, 0]

t = np.arange(nobs)
print("\nt, y, y(t|t-1) [theta1], y(t|t-1) [theta2] -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, y_cond])[:10], 3))
print(np.round(np.column_stack([t, y, y_cond])[-10:], 3))

Time-domain vs. frequency-domain (Whittle) Kalman ML — ects/kalman1c.py

"""Translated from Examples/ects/kalman1c.m -- Harvey [1990],
"Forecasting, Structural Time Series and the Kalman Filter", pages
89-90: compares time-domain (exact) maximum likelihood against
frequency-domain (Whittle) maximum likelihood estimation of the
local-level model's (sigma_epsilon, sigma_eta), on the same Purse.asc
series as kalman1b.py/panel1.py/whittle1.py.

As in the original, the time-domain likelihood is evaluated with a0 = 0,
P0 = 0 during estimation, but the final filter run used to compare the
two estimates' one-step-ahead predictions uses a0 = y[0], P0 = 0 (the
original's own, slightly inconsistent, choice -- preserved here rather
than "fixed")."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter
from quanttoolbox.econometrics.whittle import whittle_local_level

y = np.array(
    [
        10,
        15,
        10,
        10,
        12,
        10,
        7,
        17,
        10,
        14,
        8,
        17,
        14,
        18,
        3,
        9,
        11,
        10,
        6,
        12,
        14,
        10,
        25,
        29,
        33,
        33,
        12,
        19,
        16,
        19,
        19,
        12,
        34,
        15,
        36,
        29,
        26,
        21,
        17,
        19,
        13,
        20,
        24,
        12,
        6,
        14,
        6,
        12,
        9,
        11,
        17,
        12,
        8,
        14,
        14,
        12,
        5,
        8,
        10,
        3,
        16,
        8,
        8,
        7,
        12,
        6,
        10,
        8,
        10,
        5,
        7,
    ],
    dtype=float,
)
nobs = y.shape[0]


def _ssm(sigma_epsilon: float, sigma_eta: float) -> StateSpaceModel:
    return StateSpaceModel(
        z=np.array([[1.0]]),
        d=np.array([0.0]),
        h=np.array([[sigma_epsilon**2]]),
        t=np.array([[1.0]]),
        c=np.array([0.0]),
        r=np.array([[1.0]]),
        q=np.array([[sigma_eta**2]]),
    )


def ll_ml(theta: np.ndarray) -> np.ndarray:
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], np.array([0.0]), np.array([[0.0]]))
    return result.log_l


sv = np.array([3.0, 1.0])

theta1 = ml_estimation(ll_ml, sv).theta
theta2 = whittle_local_level(y, sv).theta

print("theta1 (sigma_epsilon, sigma_eta), time-domain ML:", np.round(theta1, 4))
print("theta2 (sigma_epsilon, sigma_eta), frequency-domain (Whittle) ML:", np.round(theta2, 4))

a0 = np.array([y[0]])
p0 = np.array([[0.0]])

y_cond = np.zeros((nobs, 2))
for i, theta in enumerate((theta1, theta2)):
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    y_cond[:, i] = result.y_pred[:, 0]

t = np.arange(nobs)
print("\nt, y, y(t|t-1) [time-domain ML], y(t|t-1) [Whittle ML] -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, y_cond])[:10], 3))
print(np.round(np.column_stack([t, y, y_cond])[-10:], 3))

Time-varying-beta model from a deliberately bad starting point — stats/kalman2.py

"""Translated from Examples/stats/kalman2.m -- same simulated
time-varying-coefficient regression as kalman1.py, but with the Kalman
filter's initial state a0 fixed at [-20, -20] (deliberately far from the
truth) rather than estimated, to illustrate how quickly the filter
recovers from a poor starting guess. Bounded ML estimation, same
optimizer substitution (`scipy.optimize.minimize` with
`method="L-BFGS-B"`), and the same seeded-RNG substitution as kalman1.py
all apply here too.

Note (preserved faithfully from the original): `sv`/`lb`/`ub` are still
5-dimensional and `theta[0:2]` are still optimized over, exactly as in
kalman1.py -- but since a0 is fixed externally here rather than being
set from `theta[0:2]` inside the likelihood, those two optimizer
dimensions have no effect on the objective. This looks like a leftover
from kalman1.py's script that was never trimmed down; it's kept as-is
rather than "fixed", since the goal is an exact translation."""

import numpy as np
from scipy.optimize import minimize

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

rng = np.random.default_rng(0)
nobs = 200

b0 = 10.0 + np.cumsum(1.5 * rng.standard_normal(nobs))
b1 = 4.0 + np.cumsum(0.2 * rng.standard_normal(nobs))

x0 = np.ones(nobs)
x1 = 25.0 * rng.standard_normal(nobs)
u = 2.0 * rng.standard_normal(nobs)
y = x0 * b0 + x1 * b1 + u

z = np.zeros((1, 2, nobs))
z[0, 0, :] = x0
z[0, 1, :] = x1
d = np.zeros((1, nobs))
t_mat = np.repeat(np.eye(2)[:, :, None], nobs, axis=2)
c = np.zeros((2, nobs))
r_mat = np.repeat(np.eye(2)[:, :, None], nobs, axis=2)

a0_fixed = np.array([-20.0, -20.0])


def _neg_sum_logl(theta: np.ndarray) -> float:
    th = theta**2
    p0 = np.diag(th[3:5])
    h = np.full((1, 1, nobs), th[2])
    q = np.repeat(p0[:, :, None], nobs, axis=2)
    ssm = StateSpaceModel(z=z, d=d, h=h, t=t_mat, c=c, r=r_mat, q=q)
    result = kalman_filter(ssm, y[:, None], a0_fixed, p0)
    return -float(np.sum(result.log_l))


sv = np.ones(5)
lb = np.array([-10.0, -10.0, 0.00001, 0.0001, 0.0001])
ub = np.array([10.0, 10.0, 100.0, 100.0, 100.0])

opt_result = minimize(_neg_sum_logl, sv, method="L-BFGS-B", bounds=list(zip(lb, ub, strict=True)))
theta = opt_result.x**2

print("theta (unused a0 dims, sigma_epsilon^2, sigma_b0^2, sigma_b1^2):", np.round(theta, 4))

p0 = np.diag(theta[3:5])
h = np.full((1, 1, nobs), theta[2])
q = np.repeat(p0[:, :, None], nobs, axis=2)
ssm = StateSpaceModel(z=z, d=d, h=h, t=t_mat, c=c, r=r_mat, q=q)
result = kalman_filter(ssm, y[:, None], a0_fixed, p0)
a_cond = result.a_pred

s = np.arange(1, nobs + 1)
print("\ns, true b0, filtered b0, true b1, filtered b1 -- first/last 10:")
print(np.round(np.column_stack([s, b0, a_cond[:, 0], b1, a_cond[:, 1]])[:10], 4))
print(np.round(np.column_stack([s, b0, a_cond[:, 0], b1, a_cond[:, 1]])[-10:], 4))

Time-varying-beta regression via bounded maximum likelihood — stats/kalman1.py

"""Translated from Examples/stats/kalman1.m -- simulates a 200-
observation time-varying-coefficient regression y_t = b0_t + x1_t*b1_t +
u_t, where b0_t and b1_t are independent random walks starting at 10 and
4 respectively, then recovers (a0, sigma_epsilon, sigma_b0, sigma_b1) by
*bounded* maximum likelihood (the original uses `fmincon` with explicit
lower/upper bounds) and re-runs the time-varying Kalman filter with the
estimated parameters.

`quanttoolbox.econometrics.estimation.ml_estimation` doesn't support box
constraints, so -- to preserve the original's bounded optimization
faithfully rather than silently dropping it -- this calls
`scipy.optimize.minimize(..., method="L-BFGS-B", bounds=...)` directly
on the negative summed Kalman log-likelihood instead. As in the
original, `theta[2:5]` (sigma_epsilon, sigma_b0, sigma_b1) are bounded
to `[1e-5, 100]` (`[1e-4, 100]` for the two sigma_b's) and then *also*
squared inside the likelihood -- `P0`'s diagonal is set from the same
squared theta[3]/theta[4] that also become the (constant) process-noise
covariance Q, i.e. the model assumes the state's initial uncertainty
equals its per-period innovation variance, exactly as the original does.

The original draws from MATLAB's unseeded `randn`; a fixed seed
(`np.random.default_rng(0)`) is substituted here, in the same call order
(b0 innovations, b1 innovations, x1, u) as the original."""

import numpy as np
from scipy.optimize import minimize

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

rng = np.random.default_rng(0)
nobs = 200

b0 = 10.0 + np.cumsum(1.5 * rng.standard_normal(nobs))
b1 = 4.0 + np.cumsum(0.2 * rng.standard_normal(nobs))

x0 = np.ones(nobs)
x1 = 25.0 * rng.standard_normal(nobs)
u = 2.0 * rng.standard_normal(nobs)
y = x0 * b0 + x1 * b1 + u

z = np.zeros((1, 2, nobs))
z[0, 0, :] = x0
z[0, 1, :] = x1
d = np.zeros((1, nobs))
t_mat = np.repeat(np.eye(2)[:, :, None], nobs, axis=2)
c = np.zeros((2, nobs))
r_mat = np.repeat(np.eye(2)[:, :, None], nobs, axis=2)


def _neg_sum_logl(theta: np.ndarray) -> float:
    th = theta.copy()
    th[2:5] = th[2:5] ** 2
    a0 = th[0:2]
    p0 = np.diag(th[3:5])
    h = np.full((1, 1, nobs), th[2])
    q = np.repeat(p0[:, :, None], nobs, axis=2)
    ssm = StateSpaceModel(z=z, d=d, h=h, t=t_mat, c=c, r=r_mat, q=q)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    return -float(np.sum(result.log_l))


sv = np.ones(5)
lb = np.array([-10.0, -10.0, 0.00001, 0.0001, 0.0001])
ub = np.array([10.0, 10.0, 100.0, 100.0, 100.0])

opt_result = minimize(_neg_sum_logl, sv, method="L-BFGS-B", bounds=list(zip(lb, ub, strict=True)))
theta = opt_result.x.copy()
theta[2:5] = theta[2:5] ** 2

print("theta (a0_1, a0_2, sigma_epsilon^2, sigma_b0^2, sigma_b1^2):", np.round(theta, 4))

a0 = theta[0:2]
p0 = np.diag(theta[3:5])
h = np.full((1, 1, nobs), theta[2])
q = np.repeat(p0[:, :, None], nobs, axis=2)
ssm = StateSpaceModel(z=z, d=d, h=h, t=t_mat, c=c, r=r_mat, q=q)
result = kalman_filter(ssm, y[:, None], a0, p0)
a_cond = result.a_pred

s = np.arange(1, nobs + 1)
print("\ns, true b0, filtered b0, true b1, filtered b1 -- first/last 10:")
print(np.round(np.column_stack([s, b0, a_cond[:, 0], b1, a_cond[:, 1]])[:10], 4))
print(np.round(np.column_stack([s, b0, a_cond[:, 0], b1, a_cond[:, 1]])[-10:], 4))

Time-varying-coefficient model with variances estimated by ML — ects/kalman4b.py

"""Translated from Examples/ects/kalman4b.m -- same simulated
time-varying-coefficient regression as kalman4a.py, but instead of using
the true (sigma_epsilon, sigma_beta1, sigma_beta2), estimates them by
maximum likelihood (starting from sv = [1, 1, 1]) and re-runs the
time-varying Kalman filter with the estimated variances, comparing the
recovered beta_t path against the simulated true path.

As in kalman1b.py/kalman2c.py/kalman3d.py, `theta = sqrt(theta.^2)` is
translated as `np.abs(theta)`."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

rng = np.random.default_rng(0)
n_t = 200

sigma1 = 0.5
sigma2 = 0.25
sigma = 1.0

beta1 = np.cumsum(sigma1 * rng.standard_normal(n_t))
beta2 = np.cumsum(sigma2 * rng.standard_normal(n_t))
beta = np.column_stack([beta1, beta2])
x = rng.random((n_t, 2))

y = np.sum(x * beta, axis=1) + sigma * rng.standard_normal(n_t)

a0 = np.zeros(2)
p0 = np.zeros((2, 2))

z = x[None, :, :].transpose(0, 2, 1)  # (1, 2, nT)
d = np.zeros((1, n_t))
t_mat = np.repeat(np.eye(2)[:, :, None], n_t, axis=2)
c = np.zeros((2, n_t))
r = np.repeat(np.eye(2)[:, :, None], n_t, axis=2)


def _ssm(sigma_epsilon: float, sigma_beta1: float, sigma_beta2: float) -> StateSpaceModel:
    h = np.full((1, 1, n_t), sigma_epsilon**2)
    q = np.repeat(np.diag([sigma_beta1**2, sigma_beta2**2])[:, :, None], n_t, axis=2)
    return StateSpaceModel(z=z, d=d, h=h, t=t_mat, c=c, r=r, q=q)


def ssm_logl(theta: np.ndarray) -> np.ndarray:
    sigma_epsilon, sigma_beta1, sigma_beta2 = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_beta1, sigma_beta2)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    return result.log_l


sv = np.ones(3)
ml_result = ml_estimation(ssm_logl, sv)
theta_hat = np.abs(ml_result.theta)

sigma_epsilon_hat, sigma1_hat, sigma2_hat = theta_hat
print(
    "Estimated (sigma_epsilon, sigma_beta1, sigma_beta2):",
    np.round(theta_hat, 4),
    "-- true:",
    [sigma, sigma1, sigma2],
)

ssm_hat = _ssm(sigma_epsilon_hat, sigma1_hat, sigma2_hat)
result = kalman_filter(ssm_hat, y[:, None], a0, p0)
at = result.a_filt

t = np.arange(1, n_t + 1)
print("\nt, true beta1, filtered beta1, true beta2, filtered beta2 -- first/last 10:")
print(np.round(np.column_stack([t, beta[:, 0], at[:, 0], beta[:, 1], at[:, 1]])[:10], 4))
print(np.round(np.column_stack([t, beta[:, 0], at[:, 0], beta[:, 1], at[:, 1]])[-10:], 4))

Time-varying-coefficient regression recovered by a Kalman filter — ects/kalman4a.py

"""Translated from Examples/ects/kalman4a.m -- simulates a
200-observation time-varying-coefficient regression y_t = x_t1*beta1_t +
x_t2*beta2_t + eps_t, where beta1_t and beta2_t follow independent
random walks (`recserar(sigma*randn(nT,1), 0, 1.0)` with AR coefficient
1.0 is exactly a cumulative sum starting at 0, translated as
`np.cumsum`), then recovers the beta_t path with a time-varying-Z Kalman
filter (Z_t = x_t, a random-walk state-space model with T=I) and
compares the filtered estimate against the simulated true path.

The original draws from MATLAB's unseeded `randn`/`rand`; a fixed seed
(`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter

rng = np.random.default_rng(0)
n_t = 200

sigma1 = 0.5
sigma2 = 0.25
sigma = 1.0

beta1 = np.cumsum(sigma1 * rng.standard_normal(n_t))
beta2 = np.cumsum(sigma2 * rng.standard_normal(n_t))
beta = np.column_stack([beta1, beta2])
x = rng.random((n_t, 2))

y = np.sum(x * beta, axis=1) + sigma * rng.standard_normal(n_t)

z = x[None, :, :].transpose(0, 2, 1)  # (1, 2, nT)
d = np.zeros((1, n_t))
h = np.full((1, 1, n_t), sigma**2)
t_mat = np.repeat(np.eye(2)[:, :, None], n_t, axis=2)
c = np.zeros((2, n_t))
r = np.repeat(np.eye(2)[:, :, None], n_t, axis=2)
q = np.repeat(np.diag([sigma1**2, sigma2**2])[:, :, None], n_t, axis=2)

ssm = StateSpaceModel(z=z, d=d, h=h, t=t_mat, c=c, r=r, q=q)

a0 = np.zeros(2)
p0 = np.zeros((2, 2))
result = kalman_filter(ssm, y[:, None], a0, p0)

at = result.a_filt

t = np.arange(1, n_t + 1)
print("t, true beta1, filtered beta1, true beta2, filtered beta2 -- first/last 10:")
print(np.round(np.column_stack([t, beta[:, 0], at[:, 0], beta[:, 1], at[:, 1]])[:10], 4))
print(np.round(np.column_stack([t, beta[:, 0], at[:, 0], beta[:, 1], at[:, 1]])[-10:], 4))

print("\nRMSE beta1:", round(float(np.sqrt(np.mean((beta[:, 0] - at[:, 0]) ** 2))), 4))
print("RMSE beta2:", round(float(np.sqrt(np.mean((beta[:, 1] - at[:, 1]) ** 2))), 4))

econometrics.whittle

Python alternatives

Keep — Whittle (frequency-domain) estimation isn't implemented in statsmodels, arch, or other common packages. Genuinely fills a gap.

quanttoolbox.econometrics.whittle

Whittle (frequency-domain) maximum likelihood estimation.

Ported from QuantToolBox/ects/{whittle_estimation, whittle_constrained_estimation,whittle_local_level, whittle_local_linear_trend}.m and QuantToolBox/maths/{periodogram,pdgm}.m (identical duplicates, both covered by periodogram here rather than in maths.py).

Translation notes:

• Whittle estimation fits a parametric spectral density function sdf(lambda, theta) to the sample periodogram by maximizing the Whittle approximate log-likelihood, via the same linearly-restricted optimization pattern as econometrics.estimation (theta = RR @ gamma
• r).
• MATLAB's "method of scoring" custom Newton-Raphson-with-information- matrix optimizer (used only when analytical gradients are supplied) is not ported as a separate code path; scipy.optimize.minimize (BFGS, using the analytical gradient if supplied) is used throughout instead -- both converge to the same MLE for this smooth, well-behaved objective, and BFGS is the standard choice.
• whittle_estimation.m's nargin-dispatch wrapper is not ported separately; call whittle_estimation directly with restriction=None for the unconstrained case.

periodogram(y, scaling=False)

Raw periodogram of y via FFT.

scaling=False (default): I(lambda) = |FFT(y)|^2 / n. scaling=True: additionally divided by 2*pi (spectral-density scaling).

Original: maths/{periodogram,pdgm}.m

Returns (lambda, fft_coeffs, periodogram_values).

Source code in src/quanttoolbox/econometrics/whittle.py

def periodogram(y: np.ndarray, scaling: bool = False) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Raw periodogram of y via FFT.

    scaling=False (default): I(lambda) = |FFT(y)|^2 / n.
    scaling=True: additionally divided by 2*pi (spectral-density scaling).

    Original: maths/{periodogram,pdgm}.m

    Returns (lambda, fft_coeffs, periodogram_values).
    """
    y = np.asarray(y, dtype=float).flatten()
    n = y.shape[0]
    fft_coeffs = np.fft.fft(y)
    intensity = np.abs(fft_coeffs) ** 2 / n
    lam = 2 * np.pi * np.arange(n) / n
    if scaling:
        intensity = intensity / (2 * np.pi)
    return lam, fft_coeffs, intensity

whittle_estimation(y, sdf_fn, sv, restriction=None, weights=1.0, sdf_jacobian_fn=None, cov='hessian', config=None)

Whittle (frequency-domain) MLE: fit a parametric spectral density sdf_fn(lambda, theta) to the periodogram of y by maximizing the Whittle log-likelihood sum(-log(2pi) - 0.5log(f) - 0.5*I/f), with theta = RR @ gamma + r.

Original: ects/{whittle_estimation,whittle_constrained_estimation}.m

Source code in src/quanttoolbox/econometrics/whittle.py

def whittle_estimation(
    y: np.ndarray,
    sdf_fn: Callable[[np.ndarray, np.ndarray], np.ndarray],
    sv: np.ndarray,
    restriction: tuple[np.ndarray, np.ndarray] | None = None,
    weights: np.ndarray | float = 1.0,
    sdf_jacobian_fn: Callable[[np.ndarray, np.ndarray], np.ndarray] | None = None,
    cov: str = "hessian",
    config: EstimationConfig | None = None,
) -> WhittleEstimationResult:
    """Whittle (frequency-domain) MLE: fit a parametric spectral density
    sdf_fn(lambda, theta) to the periodogram of y by maximizing the
    Whittle log-likelihood sum(-log(2*pi) - 0.5*log(f) - 0.5*I/f), with
    theta = RR @ gamma + r.

    Original: ects/{whittle_estimation,whittle_constrained_estimation}.m
    """
    if config is None:
        config = EstimationConfig()

    sv = np.asarray(sv, dtype=float).flatten()
    n_params = sv.shape[0]
    rr, r = (np.eye(n_params), np.zeros(n_params)) if restriction is None else restriction
    rr = np.asarray(rr, dtype=float)
    r = np.asarray(r, dtype=float).flatten()
    n_free = rr.shape[1]

    lam, _, intensity = periodogram(y, scaling=True)
    n_obs = intensity.shape[0]
    w = (
        np.broadcast_to(weights, (n_obs,))
        if np.isscalar(weights)
        else np.asarray(weights, dtype=float)
    )

    def _logpdf(theta: np.ndarray) -> np.ndarray:
        f_sdf = sdf_fn(lam, theta)
        return -np.log(2 * np.pi) - 0.5 * np.log(f_sdf) - 0.5 * (intensity / f_sdf)

    def _neg_sum_logl(gamma: np.ndarray) -> float:
        theta = rr @ gamma + r
        ll = w * _logpdf(theta)
        return float(-np.sum(ll[~np.isnan(ll)]))

    jac_for_opt = None
    if sdf_jacobian_fn is not None:

        def jac_for_opt(gamma: np.ndarray) -> np.ndarray:
            theta = rr @ gamma + r
            f_sdf = sdf_fn(lam, theta)
            f_sdf = np.where(f_sdf == 0, 0.01, f_sdf)
            j = sdf_jacobian_fn(lam, theta)
            j = w[:, None] * j
            inv_f = 1.0 / f_sdf
            grad_terms = 0.5 * (intensity * inv_f - 1.0)[:, None] * (j * inv_f[:, None])
            grad_terms = w[:, None] * grad_terms
            return -rr.T @ np.sum(grad_terms[~np.isnan(grad_terms).any(axis=1)], axis=0)

    result = minimize(_neg_sum_logl, sv, method="BFGS", jac=jac_for_opt)
    gamma = result.x
    theta = rr @ gamma + r

    logl = _logpdf(theta)
    valid = ~np.isnan(logl)
    n_valid = int(np.sum(valid))
    df = n_valid - n_free
    sum_logl = float(np.sum(logl[valid]))

    if cov in ("opg", "hc"):
        if sdf_jacobian_fn is not None:
            f_sdf = sdf_fn(lam, theta)
            f_sdf = np.where(f_sdf == 0, 0.01, f_sdf)
            j = sdf_jacobian_fn(lam, theta)
            inv_f = 1.0 / f_sdf
            g = 0.5 * (intensity * inv_f - 1.0)[:, None] * (j * inv_f[:, None])
        else:
            g = _numerical_jacobian(_logpdf, theta)
        g = w[:, None] * g
        g_valid = g[valid] @ rr
        gg_valid = g_valid.T @ g_valid

    if cov == "opg":
        try:
            vcv_free = np.linalg.inv(gg_valid)
        except np.linalg.LinAlgError:
            vcv_free = np.full((n_free, n_free), np.nan)
        cov_type = "opg"
    else:
        h = -_numerical_hessian(_neg_sum_logl, gamma)
        h_valid = h  # already in the gamma (restricted) parameterization
        try:
            inv_h_valid = np.linalg.inv(h_valid)
        except np.linalg.LinAlgError:
            inv_h_valid = np.full((n_free, n_free), np.nan)

        if cov == "hc":
            vcv_free = inv_h_valid @ gg_valid @ inv_h_valid
            cov_type = "hc"
        else:
            vcv_free = -inv_h_valid
            cov_type = "hessian"

    vcv = rr @ vcv_free @ rr.T
    stderr = np.sqrt(np.diag(vcv))

    return WhittleEstimationResult(
        theta=theta,
        stderr=stderr,
        vcv=vcv,
        log_l=logl,
        sum_log_l=sum_logl,
        df=df,
        n_obs=n_obs,
        n_obs_valid=n_valid,
        cov_type=cov_type,
        converged=result.success,
    )

whittle_local_level(y, sv, **kwargs)

Fit a local-level (random-walk-plus-noise) model y_t = mu_t + eps_t, mu_t = mu_{t-1} + eta_t via Whittle estimation of (sigma_epsilon, sigma_eta) on the first-differenced series.

Original: ects/whittle_local_level.m

Source code in src/quanttoolbox/econometrics/whittle.py

def whittle_local_level(y: np.ndarray, sv: np.ndarray, **kwargs) -> WhittleEstimationResult:
    """Fit a local-level (random-walk-plus-noise) model
    y_t = mu_t + eps_t, mu_t = mu_{t-1} + eta_t via Whittle estimation of
    (sigma_epsilon, sigma_eta) on the first-differenced series.

    Original: ects/whittle_local_level.m
    """
    sv = np.asarray(sv, dtype=float).flatten()
    if sv.shape[0] != 2:
        raise ValueError("whittle_local_level: sv must have length 2")

    y = np.asarray(y, dtype=float).flatten()
    dy = y[1:] - y[:-1]
    dy = dy[~np.isnan(dy)]

    return whittle_estimation(
        dy, _local_level_sdf, sv, sdf_jacobian_fn=_local_level_sdf_jacobian, **kwargs
    )

whittle_local_linear_trend(y, sv, **kwargs)

Fit a local-linear-trend model via Whittle estimation of (sigma_epsilon, sigma_eta, sigma_zeta) on the twice-differenced series.

Original: ects/whittle_local_linear_trend.m

Source code in src/quanttoolbox/econometrics/whittle.py

def whittle_local_linear_trend(y: np.ndarray, sv: np.ndarray, **kwargs) -> WhittleEstimationResult:
    """Fit a local-linear-trend model via Whittle estimation of
    (sigma_epsilon, sigma_eta, sigma_zeta) on the twice-differenced series.

    Original: ects/whittle_local_linear_trend.m
    """
    sv = np.asarray(sv, dtype=float).flatten()
    if sv.shape[0] != 3:
        raise ValueError("whittle_local_linear_trend: sv must have length 3")

    y = np.asarray(y, dtype=float).flatten()
    dy = y[1:] - y[:-1]
    dy2 = dy[1:] - dy[:-1]
    dy2 = dy2[~np.isnan(dy2)]

    return whittle_estimation(
        dy2, _local_linear_trend_sdf, sv, sdf_jacobian_fn=_local_linear_trend_sdf_jacobian, **kwargs
    )

Examples

Local Linear Trend model fit by Whittle (frequency-domain) ML — ects/kalman2b.py

"""Translated from Examples/ects/kalman2b.m -- Harvey [1990], pages
89-90: frequency-domain (Whittle) maximum-likelihood estimation of the
Local Linear Trend model's (sigma_epsilon, sigma_eta, sigma_zeta) on the
same 61-observation Gnp.asc series as kalman2a.py, then runs the Kalman
filter with the estimated parameters (numeric core only, plot dropped).

The original calls `whittle_local_linear_trend` twice -- once with
`WHITTLE_algorithm = 1` (BFGS) and once with `WHITTLE_algorithm = 3`
("scoring method") -- keeping only the second result. As established in
`whittle1.py`/`whittle.py`'s module docstring, the port only has a single
optimizer path (BFGS), so there is just one call here."""

import numpy as np

from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter
from quanttoolbox.econometrics.whittle import whittle_local_linear_trend

y = np.array(
    [
        116.8,
        120.0,
        123.2,
        130.2,
        131.4,
        125.6,
        124.5,
        134.3,
        135.2,
        151.8,
        146.4,
        139.0,
        127.8,
        147.0,
        165.9,
        165.5,
        179.4,
        190.0,
        189.8,
        190.9,
        203.6,
        183.5,
        169.3,
        144.2,
        141.5,
        154.3,
        169.5,
        193.0,
        203.2,
        192.9,
        209.4,
        227.2,
        263.7,
        297.8,
        337.1,
        361.3,
        355.2,
        312.6,
        309.9,
        323.7,
        324.1,
        255.3,
        383.4,
        395.1,
        412.8,
        406.0,
        438.0,
        446.1,
        452.5,
        447.3,
        475.9,
        487.7,
        497.2,
        529.8,
        551.0,
        581.1,
        617.8,
        658.1,
        675.2,
        706.6,
        724.7,
    ],
    dtype=float,
)
nobs = y.shape[0]

# BFGS method (very sensitive to the starting values -- kept as-is)
sv = 15.0 * np.ones(3)
result_est = whittle_local_linear_trend(y, sv)
theta = result_est.theta

sigma_epsilon, sigma_eta, sigma_zeta = theta
print("theta (sigma_epsilon, sigma_eta, sigma_zeta):", np.round(theta, 4))

ssm = StateSpaceModel(
    z=np.array([[1.0, 0.0]]),
    d=np.array([0.0]),
    h=np.array([[sigma_epsilon**2]]),
    t=np.array([[1.0, 1.0], [0.0, 1.0]]),
    c=np.array([0.0, 0.0]),
    r=np.eye(2),
    q=np.diag([sigma_eta**2, sigma_zeta**2]),
)

a0 = np.array([y[0], 0.0])
p0 = np.zeros((2, 2))
result = kalman_filter(ssm, y[:, None], a0, p0)

t = np.arange(1909, 1909 + nobs)
print("\nt, y, level a(t|t-1), slope a(t|t-1) -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, result.a_pred])[:10], 3))
print(np.round(np.column_stack([t, y, result.a_pred])[-10:], 3))

Raw periodogram of a 4-observation series — maths/pdgm1.py

"""Translated from Examples/maths/pdgm1.m -- raw periodogram of a
4-observation series, unscaled and scaled by 1/(2*pi). Despite the
tracker's earlier note, this is a small numeric example, not a plot;
`pdgm(y)`/`pdgm(y,1)` map directly to
`periodogram(y, scaling=False)`/`periodogram(y, scaling=True)`."""

import numpy as np

from quanttoolbox.econometrics.whittle import periodogram

y = np.array([0.1, 0.3, 0.4, 0.5])

lam, f, i1 = periodogram(y, scaling=False)
_, _, i2 = periodogram(y, scaling=True)

print("FFT coefficients:")
print(f)

print("\nlambda, unscaled periodogram, scaled periodogram:")
print(np.round(np.column_stack([lam, i1, i2]), 10))

Raw periodogram of an 8-observation series — maths/pdgm2.py

"""Translated from Examples/maths/pdgm2.m -- same as pdgm1.py, on an
8-observation series."""

import numpy as np

from quanttoolbox.econometrics.whittle import periodogram

y = np.array([0.01, -0.3, 1.4, 1.5, -0.7, 0.3, 0.1, 0.0])

lam, f, i1 = periodogram(y, scaling=False)
_, _, i2 = periodogram(y, scaling=True)

print("FFT coefficients:")
print(f)

print("\nlambda, unscaled periodogram, scaled periodogram:")
print(np.round(np.column_stack([lam, i1, i2]), 10))

Time-domain vs. frequency-domain (Whittle) Kalman ML — ects/kalman1c.py

"""Translated from Examples/ects/kalman1c.m -- Harvey [1990],
"Forecasting, Structural Time Series and the Kalman Filter", pages
89-90: compares time-domain (exact) maximum likelihood against
frequency-domain (Whittle) maximum likelihood estimation of the
local-level model's (sigma_epsilon, sigma_eta), on the same Purse.asc
series as kalman1b.py/panel1.py/whittle1.py.

As in the original, the time-domain likelihood is evaluated with a0 = 0,
P0 = 0 during estimation, but the final filter run used to compare the
two estimates' one-step-ahead predictions uses a0 = y[0], P0 = 0 (the
original's own, slightly inconsistent, choice -- preserved here rather
than "fixed")."""

import numpy as np

from quanttoolbox.econometrics.estimation import ml_estimation
from quanttoolbox.econometrics.kalman import StateSpaceModel, kalman_filter
from quanttoolbox.econometrics.whittle import whittle_local_level

y = np.array(
    [
        10,
        15,
        10,
        10,
        12,
        10,
        7,
        17,
        10,
        14,
        8,
        17,
        14,
        18,
        3,
        9,
        11,
        10,
        6,
        12,
        14,
        10,
        25,
        29,
        33,
        33,
        12,
        19,
        16,
        19,
        19,
        12,
        34,
        15,
        36,
        29,
        26,
        21,
        17,
        19,
        13,
        20,
        24,
        12,
        6,
        14,
        6,
        12,
        9,
        11,
        17,
        12,
        8,
        14,
        14,
        12,
        5,
        8,
        10,
        3,
        16,
        8,
        8,
        7,
        12,
        6,
        10,
        8,
        10,
        5,
        7,
    ],
    dtype=float,
)
nobs = y.shape[0]


def _ssm(sigma_epsilon: float, sigma_eta: float) -> StateSpaceModel:
    return StateSpaceModel(
        z=np.array([[1.0]]),
        d=np.array([0.0]),
        h=np.array([[sigma_epsilon**2]]),
        t=np.array([[1.0]]),
        c=np.array([0.0]),
        r=np.array([[1.0]]),
        q=np.array([[sigma_eta**2]]),
    )


def ll_ml(theta: np.ndarray) -> np.ndarray:
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], np.array([0.0]), np.array([[0.0]]))
    return result.log_l


sv = np.array([3.0, 1.0])

theta1 = ml_estimation(ll_ml, sv).theta
theta2 = whittle_local_level(y, sv).theta

print("theta1 (sigma_epsilon, sigma_eta), time-domain ML:", np.round(theta1, 4))
print("theta2 (sigma_epsilon, sigma_eta), frequency-domain (Whittle) ML:", np.round(theta2, 4))

a0 = np.array([y[0]])
p0 = np.array([[0.0]])

y_cond = np.zeros((nobs, 2))
for i, theta in enumerate((theta1, theta2)):
    sigma_epsilon, sigma_eta = np.abs(theta)
    ssm = _ssm(sigma_epsilon, sigma_eta)
    result = kalman_filter(ssm, y[:, None], a0, p0)
    y_cond[:, i] = result.y_pred[:, 0]

t = np.arange(nobs)
print("\nt, y, y(t|t-1) [time-domain ML], y(t|t-1) [Whittle ML] -- first/last 10 observations:")
print(np.round(np.column_stack([t, y, y_cond])[:10], 3))
print(np.round(np.column_stack([t, y, y_cond])[-10:], 3))

Whittle estimation with a custom Bloomfield spectral density — ects/whittle2.py

"""Translated from Examples/ects/whittle2.m -- fits a Bloomfield
exponential spectral density (Dzhaparidze [1986], p.125) of order 4 to
the 500-observation Whittle2.asc series (embedded directly below) via
Whittle (frequency-domain) maximum likelihood, using a custom `sdf_fn`
(not one of the two built-in models in `whittle.py`).

The original loops `WHITTLE_algorithm` over 1 and 3 to compare two
optimizer implementations of the same problem; as established in
whittle1.py/whittle.py's module docstring, the port only has a single
optimizer path (BFGS), so there is just one call here. The trailing `0,
0` positional arguments in `whittle_estimation(y, sdf, sv, 0, 0)`
(restriction, weights) are the GAUSS/MATLAB-toolbox "use the default"
sentinel and map to `restriction=None, weights=1.0`, i.e. the keyword
defaults -- so they're simply omitted below."""

import io

import numpy as np

from quanttoolbox.econometrics.whittle import whittle_estimation

_WHITTLE2_DATA = """
0.0000000000 0.0000000000 -0.0070836130 -0.5896696335 -0.8779192198 0.7739608676 1.8687461178 1.1272122892 0.8647355090 0.5353193167
1.4279839232 0.9644046845 1.3189917804 2.7187174736 -0.0265101591 0.6346647533 2.3763572557 1.6758165418 -0.2518899869 0.5140942167
0.1579653667 -0.4346519759 -0.1831780798 1.3166446987 0.3045805948 0.0834965505 -0.2721840893 0.0579250026 1.1119641044 1.4873721188
0.4545474039 0.1292263709 -0.2767236072 0.0256743859 0.1264964685 -0.1985822034 -0.8440399349 -1.2100443680 0.7077841221 -0.3373624656
-0.1321149645 -0.1310297612 1.4732762921 0.0401150657 1.7418622017 0.9637928035 1.4688291220 1.1283169081 0.1486218473 1.2245699016
0.3862002629 -0.1335736872 0.2834417006 -1.1079203204 -2.7073527498 -3.0159361194 -2.5576362307 -3.8092608480 -4.9122551627 -3.8686181852
-3.0861011760 -3.4524470855 -2.6502211283 -2.3989063355 -1.3435186309 -1.6064055128 -0.6907728960 -1.4407067218 -0.5739565416 -0.3463894808
-2.1216361003 -0.2147356769 -1.9100128865 -1.9747299730 -1.5910681402 -0.9187868533 -1.0914199362 -0.3385159067 -0.3647006079 1.5247997776
-0.5228135946 0.5056004160 -0.0257481548 0.2163354188 -0.2721861552 0.9405338024 1.7673421054 0.1413767482 1.6243141484 1.5587770346
2.6652868896 1.6286323647 2.8103954285 2.6524479848 2.4308940477 3.3312573035 3.3088050444 2.6008098961 1.8369500109 1.1410512106
3.0394022862 1.0523694042 2.5200472651 2.9892246724 2.0830570503 1.4955517799 2.3183629160 1.0644185033 2.7340511449 2.2515023281
3.6342729393 2.5802329592 2.7820409892 3.7279442817 2.2499036746 2.0065945505 0.8871471692 0.0330941250 0.5321053663 0.3813835930
-0.5512652790 -2.6518014663 0.6433895026 -0.7323520561 0.6973701077 -0.4084725830 2.2674168094 0.1393896738 2.3343004187 1.4625949048
0.9883464917 2.1361150179 1.1711251253 3.2753310787 0.1242448706 0.9387692773 0.5470237839 0.5273237274 -0.9745772654 0.0528308833
0.0346781921 0.4968722679 0.0796808836 0.4007851085 0.0097992552 -0.2884873429 0.2790453437 -0.0907409944 -0.4533929754 0.2809886552
0.1186961215 0.1057917070 -0.3088850630 -0.1852426195 -1.0324080925 -0.9631908464 0.0322378066 0.3797769137 1.1541253083 1.6998789937
1.4799438100 1.6883127074 0.9142055586 0.0446462834 -1.1610371583 -1.0074875444 -2.4734098135 -2.0240106582 -2.9179975860 -1.9191929931
-1.3765893671 -0.9672529077 -1.2966709606 -1.3173713811 0.7213099535 0.4906078485 -0.8528896393 -0.4324850144 -0.9786037142 0.3769233947
-0.4911307534 0.6246007893 0.6185136580 0.6790791050 -0.3252387268 1.5383808954 -0.5467010145 1.7203477314 -0.1061824069 1.0262551761
0.0708747464 0.7857255645 0.2075316074 0.3800127639 -0.0905192787 -0.8019745157 0.4156637503 0.5232173272 0.8132963581 0.1564847546
0.4342887815 1.7303198792 -0.2217797885 1.0789013218 0.3088358751 0.8889960553 0.1617417644 1.0331249799 0.4216053018 1.5648264901
1.8950547595 0.4578086794 1.6257064986 -0.2535497016 0.5416104085 2.6463790740 1.6893045547 1.9270592835 0.7287558798 0.6684989182
0.1679263462 -1.5959773398 -0.7516155784 0.5111695249 -0.3550029160 -0.3311187318 -1.8190305221 0.0930786526 1.1119247681 -0.8257627260
0.2016678520 -0.7612220943 0.8287686500 1.2063078277 -1.0166464244 0.5974603776 -0.7116803099 0.3341812880 -0.6067663371 -0.2591284789
-0.1597049420 -0.4497015391 -0.1399188820 0.3458410864 0.9474377476 0.7049700830 -0.4267943271 -1.0293282759 -0.1540129988 -1.1610944356
-0.5929523078 -0.7720332068 0.5530785848 -0.5867182079 0.8326287985 1.4451946900 1.1147540292 0.9830127334 0.0378238303 -0.6291839605
1.3424331757 0.9508200788 0.5311259848 1.7388682086 0.0715751871 0.7693526238 -0.1655581265 0.4975225875 -0.6663411229 -0.8515400234
-1.3308249865 -1.0485229241 -2.1957151355 -2.3735630056 -0.4587174437 -1.9898753152 -1.0193999556 -2.0615469770 -0.1011848265 -0.6923864072
-1.1902095441 -2.9435014484 -2.1426363241 -3.4844420966 -1.9907697329 -2.4524783217 -2.9946879812 -2.9685397226 -2.5867271246 -1.9996182278
0.0244487531 -0.7270891728 -0.6294386039 -0.9354022809 -0.5119473786 0.2353552607 -1.0381172842 0.0628670857 0.4052962105 -0.5110590221
0.9331818471 0.4968508214 1.0926383046 1.5331303247 0.3754681107 1.1784390454 0.8723818579 0.9889494801 -0.1411149874 -0.2923141929
1.7037124688 0.2971921413 2.1966852059 0.8175186668 1.1316179201 2.3725381352 1.7352804463 2.5459059704 0.8234880658 3.0121675922
2.7627588508 2.5240524966 2.4953758871 3.3476557205 1.6576147597 1.9314050375 1.1701034585 2.0773832380 1.4258449122 2.7006228848
-0.3224222055 0.7212858203 0.1805975389 -0.2921989262 0.0983673767 -0.7666920251 0.5301336616 0.4418883598 -0.1635880775 0.4437044336
0.7979835741 1.3242366416 2.0673973326 0.6387922782 2.3739658417 1.5645394816 0.2410157887 1.7106585843 0.2852274604 1.3293673648
-0.2360870111 0.2407690379 0.4182344625 0.8797050449 0.2687394013 1.5394221528 1.3449795102 1.1004288321 2.6567681580 3.0806073943
2.2835489957 2.7606598800 2.6205482816 2.0343066797 1.9808169524 1.6092639738 2.7527261064 2.3249644962 1.0105636472 2.0950495986
2.4426386353 1.6264072552 2.4974127440 -0.5119215761 0.5550011209 0.9635268173 -0.0826570943 0.1524618788 1.0529150285 -0.9660969180
-0.0110718871 0.3647659977 1.1772600940 -0.0157121637 -0.0209914320 2.9715377863 0.9752705668 1.3084549557 0.4720939864 0.3169124084
0.4499938310 1.2296193584 0.8792811568 1.5439568105 1.4828102421 0.8063835330 -0.0479898305 0.3385269847 1.5251557894 1.3606765429
-0.0648313149 -0.3907450756 -0.2351488463 0.3062231362 0.3772359222 -0.1462953655 -0.8979647086 0.8713691913 0.8265067724 0.2212269213
-1.2757641961 -1.4425053632 -0.3919020752 -0.4369329863 0.7000134311 1.7649491501 1.4967053798 0.5704227840 -0.9505335532 -1.0818904396
0.4234553970 0.1732907885 -0.4597249981 0.4382408439 1.2410708686 -0.0142690244 1.6680404889 1.8453462011 1.7015967207 1.5851986813
1.8115357863 -0.6370623397 2.1489073971 1.0415094477 1.6965608202 2.0181015453 1.9682054045 0.6645872984 1.3063893453 -0.4619972029
-0.0540595420 1.8735316299 2.3981134490 2.5395124143 2.8285664350 3.1442444357 2.5092920047 1.5098956896 0.6273532482 0.7780697737
1.4471674499 1.4609563336 2.1686345059 1.5631902956 3.0100024713 2.5086544116 1.9912642517 1.3741078654 3.3142570417 3.2495004289
2.1885111886 2.7473103630 1.4101328702 1.6804557008 -0.1848695253 1.8179271485 -0.2239986918 0.5895546129 0.8048923486 1.2632640533
0.1253163171 -1.1458201207 0.2729754981 -0.5112755016 0.3566590967 0.3234413144 -1.2866935997 0.2444798346 -0.0941590879 1.7620134253
1.1184802676 -0.4885682635 0.3333726488 -0.1491835030 -0.6291479889 -0.4930987348 -0.1209361200 1.5936256050 0.8742496016 0.9537200520
0.7066716928 0.7432255084 1.5859461053 1.1064332961 1.0494326399 1.0598341144 0.9134280771 1.8421271552 2.7617449377 2.8015870733
"""

y = np.loadtxt(io.StringIO(_WHITTLE2_DATA)).flatten()

order = 4


def bloomfield_sdf(lam: np.ndarray, theta: np.ndarray) -> np.ndarray:
    gamma_ = theta[:order]
    sigma = theta[order]
    j = np.arange(1, order + 1)
    w = np.cos(np.outer(lam, j))
    sdf = (sigma**2) * np.exp(2 * w @ gamma_)
    return sdf / (2 * np.pi)


sv = np.ones(order + 1)

result = whittle_estimation(y, bloomfield_sdf, sv)

print("theta (gamma_1..gamma_4, sigma):", np.round(result.theta, 4))
print("stderr:", np.round(result.stderr, 4))
print("sum log-likelihood:", round(result.sum_log_l, 4))

Whittle local-level MLE on real macroeconomic data — ects/whittle1.py

"""Translated from Examples/ects/whittle1.m -- Harvey [1990],
"Forecasting, Structural Time Series and the Kalman Filter", pages
89-90: estimate a local-level model's variance parameters in the
frequency domain via Whittle maximum likelihood, on the same 71-
observation "Purse" series used in panel1.py.

The original loops `WHITTLE_algorithm` over 1, 2, and 3 to compare three
different optimizer implementations of the same Whittle MLE problem;
`whittle_local_level` (see `whittle.py`'s module docstring) ports only a
single optimizer path (`scipy.optimize.minimize` with BFGS), since all
three MATLAB variants solve the identical optimization problem and
should converge to the same estimate -- so there is only one call here,
not three."""

import numpy as np

from quanttoolbox.econometrics.whittle import whittle_local_level

y = np.array(
    [
        10,
        15,
        10,
        10,
        12,
        10,
        7,
        17,
        10,
        14,
        8,
        17,
        14,
        18,
        3,
        9,
        11,
        10,
        6,
        12,
        14,
        10,
        25,
        29,
        33,
        33,
        12,
        19,
        16,
        19,
        19,
        12,
        34,
        15,
        36,
        29,
        26,
        21,
        17,
        19,
        13,
        20,
        24,
        12,
        6,
        14,
        6,
        12,
        9,
        11,
        17,
        12,
        8,
        14,
        14,
        12,
        5,
        8,
        10,
        3,
        16,
        8,
        8,
        7,
        12,
        6,
        10,
        8,
        10,
        5,
        7,
    ],
    dtype=float,
)

sv = np.array([3.0, 1.0])

result = whittle_local_level(y, sv)

print("theta (sigma_epsilon, sigma_eta):", np.round(result.theta, 4))
print("stderr:", np.round(result.stderr, 4))
print("sum log-likelihood:", round(result.sum_log_l, 4))

econometrics.tests (ADF)

Python alternatives

Already switched to statsmodels.tsa.stattools.adfuller. The arch package's arch.unitroot module has a wider family of unit-root tests (Phillips-Perron, DFGLS, KPSS, Zivot-Andrews) if more than ADF is ever needed.

quanttoolbox.econometrics.tests

Unit-root testing (Augmented Dickey-Fuller).

Ported from QuantToolBox/ects/adf_test.m

Translation notes:

• The original hand-rolls the ADF regression and hardcodes ~15x50 tables of MacKinnon critical values (interpolated for sample size and significance level) across three model specifications (no constant, a constant, and a constant + trend) -- several hundred numeric literals. statsmodels.tsa.stattools.adfuller implements the same test using the same underlying MacKinnon (1994/2010) critical-value approximation, is the standard, actively-maintained Python implementation, and is already a project dependency -- so this wraps it directly instead of transcribing the original's tables.
• The original scans lag orders 0..pLags and all three specifications ("n"/"c"/"ct") in one call, returning a (3, pLags+1) result grid; adf_test here reproduces that scan structure, calling adfuller(..., autolag=None, maxlag=p) for each (specification, lag) combination to match the original's fixed-lag (not automatically selected) behavior.
• wald_test (originally also in ects/) lives in econometrics.estimation instead, alongside the estimators it's most often used with.

adf_test(y, max_lags)

Augmented Dickey-Fuller unit-root test, scanned over lag orders 0..max_lags and all three specifications (no constant, constant, constant + trend).

Original: ects/adf_test.m

Uses statsmodels.tsa.stattools.adfuller -- see module docstring.

Source code in src/quanttoolbox/econometrics/tests.py

def adf_test(y: np.ndarray, max_lags: int) -> ADFTestResult:
    """Augmented Dickey-Fuller unit-root test, scanned over lag orders
    0..max_lags and all three specifications (no constant, constant,
    constant + trend).

    Original: ects/adf_test.m

    Uses statsmodels.tsa.stattools.adfuller -- see module docstring.
    """
    y = np.asarray(y, dtype=float).flatten()
    y = y[~np.isnan(y)]

    lags = np.arange(0, max_lags + 1)
    n_specs = len(_SPECIFICATIONS)
    n_lags = lags.shape[0]

    tau = np.full((n_specs, n_lags), np.nan)
    p_value = np.full((n_specs, n_lags), np.nan)
    critical_values = np.full((n_specs, n_lags, 3), np.nan)

    for i, spec in enumerate(_SPECIFICATIONS):
        for j, lag in enumerate(lags):
            try:
                stat, pval, _, _, crit = adfuller(
                    y, maxlag=int(lag), regression=spec, autolag=None, store=False, regresults=False
                )
            except (ValueError, np.linalg.LinAlgError):
                continue
            tau[i, j] = stat
            p_value[i, j] = pval
            critical_values[i, j] = [crit["1%"], crit["5%"], crit["10%"]]

    return ADFTestResult(
        specification=np.array(_SPECIFICATIONS),
        lags=lags,
        tau=tau,
        p_value=p_value,
        critical_values=critical_values,
    )