diff --git a/mlair/helpers/statistics.py b/mlair/helpers/statistics.py
index 8b510f704396b35cf022089e7d94e037f4c62a2b..9a0346e4f09bbbe75d2c8dd70dac2d26d1b5b146 100644
--- a/mlair/helpers/statistics.py
+++ b/mlair/helpers/statistics.py
@@ -502,3 +502,59 @@ class SkillScores:
 
         return monthly_mean
 
+
+def create_single_bootstrap_realization(data: xr.DataArray, dim_name_time: str) -> xr.DataArray:
+    """
+    Return a bootstraped realization of data
+    :param data: data from which to draw ONE bootstrap realization
+    :param dim_name_time: name of time dimension
+    :return: bootstrapped realization of data
+    """
+
+    num_of_blocks = data.coords[dim_name_time].shape[0]
+    boot_idx = np.random.choice(num_of_blocks, size=num_of_blocks, replace=True)
+    return data.isel({dim_name_time: boot_idx})
+
+
+def calculate_average(data: xr.DataArray, **kwargs) -> xr.DataArray:
+    """
+    Calculate mean of data
+    :param data: data for which to calculate mean
+    :return: mean of data
+    """
+    return data.mean(**kwargs)
+
+
+def create_n_bootstrap_realizations(data: xr.DataArray, dim_name_time, dim_name_model, n_boots: int = 1000,
+                                    dim_name_boots='boots') -> xr.DataArray:
+    """
+    Create n bootstrap realizations and calculate averages across realizations
+
+    :param data: original data from which to create bootstrap realizations
+    :param dim_name_time: name of time dimension
+    :param dim_name_model: name of model dimension
+    :param n_boots: number of bootstap realizations
+    :param dim_name_boots: name of bootstap dimension
+    :return:
+    """
+    res_dims = [dim_name_boots]
+    dims = list(data.dims)
+    coords = {dim_name_boots: range(n_boots), dim_name_model: data.coords[dim_name_model] }
+    if len(dims) > 1:
+        res_dims = res_dims + dims[1:]
+    res = xr.DataArray(np.nan, dims=res_dims, coords=coords)
+    for boot in range(n_boots):
+        res[boot] = (calculate_average(
+            create_single_bootstrap_realization(data, dim_name_time=dim_name_time),
+            dim=dim_name_time
+        ))
+    return res
+
+
+
+
+
+
+
+
+