diff --git a/fastcan/narx/_base.py b/fastcan/narx/_base.py
index 1b6d0eb..e4cf043 100644
--- a/fastcan/narx/_base.py
+++ b/fastcan/narx/_base.py
@@ -205,7 +205,7 @@ def fit(self, X, y, sample_weight=None, coef_init=None, **params):
         """
         none_inputs = False
         if X is None:  # Auto-regressive model
-            X = np.empty((1, 0), dtype=float)  # Skip validation
+            X = np.empty((1, 0), dtype=float, order="C")  # Skip validation
             none_inputs = True
         check_X_params = dict(
             dtype=float, order="C", ensure_all_finite="allow-nan", ensure_min_features=0
@@ -572,6 +572,7 @@ def predict(self, X, y_init=None):
                 y_init,
                 ensure_2d=False,
                 dtype=float,
+                order="C",
                 ensure_min_samples=0,
                 ensure_all_finite="allow-nan",
             )
diff --git a/fastcan/narx/tests/test_narx.py b/fastcan/narx/tests/test_narx.py
index a54cea8..9e18840 100644
--- a/fastcan/narx/tests/test_narx.py
+++ b/fastcan/narx/tests/test_narx.py
@@ -43,60 +43,68 @@ def test_time_ids():
 @pytest.mark.parametrize("nan", [False, True])
 def test_narx(nan, multi_output):
     """Test NARX"""
-    if multi_output:
-        rng = np.random.default_rng(12345)
-        n_samples = 1000
-        max_delay = 3
-        e0 = rng.normal(0, 0.1, n_samples)
-        e1 = rng.normal(0, 0.02, n_samples)
-        u0 = rng.uniform(0, 1, n_samples + max_delay)
-        u1 = rng.normal(0, 0.1, n_samples + max_delay)
-        y0 = np.zeros(n_samples + max_delay)
-        y1 = np.zeros(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]
-                + 0.3 * u0[i] ** 2
-                + 2 * u0[i - 1] * u0[i - 3]
-                + 1.5 * u0[i - 2] * u1[i - 3]
-                + 1
-            )
-            y1[i] = (
-                0.6 * y1[i - 1]
-                - 0.2 * y0[i - 1] * y1[i - 2]
-                + 0.3 * u1[i] ** 2
-                + 1.5 * u1[i - 2] * u0[i - 3]
-                + 0.5
-            )
-        y = np.c_[y0[max_delay:] + e0, y1[max_delay:] + e1]
-        X = np.c_[u0[max_delay:], u1[max_delay:]]
-        n_outputs = 2
-    else:
-        rng = np.random.default_rng(12345)
-        n_samples = 1000
-        max_delay = 3
-        e = rng.normal(0, 0.1, n_samples)
-        u0 = rng.uniform(0, 1, n_samples + max_delay)
-        u1 = rng.normal(0, 0.1, n_samples)
-        y = np.zeros(n_samples + max_delay)
-        for i in range(max_delay, n_samples + max_delay):
-            y[i] = (
-                0.5 * y[i - 1]
-                + 0.3 * u0[i] ** 2
-                + 2 * u0[i - 1] * u0[i - 3]
-                + 1.5 * u0[i - 2] * u1[i - max_delay]
-                + 1
-            )
-        y = y[max_delay:] + e
-        X = np.c_[u0[max_delay:], u1]
-        n_outputs = 1
 
-    if nan:
-        X_nan_ids = rng.choice(n_samples, 20, replace=False)
-        y_nan_ids = rng.choice(n_samples, 10, replace=False)
-        X[X_nan_ids] = np.nan
-        y[y_nan_ids] = np.nan
+    def make_data(multi_output, nan, rng):
+        if multi_output:
+            n_samples = 1000
+            max_delay = 3
+            e0 = rng.normal(0, 0.1, n_samples)
+            e1 = rng.normal(0, 0.02, n_samples)
+            u0 = rng.uniform(0, 1, n_samples + max_delay)
+            u1 = rng.normal(0, 0.1, n_samples + max_delay)
+            y0 = np.zeros(n_samples + max_delay)
+            y1 = np.zeros(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]
+                    + 0.3 * u0[i] ** 2
+                    + 2 * u0[i - 1] * u0[i - 3]
+                    + 1.5 * u0[i - 2] * u1[i - 3]
+                    + 1
+                )
+                y1[i] = (
+                    0.6 * y1[i - 1]
+                    - 0.2 * y0[i - 1] * y1[i - 2]
+                    + 0.3 * u1[i] ** 2
+                    + 1.5 * u1[i - 2] * u0[i - 3]
+                    + 0.5
+                )
+            y = np.c_[y0[max_delay:] + e0, y1[max_delay:] + e1]
+            X = np.c_[u0[max_delay:], u1[max_delay:]]
+            n_outputs = 2
+        else:
+            rng = np.random.default_rng(12345)
+            n_samples = 1000
+            max_delay = 3
+            e = rng.normal(0, 0.1, n_samples)
+            u0 = rng.uniform(0, 1, n_samples + max_delay)
+            u1 = rng.normal(0, 0.1, n_samples)
+            y = np.zeros(n_samples + max_delay)
+            for i in range(max_delay, n_samples + max_delay):
+                y[i] = (
+                    0.5 * y[i - 1]
+                    + 0.3 * u0[i] ** 2
+                    + 2 * u0[i - 1] * u0[i - 3]
+                    + 1.5 * u0[i - 2] * u1[i - max_delay]
+                    + 1
+                )
+            y = y[max_delay:] + e
+            X = np.c_[u0[max_delay:], u1]
+            n_outputs = 1
+
+        if nan:
+            X_nan_ids = rng.choice(n_samples, 20, replace=False)
+            y_nan_ids = rng.choice(n_samples, 10, replace=False)
+            X[X_nan_ids] = np.nan
+            y[y_nan_ids] = np.nan
+
+        X = np.asfortranarray(X)
+        y = np.asfortranarray(y)
+        return X, y, n_outputs
+
+    rng = np.random.default_rng(12345)
+    X, y, n_outputs = make_data(multi_output, nan, rng)
 
     if multi_output:
         narx_score = make_narx(
@@ -178,7 +186,7 @@ def test_narx(nan, multi_output):
     assert np.any(narx_osa_msa_coef != narx_array_init_msa.coef_)
 
     if multi_output:
-        y_init = np.ones((narx_array_init_msa.max_delay_, n_outputs))
+        y_init = np.ones((narx_array_init_msa.max_delay_, n_outputs), order="F")
     else:
         y_init = [1] * narx_array_init_msa.max_delay_
     y_hat = narx_array_init_msa.predict(X, y_init=y_init)