diff --git a/mlair/helpers/statistics.py b/mlair/helpers/statistics.py
index 87f780f9a6133edfcb2f9c71c2956b92f332e915..546c2ff71424676882c5d1b27c8ab807941d5ca8 100644
--- a/mlair/helpers/statistics.py
+++ b/mlair/helpers/statistics.py
@@ -481,3 +481,49 @@ class SkillScores:
 
         return monthly_mean
 
+
+def create_single_bootstrap_realization(data: xr.DataArray) -> xr.DataArray:
+    """
+    Return a bootstraped realization of data
+    :param data: data from which to draw ONE bootstrap realization
+    :return: bootstrapped realization of data
+    """
+
+    num_of_blocks = data.shape[0]
+    time_block_dim = data.dims[0]
+    boot_idx = np.random.choice(num_of_blocks, size=num_of_blocks, replace=True)
+    return data.isel({time_block_dim: 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
+    """
+    axis = kwargs.pop('axis', 0)
+    return data.mean(axis=axis, **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:
+    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)
+        ))
+    return res
+
+
+
+
+
+
+
+
+