Nothing
rmdcev vs Apollo: Parameter Recovery and Timing Benchmark
In rmdcev: Kuhn-Tucker and Multiple Discrete-Continuous Extreme Value Models
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval = FALSE # heavy computation — run manually
)
This vignette compares rmdcev and apollo on two dimensions:
- Parameter recovery — do both packages converge to the same estimates on synthetic data?
- Timing — how does wall-clock time scale with sample size ($n$) and number of alternatives ($J$)?
Apollo must be installed separately (install.packages("apollo")). All chunks use
eval = FALSE; run the script interactively to reproduce the results.
Setup
library(rmdcev)
library(apollo)
library(dplyr)
library(tidyr)
library(bench) # for timing; install.packages("bench") if needed
Helper functions
Wide-format conversion for Apollo
Apollo requires one row per individual (wide format). This function converts the long-format
output of GenerateMDCEVData() to the wide format Apollo expects, and computes the numeraire
quantity $x_0 = \text{income} - \sum_j p_j x_j$.
to_wide <- function(data_long, nalts) {
data_long <- as.data.frame(data_long)
# Individual-level covariates (b-cols that don't vary within id)
# b1..b_npsi_j are alt-specific; remaining b-cols are individual-specific
all_b <- grep("^b[0-9]", names(data_long), value = TRUE)
# detect individual-specific cols: same value for all rows of an id
ind_b <- all_b[sapply(all_b, function(v) {
all(tapply(data_long[[v]], data_long$id, function(x) length(unique(x)) == 1))
})]
ind_vars <- data_long %>% distinct(id, income, across(all_of(ind_b)))
alt_vars <- data_long %>%
select(id, alt, quant, price, all_of(setdiff(all_b, ind_b))) %>%
pivot_wider(names_from = alt, values_from = c(quant, price, setdiff(all_b, ind_b)))
db <- left_join(ind_vars, alt_vars, by = "id")
# Numeraire quantity
pj_xj <- Reduce("+", lapply(1:nalts, function(j) {
db[[paste0("price_", j)]] * db[[paste0("quant_", j)]]
}))
db$quant_num <- db$income - pj_xj
db
}
Apollo probability function factory
Returns a apollo_probabilities closure for the hybrid0 or gamma model with nalts
non-numeraire alternatives, npsi_j alt-varying psi covariates and one individual-level
psi covariate (b_ind).
make_apollo_prob <- function(nalts, model = c("hybrid0", "gamma"),
npsi_j = 1, ind_col = "b2") {
model <- match.arg(model)
alt_names <- paste0("alt", 1:nalts)
all_names <- c("outside", alt_names)
function(apollo_beta, apollo_inputs, functionality = "estimate") {
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P <- list()
# Utility for each non-numeraire alternative
V <- setNames(vector("list", nalts + 1), all_names)
V[["outside"]] <- 0
for (j in seq_len(nalts)) {
vj <- get(ind_col) * psi_ind # individual-specific component
for (k in seq_len(npsi_j)) {
vj <- vj + get(paste0("b", k, "_", j)) * get(paste0("psi_j", k))
}
V[[paste0("alt", j)]] <- vj
}
# Alpha: hybrid0 = 0 for all; gamma = free for outside, 0 for others
alpha_out <- if (model == "gamma") exp(ln_alpha_out) / (1 + exp(ln_alpha_out)) else 0
alpha <- setNames(
c(list(alpha_out), lapply(seq_len(nalts), function(j) 0)),
all_names
)
gamma <- setNames(
c(list(1), lapply(seq_len(nalts), function(j) exp(get(paste0("ln_g", j))))),
all_names
)
sigma <- exp(ln_sigma)
cChoice <- setNames(
c(list(quant_num), lapply(seq_len(nalts), function(j) get(paste0("quant_", j)))),
all_names
)
cost <- setNames(
c(list(1), lapply(seq_len(nalts), function(j) get(paste0("price_", j)))),
all_names
)
avail <- setNames(
c(list(1), lapply(seq_len(nalts), function(j) 1)),
all_names
)
mdcev_settings <- list(
alternatives = all_names,
budget = income,
V = V,
alpha = alpha,
gamma = gamma,
sigma = sigma,
continuousChoice = cChoice,
cost = cost,
avail = avail
)
P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
P <- apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
}
Apollo starting-values vector
make_apollo_beta <- function(nalts, model = c("hybrid0", "gamma"), npsi_j = 1) {
model <- match.arg(model)
betas <- c(psi_ind = 0)
for (k in seq_len(npsi_j)) betas[paste0("psi_j", k)] <- 0
for (j in seq_len(nalts)) betas[paste0("ln_g", j)] <- 0
if (model == "gamma") betas["ln_alpha_out"] <- 0
betas["ln_sigma"] <- 0
betas
}
Fit both packages, return comparison table
fit_both <- function(sim, model = "hybrid0", npsi_j = 1, ind_col = "b2") {
nalts <- sim$data %>% distinct(alt) %>% nrow()
# ---- rmdcev ----
rmdcev_formula <- as.formula(paste0(
"~ ", paste(paste0("b", 1:(npsi_j + 1)), collapse = " + ")
))
t_rmdcev <- system.time(
fit_r <- mdcev(
formula = rmdcev_formula,
data = sim$data,
psi_ascs = 0,
model = model,
algorithm = "MLE",
print_iterations = FALSE,
backend = "rstan"
)
)
# ---- Apollo ----
database <- to_wide(sim$data, nalts)
db_col <- paste0("b", npsi_j + 1) # individual-specific column
apollo_control <- list(
modelName = paste0("apollo_", model),
modelDescr = paste0("Apollo ", model),
indivID = "id",
outputDirectory = tempdir(),
nCores = 1
)
apollo_beta <- make_apollo_beta(nalts, model, npsi_j)
apollo_fixed <- character(0)
apollo_inputs <- apollo_validateInputs(
apollo_beta = apollo_beta,
apollo_fixed = apollo_fixed,
database = database,
apollo_control = apollo_control
)
apollo_prob <- make_apollo_prob(nalts, model, npsi_j, db_col)
t_apollo <- system.time(
fit_a <- apollo_estimate(
apollo_beta, apollo_fixed, apollo_prob, apollo_inputs,
estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
printLevel = 0L)
)
)
# ---- Extract estimates ----
r_psi <- as.numeric(fit_r$stan_fit$par$psi)
r_gamma <- as.numeric(fit_r$stan_fit$par$gamma)
r_scale <- as.numeric(fit_r$stan_fit$par$scale)
a_est <- fit_a$estimate
a_gamma <- exp(a_est[grep("^ln_g", names(a_est))])
a_scale <- exp(a_est["ln_sigma"])
# Map column names: psi_j1 = alt-varying, psi_ind = individual-varying
a_psi_j <- a_est[grep("^psi_j", names(a_est))]
a_psi_ind <- a_est["psi_ind"]
a_psi <- c(a_psi_j, a_psi_ind)
list(
ll = c(rmdcev = fit_r$log.likelihood, apollo = fit_a$LLout),
psi = rbind(rmdcev = r_psi, apollo = a_psi),
gamma = rbind(rmdcev = r_gamma[1:min(5, nalts)],
apollo = a_gamma[1:min(5, nalts)]),
scale = c(rmdcev = r_scale, apollo = as.numeric(a_scale)),
time_s = c(rmdcev = t_rmdcev["elapsed"], apollo = t_apollo["elapsed"])
)
}
1 · Parameter recovery
We generate data from a known hybrid0 DGP with $n = 500$, $J = 5$, verify that both
packages recover the true parameters, and that their log-likelihoods agree.
set.seed(42)
true_psi_j <- 0.5
true_psi_ind <- -1.5
true_gamma <- rep(5, 5)
true_scale <- 1.0
sim_h0 <- GenerateMDCEVData(
model = "hybrid0",
nobs = 500,
nalts = 5,
psi_j_parms = true_psi_j,
psi_i_parms = true_psi_ind,
gamma_parms = true_gamma,
scale_parms = true_scale
)
cmp_h0 <- fit_both(sim_h0, model = "hybrid0", npsi_j = 1)
cat("=== hybrid0 parameter comparison ===\n")
cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n",
cmp_h0$ll["rmdcev"], cmp_h0$ll["apollo"]))
cat(sprintf(" LL difference = %.6f\n",
abs(cmp_h0$ll["rmdcev"] - cmp_h0$ll["apollo"])))
cat("\nPsi (true: psi_j=0.5, psi_ind=-1.5)\n")
print(round(cmp_h0$psi, 4))
cat("\nGamma (true: all 5)\n")
print(round(cmp_h0$gamma, 4))
cat(sprintf("\nScale (true: 1.0): rmdcev=%.4f apollo=%.4f\n",
cmp_h0$scale["rmdcev"], cmp_h0$scale["apollo"]))
Repeat for the gamma model (alpha for numeraire is free):
set.seed(42)
true_alpha <- 0.5
sim_g <- GenerateMDCEVData(
model = "gamma",
nobs = 500,
nalts = 5,
psi_j_parms = true_psi_j,
psi_i_parms = true_psi_ind,
gamma_parms = true_gamma,
alpha_parms = true_alpha,
scale_parms = true_scale
)
cmp_g <- fit_both(sim_g, model = "gamma", npsi_j = 1)
cat("=== gamma parameter comparison ===\n")
cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n",
cmp_g$ll["rmdcev"], cmp_g$ll["apollo"]))
cat(sprintf(" LL difference = %.6f\n",
abs(cmp_g$ll["rmdcev"] - cmp_g$ll["apollo"])))
cat("\nPsi:\n"); print(round(cmp_g$psi, 4))
cat("\nGamma (first 5):\n"); print(round(cmp_g$gamma, 4))
cat(sprintf("\nScale: rmdcev=%.4f apollo=%.4f\n",
cmp_g$scale["rmdcev"], cmp_g$scale["apollo"]))
2 · Timing benchmark
We time MLE estimation across five (n, J) combinations with 3 replicates each.
Demand simulation timing uses the same grid.
2a · MLE estimation timing
grid <- expand.grid(
n = c(500, 1000, 5000),
J = c(50, 100),
rep = 1:3,
stringsAsFactors = FALSE
) %>%
filter(!(n == 1000 & J == 100)) # keep n=1000 for J=50 only
timing_mle <- lapply(seq_len(nrow(grid)), function(i) {
n <- grid$n[i]; J <- grid$J[i]; r <- grid$rep[i]
set.seed(r * 100 + i)
sim <- GenerateMDCEVData(
model = "hybrid0",
nobs = n,
nalts = J,
psi_j_parms = 0.5,
psi_i_parms = -1.5,
gamma_parms = rep(5, J),
scale_parms = 1.0
)
res <- fit_both(sim, model = "hybrid0", npsi_j = 1)
data.frame(n=n, J=J, rep=r,
rmdcev_s = res$time_s["rmdcev"],
apollo_s = res$time_s["apollo"],
ll_diff = abs(res$ll["rmdcev"] - res$ll["apollo"]))
})
timing_mle <- do.call(rbind, timing_mle)
summary_mle <- timing_mle %>%
group_by(n, J) %>%
summarise(
rmdcev_mean_s = round(mean(rmdcev_s), 1),
apollo_mean_s = round(mean(apollo_s), 1),
speedup = round(mean(apollo_s) / mean(rmdcev_s), 1),
ll_diff_mean = signif(mean(ll_diff), 3),
.groups = "drop"
)
cat("=== MLE estimation timing (seconds, 3 replicates) ===\n")
print(summary_mle, n = Inf)
3 · Demand simulation: 1% price increase
We compare predicted demand changes from a 1% cost increase using rmdcev's built-in
simulation pipeline and Apollo's apollo_prediction() function.
Policy specification
The price_p slot in rmdcev takes additive price changes. For a 1% increase we
add $0.01 \times \bar{p}_j$ (1% of the cross-individual mean price for each alternative).
To ensure identical policies, Apollo's database applies the same constant additive change
per alternative ($+0.01 \bar{p}_j$) rather than a proportional scale. Apollo's prediction
is then run twice — baseline and policy — and the difference gives the demand change.
3a · Apollo apollo_prediction() approach (convergence check, n=500, J=5)
apollo_prediction() does source-code inspection and cannot accept factory-generated
closures, so we write a concrete function for the J=5 convergence check.
# Concrete apollo_probabilities for J=5 (used only in this section)
apollo_probabilities_j5 <- function(apollo_beta, apollo_inputs, functionality = "estimate") {
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P <- list()
V <- list(
outside = 0,
alt1 = psi_ind * b2 + psi_j1 * b1_1,
alt2 = psi_ind * b2 + psi_j1 * b1_2,
alt3 = psi_ind * b2 + psi_j1 * b1_3,
alt4 = psi_ind * b2 + psi_j1 * b1_4,
alt5 = psi_ind * b2 + psi_j1 * b1_5
)
mdcev_settings <- list(
alternatives = c("outside","alt1","alt2","alt3","alt4","alt5"),
budget = income, V = V,
alpha = list(outside=0, alt1=0, alt2=0, alt3=0, alt4=0, alt5=0),
gamma = list(outside=1, alt1=exp(ln_g1), alt2=exp(ln_g2), alt3=exp(ln_g3),
alt4=exp(ln_g4), alt5=exp(ln_g5)),
sigma = exp(ln_sigma),
continuousChoice = list(outside=quant_num, alt1=quant_1, alt2=quant_2,
alt3=quant_3, alt4=quant_4, alt5=quant_5),
cost = list(outside=1, alt1=price_1, alt2=price_2,
alt3=price_3, alt4=price_4, alt5=price_5),
avail = list(outside=1, alt1=1, alt2=1, alt3=1, alt4=1, alt5=1)
)
P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
P <- apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}
set.seed(42)
J5 <- 5
n500 <- 500
sim5 <- GenerateMDCEVData(
model = "hybrid0",
nobs = n500,
nalts = J5,
psi_j_parms = 0.5,
psi_i_parms = -1.5,
gamma_parms = rep(5, J5),
scale_parms = 1.0
)
# ── Fit both packages ──────────────────────────────────────────────────────
fit_r5 <- mdcev(
formula = ~ b1 + b2, data = sim5$data, psi_ascs = 0,
model = "hybrid0", algorithm = "MLE",
print_iterations = FALSE, backend = "rstan"
)
database5 <- to_wide(sim5$data, J5)
ac5 <- list(modelName = "apl_sim5", modelDescr = "sim", indivID = "id",
outputDirectory = tempdir(), nCores = 1)
ab5 <- make_apollo_beta(J5, "hybrid0", npsi_j = 1)
ai5 <- suppressMessages(
apollo_validateInputs(apollo_beta = ab5, apollo_fixed = character(0),
database = database5, apollo_control = ac5))
fit_a5 <- suppressMessages(suppressWarnings(
apollo_estimate(ab5, character(0), apollo_probabilities_j5, ai5,
estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
printLevel = 0L))
))
# ── 1% price-change policy ─────────────────────────────────────────────────
baseline_prices <- matrix(as.numeric(fit_r5$stan_data$price_j), nrow = n500)
price_change <- 0.01 * colMeans(baseline_prices) # 1% of mean price per alt
# rmdcev: additive policy vector (0 for numeraire, then per-alt changes)
policies5 <- CreateBlankPolicies(npols = 2, fit_r5, price_change_only = TRUE)
policies5$price_p[[2]] <- c(0, price_change)
df_sim5 <- PrepareSimulationData(fit_r5, policies5, nsims = 1)
t_rmdcev5 <- system.time(
demand_r5 <- mdcev.sim(
df_sim5$df_indiv, df_common = df_sim5$df_common,
sim_options = df_sim5$sim_options,
cond_err = 0, nerrs = 100, sim_type = "demand"
)
)
# Extract mean demand per alt: demand[[indiv]][[sim]][policy, alt_index]
mean_demand <- function(demand_list, policy_idx, alt_idx) {
mean(sapply(demand_list, function(ind)
mean(sapply(ind, function(s) s[policy_idx, alt_idx]))))
}
d_base_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 1, a))
d_policy_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 2, a))
delta_r <- d_policy_r - d_base_r
# Apollo: predict at baseline and policy databases
# Policy database: same additive change as rmdcev (+ 0.01 * mean_price_j per alt)
database5_pol <- database5
for (j in 1:J5) {
database5_pol[[paste0("price_", j)]] <- database5[[paste0("price_", j)]] + price_change[j]
}
# Recompute quant_num at observed quantities under new prices
pj_xj_pol <- Reduce("+", lapply(1:J5, function(j)
database5_pol[[paste0("price_", j)]] * database5_pol[[paste0("quant_", j)]]))
database5_pol$quant_num <- database5_pol$income - pj_xj_pol
ai5_pol <- suppressMessages(
apollo_validateInputs(apollo_beta = fit_a5$estimate, apollo_fixed = character(0),
database = database5_pol, apollo_control = ac5))
t_apollo5_base <- system.time(
pred_base <- suppressMessages(suppressWarnings(
apollo_prediction(fit_a5, apollo_probabilities_j5, ai5,
prediction_settings = list(nRep = 100, silent = TRUE))
))
)
t_apollo5_pol <- system.time(
pred_pol <- suppressMessages(suppressWarnings(
apollo_prediction(fit_a5, apollo_probabilities_j5, ai5_pol,
prediction_settings = list(nRep = 100, silent = TRUE))
))
)
t_apollo5 <- t_apollo5_base + t_apollo5_pol
delta_a <- sapply(1:J5, function(j) {
mean(pred_pol[[paste0("alt", j, "_cont_mean")]]) -
mean(pred_base[[paste0("alt", j, "_cont_mean")]])
})
cat("=== Demand change: 1% price increase (n=500, J=5) ===\n")
cat(sprintf("%-12s %8s %8s %8s\n", "Alt", "rmdcev", "Apollo", "Diff"))
for (j in 1:J5) {
cat(sprintf("alt%-9d %8.4f %8.4f %8.4f\n", j, delta_r[j], delta_a[j],
abs(delta_r[j] - delta_a[j])))
}
cat(sprintf("\nMean abs difference: %.5f\n", mean(abs(delta_r - delta_a))))
cat(sprintf("rmdcev %% change: %s%%\n",
paste(round(100 * delta_r / d_base_r, 2), collapse = ", ")))
cat(sprintf("\nTiming: rmdcev %.2f s | Apollo (2 predictions) %.2f s\n",
t_rmdcev5["elapsed"], t_apollo5["elapsed"]))
3b · Demand simulation timing (rmdcev vs Apollo, n=500)
Apollo has no batch simulation pipeline for multiple policies. The equivalent workflow
requires calling apollo_prediction() twice per policy (baseline + policy database) with
re-validated inputs, looping over policies in R. Because apollo_prediction() requires a
concrete (non-factory) probability function, we generate one dynamically for each J using
eval(parse()).
Both packages use 100 draws (nerrs = 100 for rmdcev; nRep = 100 for Apollo). rmdcev
runs unconditional draws via its compiled C++ engine; Apollo runs its interpreted R loop.
Timings are the mean of 3 replicates.
# Dynamically generate a concrete apollo_probabilities function for any J.
# apollo_prediction() does source-code inspection and cannot use factory closures,
# so we build and eval() a function body with the right number of alternatives.
make_apollo_prob_concrete <- function(nalts, npsi_j = 1, ind_col = "b2") {
all_str <- paste0('c("outside", ', paste0('"alt', seq_len(nalts), '"', collapse = ", "), ")")
V_items <- paste(paste0("alt", seq_len(nalts), " = psi_ind * ", ind_col,
" + psi_j1 * b1_", seq_len(nalts)), collapse = ",\n ")
alpha_items <- paste(c("outside = 0", paste0("alt", seq_len(nalts), " = 0")), collapse = ", ")
gamma_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = exp(ln_g", seq_len(nalts), ")")), collapse = ", ")
choice_items <- paste(c("outside = quant_num", paste0("alt", seq_len(nalts), " = quant_", seq_len(nalts))), collapse = ", ")
cost_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = price_", seq_len(nalts))), collapse = ", ")
avail_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = 1")), collapse = ", ")
eval(parse(text = sprintf(
'function(apollo_beta, apollo_inputs, functionality = "estimate") {
apollo_attach(apollo_beta, apollo_inputs)
on.exit(apollo_detach(apollo_beta, apollo_inputs))
P <- list()
V <- list(outside = 0, %s)
mdcev_settings <- list(
alternatives = %s, budget = income, V = V,
alpha = list(%s), gamma = list(%s), sigma = exp(ln_sigma),
continuousChoice = list(%s), cost = list(%s), avail = list(%s)
)
P[["model"]] <- apollo_mdcev(mdcev_settings, functionality)
P <- apollo_prepareProb(P, apollo_inputs, functionality)
return(P)
}', V_items, all_str, alpha_items, gamma_items, choice_items, cost_items, avail_items)))
}
timing_sim <- lapply(c(50, 100), function(J) {
n <- 500
set.seed(J)
sim <- GenerateMDCEVData(
model = "hybrid0", nobs = n, nalts = J,
psi_j_parms = 0.5, psi_i_parms = -1.5,
gamma_parms = rep(5, J), scale_parms = 1.0
)
# ── rmdcev ──────────────────────────────────────────────────────────────
fit_r <- mdcev(~ b1 + b2, data = sim$data, psi_ascs = 0, model = "hybrid0",
algorithm = "MLE", print_iterations = FALSE, backend = "rstan")
bp <- matrix(as.numeric(fit_r$stan_data$price_j), nrow = n)
pc <- 0.01 * colMeans(bp)
pols <- CreateBlankPolicies(npols = 2, fit_r, price_change_only = TRUE)
pols$price_p[[2]] <- c(0, pc)
df_s <- PrepareSimulationData(fit_r, pols, nsims = 1)
rmdcev_times <- sapply(1:3, function(r)
system.time(mdcev.sim(df_s$df_indiv, df_common = df_s$df_common,
sim_options = df_s$sim_options,
cond_err = 0, nerrs = 100, sim_type = "demand"))["elapsed"])
# ── Apollo ──────────────────────────────────────────────────────────────
database <- to_wide(sim$data, J)
ac <- list(modelName = paste0("apl_", J), modelDescr = "", indivID = "id",
outputDirectory = tempdir(), nCores = 1)
ab <- make_apollo_beta(J, npsi_j = 1)
pfn <- make_apollo_prob_concrete(J)
ai <- apollo_validateInputs(apollo_beta = ab, apollo_fixed = character(0),
database = database, apollo_control = ac)
fit_a <- apollo_estimate(ab, character(0), pfn, ai,
estimate_settings = list(writeIter = FALSE, hessianRoutine = "none",
printLevel = 0L))
# Policy database: 1% additive price increase
db_pol <- database
for (j in seq_len(J)) db_pol[[paste0("price_", j)]] <- database[[paste0("price_", j)]] + pc[j]
pj_xj_pol <- Reduce("+", lapply(seq_len(J), function(j)
db_pol[[paste0("price_", j)]] * db_pol[[paste0("quant_", j)]]))
db_pol$quant_num <- db_pol$income - pj_xj_pol
ai_pol <- apollo_validateInputs(apollo_beta = fit_a$estimate, apollo_fixed = character(0),
database = db_pol, apollo_control = ac)
apollo_times <- sapply(1:3, function(r) {
t1 <- system.time(apollo_prediction(fit_a, pfn, ai,
prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
t2 <- system.time(apollo_prediction(fit_a, pfn, ai_pol,
prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"]
t1 + t2
})
data.frame(n = n, J = J,
rmdcev_mean = round(mean(rmdcev_times), 2),
rmdcev_sd = round(sd(rmdcev_times), 2),
apollo_mean = round(mean(apollo_times), 2),
apollo_sd = round(sd(apollo_times), 2),
speedup = round(mean(apollo_times) / mean(rmdcev_times), 1))
})
timing_sim <- do.call(rbind, timing_sim)
print(timing_sim, row.names = FALSE)
4 · Results summary
MLE estimation
Results from a Windows laptop (Intel i7, 16 GB RAM), hybrid0 model, 3 replicates per cell.
| n | J | rmdcev (s) | Apollo (s) | Speedup | LL diff |
|--:|--:|-----------:|-----------:|--------:|--------:|
| 500 | 50 | 1.3 | 7.2 | 5.7× | 3.5e-04 |
| 1000 | 50 | 2.6 | 9.5 | 3.6× | 2.4e-04 |
| 500 | 100 | 5.1 | 34.4 | 6.8× | 3.6e-04 |
| 5000 | 50 | 24.3 | 40.2 | 1.7× | 2.3e-03 |
| 5000 | 100 | 53.1 | 170.0 | 3.2× | 2.5e-02 |
Demand simulation (100 draws per individual, 3 replicates, n=500)
Fit once per cell; rmdcev timed 3× with 100 unconditional draws via mdcev.sim(). Apollo
timed 3× as two apollo_prediction() calls (baseline + policy), the minimum required per
policy scenario. Apollo has no batch pipeline — each policy requires a separate call with
re-validated inputs. Results from rmdcev 1.3.2 (Windows laptop, Intel i7, 16 GB RAM).
| n | J | rmdcev mean (s) | rmdcev sd (s) | Apollo mean (s) | Apollo sd (s) | Speedup |
|--:|--:|----------------:|--------------:|----------------:|--------------:|--------:|
| 500 | 50 | 1.44 | 0.06 | 15.49 | 2.36 | 10.7× |
| 500 | 100 | 3.18 | 0.15 | 30.24 | 5.55 | 9.5× |
Simulation time scales approximately linearly with $n \times J$ (rmdcev's compiled C++
loop runs at ~4–5 million individual-alternative-draw operations per second on this hardware).
Demand change convergence (n=500, J=5, 1% price increase, same additive policy)
Both packages apply an additive constant of $0.01 \times \bar{p}_j$ per alternative
and use 100 MC draws (nerrs = 100 for rmdcev; nRep = 100 for Apollo).
rmdcev completed in 0.50 s; Apollo (two apollo_prediction() calls) took 3.60 s (~7×).
| Alt | rmdcev Δ demand | Apollo Δ demand | \|Diff\| |
|----:|----------------:|----------------:|---------:|
| 1 | −1.4993 | −1.4948 | 0.0044 |
| 2 | −1.4325 | −1.4409 | 0.0085 |
| 3 | −1.1040 | −1.0960 | 0.0080 |
| 4 | −1.3146 | −1.3202 | 0.0056 |
| 5 | −1.2062 | −1.2006 | 0.0057 |
Mean absolute difference: 0.006 — MC noise from independent RNG streams, not a model discrepancy.
All alternatives show ~1.5% demand reduction from the 1% cost increase, consistent with
near-unit elasticity around the estimated satiation parameters.
Key findings:
- Demand convergence: when policies are specified identically and both packages use
100 draws, rmdcev and Apollo agree to within ~0.006 units across all alternatives.
The residual difference reflects MC noise from independent RNG streams, not a model
discrepancy.
- Simulation speed: rmdcev's compiled C++ engine is 10.7× faster at n=500, J=50
and 9.5× faster at n=500, J=100. Apollo's
apollo_prediction() requires two separate
calls per policy (baseline + policy database with re-validated inputs) and runs the
Pinjari-Bhat demand loop in interpreted R. rmdcev processes all policies in a single
compiled C++ pass, making Apollo impractical for large n × J grids or multi-policy
scenarios.
- Log-likelihood: rmdcev and Apollo agree to <0.05 LL units across all MLE cells.
- MLE speed: rmdcev is 1.7–6.8× faster than Apollo, driven by Stan analytical
gradients vs Apollo's finite-difference Jacobian. The gap grows with J because each
γ parameter requires an extra likelihood evaluation in Apollo.
Session info
sessionInfo()
Try the rmdcev package in your browser
Any scripts or data that you put into this service are public.
rmdcev documentation built on Sept. 15, 2026, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE # heavy computation — run manually )
This vignette compares rmdcev and apollo on two dimensions:
- Parameter recovery — do both packages converge to the same estimates on synthetic data?
- Timing — how does wall-clock time scale with sample size ($n$) and number of alternatives ($J$)?
Apollo must be installed separately (install.packages("apollo")). All chunks use
eval = FALSE; run the script interactively to reproduce the results.
Setup
library(rmdcev) library(apollo) library(dplyr) library(tidyr) library(bench) # for timing; install.packages("bench") if needed
Helper functions
Wide-format conversion for Apollo
Apollo requires one row per individual (wide format). This function converts the long-format
output of GenerateMDCEVData() to the wide format Apollo expects, and computes the numeraire
quantity $x_0 = \text{income} - \sum_j p_j x_j$.
to_wide <- function(data_long, nalts) { data_long <- as.data.frame(data_long) # Individual-level covariates (b-cols that don't vary within id) # b1..b_npsi_j are alt-specific; remaining b-cols are individual-specific all_b <- grep("^b[0-9]", names(data_long), value = TRUE) # detect individual-specific cols: same value for all rows of an id ind_b <- all_b[sapply(all_b, function(v) { all(tapply(data_long[[v]], data_long$id, function(x) length(unique(x)) == 1)) })] ind_vars <- data_long %>% distinct(id, income, across(all_of(ind_b))) alt_vars <- data_long %>% select(id, alt, quant, price, all_of(setdiff(all_b, ind_b))) %>% pivot_wider(names_from = alt, values_from = c(quant, price, setdiff(all_b, ind_b))) db <- left_join(ind_vars, alt_vars, by = "id") # Numeraire quantity pj_xj <- Reduce("+", lapply(1:nalts, function(j) { db[[paste0("price_", j)]] * db[[paste0("quant_", j)]] })) db$quant_num <- db$income - pj_xj db }
Apollo probability function factory
Returns a apollo_probabilities closure for the hybrid0 or gamma model with nalts
non-numeraire alternatives, npsi_j alt-varying psi covariates and one individual-level
psi covariate (b_ind).
make_apollo_prob <- function(nalts, model = c("hybrid0", "gamma"), npsi_j = 1, ind_col = "b2") { model <- match.arg(model) alt_names <- paste0("alt", 1:nalts) all_names <- c("outside", alt_names) function(apollo_beta, apollo_inputs, functionality = "estimate") { apollo_attach(apollo_beta, apollo_inputs) on.exit(apollo_detach(apollo_beta, apollo_inputs)) P <- list() # Utility for each non-numeraire alternative V <- setNames(vector("list", nalts + 1), all_names) V[["outside"]] <- 0 for (j in seq_len(nalts)) { vj <- get(ind_col) * psi_ind # individual-specific component for (k in seq_len(npsi_j)) { vj <- vj + get(paste0("b", k, "_", j)) * get(paste0("psi_j", k)) } V[[paste0("alt", j)]] <- vj } # Alpha: hybrid0 = 0 for all; gamma = free for outside, 0 for others alpha_out <- if (model == "gamma") exp(ln_alpha_out) / (1 + exp(ln_alpha_out)) else 0 alpha <- setNames( c(list(alpha_out), lapply(seq_len(nalts), function(j) 0)), all_names ) gamma <- setNames( c(list(1), lapply(seq_len(nalts), function(j) exp(get(paste0("ln_g", j))))), all_names ) sigma <- exp(ln_sigma) cChoice <- setNames( c(list(quant_num), lapply(seq_len(nalts), function(j) get(paste0("quant_", j)))), all_names ) cost <- setNames( c(list(1), lapply(seq_len(nalts), function(j) get(paste0("price_", j)))), all_names ) avail <- setNames( c(list(1), lapply(seq_len(nalts), function(j) 1)), all_names ) mdcev_settings <- list( alternatives = all_names, budget = income, V = V, alpha = alpha, gamma = gamma, sigma = sigma, continuousChoice = cChoice, cost = cost, avail = avail ) P[["model"]] <- apollo_mdcev(mdcev_settings, functionality) P <- apollo_prepareProb(P, apollo_inputs, functionality) return(P) } }
Apollo starting-values vector
make_apollo_beta <- function(nalts, model = c("hybrid0", "gamma"), npsi_j = 1) { model <- match.arg(model) betas <- c(psi_ind = 0) for (k in seq_len(npsi_j)) betas[paste0("psi_j", k)] <- 0 for (j in seq_len(nalts)) betas[paste0("ln_g", j)] <- 0 if (model == "gamma") betas["ln_alpha_out"] <- 0 betas["ln_sigma"] <- 0 betas }
Fit both packages, return comparison table
fit_both <- function(sim, model = "hybrid0", npsi_j = 1, ind_col = "b2") { nalts <- sim$data %>% distinct(alt) %>% nrow() # ---- rmdcev ---- rmdcev_formula <- as.formula(paste0( "~ ", paste(paste0("b", 1:(npsi_j + 1)), collapse = " + ") )) t_rmdcev <- system.time( fit_r <- mdcev( formula = rmdcev_formula, data = sim$data, psi_ascs = 0, model = model, algorithm = "MLE", print_iterations = FALSE, backend = "rstan" ) ) # ---- Apollo ---- database <- to_wide(sim$data, nalts) db_col <- paste0("b", npsi_j + 1) # individual-specific column apollo_control <- list( modelName = paste0("apollo_", model), modelDescr = paste0("Apollo ", model), indivID = "id", outputDirectory = tempdir(), nCores = 1 ) apollo_beta <- make_apollo_beta(nalts, model, npsi_j) apollo_fixed <- character(0) apollo_inputs <- apollo_validateInputs( apollo_beta = apollo_beta, apollo_fixed = apollo_fixed, database = database, apollo_control = apollo_control ) apollo_prob <- make_apollo_prob(nalts, model, npsi_j, db_col) t_apollo <- system.time( fit_a <- apollo_estimate( apollo_beta, apollo_fixed, apollo_prob, apollo_inputs, estimate_settings = list(writeIter = FALSE, hessianRoutine = "none", printLevel = 0L) ) ) # ---- Extract estimates ---- r_psi <- as.numeric(fit_r$stan_fit$par$psi) r_gamma <- as.numeric(fit_r$stan_fit$par$gamma) r_scale <- as.numeric(fit_r$stan_fit$par$scale) a_est <- fit_a$estimate a_gamma <- exp(a_est[grep("^ln_g", names(a_est))]) a_scale <- exp(a_est["ln_sigma"]) # Map column names: psi_j1 = alt-varying, psi_ind = individual-varying a_psi_j <- a_est[grep("^psi_j", names(a_est))] a_psi_ind <- a_est["psi_ind"] a_psi <- c(a_psi_j, a_psi_ind) list( ll = c(rmdcev = fit_r$log.likelihood, apollo = fit_a$LLout), psi = rbind(rmdcev = r_psi, apollo = a_psi), gamma = rbind(rmdcev = r_gamma[1:min(5, nalts)], apollo = a_gamma[1:min(5, nalts)]), scale = c(rmdcev = r_scale, apollo = as.numeric(a_scale)), time_s = c(rmdcev = t_rmdcev["elapsed"], apollo = t_apollo["elapsed"]) ) }
1 · Parameter recovery
We generate data from a known hybrid0 DGP with $n = 500$, $J = 5$, verify that both packages recover the true parameters, and that their log-likelihoods agree.
set.seed(42) true_psi_j <- 0.5 true_psi_ind <- -1.5 true_gamma <- rep(5, 5) true_scale <- 1.0 sim_h0 <- GenerateMDCEVData( model = "hybrid0", nobs = 500, nalts = 5, psi_j_parms = true_psi_j, psi_i_parms = true_psi_ind, gamma_parms = true_gamma, scale_parms = true_scale ) cmp_h0 <- fit_both(sim_h0, model = "hybrid0", npsi_j = 1) cat("=== hybrid0 parameter comparison ===\n") cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n", cmp_h0$ll["rmdcev"], cmp_h0$ll["apollo"])) cat(sprintf(" LL difference = %.6f\n", abs(cmp_h0$ll["rmdcev"] - cmp_h0$ll["apollo"]))) cat("\nPsi (true: psi_j=0.5, psi_ind=-1.5)\n") print(round(cmp_h0$psi, 4)) cat("\nGamma (true: all 5)\n") print(round(cmp_h0$gamma, 4)) cat(sprintf("\nScale (true: 1.0): rmdcev=%.4f apollo=%.4f\n", cmp_h0$scale["rmdcev"], cmp_h0$scale["apollo"]))
Repeat for the gamma model (alpha for numeraire is free):
set.seed(42) true_alpha <- 0.5 sim_g <- GenerateMDCEVData( model = "gamma", nobs = 500, nalts = 5, psi_j_parms = true_psi_j, psi_i_parms = true_psi_ind, gamma_parms = true_gamma, alpha_parms = true_alpha, scale_parms = true_scale ) cmp_g <- fit_both(sim_g, model = "gamma", npsi_j = 1) cat("=== gamma parameter comparison ===\n") cat(sprintf("LL: rmdcev = %.2f apollo = %.2f\n", cmp_g$ll["rmdcev"], cmp_g$ll["apollo"])) cat(sprintf(" LL difference = %.6f\n", abs(cmp_g$ll["rmdcev"] - cmp_g$ll["apollo"]))) cat("\nPsi:\n"); print(round(cmp_g$psi, 4)) cat("\nGamma (first 5):\n"); print(round(cmp_g$gamma, 4)) cat(sprintf("\nScale: rmdcev=%.4f apollo=%.4f\n", cmp_g$scale["rmdcev"], cmp_g$scale["apollo"]))
2 · Timing benchmark
We time MLE estimation across five (n, J) combinations with 3 replicates each. Demand simulation timing uses the same grid.
2a · MLE estimation timing
grid <- expand.grid( n = c(500, 1000, 5000), J = c(50, 100), rep = 1:3, stringsAsFactors = FALSE ) %>% filter(!(n == 1000 & J == 100)) # keep n=1000 for J=50 only
timing_mle <- lapply(seq_len(nrow(grid)), function(i) { n <- grid$n[i]; J <- grid$J[i]; r <- grid$rep[i] set.seed(r * 100 + i) sim <- GenerateMDCEVData( model = "hybrid0", nobs = n, nalts = J, psi_j_parms = 0.5, psi_i_parms = -1.5, gamma_parms = rep(5, J), scale_parms = 1.0 ) res <- fit_both(sim, model = "hybrid0", npsi_j = 1) data.frame(n=n, J=J, rep=r, rmdcev_s = res$time_s["rmdcev"], apollo_s = res$time_s["apollo"], ll_diff = abs(res$ll["rmdcev"] - res$ll["apollo"])) }) timing_mle <- do.call(rbind, timing_mle)
summary_mle <- timing_mle %>% group_by(n, J) %>% summarise( rmdcev_mean_s = round(mean(rmdcev_s), 1), apollo_mean_s = round(mean(apollo_s), 1), speedup = round(mean(apollo_s) / mean(rmdcev_s), 1), ll_diff_mean = signif(mean(ll_diff), 3), .groups = "drop" ) cat("=== MLE estimation timing (seconds, 3 replicates) ===\n") print(summary_mle, n = Inf)
3 · Demand simulation: 1% price increase
We compare predicted demand changes from a 1% cost increase using rmdcev's built-in
simulation pipeline and Apollo's apollo_prediction() function.
Policy specification
The price_p slot in rmdcev takes additive price changes. For a 1% increase we
add $0.01 \times \bar{p}_j$ (1% of the cross-individual mean price for each alternative).
To ensure identical policies, Apollo's database applies the same constant additive change
per alternative ($+0.01 \bar{p}_j$) rather than a proportional scale. Apollo's prediction
is then run twice — baseline and policy — and the difference gives the demand change.
3a · Apollo apollo_prediction() approach (convergence check, n=500, J=5)
apollo_prediction() does source-code inspection and cannot accept factory-generated
closures, so we write a concrete function for the J=5 convergence check.
# Concrete apollo_probabilities for J=5 (used only in this section) apollo_probabilities_j5 <- function(apollo_beta, apollo_inputs, functionality = "estimate") { apollo_attach(apollo_beta, apollo_inputs) on.exit(apollo_detach(apollo_beta, apollo_inputs)) P <- list() V <- list( outside = 0, alt1 = psi_ind * b2 + psi_j1 * b1_1, alt2 = psi_ind * b2 + psi_j1 * b1_2, alt3 = psi_ind * b2 + psi_j1 * b1_3, alt4 = psi_ind * b2 + psi_j1 * b1_4, alt5 = psi_ind * b2 + psi_j1 * b1_5 ) mdcev_settings <- list( alternatives = c("outside","alt1","alt2","alt3","alt4","alt5"), budget = income, V = V, alpha = list(outside=0, alt1=0, alt2=0, alt3=0, alt4=0, alt5=0), gamma = list(outside=1, alt1=exp(ln_g1), alt2=exp(ln_g2), alt3=exp(ln_g3), alt4=exp(ln_g4), alt5=exp(ln_g5)), sigma = exp(ln_sigma), continuousChoice = list(outside=quant_num, alt1=quant_1, alt2=quant_2, alt3=quant_3, alt4=quant_4, alt5=quant_5), cost = list(outside=1, alt1=price_1, alt2=price_2, alt3=price_3, alt4=price_4, alt5=price_5), avail = list(outside=1, alt1=1, alt2=1, alt3=1, alt4=1, alt5=1) ) P[["model"]] <- apollo_mdcev(mdcev_settings, functionality) P <- apollo_prepareProb(P, apollo_inputs, functionality) return(P) }
set.seed(42) J5 <- 5 n500 <- 500 sim5 <- GenerateMDCEVData( model = "hybrid0", nobs = n500, nalts = J5, psi_j_parms = 0.5, psi_i_parms = -1.5, gamma_parms = rep(5, J5), scale_parms = 1.0 ) # ── Fit both packages ────────────────────────────────────────────────────── fit_r5 <- mdcev( formula = ~ b1 + b2, data = sim5$data, psi_ascs = 0, model = "hybrid0", algorithm = "MLE", print_iterations = FALSE, backend = "rstan" ) database5 <- to_wide(sim5$data, J5) ac5 <- list(modelName = "apl_sim5", modelDescr = "sim", indivID = "id", outputDirectory = tempdir(), nCores = 1) ab5 <- make_apollo_beta(J5, "hybrid0", npsi_j = 1) ai5 <- suppressMessages( apollo_validateInputs(apollo_beta = ab5, apollo_fixed = character(0), database = database5, apollo_control = ac5)) fit_a5 <- suppressMessages(suppressWarnings( apollo_estimate(ab5, character(0), apollo_probabilities_j5, ai5, estimate_settings = list(writeIter = FALSE, hessianRoutine = "none", printLevel = 0L)) )) # ── 1% price-change policy ───────────────────────────────────────────────── baseline_prices <- matrix(as.numeric(fit_r5$stan_data$price_j), nrow = n500) price_change <- 0.01 * colMeans(baseline_prices) # 1% of mean price per alt # rmdcev: additive policy vector (0 for numeraire, then per-alt changes) policies5 <- CreateBlankPolicies(npols = 2, fit_r5, price_change_only = TRUE) policies5$price_p[[2]] <- c(0, price_change) df_sim5 <- PrepareSimulationData(fit_r5, policies5, nsims = 1) t_rmdcev5 <- system.time( demand_r5 <- mdcev.sim( df_sim5$df_indiv, df_common = df_sim5$df_common, sim_options = df_sim5$sim_options, cond_err = 0, nerrs = 100, sim_type = "demand" ) ) # Extract mean demand per alt: demand[[indiv]][[sim]][policy, alt_index] mean_demand <- function(demand_list, policy_idx, alt_idx) { mean(sapply(demand_list, function(ind) mean(sapply(ind, function(s) s[policy_idx, alt_idx])))) } d_base_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 1, a)) d_policy_r <- sapply(2:(J5 + 1), function(a) mean_demand(demand_r5, 2, a)) delta_r <- d_policy_r - d_base_r # Apollo: predict at baseline and policy databases # Policy database: same additive change as rmdcev (+ 0.01 * mean_price_j per alt) database5_pol <- database5 for (j in 1:J5) { database5_pol[[paste0("price_", j)]] <- database5[[paste0("price_", j)]] + price_change[j] } # Recompute quant_num at observed quantities under new prices pj_xj_pol <- Reduce("+", lapply(1:J5, function(j) database5_pol[[paste0("price_", j)]] * database5_pol[[paste0("quant_", j)]])) database5_pol$quant_num <- database5_pol$income - pj_xj_pol ai5_pol <- suppressMessages( apollo_validateInputs(apollo_beta = fit_a5$estimate, apollo_fixed = character(0), database = database5_pol, apollo_control = ac5)) t_apollo5_base <- system.time( pred_base <- suppressMessages(suppressWarnings( apollo_prediction(fit_a5, apollo_probabilities_j5, ai5, prediction_settings = list(nRep = 100, silent = TRUE)) )) ) t_apollo5_pol <- system.time( pred_pol <- suppressMessages(suppressWarnings( apollo_prediction(fit_a5, apollo_probabilities_j5, ai5_pol, prediction_settings = list(nRep = 100, silent = TRUE)) )) ) t_apollo5 <- t_apollo5_base + t_apollo5_pol delta_a <- sapply(1:J5, function(j) { mean(pred_pol[[paste0("alt", j, "_cont_mean")]]) - mean(pred_base[[paste0("alt", j, "_cont_mean")]]) }) cat("=== Demand change: 1% price increase (n=500, J=5) ===\n") cat(sprintf("%-12s %8s %8s %8s\n", "Alt", "rmdcev", "Apollo", "Diff")) for (j in 1:J5) { cat(sprintf("alt%-9d %8.4f %8.4f %8.4f\n", j, delta_r[j], delta_a[j], abs(delta_r[j] - delta_a[j]))) } cat(sprintf("\nMean abs difference: %.5f\n", mean(abs(delta_r - delta_a)))) cat(sprintf("rmdcev %% change: %s%%\n", paste(round(100 * delta_r / d_base_r, 2), collapse = ", "))) cat(sprintf("\nTiming: rmdcev %.2f s | Apollo (2 predictions) %.2f s\n", t_rmdcev5["elapsed"], t_apollo5["elapsed"]))
3b · Demand simulation timing (rmdcev vs Apollo, n=500)
Apollo has no batch simulation pipeline for multiple policies. The equivalent workflow
requires calling apollo_prediction() twice per policy (baseline + policy database) with
re-validated inputs, looping over policies in R. Because apollo_prediction() requires a
concrete (non-factory) probability function, we generate one dynamically for each J using
eval(parse()).
Both packages use 100 draws (nerrs = 100 for rmdcev; nRep = 100 for Apollo). rmdcev
runs unconditional draws via its compiled C++ engine; Apollo runs its interpreted R loop.
Timings are the mean of 3 replicates.
# Dynamically generate a concrete apollo_probabilities function for any J. # apollo_prediction() does source-code inspection and cannot use factory closures, # so we build and eval() a function body with the right number of alternatives. make_apollo_prob_concrete <- function(nalts, npsi_j = 1, ind_col = "b2") { all_str <- paste0('c("outside", ', paste0('"alt', seq_len(nalts), '"', collapse = ", "), ")") V_items <- paste(paste0("alt", seq_len(nalts), " = psi_ind * ", ind_col, " + psi_j1 * b1_", seq_len(nalts)), collapse = ",\n ") alpha_items <- paste(c("outside = 0", paste0("alt", seq_len(nalts), " = 0")), collapse = ", ") gamma_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = exp(ln_g", seq_len(nalts), ")")), collapse = ", ") choice_items <- paste(c("outside = quant_num", paste0("alt", seq_len(nalts), " = quant_", seq_len(nalts))), collapse = ", ") cost_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = price_", seq_len(nalts))), collapse = ", ") avail_items <- paste(c("outside = 1", paste0("alt", seq_len(nalts), " = 1")), collapse = ", ") eval(parse(text = sprintf( 'function(apollo_beta, apollo_inputs, functionality = "estimate") { apollo_attach(apollo_beta, apollo_inputs) on.exit(apollo_detach(apollo_beta, apollo_inputs)) P <- list() V <- list(outside = 0, %s) mdcev_settings <- list( alternatives = %s, budget = income, V = V, alpha = list(%s), gamma = list(%s), sigma = exp(ln_sigma), continuousChoice = list(%s), cost = list(%s), avail = list(%s) ) P[["model"]] <- apollo_mdcev(mdcev_settings, functionality) P <- apollo_prepareProb(P, apollo_inputs, functionality) return(P) }', V_items, all_str, alpha_items, gamma_items, choice_items, cost_items, avail_items))) }
timing_sim <- lapply(c(50, 100), function(J) { n <- 500 set.seed(J) sim <- GenerateMDCEVData( model = "hybrid0", nobs = n, nalts = J, psi_j_parms = 0.5, psi_i_parms = -1.5, gamma_parms = rep(5, J), scale_parms = 1.0 ) # ── rmdcev ────────────────────────────────────────────────────────────── fit_r <- mdcev(~ b1 + b2, data = sim$data, psi_ascs = 0, model = "hybrid0", algorithm = "MLE", print_iterations = FALSE, backend = "rstan") bp <- matrix(as.numeric(fit_r$stan_data$price_j), nrow = n) pc <- 0.01 * colMeans(bp) pols <- CreateBlankPolicies(npols = 2, fit_r, price_change_only = TRUE) pols$price_p[[2]] <- c(0, pc) df_s <- PrepareSimulationData(fit_r, pols, nsims = 1) rmdcev_times <- sapply(1:3, function(r) system.time(mdcev.sim(df_s$df_indiv, df_common = df_s$df_common, sim_options = df_s$sim_options, cond_err = 0, nerrs = 100, sim_type = "demand"))["elapsed"]) # ── Apollo ────────────────────────────────────────────────────────────── database <- to_wide(sim$data, J) ac <- list(modelName = paste0("apl_", J), modelDescr = "", indivID = "id", outputDirectory = tempdir(), nCores = 1) ab <- make_apollo_beta(J, npsi_j = 1) pfn <- make_apollo_prob_concrete(J) ai <- apollo_validateInputs(apollo_beta = ab, apollo_fixed = character(0), database = database, apollo_control = ac) fit_a <- apollo_estimate(ab, character(0), pfn, ai, estimate_settings = list(writeIter = FALSE, hessianRoutine = "none", printLevel = 0L)) # Policy database: 1% additive price increase db_pol <- database for (j in seq_len(J)) db_pol[[paste0("price_", j)]] <- database[[paste0("price_", j)]] + pc[j] pj_xj_pol <- Reduce("+", lapply(seq_len(J), function(j) db_pol[[paste0("price_", j)]] * db_pol[[paste0("quant_", j)]])) db_pol$quant_num <- db_pol$income - pj_xj_pol ai_pol <- apollo_validateInputs(apollo_beta = fit_a$estimate, apollo_fixed = character(0), database = db_pol, apollo_control = ac) apollo_times <- sapply(1:3, function(r) { t1 <- system.time(apollo_prediction(fit_a, pfn, ai, prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"] t2 <- system.time(apollo_prediction(fit_a, pfn, ai_pol, prediction_settings = list(nRep = 100, silent = TRUE)))["elapsed"] t1 + t2 }) data.frame(n = n, J = J, rmdcev_mean = round(mean(rmdcev_times), 2), rmdcev_sd = round(sd(rmdcev_times), 2), apollo_mean = round(mean(apollo_times), 2), apollo_sd = round(sd(apollo_times), 2), speedup = round(mean(apollo_times) / mean(rmdcev_times), 1)) }) timing_sim <- do.call(rbind, timing_sim) print(timing_sim, row.names = FALSE)
4 · Results summary
MLE estimation
Results from a Windows laptop (Intel i7, 16 GB RAM), hybrid0 model, 3 replicates per cell.
| n | J | rmdcev (s) | Apollo (s) | Speedup | LL diff | |--:|--:|-----------:|-----------:|--------:|--------:| | 500 | 50 | 1.3 | 7.2 | 5.7× | 3.5e-04 | | 1000 | 50 | 2.6 | 9.5 | 3.6× | 2.4e-04 | | 500 | 100 | 5.1 | 34.4 | 6.8× | 3.6e-04 | | 5000 | 50 | 24.3 | 40.2 | 1.7× | 2.3e-03 | | 5000 | 100 | 53.1 | 170.0 | 3.2× | 2.5e-02 |
Demand simulation (100 draws per individual, 3 replicates, n=500)
Fit once per cell; rmdcev timed 3× with 100 unconditional draws via mdcev.sim(). Apollo
timed 3× as two apollo_prediction() calls (baseline + policy), the minimum required per
policy scenario. Apollo has no batch pipeline — each policy requires a separate call with
re-validated inputs. Results from rmdcev 1.3.2 (Windows laptop, Intel i7, 16 GB RAM).
| n | J | rmdcev mean (s) | rmdcev sd (s) | Apollo mean (s) | Apollo sd (s) | Speedup | |--:|--:|----------------:|--------------:|----------------:|--------------:|--------:| | 500 | 50 | 1.44 | 0.06 | 15.49 | 2.36 | 10.7× | | 500 | 100 | 3.18 | 0.15 | 30.24 | 5.55 | 9.5× |
Simulation time scales approximately linearly with $n \times J$ (rmdcev's compiled C++ loop runs at ~4–5 million individual-alternative-draw operations per second on this hardware).
Demand change convergence (n=500, J=5, 1% price increase, same additive policy)
Both packages apply an additive constant of $0.01 \times \bar{p}_j$ per alternative
and use 100 MC draws (nerrs = 100 for rmdcev; nRep = 100 for Apollo).
rmdcev completed in 0.50 s; Apollo (two apollo_prediction() calls) took 3.60 s (~7×).
| Alt | rmdcev Δ demand | Apollo Δ demand | \|Diff\| | |----:|----------------:|----------------:|---------:| | 1 | −1.4993 | −1.4948 | 0.0044 | | 2 | −1.4325 | −1.4409 | 0.0085 | | 3 | −1.1040 | −1.0960 | 0.0080 | | 4 | −1.3146 | −1.3202 | 0.0056 | | 5 | −1.2062 | −1.2006 | 0.0057 |
Mean absolute difference: 0.006 — MC noise from independent RNG streams, not a model discrepancy. All alternatives show ~1.5% demand reduction from the 1% cost increase, consistent with near-unit elasticity around the estimated satiation parameters.
Key findings:
- Demand convergence: when policies are specified identically and both packages use 100 draws, rmdcev and Apollo agree to within ~0.006 units across all alternatives. The residual difference reflects MC noise from independent RNG streams, not a model discrepancy.
- Simulation speed: rmdcev's compiled C++ engine is 10.7× faster at n=500, J=50
and 9.5× faster at n=500, J=100. Apollo's
apollo_prediction()requires two separate calls per policy (baseline + policy database with re-validated inputs) and runs the Pinjari-Bhat demand loop in interpreted R. rmdcev processes all policies in a single compiled C++ pass, making Apollo impractical for large n × J grids or multi-policy scenarios. - Log-likelihood: rmdcev and Apollo agree to <0.05 LL units across all MLE cells.
- MLE speed: rmdcev is 1.7–6.8× faster than Apollo, driven by Stan analytical gradients vs Apollo's finite-difference Jacobian. The gap grows with J because each γ parameter requires an extra likelihood evaluation in Apollo.
Session info
sessionInfo()
Try the rmdcev package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.