mnl_PMI.R

### 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_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_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_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 = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)

# ################################################################# #
#### MODEL OUTPUTS                                               ####
# ################################################################# #

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN)                               ----
# ----------------------------------------------------------------- #

apollo_modelOutput(model_PMI)

# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name)               ----
# ----------------------------------------------------------------- #
summary(model_PMI)
max_PMI=max(database$PMI)
min_PMI=min(database$PMI)

x_list=seq(min_PMI,max_PMI,0.1)

df <- data.frame(Mean=double(),
                 Index=double(),
                 UB=double(),
                 LB=double())
i=1

deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold/b_profit)*100")))
apollo_deltaMethod(model_PMI, deltaMethod_settings)
deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold_PMI/b_profit)*100")))
apollo_deltaMethod(model_PMI, deltaMethod_settings)

for(j in x_list){
  deltaMethod_settings=list(expression=c(paste0(mean_rate="(b_hold/b_profit+b_hold_PMI/b_profit*",j,")*100")))
  a<-apollo_deltaMethod(model_PMI, deltaMethod_settings)
  df[i,"Mean"]<-a[["Value"]]
  df[i,"Index"]<-j
  df[i,"UB"]<-a[["Value"]]+1.96*a[["Robust s.e."]]
  df[i,"LB"]<-a[["Value"]]-1.96*a[["Robust s.e."]]
  i=i+1
}


fig_PMI<-ggplot(df, aes(x=Index, y=Mean)) + 
  geom_ribbon(aes(ymin=LB, ymax=UB),alpha=.2, fill="#00BFC4") +
  geom_line(color="#00BFC4") +
  labs(x = "PMI", y= "Money discount rate in %")
ggsave("Graphs/PMI.png")
ggsave("Graphs/PMI.tiff")