diff --git a/sysidentpy/basis_function/_polynomial.py b/sysidentpy/basis_function/_polynomial.py
index 165b70a3..9cd2a018 100644
--- a/sysidentpy/basis_function/_polynomial.py
+++ b/sysidentpy/basis_function/_polynomial.py
@@ -26,6 +26,8 @@ class Polynomial(BaseBasisFunction):
     ----------
     degree : int (max_degree), default=2
         The maximum degree of the polynomial features.
+    include_bias : bool, default=True
+        Whether to include the bias (constant) term in the output feature matrix.
 
     Notes
     -----
@@ -37,8 +39,10 @@ class Polynomial(BaseBasisFunction):
     def __init__(
         self,
         degree: int = 2,
+        include_bias: bool = True,
     ):
         self.degree = degree
+        self.include_bias = include_bias
 
     def fit(
         self,
@@ -94,6 +98,11 @@ def fit(
             ]
         )
         psi = psi[max_lag:, :]
+
+        if self.include_bias:
+            bias_column = np.ones((psi.shape[0], 1))
+            psi = np.hstack((bias_column, psi))
+
         return psi
 
     def transform(
diff --git a/sysidentpy/basis_function/tests/test_basis_functions.py b/sysidentpy/basis_function/tests/test_basis_functions.py
index 08bfe253..79b0c968 100644
--- a/sysidentpy/basis_function/tests/test_basis_functions.py
+++ b/sysidentpy/basis_function/tests/test_basis_functions.py
@@ -5,7 +5,29 @@
 
 
 def test_fit_polynomial():
-    basis_function = Polynomial(degree=2)
+    basis_function = Polynomial(degree=2, include_bias=True)
+    data = np.array(([1, 1, 1], [2, 3, 4], [3, 3, 3]))
+    max_lag = 1
+    output = np.array([[1, 4, 6, 8, 9, 12, 16], [1, 9, 9, 9, 9, 9, 9]])
+
+    r = basis_function.fit(data=data, max_lag=max_lag)
+
+    assert_array_equal(output, r)
+
+
+def test_fit_polynomial_with_bias():
+    basis_function = Polynomial(degree=2, include_bias=True)
+    data = np.array(([1, 1, 1], [2, 3, 4], [3, 3, 3]))
+    max_lag = 1
+    output = np.array([[1, 4, 6, 8, 9, 12, 16], [1, 9, 9, 9, 9, 9, 9]])
+
+    r = basis_function.fit(data=data, max_lag=max_lag)
+
+    assert_array_equal(output, r)
+
+
+def test_fit_polynomial_without_bias():
+    basis_function = Polynomial(degree=2, include_bias=False)
     data = np.array(([1, 1, 1], [2, 3, 4], [3, 3, 3]))
     max_lag = 1
     output = np.array([[4, 6, 8, 9, 12, 16], [9, 9, 9, 9, 9, 9]])
@@ -16,11 +38,11 @@ def test_fit_polynomial():
 
 
 def test_fit_polynomial_predefined():
-    basis_function = Polynomial(degree=2)
+    basis_function = Polynomial(degree=2, include_bias=True)
     data = np.array(([1, 1, 1], [2, 3, 4], [3, 3, 3]))
     max_lag = 1
     predefined_regressors = np.array([0, 2, 4])
-    output = np.array([[4, 8, 12], [9, 9, 9]])
+    output = np.array([[1, 4, 8, 12], [1, 9, 9, 9]])
 
     r = basis_function.fit(
         data=data, max_lag=max_lag, predefined_regressors=predefined_regressors
@@ -30,10 +52,10 @@ def test_fit_polynomial_predefined():
 
 
 def test_transform_polynomial():
-    basis_function = Polynomial(degree=2)
+    basis_function = Polynomial(degree=2, include_bias=True)
     data = np.array(([1, 1, 1], [2, 3, 4], [3, 3, 3]))
     max_lag = 1
-    output = np.array([[4, 6, 8, 9, 12, 16], [9, 9, 9, 9, 9, 9]])
+    output = np.array([[1, 4, 6, 8, 9, 12, 16], [1, 9, 9, 9, 9, 9, 9]])
 
     r = basis_function.transform(data=data, max_lag=max_lag)