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)
+
+