diff --git a/arch/tests/univariate/test_mean.py b/arch/tests/univariate/test_mean.py
index 4ffc6ca1e5..2e66914be6 100644
--- a/arch/tests/univariate/test_mean.py
+++ b/arch/tests/univariate/test_mean.py
@@ -1249,6 +1249,76 @@ def test_arch_lm(simulated_data):
     assert "H0: Standardized" in wald.__repr__()
 
 
+def test_li_mak(simulated_data):
+    # --- default call ---
+    zm = ZeroMean(simulated_data, volatility=GARCH())
+    res = zm.fit(disp=DISPLAY)
+    wald = res.li_mak_test()
+
+    nobs = simulated_data.shape[0]
+    default_lags = int(np.ceil(12.0 * np.power(nobs / 100.0, 1 / 4.0)))
+    # df = lags - fitdf where fitdf = GARCH(1,1).num_params - 1 = 2
+    assert wald.df == default_lags - 2
+    assert "uncorrelated" in wald.null
+    assert "Li-Mak" in wald.__repr__()
+    assert wald.stat >= 0.0
+    assert 0.0 <= wald.pval <= 1.0
+    assert_almost_equal(wald.pval, 1 - stats.chi2(wald.df).cdf(wald.stat))
+
+    # --- explicit lags ---
+    wald5 = res.li_mak_test(lags=5)
+    assert wald5.df == 5 - 2  # fitdf=2 for GARCH(1,1)
+    assert wald5.stat >= 0.0
+
+    # --- fitdf=0 gives uncorrected Ljung-Box on squared standardized residuals ---
+    wald_lb = res.li_mak_test(lags=5, fitdf=0)
+    assert wald_lb.df == 5
+
+    # manually compute Ljung-Box on squared standardized residuals
+    e = res.resid / res.conditional_volatility
+    a = e**2
+    a = a[~np.isnan(a)]
+    n = a.shape[0]
+    a_dm = a - a.mean()
+    a_var = (a_dm**2).sum()
+    stat_manual = 0.0
+    for k in range(1, 6):
+        rk = (a_dm[k:] * a_dm[:-k]).sum() / a_var
+        stat_manual += rk**2 / (n - k)
+    stat_manual *= n * (n + 2)
+    assert_almost_equal(wald_lb.stat, stat_manual, decimal=8)
+
+    # --- li_mak != arch_lm on standardized residuals (different formula + df) ---
+    wald_archlm = res.arch_lm_test(lags=5, standardized=True)
+    # ARCH-LM uses nobs*R^2 (OLS); Li-Mak uses Ljung-Box — different statistics
+    assert wald_archlm.stat != pytest.approx(wald5.stat)
+
+    # --- fitdf parameter override ---
+    wald_custom = res.li_mak_test(lags=10, fitdf=5)
+    assert wald_custom.df == 5
+
+
+def test_li_mak_constant_variance(simulated_data):
+    """With constant variance, fitdf defaults to 0 (no df correction needed)."""
+    from arch.univariate.volatility import ConstantVariance
+
+    zm = ZeroMean(simulated_data, volatility=ConstantVariance())
+    res = zm.fit(disp=DISPLAY)
+    wald = res.li_mak_test(lags=5)
+    assert wald.df == 5  # fitdf = ConstantVariance.num_params - 1 = 0
+
+
+def test_li_mak_pval_structure(simulated_data):
+    """p-value, critical values and repr are consistent."""
+    zm = ZeroMean(simulated_data, volatility=GARCH())
+    res = zm.fit(disp=DISPLAY)
+    wald = res.li_mak_test(lags=10)
+    assert len(wald.critical_values) == 3
+    assert "10%" in wald.critical_values
+    assert "Li-Mak" in repr(wald)
+    assert_almost_equal(wald.pval, 1 - stats.chi2(wald.df).cdf(wald.stat))
+
+
 def test_autoscale():
     rs = np.random.RandomState(34254321)
     dist = Normal(seed=rs)
diff --git a/arch/univariate/base.py b/arch/univariate/base.py
index b8543e8460..09089f0a4f 100644
--- a/arch/univariate/base.py
+++ b/arch/univariate/base.py
@@ -1748,6 +1748,95 @@ def arch_lm_test(
             stat, df=lags, null=null, alternative=alt, name="ARCH-LM Test"
         )
 
+    def li_mak_test(
+        self, lags: int | None = None, fitdf: int | None = None
+    ) -> WaldTestStatistic:
+        r"""
+        Li-Mak portmanteau test for remaining ARCH effects in standardized residuals
+
+        Applies the Ljung-Box statistic to the *squared* standardized residuals
+        :math:`\hat{e}_t^2 = (r_t / \hat{\sigma}_t)^2` and corrects the degrees
+        of freedom for the number of estimated GARCH variance parameters, as
+        derived in Li and Mak (1994).  Under a correctly specified GARCH model,
+        the statistic is asymptotically :math:`\chi^2(m - \text{fitdf})`, where
+        :math:`m` is the number of lags and *fitdf* is the number of estimated
+        volatility lag parameters.
+
+        Unlike the ARCH-LM test, which only tests raw residuals for heteroskedasticity,
+        the Li-Mak test directly checks whether the fitted volatility model has
+        removed *all* serial dependence in the conditional variance.
+
+        Parameters
+        ----------
+        lags : int, optional
+            Number of autocorrelation lags to include.  Defaults to
+            :math:`\lceil 1.2\,T^{1/4} \rceil`.
+        fitdf : int, optional
+            Number of estimated volatility lag parameters to subtract from the
+            degrees of freedom.  Defaults to ``volatility.num_params - 1``
+            (all variance parameters excluding the constant :math:`\omega`).
+            Pass ``fitdf=0`` to obtain the uncorrected Ljung-Box statistic.
+
+        Returns
+        -------
+        result : WaldTestStatistic
+            Test result with statistic, p-value, and null/alternative descriptions.
+
+        References
+        ----------
+        .. [1] Li, W. K. and Mak, T. K. (1994). On the Squared Residual
+               Autocorrelations in Non-Linear Time Series with Conditional
+               Heteroskedasticity. *Journal of Time Series Analysis*,
+               15(6), 627--636.
+
+        Examples
+        --------
+        >>> from arch import arch_model
+        >>> import numpy as np
+        >>> rng = np.random.default_rng(0)
+        >>> returns = rng.standard_normal(500)
+        >>> res = arch_model(returns, vol="Garch", p=1, q=1).fit(disp="off")
+        >>> li_mak = res.li_mak_test(lags=10)
+        >>> li_mak.stat  # doctest: +SKIP
+        """
+        resids = self.resid / np.asarray(self.conditional_volatility)
+        resids = resids[~np.isnan(resids)]
+        nobs = resids.shape[0]
+        if nobs < 3:
+            raise ValueError("Test requires at least 3 non-nan observations")
+
+        a = resids**2
+        lags = (
+            int(np.ceil(12.0 * np.power(nobs / 100.0, 1 / 4.0)))
+            if lags is None
+            else lags
+        )
+        lags = max(min(nobs // 2 - 1, lags), 1)
+        assert isinstance(lags, int)
+
+        if fitdf is None:
+            fitdf = max(self.model.volatility.num_params - 1, 0)
+
+        # Ljung-Box Q statistic on squared standardized residuals
+        a_mean = a.mean()
+        a_dm = a - a_mean
+        a_var = (a_dm**2).sum()
+
+        stat = 0.0
+        for k in range(1, lags + 1):
+            rk = (a_dm[k:] * a_dm[:-k]).sum() / a_var
+            stat += rk**2 / (nobs - k)
+        stat *= nobs * (nobs + 2)
+
+        df = max(lags - fitdf, 1)
+        return WaldTestStatistic(
+            stat,
+            df=df,
+            null="Squared standardized residuals are uncorrelated.",
+            alternative="Squared standardized residuals are autocorrelated (ARCH effects remain).",
+            name="Li-Mak Test",
+        )
+
 
 class ARCHModelResult(ARCHModelFixedResult):
     """