diff --git a/Scripts/mxl/mxl_wtp_space_caseD.R b/Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R
similarity index 71%
rename from Scripts/mxl/mxl_wtp_space_caseD.R
rename to Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R
index d32875f8946cbe36e65e3d2496362ca43693977c..612c5147efe4a31b13964b4ad4d2bb2687a67d70 100644
--- a/Scripts/mxl/mxl_wtp_space_caseD.R
+++ b/Scripts/mxl/mxl_wtp_space_caseD_RentINT_X.R
@@ -23,8 +23,8 @@ apollo_initialise()
 
 ### Set core controls
 apollo_control = list(
-  modelName  = "MXL_wtp_Case_D",
-  modelDescr = "MXL wtp space Case D",
+  modelName  = "MXL_wtp_Case_D_X",
+  modelDescr = "MXL wtp space Case D Interactions",
   indivID    ="id",
   mixing     = TRUE,
   HB= FALSE,
@@ -42,15 +42,27 @@ apollo_beta=c(mu_natural = 15,
               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,
@@ -100,25 +112,44 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
   # 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)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
+  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 * QFIncome)*
+                            (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)
+                            +  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 * QFIncome * Naturalness_1 + mu_walking_NR * Z_Mean_NR * WalkingDistance_1
+                            +  mu_walking_Age * Age_mean * WalkingDistance_1 + mu_walking_Income * QFIncome * 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)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 
+  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 * QFIncome)*                          
+                            (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)
+                            + 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 * QFIncome * Naturalness_2 + mu_walking_NR * Z_Mean_NR * WalkingDistance_2
+                            +  mu_walking_Age * Age_mean * WalkingDistance_2 + mu_walking_Income * QFIncome * 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)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3  
+  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 * QFIncome)*
+                            (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)
+                            +  mu_ASC_sq_no_info * Dummy_no_info 
+                            +  mu_ASC_NR * Z_Mean_NR + mu_ASC_Age * Age_mean + mu_ASC_Income * QFIncome 
+                            +  mu_nat_NR * Z_Mean_NR *Naturalness_3 + mu_nat_Age * Age_mean * Naturalness_3
+                            +  mu_nat_Income * QFIncome * Naturalness_3 + mu_walking_NR * Z_Mean_NR * WalkingDistance_3
+                            +  mu_walking_Age * Age_mean * WalkingDistance_3 + mu_walking_Income * QFIncome * WalkingDistance_3
+                            - Rent_3)
   
   
   ### Define settings for MNL model component
@@ -151,7 +182,7 @@ apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimat
 # ################################################################# #
 # estimate model with bfgs algorithm
 
-mxl_wtp_case_c = apollo_estimate(apollo_beta, apollo_fixed,
+mxl_wtp_case_d_rentX = apollo_estimate(apollo_beta, apollo_fixed,
                                  apollo_probabilities, apollo_inputs, 
                                  estimate_settings=list(maxIterations=400,
                                                         estimationRoutine="bfgs",
@@ -162,6 +193,6 @@ mxl_wtp_case_c = apollo_estimate(apollo_beta, apollo_fixed,
 # ################################################################# #
 #### MODEL OUTPUTS                                               ##
 # ################################################################# #
-apollo_saveOutput(mxl_wtp_case_c)
+apollo_saveOutput(mxl_wtp_case_d_rentX)
 
 
diff --git a/Scripts/mxl/mxl_wtp_space_caseD_rentINT.R b/Scripts/mxl/mxl_wtp_space_caseD_rentINT.R
index 7d5699c85140299317eaf49291ecf34c45a50c92..d32875f8946cbe36e65e3d2496362ca43693977c 100644
--- a/Scripts/mxl/mxl_wtp_space_caseD_rentINT.R
+++ b/Scripts/mxl/mxl_wtp_space_caseD_rentINT.R
@@ -3,147 +3,165 @@
 library(apollo) # Load apollo package 
 
 
-database <- database_full %>% filter(!is.na(Treatment_A)) %>% 
-  mutate(Dummy_Choice_Info = case_when(Treatment_new == 1 ~ 1, Treatment_new == 4 ~ 1, Treatment_new == 5 ~ 1, TRUE ~ 0),
-         Dummy_Choice_No_Info = case_when(Treatment_new == 2 ~ 1, Treatment_new == 3 ~ 1, TRUE ~ 0))
 
-  #initialize model 
-  
-  apollo_initialise()
-  
-  
-  ### Set core controls
-  apollo_control = list(
-    modelName  = "MXL_wtp Case D Rent Int",
-    modelDescr = "MXL_wtp Case D Rent Int",
-    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_I = 0,
-                mu_nat_I = 0,
-                mu_wd_I= 0,
-                mu_asc_I = 0,
-                mu_rent_NI = 0,
-                mu_nat_NI = 0,
-                mu_wd_NI= 0,
-                mu_asc_NI = 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)
-  }
+# 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",
+  modelDescr = "MXL wtp space Case D",
+  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_rent_treated = 0,
+              mu_rent_vol_treated = 0,
+              mu_rent_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
   
-  ### 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_I*Dummy_Choice_Info + mu_rent_NI*Dummy_Choice_No_Info)*(b_mu_natural * Naturalness_1 + b_mu_walking * WalkingDistance_1 +
-                                mu_nat_I * Naturalness_1 * Dummy_Choice_Info + mu_wd_I * WalkingDistance_1 * Dummy_Choice_Info+
-                                mu_nat_NI * Naturalness_1 * Dummy_Choice_No_Info + mu_wd_NI * WalkingDistance_1 * Dummy_Choice_No_Info
-                              - Rent_1)
-    
-    V[['alt2']] = -(b_mu_rent + mu_rent_I*Dummy_Choice_Info + mu_rent_NI*Dummy_Choice_No_Info)*(b_mu_natural * Naturalness_2 + b_mu_walking * WalkingDistance_2 +
-                                mu_nat_I * Naturalness_2 * Dummy_Choice_Info + mu_wd_I * WalkingDistance_2 * Dummy_Choice_Info +
-                              mu_nat_NI * Naturalness_2 * Dummy_Choice_No_Info + mu_wd_NI * WalkingDistance_2 * Dummy_Choice_No_Info
-                              - Rent_2)
-    
-    V[['alt3']] = -(b_mu_rent + mu_rent_I*Dummy_Choice_Info + mu_rent_NI*Dummy_Choice_No_Info)*(b_ASC_sq + b_mu_natural * Naturalness_3 + b_mu_walking * WalkingDistance_3 + mu_asc_I *  Dummy_Choice_Info +
-                                mu_asc_NI *  Dummy_Choice_No_Info +mu_nat_I * Naturalness_3 * Dummy_Choice_Info + mu_wd_I * WalkingDistance_3 * Dummy_Choice_Info + 
-                                mu_nat_NI * Naturalness_3 * Dummy_Choice_No_Info + mu_wd_NI * WalkingDistance_3 * Dummy_Choice_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)
+  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)*(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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info)*(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 + mu_rent_treated *Dummy_Treated + mu_rent_vol_treated * Dummy_Vol_Treated + mu_rent_no_info * Dummy_no_info)*(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
     
-    ### 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_rentINT = apollo_estimate(apollo_beta, apollo_fixed,
-                        apollo_probabilities, apollo_inputs, 
-                        estimate_settings=list(maxIterations=400,
-                                               estimationRoutine="bfgs",
-                                               hessianRoutine="analytic"))
+  ### 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)
   
-  # ################################################################# #
-  #### MODEL OUTPUTS                                               ##
-  # ################################################################# #
-  apollo_saveOutput(mxl_wtp_case_d_rentINT)
+  ### 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)