diff --git a/mlair/helpers/statistics.py b/mlair/helpers/statistics.py
index 57143c9aed0e730c81adbed33f7ba62fea39b298..b0ecadb12310a179b6446d4340804ae56d8f39c7 100644
--- a/mlair/helpers/statistics.py
+++ b/mlair/helpers/statistics.py
@@ -211,6 +211,14 @@ def mean_absolute_error(a, b, dim=None):
     return np.abs(a - b).mean(dim)
 
 
+def index_of_agreement(a, b, dim=None):
+    """Calculate index of agreement (IOA) where a is the forecast and b the reference (e.g. observation)."""
+    num = np.sum(np.square(b - a))
+    b_mean = (b * np.ones(1)).mean(dim)
+    den = np.sum(np.square(np.abs(b - b_mean) + np.abs(a - b_mean)))
+    return 1 - num / den
+
+
 def calculate_error_metrics(a, b, dim):
     """Calculate MSE, RMSE, and MAE. Additionally return number of used values for calculation."""
     mse = mean_squared_error(a, b, dim)
diff --git a/test/test_helpers/test_statistics.py b/test/test_helpers/test_statistics.py
index a3f645937258604c2dbbda07b36a58d83e879065..58432ca4a4756df17c3f9ede50b35b0ff63bf43f 100644
--- a/test/test_helpers/test_statistics.py
+++ b/test/test_helpers/test_statistics.py
@@ -6,7 +6,7 @@ import xarray as xr
 from mlair.helpers.statistics import standardise, standardise_inverse, standardise_apply, centre, centre_inverse, \
     centre_apply, apply_inverse_transformation, min_max, min_max_inverse, min_max_apply, log, log_inverse, log_apply, \
     create_single_bootstrap_realization, calculate_average, create_n_bootstrap_realizations, mean_squared_error, \
-    mean_absolute_error, calculate_error_metrics
+    mean_absolute_error, calculate_error_metrics, index_of_agreement
 from mlair.helpers.testing import check_nested_equality
 
 lazy = pytest.lazy_fixture
@@ -297,6 +297,30 @@ class TestMeanAbsoluteError:
         assert xr.testing.assert_equal(mean_absolute_error(x_array1, x_array2, "value"), expected) is None
 
 
+class TestIndexOfAgreement:
+
+    def test_index_of_agreement(self):
+        d1 = np.array([1, 2, 3, 4, 5])
+        d2 = np.array([1, 2, 3, 4, 5])
+        assert index_of_agreement(d1, d2) == 1
+        d1 = np.array([1, 2, 3, 4, 7])
+        assert np.testing.assert_almost_equal(index_of_agreement(d1, d2), 0.9333, 3) is None
+        d1 = np.array([3, 4, 5, 6, 7])
+        assert np.testing.assert_almost_equal(index_of_agreement(d1, d2), 0.687, 3) is None
+
+    def test_index_of_agreement_xarray(self):
+        d1 = np.array([np.array([1, 2, 3, 4, 5]), np.array([1, 2, 3, 4, 5]), np.array([1, 2, 3, 4, 5])])
+        d2 = np.array([np.array([2, 4, 3, 4, 6]), np.array([2, 3, 3, 4, 5]), np.array([0, 1, 3, 4, 5])])
+        shape = d1.shape
+        coords = {'index': range(shape[0]), 'value': range(shape[1])}
+        x_array1 = xr.DataArray(d1, coords=coords, dims=coords.keys())
+        x_array2 = xr.DataArray(d2, coords=coords, dims=coords.keys())
+        res = index_of_agreement(x_array1, x_array1, dim="value")
+        assert np.testing.assert_almost_equal(res, 1) is None
+        res = index_of_agreement(x_array1, x_array2, dim="value")
+        assert np.testing.assert_almost_equal(res, 0.919, 2) is None
+
+
 class TestCalculateErrorMetrics:
 
     def test_calculate_error_metrics(self):