--- title: "rmdcev vs Apollo: Parameter Recovery and Timing Benchmark" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{rmdcev vs Apollo: Parameter Recovery and Timing Benchmark} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE # heavy computation — run manually ) ``` This vignette compares **rmdcev** and **apollo** on two dimensions: 1. **Parameter recovery** — do both packages converge to the same estimates on synthetic data? 2. **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 ```{r packages} 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$. ```{r fn-to-wide} 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`). ```{r fn-apollo-prob} 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 ```{r fn-apollo-beta} 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 ```{r fn-fit-both} 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. ```{r param-recovery} 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): ```{r param-recovery-gamma} 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 ```{r timing-setup} 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 ``` ```{r timing-mle} 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) ``` ```{r timing-mle-summary} 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. ```{r sim-apollo-prob5} # 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) } ``` ```{r sim-convergence} 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. ```{r make-concrete-prob} # 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))) } ``` ```{r timing-sim-grid} 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 ```{r session} sessionInfo() ```