Repository: statsmodels/statsmodels (Python)

Issue:
linear_reset can't be set with "HAC" and maxlags, KeyErrors
#### Describe the bug
Hello, thanks for your time. We've tried several things but it still doesn't work even with the newest version of Statsmodels. We think this might be a bug. Issue is: we were trying to use "statsmodels.stats.diagnostic.linear_reset" with HAC as coviance type and maxlags as the key word args, but we got KeyErrors.

#### Code Sample, a copy-pastable example if possible
```
# Import modules
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.diagnostic import linear_reset

# Generate some random data
np.random.seed(42)
x = np.random.randn(100)
y = 2 * x + np.random.randn(100)

# Fit a linear regression model
model = sm.OLS(y, sm.add_constant(x))
# results = model.fit()

# Perform linear reset test
## opt 1
reset_test = linear_reset(results, power=2, test_type='fitted', cov_type='HAC', cov_kwargs={'maxlags': 2})

## opt 2
# test = linear_reset(results, cov_type='HAC', cov_kwds={'maxlags': 2})

## opt 3
# kwargs = {'maxlags': 2}
# test = linear_reset(results, cov_type='HAC', **kwargs)

## opt 4
# test = linear_reset(results, cov_type='HAC', maxlags=2)

```
```python
# Your code here that produces the bug
# This example should be self-contained, and so not rely on external data.
# It should run in a fresh ipython session, and so include all relevant imports.
```
<details>
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File <command-2970740754718372>:19
     16 model = sm.OLS(y, x).fit(cov_type='HAC', cov_kwds={'maxlags': 2})
     18 # Perform linear reset test
---> 19 reset_test = linear_reset(model, power=2, test_type='fitted', cov_type='HAC', cov_kwargs={'maxlags': 2})
     21 # # Perform the linear reset test with HAC covariance and maxlags=2
     22 # test = linear_reset(results, cov_type='HAC', cov_kwds={'maxlags': 2})
     23 
   (...)
     33 
     34 # Print the test results
     35 print(reset_test.summary())

File /databricks/python/lib/python3.9/site-packages/pandas/util/_decorators.py:207, in deprecate_kwarg.<locals>._deprecate_kwarg.<locals>.wrapper(*args, **kwargs)
    205     else:
    206         kwargs[new_arg_name] = new_arg_value
--> 207 return func(*args, **kwargs)

File /local_disk0/.ephemeral_nfs/cluster_libraries/python/lib/python3.9/site-packages/statsmodels/stats/diagnostic.py:1098, in linear_reset(res, power, test_type, use_f, cov_type, cov_kwargs)
   1096 mod = mod_class(res.model.data.endog, aug_exog)
   1097 cov_kwargs = {} if cov_kwargs is None else cov_kwargs
-> 1098 res = mod.fit(cov_type=cov_type, cov_kwargs=cov_kwargs)
   1099 nrestr = aug_exog.shape[1] - exog.shape[1]
   1100 nparams = aug_exog.shape[1]

File /local_disk0/.ephemeral_nfs/cluster_libraries/python/lib/python3.9/site-packages/statsmodels/regression/linear_model.py:373, in RegressionModel.fit(self, method, cov_type, cov_kwds, use_t, **kwargs)
    370     self.df_resid = self.nobs - self.rank
    372 if isinstance(self, OLS):
--> 373     lfit = OLSResults(
    374         self, beta,
    375         normalized_cov_params=self.normalized_cov_params,
    376         cov_type=cov_type, cov_kwds=cov_kwds, use_t=use_t)
    377 else:
    378     lfit = RegressionResults(
    379         self, beta,
    380         normalized_cov_params=self.normalized_cov_params,
    381         cov_type=cov_type, cov_kwds=cov_kwds, use_t=use_t,
    382         **kwargs)

File /local_disk0/.ephemeral_nfs/cluster_libraries/python/lib/python3.9/site-packages/statsmodels/regression/linear_model.py:1653, in RegressionResults.__init__(self, model, params, normalized_cov_params, scale, cov_type, cov_kwds, use_t, **kwargs)
   1651             use_t = use_t_2
   1652         # TODO: warn or not?
-> 1653     self.get_robustcov_results(cov_type=cov_type, use_self=True,
   1654                                use_t=use_t, **cov_kwds)
   1655 for key in kwargs:
   1656     setattr(self, key, kwargs[key])

File /local_disk0/.ephemeral_nfs/cluster_libraries/python/lib/python3.9/site-packages/statsmodels/regression/linear_model.py:2570, in RegressionResults.get_robustcov_results(self, cov_type, use_t, **kwargs)
   2567     res.cov_params_default = getattr(self, 'cov_' + cov_type.upper())
   2568 elif cov_type.lower() == 'hac':
   2569     # TODO: check if required, default in cov_hac_simple
-> 2570     maxlags = kwargs['maxlags']
   2571     res.cov_kwds['maxlags'] = maxlags
   2572     weights_func = kwargs.get('weights_func', sw.weights_bartlett)
KeyError: 'maxlags'

**Note**: As you can see, there are many issues on our GitHub tracker, so it is very possible that your issue has been posted before. Please check first before submitting so that we do not have to handle and close duplicates.

**Note**: Please be sure you are using the latest released version of `statsmodels`, or a recent build of `main`. If your problem has been fixed in an unreleased version, you might be able to use `main` until a new release occurs. 

**Note**: If you are using a released version, have you verified that the bug exists in the main branch of this repository? It helps the limited resources if we know problems exist in the current main branch so that they do not need to check whether the code sample produces a bug in the next release.

</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output
We expect that the code would produce the p-value from linear RESET test with "HAC" as the covariance estimator.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

[paste the output of ``import statsmodels.api as sm; sm.show_versions()`` here below this line]
0.14.1
</details>



Hints:
**update** this does not apply here, new issue #9122


robust cov_type is not implemented for, AFAIR, any of the diagnostic and specification tests.

Main reason: When I implemented those I only had references (and other packages for unit tests) for the standard version with basic OLS in the auxiliary regression.

At the time I was not sure about the theory and did not find explicit references for using `cov_type`s, although it's mentioned e.g. by Wooldridge that we can just use a specification robust cov_type in the auxiliary regression.

Also for some LM specification tests, we use uncentered rsquared as computational shortcut and not a standard wald test, i.e. it hardcodes the nonrobust hypothesis test.

For some specification/diagnostic tests we might be able to just forward the cov_type information to the auxiliary regression,
even if we cannot write unit test to verify against other packages.
In other cases it might be easier to use generic conditional moment test and newer score_test/lm_test which AFAIR already allow for cov_type (although currently with some restrictions).

I add elevated priority to issue to look at this.
I worked a lot on misspecification robust (HC, ...) CMT, score/lm test in the last years (but mostly for GLM and discrete), and should be able to get back to this. 
 
Note: using wald test with augmented `exog` is easy to perform by just defining a new model.
e.g. one simple case for poisson https://gist.github.com/josef-pkt/1fdcc5c57cd63824fb3130d013e30103

Aside: 
The older literature in econometrics was emphasizing auxiliary regression because it was relatively easy to use existing models from statistical packages for this. It also relied in many cases on an OPG version which is much easier to compute but has worse small sample properties (mainly relevant for nonlinear or non-gaussian models)
For statsmodels it's easier and provides more options to implement these tests based on generic theory and generic methods.



my mistake
the reset test in `diagnostic` was written by @bashtage and has cov_type keywords. 
I found the bug

the reset test uses cov_kwargs, but the model.fit use cov_kwds
the correct code in the reset function should be
`res = mod.fit(cov_type=cov_type, cov_kwds=cov_kwargs)`

(I don't like the word kwargs, and prefer `xxx_kwds` for named options)

Is this problem been resolved yet? I'm using the latest version of statsmodels 0.14.2 and still have the exact same issue.

Patch:
diff --git a/statsmodels/stats/diagnostic.py b/statsmodels/stats/diagnostic.py
--- a/statsmodels/stats/diagnostic.py
+++ b/statsmodels/stats/diagnostic.py
@@ -237,8 +237,9 @@ def compare_j(results_x, results_z, store=False):
     return tstat, pval
 
 
+@deprecate_kwarg("cov_kwargs", "cov_kwds")
 def compare_encompassing(results_x, results_z, cov_type="nonrobust",
-                         cov_kwargs=None):
+                         cov_kwds=None):
     r"""
     Davidson-MacKinnon encompassing test for comparing non-nested models
 
@@ -253,7 +254,7 @@ def compare_encompassing(results_x, results_z, cov_type="nonrobust",
         OLS covariance estimator. Specify one of "HC0", "HC1", "HC2", "HC3"
         to use White's covariance estimator. All covariance types supported
         by ``OLS.fit`` are accepted.
-    cov_kwargs : dict, default None
+    cov_kwds : dict, default None
         Dictionary of covariance options passed to ``OLS.fit``. See OLS.fit
         for more details.
 
@@ -317,8 +318,8 @@ def _test_nested(endog, a, b, cov_est, cov_kwds):
         df_num, df_denom = int(test.df_num), int(test.df_denom)
         return stat, pvalue, df_num, df_denom
 
-    x_nested = _test_nested(y, x, z, cov_type, cov_kwargs)
-    z_nested = _test_nested(y, z, x, cov_type, cov_kwargs)
+    x_nested = _test_nested(y, x, z, cov_type, cov_kwds)
+    z_nested = _test_nested(y, z, x, cov_type, cov_kwds)
     return pd.DataFrame([x_nested, z_nested],
                         index=["x", "z"],
                         columns=["stat", "pvalue", "df_num", "df_denom"])
@@ -479,9 +480,10 @@ def acorr_ljungbox(x, lags=None, boxpierce=False, model_df=0, period=None,
                         index=lags)
 
 
+@deprecate_kwarg("cov_kwargs", "cov_kwds")
 @deprecate_kwarg("maxlag", "nlags")
 def acorr_lm(resid, nlags=None, store=False, *, period=None,
-             ddof=0, cov_type="nonrobust", cov_kwargs=None):
+             ddof=0, cov_type="nonrobust", cov_kwds=None):
     """
     Lagrange Multiplier tests for autocorrelation.
 
@@ -510,7 +512,7 @@ def acorr_lm(resid, nlags=None, store=False, *, period=None,
         OLS covariance estimator. Specify one of "HC0", "HC1", "HC2", "HC3"
         to use White's covariance estimator. All covariance types supported
         by ``OLS.fit`` are accepted.
-    cov_kwargs : dict, default None
+    cov_kwds : dict, default None
         Dictionary of covariance options passed to ``OLS.fit``. See OLS.fit for
         more details.
 
@@ -545,8 +547,8 @@ def acorr_lm(resid, nlags=None, store=False, *, period=None,
     """
     resid = array_like(resid, "resid", ndim=1)
     cov_type = string_like(cov_type, "cov_type")
-    cov_kwargs = {} if cov_kwargs is None else cov_kwargs
-    cov_kwargs = dict_like(cov_kwargs, "cov_kwargs")
+    cov_kwds = {} if cov_kwds is None else cov_kwds
+    cov_kwds = dict_like(cov_kwds, "cov_kwds")
     nobs = resid.shape[0]
     if period is not None and nlags is None:
         maxlag = min(nobs // 5, 2 * period)
@@ -563,7 +565,7 @@ def acorr_lm(resid, nlags=None, store=False, *, period=None,
     usedlag = maxlag
 
     resols = OLS(xshort, xdall[:, :usedlag + 1]).fit(cov_type=cov_type,
-                                                     cov_kwargs=cov_kwargs)
+                                                     cov_kwds=cov_kwds)
     fval = float(resols.fvalue)
     fpval = float(resols.f_pvalue)
     if cov_type == "nonrobust":
@@ -985,9 +987,10 @@ def het_goldfeldquandt(y, x, idx=None, split=None, drop=None,
     return fval, fpval, ordering
 
 
+@deprecate_kwarg("cov_kwargs", "cov_kwds")
 @deprecate_kwarg("result", "res")
 def linear_reset(res, power=3, test_type="fitted", use_f=False,
-                 cov_type="nonrobust", cov_kwargs=None):
+                 cov_type="nonrobust", cov_kwds=None):
     r"""
     Ramsey's RESET test for neglected nonlinearity
 
@@ -1015,7 +1018,7 @@ def linear_reset(res, power=3, test_type="fitted", use_f=False,
         OLS covariance estimator. Specify one of "HC0", "HC1", "HC2", "HC3"
         to use White's covariance estimator. All covariance types supported
         by ``OLS.fit`` are accepted.
-    cov_kwargs : dict, default None
+    cov_kwds : dict, default None
         Dictionary of covariance options passed to ``OLS.fit``. See OLS.fit
         for more details.
 
@@ -1055,7 +1058,7 @@ def linear_reset(res, power=3, test_type="fitted", use_f=False,
                          "non-constant column.")
     test_type = string_like(test_type, "test_type",
                             options=("fitted", "exog", "princomp"))
-    cov_kwargs = dict_like(cov_kwargs, "cov_kwargs", optional=True)
+    cov_kwds = dict_like(cov_kwds, "cov_kwds", optional=True)
     use_f = bool_like(use_f, "use_f")
     if isinstance(power, int):
         if power < 2:
@@ -1093,8 +1096,8 @@ def linear_reset(res, power=3, test_type="fitted", use_f=False,
     aug_exog = np.hstack([exog] + [aug ** p for p in power])
     mod_class = res.model.__class__
     mod = mod_class(res.model.data.endog, aug_exog)
-    cov_kwargs = {} if cov_kwargs is None else cov_kwargs
-    res = mod.fit(cov_type=cov_type, cov_kwargs=cov_kwargs)
+    cov_kwds = {} if cov_kwds is None else cov_kwds
+    res = mod.fit(cov_type=cov_type, cov_kwds=cov_kwds)
     nrestr = aug_exog.shape[1] - exog.shape[1]
     nparams = aug_exog.shape[1]
     r_mat = np.eye(nrestr, nparams, k=nparams-nrestr)


Test patch:
diff --git a/statsmodels/stats/tests/test_diagnostic.py b/statsmodels/stats/tests/test_diagnostic.py
--- a/statsmodels/stats/tests/test_diagnostic.py
+++ b/statsmodels/stats/tests/test_diagnostic.py
@@ -1805,8 +1805,8 @@ def test_encompasing_error(reset_randomstate):
 @pytest.mark.parametrize(
     "cov",
     [
-        dict(cov_type="nonrobust", cov_kwargs={}),
-        dict(cov_type="HC0", cov_kwargs={}),
+        dict(cov_type="nonrobust", cov_kwds={}),
+        dict(cov_type="HC0", cov_kwds={}),
     ],
 )
 def test_reset_smoke(power, test_type, use_f, cov, reset_randomstate):
@@ -1826,8 +1826,8 @@ def test_reset_smoke(power, test_type, use_f, cov, reset_randomstate):
 @pytest.mark.parametrize(
     "cov",
     [
-        dict(cov_type="nonrobust", cov_kwargs={}),
-        dict(cov_type="HC0", cov_kwargs={}),
+        dict(cov_type="nonrobust", cov_kwds={}),
+        dict(cov_type="HC0", cov_kwds={}),
     ],
 )
 def test_acorr_lm_smoke(store, ddof, cov, reset_randomstate):
@@ -1996,3 +1996,36 @@ def test_diagnostics_pandas(reset_randomstate):
         res, order_by=np.arange(y.shape[0] - 1, 0 - 1, -1)
     )
     smsdia.spec_white(res.resid, x)
+
+
+def test_deprecated_argument():
+    x = np.random.randn(100)
+    y = 2 * x + np.random.randn(100)
+    result = OLS(y, add_constant(x)).fit(
+        cov_type="HAC", cov_kwds={"maxlags": 2}
+    )
+    with pytest.warns(FutureWarning, match="the "):
+        smsdia.linear_reset(
+            result,
+            power=2,
+            test_type="fitted",
+            cov_type="HAC",
+            cov_kwargs={"maxlags": 2},
+        )
+
+
+def test_diagnostics_hac(reset_randomstate):
+    x = np.random.randn(100)
+    y = 2 * x + np.random.randn(100)
+    result = OLS(y, add_constant(x)).fit(
+        cov_type="HAC", cov_kwds={"maxlags": 2}
+    )
+    reset_test = smsdia.linear_reset(
+        result,
+        power=2,
+        test_type="fitted",
+        cov_type="HAC",
+        cov_kwds={"maxlags": 2},
+    )
+    assert reset_test.statistic > 0
+    assert 0 <= reset_test.pvalue <= 1
diff --git a/statsmodels/tsa/tests/test_seasonal.py b/statsmodels/tsa/tests/test_seasonal.py
--- a/statsmodels/tsa/tests/test_seasonal.py
+++ b/statsmodels/tsa/tests/test_seasonal.py
@@ -315,17 +315,27 @@ def test_seasonal_decompose_multiple():
 
 
 @pytest.mark.matplotlib
-@pytest.mark.parametrize('model', ['additive', 'multiplicative'])
-@pytest.mark.parametrize('freq', [4, 12])
-@pytest.mark.parametrize('two_sided', [True, False])
-@pytest.mark.parametrize('extrapolate_trend', [True, False])
-def test_seasonal_decompose_plot(model, freq, two_sided, extrapolate_trend):
+@pytest.mark.parametrize("model", ["additive", "multiplicative"])
+@pytest.mark.parametrize("freq", [4, 12])
+@pytest.mark.parametrize("two_sided", [True, False])
+@pytest.mark.parametrize("extrapolate_trend", [True, False])
+def test_seasonal_decompose_plot(
+    model, freq, two_sided, extrapolate_trend, close_figures
+):
     x = np.array([-50, 175, 149, 214, 247, 237, 225, 329, 729, 809,
                   530, 489, 540, 457, 195, 176, 337, 239, 128, 102,
                   232, 429, 3, 98, 43, -141, -77, -13, 125, 361, -45, 184])
     x -= x.min() + 1
     x2 = np.r_[x[12:], x[:12]]
     x = np.c_[x, x2]
-    res = seasonal_decompose(x, period=freq, two_sided=two_sided,
-                             extrapolate_trend=extrapolate_trend)
-    res.plot()
+    res = seasonal_decompose(
+        x,
+        period=freq,
+        two_sided=two_sided,
+        extrapolate_trend=extrapolate_trend
+    )
+    fig = res.plot()
+
+    import matplotlib.pyplot as plt
+
+    plt.close(fig)