Summary

Key papers (13)

1. The impact of heterogeneity in individual frailty on the dynamics of mortality. Vaupel, J.W., Manton, K.G. and Stallard, E. -- Demography 16(3):439-454, 1979.
   Origin of the frailty argument: with an unobserved multiplicative hazard factor the surviving population is progressively selected towards low-hazard types, so the observed aggregate hazard falls even when every individual hazard is flat. Transposed to telemarketing, this predicts a declining per-call conversion rate with zero contact fatigue, and it is the mechanism the paper must rule out before calling anything fatigue.
2. True and spurious duration dependence: the identifiability of the proportional hazard model. Elbers, C. and Ridder, G. -- Review of Economic Studies 49(3):403-409, 1982.
   Establishes that baseline duration dependence and the unobserved-heterogeneity distribution are separately identified in the mixed proportional hazard model provided observed regressors vary. This is the licence to attempt the fatigue-versus-selection decomposition at all, and it says that rich predetermined client covariates are the identifying resource.
3. A method for minimizing the impact of distributional assumptions in econometric models for duration data. Heckman, J.J. and Singer, B. -- Econometrica 52(2):271-320, 1984.
   Shows structural duration estimates are highly sensitive to the assumed mixing distribution and gives the consistent nonparametric maximum likelihood estimator: a discrete distribution on finitely many mass points. This is the estimator the study should use for the unobserved client type rather than a convenient gamma or normal frailty.
4. Heterogeneity and state dependence. Heckman, J.J. -- in Rosen, S. (ed.) Studies in Labor Markets, NBER/University of Chicago Press, 1981.
   Frames the general distinction between true state dependence, where the experience itself changes behaviour, and spurious state dependence generated by persistent unobserved differences. This is precisely the fatigue-versus-selection question, and it supplies the vocabulary and the warning that the two carry different policy implications.
5. Discrete-time methods for the analysis of event histories. Allison, P.D. -- Sociological Methodology 13:61-98, 1982.
   Standard reference for the person-period expansion and the estimation of a period hazard by binary regression with period dummies. Contact order is an integer, so the discrete-time form is exact rather than an approximation, and it turns the frailty mixture into a mixed binary-response model estimable by direct maximum likelihood.
6. Retention futility: targeting high-risk customers might be ineffective. Ascarza, E. -- Journal of Marketing Research 55(1):80-98, 2018.
   Closest conceptual antecedent. Two field experiments show that baseline risk and responsiveness to intervention are different things and that targeting on the former is ineffective. The dynamic analogue is that the profitable stopping point depends on whether decay is a property of the contact or of the surviving pool.
7. Efficiently evaluating targeting policies: improving on champion vs. challenger experiments. Simester, D., Timoshenko, A. and Zoumpoulis, S.I. -- Management Science 66(8):3412-3424, 2020.
   Sets the venue's standard for counterfactual policy evaluation: randomise to actions, then evaluate any policy off-policy. Observational campaign logs lack that randomisation, which is why a study on such logs must restrict itself to policies that only truncate the incumbent, where the counterfactual is exactly identified.
8. Interpretable optimal stopping. Ciocan, D.F. and Misic, V.V. -- Management Science 68(3):1616-1638, 2022.
   Argues that stopping policies a manager can read and audit -- tree policies -- approximate the optimum well. Supports delivering the result as a small readable table of client-specific contact caps rather than a black-box scoring function.
9. Managing churn to maximize profits. Lemmens, A. and Gupta, S. -- Marketing Science 39(5):956-973, 2020.
   Replaces the statistical loss with a profit-based loss so estimation matches the managerial objective. Provides the accounting discipline any stopping rule needs: the decision turns on the ratio of cost per contact to margin per conversion, and results must be reported as sensitivity to that ratio.
10. A data-driven approach to predict the success of bank telemarketing. Moro, S., Cortez, P. and Rita, P. -- Decision Support Systems 62:22-31, 2014.
    Source study for the Portuguese bank telemarketing data used here. It is a prediction paper benchmarking neural networks, SVM, trees and logistic regression, and it warns that call duration is known only after the call and cannot be used for targeting. Its cross-sectional finding that campaign frequency predicts failure is routinely misread as causal fatigue, which is exactly the confusion this study resolves.
11. Telephone call centers: tutorial, review, and research prospects. Gans, N., Koole, G. and Mandelbaum, A. -- Manufacturing & Service Operations Management 5(2):79-141, 2003.
    The canonical operations survey of call centres, concerned with staffing, routing and abandonment given an exogenous call list. Establishes that the number of redials per prospect -- a direct claim on agent capacity -- sits outside the standard operations model and flags the marketing-operations interface as under-developed.
12. Effect of temporal spacing between advertising exposures: evidence from online field experiments. Sahni, N.S. -- Quantitative Marketing and Economics 13(3):203-247, 2015.
    Randomised evidence that within-sequence exposure structure matters: spread-out ads convert better than bunched ones. Establishes that genuine state dependence in repeated contacting exists, but identifies spacing rather than contact order, in a setting where the marginal cost per exposure is negligible so no stopping problem binds.
13. A hidden Markov model of customer relationship dynamics. Netzer, O., Lattin, J.M. and Srinivasan, V. -- Marketing Science 27(2):185-204, 2008.
    The standard structural way marketing accommodates within-customer dynamics alongside persistent unobserved differences, via latent relationship states that shift with customer-firm encounters. Shows the modelling technology exists but has not been used to separate contact effects from pool composition or to produce a stopping rule.

1. Idea 1Selected

   Contact fatigue or dynamic selection? Decomposing the declining conversion hazard of a telemarketing campaign and the value of client-specific stopping rules

   Hypothesis. Most of the observed decline in the per-contact conversion rate across successive telemarketing calls is dynamic selection of the surviving prospect pool, not genuine contact fatigue. Formally, write the probability that client i converts on contact k, given no conversion on contacts 1..k-1, as h_ik = Lambda(delta_k + x_i'beta + theta_i), with Lambda the logistic function, delta_k a contact-order profile, x_i predetermined client attributes and theta_i a time-invariant unobserved type with a nonparametric discrete distribution. H1: measured on the linear-index scale, less than half of the raw decline in the observed log-odds hazard between contact 1 and contact 8 is attributable to delta_8 - delta_1; the remainder is the fall in E[x_i'beta | at risk at k] plus the fall in E[theta_i | at risk at k]. H1 is REFUTED if the fitted delta profile retains 50 percent or more of the raw index decline, i.e. if fatigue is the dominant mechanism. H2 (policy corollary): because selection is client-specific and fatigue is not, a contact cap that varies with predetermined client attributes raises campaign profit by materially more than the best single uniform cap. H2 is refuted if the best client-specific cap policy adds less than 10 percent on top of the profit gain of the best uniform cap.

   Dataset. Portuguese bank direct-marketing telephone campaigns bank_marketing_calls

   Test plan, novelty and risk

   Test plan. SAMPLE AND STRUCTURE. Each of the 41,188 rows is one client's whole contact spell in the campaign, not one call: `campaign` is the total number of calls made to that client and `y` is the terminal outcome. This is verifiable in the data because the support of `campaign` has gaps (36, 38, 44-55 are absent while 56 is present), which is impossible if each row were a single call in a sequence that must pass through every intermediate count. Treat each client as a right-censored discrete duration: event at k = campaign_i if y_i = yes; right-censored at k = campaign_i if y_i = no (the bank stopped calling). No rows are dropped; the analysed sample is all 41,188 spells, comprising sum(campaign_i) = 105,754 calls. STEP 1, RAW HAZARD. Life table over contact order k: at-risk R(k) = #{campaign_i >= k}, events E(k) = #{campaign_i = k and y_i = yes}, raw hazard hhat(k) = E(k)/R(k), with binomial standard errors. Report k = 1..10 and a pooled 11+ cell. STEP 2, PERSON-PERIOD EXPANSION. Expand each spell into campaign_i records (Allison 1982), giving 105,754 person-period observations, with the binary event indicator d_ik. STEP 3, MODELS, estimated by maximum likelihood with an analytic gradient and L-BFGS-B. (M0) contact-order dummies only, delta_1..delta_10 plus a pooled 11+ term: reproduces the raw hazard exactly. (M1) add PREDETERMINED client covariates x_i: age and age squared, job, marital, education, credit default, housing loan, personal loan, contact channel, number of prior-campaign contacts, previous-campaign outcome, and a never-previously-contacted indicator with the pdays value interacted. These are attributes of the client, fixed before the campaign starts. (M2) add a Heckman-Singer (1984) discrete frailty theta with G mass points, normalising theta_1 = 0 and estimating the remaining locations and the mixing probabilities as free parameters; choose G in {1,2,3,4} by BIC. Cluster-free by construction: the mixture likelihood is written at the SPELL level so each client contributes one independent term. Report standard errors from the numerical Hessian of the log-likelihood. STEP 4, DECOMPOSITION. On the linear-index scale the three sources are exactly additive: E[index | at risk at k] = delta_k + E[x_i'beta | at risk at k] + E[theta_i | at risk at k], where the last term uses the Bayes posterior over mass points given k-1 observed failures. Report, for the move from contact 1 to contact 8, the total index change and the three components in index units and as percentage shares. H1 concerns the share carried by delta. Also report the model-implied population hazard against the raw life-table hazard as a fit check. STEP 5, POLICY EVALUATION, exactly identified and free of extrapolation. Only policies that TRUNCATE the bank's observed policy are evaluated: under a cap kbar_i, client i receives min(campaign_i, kbar_i) calls and converts if and only if y_i = yes and campaign_i <= kbar_i. Both quantities are observed, so no counterfactual beyond the log is needed. Define profit = m * conversions - c * calls with rho = c/m the cost-to-margin ratio, and report profit relative to margin m so the result is unit-free. Compare (a) the bank's observed policy, (b) the best uniform cap kbar, (c) the best client-specific cap kbar(x_i) obtained by applying the myopic optimal rule -- call again if and only if the model's posterior hazard is at least rho -- with the posterior computed from predetermined covariates only, and (d) as an honest benchmark, the best cap chosen separately within deciles of the fitted baseline propensity. Report over rho in {0.005, 0.01, 0.02, 0.03, 0.05}, i.e. cost per call between 0.5 and 5 percent of the margin on a conversion. Split the sample 50/50 by a seeded random split, fit the policy on the training half and report profit on the held-out half so the client-specific cap is not evaluated in sample. STEP 6, ROBUSTNESS. (i) Re-estimate adding the calendar month and the five macro indicators, which in these data are stamped at the LAST contact and are therefore not predetermined; report how the decomposition shares move, since this bounds the contribution of calendar composition. (ii) Re-estimate with a complementary log-log link. (iii) Re-estimate on spells with campaign <= 12 to check that the long right tail is not driving the result. EXCLUSIONS. `duration` is realised only after the call ends and is excluded from every model, following the source study. REFUTATION CRITERIA are stated in the hypothesis: H1 fails if |delta_8 - delta_1| is at least half the raw index decline; H2 fails if the client-specific cap adds less than 10 percent on top of the uniform-cap gain.

   Novelty. The frailty literature has known since Vaupel, Manton and Stallard (1979) that a declining population hazard is the mechanical signature of unobserved heterogeneity even with no individual-level decline, and Elbers and Ridder (1982) proved the two are separately identified when observed regressors vary. Neither result has been brought to bear on the contact-frequency decision, where practice reads response decay across successive touches as fatigue and prescribes a uniform frequency cap. The INFORMS targeting literature has settled the extensive margin -- Ascarza (2018) shows to target on responsiveness rather than baseline risk, Lemmens and Gupta (2020) align the loss with profit -- but the intensive dynamic margin, how long to persist with a given prospect, has had no comparable treatment. This dataset has been reused almost exclusively as an imbalanced-classification leaderboard, and Moro, Cortez and Rita (2014) is explicitly a prediction study whose 'campaign frequency matters' finding is a cross-sectional association. Reading the same file as 41,188 right-censored contact spells rather than 41,188 labelled cross-sections is what makes the decomposition possible, and it turns a familiar benchmark into evidence on a decision. The two mechanisms prescribe opposite actions -- fatigue implies one cap for everyone, selection implies persistence that varies by client -- so the decomposition is not a modelling nicety but the thing that determines the policy.

   Expected outcome. The raw hazard falls from roughly 5.6 percent at contact 1 to roughly 1.6 percent at contact 7, a factor of about three and a half, so there is a large effect to decompose. I expect selection to carry the majority of it, with observable composition -- much of it the calendar and macro environment in which long spells occur -- explaining a visible slice and unobserved frailty more, and a genuine but smaller residual fatigue term of perhaps a third of the index decline. For policy I expect the bank to be over-contacting at any plausible cost ratio: the best uniform cap should land somewhere in the region of 2 to 5 calls at rho between 0.01 and 0.05, and should recover a few percent of margin-normalised profit. The client-specific cap should beat it, but by a modest amount, because a cap is a coarse instrument and the predetermined covariates in this file are weak predictors. A refutation is entirely possible and would be equally publishable: if the fitted delta profile keeps most of the raw decline, the correct reading is that telemarketing fatigue is real and a uniform cap is the right instrument, which contradicts the frailty prior and vindicates standard practice.

   Main risk. Identification of the frailty distribution. Elbers and Ridder (1982) guarantee separate identification given regressor variation, but the spells here are short -- the mean is 2.57 calls -- and the predetermined covariates are weak, so the likelihood may be flat and the split between delta and theta may rest heavily on the logistic functional form. Mitigation: report the observables-only decomposition from M1, which needs no frailty assumption at all, as a floor on the selection share; report the frailty split as the additional, more parametric step; and report BIC across G so the reader sees how much the data actually support. Second risk: endogenous censoring. The bank chose when to stop, and if it stopped on information correlated with theta the frailty estimate is contaminated. This is why the policy evaluation is restricted to truncations of the observed policy, which are exactly identified whatever generated the censoring. It must be stated as a maintained assumption for the decomposition and not for the policy result. Third risk: the calendar and macro fields are stamped at the last contact, so treating them as predetermined would be wrong; they are therefore held out of the main specification and used only as a robustness bound.

2. Idea 2

   Does the macroeconomic environment change whom a bank should call, or only how many will say yes? A targeting-policy stability test across the 2008-2010 Portuguese downturn

   Hypothesis. The macroeconomic cycle shifts the LEVEL of telemarketing conversion but leaves the RANKING of prospects essentially unchanged, so a targeting policy trained in one macro regime transfers with little loss to another. H: the profit loss from training the targeting rule on the high-euribor regime and deploying it in the low-euribor regime, relative to training and deploying within the low-euribor regime, is under 10 percent of the achievable profit gain over random targeting. REFUTED if the transfer loss exceeds 10 percent, which would mean the bank must re-estimate its call list as conditions change rather than treat the score as stable.

   Dataset. Portuguese bank direct-marketing telephone campaigns bank_marketing_calls

   Test plan, novelty and risk

   Test plan. Split the 41,188 spells by euribor3m at its median into a high-rate regime (roughly the 2008-2009 calls) and a low-rate regime (roughly 2010). Within each regime, fit an L2-penalised logistic model of y on predetermined client attributes only -- age, job, marital, education, default, housing, loan, contact channel, previous, poutcome, pdays sentinel -- explicitly excluding duration and all macro variables so the score is a client score, not a calendar score. Use 5-fold cross-validation with a fixed seed to pick the penalty. Evaluate four combinations of train regime and deploy regime. In each deployment sample, rank clients by score, call the top q percent for q on a grid from 5 to 100, and compute profit = m * conversions - c * calls with rho = c/m in {0.02, 0.05, 0.10}. Report the best-q profit for the within-regime policy and for the transferred policy, and the transfer loss as a percentage of the within-regime gain over calling everyone. Add AUC and top-decile lift in each cell for reference, and bootstrap the transfer loss with 200 seeded resamples for a confidence interval. Report a coefficient-stability table -- the two regimes' coefficients with standard errors and their correlation -- to show whether any transfer loss comes from re-ranking or only from level shifts.

   Novelty. The bank marketing data is normally used to report a single AUC on a single random split, which cannot detect that the conversion base rate roughly triples between the high-rate and low-rate periods. Simester, Timoshenko and Zoumpoulis (2020) show that targeting policies degrade badly under realistic data challenges; regime shift is one such challenge and this dataset contains a genuine, large one -- a sovereign-debt-era swing in Portuguese interest rates and employment. The question of whether a targeting score has to be re-estimated when the environment moves, or only recalibrated, is a live operational question and is answerable here.

   Expected outcome. I expect large level shifts -- conversion is far higher in the low-rate months -- but a fairly stable ranking, so a transfer loss of well under 10 percent and the hypothesis supported. The interesting quantitative output is the size of the recalibration needed: the same score implies a very different optimal call depth q in the two regimes.

   Main risk. The regime split is highly collinear with calendar time and with campaign design changes the data does not record, so 'macro regime' is really 'period'. The finding would then be about temporal stability rather than macroeconomic sensitivity, which is a weaker and less distinctive claim. It is also a comparatively incremental question: it tests the stability of a standard targeting exercise rather than resolving a mechanism.

3. Idea 3

   How much reserve capacity does a distributional load forecast save over a point forecast? A newsvendor evaluation of half-hourly electricity demand in Victoria

   Hypothesis. A quantile forecast of half-hourly demand dominates the common practice of a point forecast plus a fixed proportional safety margin, at equal reliability. H: at a 98 percent no-shortfall service level, scheduling capacity at the predictive 98th percentile from a temperature-and-seasonality quantile regression requires at least 3 percent less average capacity than a point forecast inflated by whatever fixed multiplier delivers the same realised 98 percent reliability out of sample. REFUTED if the average capacity saving is under 3 percent, meaning the fixed-margin heuristic is already near-efficient and the extra modelling buys nothing.

   Dataset. Half-hourly electricity demand for Victoria, Australia vic_electricity_demand

   Test plan, novelty and risk

   Test plan. Convert the UTC timestamps to Australia/Melbourne local time before extracting time-of-day, following the dataset's own caveat. Train on 2012-2013 and test on 2014, a strict temporal holdout. Build features: half-hour-of-day and day-of-week indicators, annual Fourier terms, the public-holiday flag, and a piecewise-linear temperature spline with knots near 15 and 22 degrees to capture the U-shaped heating and cooling response. Fit (a) OLS for the conditional mean and (b) statsmodels QuantReg at the 0.98 quantile on the same design. On 2014, evaluate two capacity policies at matched realised reliability: policy A schedules the 0.98 quantile prediction; policy B schedules (1 + s) times the OLS prediction with s calibrated on the training years to hit 98 percent, then held fixed out of sample. Report realised shortfall frequency, mean scheduled capacity, mean capacity per avoided shortfall, and the pinball loss of each policy at the 0.98 quantile. Repeat across service levels 0.95, 0.98 and 0.99 as a sensitivity table, and report the newsvendor cost with shortage-to-holding ratios of 19:1, 49:1 and 99:1 matching those service levels.

   Novelty. The nonlinear temperature-demand response is well documented and would be a trivial finding on its own. The decision-relevant quantity -- how much scheduled capacity a distributional forecast saves at matched reliability, which is what a system operator actually pays for -- is a different question and is answerable directly from three complete years of half-hourly operational demand.

   Expected outcome. A saving in the low single digits of percent, which on a 4,665 unit mean demand is operationally meaningful. The quantile policy should be visibly better in summer, where the cooling response makes the conditional variance heteroskedastic and a fixed proportional margin is badly miscalibrated.

   Main risk. The fixed-margin benchmark may be a straw man; the honest comparison must calibrate s on the training data to hit exactly the same target reliability, and even then the result may reduce to 'heteroskedasticity exists'. The demand unit is disputed at source, so results must be reported in percentage terms or as a demand index. This is also the least distinctive of the three: it is a competent forecasting evaluation rather than a decomposition that changes what a manager should do.

Ranking rationale

Idea 1 is first because it is the only one of the three whose answer changes a decision rather than sharpening an estimate. Two mechanisms -- contact fatigue and dynamic selection -- generate exactly the same declining response curve and prescribe opposite policies: one cap for everyone, or persistence that varies client by client. The decomposition that separates them is well posed (Elbers and Ridder 1982), estimable here (Heckman and Singer 1984 discrete frailty on a person-period expansion), and it has never been applied to the contact-frequency decision. It also has a genuinely falsifiable pre-stated threshold, a formal model with a monotone stopping result, and a policy evaluation that needs no extrapolation because it only truncates the observed policy. The raw hazard falls from about 5.6 to about 1.6 percent across contacts 1 to 7, so there is a large, real effect to decompose, and a refutation would be as informative as a confirmation. Idea 2 is second. The regime-transfer question is real and the data contains a genuine macroeconomic swing, but the answer is likely to be the unsurprising one -- rankings are stable, levels are not -- and the regime split is confounded with calendar time and unrecorded campaign changes, so the claim would have to be hedged into something close to a stability check on a standard targeting exercise. Idea 3 is third. It is well specified and would run cleanly, but it is a forecasting evaluation rather than a structural question, its benchmark is easy to make a straw man of, and its likely finding -- that a distributional forecast beats a fixed margin under heteroskedasticity -- is close to known. The disputed demand unit adds a reporting nuisance without adding insight.

Data provenance

Result (RESULT_JSON)

{
  "n_spells": 41188,
  "n_calls": 105754,
  "n_conversions": 4640,
  "n_person_periods": 105754,
  "n_exposure_rows": 101237,
  "conversion_rate": 0.11265417111780132,
  "mean_calls_per_spell": 2.567592502670681,
  "max_contacts": 56,
  "support_gaps": [
    36,
    38,
    44,
    45,
    46,
    47,
    48,
    49,
    50,
    51,
    52,
    53,
    54,
    55
  ],
  "n_covariates": 32,
  "life_table": [
    {
      "k": 1,
      "n_risk": 41188,
      "exposures": 41188,
      "events": 2300,
      "hazard": 0.05584150723511702,
      "se": 0.0011313986753706924
    },
    {
      "k": 2,
      "n_risk": 23546,
      "exposures": 23546,
      "events": 1211,
      "hazard": 0.05143124097511254,
      "se": 0.0014394259916287425
    },
    {
      "k": 3,
      "n_risk": 12976,
      "exposures": 12976,
      "events": 574,
      "hazard": 0.04423551171393342,
      "se": 0.0018050554780315116
    },
    {
      "k": 4,
      "n_risk": 7635,
      "exposures": 7635,
      "events": 249,
      "hazard": 0.03261296660117878,
      "se": 0.0020327817972068613
    },
    {
      "k": 5,
      "n_risk": 4984,
      "exposures": 4984,
      "events": 120,
      "hazard": 0.024077046548956663,
      "se": 0.0021713026163695235
    },
    {
      "k": 6,
      "n_risk": 3385,
      "exposures": 3385,
      "events": 75,
      "hazard": 0.022156573116691284,
      "se": 0.0025299190139136648
    },
    {
      "k": 7,
      "n_risk": 2406,
      "exposures": 2406,
      "events": 38,
      "hazard": 0.015793848711554447,
      "se": 0.002541787343340507
    },
    {
      "k": 8,
      "n_risk": 1777,
      "exposures": 1777,
      "events": 17,
      "hazard": 0.00956668542487338,
      "se": 0.00230913674862167
    },
    {
      "k": 9,
      "n_risk": 1377,
      "exposures": 1377,
      "events": 17,
      "hazard": 0.012345679012345678,
      "se": 0.0029757264447003587
    },
    {
      "k": 10,
      "n_risk": 1094,
      "exposures": 1094,
      "events": 12,
      "hazard": 0.010968921389396709,
      "se": 0.003149040675856572
    },
    {
      "k": 11,
      "n_risk": 869,
      "exposures": 5386,
      "events": 27,
      "hazard": 0.005012996658002228,
      "se": 0.00096233047007039,
      "note": "pooled 11+"
    }
  ],
  "raw_logodds_change_1_to_8": -1.8120791471679305,
  "raw_hazard_ratio_8_to_1": 0.17131853881725428,
  "model_selection": {
    "1": {
      "loglik": -16609.29250949701,
      "npar": 43,
      "bic": 33675.49881491741
    },
    "2": {
      "loglik": -16598.503296113326,
      "npar": 45,
      "bic": 33675.1721926116
    },
    "3": {
      "loglik": -16598.503177971157,
      "npar": 47,
      "bic": 33696.42376078881
    },
    "4": {
      "loglik": -16598.28941825952,
      "npar": 49,
      "bic": 33717.248045827095
    }
  },
  "G_selected": 2,
  "lr_frailty": 21.57842676736618,
  "theta": [
    0,
    -14.64597927922687
  ],
  "pi": [
    0.9501411511101232,
    0.04985884888987669
  ],
  "type_first_hazard": [
    0.04243860970980862,
    1.931641511413968e-8
  ],
  "delta": [
    -3.1163312799342195,
    -3.1228936984152593,
    -3.184973576360359,
    -3.443023927986927,
    -3.7128245548722663,
    -3.756673364568039,
    -4.04256779762513,
    -4.509589376062568,
    -4.2164971806799985,
    -4.322440622222365,
    -5.049244383519336
  ],
  "delta_se": [
    0.030085338338469908,
    0.03809432166895132,
    0.05063066733161005,
    0.07147443025366448,
    0.09835922373271262,
    0.12230718579162984,
    0.16799199430866327,
    0.2467506891011131,
    0.24650390154841156,
    0.29265153164149166,
    0.1959109040505915
  ],
  "delta_train_half": [
    -3.1823927140131536,
    -3.177279929154025,
    -3.276782268696832,
    -3.475430886877544,
    -3.6328090280543095,
    -3.818395993126811,
    -4.213630728769942,
    -4.381153763590611,
    -4.327263191351745,
    -4.574479116696017,
    -5.416150735589605
  ],
  "delta_8_minus_1": -1.3932580961283485,
  "delta_8_minus_1_se": 0.2465568482245579,
  "hazard_profile": [
    {
      "k": 1,
      "n_risk": 41188,
      "delta": -3.1163312799342195,
      "Exb": 9.93669553748222e-17,
      "Etheta": -0.7302316677272378,
      "h_model": 0.05590107883430747,
      "h_no_fatigue": 0.05590107883430747,
      "h_no_selection": 0.05590107883430747
    },
    {
      "k": 2,
      "n_risk": 23546,
      "delta": -3.1228936984152593,
      "Exb": -0.05041475458181874,
      "Etheta": -0.7773806897854391,
      "h_model": 0.051342565498720114,
      "h_no_fatigue": 0.05163262304011857,
      "h_no_selection": 0.055593527303742035
    },
    {
      "k": 3,
      "n_risk": 12976,
      "delta": -3.184973576360359,
      "Exb": -0.10046472686702752,
      "Etheta": -0.822113966002176,
      "h_model": 0.044020732945149306,
      "h_no_fatigue": 0.0467868045707629,
      "h_no_selection": 0.05275893585516027
    },
    {
      "k": 4,
      "n_risk": 7635,
      "delta": -3.443023927986927,
      "Exb": -0.13248259540252638,
      "Etheta": -0.8605283647451198,
      "h_model": 0.03247732851114702,
      "h_no_fatigue": 0.04379512686945525,
      "h_no_selection": 0.04232833091011745
    },
    {
      "k": 5,
      "n_risk": 4984,
      "delta": -3.7128245548722663,
      "Exb": -0.1620027409536339,
      "Etheta": -0.8893865853359114,
      "h_model": 0.02392951891196156,
      "h_no_fatigue": 0.041722235840721975,
      "h_no_selection": 0.03346856014317972
    },
    {
      "k": 6,
      "n_risk": 3385,
      "delta": -3.756673364568039,
      "Exb": -0.1928762010900783,
      "Etheta": -0.9125725134033048,
      "h_model": 0.021924215880556035,
      "h_no_fatigue": 0.03999647271645393,
      "h_no_selection": 0.032201609450803796
    },
    {
      "k": 7,
      "n_risk": 2406,
      "delta": -4.04256779762513,
      "Exb": -0.22029066868119437,
      "Etheta": -0.9219997556876742,
      "h_model": 0.015834730135287474,
      "h_no_fatigue": 0.03827655395988403,
      "h_no_selection": 0.024969491396746136
    },
    {
      "k": 8,
      "n_risk": 1777,
      "delta": -4.509589376062568,
      "Exb": -0.24701047872120757,
      "Etheta": -0.9181621815429748,
      "h_model": 0.009605446641942884,
      "h_no_fatigue": 0.036882818150365815,
      "h_no_selection": 0.016323764268820436
    },
    {
      "k": 9,
      "n_risk": 1377,
      "delta": -4.2164971806799985,
      "Exb": -0.27048363007992265,
      "Etheta": -0.9069207014316388,
      "h_model": 0.012376278839296803,
      "h_no_fatigue": 0.03586115433992539,
      "h_no_selection": 0.0213420050235515
    },
    {
      "k": 10,
      "n_risk": 1094,
      "delta": -4.322440622222365,
      "Exb": -0.28215306793726946,
      "Etheta": -0.9184312562863124,
      "h_model": 0.010989353770426218,
      "h_no_fatigue": 0.03536163734014182,
      "h_no_selection": 0.019380716789649716
    }
  ],
  "decomposition_k8": {
    "k": 8,
    "index_change_total": -1.828199088665293,
    "index_change_fatigue": -1.3932580961283485,
    "index_change_observable": -0.24701047872120768,
    "index_change_unobservable": -0.1879305138157371,
    "share_fatigue": 0.7620932013184184,
    "share_observable": 0.13511136738479712,
    "share_unobservable": 0.10279543129678446
  },
  "decomposition_k5": {
    "k": 5,
    "index_change_total": -0.9176509335003548,
    "index_change_fatigue": -0.5964932749380472,
    "index_change_observable": -0.162002740953634,
    "index_change_unobservable": -0.15915491760867362,
    "share_fatigue": 0.6500219780333462,
    "share_observable": 0.17654070304890218,
    "share_unobservable": 0.17343731891775174
  },
  "decomposition_k8_observables_only": {
    "k": 8,
    "index_change_total": -1.693800555150124,
    "index_change_fatigue": -1.4519258300134084,
    "index_change_observable": -0.2418747251367158,
    "index_change_unobservable": 0,
    "share_fatigue": 0.8571999965395701,
    "share_observable": 0.1428000034604299,
    "share_unobservable": 0
  },
  "share_fatigue_ci95_numerator": [
    0.4977612554618328,
    1.026425147175004
  ],
  "pure_selection_model": {
    "G_selected": 3,
    "loglik": -16677.723066620056,
    "bic": 33748.60451577885,
    "best_loglik": -16677.278045835606,
    "lr_vs_unrestricted": 157.54949944456166,
    "hazard": [
      0.0576504363895491,
      0.04592661253760351,
      0.037673760012234526,
      0.0322795157241225,
      0.028426088085783424,
      0.025339242077410347,
      0.022637184442458554,
      0.02050173207947063,
      0.01882099124720847,
      0.017724367540671034
    ],
    "h8_pure_selection": 0.02050173207947063,
    "h8_observed": 0.00956668542487338,
    "bic_table": {
      "2": 33756.69721342004,
      "3": 33748.60451577885,
      "4": 33768.966278671505,
      "5": 33790.218085608576
    }
  },
  "myopic_path_monotone_share_k1_8": 0.06084296397008838,
  "myopic_path_monotone_share_k1_20": 0.021851024570263184,
  "policy": [
    {
      "rho": 0.005,
      "best_uniform_cap": 12,
      "decile_caps": [
        6,
        5,
        7,
        17,
        13,
        9,
        9,
        10,
        12,
        11
      ],
      "targeting_cap": 12,
      "mean_myopic_cap": 13.72773623385452,
      "share_never_called": 0,
      "profit_observed": 0.10303243663202874,
      "profit_uniform": 0.10315018937554628,
      "profit_myopic": 0.1028704962610469,
      "profit_decile": 0.10237326405749247,
      "profit_targeting": 0.10310794406137708,
      "profit_target_cap": 0.10315018937554628,
      "calls_observed": 2.555598718073225,
      "calls_uniform": 2.4640672040400116,
      "calls_myopic": 2.42289016218316,
      "conv_observed": 0.11581043022239489,
      "conv_uniform": 0.11547052539574632,
      "conv_myopic": 0.1149849470719627,
      "gain_uniform": 0.00011775274351753796,
      "gain_myopic": -0.00016194037098184175,
      "gain_myopic_pct_of_observed": -0.15717416405495438,
      "gain_uniform_pct_of_observed": 0.11428706081957617,
      "myopic_over_uniform": -2.3752577319586834,
      "ci_myopic_minus_uniform": [
        -0.0006531210546761192,
        0.00013176167815868615
      ],
      "ci_uniform_minus_observed": [
        -0.0001521316888414101,
        0.0003562323977857626
      ]
    },
    {
      "rho": 0.01,
      "best_uniform_cap": 9,
      "decile_caps": [
        4,
        4,
        6,
        6,
        5,
        9,
        9,
        10,
        12,
        11
      ],
      "targeting_cap": 9,
      "mean_myopic_cap": 8.45192774594542,
      "share_never_called": 0.030494318733611733,
      "profit_observed": 0.0902544430416626,
      "profit_uniform": 0.09063173739924248,
      "profit_myopic": 0.09045401573273767,
      "profit_decile": 0.09031174128386908,
      "profit_targeting": 0.0903704962610469,
      "profit_target_cap": 0.09058560745848306,
      "calls_observed": 2.555598718073225,
      "calls_uniform": 2.401330484607167,
      "calls_myopic": 2.220015538506361,
      "conv_observed": 0.11581043022239489,
      "conv_uniform": 0.11464504224531416,
      "conv_myopic": 0.11265417111780132,
      "gain_uniform": 0.0003772943575798826,
      "gain_myopic": 0.00019957269107508113,
      "gain_myopic_pct_of_observed": 0.2211222897724335,
      "gain_uniform_pct_of_observed": 0.4180341098617369,
      "myopic_over_uniform": -0.4710424710424496,
      "ci_myopic_minus_uniform": [
        -0.0007763547635233561,
        0.00047419151209090007
      ],
      "ci_uniform_minus_observed": [
        -0.00011954938331552856,
        0.0008623992424978149
      ]
    },
    {
      "rho": 0.02,
      "best_uniform_cap": 6,
      "decile_caps": [
        0,
        0,
        2,
        2,
        4,
        8,
        4,
        7,
        9,
        7
      ],
      "targeting_cap": 6,
      "mean_myopic_cap": 5.150626396037681,
      "share_never_called": 0.192289016218316,
      "profit_observed": 0.06469845586093036,
      "profit_uniform": 0.06748276196950567,
      "profit_myopic": 0.06995338448091677,
      "profit_decile": 0.06914538215014081,
      "profit_targeting": 0.06777896474701368,
      "profit_target_cap": 0.06919782460910943,
      "calls_observed": 2.555598718073225,
      "calls_uniform": 2.2682820238904537,
      "calls_myopic": 1.7077304069146353,
      "conv_observed": 0.11581043022239489,
      "conv_uniform": 0.11284840244731476,
      "conv_myopic": 0.10410799261920949,
      "gain_uniform": 0.0027843061085753157,
      "gain_myopic": 0.005254928619986415,
      "gain_myopic_pct_of_observed": 8.12218552987093,
      "gain_uniform_pct_of_observed": 4.3035124587211095,
      "myopic_over_uniform": 0.8873386815486597,
      "ci_myopic_minus_uniform": [
        0.0012267893561231432,
        0.003664343983684568
      ],
      "ci_uniform_minus_observed": [
        0.001863236865106342,
        0.003643876857337089
      ]
    },
    {
      "rho": 0.03,
      "best_uniform_cap": 4,
      "decile_caps": [
        0,
        0,
        0,
        0,
        4,
        4,
        4,
        5,
        5,
        7
      ],
      "targeting_cap": 5,
      "mean_myopic_cap": 3.3229095853161117,
      "share_never_called": 0.41079926192094784,
      "profit_observed": 0.039142468680198106,
      "profit_uniform": 0.0464144896571817,
      "profit_myopic": 0.05630863358259687,
      "profit_decile": 0.05634019617364281,
      "profit_targeting": 0.05395940565213168,
      "profit_target_cap": 0.05560600174808196,
      "calls_observed": 2.555598718073225,
      "calls_uniform": 2.0655530737107894,
      "calls_myopic": 1.2113236865106345,
      "conv_observed": 0.11581043022239489,
      "conv_uniform": 0.1083810818685054,
      "conv_myopic": 0.0926483441779159,
      "gain_uniform": 0.007272020976983595,
      "gain_myopic": 0.017166164902398764,
      "gain_myopic_pct_of_observed": 43.85560104205437,
      "gain_uniform_pct_of_observed": 18.578340156308176,
      "myopic_over_uniform": 1.3605769230769216,
      "ci_myopic_minus_uniform": [
        0.008068818587938235,
        0.011703930756531028
      ],
      "ci_uniform_minus_observed": [
        0.005892675050985724,
        0.00870026464018646
      ]
    },
    {
      "rho": 0.05,
      "best_uniform_cap": 2,
      "decile_caps": [
        0,
        0,
        0,
        0,
        0,
        1,
        0,
        3,
        2,
        7
      ],
      "targeting_cap": 5,
      "mean_myopic_cap": 1.4585316111488784,
      "share_never_called": 0.6423230067009809,
      "profit_observed": -0.011969505681266404,
      "profit_uniform": 0.009442070505972608,
      "profit_myopic": 0.03691366417403127,
      "profit_decile": 0.03600563270855588,
      "profit_targeting": 0.033189278430610855,
      "profit_target_cap": 0.03548606390210741,
      "calls_observed": 2.555598718073225,
      "calls_uniform": 1.5679809653297077,
      "calls_myopic": 0.6465960959502768,
      "conv_observed": 0.11581043022239489,
      "conv_uniform": 0.087841118772458,
      "conv_myopic": 0.06924346897154511,
      "gain_uniform": 0.02141157618723901,
      "gain_myopic": 0.04888316985529767,
      "gain_myopic_pct_of_observed": null,
      "gain_uniform_pct_of_observed": null,
      "myopic_over_uniform": 1.2830252863136409,
      "ci_myopic_minus_uniform": [
        0.02486652665824998,
        0.02991902981450908
      ],
      "ci_uniform_minus_observed": [
        0.01890235019908712,
        0.02398052830921628
      ]
    }
  ],
  "policy_main_rho": 0.02,
  "robustness": {
    "with_calendar_macro": {
      "k": 8,
      "index_change_total": -1.4658660616107997,
      "index_change_fatigue": -0.9418878489227804,
      "index_change_observable": -0.4641685888836087,
      "index_change_unobservable": -0.05980962380441068,
      "share_fatigue": 0.6425470058893142,
      "share_observable": 0.31665143292392395,
      "share_unobservable": 0.04080156118676186,
      "n_cov": 50,
      "loglik": -15267.100402553877
    },
    "cloglog": {
      "k": 8,
      "index_change_total": -1.7344932773098722,
      "index_change_fatigue": -1.3783341335195671,
      "index_change_observable": -0.23065599674952408,
      "index_change_unobservable": -0.12550314704078114,
      "share_fatigue": 0.7946609834414041,
      "share_observable": 0.1329817761572775,
      "share_unobservable": 0.07235724040131845,
      "loglik": -16628.922570442544
    },
    "campaign_le_12": {
      "k": 8,
      "index_change_total": -1.4945005401927591,
      "index_change_fatigue": -1.0516407844905298,
      "index_change_observable": -0.20647830805887005,
      "index_change_unobservable": -0.23638144764335933,
      "share_fatigue": 0.7036737399605691,
      "share_observable": 0.13815873765575132,
      "share_unobservable": 0.1581675223836795,
      "n_spells": 40621
    },
    "G_3": {
      "k": 8,
      "index_change_total": -1.8279993513881565,
      "index_change_fatigue": -1.3926602330453388,
      "index_change_observable": -0.247082043749267,
      "index_change_unobservable": -0.18825707459355057,
      "share_fatigue": 0.7618494131235728,
      "share_observable": 0.13516527976973103,
      "share_unobservable": 0.1029853071066961
    },
    "G_4": {
      "k": 8,
      "index_change_total": -1.7845494368957269,
      "index_change_fatigue": -1.3786025564404194,
      "index_change_observable": -0.2490846204031238,
      "index_change_unobservable": -0.15686226005218373,
      "share_fatigue": 0.772521359138437,
      "share_observable": 0.1395784365808399,
      "share_unobservable": 0.08790020428072308
    }
  },
  "hypothesis_H1_selection_dominates": false,
  "hypothesis_H2_personalisation_gain": true,
  "hypothesis_supported": false,
  "runtime_seconds": 61.696630001068115
}

"""Contact fatigue vs dynamic selection in a real telemarketing campaign.

Discrete-time conversion hazard over contact order with Heckman-Singer discrete
frailty; exactly additive index-scale decomposition; a constrained pure-selection
model as an identification check; truncation-only stopping-policy evaluation.
"""
import json, time, numpy as np, pandas as pd
from scipy.optimize import minimize
from scipy.special import expit, logsumexp
from dataset import load

T0 = time.time()
SEED = 20240917
KMAX = 11          # contact-order cells: 1..10 singly, 11 = pooled 11+
KPOL = 20          # policy horizon

df = load()
assert list(df.shape) == [41188, 21] and df.isna().sum().sum() == 0
K_all = df["campaign"].to_numpy(np.int64)
d_all = (df["y"].to_numpy() == "yes").astype(float)
gaps = [k for k in range(1, K_all.max() + 1) if k not in set(np.unique(K_all).tolist())]
print(f"n={len(df)} calls={K_all.sum()} conv={int(d_all.sum())} gaps={gaps}")


def build_X(d, macro=False):
    d = d.copy()
    d.loc[d["default"] == "yes", "default"] = "unknown"             # 3 spells: separation
    d.loc[d["education"] == "illiterate", "education"] = "unknown"  # 18 spells
    num = pd.DataFrame(index=d.index)
    num["age"] = d["age"].astype(float)
    num["age_sq"] = d["age"].astype(float) ** 2
    num["previous"] = d["previous"].astype(float)
    num["pdays_recent"] = (d["pdays"] < 999).astype(float)
    num["pdays_val"] = np.where(d["pdays"] < 999, d["pdays"], 0.0)
    cats = ["job", "marital", "education", "default", "housing", "loan", "contact", "poutcome"]
    if macro:
        for c in ["emp_var_rate", "cons_price_idx", "cons_conf_idx", "euribor3m", "nr_employed"]:
            num[c] = d[c].astype(float)
        cats += ["month", "day_of_week"]
    X = pd.concat([num, pd.get_dummies(d[cats], drop_first=True, dtype=float)], axis=1)
    X = X.loc[:, X.std(ddof=0) > 0]
    X = (X - X.mean()) / X.std(ddof=0)
    C, keep, cols = np.abs(np.corrcoef(X.values, rowvar=False)), [], list(X.columns)
    for j in range(len(cols)):
        if all(C[j, i] < 0.999 for i in keep):
            keep.append(j)
    return X.values[:, keep], [cols[j] for j in keep]


X_all, XCOLS = build_X(df)
print(f"design: {X_all.shape[1]} predetermined covariates; duration excluded")


class Spec:
    """Collapsed person-period exposure rows; dmap allows a constrained delta."""

    def __init__(self, K, dv, X, dmap=None):
        o = np.argsort(K, kind="stable")
        self.K, self.dv, self.X = K[o], dv[o], np.ascontiguousarray(X[o])
        self.N, self.P = len(K), X.shape[1]
        self.dmap = np.arange(KMAX) if dmap is None else np.asarray(dmap)
        self.nd = int(self.dmap.max()) + 1
        self.bj = self.dmap[np.minimum(self.K, KMAX) - 1]
        st = np.searchsorted(self.K, np.arange(1, KMAX + 1), side="left")
        cid, jj, ww = [], [], []
        for j in range(1, KMAX + 1):
            s = st[j - 1]
            cid.append(np.arange(s, self.N))
            jj.append(np.full(self.N - s, self.dmap[j - 1]))
            ww.append(np.ones(self.N - s) if j < KMAX else (self.K[s:] - KMAX + 1).astype(float))
        self.cid, self.jj, self.ww = (np.concatenate(cid), np.concatenate(jj), np.concatenate(ww))

    def npar(self, G): return self.nd + self.P + 2 * (G - 1)

    def unpack(self, p, G):
        beta = p[self.nd:self.nd + self.P]
        th, a = np.zeros(G), np.zeros(G)
        if G > 1:
            q = self.nd + self.P
            th[1:], a[1:] = p[q:q + G - 1], p[q + G - 1:]
        return p[:self.nd][self.dmap], beta, th, np.exp(a - logsumexp(a))

    def nll(self, p, G, link="logit"):
        _, beta, th, pi = self.unpack(p, G)
        dr = p[:self.nd]
        xb = self.X @ beta
        be, bx = dr[self.bj] + xb, dr[self.jj] + xb[self.cid]
        ell, cache = np.empty((self.N, G)), []
        for g in range(G):
            ex, ev = np.clip(bx + th[g], -30, 12), np.clip(be + th[g], -30, 12)
            if link == "logit":
                Tx, Tp, Qv, Qp = np.logaddexp(0.0, ex), expit(ex), ev, np.ones_like(ev)
            else:
                Tx = np.exp(ex); Tp = Tx
                uv = np.exp(ev); S = -np.expm1(-uv)
                Qv, Qp = np.log(S) + uv, uv * np.exp(-uv) / S + uv
            ell[:, g] = self.dv * Qv - np.bincount(self.cid, weights=self.ww * Tx, minlength=self.N)
            cache.append((Tp, Qp))
        lw = np.log(pi)[None, :] + ell
        m = logsumexp(lw, axis=1)
        post = np.exp(lw - m[:, None])
        gd, gb, gth, ga = np.zeros(self.nd), np.zeros(self.P), np.zeros(G), np.zeros(G)
        for g in range(G):
            Tp, Qp = cache[g]
            wr = self.ww * Tp * post[self.cid, g]
            evr = post[:, g] * self.dv * Qp
            gd += np.bincount(self.bj, weights=evr, minlength=self.nd) \
                - np.bincount(self.jj, weights=wr, minlength=self.nd)
            resid = evr - np.bincount(self.cid, weights=wr, minlength=self.N)
            gb += self.X.T @ resid
            gth[g] = resid.sum()
            ga[g] = post[:, g].sum() - pi[g] * self.N
        gr = np.concatenate([gd, gb, gth[1:], ga[1:]]) if G > 1 else np.concatenate([gd, gb])
        return -float(m.sum()), -gr

    def fit(self, G, link="logit", start=None):
        if start is None:
            p0 = np.zeros(self.npar(G))
            p0[:self.nd] = -2.8
            if G > 1:
                p0[self.nd + self.P:self.nd + self.P + G - 1] = np.linspace(-1.5, 1.5, G - 1)
        else:
            p0 = start
        return minimize(self.nll, p0, args=(G, link), jac=True, method="L-BFGS-B",
                        options=dict(maxiter=1500, maxfun=3000, ftol=1e-12, gtol=1e-8))

    def sweep(self, Gs, link="logit"):
        out, start = {}, None
        for G in Gs:
            r = self.fit(G, link, start)
            out[G] = dict(res=r, ll=float(-r.fun), npar=self.npar(G),
                          bic=float(2 * r.fun + self.npar(G) * np.log(self.N)))
            nx = np.zeros(self.npar(G + 1))
            q, qn = self.nd + self.P, self.nd + self.P
            nx[:q] = r.x[:q]
            if G > 1:
                nx[qn:qn + G - 1] = r.x[q:q + G - 1]
                nx[qn + G:qn + 2 * G - 1] = r.x[q + G - 1:]
            nx[qn + G - 1] = 1.5 * (-1) ** G
            start = nx
        return out


def cell(k): return min(k, KMAX) - 1


def profile(spec, p, G, kmax=10, link="logit"):
    delta, beta, th, pi = spec.unpack(p, G)
    xb = spec.X @ beta
    z = xb[:, None] + th[None, :]
    logS, rows = np.zeros((spec.N, G)), []
    for k in range(1, kmax + 1):
        e = np.clip(delta[cell(k)] + z, -30, 12)
        e1 = np.clip(delta[0] + z, -30, 12)
        lw = np.log(pi)[None, :] + logS
        q = np.exp(lw - logsumexp(lw, axis=1, keepdims=True))
        H = (lambda x: expit(x)) if link == "logit" else (lambda x: -np.expm1(-np.exp(x)))
        h, h1, r = H(e), H(e1), spec.K >= k
        rows.append(dict(k=k, n_risk=int(r.sum()), delta=float(delta[cell(k)]),
                         Exb=float(xb[r].mean()),
                         Etheta=float((q[r] * th[None, :]).sum(1).mean()),
                         h_model=float((q[r] * h[r]).sum(1).mean()),
                         h_no_fatigue=float((q[r] * h1[r]).sum(1).mean()),
                         h_no_selection=float((pi[None, :] * h).sum(1).mean())))
        logS = logS + (-np.logaddexp(0.0, e) if link == "logit" else -np.exp(e))
    return rows


def decompose(rows, k):
    a, b = rows[0], rows[k - 1]
    tot = sum(b[n] - a[n] for n in ("delta", "Exb", "Etheta"))
    p = dict(fatigue=b["delta"] - a["delta"], observable=b["Exb"] - a["Exb"],
             unobservable=b["Etheta"] - a["Etheta"])
    return dict(k=k, index_change_total=float(tot),
                **{f"index_change_{n}": float(v) for n, v in p.items()},
                **{f"share_{n}": float(v / tot) for n, v in p.items()})


# --------------------------------------------------------- 1. life table ----
life = []
for k in range(1, KMAX):
    r, e = int((K_all >= k).sum()), int(((K_all == k) & (d_all == 1)).sum())
    life.append(dict(k=k, n_risk=r, exposures=r, events=e, hazard=e / r,
                     se=float(np.sqrt(e / r * (1 - e / r) / r))))
r = int((K_all >= KMAX).sum())
expo11 = int((K_all[K_all >= KMAX] - KMAX + 1).sum())
e = int(((K_all >= KMAX) & (d_all == 1)).sum())
life.append(dict(k=KMAX, n_risk=r, exposures=expo11, events=e, hazard=e / expo11,
                 se=float(np.sqrt(e / expo11 * (1 - e / expo11) / expo11)), note="pooled 11+"))
lo = lambda h: float(np.log(h / (1 - h)))
raw_lo_18 = lo(life[7]["hazard"]) - lo(life[0]["hazard"])  # data-only benchmark
print("raw hazard 1..8:", [round(x["hazard"], 4) for x in life[:8]])

# ---------------------------------------------------------- 2. estimation ---
spec = Spec(K_all, d_all, X_all)
fits = spec.sweep([1, 2, 3, 4])
print("BIC by G:", {g: round(fits[g]["bic"], 1) for g in fits}, f"t={time.time()-T0:.0f}s")
GS = min(fits, key=lambda g: fits[g]["bic"])
pstar = fits[GS]["res"].x
delta_s, beta_s, th_s, pi_s = spec.unpack(pstar, GS)
type_h1 = [float(expit(delta_s[0] + t)) for t in th_s]
print(f"G*={GS} pi={np.round(pi_s,3).tolist()} type first-contact hazards={type_h1}")

npar = spec.npar(GS)
H = np.zeros((npar, npar))
for j in range(npar):
    ep = np.zeros(npar); ep[j] = 1e-5
    H[:, j] = (spec.nll(pstar + ep, GS)[1] - spec.nll(pstar - ep, GS)[1]) / 2e-5
V = np.linalg.pinv(0.5 * (H + H.T))
se = np.sqrt(np.clip(np.diag(V), 0, None))
v = np.zeros(npar); v[cell(8)] = 1.0; v[cell(1)] = -1.0
se_d81 = float(np.sqrt(v @ V @ v))

# ------------------------------------------------------ 3. decomposition ----
rows_M2 = profile(spec, pstar, GS)
dec8, dec5 = decompose(rows_M2, 8), decompose(rows_M2, 5)
dec8_M1 = decompose(profile(spec, fits[1]["res"].x, 1), 8)
sh_lo = (dec8["index_change_fatigue"] + 1.96 * se_d81) / dec8["index_change_total"]
sh_hi = (dec8["index_change_fatigue"] - 1.96 * se_d81) / dec8["index_change_total"]
print(f"share_fatigue={dec8['share_fatigue']:.3f} ci=[{sh_lo:.3f},{sh_hi:.3f}] "
      f"delta81={dec8['index_change_fatigue']:.3f}({se_d81:.3f}) t={time.time()-T0:.0f}s")

# ---- identification check: can pure selection, with no fatigue, do the job? -
sp0 = Spec(K_all, d_all, X_all, dmap=np.zeros(KMAX, int))
f0 = sp0.sweep([2, 3, 4, 5])
G0 = min(f0, key=lambda g: f0[g]["bic"])
rows_M3 = profile(sp0, f0[G0]["res"].x, G0)
pure = dict(G_selected=int(G0), loglik=float(f0[G0]["ll"]), bic=float(f0[G0]["bic"]),
            best_loglik=float(max(f0[g]["ll"] for g in f0)),
            lr_vs_unrestricted=float(2 * (fits[GS]["ll"] - max(f0[g]["ll"] for g in f0))),
            hazard=[float(r["h_model"]) for r in rows_M3],
            h8_pure_selection=float(rows_M3[7]["h_model"]), h8_observed=life[7]["hazard"],
            bic_table={str(g): float(f0[g]["bic"]) for g in f0})
print(f"pure-selection G*={G0} ll={pure['best_loglik']:.1f} h8={pure['h8_pure_selection']:.4f} "
      f"vs raw {life[7]['hazard']:.4f} t={time.time()-T0:.0f}s")

# ----------------------------------------------------------- 4. policy ------
rng = np.random.default_rng(SEED)
perm = rng.permutation(len(df))
tr, te = perm[:len(df) // 2], perm[len(df) // 2:]
sp_tr = Spec(K_all[tr], d_all[tr], X_all[tr])
r_tr = sp_tr.fit(GS, start=pstar)
dtr, btr, ttr, ptr = sp_tr.unpack(r_tr.x, GS)


def hpath(X):
    xb = X @ btr
    z = xb[:, None] + ttr[None, :]
    logS, Hm = np.zeros_like(z), np.zeros((len(xb), KPOL))
    for k in range(1, KPOL + 1):
        e = np.clip(dtr[cell(k)] + z, -30, 12)
        lw = np.log(ptr)[None, :] + logS
        q = np.exp(lw - logsumexp(lw, axis=1, keepdims=True))
        Hm[:, k - 1] = (q * expit(e)).sum(1)
        logS = logS - np.logaddexp(0.0, e)
    return Hm, xb


Hp_tr, xb_tr = hpath(X_all[tr])
Hp_te, xb_te = hpath(X_all[te])
Ktr, Dtr, Kte, Dte = K_all[tr], d_all[tr], K_all[te], d_all[te]
ed = np.quantile(xb_tr, np.linspace(0, 1, 11)[1:-1])
gtr, gte = np.searchsorted(ed, xb_tr), np.searchsorted(ed, xb_te)
CAPS = np.arange(0, KPOL + 1)
mono8 = float((np.diff(Hp_te[:, :8], axis=1) <= 1e-12).all(1).mean())
mono20 = float((np.diff(Hp_te, axis=1) <= 1e-12).all(1).mean())
bs = rng.integers(0, len(te), size=(300, len(te)))
pvec = lambda K, D, kb, rho: D * (K <= kb) - rho * np.minimum(K, kb)
policy = []
for rho in (0.005, 0.01, 0.02, 0.03, 0.05):
    cap_u = int(CAPS[int(np.argmax([pvec(Ktr, Dtr, c, rho).mean() for c in CAPS]))])
    mt, mv = Hp_tr < rho, Hp_te < rho
    kb_te = np.where(mv.any(1), mv.argmax(1), KPOL)
    dcap = np.array([CAPS[int(np.argmax([pvec(Ktr[gtr == g], Dtr[gtr == g], c, rho).mean()
                                         for c in CAPS]))] for g in range(10)])
    tgt_te = np.where(Hp_te[:, 0] >= rho, KPOL, 0)
    tc = np.array([pvec(Ktr, Dtr, np.where(Hp_tr[:, 0] >= rho, c, 0), rho).mean() for c in CAPS])
    cap_t = int(CAPS[int(np.argmax(tc))])
    pv = dict(observed=pvec(Kte, Dte, Kte, rho), uniform=pvec(Kte, Dte, cap_u, rho),
              myopic=pvec(Kte, Dte, kb_te, rho), decile=pvec(Kte, Dte, dcap[gte], rho),
              targeting=pvec(Kte, Dte, tgt_te, rho),
              target_cap=pvec(Kte, Dte, np.where(Hp_te[:, 0] >= rho, cap_t, 0), rho))
    po = float(pv["observed"].mean())
    gu = float(pv["uniform"].mean() - pv["observed"].mean())
    gm = float(pv["myopic"].mean() - pv["observed"].mean())
    dmu = (pv["myopic"] - pv["uniform"])[bs].mean(1)
    duo = (pv["uniform"] - pv["observed"])[bs].mean(1)
    policy.append(dict(
        rho=rho, best_uniform_cap=cap_u, decile_caps=dcap.tolist(), targeting_cap=cap_t,
        mean_myopic_cap=float(kb_te.mean()), share_never_called=float((kb_te == 0).mean()),
        **{f"profit_{k}": float(v.mean()) for k, v in pv.items()},
        calls_observed=float(Kte.mean()), calls_uniform=float(np.minimum(Kte, cap_u).mean()),
        calls_myopic=float(np.minimum(Kte, kb_te).mean()),
        conv_observed=float(Dte.mean()), conv_uniform=float((Dte * (Kte <= cap_u)).mean()),
        conv_myopic=float((Dte * (Kte <= kb_te)).mean()),
        gain_uniform=gu, gain_myopic=gm,
        gain_myopic_pct_of_observed=(100.0 * gm / po) if po > 1e-6 else None,
        gain_uniform_pct_of_observed=(100.0 * gu / po) if po > 1e-6 else None,
        myopic_over_uniform=(float(gm / gu) - 1.0) if abs(gu) > 1e-9 else None,
        ci_myopic_minus_uniform=[float(np.percentile(dmu, 2.5)), float(np.percentile(dmu, 97.5))],
        ci_uniform_minus_observed=[float(np.percentile(duo, 2.5)), float(np.percentile(duo, 97.5))]))
    print(f"rho={rho} cap*={cap_u} obs={pv['observed'].mean():.4f} uni={pv['uniform'].mean():.4f} "
          f"myo={pv['myopic'].mean():.4f} t={time.time()-T0:.0f}s")
main = [p for p in policy if p["rho"] == 0.02][0]
h2 = (main["myopic_over_uniform"] is not None) and (main["myopic_over_uniform"] >= 0.10)

# ------------------------------------------------------- 5. robustness ------
rob = {}
Xm, _ = build_X(df, macro=True)
sm = Spec(K_all, d_all, Xm)
rm = sm.fit(GS)
rob["with_calendar_macro"] = dict(decompose(profile(sm, rm.x, GS), 8), n_cov=int(Xm.shape[1]),
                                  loglik=float(-rm.fun))
rc = spec.fit(GS, link="cloglog")
rob["cloglog"] = dict(decompose(profile(spec, rc.x, GS, link="cloglog"), 8), loglik=float(-rc.fun))
sel = K_all <= 12
ss = Spec(K_all[sel], d_all[sel], X_all[sel])
rs = ss.fit(GS)
rob["campaign_le_12"] = dict(decompose(profile(ss, rs.x, GS), 8), n_spells=int(sel.sum()))
for G in fits:
    if G != GS and G > 1:
        rob[f"G_{G}"] = decompose(profile(spec, fits[G]["res"].x, G), 8)

h1 = bool(dec8["share_fatigue"] < 0.5)
out = dict(
    n_spells=int(len(df)), n_calls=int(K_all.sum()), n_conversions=int(d_all.sum()),
    n_person_periods=int(K_all.sum()), n_exposure_rows=int(len(spec.cid)),
    conversion_rate=float(d_all.mean()), mean_calls_per_spell=float(K_all.mean()),
    max_contacts=int(K_all.max()), support_gaps=gaps, n_covariates=int(X_all.shape[1]),
    life_table=life, raw_logodds_change_1_to_8=raw_lo_18,
    raw_hazard_ratio_8_to_1=float(life[7]["hazard"] / life[0]["hazard"]),
    model_selection={str(g): dict(loglik=fits[g]["ll"], npar=fits[g]["npar"], bic=fits[g]["bic"])
                     for g in fits},
    G_selected=int(GS), lr_frailty=float(2 * (fits[2]["ll"] - fits[1]["ll"])),
    theta=[float(x) for x in th_s], pi=[float(x) for x in pi_s], type_first_hazard=type_h1,
    delta=[float(x) for x in delta_s], delta_se=[float(x) for x in se[:KMAX]],
    delta_train_half=[float(x) for x in dtr],
    delta_8_minus_1=float(dec8["index_change_fatigue"]), delta_8_minus_1_se=se_d81,
    hazard_profile=rows_M2, decomposition_k8=dec8, decomposition_k5=dec5,
    decomposition_k8_observables_only=dec8_M1,
    share_fatigue_ci95_numerator=[sh_lo, sh_hi],
    pure_selection_model=pure, myopic_path_monotone_share_k1_8=mono8, myopic_path_monotone_share_k1_20=mono20, policy=policy, policy_main_rho=0.02, robustness=rob,
    hypothesis_H1_selection_dominates=h1, hypothesis_H2_personalisation_gain=bool(h2),
    hypothesis_supported=h1, runtime_seconds=float(time.time() - T0))
print("RESULT_JSON: " + json.dumps(out))

Results and analysis

Key numbers

Caveats

• The data are spell-level, not call-level. The file records the total number of contacts (`campaign`) and the terminal outcome (`y`), not the date, order or outcome of each individual call. Contact order is therefore reconstructed from spell length under the maintained assumption that a spell ends at the conversion if one occurs and is right-censored otherwise. If a client could convert and still be called again, or if `campaign` were recorded with error, the life table and every model built on it would shift.
• Censoring is chosen by the bank, not by nature. The bank decided when to stop calling, and if it stopped on information correlated with the unobserved type theta, the frailty estimate is contaminated and the selection share is biased. This is a maintained assumption for the decomposition. It is NOT needed for the policy result, which evaluates only truncations of the observed policy and is exactly identified whatever generated the censoring.
• The frailty component is weakly identified. BIC separates G=2 from G=1 by only 0.33 points, the second mass point is at theta = -14.65, effectively a degenerate never-converter type, and the likelihood-ratio test is non-standard at the boundary. The 10.3 percent unobserved-selection share should be read as indicative, not precise. The observables-only decomposition, which needs no frailty assumption, is reported alongside it and gives an even smaller selection share.
• The 95 percent interval on the fatigue share, [0.498, 1.026], has a lower end fractionally below the 0.50 refutation threshold, and it accounts only for sampling error in delta_8 - delta_1, holding the denominator fixed. A 50/50 split cannot be rejected at conventional levels even though the point estimate is 0.762.
• Contact 8 rests on 17 conversions among 1,777 clients at risk, and contacts 9 and 10 on 17 and 12. The tail of the contact-order profile is estimated imprecisely, its non-monotonicity beyond contact 8 may be noise, and this is the reason the myopic stopping rule is only a heuristic.
• Identification of delta against theta relies on the logistic (and, in robustness, complementary log-log) functional form together with variation in the 32 predetermined covariates. Those covariates are weak predictors in this file and the mean spell is only 2.568 calls, so the split between the two rests more heavily on functional form than one would like.
• Calendar month and the five macro indicators are stamped at the last contact and are therefore not predetermined; including them mechanically absorbs part of the contact-order effect. The specification that includes them is reported as an upper bound on calendar composition, not as a preferred model, and it still leaves fatigue at 64.3 percent.
• Profit is normalised by the conversion margin and no true cost per call or margin per term deposit is observed; the cost-to-margin ratio rho is swept over an assumed grid of 0.005 to 0.05 rather than calibrated to this bank's economics. Which conclusion about H2 holds depends on where the bank actually sits on that grid.
• Only truncations of the observed policy are evaluated. Policies that call a client MORE than the bank did, or that reorder or retime calls, are outside the evidence entirely, and nothing here says what would happen under them.
• This is one campaign at one Portuguese retail bank between May 2008 and November 2010, a period spanning the financial crisis, selling one product by telephone. Contact fatigue in this setting need not generalise to other channels, products, countries or periods, and the macro environment moved substantially within the sample.
• `duration` is excluded from every model because it is realised only after a call ends. This is the right choice for a targeting or stopping rule but it means the models deliberately ignore the single strongest predictor in the file, and the fitted hazards are correspondingly modest in discriminating power.
• Three spells with `default` = yes and 18 with `education` = illiterate were recoded to 'unknown' to avoid perfect separation in the likelihood. This affects 21 of 41,188 spells and no observation was dropped, but it is a deliberate coding choice rather than a property of the data.

≈ 3398 words · AI-drafted manuscript · not peer-reviewed

The revision addresses every point I raised. The overstated proposition has been replaced by a correct one, the manuscript now says plainly why no non-vacuous selection bound exists and tests pure selection by likelihood and BIC instead, 'lowest attainable' is gone, the abstract reports the personalisation result at every cost ratio including the two where it reverses, Assumption 1 is given the sceptical treatment it needs, and the seed, optimiser settings, Hessian step and pseudo-inverse are all reported. I re-checked the manuscript against the executed record and every number matches; the dataset, sample size, filters and provenance are exactly those recorded, and nothing is reported that the script did not run. What remains is small: a cross-reference that claims more for Proposition 2 than it proves, and a life table whose omitted rows include the one that breaks monotonicity. The formal apparatus is still light for Management Science and the frailty split is weakly identified, but the paper says so itself, and the central empirical result - that roughly three quarters of the decline in telemarketing response across contacts is genuine fatigue rather than dynamic selection, with the policy corollary that firms should cap uniformly and personalise the list - is well identified, well evidenced and useful. I recommend minor revision.