mfe.tests_stat¶

mfe.tests_stat ¶

mfe.tests_stat — Statistical tests for financial time series.

Serial correlation ljung_box Ljung-Box Q statistic (not robust to heteroskedasticity) lm_test HAC-robust LM serial correlation test (MFE lmtest1)

Conditional heteroskedasticity arch_lm Engle (1982) ARCH-LM test

Forecast evaluation mincer_zarnowitz MZ regression-based forecast evaluation diebold_mariano DM test for equal predictive accuracy (MSE/MAE/QLIKE)

MZResult `dataclass` ¶

MZResult(alpha: float, beta: float, alpha_se: float, beta_se: float, t_stat_alpha: float, t_stat_beta: float, f_stat: float, f_pvalue: float, r_squared: float, n_obs: int)

Mincer-Zarnowitz regression output.

ljung_box ¶

ljung_box(data: FloatArray, max_lags: int = 10) -> LjungBoxResult

Ljung-Box Q statistic for serial correlation.

Q_k = T(T+2) * sum_{j=1}^{k} rho_hat_j^2 / (T - j)

Under H0 of no autocorrelation, Q_k ~ chi2(k) asymptotically.

NOTE: Not appropriate for heteroskedastic data (use lm_test instead).

Parameters:

Name	Type	Description	Default
`data`	`FloatArray`		required
`max_lags`	`number of lags to test; returns one statistic per lag 1..max_lags`		`10`

Source code in src/mfe/tests_stat/serial.py

def ljung_box(
    data: FloatArray,
    max_lags: int = 10,
) -> LjungBoxResult:
    """
    Ljung-Box Q statistic for serial correlation.

    Q_k = T(T+2) * sum_{j=1}^{k} rho_hat_j^2 / (T - j)

    Under H0 of no autocorrelation, Q_k ~ chi2(k) asymptotically.

    NOTE: Not appropriate for heteroskedastic data (use lm_test instead).

    Parameters
    ----------
    data     : (T,) time series (demeaned or residuals)
    max_lags : number of lags to test; returns one statistic per lag 1..max_lags
    """
    x = np.asarray(data, dtype=np.float64)
    x = x - x.mean()
    T = len(x)

    # Sample autocorrelations
    acov0 = float(x @ x) / T
    rho = np.array([
        float(x[lag:] @ x[:T - lag]) / (T * acov0)
        for lag in range(1, max_lags + 1)
    ])

    lags_arr = np.arange(1, max_lags + 1)
    # Q_k = cumulative sum up to lag k
    q_terms = rho ** 2 / (T - lags_arr)
    Q = T * (T + 2) * np.cumsum(q_terms)

    p_vals = np.array([
        float(1 - stats.chi2.cdf(Q[k], df=k + 1))
        for k in range(max_lags)
    ])

    return LjungBoxResult(statistics=Q, p_values=p_vals, lags=lags_arr.astype(float))

lm_test ¶

lm_test(data: FloatArray, max_lags: int = 10, robust: bool = True) -> LMTestResult

LM test for serial correlation in up to max_lags lags.

The test is an LM-test for testing the null that all of the regression coefficients are zero in the auxiliary regression of y_t on lags 1..Q. The null tested is H0: phi_1 = phi_2 = ... = phi_Q = 0.

Parameters:

Name	Type	Description	Default
`data`	`FloatArray`		required
`max_lags`	`maximum lag order to test`		`10`
`robust`	`bool`	`if False, use classical homoskedastic VCV`	`True`

Notes

Equivalent to MFE toolbox lmtest1.m.

Source code in src/mfe/tests_stat/serial.py

def lm_test(
    data: FloatArray,
    max_lags: int = 10,
    robust: bool = True,
) -> LMTestResult:
    """
    LM test for serial correlation in up to max_lags lags.

    The test is an LM-test for testing the null that all of the regression
    coefficients are zero in the auxiliary regression of y_t on lags 1..Q.
    The null tested is H0: phi_1 = phi_2 = ... = phi_Q = 0.

    Parameters
    ----------
    data     : (T,) time series (typically GARCH residuals or raw returns)
    max_lags : maximum lag order to test
    robust   : if True (default), use heteroskedasticity-robust (White) VCV
               if False, use classical homoskedastic VCV

    Notes
    -----
    Equivalent to MFE toolbox lmtest1.m.
    """
    x = np.asarray(data, dtype=np.float64)
    T = len(x)
    x_dm = x - x.mean()

    stats_arr = np.empty(max_lags, dtype=np.float64)
    pvals_arr = np.empty(max_lags, dtype=np.float64)

    for q in range(1, max_lags + 1):
        # Build regressor matrix: T-q by q lags of x_dm
        n = T - q
        X = np.column_stack([x_dm[q - j - 1: T - j - 1] for j in range(q)])  # (n, q)
        eps_tilde = x_dm[q:]  # residual under null = demeaned data

        # Score: s_t = eps_tilde_t * X_t
        scores = X * eps_tilde[:, None]   # (n, q)
        s_bar = scores.mean(axis=0)       # (q,)

        if robust:
            # White sandwich: S = T^{-1} sum scores'*scores
            S = scores.T @ scores / n    # (q, q)
        else:
            # Homoskedastic: S = sigma^2 * T^{-1} X'X
            sigma2 = float(np.mean(eps_tilde ** 2))
            S = sigma2 * (X.T @ X) / n

        try:
            S_inv = np.linalg.inv(S)
        except np.linalg.LinAlgError:
            stats_arr[q - 1] = np.nan
            pvals_arr[q - 1] = np.nan
            continue

        lm = float(n * s_bar @ S_inv @ s_bar)
        stats_arr[q - 1] = lm
        pvals_arr[q - 1] = float(1 - stats.chi2.cdf(lm, df=q))

    return LMTestResult(
        statistics=stats_arr,
        p_values=pvals_arr,
        lags=np.arange(1, max_lags + 1, dtype=float),
        robust=robust,
    )

mincer_zarnowitz ¶

mincer_zarnowitz(realized: FloatArray, forecast: FloatArray, nw_lags: int = 0) -> MZResult

Mincer-Zarnowitz regression: realized = alpha + beta * forecast + eps.

Tests alpha = 0, beta = 1 (unbiased forecast), and the joint H0.

Parameters:

Name	Type	Default
`realized`	`(T,) actual realized values (e.g. RV_t)`	required
`forecast`	`(T,) model forecasts (e.g. h_{t\|t-1})`	required
`nw_lags`	`int`	`0`

Returns:

Type	Description
`MZResult`

Source code in src/mfe/tests_stat/forecast_eval.py

def mincer_zarnowitz(
    realized: FloatArray,
    forecast: FloatArray,
    nw_lags: int = 0,
) -> MZResult:
    """
    Mincer-Zarnowitz regression: realized = alpha + beta * forecast + eps.

    Tests alpha = 0, beta = 1 (unbiased forecast), and the joint H0.

    Parameters
    ----------
    realized : (T,) actual realized values (e.g. RV_t)
    forecast : (T,) model forecasts (e.g. h_{t|t-1})
    nw_lags  : Newey-West lags for HAC standard errors

    Returns
    -------
    MZResult
    """
    y = np.asarray(realized, dtype=np.float64)
    yhat = np.asarray(forecast, dtype=np.float64)
    T = len(y)

    X = np.column_stack([np.ones(T), yhat])  # (T, 2)
    XTX_inv = np.linalg.inv(X.T @ X)
    beta_hat = XTX_inv @ (X.T @ y)
    resid = y - X @ beta_hat

    # R-squared
    ss_tot = np.sum((y - np.mean(y)) ** 2)
    r2 = 1 - np.sum(resid ** 2) / ss_tot

    # Standard errors
    if nw_lags > 0:
        scores = X * resid[:, None]
        B_nw = newey_west(scores, bandwidth=nw_lags)
        vcv = XTX_inv @ B_nw @ XTX_inv
    else:
        s2 = np.sum(resid ** 2) / (T - 2)
        vcv = s2 * XTX_inv

    se = np.sqrt(np.diag(vcv))
    alpha, beta = float(beta_hat[0]), float(beta_hat[1])
    alpha_se, beta_se = float(se[0]), float(se[1])

    t_alpha = alpha / alpha_se
    t_beta = (beta - 1.0) / beta_se

    # Joint F-test: alpha=0, beta=1
    R = np.array([[1.0, 0.0], [0.0, 1.0]])
    r_vec = np.array([0.0, 1.0])
    diff = R @ beta_hat - r_vec
    try:
        f_stat = float(diff @ np.linalg.solve(R @ vcv @ R.T, diff)) / 2
        f_pvalue = float(1 - stats.f.cdf(f_stat, dfn=2, dfd=T - 2))
    except np.linalg.LinAlgError:
        f_stat = np.nan
        f_pvalue = np.nan

    return MZResult(
        alpha=alpha,
        beta=beta,
        alpha_se=alpha_se,
        beta_se=beta_se,
        t_stat_alpha=t_alpha,
        t_stat_beta=t_beta,
        f_stat=f_stat,
        f_pvalue=f_pvalue,
        r_squared=float(r2),
        n_obs=T,
    )

diebold_mariano ¶

diebold_mariano(errors1: FloatArray, errors2: FloatArray, loss: str = 'mse', nw_lags: int | None = None, alternative: str = 'two-sided') -> DMResult

Diebold-Mariano test for equal predictive accuracy.

Diebold, F.X. & Mariano, R.S. (1995): "Comparing Predictive Accuracy", JBES.

Parameters:

Name	Type	Description	Default
`errors1`	`(T,) forecast error arrays from two models`		required
`errors2`	`(T,) forecast error arrays from two models`		required
`loss`	`str`		`'mse'`
`nw_lags`	`int \| None`		`None`
`alternative`	`str`	`"greater" means model1 is less accurate (H1: model1 worse)`	`'two-sided'`

Source code in src/mfe/tests_stat/forecast_eval.py

def diebold_mariano(
    errors1: FloatArray,
    errors2: FloatArray,
    loss: str = "mse",
    nw_lags: int | None = None,
    alternative: str = "two-sided",
) -> DMResult:
    """
    Diebold-Mariano test for equal predictive accuracy.

    Diebold, F.X. & Mariano, R.S. (1995): "Comparing Predictive Accuracy",
    JBES.

    Parameters
    ----------
    errors1, errors2 : (T,) forecast error arrays from two models
    loss             : "mse" | "mae" | "qlike"
    nw_lags          : Newey-West lags; if None uses int(T^{1/3})
    alternative      : "two-sided" | "greater" | "less"
                       "greater" means model1 is less accurate (H1: model1 worse)
    """
    e1 = np.asarray(errors1, dtype=np.float64)
    e2 = np.asarray(errors2, dtype=np.float64)
    T = len(e1)

    if loss == "mse":
        d = e1 ** 2 - e2 ** 2
    elif loss == "mae":
        d = np.abs(e1) - np.abs(e2)
    elif loss == "qlike":
        # Patton (2011) QLIKE: L = sigma2/h - log(sigma2/h) - 1
        # For errors only (assuming h=1 or normalized): approximated as e^2 - log(e^2)
        d = e1 ** 2 - np.log(e1 ** 2 + 1e-30) - (e2 ** 2 - np.log(e2 ** 2 + 1e-30))
    else:
        raise ValueError(f"loss must be 'mse', 'mae', or 'qlike', got '{loss}'")

    d_bar = float(np.mean(d))

    if nw_lags is None:
        nw_lags = int(T ** (1 / 3))

    # HAC variance of d_bar
    d_centered = (d - d_bar)[:, None]
    B_nw = newey_west(d_centered, bandwidth=nw_lags)
    var_d = float(B_nw[0, 0]) / T
    se_d = np.sqrt(max(var_d, 0.0))

    dm_stat = d_bar / max(se_d, 1e-30)

    if alternative == "two-sided":
        p_val = 2 * (1 - stats.norm.cdf(abs(dm_stat)))
    elif alternative == "greater":
        p_val = 1 - stats.norm.cdf(dm_stat)
    else:
        p_val = stats.norm.cdf(dm_stat)

    return DMResult(
        statistic=float(dm_stat),
        p_value=float(p_val),
        loss_diff_mean=d_bar,
        n_obs=T,
    )

arch_lm ¶

ARCH-LM test for conditional heteroskedasticity.

Engle, R.F. (1982): "Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation", Econometrica.

The test regresses squared residuals on their lagged values: eps_t^2 = alpha_0 + alpha_1 * eps_{t-1}^2 + ... + alpha_q * eps_{t-q}^2 + u_t

H0: alpha_1 = ... = alpha_q = 0 (no ARCH effects)

Test statistic: T * R^2 ~ chi2(q) under H0.

Also computes the F-form (more reliable in small samples).

arch_lm ¶

arch_lm(residuals: FloatArray, lags: int = 5) -> ARCHLMResult

Engle ARCH-LM test for conditional heteroskedasticity.

Parameters:

Name	Type	Description	Default
`residuals`	`(T,) return or residual series`		required
`lags`	`int`		`5`

Returns:

Type	Description
`ARCHLMResult`	.lm_stat / .lm_pval : LM test (chi2 form, q df) .f_stat / .f_pval : F-test form (more reliable in small T)

Source code in src/mfe/tests_stat/arch_lm.py

def arch_lm(
    residuals: FloatArray,
    lags: int = 5,
) -> ARCHLMResult:
    """
    Engle ARCH-LM test for conditional heteroskedasticity.

    Parameters
    ----------
    residuals : (T,) return or residual series
    lags      : number of lags q in the auxiliary regression

    Returns
    -------
    ARCHLMResult
        .lm_stat / .lm_pval : LM test (chi2 form, q df)
        .f_stat  / .f_pval  : F-test form (more reliable in small T)
    """
    e = np.asarray(residuals, dtype=np.float64)
    T = len(e)
    e2 = e ** 2

    # Build auxiliary regression: e2_t on constant + e2_{t-1..t-q}
    n = T - lags
    y = e2[lags:]                            # (n,)
    X = np.column_stack(                     # (n, lags+1)
        [np.ones(n)] + [e2[lags - j - 1: T - j - 1] for j in range(lags)]
    )

    # OLS
    try:
        beta = np.linalg.lstsq(X, y, rcond=None)[0]
    except np.linalg.LinAlgError:
        return ARCHLMResult(np.nan, np.nan, np.nan, np.nan, lags, np.nan)

    y_hat = X @ beta
    resid = y - y_hat
    ss_res = float(resid @ resid)
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0

    # LM stat = T * R^2
    lm = float(n * r2)
    lm_pval = float(1 - stats.chi2.cdf(lm, df=lags))

    # F stat
    k = lags       # number of restrictions
    denom = ss_res / (n - lags - 1) if n > lags + 1 else np.nan
    f_stat = ((ss_tot - ss_res) / k) / denom if np.isfinite(denom) and denom > 0 else np.nan
    f_pval = float(1 - stats.f.cdf(f_stat, dfn=k, dfd=n - k - 1)) if np.isfinite(f_stat) else np.nan

    return ARCHLMResult(
        lm_stat=lm,
        lm_pval=lm_pval,
        f_stat=f_stat,
        f_pval=f_pval,
        lags=lags,
        r_squared=r2,
    )

forecast_eval ¶

Forecast evaluation tests for volatility models.

Mincer-Zarnowitz (1969): regression-based evaluation of forecasts. Diebold & Mariano (1995): test for equal predictive accuracy. Hansen (2005): Superior Predictive Ability test.

MZResult `dataclass` ¶

MZResult(alpha: float, beta: float, alpha_se: float, beta_se: float, t_stat_alpha: float, t_stat_beta: float, f_stat: float, f_pvalue: float, r_squared: float, n_obs: int)

Mincer-Zarnowitz regression output.

mincer_zarnowitz ¶

mincer_zarnowitz(realized: FloatArray, forecast: FloatArray, nw_lags: int = 0) -> MZResult

Mincer-Zarnowitz regression: realized = alpha + beta * forecast + eps.

Tests alpha = 0, beta = 1 (unbiased forecast), and the joint H0.

Parameters:

Name	Type	Default
`realized`	`(T,) actual realized values (e.g. RV_t)`	required
`forecast`	`(T,) model forecasts (e.g. h_{t\|t-1})`	required
`nw_lags`	`int`	`0`

Returns:

Type	Description
`MZResult`

Source code in src/mfe/tests_stat/forecast_eval.py

def mincer_zarnowitz(
    realized: FloatArray,
    forecast: FloatArray,
    nw_lags: int = 0,
) -> MZResult:
    """
    Mincer-Zarnowitz regression: realized = alpha + beta * forecast + eps.

    Tests alpha = 0, beta = 1 (unbiased forecast), and the joint H0.

    Parameters
    ----------
    realized : (T,) actual realized values (e.g. RV_t)
    forecast : (T,) model forecasts (e.g. h_{t|t-1})
    nw_lags  : Newey-West lags for HAC standard errors

    Returns
    -------
    MZResult
    """
    y = np.asarray(realized, dtype=np.float64)
    yhat = np.asarray(forecast, dtype=np.float64)
    T = len(y)

    X = np.column_stack([np.ones(T), yhat])  # (T, 2)
    XTX_inv = np.linalg.inv(X.T @ X)
    beta_hat = XTX_inv @ (X.T @ y)
    resid = y - X @ beta_hat

    # R-squared
    ss_tot = np.sum((y - np.mean(y)) ** 2)
    r2 = 1 - np.sum(resid ** 2) / ss_tot

    # Standard errors
    if nw_lags > 0:
        scores = X * resid[:, None]
        B_nw = newey_west(scores, bandwidth=nw_lags)
        vcv = XTX_inv @ B_nw @ XTX_inv
    else:
        s2 = np.sum(resid ** 2) / (T - 2)
        vcv = s2 * XTX_inv

    se = np.sqrt(np.diag(vcv))
    alpha, beta = float(beta_hat[0]), float(beta_hat[1])
    alpha_se, beta_se = float(se[0]), float(se[1])

    t_alpha = alpha / alpha_se
    t_beta = (beta - 1.0) / beta_se

    # Joint F-test: alpha=0, beta=1
    R = np.array([[1.0, 0.0], [0.0, 1.0]])
    r_vec = np.array([0.0, 1.0])
    diff = R @ beta_hat - r_vec
    try:
        f_stat = float(diff @ np.linalg.solve(R @ vcv @ R.T, diff)) / 2
        f_pvalue = float(1 - stats.f.cdf(f_stat, dfn=2, dfd=T - 2))
    except np.linalg.LinAlgError:
        f_stat = np.nan
        f_pvalue = np.nan

    return MZResult(
        alpha=alpha,
        beta=beta,
        alpha_se=alpha_se,
        beta_se=beta_se,
        t_stat_alpha=t_alpha,
        t_stat_beta=t_beta,
        f_stat=f_stat,
        f_pvalue=f_pvalue,
        r_squared=float(r2),
        n_obs=T,
    )

diebold_mariano ¶

diebold_mariano(errors1: FloatArray, errors2: FloatArray, loss: str = 'mse', nw_lags: int | None = None, alternative: str = 'two-sided') -> DMResult

Diebold-Mariano test for equal predictive accuracy.

Diebold, F.X. & Mariano, R.S. (1995): "Comparing Predictive Accuracy", JBES.

Parameters:

Name	Type	Description	Default
`errors1`	`(T,) forecast error arrays from two models`		required
`errors2`	`(T,) forecast error arrays from two models`		required
`loss`	`str`		`'mse'`
`nw_lags`	`int \| None`		`None`
`alternative`	`str`	`"greater" means model1 is less accurate (H1: model1 worse)`	`'two-sided'`

Source code in src/mfe/tests_stat/forecast_eval.py

def diebold_mariano(
    errors1: FloatArray,
    errors2: FloatArray,
    loss: str = "mse",
    nw_lags: int | None = None,
    alternative: str = "two-sided",
) -> DMResult:
    """
    Diebold-Mariano test for equal predictive accuracy.

    Diebold, F.X. & Mariano, R.S. (1995): "Comparing Predictive Accuracy",
    JBES.

    Parameters
    ----------
    errors1, errors2 : (T,) forecast error arrays from two models
    loss             : "mse" | "mae" | "qlike"
    nw_lags          : Newey-West lags; if None uses int(T^{1/3})
    alternative      : "two-sided" | "greater" | "less"
                       "greater" means model1 is less accurate (H1: model1 worse)
    """
    e1 = np.asarray(errors1, dtype=np.float64)
    e2 = np.asarray(errors2, dtype=np.float64)
    T = len(e1)

    if loss == "mse":
        d = e1 ** 2 - e2 ** 2
    elif loss == "mae":
        d = np.abs(e1) - np.abs(e2)
    elif loss == "qlike":
        # Patton (2011) QLIKE: L = sigma2/h - log(sigma2/h) - 1
        # For errors only (assuming h=1 or normalized): approximated as e^2 - log(e^2)
        d = e1 ** 2 - np.log(e1 ** 2 + 1e-30) - (e2 ** 2 - np.log(e2 ** 2 + 1e-30))
    else:
        raise ValueError(f"loss must be 'mse', 'mae', or 'qlike', got '{loss}'")

    d_bar = float(np.mean(d))

    if nw_lags is None:
        nw_lags = int(T ** (1 / 3))

    # HAC variance of d_bar
    d_centered = (d - d_bar)[:, None]
    B_nw = newey_west(d_centered, bandwidth=nw_lags)
    var_d = float(B_nw[0, 0]) / T
    se_d = np.sqrt(max(var_d, 0.0))

    dm_stat = d_bar / max(se_d, 1e-30)

    if alternative == "two-sided":
        p_val = 2 * (1 - stats.norm.cdf(abs(dm_stat)))
    elif alternative == "greater":
        p_val = 1 - stats.norm.cdf(dm_stat)
    else:
        p_val = stats.norm.cdf(dm_stat)

    return DMResult(
        statistic=float(dm_stat),
        p_value=float(p_val),
        loss_diff_mean=d_bar,
        n_obs=T,
    )

serial ¶

Serial correlation tests for financial time series.

Ljung & Box (1978): Q-statistic for autocorrelation up to lag K. Godfrey (1978) / Breusch (1978): LM test for serial correlation.

Key difference from statsmodels: statsmodels.stats.diagnostic.acorr_ljungbox only provides the standard LB test. lmtest here provides a HAC-robust LM variant (lmtest1 from MFE toolbox) which is appropriate for heteroskedastic series — the standard LB test is not.

The heteroskedasticity-robust LM test is essentially an LR-class test: LM = T * s_hat' * S_hat^{-1} * s_hat where s_hat = T^{-1} X'eps_tilde (gradient under null) and S_hat is estimated under the alternative using the White sandwich.

ljung_box ¶

ljung_box(data: FloatArray, max_lags: int = 10) -> LjungBoxResult

Ljung-Box Q statistic for serial correlation.

Q_k = T(T+2) * sum_{j=1}^{k} rho_hat_j^2 / (T - j)

Under H0 of no autocorrelation, Q_k ~ chi2(k) asymptotically.

NOTE: Not appropriate for heteroskedastic data (use lm_test instead).

Parameters:

Name	Type	Description	Default
`data`	`FloatArray`		required
`max_lags`	`number of lags to test; returns one statistic per lag 1..max_lags`		`10`

Source code in src/mfe/tests_stat/serial.py

def ljung_box(
    data: FloatArray,
    max_lags: int = 10,
) -> LjungBoxResult:
    """
    Ljung-Box Q statistic for serial correlation.

    Q_k = T(T+2) * sum_{j=1}^{k} rho_hat_j^2 / (T - j)

    Under H0 of no autocorrelation, Q_k ~ chi2(k) asymptotically.

    NOTE: Not appropriate for heteroskedastic data (use lm_test instead).

    Parameters
    ----------
    data     : (T,) time series (demeaned or residuals)
    max_lags : number of lags to test; returns one statistic per lag 1..max_lags
    """
    x = np.asarray(data, dtype=np.float64)
    x = x - x.mean()
    T = len(x)

    # Sample autocorrelations
    acov0 = float(x @ x) / T
    rho = np.array([
        float(x[lag:] @ x[:T - lag]) / (T * acov0)
        for lag in range(1, max_lags + 1)
    ])

    lags_arr = np.arange(1, max_lags + 1)
    # Q_k = cumulative sum up to lag k
    q_terms = rho ** 2 / (T - lags_arr)
    Q = T * (T + 2) * np.cumsum(q_terms)

    p_vals = np.array([
        float(1 - stats.chi2.cdf(Q[k], df=k + 1))
        for k in range(max_lags)
    ])

    return LjungBoxResult(statistics=Q, p_values=p_vals, lags=lags_arr.astype(float))

lm_test ¶

lm_test(data: FloatArray, max_lags: int = 10, robust: bool = True) -> LMTestResult

LM test for serial correlation in up to max_lags lags.

The test is an LM-test for testing the null that all of the regression coefficients are zero in the auxiliary regression of y_t on lags 1..Q. The null tested is H0: phi_1 = phi_2 = ... = phi_Q = 0.

Parameters:

Name	Type	Description	Default
`data`	`FloatArray`		required
`max_lags`	`maximum lag order to test`		`10`
`robust`	`bool`	`if False, use classical homoskedastic VCV`	`True`

Notes

Equivalent to MFE toolbox lmtest1.m.

Source code in src/mfe/tests_stat/serial.py

def lm_test(
    data: FloatArray,
    max_lags: int = 10,
    robust: bool = True,
) -> LMTestResult:
    """
    LM test for serial correlation in up to max_lags lags.

    The test is an LM-test for testing the null that all of the regression
    coefficients are zero in the auxiliary regression of y_t on lags 1..Q.
    The null tested is H0: phi_1 = phi_2 = ... = phi_Q = 0.

    Parameters
    ----------
    data     : (T,) time series (typically GARCH residuals or raw returns)
    max_lags : maximum lag order to test
    robust   : if True (default), use heteroskedasticity-robust (White) VCV
               if False, use classical homoskedastic VCV

    Notes
    -----
    Equivalent to MFE toolbox lmtest1.m.
    """
    x = np.asarray(data, dtype=np.float64)
    T = len(x)
    x_dm = x - x.mean()

    stats_arr = np.empty(max_lags, dtype=np.float64)
    pvals_arr = np.empty(max_lags, dtype=np.float64)

    for q in range(1, max_lags + 1):
        # Build regressor matrix: T-q by q lags of x_dm
        n = T - q
        X = np.column_stack([x_dm[q - j - 1: T - j - 1] for j in range(q)])  # (n, q)
        eps_tilde = x_dm[q:]  # residual under null = demeaned data

        # Score: s_t = eps_tilde_t * X_t
        scores = X * eps_tilde[:, None]   # (n, q)
        s_bar = scores.mean(axis=0)       # (q,)

        if robust:
            # White sandwich: S = T^{-1} sum scores'*scores
            S = scores.T @ scores / n    # (q, q)
        else:
            # Homoskedastic: S = sigma^2 * T^{-1} X'X
            sigma2 = float(np.mean(eps_tilde ** 2))
            S = sigma2 * (X.T @ X) / n

        try:
            S_inv = np.linalg.inv(S)
        except np.linalg.LinAlgError:
            stats_arr[q - 1] = np.nan
            pvals_arr[q - 1] = np.nan
            continue

        lm = float(n * s_bar @ S_inv @ s_bar)
        stats_arr[q - 1] = lm
        pvals_arr[q - 1] = float(1 - stats.chi2.cdf(lm, df=q))

    return LMTestResult(
        statistics=stats_arr,
        p_values=pvals_arr,
        lags=np.arange(1, max_lags + 1, dtype=float),
        robust=robust,
    )

mfe.tests_stat¶

mfe.tests_stat ¶

MZResult dataclass ¶

ljung_box ¶

lm_test ¶

mincer_zarnowitz ¶

diebold_mariano ¶

arch_lm ¶

arch_lm ¶

forecast_eval ¶

MZResult dataclass ¶

mincer_zarnowitz ¶

diebold_mariano ¶

serial ¶

ljung_box ¶

lm_test ¶

MZResult `dataclass` ¶

MZResult `dataclass` ¶