#### Apollo standard script #####

library(apollo) # Load apollo package 



# Test treatment effect

database <- database_full %>%
  filter(!is.na(Treatment_new)) %>%
  mutate(Dummy_Treated = case_when(Treatment_new == 1|Treatment_new == 2  ~ 1, TRUE ~ 0),
         Dummy_Vol_Treated = case_when(Treatment_new == 5 |Treatment_new == 4 ~ 1, TRUE ~ 0),
         Dummy_no_info = case_when(Treatment_new == 3 ~ 1, TRUE~0))

table(database$Dummy_Treated)
table(database$Dummy_Vol_Treated)
table(database$Dummy_no_info)

#initialize model 

apollo_initialise()


### Set core controls
apollo_control = list(
  modelName  = "MXL_wtp_Case_D_X",
  modelDescr = "MXL wtp space Case D Interactions",
  indivID    ="id",
  mixing     = TRUE,
  HB= FALSE,
  nCores     = n_cores, 
  outputDirectory = "Estimation_results/mxl"
)

##### Define model parameters depending on your attributes and model specification! ####
# set values to 0 for conditional logit model

apollo_beta=c(mu_natural = 15,
              mu_walking = -1,
              mu_rent = -2,
              ASC_sq = 0,
              mu_ASC_sq_treated = 0,
              mu_ASC_sq_vol_treated = 0,
              mu_ASC_sq_no_info = 0,
              mu_ASC_NR = 0,
              mu_ASC_Age = 0,
              mu_ASC_Income = 0,
              mu_rent_treated = 0,
              mu_rent_vol_treated = 0,
              mu_rent_no_info = 0,
              mu_rent_NR = 0,
              mu_rent_Age = 0,
              mu_rent_Income = 0,
              mu_nat_treated =0,
              mu_nat_vol_treated = 0,
              mu_nat_no_info = 0,
              mu_nat_NR = 0,
              mu_nat_Age = 0,
              mu_nat_Income = 0,
              mu_walking_treated =0,
              mu_walking_vol_treated = 0,
              mu_walking_no_info = 0,
              mu_walking_NR = 0,
              mu_walking_Age = 0,
              mu_walking_Income = 0,
              sig_natural = 15,
              sig_walking = 2,
              sig_rent = 2,
              sig_ASC_sq = 2)

### specify parameters that should be kept fixed, here = none
apollo_fixed = c()

### Set parameters for generating draws, use 2000 sobol draws
apollo_draws = list(
  interDrawsType = "sobol",
  interNDraws    = n_draws,
  interUnifDraws = c(),
  interNormDraws = c("draws_natural", "draws_walking", "draws_rent", "draws_asc"),
  intraDrawsType = "halton",
  intraNDraws    = 0,
  intraUnifDraws = c(),
  intraNormDraws = c()
)

### Create random parameters, define distribution of the parameters
apollo_randCoeff = function(apollo_beta, apollo_inputs){
  randcoeff = list()
  
  randcoeff[["b_mu_natural"]] = mu_natural + sig_natural * draws_natural
  randcoeff[["b_mu_walking"]] = mu_walking + sig_walking * draws_walking
  randcoeff[["b_mu_rent"]] = -exp(mu_rent + sig_rent * draws_rent)
  randcoeff[["b_ASC_sq"]] = ASC_sq + sig_ASC_sq * draws_asc
  
  return(randcoeff)
}


### validate 
apollo_inputs = apollo_validateInputs()
apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
  
  ### Function initialisation: do not change the following three commands
  ### Attach inputs and detach after function exit
  apollo_attach(apollo_beta, apollo_inputs)
  on.exit(apollo_detach(apollo_beta, apollo_inputs))
  
  ### Create list of probabilities P
  P = list()
  
  #### List of utilities (later integrated in mnl_settings below)  ####
  # Define utility functions here:
  
  V = list()
  V[['alt1']] = -(b_mu_rent + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info 
                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*
                            (b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
                            +  mu_nat_treated * Naturalness_1 *Dummy_Treated + mu_nat_no_info * Naturalness_1 * Dummy_no_info
                            +  mu_nat_vol_treated * Naturalness_1 * Dummy_Vol_Treated
                            +  mu_walking_treated * WalkingDistance_1 *Dummy_Treated + mu_walking_no_info * WalkingDistance_1 * Dummy_no_info
                            +  mu_walking_vol_treated * WalkingDistance_1 * Dummy_Vol_Treated
                            +  mu_nat_NR * Z_Mean_NR *Naturalness_1 + mu_nat_Age * Age_mean * Naturalness_1
                            +  mu_nat_Income * Income_mean * Naturalness_1 + mu_walking_NR * Z_Mean_NR * WalkingDistance_1
                            +  mu_walking_Age * Age_mean * WalkingDistance_1 + mu_walking_Income * Income_mean * WalkingDistance_1
                            - Rent_1)
  
  V[['alt2']] = -(b_mu_rent + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*                          
                            (b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
                            + mu_nat_treated * Naturalness_2 *Dummy_Treated + mu_nat_no_info * Naturalness_2 * Dummy_no_info
                            + mu_nat_vol_treated * Naturalness_2 * Dummy_Vol_Treated
                            + mu_walking_treated * WalkingDistance_2 *Dummy_Treated + mu_walking_no_info * WalkingDistance_2 * Dummy_no_info
                            + mu_walking_vol_treated * WalkingDistance_2 * Dummy_Vol_Treated
                            +  mu_nat_NR * Z_Mean_NR *Naturalness_2 + mu_nat_Age * Age_mean * Naturalness_2
                            +  mu_nat_Income * Income_mean * Naturalness_2 + mu_walking_NR * Z_Mean_NR * WalkingDistance_2
                            +  mu_walking_Age * Age_mean * WalkingDistance_2 + mu_walking_Income * Income_mean * WalkingDistance_2
                            - Rent_2)
  
  V[['alt3']] = -(b_mu_rent + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info
                  + mu_rent_NR * Z_Mean_NR + mu_rent_Age * Age_mean + mu_rent_Income * Income_mean)*
                            (b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
                            +  mu_nat_treated * Naturalness_3 *Dummy_Treated + mu_nat_no_info * Naturalness_3 * Dummy_no_info
                            +  mu_nat_vol_treated * Naturalness_3 * Dummy_Vol_Treated 
                            +  mu_walking_treated * WalkingDistance_3 *Dummy_Treated + mu_walking_no_info * WalkingDistance_3 * Dummy_no_info
                            +  mu_walking_vol_treated * WalkingDistance_3 * Dummy_Vol_Treated
                            +  mu_ASC_sq_treated * Dummy_Treated + mu_ASC_sq_vol_treated * Dummy_Vol_Treated
                            +  mu_ASC_sq_no_info * Dummy_no_info 
                            +  mu_ASC_NR * Z_Mean_NR + mu_ASC_Age * Age_mean + mu_ASC_Income * Income_mean 
                            +  mu_nat_NR * Z_Mean_NR *Naturalness_3 + mu_nat_Age * Age_mean * Naturalness_3
                            +  mu_nat_Income * Income_mean * Naturalness_3 + mu_walking_NR * Z_Mean_NR * WalkingDistance_3
                            +  mu_walking_Age * Age_mean * WalkingDistance_3 + mu_walking_Income * Income_mean * WalkingDistance_3
                            - Rent_3)
  
  
  ### Define settings for MNL model component
  mnl_settings = list(
    alternatives  = c(alt1=1, alt2=2, alt3=3),
    avail         = 1, # all alternatives are available in every choice
    choiceVar     = choice,
    V             = V#,  # tell function to use list vector defined above
    
  )
  
  ### Compute probabilities using MNL model
  P[['model']] = apollo_mnl(mnl_settings, functionality)
  
  ### Take product across observation for same individual
  P = apollo_panelProd(P, apollo_inputs, functionality)
  
  ### Average across inter-individual draws - nur bei Mixed Logit!
  P = apollo_avgInterDraws(P, apollo_inputs, functionality)
  
  ### Prepare and return outputs of function
  P = apollo_prepareProb(P, apollo_inputs, functionality)
  return(P)
}



# ################################################################# #
#### MODEL ESTIMATION                                            ##
# ################################################################# #
# estimate model with bfgs algorithm

mxl_wtp_case_d_rentX = apollo_estimate(apollo_beta, apollo_fixed,
                                 apollo_probabilities, apollo_inputs, 
                                 estimate_settings=list(maxIterations=400,
                                                        estimationRoutine="bfgs",
                                                        hessianRoutine="analytic"))



# ################################################################# #
#### MODEL OUTPUTS                                               ##
# ################################################################# #
apollo_saveOutput(mxl_wtp_case_d_rentX)