# Copyright (c) 2026 INNOVATIO SAS
# SPDX-License-Identifier: MIT
"""MSC-P-045 (marketing-science-center.com): value a firm's current and future customers.

Standard library only. The customer-based valuation of Gupta, Lehmann and Stuart (2004) is checked first on
the inputs they publish for five firms (their Tables 1 and 2): the program recomputes their customer values of
March 2002 and their Table 4 elasticities, and prints what it reproduces and what it does not. It is then
applied to a synthetic online retailer, declared as such: the cumulative number of customers who ever bought
follows the S-shaped curve alpha / (1 + exp(-beta - gamma t)); quarterly acquisitions are observed with
multiplicative noise; active customers decay at a constant annual retention rate. The analyst sees the active
customers only, rebuilds acquisitions from an assumed retention rate, fits the curve by nonlinear least
squares and values current and future customers. Three shortcuts are written out in full: current customers
only, the expected-lifetime shortcut, and a curve fitted before its inflection point. The R reference prints
exactly the same lines.
"""
import math

LCG_SEED = 20261006
LCG_MULTIPLIER = 1103515245
LCG_INCREMENT = 12345
LCG_MODULUS = 2147483648

# Synthetic online retailer, quarterly data; t = 1 is the first observed quarter, t = 0 the stock before it.
TRUE_ALPHA = 2000000.0
TRUE_BETA = -4.0
TRUE_GAMMA = 0.30
QUARTERS = 16
EARLY_QUARTERS = 8
NOISE_SD = 0.10
RETENTION = 0.80          # annual
MARGIN = 12.0             # per active customer and quarter, euros
ACQUISITION_COST = 40.0   # per acquired customer, euros
DISCOUNT = 0.10           # annual
BOOTSTRAP = 1000
REPLICATIONS = 100
COVERAGE_BOOTSTRAP = 200  # bootstrap draws per replication, to measure the interval's actual coverage
STEP = 0.05               # Simpson step, quarters
HORIZON = 400.0           # quarters integrated after the valuation date
NM_MAX_ITER = 20000
NM_TOL = 1e-12
# The fitted market size alpha is bounded: above 1,000 times the customers already counted the simplex is
# told "much worse". A fit that ends on the bound is reported as such.
ALPHA_BOUND_FACTOR = 1000.0
PENALTY = 1e300

# Gupta, Lehmann and Stuart: Table 1 (customers in March 2002, quarterly margin, acquisition cost, annual
# retention), Table 2 (alpha in millions, beta, gamma), quarters of data with t = 1 for the first quarter,
# Table 3 (customer value, $ billion) and Table 4 (% increase for a 1 % improvement).
GLS_TAX = 0.38
GLS_DISCOUNT = 0.12
GLS = [
    # name,        customers,   margin, cost,   ret,  alpha,   beta,   gamma, T,  value, ret%,  acq%, mar%, disc%
    ("Amazon",      33800000.0,   3.87,   7.70, 0.70,  67.045, -4.114, 0.265, 21, 0.82, 2.45, 0.07, 1.07, 0.46),
    ("Ameritrade",   1877000.0,  50.39, 203.44, 0.95,   2.482, -3.345, 0.263, 19, 1.62, 6.75, 0.03, 1.03, 1.17),
    ("Capital One", 46600000.0,  13.71,  75.49, 0.85, 171.200, -3.052, 0.149, 22, 11.00, 5.12, 0.32, 1.32, 1.11),
    ("Ebay",        46100000.0,   4.31,  11.26, 0.80,  81.945, -6.009, 0.317, 22, 1.89, 3.42, 0.08, 1.08, 0.63),
    ("E*Trade",      4117370.0,  43.02, 391.00, 0.95,   4.719, -3.441, 0.365, 18, 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 = {"Amazon": (1997, 3), "Ameritrade": (1997, 9), "Capital One": (1996, 12), "Ebay": (1996, 12), "E*Trade": (1997, 12)}
GLS_PEAK = {"Amazon": (2000, 12), "Ameritrade": (2000, 9), "Capital One": (2001, 12), "Ebay": (2001, 6), "E*Trade": (2000, 3)}
# Table 5: customer value ($ billion) at discount 8, 12, 16 % (rows) and retention 70, 80, 90 % (columns).
GLS_TABLE5 = {
    "Amazon": ((0.90, 1.25, 1.97), (0.82, 1.10, 1.62), (0.75, 0.98, 1.38)),
    "Ameritrade": ((0.71, 0.97, 1.52), (0.65, 0.85, 1.25), (0.59, 0.76, 1.06)),
    "Capital One": ((7.33, 10.95, 18.94), (6.21, 8.88, 14.14), (5.35, 7.39, 11.04)),
    "Ebay": ((1.56, 2.18, 3.47), (1.39, 1.89, 2.82), (1.29, 1.70, 2.41)),
    "E*Trade": ((1.07, 1.52, 2.46), (0.98, 1.34, 2.03), (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 = {"Amazon": 3.0, "Ebay": 5.3, "E*Trade": 3.89}


class Declared:
    def __init__(self, seed):
        self.state = seed

    def uniform(self):
        self.state = (LCG_MULTIPLIER * self.state + LCG_INCREMENT) % LCG_MODULUS
        return (self.state + 0.5) / LCG_MODULUS

    def normal(self):
        u1 = self.uniform()
        u2 = self.uniform()
        return math.sqrt(-2.0 * math.log(u1)) * math.cos(2.0 * math.pi * u2)


def plain_sum(values):
    """Left-to-right double sum; Python 3.12+ sum() compensates, R's sum() uses long double."""
    total = 0.0
    for v in values:
        total += v
    return total


def fmt(x, digits):
    text = f"{x:.{digits}f}"
    return text[1:] if text.startswith("-") and float(text) == 0.0 else text


def pct(x, digits=1):
    return fmt(100.0 * x, digits)


def cumulative(par, t):
    a, b, g = par
    return a / (1.0 + math.exp(-b - g * t))


def acquisition_rate(par, t):
    a, b, g = par
    e = math.exp(-b - g * t)
    return a * g * e / ((1.0 + e) * (1.0 + e))


def quarterly(annual_retention, annual_discount):
    return annual_retention ** 0.25, (1.0 + annual_discount) ** 0.25 - 1.0


def lifetime_value(margin, rq, iq):
    """Margins of one customer from the next quarter on, kept with probability rq per quarter."""
    return margin * rq / (1.0 + iq - rq)


def future_cohorts(par, value_per_customer, iq, start):
    """Integral over k > start of n(k) * value_per_customer * (1 + iq)^-(k - start), by Simpson's rule."""
    steps = int(round(HORIZON / STEP))
    log_d = math.log(1.0 + iq)
    terms = []
    for j in range(steps + 1):
        k = start + j * STEP
        w = 1.0 if j == 0 or j == steps else (4.0 if j % 2 == 1 else 2.0)
        terms.append(w * acquisition_rate(par, k) * math.exp(-log_d * (k - start)))
    return value_per_customer * plain_sum(terms) * STEP / 3.0


def base_value(par, current, start, margin, cost, annual_retention, annual_discount):
    """(value of current customers, value of future customers), before tax."""
    rq, iq = quarterly(annual_retention, annual_discount)
    lv = lifetime_value(margin, rq, iq)
    return current * lv, future_cohorts(par, lv - cost, iq, start)


def gls_value(row, margin_mult=1.0, cost_mult=1.0, retention=None, discount=GLS_DISCOUNT, future_mult=1.0):
    name, n, m, c, r, a, b, g, t, *_ = row
    r = r if retention is None else retention
    cur, fut = base_value((a * 1e6, b, g), n, t, m * margin_mult, c * cost_mult, r, discount)
    return (cur + future_mult * fut) * (1.0 - GLS_TAX) / 1e9


def quarter_after(start, quarters):
    """Calendar (year, month) of the quarter that lies `quarters` quarters after `start`."""
    months = start[0] * 12 + (start[1] - 1) + 3 * quarters
    return months // 12, months % 12 + 1


# ---- nonlinear least squares ------------------------------------------------------------------------------

def nelder_mead(f, start, step):
    k = len(start)
    pts = [list(start)]
    for j in range(k):
        p = list(start)
        p[j] += step
        pts.append(p)
    vals = [f(p) for p in pts]
    converged = False
    for _ in range(NM_MAX_ITER):
        order = sorted(range(k + 1), key=lambda i: (vals[i], i))
        pts = [pts[i] for i in order]
        vals = [vals[i] for i in order]
        size = max(abs(pts[i][j] - pts[0][j]) for i in range(1, k + 1) for j in range(k))
        if vals[k] - vals[0] < NM_TOL and size < 1e-8:
            converged = True
            break
        centroid = []
        for j in range(k):
            c = 0.0
            for i in range(k):
                c += pts[i][j]
            centroid.append(c / k)
        refl = [centroid[j] + (centroid[j] - pts[k][j]) for j in range(k)]
        fr = f(refl)
        if fr < vals[0]:
            exp_pt = [centroid[j] + 2.0 * (centroid[j] - pts[k][j]) for j in range(k)]
            fe = f(exp_pt)
            if fe < fr:
                pts[k], vals[k] = exp_pt, fe
            else:
                pts[k], vals[k] = refl, fr
        elif fr < vals[k - 1]:
            pts[k], vals[k] = refl, fr
        else:
            if fr < vals[k]:
                con = [centroid[j] + 0.5 * (refl[j] - centroid[j]) for j in range(k)]
            else:
                con = [centroid[j] + 0.5 * (pts[k][j] - centroid[j]) for j in range(k)]
            fc = f(con)
            if fc < min(fr, vals[k]):
                pts[k], vals[k] = con, fc
            else:
                for i in range(1, k + 1):
                    pts[i] = [pts[0][j] + 0.5 * (pts[i][j] - pts[0][j]) for j in range(k)]
                    vals[i] = f(pts[i])
    best = min(range(k + 1), key=lambda i: (vals[i], i))
    return pts[best], vals[best], converged


def triers(active, annual_retention):
    """Customers who ever bought, rebuilt from active customers: n_t = A_t - rq * A_(t-1), N_0 = A_0."""
    rq = annual_retention ** 0.25
    out = [active[0]]
    for t in range(1, len(active)):
        out.append(out[-1] + active[t] - rq * active[t - 1])
    return out


def fit_curve(cum):
    """Least squares of the S-curve on t = 0 .. len(cum) - 1, in millions; z = (log alpha, beta, log gamma)."""
    scale = 1e6
    ys = [v / scale for v in cum]
    log_bound = math.log(ALPHA_BOUND_FACTOR * ys[-1])

    def sse(z):
        if z[0] > log_bound or z[2] > 3.0:
            return PENALTY
        a, b, g = math.exp(z[0]), z[1], math.exp(z[2])
        return plain_sum([(y - a / (1.0 + math.exp(-b - g * t))) ** 2 for t, y in enumerate(ys)])

    z = [math.log(2.0 * ys[-1]), -3.0, math.log(0.2)]
    value = sse(z)
    for _ in range(50):
        z_new, v_new, converged = nelder_mead(sse, z, 0.5)
        if not converged:
            raise SystemExit("error: the curve fit did not converge within NM_MAX_ITER iterations")
        moved = value - v_new
        z, value = z_new, v_new
        if moved < 1e-12:
            par = (math.exp(z[0]) * scale, z[1], math.exp(z[2]))
            return par, value, z[0] > log_bound - 1e-3
    raise SystemExit("error: the curve fit kept moving after 50 restarts")


# ---- synthetic data ---------------------------------------------------------------------------------------

def simulate(stream, par, quarters, annual_retention, noise):
    """Active customers A_0 .. A_quarters: A_0 = N(0); A_t = rq A_(t-1) + n_t exp(noise e - noise^2/2)."""
    rq = annual_retention ** 0.25
    active = [cumulative(par, 0)]
    for t in range(1, quarters + 1):
        n_t = cumulative(par, t) - cumulative(par, t - 1)
        shock = math.exp(noise * stream.normal() - 0.5 * noise * noise)
        active.append(math.floor(rq * active[-1] + n_t * shock + 0.5))
    active[0] = math.floor(active[0] + 0.5)
    return active


def lifetime_shortcut(annual_margin, annual_retention, annual_discount):
    """Expected lifetime 1 / (1 - r) years, margins discounted over that many years (Gupta et al., note 3)."""
    years = 1.0 / (1.0 - annual_retention)
    return annual_margin * (1.0 - (1.0 + annual_discount) ** (-years)) / annual_discount


def value_all(par, current, start, retention=RETENTION, discount=DISCOUNT, margin=MARGIN, cost=ACQUISITION_COST):
    cur, fut = base_value(par, current, start, margin, cost, retention, discount)
    return cur + fut


def quantile_pair(xs):
    """2.5th and 97.5th sorted values (the 25th and 975th of 1,000)."""
    s = sorted(xs)
    return s[round(0.025 * len(s)) - 1], s[round(0.975 * len(s)) - 1]


def noise_sd(resid, parameters):
    """Standard deviation of residuals, with one degree of freedom per fitted parameter."""
    m = plain_sum(resid) / len(resid)
    return math.sqrt(plain_sum([(x - m) ** 2 for x in resid]) / (len(resid) - parameters))


def log_residuals(cum, par):
    out = []
    for t in range(1, len(cum)):
        out.append(math.log((cum[t] - cum[t - 1]) / (cumulative(par, t) - cumulative(par, t - 1))))
    return out


def mean_sd(xs):
    m = plain_sum(xs) / len(xs)
    return m, math.sqrt(plain_sum([(x - m) ** 2 for x in xs]) / (len(xs) - 1))


# ---- invariants -------------------------------------------------------------------------------------------

def check_invariants():
    problems = []
    checked = 0
    par = (TRUE_ALPHA, TRUE_BETA, TRUE_GAMMA)
    # Simpson integral of the acquisition rate equals the remaining market, alpha - N(T).
    rq, iq = quarterly(RETENTION, 0.0)
    remaining = future_cohorts(par, 1.0, 0.0, QUARTERS)
    if abs(remaining - (TRUE_ALPHA - cumulative(par, QUARTERS))) > 1e-3:
        problems.append("Simpson integral of n(t)")
    checked += 1
    # Lifetime value closed form equals its series.
    rq, iq = quarterly(RETENTION, DISCOUNT)
    series = plain_sum([MARGIN * rq ** t / (1.0 + iq) ** t for t in range(1, 3001)])
    if abs(series - lifetime_value(MARGIN, rq, iq)) > 1e-9:
        problems.append("lifetime value series")
    checked += 1
    # Gupta et al., note 3: margin 100, retention 80 %, discount 12 % -> 250; five-year annuity -> 360.
    if abs(lifetime_value(100.0, 0.8, 0.12) - 250.0) > 1e-9:
        problems.append("note 3 lifetime value")
    checked += 1
    if abs(lifetime_shortcut(100.0, 0.8, 0.12) - 360.4776) > 1e-4:
        problems.append("note 3 expected-lifetime value")
    checked += 1
    # Their acquisition example: 100,000 then 130,000 customers at 80 % retention -> 50,000 acquired.
    if abs(triers([100000.0, 130000.0], 0.8 ** 4)[1] - 150000.0) > 1e-6:
        problems.append("acquisitions rebuilt from active customers")
    checked += 1
    # The acquisition peak is at t = -beta / gamma.
    peak = -TRUE_BETA / TRUE_GAMMA
    if not (acquisition_rate(par, peak) > acquisition_rate(par, peak - 0.01) and
            acquisition_rate(par, peak) > acquisition_rate(par, peak + 0.01)):
        problems.append("acquisition peak")
    checked += 1
    # Without noise and with the true retention, the rebuilt triers are the true curve.
    clean = simulate(Declared(1), par, QUARTERS, RETENTION, 0.0)
    rebuilt = triers(clean, RETENTION)
    if max(abs(rebuilt[t] - cumulative(par, t)) for t in range(QUARTERS + 1)) > 2.0 * QUARTERS:
        problems.append("triers rebuilt without noise")
    checked += 1
    if problems:
        raise SystemExit("error: invariant failed: " + ", ".join(problems))
    return checked


def published_checks():
    print("== Check on the published case: Gupta, Lehmann and Stuart, Tables 1 to 5 ==")
    print("convention: quarterly retention r^(1/4), quarterly discount 1.12^(1/4) - 1, t = 1 for the first quarter of data, tax 38 %")
    print(f"note 3: lifetime value {fmt(lifetime_value(100.0, 0.8, 0.12), 2)} (first margin after one period); "
          f"equation (2) summed from t = 0 would give {fmt(100.0 * 1.12 / (1.12 - 0.8), 2)}; "
          f"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; "
          f"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 = 0
    for row in GLS:
        name = row[0]
        v = gls_value(row)
        e_ret = gls_value(row, retention=row[4] * 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
        published = row[9:]
        ok_value = abs(v - published[0]) <= 0.0151
        ok_el = all(abs(100.0 * e - p) <= 0.0251 for e, p in zip((e_ret, e_acq, e_mar), published[1:4]))
        reproduced += ok_value and ok_el
        print(f"{name}: value {fmt(v, 3)} (published {fmt(published[0], 2)}, {'reproduced' if ok_value else 'NOT reproduced'}); "
              f"1 % better retention {pct(e_ret, 2)} % (published {fmt(published[1], 2)}), "
              f"acquisition cost {pct(e_acq, 2)} % ({fmt(published[2], 2)}), margin {pct(e_mar, 2)} % ({fmt(published[3], 2)}), "
              f"margin minus acquisition {fmt(100.0 * (e_mar - e_acq), 3)}; "
              f"discount rate x 0.99 {pct(e_dis, 2)} %, 12 % to 11.80 % {pct(e_dis_118, 2)} % ({fmt(published[4], 2)}); "
              f"retention to discount ratio {fmt(e_ret / e_dis, 1)}")
    print(f"firms whose value and three elasticities are reproduced: {reproduced} of 5")
    if reproduced != 4:
        raise SystemExit("error: the published check no longer reproduces four firms")
    cap = GLS[2]
    v_cap = gls_value(cap, future_mult=1.25)
    print(f"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)} %, "
          f"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 = []
    for row in GLS:
        name, peak = row[0], -row[6] / row[7]
        cells = []
        for label, offset in (("t = 1", 1), ("t = 0", 0)):
            y, mth = quarter_after(GLS_START[name], math.ceil(peak) - offset)
            cells.append(f"{label} {mth:02d}/{y}")
        y, mth = GLS_PEAK[name]
        peaks.append(f"{name} {fmt(peak, 2)}: {', '.join(cells)} (published {mth:02d}/{y})")
    print("acquisition peak, calendar quarter: " + "; ".join(peaks))
    print("Table 5 recomputed with the curve held fixed (published), discount 12 %:")
    lower = higher = same = 0
    for row in GLS:
        name, base = row[0], row[4]
        cells = []
        for d_i, d in enumerate((0.08, 0.12, 0.16)):
            for r_i, r in enumerate((0.70, 0.80, 0.90)):
                v = gls_value(row, retention=r, discount=d)
                pub = GLS_TABLE5[name][d_i][r_i]
                if d == 0.12:
                    cells.append(f"r {pct(r, 0)} % {fmt(v, 2)} ({fmt(pub, 2)})")
                if abs(r - base) < 1e-9:
                    same += abs(v - pub) <= 0.0151
                elif (r > base and pub < v) or (r < base and pub > v):
                    lower += 1
                else:
                    higher += 1
        full = GLS_FULL_RETENTION.get(name)
        extra = f"; retention 100 % {fmt(gls_value(row, retention=1.0), 2)} (text {fmt(full, 2)})" if full else ""
        print(f"  {name}: " + "; ".join(cells) + extra)
    print(f"Table 5 cells away from the firm's own retention: {lower} of {lower + higher} move the way a refitted curve would; "
          f"cells at the firm's own retention within 0.015: {same} of {45 - lower - higher}")
    print()


def main():
    print("MSC-P-045 reference output")
    print(f"invariants checked: {check_invariants()}")
    print()
    published_checks()

    par = (TRUE_ALPHA, TRUE_BETA, TRUE_GAMMA)
    stream = Declared(LCG_SEED)
    active = simulate(stream, par, QUARTERS, RETENTION, NOISE_SD)
    print("== Synthetic online retailer ==")
    print(f"design: seed {LCG_SEED}, acquisition noise sd {fmt(NOISE_SD, 2)}, Simpson step {fmt(STEP, 3)} quarter over {fmt(HORIZON, 0)} quarters, "
          f"bootstrap {BOOTSTRAP}, replications {REPLICATIONS}")
    print(f"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)}")
    print(f"annual retention {pct(RETENTION, 0)} %, quarterly margin {fmt(MARGIN, 2)}, acquisition cost {fmt(ACQUISITION_COST, 2)}, annual discount {pct(DISCOUNT, 0)} %")
    print("active customers by quarter: " + " ".join(fmt(a, 0) for a in active))
    rq, iq = quarterly(RETENTION, DISCOUNT)
    lv = lifetime_value(MARGIN, rq, iq)
    print(f"quarterly retention {fmt(rq, 4)}, quarterly discount {fmt(100.0 * iq, 3)} %, value of one customer {fmt(lv, 2)}")
    cum = triers(active, RETENTION)
    print(f"customers who ever bought, rebuilt: {fmt(cum[-1], 0)} (true {fmt(cumulative(par, QUARTERS), 0)})")
    est, sse, bound = fit_curve(cum)
    print(f"fitted curve: alpha {fmt(est[0], 0)}, beta {fmt(est[1], 3)}, gamma {fmt(est[2], 4)}; peak {fmt(-est[1] / est[2], 2)}; "
          f"SSE {fmt(sse, 6)} (millions squared); bound reached: {'yes' if bound else 'no'}")
    current = active[-1]
    cur_e, fut_e = base_value(est, current, QUARTERS, MARGIN, ACQUISITION_COST, RETENTION, DISCOUNT)
    cur_t, fut_t = base_value(par, current, QUARTERS, MARGIN, ACQUISITION_COST, RETENTION, DISCOUNT)
    ref = cur_t + fut_t
    value = cur_e + fut_e
    print(f"current customers {fmt(current, 0)}: value {fmt(cur_e / 1e6, 2)} M")
    print(f"future customers: estimated {fmt(fut_e / 1e6, 2)} M, true curve {fmt(fut_t / 1e6, 2)} M; "
          f"to acquire: estimated {fmt(est[0] - cumulative(est, QUARTERS), 0)}, true {fmt(TRUE_ALPHA - cumulative(par, QUARTERS), 0)}")
    print(f"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)} %); "
          f"future customers {pct(fut_e / value)} % of it")
    print()

    print("== Shortcuts ==")
    print(f"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)
    print(f"expected-lifetime shortcut: {fmt(1.0 / (1.0 - RETENTION), 1)} years, {fmt(short, 2)} per customer instead of {fmt(lv, 2)} "
          f"({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)
    print(f"same shortcut against annual values: {fmt(annual_next, 2)} with the first margin after one year ({pct(short / annual_next - 1.0)} %), "
          f"{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[:EARLY_QUARTERS + 1], RETENTION)
    est_e, sse_e, bound_e = fit_curve(early)
    v_early = value_all(est_e, active[EARLY_QUARTERS], EARLY_QUARTERS)
    ref_early = value_all(par, active[EARLY_QUARTERS], EARLY_QUARTERS)
    print(f"curve fitted on quarters 0 to {EARLY_QUARTERS}, before the peak: alpha {fmt(est_e[0], 0)}, beta {fmt(est_e[1], 3)}, gamma {fmt(est_e[2], 4)}; "
          f"bound reached: {'yes' if bound_e else 'no'}; value at quarter {EARLY_QUARTERS} {fmt(v_early / 1e6, 2)} M "
          f"(reference {fmt(ref_early / 1e6, 2)} M, {pct(v_early / ref_early - 1.0)} %)")
    print()

    print("== Elasticities (curve held fixed, as in Table 4) ==")
    print(f"improved retention: {pct(RETENTION * 1.01)} % instead of {pct(RETENTION, 0)} %")
    for label, kwargs in (("retention x 1.01", dict(retention=RETENTION * 1.01)),
                          ("margin x 1.01", dict(margin=MARGIN * 1.01)),
                          ("acquisition cost x 0.99", dict(cost=ACQUISITION_COST * 0.99)),
                          ("discount rate x 0.99", dict(discount=DISCOUNT * 0.99))):
        print(f"{label}: {pct(value_all(est, current, QUARTERS, **kwargs) / 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))
    print(f"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)}; "
          f"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)} %")
    print()

    print("== Retention and discount grid, value in M (curve held fixed / acquisitions rebuilt and curve refitted) ==")
    for d in (0.08, 0.10, 0.12):
        cells = []
        for r in (0.70, 0.80, 0.90):
            fixed = value_all(est, current, QUARTERS, retention=r, discount=d)
            refit, _, _ = fit_curve(triers(active, r))
            moved = value_all(refit, current, QUARTERS, retention=r, discount=d)
            cells.append(f"r {pct(r, 0)} %: {fmt(fixed / 1e6, 2)} / {fmt(moved / 1e6, 2)}")
        print(f"discount {pct(d, 0)} %: " + "; ".join(cells))
    for r in (0.70, 0.90):
        refit, _, _ = fit_curve(triers(active, r))
        print(f"assumed retention {pct(r, 0)} %: rebuilt customers who ever bought {fmt(triers(active, r)[-1], 0)}, fitted alpha {fmt(refit[0], 0)}")
    print()

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

    print("== 100 independent retailers drawn from the true curve ==")
    full_err, early_err, bound_hits, current_only, covered, early_hit = [], [], 0, [], 0, []
    for _ in range(REPLICATIONS):
        sim = simulate(stream, par, QUARTERS, RETENTION, NOISE_SD)
        sim_cum = triers(sim, RETENTION)
        f_par, _, _ = fit_curve(sim_cum)
        ref_s = value_all(par, sim[-1], QUARTERS)
        f_sd = noise_sd(log_residuals(sim_cum, f_par), 3)
        draws = []
        for _ in range(COVERAGE_BOOTSTRAP):
            b_par, _, _ = fit_curve(triers(simulate(stream, f_par, QUARTERS, RETENTION, f_sd), RETENTION))
            draws.append(value_all(b_par, sim[-1], QUARTERS))
        lo_b, hi_b = quantile_pair(draws)
        covered += lo_b <= ref_s <= hi_b
        full_err.append(value_all(f_par, sim[-1], QUARTERS) / ref_s - 1.0)
        cur_only, _ = base_value(f_par, sim[-1], QUARTERS, MARGIN, ACQUISITION_COST, RETENTION, DISCOUNT)
        current_only.append(cur_only / ref_s - 1.0)
        e_par, _, hit = fit_curve(triers(sim[:EARLY_QUARTERS + 1], RETENTION))
        bound_hits += hit
        early_hit.append(hit)
        early_err.append(value_all(e_par, sim[EARLY_QUARTERS], EARLY_QUARTERS) / value_all(par, sim[EARLY_QUARTERS], EARLY_QUARTERS) - 1.0)
    share = covered / REPLICATIONS
    print(f"bootstrap interval (nominal 95 %, {COVERAGE_BOOTSTRAP} draws each) contains the reference in {covered} of {REPLICATIONS} retailers "
          f"(simulation error {pct(math.sqrt(share * (1.0 - share) / REPLICATIONS))} %)")
    m, sd = mean_sd(full_err)
    print(f"full window: mean gap {pct(m)} % (sd {pct(sd)} %, simulation error {pct(sd / math.sqrt(REPLICATIONS))} %), largest {pct(max(abs(x) for x in full_err))} %")
    m, sd = mean_sd(current_only)
    print(f"current customers only: mean gap {pct(m)} % (sd {pct(sd)} %, simulation error {pct(sd / math.sqrt(REPLICATIONS))} %)")
    s = sorted(early_err)
    half = len(s) // 2
    within = sum(1 for x in early_err if abs(x) <= 0.20)
    free = sorted(x for x, h in zip(early_err, early_hit) if not h)
    print(f"curve fitted before the peak: median gap {pct((s[half - 1] + s[half]) / 2.0)} %; within 20 %: {within} of {REPLICATIONS}; "
          f"ceiling at its bound (1,000 times the customers counted): {bound_hits}; "
          f"other fits from {pct(free[0])} % to {pct(free[-1])} %")


if __name__ == "__main__":
    main()