diff --git a/Scripts/MAKE_FILE.R b/Scripts/MAKE_FILE.R
index be870b69e147a6e6cef2f88f287922bc7f55f260..a165a46abbee8521fdece9f5dbd32a626aab29cc 100644
--- a/Scripts/MAKE_FILE.R
+++ b/Scripts/MAKE_FILE.R
@@ -1,100 +1,105 @@
-rm(list=ls())
-library(tidyverse)
-library(tidylog)
-library(apollo)
-library(reshape2)
-library(xtable)
-library(stargazer)
-library(texreg)
-#test
-# Set values for estimation in Apollo
-n_draws <- 2000
-n_cores <- min(parallel::detectCores()-1, 25)
-
-# Load data
-load("Data/database_full.RData")
-load("Data/database.RData")
-
-# Data preparation
-source("Scripts/data_prep.R")
-source("Scripts/treatment.R")
-
-####### Estimate models ######
-
-### Logit
-source("Scripts/logit/chr_vol_treat.R")
-source("Scripts/logit/protesters.R")
-       
-### OLS
-source("Scripts/ols/ols_time_spent.R")
-source("Scripts/ols/ols_quiz.R")
-source("Scripts/ols/ols_opt_out.R")
-source("Scripts/ols/ols_nr.R")
-source("Scripts/ols/ols_consequentiality.R")
-
-##### Conditional Logits #####
-#source("Scripts/clogit.R")
-#source("Scripts/clogit_wtp.R")
-
-##### Mixed Logit Models ######
-
-#source("Scripts/mxl/mxl_wtp_space.R")
-#source("Scripts/mxl/mxl_wtp_space_4d.R")
-#source("Scripts/mxl/mxl_wtp_space_interact_all.R")
-#source("Scripts/mxl/mxl_socio_int.R")
-#source("Scripts/mxl/mxl_treatment_time.R")
-#source("Scripts/mxl/mxl_treatment_time_interaction.R")
-#source("Scripts/mxl/mxl_treatment_time_Dummies.R")
-#source("Scripts/mxl/mxl_wtp_space_interact_everything.R")
-#source("Scripts/mxl/mxl_wtp_space_4d_interact_everything.R")
-#############################
-
-##### Load models ############
-
-mxl_wtp <- apollo_loadModel("Estimation_results/mxl/MXL_wtp")
-mxl_wtp_4d <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_4d")
-mxl_wtp_all_int <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_interact_all")
-mxl_wtp_socio <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_socio_int")
-mxl_wtp_tt <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_treatment_time")
-mxl_wtp_tt_interaction <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_treatment_time_interaction")
-mxl_wtp_time_groups <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_time_groups")
-mxl_wtp_everything <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_interact_everything")
-mxl_wtp_case_a <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case A")
-mxl_wtp_case_a_NR <- apollo_loadModel("Estimation_results/mxl/MXL_wtp NR A")
-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")
-
-# rent interactions models
-mxl_wtp_case_a_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case A Rent Int")
-mxl_wtp_case_b_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case B Rent Int")
-mxl_wtp_case_c_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_Case_C Rent INT")
-mxl_wtp_case_d_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_Case_D")
-
-# rent interactions models NR
-mxl_wtp_NR_case_a_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp NR A Rent INT")
-mxl_wtp_NR_case_c_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_NR_Case_C RENT INT X")
-
-# Alternative case
-case_d <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case D Rent Int")
-
-##############################
-
-# Model analysis
-source("Scripts/visualize_models.R")
-
-source("Scripts/compare_split_samples.R")
-
-source("Scripts/create_tables.R")
-
-source("Scripts/interaction_plots_presi.R")
-
-
-### Old models ###
-
-
-# # without protesters
-# mxl_wtp_case_a_prot <- apollo_loadModel("Estimation_results/mxl/without_protesters/MXL_wtp Case A prot")
-# mxl_wtp_case_b_prot <- apollo_loadModel("Estimation_results/mxl/without_protesters/MXL_wtp Case B prot")
+rm(list=ls())
+library(tidyverse)
+library(tidylog)
+library(apollo)
+library(reshape2)
+library(xtable)
+library(stargazer)
+library(texreg)
+#test
+# Set values for estimation in Apollo
+n_draws <- 2000
+n_cores <- min(parallel::detectCores()-1, 25)
+
+# Load data
+load("Data/database_full.RData")
+load("Data/database.RData")
+
+# Data preparation
+source("Scripts/data_prep.R")
+source("Scripts/treatment.R")
+
+####### Estimate models ######
+
+### Logit
+source("Scripts/logit/chr_vol_treat.R")
+source("Scripts/logit/protesters.R")
+       
+### OLS
+source("Scripts/ols/ols_time_spent.R")
+source("Scripts/ols/ols_quiz.R")
+source("Scripts/ols/ols_opt_out.R")
+source("Scripts/ols/ols_nr.R")
+source("Scripts/ols/ols_consequentiality.R")
+
+##### Conditional Logits #####
+#source("Scripts/clogit.R")
+#source("Scripts/clogit_wtp.R")
+
+##### Mixed Logit Models ######
+
+#source("Scripts/mxl/mxl_wtp_space.R")
+#source("Scripts/mxl/mxl_wtp_space_4d.R")
+#source("Scripts/mxl/mxl_wtp_space_interact_all.R")
+#source("Scripts/mxl/mxl_socio_int.R")
+#source("Scripts/mxl/mxl_treatment_time.R")
+#source("Scripts/mxl/mxl_treatment_time_interaction.R")
+#source("Scripts/mxl/mxl_treatment_time_Dummies.R")
+#source("Scripts/mxl/mxl_wtp_space_interact_everything.R")
+#source("Scripts/mxl/mxl_wtp_space_4d_interact_everything.R")
+#############################
+
+##### Load models ############
+
+mxl_wtp <- apollo_loadModel("Estimation_results/mxl/MXL_wtp")
+mxl_wtp_4d <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_4d")
+mxl_wtp_all_int <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_interact_all")
+mxl_wtp_socio <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_socio_int")
+mxl_wtp_tt <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_treatment_time")
+mxl_wtp_tt_interaction <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_treatment_time_interaction")
+mxl_wtp_time_groups <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_time_groups")
+mxl_wtp_everything <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_interact_everything")
+mxl_wtp_case_a <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case A")
+mxl_wtp_case_a_NR <- apollo_loadModel("Estimation_results/mxl/MXL_wtp NR A")
+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")
+
+# rent interactions models
+mxl_wtp_case_a_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case A Rent Int")
+mxl_wtp_case_b_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case B Rent Int")
+mxl_wtp_case_c_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_Case_C Rent INT")
+mxl_wtp_case_d_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_Case_D")
+
+# rent interactions models NR
+mxl_wtp_NR_case_a_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp NR A Rent INT")
+mxl_wtp_NR_case_c_rentINT <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_NR_Case_C RENT INT X")
+
+# Alternative case
+case_d <- apollo_loadModel("Estimation_results/mxl/MXL_wtp Case D Rent Int")
+
+# New Case Text Video merged
+new_case_b <- apollo_loadModel("Estimation_results/mxl/MXL_wtp_Case_D")
+
+new_case_b_NR <- apollo_loadModel("Estimation_results/mxl/D_NR")
+
+##############################
+
+# Model analysis
+source("Scripts/visualize_models.R")
+
+source("Scripts/compare_split_samples.R")
+
+source("Scripts/create_tables.R")
+
+source("Scripts/interaction_plots_presi.R")
+
+
+### Old models ###
+
+
+# # without protesters
+# mxl_wtp_case_a_prot <- apollo_loadModel("Estimation_results/mxl/without_protesters/MXL_wtp Case A prot")
+# mxl_wtp_case_b_prot <- apollo_loadModel("Estimation_results/mxl/without_protesters/MXL_wtp Case B prot")
 # mxl_wtp_case_c_prot <- apollo_loadModel("Estimation_results/mxl/without_protesters/MXL_wtp_Case_C prot")
\ No newline at end of file
diff --git a/Scripts/create_tables.R b/Scripts/create_tables.R
index 45600ebe9b2d13d0c9d73a73b421b557e1109bb7..266d226f2f182fa92e9d83fe612bec792d40effe 100644
--- a/Scripts/create_tables.R
+++ b/Scripts/create_tables.R
@@ -1,203 +1,252 @@
-library(choiceTools)
-
-dir.create("Tables/mxl")
-dir.create("Tables/logit")
-dir.create("Tables/ols/")
-
-
-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")
-# OLS A
-texreg(l=list(ols_percentage_correct_control_A, ols_time_spent_control_A, ols_time_cc_control_A, conseq_model_control_A),
-       custom.model.names = c("Quiz", "Interview Time",  "CC Time",  "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 represents the percentage of correct quiz questions, Model 2A refers to net interview time, 
-       Model 3A denotes choice card time, 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",
-       caption = "Results of OLS regressions for Scenario Case A.",
-       file="Tables/ols/ols_A.tex")
-
-# OLS D
-texreg(l=list(ols_percentage_correct_control_D, ols_time_spent_control_D, ols_time_cc_control_D, conseq_model_control_D),
-       custom.model.names = c("Quiz", "Interview Time",  "CC Time",  "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 represents the percentage of correct quiz questions, Model 2B refers to net interview time, 
-       Model 3B denotes choice card time, 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.",
-       file="Tables/ols/ols_D.tex")
-
-# Manipulation check
-texreg(l=list(ols_percentage_correct_A,  ols_percentage_correct_control_A, ols_percentage_correct_C, ols_percentage_correct_control_C),
-       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
-       custom.header = list("Dependent Variable: Percentage of correct quiz statements" = 1:4),
-       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb", 
-       custom.note = "%stars. Standard errors in parentheses.",
-       label = "tab:mani",
-       caption = "Results of OLS on percentage of correct quiz statements.",
-       file="Tables/ols/manipulation.tex")
-
-
-# Net interview time 
-texreg(l=list(ols_time_spent_A,  ols_time_spent_control_A, ols_time_spent_C,  ols_time_spent_control_C),
-       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
-       custom.header = list("Dependent variable: Net interview time" = 1:4),
-       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
-       custom.note = "%stars. Standard errors in parentheses.",
-       label = "tab:net_int",
-       caption = "Results of OLS on net interview time.",
-       file="Tables/ols/interviewtime.tex")
-
-
-# CC Time
-texreg(l=list(ols_time_cc_A,  ols_time_cc_control_A, ols_time_cc_C,  ols_time_cc_control_C),
-       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
-       custom.header = list("Dependent variable: Mean choice card time" = 1:4),
-       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
-       custom.note = "%stars. Standard errors in parentheses.",
-       label = "tab:cctime",
-       caption = "Results of OLS on mean choice card time.",
-       file="Tables/ols/cctime.tex")
-
-# Consequentiality
-texreg(l=list(conseq_model_A, conseq_model_control_A, conseq_model_C, conseq_model_control_C),
-       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
-       custom.header = list("Dependent variable: Consequentiality score" = 1:4),
-       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
-       custom.note = "%stars. Standard errors in parentheses.",
-       label = "tab:conseq",
-       caption = "Results of OLS on consequentiality score.",
-       file="Tables/ols/consequentiality.tex")
-
-# Opt Out
-texreg(l=list(ols_opt_out_A, ols_opt_out_control_A, ols_opt_out_C, ols_opt_out_control_C),
-       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
-       custom.header = list("Dependent variable: Number of opt-out choices" = 1:4),
-       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
-       custom.note = "%stars. Standard errors in parentheses.",
-       label = "tab:optout",
-       caption = "Results of OLS on number of opt-out choices.",
-       file="Tables/ols/optout.tex")
-
-# NR
-texreg(l=list(nr_model_treat_A),
-       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.",
-       file="Tables/ols/nr_ols.tex")
-
-#### Logit #####
-
-texreg(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",
-       caption = "Results of logit regression on the access of optional information.",
-       file="Tables/logit/chose_treatment.tex")
-
-
-##### MXL #######
-
-### Baseline case A
-case_A <- quicktexregapollo(mxl_wtp_case_a_rentINT)
-
-coef_names <- case_A@coef.names 
-coef_names <- sub("^(mu_)(.*)(_T|_VT)$", "\\2\\3", coef_names)
-coef_names[4] <- "mu_ASC_sq"
-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",
-                              "wd" = "Walking Distance", "asc" = "ASC SQ"),
-       custom.model.names = c("Mean", "SD", "Treated", "Voluntary Treated"), custom.note = "%stars. 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")
-
-### Baseline case C
-case_C <- quicktexregapollo(mxl_wtp_case_c_rentINT)
-
-coef_names <- case_C@coef.names 
-coef_names <- sub("^(mu_)(.*)(1|2|info)$", "\\2\\3", coef_names)
-coef_names[4] <- "mu_ASC_sq"
-case_C@coef.names <- coef_names
-
-
-case_C_cols <- map(c("^mu_", "^sig_", "_vid1$", "_vid2$", "_nv1$", "_nv2$", "_no_info$"), subcoef, case_C)
-
-texreg(c(case_C_cols[1], remGOF(case_C_cols[2:7])),
-       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", "Video 1", "Video 2", "Text 1", "Text 2", "No Info"), custom.note = "%stars. 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.",
-       file="Tables/mxl/case_C_rent_INT.tex")
-
-### Rent NR model case C
-case_C_NR <- quicktexregapollo(mxl_wtp_NR_case_c_rentINT)
-
-coef_names <- case_C_NR@coef.names 
-coef_names <- sub("^(mu_)(.*)(1|2|info|NR)$", "\\2\\3", coef_names)
-coef_names[4] <- "mu_ASC_sq"
-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",
-                              "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.",
-       stars = c(0.01, 0.05, 0.1), float.pos="tb",
-       label = "tab:mxl_NR",
-       caption = "Results of mixed logit model with treatment and NR-index interactions for Case B.",
-       file="Tables/mxl/case_C_rent_INT_NR.tex")
-# Main model
-# texreg(l=list(mxl_wtp_case_a_rentINT),
-#        custom.coef.map = list("mu_natural" = "Naturalness", "mu_walking" = "Walking Distance", "mu_rent" = "Rent",
-#                               "ASC_sq" = "ASC SQ", "sig_natural" = "Naturalness SD", "sig_walking" = "Walking Distance SD",
-#                               "sig_rent" = "Rent SD", "sig_ASC_sq" = "ASC SD",
-#                               "mu_nat_T" = "Naturalness X Treated", "mu_wd_T" = "Walking Distance X Treated", "mu_rent_T" = "Rent X Treated",
-#                               "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")
+library(choiceTools)
+
+dir.create("Tables/mxl")
+dir.create("Tables/logit")
+dir.create("Tables/ols/")
+
+
+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")
+# OLS A
+texreg(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",
+       caption = "Results of OLS regressions for Scenario Case A.",
+       file="Tables/ols/ols_A.tex")
+
+# OLS D
+texreg(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.",
+       file="Tables/ols/ols_D.tex")
+
+# Manipulation check
+texreg(l=list(ols_percentage_correct_A,  ols_percentage_correct_control_A, ols_percentage_correct_C, ols_percentage_correct_control_C),
+       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
+       custom.header = list("Dependent Variable: Percentage of correct quiz statements" = 1:4),
+       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb", 
+       custom.note = "%stars. Standard errors in parentheses.",
+       label = "tab:mani",
+       caption = "Results of OLS on percentage of correct quiz statements.",
+       file="Tables/ols/manipulation.tex")
+
+
+# Net interview time 
+texreg(l=list(ols_time_spent_A,  ols_time_spent_control_A, ols_time_spent_C,  ols_time_spent_control_C),
+       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
+       custom.header = list("Dependent variable: Net interview time" = 1:4),
+       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       custom.note = "%stars. Standard errors in parentheses.",
+       label = "tab:net_int",
+       caption = "Results of OLS on net interview time.",
+       file="Tables/ols/interviewtime.tex")
+
+
+# CC Time
+texreg(l=list(ols_time_cc_A,  ols_time_cc_control_A, ols_time_cc_C,  ols_time_cc_control_C),
+       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
+       custom.header = list("Dependent variable: Mean choice card time" = 1:4),
+       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       custom.note = "%stars. Standard errors in parentheses.",
+       label = "tab:cctime",
+       caption = "Results of OLS on mean choice card time.",
+       file="Tables/ols/cctime.tex")
+
+# Consequentiality
+texreg(l=list(conseq_model_A, conseq_model_control_A, conseq_model_C, conseq_model_control_C),
+       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
+       custom.header = list("Dependent variable: Consequentiality score" = 1:4),
+       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       custom.note = "%stars. Standard errors in parentheses.",
+       label = "tab:conseq",
+       caption = "Results of OLS on consequentiality score.",
+       file="Tables/ols/consequentiality.tex")
+
+# Opt Out
+texreg(l=list(ols_opt_out_A, ols_opt_out_control_A, ols_opt_out_C, ols_opt_out_control_C),
+       custom.model.names = c("Case A", "with Controls", "Case B", "with Controls"),
+       custom.header = list("Dependent variable: Number of opt-out choices" = 1:4),
+       custom.coef.map = list_ols,  stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       custom.note = "%stars. Standard errors in parentheses.",
+       label = "tab:optout",
+       caption = "Results of OLS on number of opt-out choices.",
+       file="Tables/ols/optout.tex")
+
+# NR
+texreg(l=list(nr_model_treat_A),
+       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.",
+       file="Tables/ols/nr_ols.tex")
+
+#### Logit #####
+
+texreg(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",
+       caption = "Results of logit regression on the access of optional information.",
+       file="Tables/logit/chose_treatment.tex")
+
+
+##### MXL #######
+
+### Baseline case A
+case_A <- quicktexregapollo(mxl_wtp_case_a_rentINT)
+
+coef_names <- case_A@coef.names 
+coef_names <- sub("^(mu_)(.*)(_T|_VT)$", "\\2\\3", coef_names)
+coef_names[4] <- "mu_ASC_sq"
+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",
+                              "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")
+
+### Baseline case C
+case_C <- quicktexregapollo(mxl_wtp_case_c_rentINT)
+
+coef_names <- case_C@coef.names 
+coef_names <- sub("^(mu_)(.*)(1|2|info)$", "\\2\\3", coef_names)
+coef_names[4] <- "mu_ASC_sq"
+case_C@coef.names <- coef_names
+
+
+case_C_cols <- map(c("^mu_", "^sig_", "_vid1$", "_vid2$", "_nv1$", "_nv2$", "_no_info$"), subcoef, case_C)
+
+texreg(c(case_C_cols[1], remGOF(case_C_cols[2:7])),
+       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", "Video 1", "Video 2", "Text 1", "Text 2", "No Info"), custom.note = "%stars. 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.",
+       file="Tables/mxl/case_C_rent_INT.tex")
+
+### Rent NR model case C
+case_C_NR <- quicktexregapollo(mxl_wtp_NR_case_c_rentINT)
+
+coef_names <- case_C_NR@coef.names 
+coef_names <- sub("^(mu_)(.*)(1|2|info|NR)$", "\\2\\3", coef_names)
+coef_names[4] <- "mu_ASC_sq"
+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",
+                              "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.",
+       stars = c(0.01, 0.05, 0.1), float.pos="tb",
+       label = "tab:mxl_NR",
+       caption = "Results of mixed logit model with treatment and NR-index interactions for Case B.",
+       file="Tables/mxl/case_C_rent_INT_NR.tex")
+
+
+### New Case B
+case_B <- quicktexregapollo(mxl_wtp_case_d_rentINT)
+
+coef_names <- case_B@coef.names 
+coef_names <- sub("^(mu_)(.*)(vol_treated|_treated|_no_info)$", "\\2\\3", coef_names)
+coef_names[4] <- "mu_ASC_sq"
+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",
+                              "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.",
+       file="Tables/mxl/case_B_rent_INT.tex")
+
+
+
+### New Case B NR
+case_B_NR <- quicktexregapollo(new_case_b_NR)
+
+coef_names <- case_B_NR@coef.names 
+coef_names <- sub("^(mu_)(.*)(vol_treated|_treated|_no_info|NR)$", "\\2\\3", coef_names)
+coef_names[4] <- "mu_ASC_sq"
+case_B_NR@coef.names <- coef_names
+
+
+case_B_cols_NR <- map(c("^mu_", "^sig_", "_treated$", "_vol_treated$","_no_info$", "_NR$"), subcoef, case_B_NR)
+
+texreg(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.",
+       file="Tables/mxl/case_B_rent_INT_NR.tex")
+
+
+# Main model
+# texreg(l=list(mxl_wtp_case_a_rentINT),
+#        custom.coef.map = list("mu_natural" = "Naturalness", "mu_walking" = "Walking Distance", "mu_rent" = "Rent",
+#                               "ASC_sq" = "ASC SQ", "sig_natural" = "Naturalness SD", "sig_walking" = "Walking Distance SD",
+#                               "sig_rent" = "Rent SD", "sig_ASC_sq" = "ASC SD",
+#                               "mu_nat_T" = "Naturalness X Treated", "mu_wd_T" = "Walking Distance X Treated", "mu_rent_T" = "Rent X Treated",
+#                               "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/mxl/mxl_wtp_space_caseD_RentINT_NR.R b/Scripts/mxl/mxl_wtp_space_caseD_RentINT_NR.R
index d0d686c41cd9b909b7619fb4928e056dfb489444..7e7b0be580ee84db7b993ce7a24b29854b631e3a 100644
--- a/Scripts/mxl/mxl_wtp_space_caseD_RentINT_NR.R
+++ b/Scripts/mxl/mxl_wtp_space_caseD_RentINT_NR.R
@@ -1,187 +1,187 @@
-#### 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",
-  modelDescr = "MXL wtp space Case D NR",
-  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_rent_treated = 0,
-              mu_rent_vol_treated = 0,
-              mu_rent_no_info = 0,
-              mu_rent_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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info 
-                  + mu_rent_NR * Z_Mean_NR)*
-                            (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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
-                  + mu_rent_NR * Z_Mean_NR)*                          
-                            (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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
-                  + mu_rent_NR * Z_Mean_NR)*
-                            (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_rent_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_rent_NR)
-
-
+#### 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",
+  modelDescr = "MXL wtp space Case D NR",
+  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_rent_treated = 0,
+              mu_rent_vol_treated = 0,
+              mu_rent_no_info = 0,
+              mu_rent_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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info 
+                  + mu_rent_NR * Z_Mean_NR)*
+                            (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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
+                  + mu_rent_NR * Z_Mean_NR)*                          
+                            (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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
+                  + mu_rent_NR * Z_Mean_NR)*
+                            (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_rent_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_rent_NR)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R b/Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R
index 2b16bb087a277655f659d410f95011ad3355c802..4682a96b61ce61af30f6b4756c76cd1fb295286c 100644
--- a/Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R
+++ b/Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R
@@ -1,198 +1,198 @@
-#### 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_X",
-  modelDescr = "MXL wtp space Case D Interactions",
-  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_ASC_Age = 0,
-              mu_ASC_Income = 0,
-              mu_rent_treated = 0,
-              mu_rent_vol_treated = 0,
-              mu_rent_no_info = 0,
-              mu_rent_NR = 0,
-              mu_rent_Age = 0,
-              mu_rent_Income = 0,
-              mu_nat_treated =0,
-              mu_nat_vol_treated = 0,
-              mu_nat_no_info = 0,
-              mu_nat_NR = 0,
-              mu_nat_Age = 0,
-              mu_nat_Income = 0,
-              mu_walking_treated =0,
-              mu_walking_vol_treated = 0,
-              mu_walking_no_info = 0,
-              mu_walking_NR = 0,
-              mu_walking_Age = 0,
-              mu_walking_Income = 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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info 
-                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*
-                            (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_nat_Age * Age_mean * Naturalness_1
-                            +  mu_nat_Income * Income_mean * Naturalness_1 + mu_walking_NR * Z_Mean_NR * WalkingDistance_1
-                            +  mu_walking_Age * Age_mean * WalkingDistance_1 + mu_walking_Income * Income_mean * WalkingDistance_1
-                            - Rent_1)
-  
-  V[['alt2']] = -(b_mu_rent + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
-                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*                          
-                            (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_nat_Age * Age_mean * Naturalness_2
-                            +  mu_nat_Income * Income_mean * Naturalness_2 + mu_walking_NR * Z_Mean_NR * WalkingDistance_2
-                            +  mu_walking_Age * Age_mean * WalkingDistance_2 + mu_walking_Income * Income_mean * WalkingDistance_2
-                            - Rent_2)
-  
-  V[['alt3']] = -(b_mu_rent + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
-                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*
-                            (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_ASC_Age * Age_mean + mu_ASC_Income * Income_mean 
-                            +  mu_nat_NR * Z_Mean_NR *Naturalness_3 + mu_nat_Age * Age_mean * Naturalness_3
-                            +  mu_nat_Income * Income_mean * Naturalness_3 + mu_walking_NR * Z_Mean_NR * WalkingDistance_3
-                            +  mu_walking_Age * Age_mean * WalkingDistance_3 + mu_walking_Income * Income_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_rentX = 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_rentX)
-
-
+#### 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_X",
+  modelDescr = "MXL wtp space Case D Interactions",
+  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_ASC_Age = 0,
+              mu_ASC_Income = 0,
+              mu_rent_treated = 0,
+              mu_rent_vol_treated = 0,
+              mu_rent_no_info = 0,
+              mu_rent_NR = 0,
+              mu_rent_Age = 0,
+              mu_rent_Income = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_NR = 0,
+              mu_nat_Age = 0,
+              mu_nat_Income = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_NR = 0,
+              mu_walking_Age = 0,
+              mu_walking_Income = 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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info 
+                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*
+                            (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_nat_Age * Age_mean * Naturalness_1
+                            +  mu_nat_Income * Income_mean * Naturalness_1 + mu_walking_NR * Z_Mean_NR * WalkingDistance_1
+                            +  mu_walking_Age * Age_mean * WalkingDistance_1 + mu_walking_Income * Income_mean * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
+                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*                          
+                            (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_nat_Age * Age_mean * Naturalness_2
+                            +  mu_nat_Income * Income_mean * Naturalness_2 + mu_walking_NR * Z_Mean_NR * WalkingDistance_2
+                            +  mu_walking_Age * Age_mean * WalkingDistance_2 + mu_walking_Income * Income_mean * WalkingDistance_2
+                            - Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
+                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*
+                            (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_ASC_Age * Age_mean + mu_ASC_Income * Income_mean 
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_3 + mu_nat_Age * Age_mean * Naturalness_3
+                            +  mu_nat_Income * Income_mean * Naturalness_3 + mu_walking_NR * Z_Mean_NR * WalkingDistance_3
+                            +  mu_walking_Age * Age_mean * WalkingDistance_3 + mu_walking_Income * Income_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_rentX = 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_rentX)
+
+