diff --git a/pyfacts/statistics.py b/pyfacts/statistics.py
index 0b133f5..7de44b7 100644
--- a/pyfacts/statistics.py
+++ b/pyfacts/statistics.py
@@ -8,7 +8,7 @@ from typing import Literal
 from pyfacts.core import date_parser
 
 from .pyfacts import TimeSeries
-from .utils import _interval_to_years
+from .utils import _interval_to_years, covariance
 
 
 @date_parser(3, 4)
@@ -212,7 +212,7 @@ def beta(
     asset_rr = asset_data.calculate_rolling_returns(**common_params)
     market_rr = market_data.calculate_rolling_returns(**common_params)
 
-    cov = statistics.covariance(asset_rr.values, market_rr.values)
+    cov = covariance(asset_rr.values, market_rr.values)
     market_var = statistics.variance(market_rr.values)
 
     beta = cov / market_var
diff --git a/pyfacts/utils.py b/pyfacts/utils.py
index 857d217..11522f0 100644
--- a/pyfacts/utils.py
+++ b/pyfacts/utils.py
@@ -1,4 +1,7 @@
+from __future__ import annotations
+
 import datetime
+import statistics
 from dataclasses import dataclass
 from typing import List, Literal, Mapping, Sequence, Tuple
 
@@ -187,3 +190,36 @@ def _is_eomonth(dates: Sequence[datetime.datetime], threshold: float = 0.7):
     eomonth_dates = [date.month != (date + relativedelta(days=1)).month for date in dates]
     eomonth_proportion = sum(eomonth_dates) / len(dates)
     return eomonth_proportion > threshold
+
+
+def covariance(series1: list, series2: list) -> float:
+    """Returns the covariance of two series
+
+        This is a compatibility function for Python versions prior to 3.10.
+        It will be replaced with statistics.covariance when support is dropped for versions <3.10.
+
+    Parameters
+    ----------
+    series1 : List
+        A list of numbers
+    series2 : list
+        A list of numbers
+
+    Returns
+    -------
+    float
+        Returns the covariance as a float value
+    """
+
+    n = len(series1)
+    if len(series2) != n:
+        raise ValueError("Lenght of both series must be same for covariance calcualtion.")
+    if n < 2:
+        raise ValueError("At least two data poitns are required for covariance calculation.")
+
+    mean1 = statistics.mean(series1)
+    mean2 = statistics.mean(series2)
+
+    xy = sum([(x - mean1) * (y - mean2) for x, y in zip(series1, series2)])
+
+    return xy / n