diff --git a/fastcan/narx.py b/fastcan/narx.py
index 2883c7a..4ee5ccc 100644
--- a/fastcan/narx.py
+++ b/fastcan/narx.py
@@ -669,7 +669,7 @@ def fit(self, X, y, sample_weight=None, coef_init=None, **params):
 
         Parameters
         ----------
-        X : {array-like, sparse matrix} of shape (n_samples, `n_features_in_`)
+        X : {array-like, sparse matrix} of shape (n_samples, `n_features_in_`) or None
             Training data.
 
         y : array-like of shape (n_samples,) or (n_samples, `n_outputs_`)
@@ -700,7 +700,9 @@ def fit(self, X, y, sample_weight=None, coef_init=None, **params):
         self : object
             Fitted Estimator.
         """
-        check_X_params = dict(dtype=float, order="C", ensure_all_finite="allow-nan")
+        check_X_params = dict(
+            dtype=float, order="C", ensure_all_finite="allow-nan", ensure_min_features=0
+        )
         check_y_params = dict(
             ensure_2d=False, dtype=float, order="C", ensure_all_finite="allow-nan"
         )
@@ -717,7 +719,10 @@ def fit(self, X, y, sample_weight=None, coef_init=None, **params):
         n_samples, n_features = X.shape
 
         if self.feat_ids is None:
-            feat_ids_ = make_poly_ids(n_features, 1) - 1
+            if n_features == 0:
+                feat_ids_ = make_poly_ids(self.n_outputs_, 1) - 1
+            else:
+                feat_ids_ = make_poly_ids(n_features, 1) - 1
         else:
             feat_ids_ = self.feat_ids
 
@@ -1152,6 +1157,7 @@ def predict(self, X, y_init=None):
             order="C",
             reset=False,
             ensure_all_finite="allow-nan",
+            ensure_min_features=0,
         )
         if y_init is None:
             y_init = np.zeros((self.max_delay_, self.n_outputs_))
@@ -1419,7 +1425,13 @@ def make_narx(
     |  0  | X[k-1,0]*X[k-3,0]  |  2.000   |
     |  0  |  X[k-2,0]*X[k,1]   |  1.528   |
     """
-    X = check_array(X, dtype=float, ensure_2d=True, ensure_all_finite="allow-nan")
+    X = check_array(
+        X,
+        dtype=float,
+        ensure_2d=True,
+        ensure_all_finite="allow-nan",
+        ensure_min_features=0,
+    )
     y = check_array(y, dtype=float, ensure_2d=False, ensure_all_finite="allow-nan")
     check_consistent_length(X, y)
     if y.ndim == 1:
diff --git a/tests/test_narx.py b/tests/test_narx.py
index 38bca0c..c72af2d 100644
--- a/tests/test_narx.py
+++ b/tests/test_narx.py
@@ -21,8 +21,12 @@
 
 
 def test_narx_is_sklearn_estimator():
+    # Skip 0 feature check for NARX, as AR models have no features
+    expected_failures = {
+        "check_estimators_empty_data_messages": ("NARX can handle 0 feature."),
+    }
     with pytest.warns(UserWarning, match="output_ids got"):
-        check_estimator(NARX())
+        check_estimator(NARX(), expected_failed_checks=expected_failures)
 
 
 def test_poly_ids():
@@ -576,3 +580,38 @@ def test_nan_split(max_delay):
     assert poly_terms_masked.shape[0] == n_sessions * (
         n_samples_per_session - narx.max_delay_
     )
+
+
+def test_default_narx_handles_zero_features():
+    """Check that default NARX handles X with 0 features without error."""
+    X = np.empty((10, 0))
+    y = np.random.rand(10, 1)
+    NARX().fit(X, y)
+
+
+def test_auto_reg():
+    """Test auto-regression with NARX"""
+    rng = np.random.default_rng(12345)
+    n_samples = 100
+    max_delay = 2
+    e0 = rng.normal(0, 0.01, n_samples)
+    e1 = rng.normal(0, 0.01, n_samples)
+    y0 = np.ones(n_samples + max_delay)
+    y1 = np.ones(n_samples + max_delay)
+    for i in range(max_delay, n_samples + max_delay):
+        y0[i] = 0.5 * y0[i - 1] + 0.8 * y1[i - 1] + 1
+        y1[i] = 0.6 * y1[i - 1] - 0.2 * y0[i - 1] * y1[i - 2] + 0.5
+    y = np.c_[y0[max_delay:] + e0, y1[max_delay:] + e1]
+    X = np.empty((n_samples, 0))  # No features, only auto-regression
+
+    model = make_narx(
+        X,
+        y,
+        n_terms_to_select=2,
+        max_delay=max_delay,
+        poly_degree=2,
+        verbose=0,
+    )
+    model.fit(X, y)
+    y_pred = model.predict(X, y_init=y[: model.max_delay_])
+    assert r2_score(y, model.predict(X, y_init=y)) > 0.5