diff --git a/example_scripts/exec_timeagg.R b/example_scripts/exec_timeagg.R
index 66df62b868b25564f863970c97473c047d3534b3..d383745a96af8fdd884934827b3201d19dd58d21 100644
--- a/example_scripts/exec_timeagg.R
+++ b/example_scripts/exec_timeagg.R
@@ -26,7 +26,7 @@ data <- Loading(recipe)
 # Units transformation
 data <- Units(recipe, data)
 # Temporal aggregation
-data <- Aggregation(recipe = recipe, data = data)
+data_agg <- Aggregation(recipe = recipe, data = data)
 # Calibrate datasets
 data <- Calibration(recipe, data)
 # Compute skill metrics 
diff --git a/modules/Aggregation/R/agg_ini_end.R b/modules/Aggregation/R/agg_ini_end.R
index 17c6940a7a77ccf65f42c1fc110856d171f9bee1..400af417de3e291bab70a60350dba8aa4c441f69 100644
--- a/modules/Aggregation/R/agg_ini_end.R
+++ b/modules/Aggregation/R/agg_ini_end.R
@@ -9,6 +9,8 @@ agg_ini_end <- function(x, ini, end, indices = NULL, method, na.rm ,ncores) {
     # take the first and last element of each indices list for time_bounds saving
     ini <- unlist(lapply(indices, function(x){x[1]}))
     end <- unlist(lapply(indices, function(x){x[length(x)]}))
+    indices <- lapply(1:length(ini), function(x) {
+                    ini[x]:end[x]})
     plotting_attr <- list(names(indices))
   }
   original_dims <- names(dim(x[[1]]$data))
@@ -16,8 +18,13 @@ agg_ini_end <- function(x, ini, end, indices = NULL, method, na.rm ,ncores) {
     x[[1]]$data <- Apply(x[[1]]$data,
                          target_dim = 'time',
                          function(y, ind) {
-                           sapply(1:length(indices), 
-                                  function(z) {mean(y[indices[[z]]], na.rm = na.rm)})
+                           res <- sapply(1:length(indices), 
+                                  function(z) {
+                                  mean(y[indices[[z]]], na.rm = na.rm)})
+                           if (is.null(dim(res))) {
+                             dim(res) <- c(length(res))
+                           }
+                           return(res)
                          },
                          ind = indices,
                          output_dims = 'time',
@@ -26,8 +33,13 @@ agg_ini_end <- function(x, ini, end, indices = NULL, method, na.rm ,ncores) {
     x[[1]]$data <- Apply(x[[1]]$data,
                          target_dim = 'time',
                          function(y, ind) {
-                           sapply(1:length(indices), 
-                                  function(z) {sum(y[indices[[z]]], na.rm = na.rm)})
+                           res <- sapply(1:length(indices), 
+                                  function(z) {
+                                  sum(y[indices[[z]]], na.rm = na.rm)})
+                           if (is.null(dim(res))) {
+                             dim(res) <- c(length(res))
+                           }
+                           return(res)
                          },
                          ind = indices,
                          output_dims = 'time',