diff --git a/inst/Regression/stepwisefit.m b/inst/Regression/stepwisefit.m
index 67c2359a8..1f8400169 100644
--- a/inst/Regression/stepwisefit.m
+++ b/inst/Regression/stepwisefit.m
@@ -139,15 +139,21 @@
     error ("stepwisefit: at least two input arguments required");
   endif
 
-  if (! ismatrix (X) || ! isvector (y))
-    error ("stepwisefit: X must be a matrix and y a vector");
+  if (! ismatrix (X))
+    error ("stepwisefit: X must be a matrix.");
+  endif
+  if (! (isvector (y) || isempty (y)))
+    error ("stepwisefit: Y must be a column vector.");
+  endif
+  if (! iscolumn (y))
+    error ("stepwisefit: Y must be a column vector.");
   endif
 
   y = y(:);
 
   ## Validate row compatibility BEFORE any concatenation
 if (rows (X) != rows (y))
-  error ("stepwisefit: X must be a matrix and y a vector");
+    error ("stepwisefit: X and Y must have the same number of rows.");
 endif
 
   ## Parse Name-Value pairs
@@ -189,7 +195,8 @@
 
   ## Validate PRemove (if provided)
   if (! isempty (PRemove))
-    if (! (isscalar (PRemove) && isnumeric (PRemove) && PRemove > 0 && PRemove < 1))
+    if (! (isscalar (PRemove) && isnumeric (PRemove) && PRemove > 0
+           && PRemove < 1))
       error ("stepwisefit: PRemove must be a scalar strictly between 0 and 1");
     endif
     if (PRemove < PEnter)
@@ -198,7 +205,8 @@
   endif
 
   ## Validate MaxIter
-  if (! (isscalar (MaxIter) && isnumeric (MaxIter) && MaxIter > 0 && fix (MaxIter) == MaxIter))
+  if (! (isscalar (MaxIter) && isnumeric (MaxIter) && MaxIter > 0
+         && fix (MaxIter) == MaxIter))
     error ("stepwisefit: MaxIter must be a positive integer");
   endif
 
@@ -210,6 +218,24 @@
   n = rows (Xc);
   p = columns (Xc);
 
+  if (n == 0)
+    b = NaN (p, 1);
+    se = NaN (p, 1);
+    pval = NaN (p, 1);
+    finalmodel = false (1, p);
+    stats = struct ('source', 'stepwisefit', 'df0', -1, 'dfe', 0, ...
+                    'SStotal', 0, 'SSresid', 0, 'fstat', NaN, 'pval', NaN, ...
+                    'rmse', NaN, 'intercept', 0, 'wasnan', wasnan, ...
+                    'xr', zeros(0, p), 'yr', zeros(0, 1), 'B', NaN (p, 1), ...
+                    'SE', NaN (p, 1), 'TSTAT', NaN (p, 1), ...
+                    'PVAL', NaN (p, 1), ...
+                    'covb', NaN (p, p));
+    history = struct ('B', zeros (p, 0), 'rmse', zeros (0, 0), ...
+                      'df0', zeros (0, 0), 'in', false (0, p));
+    nextstep = 0;
+    return;
+  endif
+
   ## Validate Keep and InModel type (if provided)
   if (! isempty (Keep) && ! islogical (Keep))
     error ("stepwisefit: Keep must be a logical vector");
@@ -649,22 +675,29 @@
 %! [b1] = stepwisefit (X, y);
 %! [b2] = stepwisefit (X, y, 'Scale', 'on');
 %! assert_equal (rows (b1) == rows (b2), true);
-%!test
-%! X = randn (20,4);
-%! y = randn (20,1);
-%! fail ('stepwisefit (X,y,''Keep'',[true false])');
+%!error <stepwisefit: Keep length must match number of predictors> ...
+%!       stepwisefit (randn (20,4), randn (20,1), 'Keep', [true false])
 
 ## Test input validation
 %!error <stepwisefit: at least two input arguments required> ...
 %!       stepwisefit ()
-%!error <stepwisefit: X must be a matrix and y a vector> ...
+%!error <stepwisefit: X must be a matrix.> ...
 %!       stepwisefit (ones (2,2,2), [1;2])
-%!error <stepwisefit: X must be a matrix and y a vector> ...
+%!error <stepwisefit: X and Y must have the same number of rows.> ...
 %!       stepwisefit (ones (3,2), ones (2,1))
 %!error <stepwisefit: unrecognized input arguments> ...
 %!       stepwisefit (randn (10,2), randn (10,1), 'UnknownOpt', 5)
 %!error <stepwisefit: Display must be 'on' or 'off'> ...
 %!       stepwisefit (randn (10,2), randn (10,1), 'Display', 'maybe')
+
+%!test
+%! ## Edge cases with empty arrays
+%! b = stepwisefit (zeros (0, 3), zeros (0, 1), 'Display', 'off');
+%! assert_equal (size (b), [3 1]);
+%! assert_equal (isnan (b), true (3, 1));
+%!
+%!error <stepwisefit: Y must be a column vector.> ...
+%!       stepwisefit ([], [], 'Display', 'off')
 %!error <stepwisefit: Scale must be 'on' or 'off'> ...
 %!       stepwisefit (randn (10,2), randn (10,1), 'Scale', 123)
 %!error <stepwisefit: PEnter must be a scalar strictly between 0 and 1> ...