diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_lognormRent_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_lognormRent_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..a992b8ee23fc6628062a9e8ecf375279038dab9d
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_lognormRent_scale.R
@@ -0,0 +1,153 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A lognorm Rent Scale",
+    modelDescr = "MXL_pref Case A lognorm Rent Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT = 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_rent_T * Rent_1 * Dummy_Treated +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated +
+                                mu_rent_VT * Rent_1 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                               mu_rent_T * Rent_2 * Dummy_Treated +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                                mu_rent_VT * Rent_2 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_rent_T * Rent_3 * Dummy_Treated + mu_asc_VT * Dummy_Vol_Treated +
+                                mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated +
+                                mu_rent_VT * Rent_3 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_lognorm_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_lognorm_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_normRent_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_normRent_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..ed63cc7ec293593043abd7d0936fe9e7bbd18ce4
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_normRent_scale.R
@@ -0,0 +1,153 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A norm Rent Scale",
+    modelDescr = "MXL_pref Case A norm Rent Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT = 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = mu_rent + sig_rent * draws_rent
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_rent_T * Rent_1 * Dummy_Treated +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated +
+                                mu_rent_VT * Rent_1 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                               mu_rent_T * Rent_2 * Dummy_Treated +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                                mu_rent_VT * Rent_2 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_rent_T * Rent_3 * Dummy_Treated + mu_asc_VT * Dummy_Vol_Treated +
+                                mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated +
+                                mu_rent_VT * Rent_3 * Dummy_Vol_Treated + mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_norm_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_norm_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_only_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_only_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..477de7cb5c435eebd92ea763ba355b80321faa77
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_only_scale.R
@@ -0,0 +1,136 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A only  Scale",
+    modelDescr = "MXL_pref Case A only  Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + 
+                                                                                   b_mu_walking * WalkingDistance_1) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + 
+                                                                                   b_mu_walking * WalkingDistance_2) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_mu_natural * Naturalness_3 + 
+                                                                                   b_mu_walking * WalkingDistance_3 + b_ASC_sq) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_only_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_only_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseA_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..16c49729bedc1b761810dcaa196bf0acc091e6a2
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseA_scale.R
@@ -0,0 +1,149 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_pref Case A  Scale",
+    modelDescr = "MXL_pref Case A  Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl/prefspace"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT = 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_1 + b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated +
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated) 
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_2 + b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                                mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated) 
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(b_mu_rent*Rent_3 + b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT * Dummy_Vol_Treated +
+                                mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated +
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated +
+                                mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated) 
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_pref_case_a_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_pref_case_a_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_pref_space_caseD_only_scale.R b/Scripts/mxl/Scale_models/mxl_pref_space_caseD_only_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..9db5aae3ad5c92bc134d21e60ac3ea20985a72cd
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_pref_space_caseD_only_scale.R
@@ -0,0 +1,150 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D scale",
+  modelDescr = "MXL wtp space Case D scale",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2,
+              lambda_T = 0,
+              lambda_VT =0,
+              lambda_NO = 0)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated + lambda_NO * Dummy_no_info)*
+                            (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                             b_mu_rent * Rent_1)
+  
+  V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated + lambda_NO * Dummy_no_info)*
+                            (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                             b_mu_rent * Rent_2)
+  
+  V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated + lambda_NO * Dummy_no_info)*
+                            (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 +
+                             b_mu_rent * Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_d_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_d_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_base_scale.R b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_base_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..c6ded1ef21d2e6bde1d10e7f1452e14d1575f082
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_base_scale.R
@@ -0,0 +1,149 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Case A Scale",
+    modelDescr = "MXL_wtp Case A Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT= 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*-(b_mu_rent)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated+
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated
+                              - Rent_1)
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*-(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                              mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated
+                              - Rent_2)
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*-(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT *  Dummy_Vol_Treated +mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated + 
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated
+                                - Rent_3)
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_case_a_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_case_a_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_exp_scale.R b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_exp_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..74574cb8652e858085e980d6cec7823731c889c9
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_exp_scale.R
@@ -0,0 +1,151 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Case A Rent Int Exp Scale",
+    modelDescr = "MXL_wtp Case A Rent Int Exp Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT= 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + exp(mu_rent_T)*Dummy_Treated + exp(mu_rent_VT)*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated+
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated
+                              - Rent_1))
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + exp(mu_rent_T)*Dummy_Treated + exp(mu_rent_VT)*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                              mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated
+                              - Rent_2))
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + exp(mu_rent_T)*Dummy_Treated + exp(mu_rent_VT)*Dummy_Vol_Treated)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT *  Dummy_Vol_Treated +mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated + 
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated
+                                - Rent_3))
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_case_a_rentINT_exp_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_case_a_rentINT_exp_scale)
+
+
diff --git a/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_scale.R b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_scale.R
new file mode 100644
index 0000000000000000000000000000000000000000..e11aaf4cfaa2ddbb4f4031005f939a9cf78e6ba4
--- /dev/null
+++ b/Scripts/mxl/Scale_models/mxl_wtp_space_caseA_rentINT_scale.R
@@ -0,0 +1,151 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
+  mutate(Dummy_Treated = case_when(Treatment_A == "Treated" ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_A == "Vol_Treated" ~ 1, TRUE ~ 0))
+
+  #initialize model 
+  
+  apollo_initialise()
+  
+  
+  ### Set core controls
+  apollo_control = list(
+    modelName  = "MXL_wtp Case A Rent Int Scale",
+    modelDescr = "MXL_wtp Case A Rent Int Scale",
+    indivID    ="id",
+    mixing     = TRUE,
+    HB= FALSE,
+    nCores     = n_cores, 
+    outputDirectory = "Estimation_results/mxl"
+  )
+  
+  ##### Define model parameters depending on your attributes and model specification! ####
+  # set values to 0 for conditional logit model
+  
+  apollo_beta=c(mu_natural = 15,
+                mu_walking = -1,
+                mu_rent = -2,
+                ASC_sq = 0,
+                mu_rent_T = 0,
+                mu_nat_T = 0,
+                mu_wd_T= 0,
+                mu_asc_T = 0,
+                mu_rent_VT = 0,
+                mu_nat_VT = 0,
+                mu_wd_VT= 0,
+                mu_asc_VT = 0,
+                sig_natural = 15,
+                sig_walking = 2,
+                sig_rent = 2,
+                sig_ASC_sq = 2,
+                lambda_T = 0,
+                lambda_VT =0)
+  
+  ### specify parameters that should be kept fixed, here = none
+  apollo_fixed = c()
+  
+  ### Set parameters for generating draws, use 2000 sobol draws
+  apollo_draws = list(
+    interDrawsType = "sobol",
+    interNDraws    = n_draws,
+    interUnifDraws = c(),
+    interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+    intraDrawsType = "halton",
+    intraNDraws    = 0,
+    intraUnifDraws = c(),
+    intraNormDraws = c()
+  )
+  
+  ### Create random parameters, define distribution of the parameters
+  apollo_randCoeff = function(apollo_beta, apollo_inputs){
+    randcoeff = list()
+    
+    randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+    randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+    randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+    randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+    
+    return(randcoeff)
+  }
+  
+  
+  ### validate 
+  apollo_inputs = apollo_validateInputs()
+  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+    
+    ### Function initialisation: do not change the following three commands
+    ### Attach inputs and detach after function exit
+    apollo_attach(apollo_beta, apollo_inputs)
+    on.exit(apollo_detach(apollo_beta, apollo_inputs))
+    
+    ### Create list of probabilities P
+    P = list()
+    
+    #### List of utilities (later integrated in mnl_settings below)  ####
+    # Define utility functions here:
+     
+    V = list()
+    V[['alt1']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + mu_rent_T*Dummy_Treated + mu_rent_VT*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                                mu_nat_T * Naturalness_1 * Dummy_Treated + mu_wd_T * WalkingDistance_1 * Dummy_Treated+
+                                mu_nat_VT * Naturalness_1 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_1 * Dummy_Vol_Treated
+                              - Rent_1))
+    
+    V[['alt2']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + mu_rent_T*Dummy_Treated + mu_rent_VT*Dummy_Vol_Treated)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
+                                mu_nat_T * Naturalness_2 * Dummy_Treated + mu_wd_T * WalkingDistance_2 * Dummy_Treated +
+                              mu_nat_VT * Naturalness_2 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_2 * Dummy_Vol_Treated
+                              - Rent_2))
+    
+    V[['alt3']] = exp(lambda_T * Dummy_Treated + lambda_VT * Dummy_Vol_Treated)*(-(b_mu_rent + mu_rent_T*Dummy_Treated + mu_rent_VT*Dummy_Vol_Treated)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_T *  Dummy_Treated +
+                                mu_asc_VT *  Dummy_Vol_Treated +mu_nat_T * Naturalness_3 * Dummy_Treated + mu_wd_T * WalkingDistance_3 * Dummy_Treated + 
+                                mu_nat_VT * Naturalness_3 * Dummy_Vol_Treated + mu_wd_VT * WalkingDistance_3 * Dummy_Vol_Treated
+                                - Rent_3))
+    
+    
+    ### Define settings for MNL model component
+    mnl_settings = list(
+      alternatives  = c(alt1=1, alt2=2, alt3=3),
+      avail         = 1, # all alternatives are available in every choice
+      choiceVar     = choice,
+      V             = V#,  # tell function to use list vector defined above
+      
+    )
+    
+    ### Compute probabilities using MNL model
+    P[['model']] = apollo_mnl(mnl_settings, functionality)
+    
+    ### Take product across observation for same individual
+    P = apollo_panelProd(P, apollo_inputs, functionality)
+    
+    ### Average across inter-individual draws - nur bei Mixed Logit!
+    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+    
+    ### Prepare and return outputs of function
+    P = apollo_prepareProb(P, apollo_inputs, functionality)
+    return(P)
+  }
+  
+  
+  
+  # ################################################################# #
+  #### MODEL ESTIMATION                                            ##
+  # ################################################################# #
+  # estimate model with bfgs algorithm
+  
+  mxl_wtp_case_a_rentINT_scale = apollo_estimate(apollo_beta, apollo_fixed,
+                        apollo_probabilities, apollo_inputs, 
+                        estimate_settings=list(maxIterations=400,
+                                               estimationRoutine="bfgs",
+                                               hessianRoutine="analytic"))
+  
+  
+  
+  # ################################################################# #
+  #### MODEL OUTPUTS                                               ##
+  # ################################################################# #
+  apollo_saveOutput(mxl_wtp_case_a_rentINT_scale)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_caseD.R b/Scripts/mxl/mxl_wtp_space_caseD.R
new file mode 100644
index 0000000000000000000000000000000000000000..d51b52d3acd20fe9293fb056b77e23ef123b495c
--- /dev/null
+++ b/Scripts/mxl/mxl_wtp_space_caseD.R
@@ -0,0 +1,164 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D base",
+  modelDescr = "MXL wtp space Case D base",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
+                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
+                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated- Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
+                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                            +  mu_ASC_sq_no_info * Dummy_no_info - Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_c = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_c)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_NR.R b/Scripts/mxl/mxl_wtp_space_caseD_NR.R
new file mode 100644
index 0000000000000000000000000000000000000000..f14a8bfda1014b46ebde6859435c3e8b77fd33cf
--- /dev/null
+++ b/Scripts/mxl/mxl_wtp_space_caseD_NR.R
@@ -0,0 +1,180 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D_NR base",
+  modelDescr = "MXL wtp space Case D NR base",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_ASC_NR = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_NR = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_NR = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*
+                            (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
+                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_1 +
+                            + mu_walking_NR * Z_Mean_NR * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent)*                          
+                            (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
+                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_2 + 
+                            + mu_walking_NR * Z_Mean_NR * WalkingDistance_2
+                            - Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent)*
+                            (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
+                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                            +  mu_ASC_sq_no_info * Dummy_no_info 
+                            +  mu_ASC_NR * Z_Mean_NR 
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_3 
+                            +  mu_walking_NR * Z_Mean_NR * WalkingDistance_3
+                            - Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_d_NR = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_d_NR)
+
+
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_age.R b/Scripts/mxl/mxl_wtp_space_caseD_age.R
new file mode 100644
index 0000000000000000000000000000000000000000..9409124d1047a9d7704cae84807a9e278d289e92
--- /dev/null
+++ b/Scripts/mxl/mxl_wtp_space_caseD_age.R
@@ -0,0 +1,180 @@
+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+
+# Test treatment effect
+
+database <- database_full %>%
+  filter(!is.na(Treatment_new)) %>%
+  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
+         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
+         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))
+
+table(database$Dummy_Treated)
+table(database$Dummy_Vol_Treated)
+table(database$Dummy_no_info)
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Case_D_Age",
+  modelDescr = "MXL wtp space Case D Age",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl"
+)
+
+##### Define model parameters depending on your attributes and model specification! ####
+# set values to 0 for conditional logit model
+
+apollo_beta=c(mu_natural = 15,
+              mu_walking = -1,
+              mu_rent = -2,
+              ASC_sq = 0,
+              mu_ASC_sq_treated = 0,
+              mu_ASC_sq_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_asc_age = 0,
+              mu_nat_treated =0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_age = 0,
+              mu_walking_treated =0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_age = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent)*
+                            (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
+                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated
+                            +  mu_nat_age * Age_mean *Naturalness_1 +
+                            + mu_walking_age * Age_mean * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] = -(b_mu_rent)*                          
+                            (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
+                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
+                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
+                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated
+                            +  mu_nat_age * Age_mean *Naturalness_2 + 
+                            + mu_walking_age * Age_mean * WalkingDistance_2
+                            - Rent_2)
+  
+  V[['alt3']] = -(b_mu_rent)*
+                            (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
+                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
+                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
+                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
+                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
+                            +  mu_ASC_sq_no_info * Dummy_no_info 
+                            +  mu_asc_age * Age_mean
+                            +  mu_nat_age * Age_mean *Naturalness_3 
+                            +  mu_walking_age * Age_mean * WalkingDistance_3
+                            - Rent_3)
+  
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3),
+    avail         = 1, # all alternatives are available in every choice
+    choiceVar     = choice,
+    V             = V#,  # tell function to use list vector defined above
+    
+  )
+  
+  ### Compute probabilities using MNL model
+  P[['model']] = apollo_mnl(mnl_settings, functionality)
+  
+  ### Take product across observation for same individual
+  P = apollo_panelProd(P, apollo_inputs, functionality)
+  
+  ### Average across inter-individual draws - nur bei Mixed Logit!
+  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ##
+# ################################################################# #
+# estimate model with bfgs algorithm
+
+mxl_wtp_case_d_age = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_d_age)
+
+
diff --git a/project_start.qmd b/project_start.qmd
index 776eb1b02e8b77e3573ac08e88cc582b3d88201f..aebd3fa49f4b97931a54f8dbd2cebbb0758958a1 100644
--- a/project_start.qmd
+++ b/project_start.qmd
@@ -155,7 +155,7 @@ list_ols <- list("(Intercept)" = "Intercept", "as.factor(Treatment_A)Treated" =
 
 +   **Socio-demographics**: Age, Gender, Income, Education
 
-+   **Attitudinal variable**: Measure derived from 21 items on **nature relatedness** [@nisbet2009nature]
++   **Attitudinal variable**: Measure derived from 21 items on **nature relatedness** (NR-Index) [@nisbet2009nature]
 
 :::
 
diff --git a/ugs_data-main.Rproj b/ugs_data-main.Rproj
index 3af27f6a64e3e141865e25d4c4e719c13f9ed422..8e3c2ebc99e2e337f7d69948b93529a437590b27 100644
--- a/ugs_data-main.Rproj
+++ b/ugs_data-main.Rproj
@@ -1,13 +1,13 @@
-Version: 1.0
-
-RestoreWorkspace: Default
-SaveWorkspace: Default
-AlwaysSaveHistory: Default
-
-EnableCodeIndexing: Yes
-UseSpacesForTab: Yes
-NumSpacesForTab: 2
-Encoding: UTF-8
-
-RnwWeave: Sweave
-LaTeX: pdfLaTeX
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX