Comparisons between MDCEV model and EMDC model
Posted: 11 Dec 2023, 08:25
Hi authors,
I have been thinking about the MDCEV model and the EMDC model. If I understand correctly, the EMDC model is a special case of the MDCEV model (or should I say it is MDCEV model with complementarity/substitution). That is to say, if I restrict the complementarity/substitution parameters delta's in an EMDC model, the estimation results/prediction results would be the same as the MDCEV model. Thus, I tried to make some changes to the sample code provided on the website and compare. However, I found that they actually generated very different estimates. I provided my code here.
I am wondering:
1. Why might these two models generate different estimates even when the delta's are fixed in the emdc model? I also tried to do a prediction using the same dataset to see whether the predicted outcome is the same as the observation. Seems that these two models give different prediction outcomes too.
2. I did realize that the MDCEV model needs to specify alpha, while the EMDC model does not need alpha information. Why is that? Could this cause the inconsistency?
3. Related to the above point, I realized that in the MDCEV model, if I do not fix alpha parameter, it does not produce s.e.
MDCEV model code:
Here is the EMDC model codeļ¼
Thank you in advance!
Best,
TW
I have been thinking about the MDCEV model and the EMDC model. If I understand correctly, the EMDC model is a special case of the MDCEV model (or should I say it is MDCEV model with complementarity/substitution). That is to say, if I restrict the complementarity/substitution parameters delta's in an EMDC model, the estimation results/prediction results would be the same as the MDCEV model. Thus, I tried to make some changes to the sample code provided on the website and compare. However, I found that they actually generated very different estimates. I provided my code here.
I am wondering:
1. Why might these two models generate different estimates even when the delta's are fixed in the emdc model? I also tried to do a prediction using the same dataset to see whether the predicted outcome is the same as the observation. Seems that these two models give different prediction outcomes too.
2. I did realize that the MDCEV model needs to specify alpha, while the EMDC model does not need alpha information. Why is that? Could this cause the inconsistency?
3. Related to the above point, I realized that in the MDCEV model, if I do not fix alpha parameter, it does not produce s.e.
MDCEV model code:
Code: Select all
#################################MDCEV model with outside good
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory
rm(list = ls())
### Load Apollo library
library(apollo)
### Initialise code
apollo_initialise()
### Set core controls
apollo_control = list(
modelName ="MDCEV_with_outside_good",
modelDescr ="MDCEV model on time use data, alpha-gamma profile with outside good and socio-demographics",
indivID ="indivID",
outputDirectory = "output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Loading data from package
### if data is to be loaded from a file (e.g. called data.csv),
### the code would be: database = read.csv("data.csv",header=TRUE)
database = apollo_timeUseData
### for data dictionary, use ?apollo_timeUseData
### Create consumption variables for combined activities
database$t_outside = rowSums(database[,c("t_a01", "t_a06", "t_a10", "t_a11", "t_a12")]) # outside good: time spent at home and travelling
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 = -50,
gamma_work = 10.425,
gamma_school = 9.955,
gamma_shopping = 2.065,
gamma_private = 4.419,
gamma_leisure = 6.927,
delta_work = -3.372,
delta_school = -5.250,
delta_shopping = -3.663,
delta_private = -3.973,
delta_leisure = -3.3,
delta_work_FT = 0.708,
delta_work_wknd = -1.6,
delta_school_young = 0.948,
delta_leisure_wknd = 0.155,
sig = 1.944)
### 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("sig" , "alpha_base")
# ################################################################# #
#### 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)
### 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 ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs)
# ################################################################# #
#### MODEL OUTPUTS ####
# ################################################################# #
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO SCREEN) ----
# ----------------------------------------------------------------- #
apollo_modelOutput(model)
# ----------------------------------------------------------------- #
#---- FORMATTED OUTPUT (TO FILE, using model name) ----
# ----------------------------------------------------------------- #
apollo_saveOutput(model)
# ################################################################# #
#### PREDICTION ####
# ################################################################# #
#############################Use the model estimates to predict and see whether it is close to what is observed
model <- apollo_loadModel(apollo_control$modelName)
apollo_inputs <- apollo_validateInputs(database=database)
apollo_inputs$apollo_control$nCores <- 4
pred <- apollo_prediction(model, apollo_probabilities, apollo_inputs)
XObs <- cbind(work = database$t_a02/60,
school = database$t_a03/60,
shopping = database$t_a04/60,
private = database$t_a05/60,
leisure = database$t_leisure/60)
XPre <- pred[,4:8]
round(sqrt(colMeans((XObs - XPre)^2)),2) # RMSE per product: 3.61 0.81 1.97 1.85 3.73
round(sqrt(mean((colSums(XObs) - colSums(XPre))^2)),2) # RMSE for sample: 975.87
Code: Select all
# ################################################################# #
#### LOAD LIBRARY AND DEFINE CORE SETTINGS ####
# ################################################################# #
### Clear memory and initialise
rm(list = ls())
library(apollo)
apollo_initialise()
### Set core controls
apollo_control = list(
modelName ="eMDC_with_budget",
modelDescr ="Extended MDC with complementarity and substitution, with observed budget and socio-demographics",
indivID ="indivID",
outputDirectory="output"
)
# ################################################################# #
#### LOAD DATA AND APPLY ANY TRANSFORMATIONS ####
# ################################################################# #
### Load data from within the Apollo package
database = 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")])
# ### Randomly split dataset into estimation (70%) and validation (30%)
# set.seed(1)
# database$validation <- runif(nrow(database))>0.7
# dbVal <- database[ database$validation,] # validation sample
# database <- database[!database$validation,] # estimation sample
# ################################################################# #
#### DEFINE MODEL PARAMETERS ####
# ################################################################# #
### Parameters starting values c(name1=value1, name2=value2, ...)
apollo_beta = c( sigma = 0.71, ###Here I remove aFemale
# Satiation
gWork = 10.425,
gSchool = 9.955,
gShopping= 2.065,
gPrivate = 4.419,
gLeisure = 6.927,
# Base utility
bWork =-3.372, bSchool =-5.250,
bShopping =-3.663, bPrivate =-3.973,
bLeisure =-3.300, bWork_FT = 0.708,
bWork_wknd=-1.600, bSchool_young= 0.948,
bLeisure_wknd= 0.155,
# Compl/subst
# dWorkScho=-0.008, dWorkShop= 0.000,
# dWorkPriv= 0.000, dWorkLeis= 0.000,
# dSchoShop= 0.000, dSchoPriv= 0.000,
# dSchoLeis= 0.000, dShopPriv= 0.010,
# dShopLeis= 0.012, dPrivLeis= 0.012
dWorkScho=-0.000, dWorkShop= 0.000,
dWorkPriv= 0.000, dWorkLeis= 0.000,
dSchoShop= 0.000, dSchoPriv= 0.000,
dSchoLeis= 0.000, dShopPriv= 0.000,
dShopLeis= 0.000, dPrivLeis= 0.000)
### Names of fixed parameters
apollo_fixed = c('dWorkShop', 'dWorkPriv', 'dSchoShop',
'dSchoPriv', 'dSchoLeis', 'dWorkLeis' ,
"dWorkScho" , "dShopPriv" , "dShopLeis" , "dPrivLeis") #############Here I also fix sigma to the true value 0.71
# ################################################################# #
#### GROUP AND VALIDATE INPUTS ####
# ################################################################# #
apollo_inputs = apollo_validateInputs()
# ################################################################# #
#### DEFINE MODEL AND LIKELIHOOD FUNCTION ####
# ################################################################# #
apollo_probabilities=function(apollo_beta, apollo_inputs,
functionality="estimate"){
### Initialise
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P = list()
### Prepare Inputs
alts = c("work", "school", "shopping", "private", "leisure")
nAlt = length(alts)
ones = setNames(as.list(rep(1, nAlt)), alts)
continuousChoice = list(work = t_a02/60,
school = t_a03/60,
shopping = t_a04/60,
private = t_a05/60,
leisure = t_leisure/60)
utilities = list(
work = bWork + bWork_FT*occ_full_time + bWork_wknd*weekend,
school = bSchool + bSchool_young*(age<=30),
shopping = bShopping,
private = bPrivate,
leisure = bLeisure + bLeisure_wknd*weekend
)
gamma = list(work = gWork,
school = gSchool,
shopping = gShopping,
private = gPrivate,
leisure = gLeisure)
delta <- c(0, 0, 0, 0, 0,
dWorkScho, 0, 0, 0, 0,
dWorkShop, dSchoShop, 0, 0, 0,
dWorkPriv, dSchoPriv, dShopPriv, 0, 0,
dWorkLeis, dSchoLeis, dShopLeis, dPrivLeis, 0)
delta <- matrix(delta, nrow=nAlt, ncol=nAlt, byrow=TRUE)
emdc_settings <- list(continuousChoice = continuousChoice,
avail = ones,
utilityOutside = 0,
utilities = utilities,
budget = 24,
sigma = sigma,
gamma = gamma,
delta = delta,
cost = ones)
P[["model"]] = apollo_emdc(emdc_settings, functionality)
### Comment out as necessary
P = apollo_panelProd(P, apollo_inputs, functionality)
P = apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
# ################################################################# #
#### MODEL ESTIMATION & OUTPUT ####
# ################################################################# #
model = apollo_estimate(apollo_beta, apollo_fixed,
apollo_probabilities, apollo_inputs)
apollo_modelOutput(model)
apollo_saveOutput(model)
# ################################################################# #
#### PREDICTION ####
# ################################################################# #
#############################Use the model estimates to predict and see whether it is close to what is observed
model <- apollo_loadModel(apollo_control$modelName)
apollo_inputs <- apollo_validateInputs(database=database)
apollo_inputs$apollo_control$nCores <- 4
pred <- apollo_prediction(model, apollo_probabilities, apollo_inputs)
XObs <- cbind(work = dbVal$t_a02/60,
school = dbVal$t_a03/60,
shopping = dbVal$t_a04/60,
private = dbVal$t_a05/60,
leisure = dbVal$t_leisure/60)
XPre <- pred[,4:8]
round(sqrt(colMeans((XObs - XPre)^2)),2) #
round(sqrt(mean((colSums(XObs) - colSums(XPre))^2)),2) #
Best,
TW