# Copyright (c) 2026 INNOVATIO SAS
# SPDX-License-Identifier: MIT
# MSC-P-045 (marketing-science-center.com): value a firm's current and future customers.
# Base R only. Same computations, same order and same printed lines as msc-p045-reference.py: the
# customer-based valuation of Gupta, Lehmann and Stuart (2004) checked on their published inputs, then
# applied to a synthetic online retailer declared as such.

LCG_SEED <- 20261006
LCG_MODULUS <- 2147483648
LCG_INCREMENT <- 12345
LCG_HIGH <- 16838   # 1103515245 = 16838 * 65536 + 20077, split to keep products exact
LCG_LOW <- 20077

TRUE_ALPHA <- 2000000
TRUE_BETA <- -4.0
TRUE_GAMMA <- 0.30
QUARTERS <- 16L
EARLY_QUARTERS <- 8L
NOISE_SD <- 0.10
RETENTION <- 0.80
MARGIN <- 12.0
ACQUISITION_COST <- 40.0
DISCOUNT <- 0.10
BOOTSTRAP <- 1000L
REPLICATIONS <- 100L
COVERAGE_BOOTSTRAP <- 200L  # bootstrap draws per replication, to measure the interval's actual coverage
STEP <- 0.05
HORIZON <- 400.0
NM_MAX_ITER <- 20000L
NM_TOL <- 1e-12
ALPHA_BOUND_FACTOR <- 1000.0
PENALTY <- 1e300

GLS_TAX <- 0.38
GLS_DISCOUNT <- 0.12
GLS <- list(
  list(name = "Amazon", n = 33800000, m = 3.87, c = 7.70, r = 0.70, a = 67.045, b = -4.114, g = 0.265, t = 21, pub = c(0.82, 2.45, 0.07, 1.07, 0.46)),
  list(name = "Ameritrade", n = 1877000, m = 50.39, c = 203.44, r = 0.95, a = 2.482, b = -3.345, g = 0.263, t = 19, pub = c(1.62, 6.75, 0.03, 1.03, 1.17)),
  list(name = "Capital One", n = 46600000, m = 13.71, c = 75.49, r = 0.85, a = 171.200, b = -3.052, g = 0.149, t = 22, pub = c(11.00, 5.12, 0.32, 1.32, 1.11)),
  list(name = "Ebay", n = 46100000, m = 4.31, c = 11.26, r = 0.80, a = 81.945, b = -6.009, g = 0.317, t = 22, pub = c(1.89, 3.42, 0.08, 1.08, 0.63)),
  list(name = "E*Trade", n = 4117370, m = 43.02, c = 391.00, r = 0.95, a = 4.719, b = -3.441, g = 0.365, t = 18, pub = c(2.69, 6.67, 0.02, 1.02, 1.14))
)

# First quarter of each firm's data (Table 1) and calendar quarter of its acquisition peak (Table 2).
GLS_START <- list("Amazon" = c(1997, 3), "Ameritrade" = c(1997, 9), "Capital One" = c(1996, 12), "Ebay" = c(1996, 12), "E*Trade" = c(1997, 12))
GLS_PEAK <- list("Amazon" = c(2000, 12), "Ameritrade" = c(2000, 9), "Capital One" = c(2001, 12), "Ebay" = c(2001, 6), "E*Trade" = c(2000, 3))
# Table 5: customer value ($ billion) at discount 8, 12, 16 % (rows) and retention 70, 80, 90 % (columns).
GLS_TABLE5 <- list(
  "Amazon" = rbind(c(0.90, 1.25, 1.97), c(0.82, 1.10, 1.62), c(0.75, 0.98, 1.38)),
  "Ameritrade" = rbind(c(0.71, 0.97, 1.52), c(0.65, 0.85, 1.25), c(0.59, 0.76, 1.06)),
  "Capital One" = rbind(c(7.33, 10.95, 18.94), c(6.21, 8.88, 14.14), c(5.35, 7.39, 11.04)),
  "Ebay" = rbind(c(1.56, 2.18, 3.47), c(1.39, 1.89, 2.82), c(1.29, 1.70, 2.41)),
  "E*Trade" = rbind(c(1.07, 1.52, 2.46), c(0.98, 1.34, 2.03), c(0.90, 1.20, 1.73))
)
# Values at 100 % retention given in the text: Amazon "about $3 billion", Ebay "$5.3 billion", E*Trade "$3.89".
GLS_FULL_RETENTION <- list("Amazon" = 3.0, "Ebay" = 5.3, "E*Trade" = 3.89)

new_stream <- function(seed) {
  env <- new.env()
  env$state <- seed
  env
}
uniform <- function(s) {
  high <- (LCG_HIGH * s$state) %% LCG_MODULUS
  s$state <- (high * 65536 + LCG_LOW * s$state + LCG_INCREMENT) %% LCG_MODULUS
  (s$state + 0.5) / LCG_MODULUS
}
normal <- function(s) {
  u1 <- uniform(s)
  u2 <- uniform(s)
  sqrt(-2.0 * log(u1)) * cos(2.0 * pi * u2)
}

plain_sum <- function(values) {
  total <- 0
  for (v in values) total <- total + v
  total
}

fail <- function(message) {
  cat(paste0("error: ", message, "\n"), file = stderr())
  quit(status = 1)
}

out <- function(...) cat(paste0(..., "\n"), sep = "")

fmt <- function(x, digits) {
  text <- sprintf(paste0("%.", digits, "f"), x)
  if (startsWith(text, "-") && as.numeric(text) == 0) substring(text, 2) else text
}

pct <- function(x, digits = 1) fmt(100.0 * x, digits)

cumulative <- function(par, t) par[1] / (1.0 + exp(-par[2] - par[3] * t))

acquisition_rate <- function(par, t) {
  e <- exp(-par[2] - par[3] * t)
  par[1] * par[3] * e / ((1.0 + e) * (1.0 + e))
}

quarterly <- function(annual_retention, annual_discount) c(annual_retention^0.25, (1.0 + annual_discount)^0.25 - 1.0)

lifetime_value <- function(margin, rq, iq) margin * rq / (1.0 + iq - rq)

future_cohorts <- function(par, value_per_customer, iq, start) {
  steps <- as.integer(round(HORIZON / STEP))
  log_d <- log(1.0 + iq)
  # Vectorised over the nodes: element by element the same operations, in the same order, as the Python loop;
  # only the sum is kept left to right (plain_sum) so that both languages add the terms identically.
  j <- 0:steps
  k <- start + j * STEP
  w <- ifelse(j == 0 | j == steps, 1.0, ifelse(j %% 2 == 1, 4.0, 2.0))
  terms <- w * acquisition_rate(par, k) * exp(-log_d * (k - start))
  value_per_customer * plain_sum(terms) * STEP / 3.0
}

base_value <- function(par, current, start, margin, cost, annual_retention, annual_discount) {
  q <- quarterly(annual_retention, annual_discount)
  lv <- lifetime_value(margin, q[1], q[2])
  c(current * lv, future_cohorts(par, lv - cost, q[2], start))
}

gls_value <- function(row, margin_mult = 1.0, cost_mult = 1.0, retention = NULL, discount = GLS_DISCOUNT, future_mult = 1.0) {
  r <- if (is.null(retention)) row$r else retention
  v <- base_value(c(row$a * 1e6, row$b, row$g), row$n, row$t, row$m * margin_mult, row$c * cost_mult, r, discount)
  (v[1] + future_mult * v[2]) * (1.0 - GLS_TAX) / 1e9
}

quarter_after <- function(start, quarters) {
  months <- start[1] * 12 + (start[2] - 1) + 3 * quarters
  c(months %/% 12, months %% 12 + 1)
}

nelder_mead <- function(f, start, step) {
  k <- length(start)
  pts <- list(start)
  for (j in seq_len(k)) {
    p <- start
    p[j] <- p[j] + step
    pts[[j + 1]] <- p
  }
  vals <- vapply(pts, f, numeric(1))
  converged <- FALSE
  for (iter in seq_len(NM_MAX_ITER)) {
    ord <- order(vals, seq_along(vals))
    pts <- pts[ord]
    vals <- vals[ord]
    size <- 0
    for (i in 2:(k + 1)) for (j in seq_len(k)) size <- max(size, abs(pts[[i]][j] - pts[[1]][j]))
    if (vals[k + 1] - vals[1] < NM_TOL && size < 1e-8) {
      converged <- TRUE
      break
    }
    centroid <- numeric(k)
    for (j in seq_len(k)) {
      cc <- 0.0
      for (i in seq_len(k)) cc <- cc + pts[[i]][j]
      centroid[j] <- cc / k
    }
    refl <- centroid + (centroid - pts[[k + 1]])
    fr <- f(refl)
    if (fr < vals[1]) {
      exp_pt <- centroid + 2.0 * (centroid - pts[[k + 1]])
      fe <- f(exp_pt)
      if (fe < fr) {
        pts[[k + 1]] <- exp_pt
        vals[k + 1] <- fe
      } else {
        pts[[k + 1]] <- refl
        vals[k + 1] <- fr
      }
    } else if (fr < vals[k]) {
      pts[[k + 1]] <- refl
      vals[k + 1] <- fr
    } else {
      con <- if (fr < vals[k + 1]) centroid + 0.5 * (refl - centroid) else centroid + 0.5 * (pts[[k + 1]] - centroid)
      fc <- f(con)
      if (fc < min(fr, vals[k + 1])) {
        pts[[k + 1]] <- con
        vals[k + 1] <- fc
      } else {
        for (i in 2:(k + 1)) {
          pts[[i]] <- pts[[1]] + 0.5 * (pts[[i]] - pts[[1]])
          vals[i] <- f(pts[[i]])
        }
      }
    }
  }
  best <- which.min(vals)
  list(point = pts[[best]], value = vals[best], converged = converged)
}

triers <- function(active, annual_retention) {
  rq <- annual_retention^0.25
  res <- numeric(length(active))
  res[1] <- active[1]
  if (length(active) > 1) for (t in 2:length(active)) res[t] <- res[t - 1] + active[t] - rq * active[t - 1]
  res
}

fit_curve <- function(cum) {
  scale <- 1e6
  ys <- cum / scale
  log_bound <- log(ALPHA_BOUND_FACTOR * ys[length(ys)])
  ts <- seq_along(ys) - 1
  sse <- function(z) {
    if (z[1] > log_bound || z[3] > 3.0) return(PENALTY)
    a <- exp(z[1]); b <- z[2]; g <- exp(z[3])
    plain_sum((ys - a / (1.0 + exp(-b - g * ts)))^2)
  }
  z <- c(log(2.0 * ys[length(ys)]), -3.0, log(0.2))
  value <- sse(z)
  for (restart in 1:50) {
    res <- nelder_mead(sse, z, 0.5)
    if (!res$converged) fail("the curve fit did not converge within NM_MAX_ITER iterations")
    moved <- value - res$value
    z <- res$point
    value <- res$value
    if (moved < 1e-12) {
      return(list(par = c(exp(z[1]) * scale, z[2], exp(z[3])), sse = value, bound = z[1] > log_bound - 1e-3))
    }
  }
  fail("the curve fit kept moving after 50 restarts")
}

simulate <- function(stream, par, quarters, annual_retention, noise) {
  rq <- annual_retention^0.25
  active <- numeric(quarters + 1)
  active[1] <- cumulative(par, 0)
  for (t in 1:quarters) {
    n_t <- cumulative(par, t) - cumulative(par, t - 1)
    shock <- exp(noise * normal(stream) - 0.5 * noise * noise)
    active[t + 1] <- floor(rq * active[t] + n_t * shock + 0.5)
  }
  active[1] <- floor(active[1] + 0.5)
  active
}

lifetime_shortcut <- function(annual_margin, annual_retention, annual_discount) {
  years <- 1.0 / (1.0 - annual_retention)
  annual_margin * (1.0 - (1.0 + annual_discount)^(-years)) / annual_discount
}

value_all <- function(par, current, start, retention = RETENTION, discount = DISCOUNT, margin = MARGIN, cost = ACQUISITION_COST) {
  v <- base_value(par, current, start, margin, cost, retention, discount)
  v[1] + v[2]
}

quantile_pair <- function(xs) {
  s <- sort(xs)
  c(s[round(0.025 * length(s))], s[round(0.975 * length(s))])
}

noise_sd <- function(resid, parameters) {
  m <- plain_sum(resid) / length(resid)
  sqrt(plain_sum((resid - m)^2) / (length(resid) - parameters))
}

log_residuals <- function(cum, par) {
  out <- numeric(length(cum) - 1)
  for (t in 1:(length(cum) - 1)) out[t] <- log((cum[t + 1] - cum[t]) / (cumulative(par, t) - cumulative(par, t - 1)))
  out
}

mean_sd <- function(xs) {
  m <- plain_sum(xs) / length(xs)
  c(m, sqrt(plain_sum((xs - m)^2) / (length(xs) - 1)))
}

check_invariants <- function() {
  problems <- character(0)
  checked <- 0L
  par <- c(TRUE_ALPHA, TRUE_BETA, TRUE_GAMMA)
  remaining <- future_cohorts(par, 1.0, 0.0, QUARTERS)
  if (abs(remaining - (TRUE_ALPHA - cumulative(par, QUARTERS))) > 1e-3) problems <- c(problems, "Simpson integral of n(t)")
  checked <- checked + 1L
  q <- quarterly(RETENTION, DISCOUNT)
  series <- plain_sum(MARGIN * q[1]^(1:3000) / (1.0 + q[2])^(1:3000))
  if (abs(series - lifetime_value(MARGIN, q[1], q[2])) > 1e-9) problems <- c(problems, "lifetime value series")
  checked <- checked + 1L
  if (abs(lifetime_value(100.0, 0.8, 0.12) - 250.0) > 1e-9) problems <- c(problems, "note 3 lifetime value")
  checked <- checked + 1L
  if (abs(lifetime_shortcut(100.0, 0.8, 0.12) - 360.4776) > 1e-4) problems <- c(problems, "note 3 expected-lifetime value")
  checked <- checked + 1L
  if (abs(triers(c(100000.0, 130000.0), 0.8^4)[2] - 150000.0) > 1e-6) problems <- c(problems, "acquisitions rebuilt from active customers")
  checked <- checked + 1L
  peak <- -TRUE_BETA / TRUE_GAMMA
  if (!(acquisition_rate(par, peak) > acquisition_rate(par, peak - 0.01) &&
        acquisition_rate(par, peak) > acquisition_rate(par, peak + 0.01))) problems <- c(problems, "acquisition peak")
  checked <- checked + 1L
  clean <- simulate(new_stream(1), par, QUARTERS, RETENTION, 0.0)
  rebuilt <- triers(clean, RETENTION)
  if (max(abs(rebuilt - cumulative(par, 0:QUARTERS))) > 2.0 * QUARTERS) problems <- c(problems, "triers rebuilt without noise")
  checked <- checked + 1L
  if (length(problems)) fail(paste0("invariant failed: ", paste(problems, collapse = ", ")))
  checked
}

published_checks <- function() {
  out("== Check on the published case: Gupta, Lehmann and Stuart, Tables 1 to 5 ==")
  out("convention: quarterly retention r^(1/4), quarterly discount 1.12^(1/4) - 1, t = 1 for the first quarter of data, tax 38 %")
  out("note 3: lifetime value ", fmt(lifetime_value(100.0, 0.8, 0.12), 2), " (first margin after one period); ",
      "equation (2) summed from t = 0 would give ", fmt(100.0 * 1.12 / (1.12 - 0.8), 2), "; ",
      "expected-lifetime value ", fmt(lifetime_shortcut(100.0, 0.8, 0.12), 2), ", ", pct(lifetime_shortcut(100.0, 0.8, 0.12) / 250.0 - 1.0), " % above 250; ",
      "with the first margin also counted at the end of year 1, ", fmt(100.0 / (1.12 - 0.8), 2), " and ", pct(lifetime_shortcut(100.0, 0.8, 0.12) / (100.0 / (1.12 - 0.8)) - 1.0), " %")
  reproduced <- 0L
  for (row in GLS) {
    v <- gls_value(row)
    e_ret <- gls_value(row, retention = row$r * 1.01) / v - 1.0
    e_acq <- gls_value(row, cost_mult = 0.99) / v - 1.0
    e_mar <- gls_value(row, margin_mult = 1.01) / v - 1.0
    e_dis <- gls_value(row, discount = GLS_DISCOUNT * 0.99) / v - 1.0
    e_dis_118 <- gls_value(row, discount = 0.118) / v - 1.0
    pub <- row$pub
    ok_value <- abs(v - pub[1]) <= 0.0151
    ok_el <- all(abs(100.0 * c(e_ret, e_acq, e_mar) - pub[2:4]) <= 0.0251)
    if (ok_value && ok_el) reproduced <- reproduced + 1L
    out(row$name, ": value ", fmt(v, 3), " (published ", fmt(pub[1], 2), ", ", if (ok_value) "reproduced" else "NOT reproduced", "); ",
        "1 % better retention ", pct(e_ret, 2), " % (published ", fmt(pub[2], 2), "), ",
        "acquisition cost ", pct(e_acq, 2), " % (", fmt(pub[3], 2), "), margin ", pct(e_mar, 2), " % (", fmt(pub[4], 2), "), ",
        "margin minus acquisition ", fmt(100.0 * (e_mar - e_acq), 3), "; ",
        "discount rate x 0.99 ", pct(e_dis, 2), " %, 12 % to 11.80 % ", pct(e_dis_118, 2), " % (", fmt(pub[5], 2), "); ",
        "retention to discount ratio ", fmt(e_ret / e_dis, 1))
  }
  out("firms whose value and three elasticities are reproduced: ", reproduced, " of 5")
  if (reproduced != 4L) fail("the published check no longer reproduces four firms")
  cap <- GLS[[3]]
  v_cap <- gls_value(cap, future_mult = 1.25)
  out("Capital One with future acquisitions x 1.25: value ", fmt(v_cap, 2), ", retention ", pct(gls_value(cap, retention = 0.85 * 1.01, future_mult = 1.25) / v_cap - 1.0, 2), " %, ",
      "acquisition cost ", pct(gls_value(cap, cost_mult = 0.99, future_mult = 1.25) / v_cap - 1.0, 2), " %, margin ", pct(gls_value(cap, margin_mult = 1.01, future_mult = 1.25) / v_cap - 1.0, 2), " %")
  peaks <- character(0)
  for (row in GLS) {
    peak <- -row$b / row$g
    cells <- character(0)
    for (k in 1:2) {
      label <- c("t = 1", "t = 0")[k]
      offset <- c(1, 0)[k]
      ym <- quarter_after(GLS_START[[row$name]], ceiling(peak) - offset)
      cells <- c(cells, paste0(label, " ", sprintf("%02d", as.integer(ym[2])), "/", as.integer(ym[1])))
    }
    ym <- GLS_PEAK[[row$name]]
    peaks <- c(peaks, paste0(row$name, " ", fmt(peak, 2), ": ", paste(cells, collapse = ", "), " (published ", sprintf("%02d", as.integer(ym[2])), "/", as.integer(ym[1]), ")"))
  }
  out("acquisition peak, calendar quarter: ", paste(peaks, collapse = "; "))
  out("Table 5 recomputed with the curve held fixed (published), discount 12 %:")
  lower <- 0L; higher <- 0L; same <- 0L
  discounts <- c(0.08, 0.12, 0.16)
  retentions <- c(0.70, 0.80, 0.90)
  for (row in GLS) {
    base <- row$r
    cells <- character(0)
    for (d_i in 1:3) {
      for (r_i in 1:3) {
        d <- discounts[d_i]
        r <- retentions[r_i]
        v <- gls_value(row, retention = r, discount = d)
        pub <- GLS_TABLE5[[row$name]][d_i, r_i]
        if (d == 0.12) cells <- c(cells, paste0("r ", pct(r, 0), " % ", fmt(v, 2), " (", fmt(pub, 2), ")"))
        if (abs(r - base) < 1e-9) {
          if (abs(v - pub) <= 0.0151) same <- same + 1L
        } else if ((r > base && pub < v) || (r < base && pub > v)) {
          lower <- lower + 1L
        } else {
          higher <- higher + 1L
        }
      }
    }
    full <- GLS_FULL_RETENTION[[row$name]]
    extra <- if (!is.null(full)) paste0("; retention 100 % ", fmt(gls_value(row, retention = 1.0), 2), " (text ", fmt(full, 2), ")") else ""
    out("  ", row$name, ": ", paste(cells, collapse = "; "), extra)
  }
  out("Table 5 cells away from the firm's own retention: ", lower, " of ", lower + higher, " move the way a refitted curve would; ",
      "cells at the firm's own retention within 0.015: ", same, " of ", 45L - lower - higher)
  out()
}

out("MSC-P-045 reference output")
out("invariants checked: ", check_invariants())
out()
published_checks()

par <- c(TRUE_ALPHA, TRUE_BETA, TRUE_GAMMA)
stream <- new_stream(LCG_SEED)
active <- simulate(stream, par, QUARTERS, RETENTION, NOISE_SD)
out("== Synthetic online retailer ==")
out("design: seed ", LCG_SEED, ", acquisition noise sd ", fmt(NOISE_SD, 2), ", Simpson step ", fmt(STEP, 3), " quarter over ", fmt(HORIZON, 0), " quarters, ",
    "bootstrap ", BOOTSTRAP, ", replications ", REPLICATIONS)
out("true curve: alpha ", fmt(TRUE_ALPHA, 0), ", beta ", fmt(TRUE_BETA, 2), ", gamma ", fmt(TRUE_GAMMA, 2), "; acquisition peak at quarter ", fmt(-TRUE_BETA / TRUE_GAMMA, 2))
out("annual retention ", pct(RETENTION, 0), " %, quarterly margin ", fmt(MARGIN, 2), ", acquisition cost ", fmt(ACQUISITION_COST, 2), ", annual discount ", pct(DISCOUNT, 0), " %")
out("active customers by quarter: ", paste(vapply(active, fmt, character(1), digits = 0), collapse = " "))
q <- quarterly(RETENTION, DISCOUNT)
rq <- q[1]; iq <- q[2]
lv <- lifetime_value(MARGIN, rq, iq)
out("quarterly retention ", fmt(rq, 4), ", quarterly discount ", fmt(100.0 * iq, 3), " %, value of one customer ", fmt(lv, 2))
cum <- triers(active, RETENTION)
out("customers who ever bought, rebuilt: ", fmt(cum[length(cum)], 0), " (true ", fmt(cumulative(par, QUARTERS), 0), ")")
fit <- fit_curve(cum)
est <- fit$par
out("fitted curve: alpha ", fmt(est[1], 0), ", beta ", fmt(est[2], 3), ", gamma ", fmt(est[3], 4), "; peak ", fmt(-est[2] / est[3], 2), "; ",
    "SSE ", fmt(fit$sse, 6), " (millions squared); bound reached: ", if (fit$bound) "yes" else "no")
current <- active[length(active)]
ve <- base_value(est, current, QUARTERS, MARGIN, ACQUISITION_COST, RETENTION, DISCOUNT)
vt <- base_value(par, current, QUARTERS, MARGIN, ACQUISITION_COST, RETENTION, DISCOUNT)
cur_e <- ve[1]; fut_e <- ve[2]
ref <- vt[1] + vt[2]
value <- cur_e + fut_e
out("current customers ", fmt(current, 0), ": value ", fmt(cur_e / 1e6, 2), " M")
out("future customers: estimated ", fmt(fut_e / 1e6, 2), " M, true curve ", fmt(vt[2] / 1e6, 2), " M; ",
    "to acquire: estimated ", fmt(est[1] - cumulative(est, QUARTERS), 0), ", true ", fmt(TRUE_ALPHA - cumulative(par, QUARTERS), 0))
out("value of the customer base: ", fmt(value / 1e6, 2), " M (reference with the true curve ", fmt(ref / 1e6, 2), " M, gap ", pct(value / ref - 1.0), " %); ",
    "future customers ", pct(fut_e / value), " % of it")
out()

out("== Shortcuts ==")
out("current customers only: ", fmt(cur_e / 1e6, 2), " M (", pct(cur_e / ref - 1.0), " %)")
short <- lifetime_shortcut(4.0 * MARGIN, RETENTION, DISCOUNT)
cur_s <- current * short
fut_s <- future_cohorts(est, short - ACQUISITION_COST, iq, QUARTERS)
out("expected-lifetime shortcut: ", fmt(1.0 / (1.0 - RETENTION), 1), " years, ", fmt(short, 2), " per customer instead of ", fmt(lv, 2), " ",
    "(", pct(short / lv - 1.0), " %); base ", fmt((cur_s + fut_s) / 1e6, 2), " M (", pct((cur_s + fut_s) / ref - 1.0), " %)")
annual_next <- 4.0 * MARGIN * RETENTION / (1.0 + DISCOUNT - RETENTION)
annual_first <- 4.0 * MARGIN / (1.0 + DISCOUNT - RETENTION)
out("same shortcut against annual values: ", fmt(annual_next, 2), " with the first margin after one year (", pct(short / annual_next - 1.0), " %), ",
    fmt(annual_first, 2), " with the first margin at the end of year 1 like the shortcut (", pct(short / annual_first - 1.0), " %)")
early <- triers(active[1:(EARLY_QUARTERS + 1)], RETENTION)
fe <- fit_curve(early)
est_e <- fe$par
v_early <- value_all(est_e, active[EARLY_QUARTERS + 1], EARLY_QUARTERS)
ref_early <- value_all(par, active[EARLY_QUARTERS + 1], EARLY_QUARTERS)
out("curve fitted on quarters 0 to ", EARLY_QUARTERS, ", before the peak: alpha ", fmt(est_e[1], 0), ", beta ", fmt(est_e[2], 3), ", gamma ", fmt(est_e[3], 4), "; ",
    "bound reached: ", if (fe$bound) "yes" else "no", "; value at quarter ", EARLY_QUARTERS, " ", fmt(v_early / 1e6, 2), " M ",
    "(reference ", fmt(ref_early / 1e6, 2), " M, ", pct(v_early / ref_early - 1.0), " %)")
out()

out("== Elasticities (curve held fixed, as in Table 4) ==")
out("improved retention: ", pct(RETENTION * 1.01), " % instead of ", pct(RETENTION, 0), " %")
out("retention x 1.01: ", pct(value_all(est, current, QUARTERS, retention = RETENTION * 1.01) / value - 1.0, 2), " %")
out("margin x 1.01: ", pct(value_all(est, current, QUARTERS, margin = MARGIN * 1.01) / value - 1.0, 2), " %")
out("acquisition cost x 0.99: ", pct(value_all(est, current, QUARTERS, cost = ACQUISITION_COST * 0.99) / value - 1.0, 2), " %")
out("discount rate x 0.99: ", pct(value_all(est, current, QUARTERS, discount = DISCOUNT * 0.99) / value - 1.0, 2), " %")
e_r <- value_all(est, current, QUARTERS, retention = RETENTION * 1.01) / value - 1.0
e_m <- value_all(est, current, QUARTERS, margin = MARGIN * 1.01) / value - 1.0
e_c <- value_all(est, current, QUARTERS, cost = ACQUISITION_COST * 0.99) / value - 1.0
e_d <- value_all(est, current, QUARTERS, discount = DISCOUNT * 0.99) / value - 1.0
refit_r <- fit_curve(triers(active, RETENTION * 1.01))$par
out("margin minus acquisition cost: ", fmt(100.0 * (e_m - e_c), 3), " (an identity of the model); retention to discount ratio ", fmt(e_r / e_d, 1), "; ",
    "retention x 1.01 with acquisitions rebuilt and the curve refitted: ", pct(value_all(refit_r, current, QUARTERS, retention = RETENTION * 1.01) / value - 1.0, 2), " %")
out()

out("== Retention and discount grid, value in M (curve held fixed / acquisitions rebuilt and curve refitted) ==")
for (d in c(0.08, 0.10, 0.12)) {
  cells <- character(0)
  for (r in c(0.70, 0.80, 0.90)) {
    fixed <- value_all(est, current, QUARTERS, retention = r, discount = d)
    refit <- fit_curve(triers(active, r))$par
    moved <- value_all(refit, current, QUARTERS, retention = r, discount = d)
    cells <- c(cells, paste0("r ", pct(r, 0), " %: ", fmt(fixed / 1e6, 2), " / ", fmt(moved / 1e6, 2)))
  }
  out("discount ", pct(d, 0), " %: ", paste(cells, collapse = "; "))
}
for (r in c(0.70, 0.90)) {
  refit <- fit_curve(triers(active, r))$par
  rebuilt <- triers(active, r)
  out("assumed retention ", pct(r, 0), " %: rebuilt customers who ever bought ", fmt(rebuilt[length(rebuilt)], 0), ", fitted alpha ", fmt(refit[1], 0))
}
out()

out("== Parametric bootstrap, 1,000 retailers drawn from the fitted curve ==")
sd_r <- noise_sd(log_residuals(cum, est), 3)
out("noise of quarterly acquisitions, estimated with three fitted parameters: ", fmt(sd_r, 4), " (true ", fmt(NOISE_SD, 2), ")")
boot_alpha <- numeric(BOOTSTRAP)
boot_value <- numeric(BOOTSTRAP)
for (b in seq_len(BOOTSTRAP)) {
  sim <- simulate(stream, est, QUARTERS, RETENTION, sd_r)
  b_par <- fit_curve(triers(sim, RETENTION))$par
  boot_alpha[b] <- b_par[1]
  boot_value[b] <- value_all(b_par, current, QUARTERS)
}
qa <- quantile_pair(boot_alpha)
out("alpha: 95 % interval ", fmt(qa[1] / 1000.0, 0), " to ", fmt(qa[2] / 1000.0, 0), " thousand")
qv <- quantile_pair(boot_value)
out("value of the customer base (current customers as observed): 95 % interval ", fmt(qv[1] / 1e6, 2), " to ", fmt(qv[2] / 1e6, 2), " M")
out()

out("== 100 independent retailers drawn from the true curve ==")
full_err <- numeric(REPLICATIONS)
early_err <- numeric(REPLICATIONS)
current_only <- numeric(REPLICATIONS)
early_hit <- logical(REPLICATIONS)
bound_hits <- 0L
covered <- 0L
for (k in seq_len(REPLICATIONS)) {
  sim <- simulate(stream, par, QUARTERS, RETENTION, NOISE_SD)
  sim_cum <- triers(sim, RETENTION)
  f_par <- fit_curve(sim_cum)$par
  ref_s <- value_all(par, sim[length(sim)], QUARTERS)
  f_sd <- noise_sd(log_residuals(sim_cum, f_par), 3)
  draws <- numeric(COVERAGE_BOOTSTRAP)
  for (b in seq_len(COVERAGE_BOOTSTRAP)) {
    b_par <- fit_curve(triers(simulate(stream, f_par, QUARTERS, RETENTION, f_sd), RETENTION))$par
    draws[b] <- value_all(b_par, sim[length(sim)], QUARTERS)
  }
  qb <- quantile_pair(draws)
  if (qb[1] <= ref_s && ref_s <= qb[2]) covered <- covered + 1L
  full_err[k] <- value_all(f_par, sim[length(sim)], QUARTERS) / ref_s - 1.0
  cur_only <- base_value(f_par, sim[length(sim)], QUARTERS, MARGIN, ACQUISITION_COST, RETENTION, DISCOUNT)[1]
  current_only[k] <- cur_only / ref_s - 1.0
  e_fit <- fit_curve(triers(sim[1:(EARLY_QUARTERS + 1)], RETENTION))
  if (e_fit$bound) bound_hits <- bound_hits + 1L
  early_hit[k] <- e_fit$bound
  early_err[k] <- value_all(e_fit$par, sim[EARLY_QUARTERS + 1], EARLY_QUARTERS) / value_all(par, sim[EARLY_QUARTERS + 1], EARLY_QUARTERS) - 1.0
}
share <- covered / REPLICATIONS
out("bootstrap interval (nominal 95 %, ", COVERAGE_BOOTSTRAP, " draws each) contains the reference in ", covered, " of ", REPLICATIONS, " retailers ",
    "(simulation error ", pct(sqrt(share * (1.0 - share) / REPLICATIONS)), " %)")
ms <- mean_sd(full_err)
out("full window: mean gap ", pct(ms[1]), " % (sd ", pct(ms[2]), " %, simulation error ", pct(ms[2] / sqrt(REPLICATIONS)), " %), largest ", pct(max(abs(full_err))), " %")
ms <- mean_sd(current_only)
out("current customers only: mean gap ", pct(ms[1]), " % (sd ", pct(ms[2]), " %, simulation error ", pct(ms[2] / sqrt(REPLICATIONS)), " %)")
s <- sort(early_err)
half <- length(s) %/% 2
within <- sum(abs(early_err) <= 0.20)
free <- sort(early_err[!early_hit])
out("curve fitted before the peak: median gap ", pct((s[half] + s[half + 1]) / 2.0), " %; within 20 %: ", within, " of ", REPLICATIONS, "; ",
    "ceiling at its bound (1,000 times the customers counted): ", bound_hits, "; ",
    "other fits from ", pct(free[1]), " % to ", pct(free[length(free)]), " %")