diff --git a/ar-python.cpp b/ar-python.cpp
index d9586a9..6cee86a 100644
--- a/ar-python.cpp
+++ b/ar-python.cpp
@@ -107,7 +107,7 @@ static const char ar_arsel_docstring[] =
 // Return type for ar_arsel initialized within initar
 static PyTypeObject *ar_ArselType = NULL;
 
-extern "C" PyObject *ar_arsel(PyObject *self, PyObject *args)
+extern "C" PyObject *ar_arsel(PyObject *self, PyObject *args, PyObject *kwargs)
 {
     PyObject *ret_args = NULL, *ret = NULL;
 
@@ -126,12 +126,24 @@ extern "C" PyObject *ar_arsel(PyObject *self, PyObject *args)
     std::size_t minorder  = DEFAULT_MINORDER;
     std::size_t maxorder  = DEFAULT_MAXORDER;
 
-    // Parse input tuple with second and subsequent arguments optional
+    // Define keyword argument names
+    static char *kwlist[] = {
+        (char *)"data",
+        (char *)"submean",
+        (char *)"absrho",
+        (char *)"criterion",
+        (char *)"minorder",
+        (char *)"maxorder",
+        NULL
+    };
+
+    // Parse input tuple and keywords with second and subsequent arguments optional
     {
         unsigned long ul_minorder = minorder, ul_maxorder = maxorder;
-        if (!PyArg_ParseTuple(args, "O|iiskk", &data_obj, &submean, &absrho,
-                                               &criterion,
-                                               &ul_minorder, &ul_maxorder)) {
+        if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O|iiskk", kwlist,
+                                         &data_obj, &submean, &absrho,
+                                         &criterion,
+                                         &ul_minorder, &ul_maxorder)) {
             return NULL;
         }
         minorder = ul_minorder;
@@ -334,7 +346,7 @@ extern "C" PyObject *ar_arsel(PyObject *self, PyObject *args)
 
 // Specification of methods available in the module
 static PyMethodDef ar_methods[] = {
-    {"arsel", ar_arsel, METH_VARARGS, ar_arsel_docstring},
+    {"arsel", (PyCFunction)ar_arsel, METH_VARARGS | METH_KEYWORDS, ar_arsel_docstring},
     {NULL, NULL, 0, NULL}
 };
 
diff --git a/test_ar.py b/test_ar.py
index 295bc00..eb444dc 100755
--- a/test_ar.py
+++ b/test_ar.py
@@ -18,8 +18,8 @@ def main():
     expected = [float(x) for x in open("test0.coeff").read().split()]
     data = [float(x) for x in open("test0.dat").read().split()]
 
-    # Fit AR model using ar.arsel
-    result = ar.arsel(
+    # Fit AR model using ar.arsel with POSITIONAL arguments
+    result_positional = ar.arsel(
         data,
         False,
         True,
@@ -27,6 +27,48 @@ def main():
         len(expected),
         len(expected),
     )
+
+    # Fit AR model using ar.arsel with KEYWORD arguments
+    result_keyword = ar.arsel(
+        data=data,
+        submean=False,
+        absrho=True,
+        criterion="CIC",
+        minorder=len(expected),
+        maxorder=len(expected),
+    )
+
+    # Verify that positional and keyword calls produce identical results
+    print("Verifying positional vs keyword arguments produce identical results...")
+
+    def values_equal(a, b):
+        """Check if two values are equal, handling NaN and inf correctly."""
+        # If both are NaN, consider them equal
+        if np.isnan(a) and np.isnan(b):
+            return True
+        # If both are inf with same sign, consider them equal
+        if np.isinf(a) and np.isinf(b) and (a > 0) == (b > 0):
+            return True
+        # Otherwise use normal equality
+        return a == b
+
+    assert len(result_positional.AR[0]) == len(result_keyword.AR[0]), \
+        "AR length mismatch between positional and keyword calls"
+    for i, (pos_val, kw_val) in enumerate(zip(result_positional.AR[0], result_keyword.AR[0])):
+        assert values_equal(pos_val, kw_val), \
+            f"AR[0][{i}] mismatch: positional={pos_val}, keyword={kw_val}"
+    assert values_equal(result_positional.mu[0], result_keyword.mu[0]), \
+        f"mu mismatch: positional={result_positional.mu[0]}, keyword={result_keyword.mu[0]}"
+    assert values_equal(result_positional.sigma2eps[0], result_keyword.sigma2eps[0]), \
+        f"sigma2eps mismatch: positional={result_positional.sigma2eps[0]}, keyword={result_keyword.sigma2eps[0]}"
+    assert values_equal(result_positional.gain[0], result_keyword.gain[0]), \
+        f"gain mismatch: positional={result_positional.gain[0]}, keyword={result_keyword.gain[0]}"
+    assert values_equal(result_positional.sigma2x[0], result_keyword.sigma2x[0]), \
+        f"sigma2x mismatch: positional={result_positional.sigma2x[0]}, keyword={result_keyword.sigma2x[0]}"
+    print("✓ Positional and keyword arguments produce IDENTICAL results\n")
+
+    # Use positional result for the rest of the test
+    result = result_positional
     estimated = result.AR[0][1:]
 
     if len(estimated) != len(expected):