diff --git a/R_Scripts/mnl_PMI_CMI.R b/R_Scripts/mnl_PMI_CMI.R
new file mode 100644
index 0000000000000000000000000000000000000000..d8b2c6b0bb9b051209ec3fd0ff364c9a85a5fc03
--- /dev/null
+++ b/R_Scripts/mnl_PMI_CMI.R
@@ -0,0 +1,178 @@
+
+### Initialise code
+apollo_initialise()
+
+### Set core controls
+apollo_control = list(
+  modelName       = "MNL_SP",
+  modelDescr      = "Simple MNL model on mode choice SP data",
+  indivID         = "id",
+  outputDirectory = "Output",
+  panelData       = TRUE
+)
+
+
+
+# ################################################################# #
+#### DEFINE MODEL PARAMETERS                                     ####
+# ################################################################# #
+
+### Vector of parameters, including any that are kept fixed in estimation
+apollo_beta=c(asc_a      = 0,
+              asc_b      = 0,
+              b_profit    =1,
+              b_hold      = 0,
+              b_loc       = 0,
+              b_admin1     = 0,
+              b_admin2     = 0,
+              b_hold_CMI = 0,
+              b_hold_PMI = 0
+              
+              
+              
+              
+)
+
+### Vector with names (in quotes) of parameters to be kept fixed at their starting value in apollo_beta, use apollo_beta_fixed = c() if none
+apollo_fixed = c()
+
+# ################################################################# #
+#### GROUP AND VALIDATE INPUTS                                   ####
+# ################################################################# #
+
+apollo_inputs = apollo_validateInputs()
+
+# ################################################################# #
+#### DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
+# ################################################################# #
+
+apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
+  
+  ### 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: these must use the same names as in mnl_settings, order is irrelevant
+  V = list()
+  V[["alt1"]]  = asc_a  +
+    b_profit  * profit1 - 
+    b_hold * hold1 - 
+    b_hold_PMI*hold1*PMI-
+    b_hold_CMI*hold1*CMI+
+    b_loc * loc1 + 
+    b_admin1 *admin_11+
+    b_admin2 *admin_21
+  
+  
+  
+  V[["alt2"]]  =asc_b  +
+    b_profit  * profit2 - 
+    b_hold * hold2 - 
+    b_hold_PMI*hold2*PMI -
+    b_hold_CMI*hold2*CMI+
+    b_loc * loc2 + 
+    b_admin1 *admin_12+
+    b_admin2 *admin_22 
+  
+  
+  
+  V[["alt3"]]  =       b_profit  * profit3
+  
+  ### Define settings for MNL model component
+  mnl_settings = list(
+    alternatives  = c(alt1=1, alt2=2, alt3=3), 
+    avail         = 1, 
+    choiceVar     = choosen,
+    utilities     = V
+  )
+  
+  ### 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)
+  
+  ### Prepare and return outputs of function
+  P = apollo_prepareProb(P, apollo_inputs, functionality)
+  return(P)
+}
+
+# ################################################################# #
+#### MODEL ESTIMATION                                            ####
+# ################################################################# #
+
+model_PMI_CMI = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
+
+# ################################################################# #
+#### MODEL OUTPUTS                                               ####
+# ################################################################# #
+
+# ----------------------------------------------------------------- #
+#---- FORMATTED OUTPUT (TO SCREEN)                               ----
+# ----------------------------------------------------------------- #
+
+apollo_modelOutput(model_PMI_CMI)
+
+# ----------------------------------------------------------------- #
+#---- FORMATTED OUTPUT (TO FILE, using model name)               ----
+# ----------------------------------------------------------------- #
+summary(model_PMI_CMI)
+deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold/b_profit)*100")))
+apollo_deltaMethod(model_PMI_CMI, deltaMethod_settings)
+deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold_PMI/b_profit)*100")))
+apollo_deltaMethod(model_PMI_CMI, deltaMethod_settings)
+deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold_CMI/b_profit)*100")))
+apollo_deltaMethod(model_PMI_CMI, deltaMethod_settings)
+
+
+
+max_PMI=max(database$PMI)
+min_PMI=min(database$PMI)
+
+x_list_PMI=seq(min_PMI,max_PMI,0.1)
+
+max_CMI=max(database$CMI)
+min_CMI=min(database$CMI)
+
+x_list_CMI=seq(min_CMI,max_CMI,0.1)
+
+df <- data.frame(Mean=double(),
+                 Index1=double(),
+                 UB=double(),
+                 LB=double(),
+                 Index=character())
+i=1
+for(j in x_list_CMI){
+  deltaMethod_settings=list(expression=c(paste0(mean_rate="100*(b_hold/b_profit+b_hold_CMI/b_profit*",j,")")))
+  a<-apollo_deltaMethod(model_PMI_CMI, deltaMethod_settings)
+  df[i,"Mean"]<-a[["Value"]]
+  df[i,"Index1"]<-j
+  df[i,"UB"]<-a[["Value"]]+1.96*a[["Robust s.e."]]
+  df[i,"LB"]<-a[["Value"]]-1.96*a[["Robust s.e."]]
+  df[i,"Index"]<-"CMI"
+  i=i+1
+}
+
+for(j in x_list_PMI){
+  deltaMethod_settings=list(expression=c(paste0(mean_rate="100*(b_hold/b_profit+b_hold_PMI/b_profit*",j,")")))
+  a<-apollo_deltaMethod(model_PMI_CMI, deltaMethod_settings)
+  df[i,"Mean"]<-a[["Value"]]
+  df[i,"Index1"]<-j
+  df[i,"UB"]<-a[["Value"]]+1.96*a[["Robust s.e."]]
+  df[i,"LB"]<-a[["Value"]]-1.96*a[["Robust s.e."]]
+  df[i,"Index"]<-"PMI"
+  i=i+1
+}
+
+fig_PMI_CMI<-ggplot(df, aes(x=Index1, y=Mean, color=Index, fill=Index)) + 
+  geom_ribbon(aes(ymin=LB, ymax=UB ),alpha=.2 ) +
+  geom_line() +
+  labs(x = "Motivation index", y= "Money discount rate in %")
+ggsave("Graphs/PMI_CMI.png")
+
+ggsave("Graphs/PMI_CMI.tiff")
+
+