Fix small sample correction handling for HAC (#1045)

* new clean branch for the draft PR

* clarify origin of default lags

* multiple smaller fixes

* commit to allow merge

* updated NW Meat function to match fixest implementation

* Newey-West complete. Tests and documentation pending

* pass function attributes along, force _nw_meat to work on arrays only (else numba does not compile)

* add basuc unit tests for HAC

* fix test bug?

* one more fix

* fix more bugs

* fix tests

* pacify mypy'

* test for errors

* add test for update

* update type hints for lags

* rename arg from lags to lag for compatibility with fixest

* require lag, time_id for HAC

* panel HAC first PR

* added balance function

* allow for non-balances panel and non-consecutive data

* _dk_meat added. documentation and tests are pending

* fix indentation error & njit & linter & mypy

* prepare tests for non-consecutive implementation

* good speed ups of time series hac

* fixed everything and added fixef_maxiter everywhere

* more numba for DK

* code reorg + check for unique time and panel

* more adjustment to new syntax

* update get_panel_idx

* tests for fepois, glm

* ar 1 errors

* add unrelated test error

* fix test bugs

* compure panel hac on scores, fix bugs for glms, support glms

* update readme

* update tests

* fix linter

* some progress ssc

* pre commit

* fix HAC ssc

* fix

* delete print statements

---------

Co-authored-by: Daman Dhaliwal <[email protected]>

Author: s3alfisc

pyfixest/estimation/feols_.py

@@ -327,7 +327,7 @@ def __init__(
         self._demean_func = impl["demean"]
         self._find_collinear_variables_func = impl["collinear"]
         self._crv1_meat_func = impl["crv1_meat"]
-        self._cound_nested_fixef_func = impl["nonnested"]
+        self._count_nested_fixef_func = impl["nonnested"]
 
         # set in get_fit()
         self._tZX = np.array([])
@@ -690,32 +690,18 @@ def vcov(
             self._vcov = self._ssc * self._vcov_hetero()
 
         elif self._vcov_type == "HAC":
-            if self._ssc_dict["k_adj"]:
-                self._ssc_dict["k_adj"] = False
-                warnings.warn(
-                    "k_adj was set to False for HAC inference because that's currently the only supported option."
-                )
-            if self._ssc_dict["G_adj"]:
-                self._ssc_dict["G_adj"] = False
-                warnings.warn(
-                    "G_adj was set to False for HAC inference because that's currently the only supported option."
-                )
-            if self._ssc_dict["k_fixef"] == "nonnested":
-                self._ssc_dict["k_fixef"] = "none"
-                warnings.warn(
-                    "k_fixef was set to none for HAC inference because nonnested is currently not supported."
-                )
-
             ssc_kwargs_hac = {
-                "k_fe_nested": 0,
-                "n_fe_fully_nested": 0,
+                "k_fe_nested": 0,  # nesting ignored / irrelevant for HAC SEs
+                "n_fe_fully_nested": 0,  # nesting ignored / irrelevant for HAC SEs
                 "vcov_sign": 1,
                 "vcov_type": "HAC",
-                "G": self._N,
+                "G": np.unique(self._data[self._time_id]).shape[
+                    0
+                ],  # number of unique time periods T used
             }
 
             all_kwargs = {**ssc_kwargs, **ssc_kwargs_hac}
-            self._ssc, self._dof_k, self._df_t = get_ssc(**all_kwargs)
+            self._ssc, self._df_k, self._df_t = get_ssc(**all_kwargs)
 
             self._vcov = self._ssc * self._vcov_hac()
 
@@ -770,7 +756,7 @@ def vcov(
             k_fe_nested = 0
             n_fe_fully_nested = 0
             if self._has_fixef and self._ssc_dict["k_fixef"] == "nonnested":
-                k_fe_nested_flag, n_fe_fully_nested = self._cound_nested_fixef_func(
+                k_fe_nested_flag, n_fe_fully_nested = self._count_nested_fixef_func(
                     all_fixef_array=np.array(
                         self._fixef.replace("^", "_").split("+"), dtype=str
                     ),

pyfixest/utils/utils.py

@@ -218,7 +218,7 @@ def get_ssc(
         adj_value = (N - 1) / (N - df_k) if vcov_type != "hetero" else N / (N - df_k)
 
     # G_adj applied with G = N for hetero but not for iid
-    if vcov_type in ["CRV"] and G_adj:
+    if vcov_type in ["CRV", "HAC"] and G_adj:
         if G_df == "conventional":
             G_adj_value = G / (G - 1)
         elif G_df == "min":
@@ -227,8 +227,7 @@ def get_ssc(
         else:
             raise ValueError("G_df is neither conventional nor min.")
 
-    df_t = N - df_k if vcov_type in ["iid", "hetero"] else G - 1
-
+    df_t = N - df_k if vcov_type in ["iid", "hetero", "HAC-TS"] else G - 1
     return np.array([adj_value * G_adj_value * vcov_sign]), df_k, df_t
 
 

tests/test_errors.py

@@ -1069,13 +1069,3 @@ def test_errors_hac():
                 vcov=vcov,
                 vcov_kwargs={"time_id": "time3", "panel_id": "panel", "lag": 5},
             )
-
-
-def test_errors_hac_inference():
-    data = pf.get_data()
-    data["Y"] = np.where(data["Y"] > 0, 1, 0)
-    with pytest.raises(
-        NotImplementedError,
-        match=r"HAC inference is not supported for this model type\.",
-    ):
-        pf.quantreg("Y ~ X1", data=data, vcov="NW")

tests/test_hac_vs_fixest.py

@@ -271,6 +271,9 @@ def _get_r_panel_kwargs(time_id, panel_id, lag, inference):
     [None, "weights"],
 )
 @pytest.mark.parametrize("fml", ols_fmls)
+@pytest.mark.parametrize("k_adj", [True, False])
+@pytest.mark.parametrize("G_adj", [True, False])
+@pytest.mark.parametrize("k_fixef", ["none", "nonnested", "full"])
 def test_single_fit_feols_hac_panel(
     data_panel,
     data_time,
@@ -279,11 +282,10 @@ def test_single_fit_feols_hac_panel(
     weights,
     fml,
     balanced,
+    k_adj,
+    G_adj,
+    k_fixef,
 ):
-    k_adj = False
-    G_adj = False
-    ssc_ = ssc(k_adj=k_adj, G_adj=G_adj)
-
     lag = vcov_kwargs.get("lag", None)
     time_id = vcov_kwargs.get("time_id", None)
     panel_id = vcov_kwargs.get("panel_id", None)
@@ -310,9 +312,9 @@ def test_single_fit_feols_hac_panel(
         if inference == "NW"
         else fixest.vcov_DK(**r_panel_kwars),
         data=data,
-        ssc=fixest.ssc(k_adj, "nested", False, G_adj, "min", "min"),
         **({"weights": ro.Formula(f"~{weights}")} if weights is not None else {}),
         panel_time_step=1,
+        ssc=fixest.ssc(k_adj, k_fixef, False, G_adj, "min", "min"),
     )
 
     mod = pf.feols(
@@ -321,7 +323,7 @@ def test_single_fit_feols_hac_panel(
         vcov=inference,
         vcov_kwargs=vcov_kwargs,
         weights=weights,
-        ssc=ssc_,
+        ssc=pf.ssc(k_adj=k_adj, k_fixef=k_fixef, G_adj=G_adj),
     )
 
     # r_fixest to global r env, needed for