diff --git a/.github/workflows/unittests.yml b/.github/workflows/unittests.yml
index 9bdf50e0..253f72b1 100644
--- a/.github/workflows/unittests.yml
+++ b/.github/workflows/unittests.yml
@@ -21,3 +21,13 @@ jobs:
           pip install smefit
       - name: Run tests
         run: pytest
+      - name: Install smefit
+        run: |
+          pip install smefit
+      - name: Test with pytest
+        run: |
+          # add test packages
+          pip install pytest
+          pip install pytest-cov
+          pip install pytest-env
+          pytest
diff --git a/tests/test_covariance.py b/tests/test_covariance.py
new file mode 100644
index 00000000..2ea8091d
--- /dev/null
+++ b/tests/test_covariance.py
@@ -0,0 +1,51 @@
+# python
+import pathlib, os
+import numpy as np
+import pytest
+from smefit.loader import Loader
+from smefit.covmat import covmat_from_systematics
+
+BASE_DIR = os.path.dirname(__file__)
+COMMONDATA_DIRS = [
+    pathlib.Path(BASE_DIR + "/../commondata"),
+    pathlib.Path(BASE_DIR + "/../commondata_projections_L0"),
+]
+
+
+def get_yaml_files(directories):
+    files = []
+    for directory in directories:
+        if pathlib.Path(directory).is_dir():
+            files.extend(
+                (f.stem, directory)
+                for f in pathlib.Path(directory).iterdir()
+                if f.suffix == ".yaml"
+            )
+    return files
+
+
+@pytest.mark.parametrize("dataset,commondata_dir", get_yaml_files(COMMONDATA_DIRS))
+def test_experimental_covmat(dataset, commondata_dir):
+    Loader.commondata_path = commondata_dir
+    Loader.theory_path = pathlib.Path(BASE_DIR + "/../theory")
+
+    loaded_dataset = Loader(
+        setname=dataset,
+        operators_to_keep=[],
+        order="LO",
+        use_quad="False",
+        use_theory_covmat="False",
+        use_multiplicative_prescription="False",
+        rot_to_fit_basis=None,
+        cutoff_scale=None,
+    )
+
+    stat_error = loaded_dataset.stat_error
+    sys_error = loaded_dataset.sys_error
+
+    exp_covmat = covmat_from_systematics([stat_error], [sys_error])
+
+    eigvals = np.linalg.eigvalsh(exp_covmat)
+    assert np.all(
+        eigvals > 0
+    ), f"Experimental covariance matrix for {dataset} is not positive definite."