The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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.These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.