diff --git a/yellowstone_analysis/analysis.R b/yellowstone_analysis/analysis.R
index e0957c9..2903c57 100755
--- a/yellowstone_analysis/analysis.R
+++ b/yellowstone_analysis/analysis.R
@@ -21,6 +21,8 @@ df_raw = read_csv("bird_event_data_processed_vars.csv")
 # Set global variables #
 ########################
 
+run_analysis <- function(TIME_MAX){
+
 # Variables for the analysis
 TIME_MAX=240 # We take the yellowstone data from 0 - 180 secondes before / after passage of snowmobile
 
@@ -321,4 +323,70 @@ fig_trends = (boxplot_before | boxplot_after) +
   theme(axis.title.x = element_text(hjust=1.2)) 
 fig_trends
 
-ggsave("Figures_MS/Figure3.png" ,fig_trends ,dpi=300, width = 9, height = 7.5)
+ggsave(paste0("Figures_MS/Figure3-", TIME_MAX, ".png"),fig_trends ,dpi=300, width = 9, height = 7.5)
+
+return(list("model_before" = m_before, "model_after" = m_after))
+
+}
+
+#########################################################################
+########################## Sensitivity Analysis #########################
+#########################################################################
+time_min <- 120
+time_max <- 300
+step_size <- 20 
+
+# Initialize lists to store results
+all_results_before <- list()
+all_results_after <- list()
+
+# Run the function for various TIME_MAX values
+for (TIME_MAX in seq(time_min, time_max, by = step_size)) {
+  
+  results <- run_analysis(TIME_MAX)
+  
+  # Assuming 'results' contains two models, one for 'before' and one for 'after'
+  all_results_before[[as.character(TIME_MAX)]] <- summary(results$model_before)
+  all_results_after[[as.character(TIME_MAX)]] <- summary(results$model_after)
+  
+}
+
+# Extracting the coefficient for model_before & model_after
+get_coefficient <- function(results_list) {
+  sapply(results_list, function(res) coef(res)$cond[2]) # [variable_name, "Estimate"]
+}
+
+coeff_before <- get_coefficient(all_results_before)
+coeff_after <- get_coefficient(all_results_after)
+
+#########################
+# Correlation analysis #
+########################
+
+# Correlation
+cor_test <- cor.test(seq(time_min, time_max, by = step_size), coeff_before)
+# Linear Regression
+lm_fit_before <- lm(coeff_before ~ seq(time_min, time_max, by = step_size))
+rb <- summary(lm_fit_before)[4]$coefficients[2,4]
+
+# Correlation 
+cor_test <- cor.test(seq(time_min, time_max, by = step_size), coeff_after)
+
+# Linear Regression
+lm_fit_after <- lm(coeff_after ~ seq(time_min, time_max, by = step_size))
+ra <- summary(lm_fit_after)[4]$coefficients[2,4]
+
+#############################################
+# Visualization of the sensitivity analysis #
+#############################################
+x_coordinate = 270
+y_coordinate_a = 0.023
+y_coordinate_b = 0.020
+
+plot(seq(time_min, time_max, by = step_size), coeff_before, type = "l", 
+     xlab = "TIME_MAX", ylab = "Coefficient",
+     main = "Sensitivity Analysis for Coefficient of before and after passage of snowmobile")
+lines(seq(time_min, time_max, by = step_size), coeff_after, col = "blue")
+legend("topright", legend = c("closest_ss_before", "closest_ss_after"), fill = c("black", "blue"))
+text(x = x_coordinate, y = y_coordinate_a, labels = paste("p-value before: ", round(rb, digits = 7)), pos = 4, cex = 0.8)
+text(x = x_coordinate, y = y_coordinate_b, labels = paste("p-value after: ", round(ra, digits = 2)), pos = 4, cex = 0.8, col="blue")