diff --git a/Scripts/compare_coef.R b/Scripts/compare_coef.R
new file mode 100644
index 0000000000000000000000000000000000000000..cd844db9832780b2e2f101e591acf378d038f2b6
--- /dev/null
+++ b/Scripts/compare_coef.R
@@ -0,0 +1,38 @@
+
+#load model
+model <- readRDS("Estimation_results/mxl/prediction/MXL_wtp_Prediction matching all cp_model.rds")
+
+compare_coefficients <- function(model, coef_name1, coef_name2) {
+  # Extract coefficients and standard errors
+  coef <- summary(model)$estimate
+  beta1 <- coef[coef_name1, "Estimate"]
+  beta2 <- coef[coef_name2, "Estimate"]
+  
+  
+  # Assuming model$robse contains robust standard errors
+  se1 <- model$se[coef_name1]
+  se2 <- model$se[coef_name2]
+  # Assuming model$robse contains robust standard errors
+  robse1 <- model$robse[coef_name1]
+  robse2 <- model$robse[coef_name2]
+  #Caluclate Difference
+  difference <- beta1-beta2
+  # Calculate Z-score
+  z_value <- (beta1 - beta2) / sqrt(se1^2 + se2^2)
+  
+  # Calculate p-value
+  p_value <- 2 * (1 - pnorm(abs(z_value)))
+  # Calculate Z-score
+  rob_z_value <-(beta1 - beta2) / sqrt(robse1^2 + robse2^2)
+  
+  # Calculate p-value
+  rob_p_value <- 2* (1 - pnorm(abs(rob_z_value)))
+  
+  # Return results as a list
+  return(list(difference = difference ,z_value = z_value, p_value = p_value ,rob_z_value = rob_z_value, rob_p_value = rob_p_value))
+}
+
+# Example usage
+compare_coefficients(model, "mu_nat_vol_treated", "mu_nat_treat_pred")
+compare_coefficients(model, "mu_nat_no_info", "mu_nat_treat_not_pred")
+