This vignette compares rmdcev and apollo on two dimensions:
Apollo must be installed separately
(install.packages("apollo")). All chunks use
eval = FALSE; run the script interactively to reproduce the
results.
library(rmdcev)
library(apollo)
library(dplyr)
library(tidyr)
library(bench) # for timing; install.packages("bench") if neededApollo 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
}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)
}
}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 <- 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"])
)
}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"]))We time MLE estimation across five (n, J) combinations with 3 replicates each. Demand simulation timing uses the same grid.
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 onlytiming_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)We compare predicted demand changes from a 1% cost increase using
rmdcev’s built-in simulation pipeline and Apollo’s
apollo_prediction() function.
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.
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"]))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)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 |
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).
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:
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.