diff --git a/Scripts/MAKE_FILE.R b/Scripts/MAKE_FILE.R
index 98189b39eaa6388ce2c279650f739801e10b2b5a..dd201e5a3475b9421e8c00a235ed162501446440 100644
--- a/Scripts/MAKE_FILE.R
+++ b/Scripts/MAKE_FILE.R
@@ -67,6 +67,7 @@ mxl_wtp_case_b <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case B")
 mxl_wtp_case_b_NR <- apollo_loadModel("Estimation_results/mxl/MXL_wtp NR B")
 mxl_wtp_case_c <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_Case_C")
 mxl_wtp_case_c_NR <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_NR_Case_C")
+mxl_wtp_case_d <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_Case_D base")
 
 # rent interactions models
 mxl_wtp_case_a_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case A Rent Int")
diff --git a/Scripts/create_tables.R b/Scripts/create_tables.R
index 9b14fa9f3a71fe1a76d02e0b12810bd5f1f2dbad..cceb70ddc0880cf762c5b7e9843e83433514ba4c 100644
--- a/Scripts/create_tables.R
+++ b/Scripts/create_tables.R
@@ -1,5 +1,5 @@
-#library(choiceTools, lib.loc = "/home/nc71qaxa/r-packages")
-library(choiceTools)
+library(choiceTools, lib.loc = "/home/nc71qaxa/r-packages")
+#library(choiceTools)
 
 dir.create("Tables/mxl")
 dir.create("Tables/logit")
@@ -131,7 +131,7 @@ texreg(l=list(logit_choice_treat_uni), stars = c(0.01, 0.05, 0.1), float.pos="tb
 ##### MXL #######
 
 ### Baseline case A
-case_A <- quicktexregapollo(mxl_wtp_case_a_rentINT)
+case_A <- quicktexregapollo(mxl_wtp_case_a)
 
 coef_names <- case_A@coef.names 
 coef_names <- sub("^(mu_)(.*)(_T|_VT)$", "\\2\\3", coef_names)
@@ -141,14 +141,14 @@ case_A@coef.names <- coef_names
 case_A_cols <- map(c("^mu_", "^sig_", "_T$", "_VT$"), subcoef, case_A)
 
 texreg(c(case_A_cols[1], remGOF(case_A_cols[2:4])),
-       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
-                              "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+       custom.coef.map = list("ASC_sq" = "ASC SQ", "natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
+                              "_natural" = "Naturalness", "nat" = "Naturalness",
                               "wd" = "Walking Distance", "asc" = "ASC SQ"),
        custom.model.names = c("Mean", "SD", "Treated", "Optional Treatment"), custom.note = "%stars (one-sided). Robust standard errors in parentheses.",
        stars = c(0.01, 0.05, 0.1), float.pos="tb", 
        label = "tab:mxl_A",
        caption = "Results of mixed logit model with treatment interactions for Case A.",
-       file="Tables/mxl/case_A_rent_INT.tex")
+       file="Tables/mxl/case_A.tex")
 
 ### Baseline case C
 case_C <- quicktexregapollo(mxl_wtp_case_c_rentINT)
@@ -184,8 +184,8 @@ case_C_NR@coef.names <- coef_names
 case_C_cols_NR <- map(c("^mu_", "^sig_", "_vid1$", "_vid2$", "_nv1$", "_nv2$", "_no_info$", "_NR$"), subcoef, case_C_NR)
 
 texreg(c(case_C_cols_NR[1], remGOF(case_C_cols_NR[2:8])),
-       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
-                              "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "ASC_sq" = "ASC SQ", "rent" = "Rent",
+                               "_natural" = "Naturalness", "nat" = "Naturalness",
                               "wd" = "Walking Distance", "asc" = "ASC SQ",
                               "ASC_sq_info" = "ASC SQ", "rent_info" = "Rent", "nat_info" = "Naturalness", "walking_info" = "Walking Distance"),
        custom.model.names = c("Mean", "SD", "Video 1", "Video 2", "Text 1", "Text 2", "No Info", "NR"), custom.note = "%stars. Robust standard errors in parentheses.",
@@ -196,7 +196,7 @@ texreg(c(case_C_cols_NR[1], remGOF(case_C_cols_NR[2:8])),
 
 
 ### New Case B
-case_B <- quicktexregapollo(mxl_wtp_case_d_rentINT)
+case_B <- quicktexregapollo(mxl_wtp_case_d)
 
 coef_names <- case_B@coef.names 
 coef_names <- sub("^(mu_)(.*)(vol_treated|_treated|_no_info)$", "\\2\\3", coef_names)
@@ -207,8 +207,8 @@ case_B@coef.names <- coef_names
 case_B_cols <- map(c("^mu_", "^sig_", "_treated$", "_vol_treated$","_no_info$"), subcoef, case_B)
 
 texreg(c(case_B_cols[1], remGOF(case_B_cols[2:5])),
-       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
-                              "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "ASC_sq" = "ASC SQ",  "rent" = "Rent",
+                              "_natural" = "Naturalness", "nat" = "Naturalness",
                               "wd" = "Walking Distance", "asc" = "ASC SQ",
                               "ASC_sq_info" = "ASC SQ", "rent_info" = "Rent", "nat_info" = "Naturalness", "walking_info" = "Walking Distance"),
        custom.model.names = c("Mean", "SD", "Treated", "Vol. Treated", "No Info"), custom.note = "%stars (one-sided). Robust standard errors in parentheses.",
@@ -251,3 +251,4 @@ texreg(c(case_B_cols_NR[1], remGOF(case_B_cols_NR[2:6])),
 #                               "mu_asc_T" = "ASC X Treated", "mu_nat_VT" = "Naturalness X Vol. Treated",  "mu_wd_VT" = "Walking Distance X Vol. Treated",  
 #                               "mu_rent_VT" = "Rent X Vol. Treated", "mu_asc_VT" = "ASC X Vol. Treated"), 
 #        stars = c(0.01, 0.05, 0.1), override.se = mxl_wtp_case_a_rentINT$robse, file="Tables/mxl/case_A_rent_INT.tex")
+
diff --git a/Scripts/data_prep.R b/Scripts/data_prep.R
index 5627a1cf170e9f75fbb5910c57f0da1aae6f28fe..90942782bd9f39bfe9cded11ab5e914f308ae3a8 100644
--- a/Scripts/data_prep.R
+++ b/Scripts/data_prep.R
@@ -1,134 +1,136 @@
-# Prepare variables that we want to use
-
-database_full <- database_full %>% rename(Gender = "Q03W123", Education = "Q06W123", HHSize = "Q41W123",
-                                  WorkingTime = "Q44W123", Birthyear = "Q01W123", Rent_net = "Q07W123",
-                                  Number_Kids = "Q42W123", Employment_type = "Q43W123", Conseq_UGS = "Q28W3",
-                                  Conseq_Money = "Q29W3") 
-
-
-database_full <- database_full %>% mutate(Gender = dplyr::recode(Gender, "A1" = 1, "A2" = 2, "A3"=3),
-                                Education = dplyr::recode(Education, "A1" = 1, "A2" = 2, "A3"=3, "A4" = 4, "A5" = 5),
-                                Employment_type = dplyr::recode(Employment_type, "A1" = 1, "A2" = 2, "A3"=3, "A4" = 4, 
-                                                         "A5" = 5, "A6" = 6),
-                                Conseq_UGS = dplyr::recode(Conseq_UGS, "A1" = 5, "A2" = 4, "A3"=3, "A4" = 2, "A5" = 1, "A6" = NA_real_),
-                                Conseq_Money = dplyr::recode(Conseq_Money, "A1" = 5, "A2" = 4, "A3"=3, "A4" = 2, "A5" = 1, "A6" = NA_real_))
-
-database_full <- database_full %>% mutate(Gender_female = case_when(Gender == 2 ~1, TRUE~0),
-                                          Age = 2023-Birthyear,
-                                          Uni_degree = case_when(Education == 5 ~1, TRUE~0),
-                                          Kids_Dummy = case_when(Number_Kids > 0 ~ 1, TRUE ~0),
-                                          Employment_full = case_when(Employment_type == 1 ~ 1, TRUE~0),
-                                          Pensioner = case_when(Employment_type == 6 ~ 1, TRUE~0),
-                                          Age_mean = Age - mean(Age),
-                                          Income_mean = QFIncome - mean(QFIncome))
-
-# Data cleaning 
-
-
-
-database_full <- database_full %>% filter(Rent_SQ <= 10000 & Rent_SQ >=50) %>%
-  filter(WalkingDistance_SQ > 0 & WalkingDistance_SQ <= 300) %>% 
-  filter(Gender!=3)
-
-database <- database %>% filter(Rent_SQ <= 10000 & Rent_SQ >=50) %>% 
-  filter(WalkingDistance_SQ > 0 & WalkingDistance_SQ <= 300)
-
-
-summary(database_full$interviewtime)
-
-database_full <- database_full %>% filter(interviewtime >= 300) # make change time to 10 seconds?
-
-
-database_full <- database_full %>%
-  filter(!is.na(Treatment_new)) %>%
-  mutate(Treatment_A = case_when(
-    Treatment == 1 ~ "Treated",
-    Treatment == 2 ~ "Vol_Treated",
-    Treatment == 3 ~ "Not_Treated",
-    TRUE ~ NA_character_
-  )) %>% 
-  mutate(Treatment_B = case_when(
-    Treatment_new == 1 | Treatment_new == 2 | Treatment_new == 4 | Treatment_new == 5 ~ "Treated",
-    Treatment_new == 3 | Treatment_new == 6 ~ "Not_Treated"
-  )) %>% 
-  mutate(Treatment_C = case_when(
-    Treatment_new == 1 ~ 'Video 1',
-    Treatment_new == 2 ~ 'No Video 1',
-    Treatment_new == 3 ~ 'No Info 2',
-    Treatment_new == 4 ~ 'No Video 2',
-    Treatment_new == 5 ~ 'Video 2',
-    Treatment_new == 6 ~ 'No Treatment 3',
-    TRUE ~ NA_character_
-  )) %>% 
-  mutate(Treatment_D = case_when(
-    Treatment_new == 1 | Treatment_new == 2 ~ 'Treated',
-    Treatment_new == 3 ~ 'No Info 2',
-    Treatment_new == 4 | Treatment_new == 5 ~ 'Vol. Treated',
-    Treatment_new == 6 ~ 'No Treatment 3',
-    TRUE ~ NA_character_
-  ))
-
-id_list <- unique(database_full$id)
-
-# Do we sill want to use this? or only database full?
-database <- database %>% filter(id %in% id_list) %>% filter(!is.na(Treatment_new))
-# Building NR Index
-
-for (i in 1:21) {
-  variable_name <- paste0("Q38S", sprintf("%02d", i), "W3")  # Generate variable name
-  cat("Table for", variable_name, ":\n")
-  print(table(database_full[[variable_name]]))
-  cat("\n")
-  database_full[[variable_name]] <- as.numeric(factor(database_full[[variable_name]], levels = c("A1", "A2", "A3", "A4", "A5")))
-  cat("Table for", variable_name, ":\n")
-  print(table(database_full[[variable_name]]))
-  cat("\n")
-}
-
-variables_to_reverse <- c("Q38S02W3", "Q38S03W3", "Q38S10W3", "Q38S11W3", "Q38S13W3", "Q38S14W3", "Q38S15W3", "Q38S18W3")
-for (variable_name in variables_to_reverse) {
-  cat("Table for", variable_name, ":\n")
-  
-  # Convert the variable to a factor with numerical levels and reverse the scores
-  database_full[[variable_name]] <- 6 - as.numeric(database_full[[variable_name]])
-  
-  # Print the table
-  print(table(database_full[[variable_name]]))
-  cat("\n")
-}
-q38_variables <- grep("^Q38", names(database_full), value = TRUE)
-database_full$Total_NR <- rowSums(database_full[q38_variables])
-hist(database_full$Total_NR)
-database_full <- database_full %>% 
-  mutate(Mean_NR=Total_NR/21) 
-mean_nr<-mean(database_full$Mean_NR, na.rm = TRUE)
-sd_nr<-sd(database_full$Mean_NR, na.rm = TRUE)
-database_full <- database_full %>% 
-  mutate(Z_Mean_NR=(Mean_NR-mean_nr)/sd_nr)
-database$Z_Mean_NR<- database_full$Z_Mean_NR
-summary(database$Z_Mean_NR)
-
-#Self-Reference Index
-
-for (i in 8:10) {
-  variable_name <- paste0("TV", sprintf("%02d", i), "W3")  # Generate variable name
-  cat("Table for", variable_name, ":\n")
-  print(table(database_full[[variable_name]]))
-  cat("\n")
-  database_full[[variable_name]] <- as.numeric(factor(database_full[[variable_name]], levels = c("A1", "A2", "A3", "A4", "A5")))
-  cat("Table for", variable_name, ":\n")
-  print(table(database_full[[variable_name]]))
-  cat("\n")
-}
-
-
-database_full$Total_SR <- database_full$TV08W3+database_full$TV09W3+database_full$TV10W3
-hist(database_full$Total_SR)
-database_full <- database_full %>% 
-  mutate(Mean_SR=Total_SR/3) 
-mean_sr<-mean(database_full$Mean_SR, na.rm = TRUE)
-sd_sr<-sd(database_full$Mean_SR, na.rm = TRUE)
-database_full <- database_full %>% 
-  mutate(Z_Mean_SR=(Mean_SR-mean_sr)/sd_sr)
-database$Z_Mean_SR<- database_full$Z_Mean_SR
-summary(database$Z_Mean_SR)
+# Prepare variables that we want to use
+
+database_full <- database_full %>% rename(Gender = "Q03W123", Education = "Q06W123", HHSize = "Q41W123",
+                                  WorkingTime = "Q44W123", Birthyear = "Q01W123", Rent_net = "Q07W123",
+                                  Number_Kids = "Q42W123", Employment_type = "Q43W123", Conseq_UGS = "Q28W3",
+                                  Conseq_Money = "Q29W3", UGS_visits = "Q23W123") 
+
+
+database_full <- database_full %>% mutate(Gender = dplyr::recode(Gender, "A1" = 1, "A2" = 2, "A3"=3),
+                                Education = dplyr::recode(Education, "A1" = 1, "A2" = 2, "A3"=3, "A4" = 4, "A5" = 5),
+                                Employment_type = dplyr::recode(Employment_type, "A1" = 1, "A2" = 2, "A3"=3, "A4" = 4, 
+                                                         "A5" = 5, "A6" = 6),
+                                Conseq_UGS = dplyr::recode(Conseq_UGS, "A1" = 5, "A2" = 4, "A3"=3, "A4" = 2, "A5" = 1, "A6" = NA_real_),
+                                Conseq_Money = dplyr::recode(Conseq_Money, "A1" = 5, "A2" = 4, "A3"=3, "A4" = 2, "A5" = 1, "A6" = NA_real_),
+                                UGS_visits = dplyr::recode(UGS_visits, "A1" = 1, "A2" = 2, "A3"=3, "A4" = 4, "A5" = 5, "A6" = 6,
+                                                           "A7" = 7, "A8" = 8, "A9" = 9, "A10"= 10))
+
+database_full <- database_full %>% mutate(Gender_female = case_when(Gender == 2 ~1, TRUE~0),
+                                          Age = 2023-Birthyear,
+                                          Uni_degree = case_when(Education == 5 ~1, TRUE~0),
+                                          Kids_Dummy = case_when(Number_Kids > 0 ~ 1, TRUE ~0),
+                                          Employment_full = case_when(Employment_type == 1 ~ 1, TRUE~0),
+                                          Pensioner = case_when(Employment_type == 6 ~ 1, TRUE~0),
+                                          Age_mean = Age - mean(Age),
+                                          Income_mean = QFIncome - mean(QFIncome))
+
+# Data cleaning 
+
+
+
+database_full <- database_full %>% filter(Rent_SQ <= 10000 & Rent_SQ >=50) %>%
+  filter(WalkingDistance_SQ > 0 & WalkingDistance_SQ <= 300) %>% 
+  filter(Gender!=3)
+
+database <- database %>% filter(Rent_SQ <= 10000 & Rent_SQ >=50) %>% 
+  filter(WalkingDistance_SQ > 0 & WalkingDistance_SQ <= 300)
+
+
+summary(database_full$interviewtime)
+
+database_full <- database_full %>% filter(interviewtime >= 300) # make change time to 10 seconds?
+
+
+database_full <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Treatment_A = case_when(
+    Treatment == 1 ~ "Treated",
+    Treatment == 2 ~ "Vol_Treated",
+    Treatment == 3 ~ "Not_Treated",
+    TRUE ~ NA_character_
+  )) %>% 
+  mutate(Treatment_B = case_when(
+    Treatment_new == 1 | Treatment_new == 2 | Treatment_new == 4 | Treatment_new == 5 ~ "Treated",
+    Treatment_new == 3 | Treatment_new == 6 ~ "Not_Treated"
+  )) %>% 
+  mutate(Treatment_C = case_when(
+    Treatment_new == 1 ~ 'Video 1',
+    Treatment_new == 2 ~ 'No Video 1',
+    Treatment_new == 3 ~ 'No Info 2',
+    Treatment_new == 4 ~ 'No Video 2',
+    Treatment_new == 5 ~ 'Video 2',
+    Treatment_new == 6 ~ 'No Treatment 3',
+    TRUE ~ NA_character_
+  )) %>% 
+  mutate(Treatment_D = case_when(
+    Treatment_new == 1 | Treatment_new == 2 ~ 'Treated',
+    Treatment_new == 3 ~ 'No Info 2',
+    Treatment_new == 4 | Treatment_new == 5 ~ 'Vol. Treated',
+    Treatment_new == 6 ~ 'No Treatment 3',
+    TRUE ~ NA_character_
+  ))
+
+id_list <- unique(database_full$id)
+
+# Do we sill want to use this? or only database full?
+database <- database %>% filter(id %in% id_list) %>% filter(!is.na(Treatment_new))
+# Building NR Index
+
+for (i in 1:21) {
+  variable_name <- paste0("Q38S", sprintf("%02d", i), "W3")  # Generate variable name
+  cat("Table for", variable_name, ":\n")
+  print(table(database_full[[variable_name]]))
+  cat("\n")
+  database_full[[variable_name]] <- as.numeric(factor(database_full[[variable_name]], levels = c("A1", "A2", "A3", "A4", "A5")))
+  cat("Table for", variable_name, ":\n")
+  print(table(database_full[[variable_name]]))
+  cat("\n")
+}
+
+variables_to_reverse <- c("Q38S02W3", "Q38S03W3", "Q38S10W3", "Q38S11W3", "Q38S13W3", "Q38S14W3", "Q38S15W3", "Q38S18W3")
+for (variable_name in variables_to_reverse) {
+  cat("Table for", variable_name, ":\n")
+  
+  # Convert the variable to a factor with numerical levels and reverse the scores
+  database_full[[variable_name]] <- 6 - as.numeric(database_full[[variable_name]])
+  
+  # Print the table
+  print(table(database_full[[variable_name]]))
+  cat("\n")
+}
+q38_variables <- grep("^Q38", names(database_full), value = TRUE)
+database_full$Total_NR <- rowSums(database_full[q38_variables])
+hist(database_full$Total_NR)
+database_full <- database_full %>% 
+  mutate(Mean_NR=Total_NR/21) 
+mean_nr<-mean(database_full$Mean_NR, na.rm = TRUE)
+sd_nr<-sd(database_full$Mean_NR, na.rm = TRUE)
+database_full <- database_full %>% 
+  mutate(Z_Mean_NR=(Mean_NR-mean_nr)/sd_nr)
+database$Z_Mean_NR<- database_full$Z_Mean_NR
+summary(database$Z_Mean_NR)
+
+#Self-Reference Index
+
+for (i in 8:10) {
+  variable_name <- paste0("TV", sprintf("%02d", i), "W3")  # Generate variable name
+  cat("Table for", variable_name, ":\n")
+  print(table(database_full[[variable_name]]))
+  cat("\n")
+  database_full[[variable_name]] <- as.numeric(factor(database_full[[variable_name]], levels = c("A1", "A2", "A3", "A4", "A5")))
+  cat("Table for", variable_name, ":\n")
+  print(table(database_full[[variable_name]]))
+  cat("\n")
+}
+
+
+database_full$Total_SR <- database_full$TV08W3+database_full$TV09W3+database_full$TV10W3
+hist(database_full$Total_SR)
+database_full <- database_full %>% 
+  mutate(Mean_SR=Total_SR/3) 
+mean_sr<-mean(database_full$Mean_SR, na.rm = TRUE)
+sd_sr<-sd(database_full$Mean_SR, na.rm = TRUE)
+database_full <- database_full %>% 
+  mutate(Z_Mean_SR=(Mean_SR-mean_sr)/sd_sr)
+database$Z_Mean_SR<- database_full$Z_Mean_SR
+summary(database$Z_Mean_SR)
diff --git a/Scripts/logit/chr_vol_treat.R b/Scripts/logit/chr_vol_treat.R
index 004568b0b44e6da61843e71e2f50fe268c96859c..6a805fe3d477d91193e8e1b2d441d60bc4809696 100644
--- a/Scripts/logit/chr_vol_treat.R
+++ b/Scripts/logit/chr_vol_treat.R
@@ -1,9 +1,13 @@
 library(dplyr)
 library(tidyr)
 library(margins)
+<<<<<<< HEAD
 library(lubridate)
 library(caret)
 library(randomForest)
+=======
+library(pROC)
+>>>>>>> e410a7a5c52ecf34b454cff74f0b641cd5615679
 
 library(xgboost)
 # Test treatment effect
@@ -14,7 +18,9 @@ database_full <- database_full %>%
          Dummy_Video_2 = case_when(Treatment_new == 5 ~ 1, TRUE ~ 0),
          Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0),
          Dummy_Info_nv1 = case_when(Treatment_new == 2 ~1, TRUE~0),
-         Dummy_Info_nv2 = case_when(Treatment_new == 4 ~1 , TRUE~0))
+         Dummy_Info_nv2 = case_when(Treatment_new == 4 ~1 , TRUE~0)) %>% 
+  mutate(Engagement_ugs = groupTime1731 + groupTime1726 + groupTime1769 + groupTime1727 + groupTime1770 + groupTime1771 + groupTime1768 +
+           groupTime1738)
 
 
 
@@ -39,9 +45,12 @@ summary(logit_choice_treat)
 
 
 logit_choice_treat_uni<-glm(Choice_Treat ~  as.factor(Gender)+Z_Mean_NR+Age_mean + QFIncome +
-                          Uni_degree + Kids_Dummy+ as.factor(Q27W123), data, family=binomial)
+                          Uni_degree + Kids_Dummy + Engagement_ugs + UGS_visits, data, family=binomial)
 summary(logit_choice_treat_uni)
 
+
+logit_sig <- glm(Choice_Treat ~  Z_Mean_NR + Age_mean + QFIncome, data, family=binomial)
+
 # Calculate marginal effects
 marginal_effects <- margins(logit_choice_treat_uni)
 
@@ -49,14 +58,16 @@ marginal_effects <- margins(logit_choice_treat_uni)
 summary(marginal_effects)               
 
 
+cut_off <- 0.45
 
 
 # Make predictions
-data$probabilities <- logit_choice_treat_uni %>% predict(data, type = "response")
-data$predicted.classes <- ifelse(data$probabilities  >= 0.5, 1, 0)
+data$probabilities <- logit_sig %>% predict(data, type = "response")
+data$predicted.classes <- ifelse(data$probabilities  >= cut_off, 1, 0)
 # Model accuracy
 mean(data$predicted.classes == data$Choice_Treat, na.rm = T)
 
+<<<<<<< HEAD
 
 
 
@@ -134,3 +145,38 @@ table(labeled_data$Choice_Treat, labeled_data$PredictedGroup)
 unlabeled_predictions <- predict(model3, newdata = unlabeled_data)
 unlabeled_data$PredictedGroup <- unlabeled_predictions
 print(model3$bestTune)
+=======
+table(data$predicted.classes ,data$Choice_Treat)
+
+calculate_metrics <- function(confusion_matrix) {
+  # Extract values from the confusion matrix
+  TN <- confusion_matrix[1, 1]
+  FP <- confusion_matrix[2, 1]
+  FN <- confusion_matrix[1, 2]
+  TP <- confusion_matrix[2, 2]
+  
+  # Calculate sensitivity
+  sensitivity <- round(TP / (TP + FN), 2)
+  
+  # Calculate specificity
+  specificity <- round(TN / (TN + FP), 2)
+  
+  # Return the results as a list
+  return(list(sensitivity = sensitivity, specificity = specificity))
+}
+
+# Calculate metrics
+metrics <- calculate_metrics(table(data$predicted.classes ,data$Choice_Treat))
+
+# Print results
+print(paste("Sensitivity:", metrics$sensitivity))
+print(paste("Specificity:", metrics$specificity))
+roc_obj <- roc(data$Choice_Treat, data$probabilities)
+plot(roc_obj, main = "ROC Curve", col = "blue", lwd = 2)
+auc_value <- auc(roc_obj)
+best_coords <- coords(roc_obj, "best", best.method="youden")
+
+
+cut_off <- best_coords$threshold
+
+>>>>>>> e410a7a5c52ecf34b454cff74f0b641cd5615679
diff --git a/Scripts/mxl/Prediction models/mxl_wtp_space_caseE.R b/Scripts/mxl/Prediction models/mxl_wtp_space_caseE.R
new file mode 100644
index 0000000000000000000000000000000000000000..d7b40621d7bc35b4bbd35c5ecf815104de135ca7
--- /dev/null
+++ b/Scripts/mxl/Prediction models/mxl_wtp_space_caseE.R	
@@ -0,0 +1,175 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database_full$probabilities <- logit_sig %>% predict(database_full, type = "response")
+database_full$predicted.classes <- ifelse(database_full$probabilities  >= cut_off, 1, 0)
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0),
+         Dummy_pred_Treated = case_when(Treatment_new == 6 & predicted.classes == 1 ~ 1, TRUE~0)) 
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_E base",
+  modelDescr = "MXL wtp space Case E base",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_ASC_sq_pred_treat = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_pred_treat = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_pred_treat = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
+                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated - Rent_1
+                            +  mu_nat_pred_treat * Naturalness_1 * Dummy_pred_Treated + mu_walking_pred_treat * WalkingDistance_1 *Dummy_pred_Treated)
+  
+  V[['alt2']] = -(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
+                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated- Rent_2
+                            + mu_nat_pred_treat * Naturalness_2 * Dummy_pred_Treated + mu_walking_pred_treat * WalkingDistance_2 *Dummy_pred_Treated)
+  
+  V[['alt3']] = -(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
+                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                            +  mu_ASC_sq_no_info * Dummy_no_info - Rent_3
+                            +  mu_nat_pred_treat * Naturalness_3 * Dummy_pred_Treated + mu_walking_pred_treat * WalkingDistance_3 *Dummy_pred_Treated
+                            +  mu_ASC_sq_pred_treat * Dummy_pred_Treated)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_e = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_e)
+
+
diff --git a/Scripts/mxl/Prediction models/mxl_wtp_space_caseF.R b/Scripts/mxl/Prediction models/mxl_wtp_space_caseF.R
new file mode 100644
index 0000000000000000000000000000000000000000..1bc9c4e08f19e756f7091baa5b973ac747c6a2fe
--- /dev/null
+++ b/Scripts/mxl/Prediction models/mxl_wtp_space_caseF.R	
@@ -0,0 +1,183 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database_full$probabilities <- logit_sig %>% predict(database_full, type = "response")
+database_full$predicted.classes <- ifelse(database_full$probabilities  >= cut_off, 1, 0)
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0),
+         Dummy_pred_Treated = case_when(Treatment_new == 6 & predicted.classes == 1 ~ 1, TRUE~0),
+         Dummy_Treated_predicted = case_when(Treatment_new == 1  & predicted.classes == 1|Treatment_new == 2 & predicted.classes == 1 ~ 1, TRUE ~ 0),
+         Dummy_Treated_notpredicted = case_when(Treatment_new == 1  & predicted.classes == 0|Treatment_new == 2 & predicted.classes == 0 ~ 1, TRUE ~ 0)) 
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_F base",
+  modelDescr = "MXL wtp space Case F base",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated_predicted = 0,
+              mu_ASC_sq_treated_notpredicted = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_ASC_sq_pred_treat = 0,
+              mu_nat_treated_predicted =0,
+              mu_nat_treated_notpredicted = 0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_pred_treat = 0,
+              mu_walking_treated_predicted =0,
+              mu_walking_treated_notpredicted =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_pred_treat = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                +  mu_nat_treated_predicted * Naturalness_1 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                              +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_1 *Dummy_Treated_notpredicted
+                              +  mu_walking_treated_predicted * WalkingDistance_1 *Dummy_Treated_predicted +  mu_walking_treated_notpredicted * WalkingDistance_1 *Dummy_Treated_notpredicted
+                              + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                              +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated - Rent_1
+                              +  mu_nat_pred_treat * Naturalness_1 * Dummy_pred_Treated + mu_walking_pred_treat * WalkingDistance_1 * Dummy_pred_Treated)
+  
+  V[['alt2']] = -(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                              + mu_nat_treated_predicted * Naturalness_2 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                              + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_2 *Dummy_Treated_notpredicted
+                              + mu_walking_treated_predicted * WalkingDistance_2 *Dummy_Treated_predicted + mu_walking_treated_notpredicted * WalkingDistance_2 *Dummy_Treated_notpredicted
+                              + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                              + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated- Rent_2
+                              + mu_nat_pred_treat * Naturalness_2 * Dummy_pred_Treated + mu_walking_pred_treat * WalkingDistance_2 * Dummy_pred_Treated)
+  
+  V[['alt3']] = -(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                              +  mu_nat_treated_predicted * Naturalness_3 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                              +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_3 *Dummy_Treated_notpredicted
+                              +  mu_walking_treated_predicted * WalkingDistance_3 *Dummy_Treated_predicted + mu_walking_treated_notpredicted * WalkingDistance_3 *Dummy_Treated_notpredicted
+                              +  mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                              +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                              +  mu_ASC_sq_treated_predicted * Dummy_Treated_predicted + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                              +  mu_ASC_sq_no_info * Dummy_no_info - Rent_3 + mu_ASC_sq_treated_notpredicted * Dummy_Treated_notpredicted
+                              +  mu_nat_pred_treat * Naturalness_3 * Dummy_pred_Treated + mu_walking_pred_treat * WalkingDistance_3 * Dummy_pred_Treated
+                              +  mu_ASC_sq_pred_treat * Dummy_pred_Treated)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_f = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_f)
+
+
diff --git a/Scripts/mxl/Prediction models/mxl_wtp_space_caseG.R b/Scripts/mxl/Prediction models/mxl_wtp_space_caseG.R
new file mode 100644
index 0000000000000000000000000000000000000000..7a8271185cc81508c4d0856d815959684fba0e3d
--- /dev/null
+++ b/Scripts/mxl/Prediction models/mxl_wtp_space_caseG.R	
@@ -0,0 +1,176 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database_full$probabilities <- logit_sig %>% predict(database_full, type = "response")
+database_full$predicted.classes <- ifelse(database_full$probabilities  >= cut_off, 1, 0)
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0),
+         Dummy_pred_Treated = case_when(Treatment_new == 6 & predicted.classes == 1 ~ 1, TRUE~0),
+         Dummy_Treated_predicted = case_when(Treatment_new == 1  & predicted.classes == 1|Treatment_new == 2 & predicted.classes == 1 ~ 1, TRUE ~ 0),
+         Dummy_Treated_notpredicted = case_when(Treatment_new == 1  & predicted.classes == 0|Treatment_new == 2 & predicted.classes == 0 ~ 1, TRUE ~ 0)) 
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_G base",
+  modelDescr = "MXL wtp space Case G base",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated_predicted = 0,
+              mu_ASC_sq_treated_notpredicted = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_nat_treated_predicted =0,
+              mu_nat_treated_notpredicted = 0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_walking_treated_predicted =0,
+              mu_walking_treated_notpredicted =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                              +  mu_nat_treated_predicted * Naturalness_1 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                              +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_1 *Dummy_Treated_notpredicted
+                              +  mu_walking_treated_predicted * WalkingDistance_1 *Dummy_Treated_predicted +  mu_walking_treated_notpredicted * WalkingDistance_1 *Dummy_Treated_notpredicted
+                              +  mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                              +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                              + mu_nat_treated_predicted * Naturalness_2 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                              + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_2 *Dummy_Treated_notpredicted
+                              + mu_walking_treated_predicted * WalkingDistance_2 *Dummy_Treated_predicted + mu_walking_treated_notpredicted * WalkingDistance_2 *Dummy_Treated_notpredicted
+                              + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                              + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated- Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                              +  mu_nat_treated_predicted * Naturalness_3 *Dummy_Treated_predicted + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                              +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated + mu_nat_treated_notpredicted * Naturalness_3 *Dummy_Treated_notpredicted
+                              +  mu_walking_treated_predicted * WalkingDistance_3 *Dummy_Treated_predicted + mu_walking_treated_notpredicted * WalkingDistance_3 *Dummy_Treated_notpredicted
+                              +  mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                              +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                              +  mu_ASC_sq_treated_predicted * Dummy_Treated_predicted + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                              +  mu_ASC_sq_no_info * Dummy_no_info - Rent_3 + mu_ASC_sq_treated_notpredicted * Dummy_Treated_notpredicted)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_g = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_g)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_NR_caseA_rentINT.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_NR_caseA_rentINT.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_NR_caseA_rentINT.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_NR_caseA_rentINT.R
diff --git a/Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_NR_caseC_RentINT.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_NR_caseC_RentINT.R
diff --git a/Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT_X.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_NR_caseC_RentINT_X.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_NR_caseC_RentINT_X.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_NR_caseC_RentINT_X.R
diff --git a/Scripts/mxl/mxl_wtp_space_caseA_rentINT.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseA_rentINT.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_caseA_rentINT.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseA_rentINT.R
diff --git a/Scripts/mxl/mxl_wtp_space_caseB_rentINT.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseB_rentINT.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_caseB_rentINT.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseB_rentINT.R
diff --git a/Scripts/mxl/mxl_wtp_space_caseC_rentINT.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseC_rentINT.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_caseC_rentINT.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseC_rentINT.R
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_rentINT.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseD_RentINT.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_caseD_rentINT.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseD_RentINT.R
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_RentINT_NR.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseD_RentINT_NR.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_caseD_RentINT_NR.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseD_RentINT_NR.R
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R b/Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseD_RentINT_X.R
similarity index 100%
rename from Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R
rename to Scripts/mxl/Rent Interaction models/mxl_wtp_space_caseD_RentINT_X.R
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_lognormRent_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_lognormRent_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..a992b8ee23fc6628062a9e8ecf375279038dab9d
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_lognormRent_scale.R
@@ -0,0 +1,153 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A lognorm Rent Scale",
+    modelDescr = "MXL_pref Case A lognorm Rent Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT = 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_rent_T * Rent_1 * Dummy_Treated +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated +
+                                mu_rent_VT * Rent_1 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                               mu_rent_T * Rent_2 * Dummy_Treated +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                                mu_rent_VT * Rent_2 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_rent_T * Rent_3 * Dummy_Treated + mu_asc_VT * Dummy_Vol_Treated +
+                                mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated +
+                                mu_rent_VT * Rent_3 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_lognorm_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_lognorm_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_normRent_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_normRent_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..ed63cc7ec293593043abd7d0936fe9e7bbd18ce4
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_normRent_scale.R
@@ -0,0 +1,153 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A norm Rent Scale",
+    modelDescr = "MXL_pref Case A norm Rent Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT = 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = mu_rent + sig_rent * draws_rent
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_rent_T * Rent_1 * Dummy_Treated +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated +
+                                mu_rent_VT * Rent_1 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                               mu_rent_T * Rent_2 * Dummy_Treated +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                                mu_rent_VT * Rent_2 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_rent_T * Rent_3 * Dummy_Treated + mu_asc_VT * Dummy_Vol_Treated +
+                                mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated +
+                                mu_rent_VT * Rent_3 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_norm_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_norm_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_only_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_only_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..477de7cb5c435eebd92ea763ba355b80321faa77
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_only_scale.R
@@ -0,0 +1,136 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A only  Scale",
+    modelDescr = "MXL_pref Case A only  Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + 
+                                                                                   b_mu_walking * WalkingDistance_1) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + 
+                                                                                   b_mu_walking * WalkingDistance_2) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_mu_natural * Naturalness_3 + 
+                                                                                   b_mu_walking * WalkingDistance_3 + b_ASC_sq) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_only_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_only_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..16c49729bedc1b761810dcaa196bf0acc091e6a2
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_scale.R
@@ -0,0 +1,149 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A  Scale",
+    modelDescr = "MXL_pref Case A  Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT = 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated +
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                                mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT * Dummy_Vol_Treated +
+                                mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated +
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseD_only_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseD_only_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..9db5aae3ad5c92bc134d21e60ac3ea20985a72cd
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseD_only_scale.R
@@ -0,0 +1,150 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D scale",
+  modelDescr = "MXL wtp space Case D scale",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2,
+              lambda_T = 0,
+              lambda_VT =0,
+              lambda_NO = 0)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated + lambda_NO * Dummy_no_info)*
+                            (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                             b_mu_rent * Rent_1)
+  
+  V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated + lambda_NO * Dummy_no_info)*
+                            (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                             b_mu_rent * Rent_2)
+  
+  V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated + lambda_NO * Dummy_no_info)*
+                            (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 +
+                             b_mu_rent * Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_d_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_d_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_base_scale.R b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_base_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..c6ded1ef21d2e6bde1d10e7f1452e14d1575f082
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_base_scale.R
@@ -0,0 +1,149 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Case A Scale",
+    modelDescr = "MXL_wtp Case A Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT= 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*-(b_mu_rent)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated+
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated
+                              - Rent_1)
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*-(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                              mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated
+                              - Rent_2)
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*-(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT *  Dummy_Vol_Treated +mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated + 
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated
+                                - Rent_3)
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_case_a_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_case_a_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_exp_scale.R b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_exp_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..74574cb8652e858085e980d6cec7823731c889c9
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_exp_scale.R
@@ -0,0 +1,151 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Case A Rent Int Exp Scale",
+    modelDescr = "MXL_wtp Case A Rent Int Exp Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT= 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + exp(mu_rent_T)*Dummy_Treated + exp(mu_rent_VT)*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated+
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated
+                              - Rent_1))
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + exp(mu_rent_T)*Dummy_Treated + exp(mu_rent_VT)*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                              mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated
+                              - Rent_2))
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + exp(mu_rent_T)*Dummy_Treated + exp(mu_rent_VT)*Dummy_Vol_Treated)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT *  Dummy_Vol_Treated +mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated + 
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated
+                                - Rent_3))
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_case_a_rentINT_exp_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_case_a_rentINT_exp_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_scale.R b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..e11aaf4cfaa2ddbb4f4031005f939a9cf78e6ba4
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_scale.R
@@ -0,0 +1,151 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Case A Rent Int Scale",
+    modelDescr = "MXL_wtp Case A Rent Int Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT= 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + mu_rent_T*Dummy_Treated + mu_rent_VT*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated+
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated
+                              - Rent_1))
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + mu_rent_T*Dummy_Treated + mu_rent_VT*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                              mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated
+                              - Rent_2))
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + mu_rent_T*Dummy_Treated + mu_rent_VT*Dummy_Vol_Treated)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT *  Dummy_Vol_Treated +mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated + 
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated
+                                - Rent_3))
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_case_a_rentINT_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_case_a_rentINT_scale)
+
+
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A_NR.R b/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A_NR.R
new file mode 100644
index 0000000000000000000000000000000000000000..b96c49bd26358c82387df6668f49398964e49a9a
--- /dev/null
+++ b/Scripts/mxl/Split_samples/mxl_wtp_space_not_tr_A_NR.R
@@ -0,0 +1,148 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Video_1 = case_when(Treatment_new == 1 ~ 1, TRUE ~ 0),
+         Dummy_Video_2 = case_when(Treatment_new == 5 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0),
+         Dummy_Info_nv1 = case_when(Treatment_new == 2 ~1, TRUE~0),
+         Dummy_Info_nv2 = case_when(Treatment_new == 4 ~1 , TRUE~0)) %>% 
+  filter(Treatment_A == "Not_Treated")
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Not_Treated A NR",
+    modelDescr = "MXL wtp space Not_Treated A NR",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/Split_samples"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                mu_nat_NR = 0,
+                mu_walking_NR = 0,
+                mu_asc_NR = 0,
+                ASC_sq = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                              mu_nat_NR * Naturalness_1 * Z_Mean_NR + mu_walking_NR * WalkingDistance_1 * Z_Mean_NR -
+                              Rent_1)
+    
+    V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                              mu_nat_NR * Naturalness_2 * Z_Mean_NR + mu_walking_NR * WalkingDistance_2 * Z_Mean_NR -
+                              Rent_2)
+    
+    V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 +
+                              mu_nat_NR * Naturalness_3 * Z_Mean_NR + mu_walking_NR * WalkingDistance_3 * Z_Mean_NR -
+                              Rent_3 + mu_asc_NR * Z_Mean_NR)
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_not_tr_a_nr = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_not_tr_a_nr)
+
+
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A_NR.R b/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A_NR.R
new file mode 100644
index 0000000000000000000000000000000000000000..28870ec5b24710f90a3cbd184534cd9b597873dc
--- /dev/null
+++ b/Scripts/mxl/Split_samples/mxl_wtp_space_tr_A_NR.R
@@ -0,0 +1,149 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Video_1 = case_when(Treatment_new == 1 ~ 1, TRUE ~ 0),
+         Dummy_Video_2 = case_when(Treatment_new == 5 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0),
+         Dummy_Info_nv1 = case_when(Treatment_new == 2 ~1, TRUE~0),
+         Dummy_Info_nv2 = case_when(Treatment_new == 4 ~1 , TRUE~0)) %>% 
+  filter(Treatment_A == "Treated")
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Treated A NR",
+    modelDescr = "MXL wtp space Treated A NR",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/Split_samples"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_nat_NR = 0,
+                mu_walking_NR = 0,
+                mu_asc_NR = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_NR * Naturalness_1 * Z_Mean_NR + mu_walking_NR * WalkingDistance_1 * Z_Mean_NR -
+                                Rent_1)
+    
+    V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_NR * Naturalness_2 * Z_Mean_NR + mu_walking_NR * WalkingDistance_2 * Z_Mean_NR -
+                                Rent_2)
+    
+    V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 +
+                                mu_nat_NR * Naturalness_3 * Z_Mean_NR + mu_walking_NR * WalkingDistance_3 * Z_Mean_NR -
+                                Rent_3 + mu_asc_NR * Z_Mean_NR)
+    
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_tr_a_nr = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_tr_a_nr)
+
+
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat.R b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat.R
new file mode 100644
index 0000000000000000000000000000000000000000..8177fbe1d15e9860278bf60372ff976de482cb67
--- /dev/null
+++ b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat.R
@@ -0,0 +1,136 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(Treatment_D == "Vol. Treated")
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Vol_Treat",
+    modelDescr = "MXL wtp space Vol_Treat",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/Split_samples"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 -
+                              Rent_1)
+    
+    V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 -
+                              Rent_2)
+    
+    V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 -
+                              Rent_3)
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_Vol_Treat = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_Vol_Treat)
+
+
diff --git a/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat_NR.R b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat_NR.R
new file mode 100644
index 0000000000000000000000000000000000000000..dfac0a3161ac0c63b5ced965063899dd063c7d1b
--- /dev/null
+++ b/Scripts/mxl/Split_samples/mxl_wtp_space_vol_treat_NR.R
@@ -0,0 +1,143 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  filter(Treatment_D == "Vol. Treated")
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Vol_Treat NR",
+    modelDescr = "MXL wtp space Vol_Treat NR",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/Split_samples"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                mu_nat_NR = 0,
+                mu_walking_NR = 0,
+                mu_asc_NR = 0,
+                ASC_sq = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = -b_mu_rent*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                              mu_nat_NR * Naturalness_1 * Z_Mean_NR + mu_walking_NR * WalkingDistance_1 * Z_Mean_NR -
+                              Rent_1)
+    
+    V[['alt2']] = -b_mu_rent*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                              mu_nat_NR * Naturalness_2 * Z_Mean_NR + mu_walking_NR * WalkingDistance_2 * Z_Mean_NR -
+                              Rent_2)
+    
+    V[['alt3']] = -b_mu_rent*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 +
+                              mu_nat_NR * Naturalness_3 * Z_Mean_NR + mu_walking_NR * WalkingDistance_3 * Z_Mean_NR -
+                              Rent_3 + mu_asc_NR * Z_Mean_NR)
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_Vol_Treat_NR = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_Vol_Treat_NR)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_caseD.R b/Scripts/mxl/mxl_wtp_space_caseD.R
new file mode 100644
index 0000000000000000000000000000000000000000..d51b52d3acd20fe9293fb056b77e23ef123b495c
--- /dev/null
+++ b/Scripts/mxl/mxl_wtp_space_caseD.R
@@ -0,0 +1,164 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D base",
+  modelDescr = "MXL wtp space Case D base",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
+                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
+                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated- Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
+                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                            +  mu_ASC_sq_no_info * Dummy_no_info - Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_c = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_c)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_NR.R b/Scripts/mxl/mxl_wtp_space_caseD_NR.R
new file mode 100644
index 0000000000000000000000000000000000000000..f14a8bfda1014b46ebde6859435c3e8b77fd33cf
--- /dev/null
+++ b/Scripts/mxl/mxl_wtp_space_caseD_NR.R
@@ -0,0 +1,180 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D_NR base",
+  modelDescr = "MXL wtp space Case D NR base",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_ASC_NR = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_NR = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_NR = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*
+                            (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
+                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_1 +
+                            + mu_walking_NR * Z_Mean_NR * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent)*                          
+                            (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
+                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_2 + 
+                            + mu_walking_NR * Z_Mean_NR * WalkingDistance_2
+                            - Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent)*
+                            (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
+                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                            +  mu_ASC_sq_no_info * Dummy_no_info 
+                            +  mu_ASC_NR * Z_Mean_NR 
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_3 
+                            +  mu_walking_NR * Z_Mean_NR * WalkingDistance_3
+                            - Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_d_NR = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_d_NR)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_age.R b/Scripts/mxl/mxl_wtp_space_caseD_age.R
new file mode 100644
index 0000000000000000000000000000000000000000..9409124d1047a9d7704cae84807a9e278d289e92
--- /dev/null
+++ b/Scripts/mxl/mxl_wtp_space_caseD_age.R
@@ -0,0 +1,180 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D_Age",
+  modelDescr = "MXL wtp space Case D Age",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_asc_age = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_age = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_age = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*
+                            (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
+                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated
+                            +  mu_nat_age * Age_mean *Naturalness_1 +
+                            + mu_walking_age * Age_mean * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent)*                          
+                            (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
+                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated
+                            +  mu_nat_age * Age_mean *Naturalness_2 + 
+                            + mu_walking_age * Age_mean * WalkingDistance_2
+                            - Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent)*
+                            (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
+                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                            +  mu_ASC_sq_no_info * Dummy_no_info 
+                            +  mu_asc_age * Age_mean
+                            +  mu_nat_age * Age_mean *Naturalness_3 
+                            +  mu_walking_age * Age_mean * WalkingDistance_3
+                            - Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_d_age = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_d_age)
+
+
diff --git a/project_start.qmd b/project_start.qmd
index f601ceff805136ce30d01764911ae85032e5caef..aff9ad0ad4b3edeb16090aea063c165975f791fb 100644
--- a/project_start.qmd
+++ b/project_start.qmd
@@ -294,9 +294,15 @@ ggpubr::ggarrange(plot_mani_A, plot_cons_A)
 
 
 
+<<<<<<< HEAD
 ```{=tex}
 \begin{equation}
     U_i = -(\beta_{C_i} + \beta_{TreatC_i} \cdot v_{Treat}) \cdot (\beta_{X_i} \cdot v_{X_i} + \beta_{TreatX_i} \cdot v_{X_i} \cdot v_{Treat} - C_i) + \epsilon_i
+=======
+```{=latex}
+\begin{equation}
+    U_i = -\beta_{C_i} \cdot (\beta_{X_i} \cdot v_{X_i} + \beta_{TreatX_i} \cdot v_{X_i} \cdot v_{Treat} - C_i) + \epsilon_i
+>>>>>>> e410a7a5c52ecf34b454cff74f0b641cd5615679
     \label{mxl_base}
 \end{equation}
 ```
@@ -324,8 +330,13 @@ with
 ::: {style="font-size: 80%;"}
 ```{r, results='asis'}
 htmlreg(c(case_A_cols[1], remGOF(case_A_cols[2:4])), single.row = TRUE,
+<<<<<<< HEAD
        custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
                               "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+=======
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "ASC_sq" = "ASC SQ",  "rent" = "Rent",
+                              "_natural" = "Naturalness", "nat" = "Naturalness",
+>>>>>>> e410a7a5c52ecf34b454cff74f0b641cd5615679
                               "wd" = "Walking Distance", "asc" = "ASC SQ"),
        custom.model.names = c("Mean", "SD", "Treated", "Optional Treatment"), custom.note = "%stars. Standard errors in parentheses.",
        stars = c(0.01, 0.05, 0.1), float.pos="tb", 
@@ -395,8 +406,13 @@ ggpubr::ggarrange(plot_mani_D, plot_cons_D)
 ::: {style="font-size: 80%;"}
 ```{r, results='asis'}
 htmlreg(c(case_B_cols[1], remGOF(case_B_cols[2:5])),
+<<<<<<< HEAD
        custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
                               "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+=======
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "ASC_sq" = "ASC SQ", "rent" = "Rent",
+                              "_natural" = "Naturalness", "nat" = "Naturalness",
+>>>>>>> e410a7a5c52ecf34b454cff74f0b641cd5615679
                               "wd" = "Walking Distance", "asc" = "ASC SQ",
                               "ASC_sq_info" = "ASC SQ", "rent_info" = "Rent", "nat_info" = "Naturalness", "walking_info" = "Walking Distance"),
        custom.model.names = c("Mean", "SD", "Treated", "Vol. Treated", "No Info"), custom.note = "%stars (one-sided). Robust standard errors in parentheses.",
@@ -411,8 +427,13 @@ htmlreg(c(case_B_cols[1], remGOF(case_B_cols[2:5])),
 ::: {style="font-size: 80%;"}
 ```{r, results='asis'}
 htmlreg(c(case_B_cols_NR[1], remGOF(case_B_cols_NR[2:6])),
+<<<<<<< HEAD
        custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
                               "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+=======
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "ASC_sq" = "ASC SQ", "rent" = "Rent",
+                              "_natural" = "Naturalness", "nat" = "Naturalness",
+>>>>>>> e410a7a5c52ecf34b454cff74f0b641cd5615679
                               "wd" = "Walking Distance", "asc" = "ASC SQ", "ASC" ="ASC SQ",
                               "ASC_sq_info" = "ASC SQ", "rent_info" = "Rent", "nat_info" = "Naturalness", "walking_info" = "Walking Distance"),
        custom.model.names = c("Mean", "SD", "Treated", "Vol. Treated", "No Info", "NR-Index"), custom.note = "%stars (one-sided). Robust standard errors in parentheses.",
diff --git a/project_start_wonv.qmd b/project_start_wonv.qmd
new file mode 100644
index 0000000000000000000000000000000000000000..aebd3fa49f4b97931a54f8dbd2cebbb0758958a1
--- /dev/null
+++ b/project_start_wonv.qmd
@@ -0,0 +1,598 @@
+---
+title: "Hot Cities, Cool Choices:" 
+subtitle: "The Effect of Optional and Obligatory Information on Stated Preferences for Urban Green Spaces"
+title-slide-attributes:
+  data-background-image: Grafics/iDiv_logo_item.png
+  data-background-size: contain
+  data-background-opacity: "0.2"
+author: "Nino Cavallaro, Fabian Marder, Julian Sagebiel"
+institute: 
+- German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig
+- Leipzig University 
+date: today
+date-format: long
+bibliography: references.bib
+filters:
+  - parse-latex
+  - custom_app.lua
+format: 
+  revealjs:
+    slide-number: true
+    smaller: true
+    logo: Grafics/iDiv_logo_item.PNG
+    footer: "WONV 2024: Hot Cities, Cool Choices"
+    scrollable: true
+    embed-resources: true
+---
+
+```{r, include=FALSE, cache=FALSE}
+source("Scripts/MAKE_FILE.R")
+```
+
+```{r loadlibs, include=FALSE}
+library(tidyverse)
+library(apollo)
+library(texreg)
+
+list_ols <- list("(Intercept)" = "Intercept", "as.factor(Treatment_A)Treated" = "Treated", "as.factor(Treatment_A)Vol_Treated" = "Optional Treatment",
+                 "as.factor(Treatment_C)No Info 2" = "No Info 2", "as.factor(Treatment_C)No Video 1" = "Text 1",
+                 "as.factor(Treatment_C)No Video 2" = "Text 2", "as.factor(Treatment_C)Video 1" = "Video 1",
+                 "as.factor(Treatment_C)Video 2" = "Video 2", "as.factor(Treatment_D)Treated" = "Treated", "as.factor(Treatment_D)Vol. Treated" = "Vol. Treated",
+                 "as.factor(Treatment_D)No Info 2" = "No Info", 
+                 "Z_Mean_NR" = "NR-Index", "as.factor(Gender)2" = "Female",
+                 "Age_mean" = "Age", "QFIncome" = "Income", "Uni_degree" = "University Degree")
+```
+
+# Motivation & Research Contribution
+
+## Motivation (1)
+
+::: incremental
+-   **Stated preference** methods are frequently applied in **environmental valuation** to estimate economic values of policies, goods, and services that cannot be valued otherwise
+-   Stated preference methods face **validity challenges**
+-   Valid value estimation requires **sufficient information** provision about the good being valued
+-   Still unclear **what formats of information provision** and **how much information** are optimal for valid preference elicitation
+:::
+
+## Motivation (2)
+
+::: incremental
+-   Too **much information** may increase survey **complexity**, leading to respondents being overburdened with it and producing less consistent choices
+-   Too **little information** may lead respondents to **not** being able to make an **informed choice**
+-   Valid preference elicitation depends not only on the provision of information, but also on the **appropriate processing and recall** of the information by the respondent
+-   **Optional information** allows the respondents to gather required information if needed and might increase efficiency of information provision
+
+-  Providing optional information should enhance optimal information seeking leading to less heterogeneity in good-specific knowledge between the respondents 
+:::
+
+## Literature
+
+::: incremental
+-   There is **little research** on the effects of **optional information provision** on choice behavior and information recall
+-   In their study, @tienhaara2022information surveyed preferences for agricultural genetic resources, allowing respondents the option to access detailed information on the valued goods prior to preference elicitation
+-   Similarly, @hu2009consumers offered respondents the opportunity to access optional information about genetic modified food before participating in a choice experiment
+-   Both studies conclude that, on average, respondents who voluntary retrieve information have **larger willingness to pay** for the good to be valued
+-   Their study design, however, does not allow comparing the optional information retrieval to a version where the additional information was shown obligatory
+:::
+
+## Research Contribution
+
+::: incremental
+-   Our study explores the impact of additional obligatory and optional information on stated preferences using an exogenous split sample approach with three treatments
+-   We investigate the effects of information treatments on survey engagement, information recall, consequentiality, and stated preferences, similar to @welling2023information, expanding our understanding of treatment effects
+-   We test who chooses additional information and to what extent their preferences differ from those who do not choose additional information
+:::
+
+## Research Questions
+
+::: incremental
+1.  Do obligatory and optional information provision affect **survey engagement**, **information recall**, **consequentiality**, and **stated preferences**?
+
+- OLS and MXL with interactions
+
+2.  Do **socio-demographic** or **attitudinal** variables influence the decision to **access optional information**?
+
+- Logit regression
+
+3.  Do **survey engagement**, **information recall**, **consequentiality**, and **stated preferences** differ between respondents who **voluntary access information** and those who do not?
+
+- OLS and MXL with interactions
+:::
+
+# Survey & Data
+
+## Discrete Choice Experiment
+
+::: incremental
+-   To investigate the research questions, we use data from a **discrete choice experiment (DCE)** on naturalness of urban green spaces in the 14 largest German cities
+-   The survey is an exact **replication** of the choice experiment of @Bronnmann062321-0072R1 and differs only in the information provided to the respondents
+-   In the DCE, respondents were asked to imagine possible **changes** to their **most frequently used UGS**
+-   This **restructuring** involved adjustments to the UGS's **naturalness** and changes to the **walking distance**
+-   The associated **costs** of this restructuring were intended to be integrated into monthly **rental payments**
+<!-- -   Participants in the DCE were presented **ten** randomly assigned **choice cards** with a choice between **two alternative programs** for the renovation of the UGS and the **current status quo**. -->
+:::
+
+## Choice Card
+
+![](images/Figure%202.PNG){width="300"}
+
+## Treatment (Information Provision)
+
+-   Info text about the effect of **natural urban green spaces** on urban **heat islands**
+-   **Optional video** (2 minutes) with the almost the same information
+
+![](images/waermeinsel.png){width="200"}
+
+## Treatment (Quiz & Self Reference Questions)
+
+::: incremental
+-   **Seven quiz questions** with strict reference to the previously provided information
+
+-   Example: *The temperature difference between the city and the surrounding area can be up to 10 degrees Celsius. (true/false)*
+
+-   **Two self reference questions**
+
+-   Example: *The city should do more to avoid heat islands. (Strongly agree - Strongly disagree)*
+:::
+
+
+
+
+
+
+
+## Additional Data
+
+::: incremental
+
++   Quiz: Evaluation of the quiz we gave to everyone after the DCE.**-\>Information recall**
+
++   Timings: We saved the net interview time and the mean Choice Card time.-\> **Survey engagement**
+
++   Two questions on **consequentiality **
+
+- Example: *To what extent do you believe that the decisions you make will have an impact on how the green spaces in your neighborhood are designed in the future? (I believe in it very much - I don’t believe in it at all)*
+
++   **Socio-demographics**: Age, Gender, Income, Education
+
++   **Attitudinal variable**: Measure derived from 21 items on **nature relatedness** (NR-Index) [@nisbet2009nature]
+
+:::
+
+
+## Experimental Setting
+
+![](Grafics/FlowChart_4_groups.png){width="300"}
+
+
+<!-- ## Methods (1) {auto-animate="true"} -->
+
+
+<!-- -   OLS regression (survey engagement, information recall, consequentiality): -->
+
+<!-- ```{=tex} -->
+<!-- \begin{equation} -->
+<!--     Y = \beta_0 +  \beta_{Treat} \cdot v_{Treat} + \beta_{Control} \cdot v_{Control} + \epsilon -->
+<!--     \label{ols} -->
+<!-- \end{equation} -->
+<!-- ``` -->
+<!-- -   Mixed logit model with interactions in WTP space: -->
+
+
+<!-- ## Methods (2) {auto-animate="true"} -->
+
+<!-- -   Mixed logit model with interactions in WTP space: -->
+
+<!-- ```{=tex} -->
+<!-- \begin{equation} -->
+<!--     U_i = -(\beta_{C_i} + \beta_{TreatC_i} \cdot v_{Treat}) \cdot (\beta_{X_i} \cdot v_{X_i} + \beta_{TreatX_i} \cdot v_{X_i} \cdot v_{Treat} - C_i) + \epsilon_i -->
+<!-- \end{equation} -->
+<!-- ``` -->
+<!-- with -->
+
+<!-- ```{=tex} -->
+<!-- \begin{equation} -->
+<!--   v_{X_i} = \{ASC_{sq_i}, Nat_i, WD_i\} -->
+<!-- \end{equation} -->
+<!-- ``` -->
+<!-- and -->
+
+<!-- ```{=tex} -->
+<!-- \begin{equation} -->
+<!--   v_{Treat_A} = \{Treated, Optional Treatment\} -->
+<!-- \end{equation} -->
+<!-- ``` -->
+<!-- ```{=tex} -->
+<!-- \begin{equation} -->
+<!--   v_{Treat_B} = \{Treated, Vol. Treated, No Info\} -->
+<!-- \end{equation} -->
+<!-- ``` -->
+
+# Case A: Obligatory vs. Optional Information
+
+
+## Case A
+
+1.  Do obligatory and optional information provision affect **survey engagement**, **information recall**, **consequentiality**, and **stated preferences**?
+
+![](Grafics/FlowChart_4_groups_A.png){width="300"}
+
+<!-- ::: {style="font-size: 45%;"} -->
+
+<!-- ```{r, results='asis'} -->
+
+<!-- htmlreg(l=list(ols_time_spent_control_A, ols_time_cc_control_A, ols_percentage_correct_control_A, conseq_model_control_A), -->
+
+<!--        custom.model.names = c("Interview Time",  "CC Time", "Quiz",  "Cons. Score"), -->
+
+<!--        custom.header = list("Model 1A" = 1:1, "Model 2A" = 2:2, "Model 3A" = 3:3, "Model 4A" = 4:4), -->
+
+<!--        custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",  -->
+
+<!--        label = "tab:olsA", -->
+
+<!--        caption = "Results of OLS regressions for Scenario Case A.") -->
+
+<!-- ``` -->
+
+<!-- ::: -->
+
+## OLS: Engagement
+
+::: panel-tabset
+
+### OLS Specification
+
+```{=tex}
+\begin{equation}
+    Y = \beta_0 +  \beta_{Treat} \cdot v_{Treat} + \beta_{SocDem} \cdot v_{SocDem} + \epsilon
+    \label{olss}
+\end{equation}
+```
+
+with
+
+```{=tex}
+\begin{equation}
+  Y = \left\{
+  \begin{aligned}
+    &\text{Net Interview Time, Mean Choice Card Time} \\
+    &\text{Percentage of correct quiz statements, Consequentiality Score} 
+  \end{aligned}
+  \right\}
+\end{equation}
+```
+
+
+```{=tex}
+\begin{equation}
+  v_{Treat} = \{\text{Treated, Optional Treatment}\}
+\end{equation}
+```
+
+```{=tex}
+\begin{equation}
+  v_{SocDem} = \{\text{Age, Gender, Income, Education, NR-Index}\}
+\end{equation}
+```
+
+
+### Results
+
+```{r}
+ggpubr::ggarrange(plot_interview_A, plot_cc_A)
+```
+
+:::
+
+## OLS: Information Recall & Consequentiality
+
+```{r}
+ggpubr::ggarrange(plot_mani_A, plot_cons_A)
+```
+
+## MXL: Effects on Stated Preferences
+
+::: panel-tabset
+
+### MXL Specification
+
+```{=tex}
+\begin{equation}
+    U_i = -(\beta_{C_i} + \beta_{TreatC_i} \cdot v_{Treat}) \cdot (\beta_{X_i} \cdot v_{X_i} + \beta_{TreatX_i} \cdot v_{X_i} \cdot v_{Treat} - C_i) + \epsilon_i
+    \label{mxl_base}
+\end{equation}
+```
+with
+
+```{=tex}
+\begin{equation}
+  v_{X_i} = \{ASC_{sq_i}, Nat_i, WD_i\}
+\end{equation}
+```
+
+```{=tex}
+\begin{equation}
+  v_{Treat} = \{\text{Treated, Optional Treatment}\}
+\end{equation}
+```
+```{=tex}
+\begin{equation}
+  C_i = \{Rent_i\}
+\end{equation}
+```
+
+### Results 
+
+::: {style="font-size: 90%;"}
+```{r, results='asis'}
+htmlreg(c(case_A_cols[1], remGOF(case_A_cols[2:4])), single.row = TRUE,
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
+                              "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+                              "wd" = "Walking Distance", "asc" = "ASC SQ"),
+       custom.model.names = c("Mean", "SD", "Treated", "Optional Treatment"), custom.note = "%stars. Standard errors in parentheses.",
+       stars = c(0.01, 0.05, 0.1), float.pos="tb", 
+       label = "tab:mxl_A",
+       caption = "Results of mixed logit model with treatment interactions for Case A.")
+```
+:::
+
+
+:::
+
+# Case B: Voluntary Information Access
+
+## Voluntary Information Access
+
+2.  Do **socio-demographic** or **attitudinal** variables influence the decision to **access optional information**?
+
+![](Grafics/FlowChart_Optional_only.png){width="300"}
+
+
+
+## Logit Regression: Who chooses Optional Information?
+
+```{=tex}
+\begin{equation}
+    Y = \beta_0 + \beta_{SocDem} \cdot v_{SocDem} + \epsilon
+    \label{simple_logit}
+\end{equation}
+```
+
+::: {style="font-size: 58%;"}
+```{r, results='asis'}
+
+
+htmlreg(l=list(logit_choice_treat_uni), stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       custom.model.names = c("Logit regression"),
+       custom.header = list("Dependent variable: Voluntary Information Access" = 1),
+       custom.coef.map = list_ols, custom.note = "%stars. Standard errors in parentheses.",
+       label = "tab:logit_vt", single.row = TRUE,
+       caption = "Results of logit regression on access of optional information.")
+
+
+```
+:::
+
+
+## Case B
+
+3.  Do **survey engagement**, **information recall**, **consequentiality**, and **stated preferences** differ between respondents who **voluntary access information** and those who do not?
+
+![](Grafics/FlowChart_4_groups_B.png){width="300"}
+
+## OLS Engagement: Interview & Choice Card Time
+
+```{r}
+ggpubr::ggarrange(plot_interview_D, plot_cc_D)
+```
+
+## OLS: Information Recall & Consequentiality
+
+```{r}
+ggpubr::ggarrange(plot_mani_D, plot_cons_D)
+```
+
+## MXL: Case B
+
+::: {style="font-size: 80%;"}
+```{r, results='asis'}
+htmlreg(c(case_B_cols[1], remGOF(case_B_cols[2:5])),
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
+                              "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+                              "wd" = "Walking Distance", "asc" = "ASC SQ",
+                              "ASC_sq_info" = "ASC SQ", "rent_info" = "Rent", "nat_info" = "Naturalness", "walking_info" = "Walking Distance"),
+       custom.model.names = c("Mean", "SD", "Treated", "Vol. Treated", "No Info"), custom.note = "%stars (one-sided). Robust standard errors in parentheses.",
+       stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       label = "tab:mxl_C",
+       caption = "Results of mixed logit model with treatment interactions for Case B.")
+```
+:::
+
+## MXL: Case B with NR-index Interaction
+
+::: {style="font-size: 80%;"}
+```{r, results='asis'}
+htmlreg(c(case_B_cols_NR[1], remGOF(case_B_cols_NR[2:6])),
+       custom.coef.map = list("natural" = "Naturalness", "walking" = "Walking Distance", "rent" = "Rent",
+                              "ASC_sq" = "ASC SQ", "_natural" = "Naturalness", "nat" = "Naturalness",
+                              "wd" = "Walking Distance", "asc" = "ASC SQ", "ASC" ="ASC SQ",
+                              "ASC_sq_info" = "ASC SQ", "rent_info" = "Rent", "nat_info" = "Naturalness", "walking_info" = "Walking Distance"),
+       custom.model.names = c("Mean", "SD", "Treated", "Vol. Treated", "No Info", "NR-Index"), custom.note = "%stars (one-sided). Robust standard errors in parentheses.",
+       stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       label = "tab:mxl_C_NR",
+       caption = "Results of mixed logit model with treatment interactions for Case B including NR-index.")
+```
+:::
+
+# Discussion
+
+## Discussion (1)
+
+1.  Do obligatory and optional information provision affect survey engagement, information recall, consequentiality, and stated preferences?
+
+::: incremental
+-   Obligatory and optional treatments do not increase survey engagement measured via time spend on the survey
+
+-   Both treatments increase information recall, stronger effect for obligatory treatment
+
+-   No effects on consequentiality
+
+-   Strong effects on stated preferences for both treatments, more pronounced effect for the obligatory treatment
+:::
+
+## Discussion (2)
+
+2.  Do socio-demographic or attitudinal variables influence the decision to access optional information?
+
+::: incremental
+-   Respondents that voluntary access information are younger, richer and have a higher natural relatedness index
+
+-   No effects of gender and education
+
+-   Respondents' preferences for the good to be valued (NR-index) influence the likelihood of accessing additional information
+:::
+
+## Discussion (3)
+
+3.  Do survey engagement, information recall, consequentiality, and stated preferences differ between respondents who voluntary access information and those who do not?
+
+::: incremental
+-   Respondents that voluntary access information do not engage more in the survey, but perform best in the quiz
+
+-   Respondents that decide to not access additional information engage less in the survey, have a lower consequentiality score and do not perform different in the quiz than the non-treated respondents
+
+-   Highest willingness to pay values in the group that voluntary accesses information
+:::
+
+## Conclusion
+
+::: incremental
+-   Providing optional information does not lead to optional information seeking
+
+-   Optional information is mostly accessed by people that are interested in the good to be valued
+
+-   We recommend to use obligatory information provision rather than optional one
+:::
+
+## References
+
+::: {#refs}
+:::
+
+# Appendix {.appendix}
+
+
+Information provision (Video) Link to the video: <https://idiv.limequery.com/upload/surveys/682191/files/urban-heat-island-effekt.mp4>
+
+## Summary Statistics
+
+::: {style="font-size: 55%;"}
+```{=latex}
+        \begin{table}[htbp]
+        \centering
+        \begin{tabular}{lcrcclc} \hline 
+         Characteristic& Sample & Treated & Opt. Treatment & Non-Treated & Vol. Treated &  No Info   \\ \hline
+         Gender & &&& &   \\
+          \hspace{0.3cm} Proportion Females&50.67\%&50.15\%&47.98\%&49.76\% &50.48\%  &  48.26\%  \\ \hline
+         Age &44.60&44.91 &44.24&44.63 & 42.51 &  44.97 \\ \hline
+         Net Household Income&&&&  & & \\
+          \hspace{0.3cm} Less than 1500\euro  &19.06\%&20.43\%&18.86\%&17.85\%&  22.12\% &  19.77\% \\ 
+          \hspace{0.3cm} 1500\euro -3000\euro &23.39\%&20.59\%&24.75\%&24.96\%&  15.87\% &  26.16\% \\
+          \hspace{0.3cm} 3000\euro -4000\euro &32.62\%&35.14\%&29.63\%&32.86\%&  34.13\% & 27.61\% \\ 
+          \hspace{0.3cm} More than 4000\euro  &19.65\%&19.04\%&20.70\%&19.27\%&   24.52\% & 18.31\% \\ \hline
+                   Education & &&&&  &   \\
+          \hspace{0.3cm} University degree &34.54\%&32.20\%&34.68\%&36.81\%&  33.17\% & 33.72\% \\ 
+          \hspace{0.3cm} No University degree &64.46\%&67.80\%&65.32\%&63.19\%&  66.83\%  &66.28\% \\ \hline
+          NR-index &0.00&0.00&-0.02&0.02&  0.15  &-0.19 \\ \hline
+           Correct quiz answers   &66.97\%&69.86\%&66.91\%&64.07\%&69.44\%&63.08\% \\ \hline
+           Consequentiality  score  &6.15&6.28&6.02&6.12&6.75   &5.76 \\ \hline
+          Net interview time   &1396&1340&1424&1426& 1441  &1350  \\
+          Mean Choice Card time &18.86&19.02&19.10&18.51 &19.03&19.24 \\ \hline
+                   Number of respondents &1873&646&594&633&  250  &344  \\ \hline
+                   \multicolumn{7}{p{\textwidth-2\tabcolsep}}{\scriptsize Notes: Age is measured in years; NR-index is z-standardized with mean=0 and standard deviation=1; Correct quiz answers measure the percentage of correctly answered quiz questions; Consequentiality score is the sum of the two Likert-type questions on consequentiality; Net interview time and Mean Choice card time are measured in seconds; Net interview time is measured as overall interview time minus treatment time; To avoid bias, we excluded the fastest 1\% and the slowest 1\% of net interview time and mean Choice Card time from the analysis.}
+    \end{tabular}
+        \caption{Summary statistics of the sample and by scenarios.}
+        \label{tab:summary}
+        \end{table}
+```
+:::
+
+
+## NR OLS
+
+::: {style="font-size: 65%;"}
+```{r, results='asis'}
+htmlreg(l=list(nr_model_treat_A), single.row = TRUE,
+       custom.model.names = c("OLS regression"),
+       custom.header = list("Dependent variable: NR-Index" = 1),
+       custom.coef.map = list("(Intercept)" = "Intercept", "as.factor(Treatment_A)Treated" = "Treated", "as.factor(Treatment_A)Vol_Treated" = "Vol. Treated",
+                              "as.factor(Treatment_C)No Info 2" = "No Info 2", "as.factor(Treatment_C)No Video 1" = "Text 1",
+                              "as.factor(Treatment_C)No Video 2" = "Text 2", "as.factor(Treatment_C)Video 1" = "Video 1",
+                              "as.factor(Treatment_C)Video 2" = "Video 2", "Z_Mean_NR" = "NR-Index", "as.factor(Gender)2" = "Female",
+                              "Age_mean" = "Age", "QFIncome" = "Income", "Uni_degree" = "University Degree", "Kids_Dummy" = "Children",
+                              "Naturalness_SQ" = "Naturalness SQ", "WalkingDistance_SQ" = "Walking Distance SQ"),
+       stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       custom.note = "%stars. Standard errors in parentheses.",
+       label = "tab:nr_ols",
+        caption = "Results of OLS on the NR-index.")
+```
+:::
+
+## OLS Models: Case A
+
+
+::: {style="font-size: 65%;"}
+```{r, results='asis'}
+htmlreg(l=list(ols_time_spent_control_A, ols_time_cc_control_A, ols_percentage_correct_control_A, conseq_model_control_A),
+       custom.model.names = c("Interview Time",  "CC Time", "Quiz",  "Cons. Score"),
+       custom.header = list("Model 1A" = 1:1, "Model 2A" = 2:2, "Model 3A" = 3:3, "Model 4A" = 4:4),
+       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb", 
+       custom.note = "Notes: (i) The dependent variables examined in this regression analysis are as follows: 
+       Model 1A refers to net interview time, 
+       Model 2A denotes choice card time, Model 3A represents the percentage of correct quiz questions, and Model 4A represents consequentiality score.
+       (ii) The variables included in the analysis are as follows: Treated is a dummy variable indicating membership 
+       in the obligatory treated group, Optional Treatment is a dummy variable indicating membership in the group with 
+       the choice to receive treatment or not, with the reference group being Non-Treated. NR-Index represents the z-standardized
+       natural relatedness index. Female is a dummy variable denoting gender, Age has been mean-centered and measured in years, 
+       Income is a continuous variable indicating a transition from one income group to the next higher, and University Degree is 
+       a dummy variable indicating whether an individual holds a university degree; (iii) %stars and standard errors in parentheses.",
+       label = "tab:olsA", single.row = TRUE,
+       caption = "Results of OLS regressions for Scenario Case A.")
+```
+:::
+
+## OLS Models: Case B
+
+::: {style="font-size: 60%;"}
+```{r, results='asis'}
+htmlreg(l=list(ols_time_spent_control_D, ols_time_cc_control_D, ols_percentage_correct_control_D,  conseq_model_control_D),
+       custom.model.names = c("Interview Time",  "CC Time", "Quiz",  "Cons. Score"),
+       custom.header = list("Model 1B" = 1:1, "Model 2B" = 2:2, "Model 3B" = 3:3, "Model 4B" = 4:4),
+       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb", 
+       custom.note = "Notes: (i) The dependent variables examined in this regression analysis are as follows: 
+       Model 1B refers to net interview time, 
+       Model 2B denotes choice card time, Model 1B represents the percentage of correct quiz questions, and Model 4B represents consequentiality score.
+       (ii) The variables included in the analysis are as follows: Treated is a dummy variable indicating membership 
+       in the obligatory treated group, Vol. Treated is a dummy variable indicating the group that voluntarily chose the optional treatment, 
+       while No Info indicates the group that did not opt for the treatment, with the reference group being Non-Treated. NR-Index represents the z-standardized
+       natural relatedness index. Female is a dummy variable denoting gender, Age has been mean-centered and measured in years, 
+       Income is a continuous variable indicating a transition from one income group to the next higher, and University Degree is 
+       a dummy variable indicating whether an individual holds a university degree; (iii) %stars and standard errors in parentheses.",
+       label = "tab:olsD",
+       caption = "Results of OLS regressions for Scenario Case B.", single.row = TRUE)
+```
+:::
+
+## MXL: Split Samples
+
+```{r}
+ggplot(data=mxl_melt_info, aes(x=Coefficent, y=abs(value), fill=variable)) +
+  geom_bar(stat="identity",  position='dodge', width = 0.9) +
+  geom_errorbar(aes(x=Coefficent, ymin=abs(value)-ME, ymax=abs(value)+ME), width=0.3, position=position_dodge(0.8)) +
+  ylab("Absolute Value") +
+  xlab("Coefficient") +
+  scale_x_discrete(guide = guide_axis(angle = 45)) +
+  scale_fill_brewer(palette = "Set2", labels = c("Treated", "Optional Treatment", "Not Treated"), name="Treatment") +
+  theme(legend.position = c(0.85, 0.8)) 
+```