matching mxl models

Scripts/logit/chr_vol_treat.R

@@ -76,7 +76,7 @@ data<-select(data, Choice_Treat, id, Age, Q02W123, Uni_degree, QFIncome,Z_Mean_N
              Uni_degree, Q02W123,Q04W123,Q05W123, Q08W123,Q09W123,Q10W123,Q11W3,Q12W123,Q13W23,Q14S01W123,Q14S02W23,
              Q14S03W123,Q14S04W123,Q14S05W123,Q14S06W23,Q14S07W3,Q14S08W2,Q14W23,
              Q15S01W3,Q15S02W3,Q16W3,Q17W13,Q18W123,Q19W3,C02W23,Q20W23,Q21W23,
-             Q22S01W123,Q22S02W23,Q23W123,Q24S01W123,
+             Q22S01W123,Q22S02W23,UGS_visits,Q24S01W123,
              Q24S02W123,Q24S03W123,Q24S04W23,Q24S05W123,Q25W23,Q26S01W123,Q26S02W123,
              Q26S03W23,Q26S04W123,Q26S05W123,Q26S06W123,Q26S07W23,Q26S08W23,Q26S99W23,
              Q27W123,Q30W23,Q31S01W23,Q31S02W23,Q31S03W23,Q31S04W23,
@@ -177,4 +177,4 @@ best_coords <- coords(roc_obj, "best", best.method="youden")
 
 cut_off <- best_coords$threshold
 
->>>>>>> e410a7a5c52ecf34b454cff74f0b641cd5615679
+

Scripts/mxl/Prediction models/mxl_wtp_space_pred_matching.R

+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+data_predictions <- readRDS("Data/predictions.RDS")
+
+database <- left_join(database_full, data_predictions, by="id")
+
+
+
+database <- database %>% 
+  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)) %>% 
+  mutate(Dummy_Treated_Pred = case_when(Dummy_Treated == 1 & PredictedGroup == 1 ~1, TRUE~0),
+         Dummy_Treated_Not_Pred = case_when(Dummy_Treated == 1 & PredictedGroup == 0 ~1, TRUE~0))
+
+
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Prediction matching",
+  modelDescr = "MXL wtp space Prediction matching",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl/prediction"
+)
+
+##### 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_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_ASC_sq_treat_pred = 0,
+              mu_ASC_sq_treat_not_pred = 0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_treat_pred = 0,
+              mu_nat_treat_not_pred = 0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_treat_pred = 0,
+              mu_walking_treat_not_pred = 0,
+              mu_rent_vol_treated = 0,
+              mu_rent_no_info = 0,
+              mu_rent_treat_pred = 0,
+              mu_rent_treat_not_pred = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                    mu_rent_treat_not_pred * Dummy_Treated_Not_Pred)*
+                            (b_mu_natural*Naturalness_1 + b_mu_walking*WalkingDistance_1 
+                            + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_1 + mu_nat_no_info * Dummy_no_info * Naturalness_1
+                            + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_1 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_1
+                            + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_1 + mu_walking_no_info* Dummy_no_info * WalkingDistance_1
+                            + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_1 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] =  -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                     mu_rent_treat_not_pred * Dummy_Treated_Not_Pred)*
+                     (b_mu_natural*Naturalness_2 + b_mu_walking*WalkingDistance_2 
+                     + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_2 + mu_nat_no_info * Dummy_no_info * Naturalness_2
+                     + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_2 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_2
+                     + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_2 + mu_walking_no_info* Dummy_no_info * WalkingDistance_2
+                     + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_2 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_2
+                    - Rent_2)
+  
+  V[['alt3']] =  -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                     mu_rent_treat_not_pred * Dummy_Treated_Not_Pred)*
+                   (b_mu_natural*Naturalness_3 + b_mu_walking*WalkingDistance_3 
+                    + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_3 + mu_nat_no_info * Dummy_no_info * Naturalness_3
+                    + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_3 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_3
+                    + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_3 + mu_walking_no_info* Dummy_no_info * WalkingDistance_3
+                    + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_3 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_3
+                    + b_ASC_sq + mu_ASC_sq_vol_treated * Dummy_Vol_Treated + mu_ASC_sq_no_info * Dummy_no_info
+                    + mu_ASC_sq_treat_pred * Dummy_Treated_Pred + mu_ASC_sq_treat_not_pred * Dummy_Treated_Not_Pred- 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_e = apollo_estimate(apollo_beta, apollo_fixed,
+                                 apollo_probabilities, apollo_inputs, 
+                                 estimate_settings=list(maxIterations=400,
+                                                        estimationRoutine="bfgs",
+                                                        hessianRoutine="analytic"))
+
+
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ##
+# ################################################################# #
+apollo_saveOutput(mxl_wtp_case_e)
+
+

Scripts/mxl/Prediction models/mxl_wtp_space_pred_matching_all.R

+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+data_predictions <- readRDS("Data/predictions.RDS")
+
+database <- left_join(database_full, data_predictions, by="id")
+
+
+
+database <- database %>% 
+  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)) %>% 
+  mutate(Dummy_Treated_Pred = case_when(Dummy_Treated == 1 & PredictedGroup == 1 ~1, TRUE~0),
+         Dummy_Treated_Not_Pred = case_when(Dummy_Treated == 1 & PredictedGroup == 0 ~1, TRUE~0)) %>% 
+  mutate(Dummy_Control_Pred = case_when(Treatment_new == 6 & PredictedGroup == 1 ~1, TRUE~0))
+
+
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Prediction matching all",
+  modelDescr = "MXL wtp space Prediction matching all",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl/prediction"
+)
+
+##### 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_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_ASC_sq_treat_pred = 0,
+              mu_ASC_sq_treat_not_pred = 0,
+              mu_ASC_sq_control_pred = 0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_treat_pred = 0,
+              mu_nat_treat_not_pred = 0,
+              mu_nat_control_pred = 0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_treat_pred = 0,
+              mu_walking_treat_not_pred = 0,
+              mu_walking_control_pred = 0,
+              mu_rent_vol_treated = 0,
+              mu_rent_no_info = 0,
+              mu_rent_treat_pred = 0,
+              mu_rent_treat_not_pred = 0,
+              mu_rent_control_pred = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                    mu_rent_treat_not_pred * Dummy_Treated_Not_Pred + mu_rent_control_pred * Dummy_Control_Pred)*
+                            (b_mu_natural*Naturalness_1 + b_mu_walking*WalkingDistance_1 
+                            + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_1 + mu_nat_no_info * Dummy_no_info * Naturalness_1
+                            + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_1 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_1 + mu_nat_control_pred * Dummy_Control_Pred * Naturalness_1
+                            + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_1 + mu_walking_no_info* Dummy_no_info * WalkingDistance_1
+                            + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_1 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_1 + mu_walking_control_pred * Dummy_Control_Pred * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] =  -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                     mu_rent_treat_not_pred * Dummy_Treated_Not_Pred + mu_rent_control_pred * Dummy_Control_Pred)*
+                     (b_mu_natural*Naturalness_2 + b_mu_walking*WalkingDistance_2 
+                     + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_2 + mu_nat_no_info * Dummy_no_info * Naturalness_2
+                     + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_2 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_2 + mu_nat_control_pred * Dummy_Control_Pred * Naturalness_2
+                     + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_2 + mu_walking_no_info* Dummy_no_info * WalkingDistance_2
+                     + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_2 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_2 + mu_walking_control_pred * Dummy_Control_Pred * WalkingDistance_2
+                    - Rent_2)
+  
+  V[['alt3']] =  -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                     mu_rent_treat_not_pred * Dummy_Treated_Not_Pred + mu_rent_control_pred * Dummy_Control_Pred)*
+                   (b_mu_natural*Naturalness_3 + b_mu_walking*WalkingDistance_3 
+                    + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_3 + mu_nat_no_info * Dummy_no_info * Naturalness_3
+                    + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_3 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_3 + mu_nat_control_pred * Dummy_Control_Pred * Naturalness_3
+                    + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_3 + mu_walking_no_info* Dummy_no_info * WalkingDistance_3
+                    + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_3 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_3 + mu_walking_control_pred * Dummy_Control_Pred * WalkingDistance_3
+                    + b_ASC_sq + mu_ASC_sq_vol_treated * Dummy_Vol_Treated + mu_ASC_sq_no_info * Dummy_no_info
+                    + mu_ASC_sq_treat_pred * Dummy_Treated_Pred + mu_ASC_sq_treat_not_pred * Dummy_Treated_Not_Pred + mu_ASC_sq_control_pred * Dummy_Control_Pred - 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_matching_all = 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_matching_all)
+
+

Scripts/mxl/Prediction models/mxl_wtp_space_pred_matching_all_ref_control_pred.R

+#### Apollo standard script #####
+
+library(apollo) # Load apollo package 
+
+
+data_predictions <- readRDS("Data/predictions.RDS")
+
+database <- left_join(database_full, data_predictions, by="id")
+
+
+
+database <- database %>% 
+  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)) %>% 
+  mutate(Dummy_Treated_Pred = case_when(Dummy_Treated == 1 & PredictedGroup == 1 ~1, TRUE~0),
+         Dummy_Treated_Not_Pred = case_when(Dummy_Treated == 1 & PredictedGroup == 0 ~1, TRUE~0)) %>% 
+  mutate(Dummy_Control_Not_Pred = case_when(Treatment_new == 6 & PredictedGroup == 0 ~1, TRUE~0))
+
+
+
+#initialize model 
+
+apollo_initialise()
+
+
+### Set core controls
+apollo_control = list(
+  modelName  = "MXL_wtp_Prediction matching all cp",
+  modelDescr = "MXL wtp space Prediction matching all cp",
+  indivID    ="id",
+  mixing     = TRUE,
+  HB= FALSE,
+  nCores     = n_cores, 
+  outputDirectory = "Estimation_results/mxl/prediction"
+)
+
+##### 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_vol_treated = 0,
+              mu_ASC_sq_no_info = 0,
+              mu_ASC_sq_treat_pred = 0,
+              mu_ASC_sq_treat_not_pred = 0,
+              mu_ASC_sq_control_not_pred = 0,
+              mu_nat_vol_treated = 0,
+              mu_nat_no_info = 0,
+              mu_nat_treat_pred = 0,
+              mu_nat_treat_not_pred = 0,
+              mu_nat_control_not_pred = 0,
+              mu_walking_vol_treated = 0,
+              mu_walking_no_info = 0,
+              mu_walking_treat_pred = 0,
+              mu_walking_treat_not_pred = 0,
+              mu_walking_control_not_pred = 0,
+              mu_rent_vol_treated = 0,
+              mu_rent_no_info = 0,
+              mu_rent_treat_pred = 0,
+              mu_rent_treat_not_pred = 0,
+              mu_rent_control_not_pred = 0,
+              sig_natural = 15,
+              sig_walking = 2,
+              sig_rent = 2,
+              sig_ASC_sq = 2)
+
+### specify parameters that should be kept fixed, here = none
+apollo_fixed = c()
+
+### Set parameters for generating draws, use 2000 sobol draws
+apollo_draws = list(
+  interDrawsType = "sobol",
+  interNDraws    = n_draws,
+  interUnifDraws = c(),
+  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
+  intraDrawsType = "halton",
+  intraNDraws    = 0,
+  intraUnifDraws = c(),
+  intraNormDraws = c()
+)
+
+### Create random parameters, define distribution of the parameters
+apollo_randCoeff = function(apollo_beta, apollo_inputs){
+  randcoeff = list()
+  
+  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
+  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
+  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
+  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
+  
+  return(randcoeff)
+}
+
+
+### validate 
+apollo_inputs = apollo_validateInputs()
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### Function initialisation: do not change the following three commands
+  ### Attach inputs and detach after function exit
+  apollo_attach(apollo_beta, apollo_inputs)
+  on.exit(apollo_detach(apollo_beta, apollo_inputs))
+  
+  ### Create list of probabilities P
+  P = list()
+  
+  #### List of utilities (later integrated in mnl_settings below)  ####
+  # Define utility functions here:
+  
+  V = list()
+  V[['alt1']] = -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                    mu_rent_treat_not_pred * Dummy_Treated_Not_Pred + mu_rent_control_not_pred * Dummy_Control_Not_Pred)*
+                            (b_mu_natural*Naturalness_1 + b_mu_walking*WalkingDistance_1 
+                            + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_1 + mu_nat_no_info * Dummy_no_info * Naturalness_1
+                            + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_1 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_1 + mu_nat_control_not_pred * Dummy_Control_Not_Pred * Naturalness_1
+                            + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_1 + mu_walking_no_info* Dummy_no_info * WalkingDistance_1
+                            + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_1 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_1 + mu_walking_control_not_pred * Dummy_Control_Not_Pred * WalkingDistance_1
+                            - Rent_1)
+  
+  V[['alt2']] =  -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                     mu_rent_treat_not_pred * Dummy_Treated_Not_Pred + mu_rent_control_not_pred * Dummy_Control_Not_Pred)*
+                     (b_mu_natural*Naturalness_2 + b_mu_walking*WalkingDistance_2 
+                     + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_2 + mu_nat_no_info * Dummy_no_info * Naturalness_2
+                     + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_2 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_2 + mu_nat_control_not_pred * Dummy_Control_Not_Pred * Naturalness_2
+                     + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_2 + mu_walking_no_info* Dummy_no_info * WalkingDistance_2
+                     + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_2 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_2 + mu_walking_control_not_pred * Dummy_Control_Not_Pred * WalkingDistance_2
+                    - Rent_2)
+  
+  V[['alt3']] =  -(b_mu_rent + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info + mu_rent_treat_pred * Dummy_Treated_Pred +
+                     mu_rent_treat_not_pred * Dummy_Treated_Not_Pred + mu_rent_control_not_pred * Dummy_Control_Not_Pred)*
+                   (b_mu_natural*Naturalness_3 + b_mu_walking*WalkingDistance_3 
+                    + mu_nat_vol_treated * Dummy_Vol_Treated * Naturalness_3 + mu_nat_no_info * Dummy_no_info * Naturalness_3
+                    + mu_nat_treat_pred * Dummy_Treated_Pred * Naturalness_3 + mu_nat_treat_not_pred * Dummy_Treated_Not_Pred * Naturalness_3 + mu_nat_control_not_pred * Dummy_Control_Not_Pred * Naturalness_3
+                    + mu_walking_vol_treated * Dummy_Vol_Treated * WalkingDistance_3 + mu_walking_no_info* Dummy_no_info * WalkingDistance_3
+                    + mu_walking_treat_pred * Dummy_Treated_Pred * WalkingDistance_3 + mu_walking_treat_not_pred * Dummy_Treated_Not_Pred * WalkingDistance_3 + mu_walking_control_not_pred * Dummy_Control_Not_Pred * WalkingDistance_3
+                    + b_ASC_sq + mu_ASC_sq_vol_treated * Dummy_Vol_Treated + mu_ASC_sq_no_info * Dummy_no_info
+                    + mu_ASC_sq_treat_pred * Dummy_Treated_Pred + mu_ASC_sq_treat_not_pred * Dummy_Treated_Not_Pred + mu_ASC_sq_control_not_pred * Dummy_Control_Not_Pred - 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_matching_all_cp = 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_matching_all_cp)
+
+