how to simulate the MDCEV probability?

[NiluCapgi]

How to simulate the MDCEV probability?

[dpalma]

Hi,

The function apollo_mdcev has an option called rawPrediction you can use to simulate choices from an MDCEV model.

The way the forecasting of MDCEV works, is that -for each individual- it simulates data from the model several times (100 times by default), and then averages the forecasting across all those simulations. By setting mdcev_settings$rawPrediction=TRUE, instead of getting the average across all those simulations, you will get each of them.

Below is the example. Notice that you will need Apollo v0.3.6 (the latest version) for it to work.

Best wishes,
David

Code: Select all

# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS                       ####
# ################################################################# #

### Initialise
rm(list = ls())
library(apollo)
apollo_setWorkDir()
apollo_initialise()

### Set core controls
apollo_control = list(
  modelName  = "mdcevSimulation",  
  modelDescr = "MDCEV on time use, alpha-gamma profile, with outside good",
  indivID    = "indivID"
)

# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
# ################################################################# #

### Loading data from within package
database = apollo_timeUseData
### for data dictionary, use ?apollo_timeUseData

### Create consumption variables for combined activities
# outside good: time spent at home and travelling
database$t_outside = rowSums(database[,c("t_a01", "t_a06", "t_a10", 
                                         "t_a11", "t_a12")]) 
database$t_leisure = rowSums(database[,c("t_a07", "t_a08", "t_a09")])

# ################################################################# #
#### DEFINE MODEL PARAMETERS                                     ####
# ################################################################# #

### Vector of parameters, including any that are kept fixed in estimation
apollo_beta = c(alpha_base         = -19,
                gamma_work         = 1,
                gamma_school       = 1,
                gamma_shopping     = 1,
                gamma_private      = 1,
                gamma_leisure      = 1,
                delta_work         = 0,
                delta_school       = 0,
                delta_shopping     = 0,
                delta_private      = 0,
                delta_leisure      = 0,
                delta_work_FT      = 0,
                delta_work_wknd    = 0,
                delta_school_young = 0,
                delta_leisure_wknd = 0,
                sig                = 1)

### Names of parameters to keep their values fixed throughout estimation
apollo_fixed = c("alpha_base", "sig")

# ################################################################# #
#### 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()
  
  ### Define individual alternatives
  alternatives  = c("outside",
                    "work",
                    "school",
                    "shopping",
                    "private", 
                    "leisure")
  
  ### Define availabilities
  avail = list(outside  = 1, 
               work     = 1, 
               school   = 1, 
               shopping = 1, 
               private  = 1,
               leisure  = 1)
  
  ### Define continuous consumption for individual alternatives
  continuousChoice = list(outside  = t_outside/60,
                          work     = t_a02/60,
                          school   = t_a03/60,
                          shopping = t_a04/60,
                          private  = t_a05/60,
                          leisure  = t_leisure/60)
  
  ### Define utilities for individual alternatives
  V = list()
  V[["outside"]]  = 0
  V[["work"]]     = delta_work     + delta_work_FT * occ_full_time + delta_work_wknd * weekend
  V[["school"]]   = delta_school   + delta_school_young * (age<=30)
  V[["shopping"]] = delta_shopping
  V[["private"]]  = delta_private
  V[["leisure"]]  = delta_leisure  + delta_leisure_wknd*weekend
  
  ### Define alpha parameters
  alpha = list(outside  = 1 /(1 + exp(-alpha_base)), 
               work     = 1 /(1 + exp(-alpha_base)), 
               school   = 1 /(1 + exp(-alpha_base)), 
               shopping = 1 /(1 + exp(-alpha_base)), 
               private  = 1 /(1 + exp(-alpha_base)),
               leisure  = 1 /(1 + exp(-alpha_base)))
  
  ### Define gamma parameters
  gamma = list(work     = gamma_work,    
               school   = gamma_school,
               shopping = gamma_shopping,
               private  = gamma_private,
               leisure  = gamma_leisure)
  
  ### Define costs for individual alternatives
  cost = list(outside  = 1, 
              work     = 1, 
              school   = 1, 
              shopping = 1, 
              private  = 1,
              leisure  = 1)
  
  ### Define settings for MDCEV model
  mdcev_settings = list(alternatives      = alternatives,
                        avail             = avail,
                        continuousChoice  = continuousChoice,
                        utilities         = V,
                        alpha             = alpha,
                        gamma             = gamma, 
                        sigma             = sig, 
                        cost              = cost,
                        budget            = 24, 
                        rawPrediction     = TRUE) #### THIS LINE IS NEW
  
  ### Compute probabilities using MDCEV model
  P[["model"]] = apollo_mdcev(mdcev_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 AND OUTPUT                                 ####
# ################################################################# #

model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, 
                        apollo_inputs, list(writeIter=FALSE))

### Output to screen
apollo_modelOutput(model)


# ################################################################# #
#### SIMULATION.                                                 ####
# ################################################################# #

pred <- apollo_prediction(model, apollo_probabilities, apollo_inputs)
# First simulation, one row per observation
simulation1 <- pred[,,1]