diff --git a/METHODOLOGY_REVIEW.md b/METHODOLOGY_REVIEW.md index 1cf8669a..a6e21e6d 100644 --- a/METHODOLOGY_REVIEW.md +++ b/METHODOLOGY_REVIEW.md @@ -30,7 +30,7 @@ Each estimator in diff-diff should be periodically reviewed to ensure: | StackedDiD | `stacked_did.py` | `stacked-did-weights` | **Complete** | 2026-02-19 | | TROP | `trop.py` | (forthcoming) | Not Started | - | | BaconDecomposition | `bacon.py` | `bacondecomp::bacon()` | Not Started | - | -| HonestDiD | `honest_did.py` | `HonestDiD` package | Not Started | - | +| HonestDiD | `honest_did.py` | `HonestDiD` package | **Complete** | 2026-03-31 | | PreTrendsPower | `pretrends.py` | `pretrends` package | Not Started | - | | PowerAnalysis | `power.py` | `pwr` / `DeclareDesign` | Not Started | - | @@ -618,14 +618,84 @@ variables appear to the left of the `|` separator. | Module | `honest_did.py` | | Primary Reference | Rambachan & Roth (2023) | | R Reference | `HonestDiD` package | -| Status | Not Started | -| Last Review | - | +| Status | **Complete** (pending R comparison) | +| Last Review | 2026-04-01 | + +**Verified Components:** +- [x] Delta^SD: second-difference constraints [1,-2,1] with delta_0=0 boundary handling +- [x] Delta^SD: T+Tbar-1 constraint rows (bridge constraint at t=0) +- [x] Delta^RM: constrains first differences (not levels), union of polyhedra per Lemma 2.2 +- [x] Identified set LP: pins delta_pre = beta_pre via equality constraints (Equations 5-6) +- [x] M=0 for Delta^SD: linear extrapolation gives finite point-identified bounds +- [x] Mbar=0 for Delta^RM: point identification (all post first-diffs = 0) +- [x] Optimal FLCI for Delta^SD: folded normal cv_alpha, Nelder-Mead over pre-period weights +- [x] Sensitivity grid: bounds computed for each M in grid, breakdown value via binary search +- [x] Survey variance (RM, M=0 smoothness): t-distribution critical values from df_survey +- [ ] Survey variance (M>0 smoothness): optimal FLCI uses asymptotic normal only; df_survey=0 → NaN +- [x] CallawaySantAnna integration: universal base period, reference period filtering +- [x] Three-period analytical case matches paper Section 2.3 +- [ ] ARP hybrid for Delta^RM: infrastructure implemented, moment inequality transformation needs calibration +- [ ] R comparison: pending (benchmark scripts need updating) + +**Test Coverage:** +- 63 existing tests in `tests/test_honest_did.py` (14 classes) — all passing +- 17 new methodology verification tests in `tests/test_methodology_honest_did.py` +- R benchmark tests (pending) **Corrections Made:** -- (None yet) +1. **DeltaRM: first differences, not levels** (`honest_did.py`, `_construct_constraints_rm_component`): + The paper's Delta^RM constrains `|delta_{t+1} - delta_t|` (consecutive first differences) + bounded by Mbar × max pre-treatment first difference. The code constrained `|delta_post|` + (absolute levels) bounded by Mbar × max `|beta_pre|`. Completely rewritten using + union-of-polyhedra decomposition per Lemma 2.2. + +2. **LP pins delta_pre = beta_pre** (`honest_did.py`, `_solve_bounds_lp`): + The paper's identified set LP (Equations 5-6) fixes pre-treatment violations to the observed + pre-treatment coefficients. The code had no equality constraint — delta_pre was unconstrained. + For Delta^SD(M=0), this made the LP unbounded. Added A_eq/b_eq equality constraints. + +3. **DeltaSD constraint matrix: delta_0=0 boundary** (`honest_did.py`, `_construct_A_sd`): + The code built second-difference matrices treating [delta_{-T},...,delta_{-1},delta_1,...,delta_{Tbar}] + as consecutive, missing delta_0=0 at the boundary. Three boundary rows were wrong: + - t=-1: `d_{-2} - 2*d_{-1} + 0` (uses delta_0=0) + - t=0: `d_{-1} + d_1` (bridge constraint, was missing) + - t=1: `0 - 2*d_1 + d_2` (uses delta_0=0) + Now produces T+Tbar-1 rows (was T+Tbar-2). + +4. **Optimal FLCI for Delta^SD** (`honest_did.py`, `_compute_optimal_flci`): + Replaced naive FLCI (`lb - z*se, ub + z*se`) with the paper's optimal FLCI (Section 4.1): + jointly optimizes affine estimator direction v and half-length chi using folded normal + critical values cv_alpha(bias/se). Significantly narrower CIs. + +5. **REGISTRY.md equations** (`docs/methodology/REGISTRY.md`): + DeltaSD equation was first differences (should be second differences). DeltaRM equation + was absolute levels (should be first differences). Both corrected with full formulations. + +6. **Performance** (`honest_did.py`): + Sensitivity grid reduced from ~9 minutes to 0.1 seconds via: Newton's method for cv_alpha + (5 iterations vs 100), centrosymmetric bias LP (1 solve vs 2), M=0 short-circuit, + looser Nelder-Mead tolerances. **Outstanding Concerns:** -- (None yet) +- **Delta^RM CI**: uses naive FLCI (conservative) instead of the paper's ARP conditional/hybrid + confidence sets. ARP infrastructure exists but moment inequality transformation needs + calibration. Tracked in TODO.md. +- R benchmark comparison not yet run (Python benchmark needs API update) +- Combined method uses single M for both SD and RM (DeltaSDRM dataclass has separate M/Mbar) + +**Deviations from R's HonestDiD:** +1. **Deviation from R:** Delta^RM CIs use naive FLCI (`lb - z*se, ub + z*se`) instead of ARP + conditional/hybrid. Conservative (wider CIs, valid coverage). ARP deferred. +2. **Note:** Delta^SD optimal FLCI matches the paper's Section 4.1 methodology: first-difference + reparameterization, slope weights with sum(w)=sum_j j*l_j constraint (Eq. 17), bias LP in + fd-space, folded normal (or folded non-central t for survey df). Nelder-Mead optimizer vs + R's custom solver may produce numerical differences at tolerance level. +3. **Note:** `method="combined"` (Delta^SDRM) uses naive FLCI on the intersection of SD and RM + bounds. The paper proves FLCI is not consistent for Delta^SDRM (Proposition 4.2). A runtime + UserWarning is emitted. Use `method="smoothness"` or `method="relative_magnitude"` separately + for paper-supported inference. +4. **Note (deviation from R):** Python warns (doesn't error) when CallawaySantAnna results use + `base_period != "universal"`. R's HonestDiD requires universal base period. --- diff --git a/TODO.md b/TODO.md index 7da9a30b..aafe0f4e 100644 --- a/TODO.md +++ b/TODO.md @@ -68,6 +68,7 @@ Deferred items from PR reviews that were not addressed before merge. | StaggeredTripleDifference R cross-validation: CSV fixtures not committed (gitignored); tests skip without local R + triplediff. Commit fixtures or generate deterministically. | `tests/test_methodology_staggered_triple_diff.py` | #245 | Medium | | StaggeredTripleDifference R parity: benchmark only tests no-covariate path (xformla=~1). Add covariate-adjusted scenarios and aggregation SE parity assertions. | `benchmarks/R/benchmark_staggered_triplediff.R` | #245 | Medium | | StaggeredTripleDifference: per-cohort group-effect SEs include WIF (conservative vs R's wif=NULL). Documented in REGISTRY. Could override mixin for exact R match. | `staggered_triple_diff.py` | #245 | Low | +| HonestDiD Delta^RM: uses naive FLCI instead of paper's ARP conditional/hybrid confidence sets (Sections 3.2.1-3.2.2). ARP infrastructure exists but moment inequality transformation needs calibration. CIs are conservative (wider, valid coverage). | `honest_did.py` | #248 | Medium | #### Performance diff --git a/diff_diff/honest_did.py b/diff_diff/honest_did.py index ae9d660e..02f84a83 100644 --- a/diff_diff/honest_did.py +++ b/diff_diff/honest_did.py @@ -73,18 +73,18 @@ class DeltaRM: """ Relative magnitudes restriction on trend violations (Delta^{RM}). - Post-treatment violations are bounded by Mbar times the maximum - absolute pre-treatment violation: - |delta_post| <= Mbar * max(|delta_pre|) + Post-treatment consecutive first differences are bounded by Mbar + times the maximum pre-treatment first difference: + |delta_{t+1} - delta_t| <= Mbar * max_{s<0} |delta_{s+1} - delta_s| - When Mbar=0, this enforces exact parallel trends post-treatment. - Mbar=1 means post-period violations can be as large as the worst - observed pre-period violation. + When Mbar=0, this enforces zero post-treatment first differences. + Mbar=1 means post-period first differences can be as large as the + worst observed pre-period first difference. Parameters ---------- Mbar : float - Scaling factor for maximum pre-period violation. + Scaling factor for maximum pre-period first difference. Examples -------- @@ -110,7 +110,7 @@ class DeltaSDRM: Imposes both: 1. Smoothness: |delta_{t+1} - 2*delta_t + delta_{t-1}| <= M - 2. Relative magnitudes: |delta_post| <= Mbar * max(|delta_pre|) + 2. Relative magnitudes: |delta_{t+1} - delta_t| <= Mbar * max_{s<0} |delta_{s+1} - delta_s| This is more restrictive than either constraint alone. @@ -119,7 +119,7 @@ class DeltaSDRM: M : float Maximum allowed second difference (smoothness). Mbar : float - Scaling factor for maximum pre-period violation (relative magnitudes). + Scaling factor for maximum pre-period first difference (relative magnitudes). Examples -------- @@ -198,7 +198,17 @@ class HonestDiDResults: survey_metadata: Optional[Any] = field(default=None, repr=False) df_survey: Optional[int] = field(default=None, repr=False) + def _ci_is_finite(self) -> bool: + """Check if CI endpoints are finite (not NaN/inf).""" + return np.isfinite(self.ci_lb) and np.isfinite(self.ci_ub) + def __repr__(self) -> str: + if not self._ci_is_finite(): + return ( + f"HonestDiDResults(bounds=[{self.lb}, {self.ub}], " + f"CI=[{self.ci_lb}, {self.ci_ub}] (undefined), " + f"M={self.M})" + ) sig = "" if self.ci_lb <= 0 <= self.ci_ub else "*" return ( f"HonestDiDResults(bounds=[{self.lb:.4f}, {self.ub:.4f}], " @@ -208,7 +218,12 @@ def __repr__(self) -> str: @property def is_significant(self) -> bool: - """Check if CI excludes zero (effect is robust to violations).""" + """Check if CI excludes zero (effect is robust to violations). + + Returns False for undefined (NaN) CIs. + """ + if not self._ci_is_finite(): + return False return not (self.ci_lb <= 0 <= self.ci_ub) @property @@ -288,7 +303,7 @@ def summary(self) -> str: if self.method == "relative_magnitude": lines.append( - f"Post-treatment violations bounded at {self.M:.1f}x max pre-period violation." + f"Post-treatment first differences bounded at {self.M:.1f}x max pre-period first difference." ) elif self.method == "smoothness": if self.M == 0: @@ -441,7 +456,7 @@ def to_dataframe(self) -> pd.DataFrame: "ub": ub, "ci_lb": ci_lb, "ci_ub": ci_ub, - "is_significant": not (ci_lb <= 0 <= ci_ub), + "is_significant": (np.isfinite(ci_lb) and np.isfinite(ci_ub) and not (ci_lb <= 0 <= ci_ub)), } ) return pd.DataFrame(rows) @@ -808,45 +823,100 @@ def _extract_event_study_params( ) -def _construct_A_sd(num_periods: int) -> np.ndarray: +def _construct_A_sd(num_pre_periods: int, num_post_periods: int) -> np.ndarray: """ Construct constraint matrix for smoothness (second differences). - For T periods, creates matrix A such that: - A @ delta gives the second differences. + Builds the matrix A such that A @ delta gives the second differences, + accounting for the normalization delta_0 = 0 at the pre-post boundary. + + The delta vector is [delta_{-T}, ..., delta_{-1}, delta_1, ..., delta_{Tbar}] + (delta_0 = 0 is omitted). Second differences at the boundary use delta_0 = 0: + t=-1: delta_{-2} - 2*delta_{-1} + 0 (if num_pre >= 2) + t= 0: delta_{-1} + delta_1 (bridge constraint, always present) + t= 1: 0 - 2*delta_1 + delta_2 (if num_post >= 2) Parameters ---------- - num_periods : int - Number of time periods. + num_pre_periods : int + Number of pre-treatment periods (T). + num_post_periods : int + Number of post-treatment periods (Tbar). Returns ------- A : np.ndarray - Constraint matrix of shape (num_periods - 2, num_periods). + Constraint matrix of shape (n_constraints, num_pre + num_post). + n_constraints = num_pre + num_post - 1 for sufficient periods, + accounting for the delta_0 = 0 boundary. """ - if num_periods < 3: - return np.zeros((0, num_periods)) - - n_constraints = num_periods - 2 - A = np.zeros((n_constraints, num_periods)) - - for i in range(n_constraints): - # Second difference: delta_{t+1} - 2*delta_t + delta_{t-1} - A[i, i] = 1 # delta_{t-1} - A[i, i + 1] = -2 # delta_t - A[i, i + 2] = 1 # delta_{t+1} - - return A + T = num_pre_periods + Tbar = num_post_periods + total = T + Tbar + + if total < 2: + return np.zeros((0, total)) + + rows = [] + + # Pure pre-period second differences: t = -T+1, ..., -2 + # These involve delta[i-1], delta[i], delta[i+1] all in the pre-period block + # Row i corresponds to: delta_{-(T-i)} - 2*delta_{-(T-i-1)} + delta_{-(T-i-2)} + for i in range(T - 2): + row = np.zeros(total) + row[i] = 1 # delta_{t-1} + row[i + 1] = -2 # delta_t + row[i + 2] = 1 # delta_{t+1} + rows.append(row) + + # Boundary constraint at t = -1: delta_{-2} - 2*delta_{-1} + delta_0 + # With delta_0 = 0: delta_{-2} - 2*delta_{-1} + if T >= 2: + row = np.zeros(total) + row[T - 2] = 1 # delta_{-2} + row[T - 1] = -2 # delta_{-1} + # delta_0 = 0, no entry needed + rows.append(row) + + # Bridge constraint at t = 0: delta_{-1} - 2*delta_0 + delta_1 + # With delta_0 = 0: delta_{-1} + delta_1 + if T >= 1 and Tbar >= 1: + row = np.zeros(total) + row[T - 1] = 1 # delta_{-1} + row[T] = 1 # delta_1 + rows.append(row) + + # Boundary constraint at t = 1: delta_0 - 2*delta_1 + delta_2 + # With delta_0 = 0: -2*delta_1 + delta_2 + if Tbar >= 2: + row = np.zeros(total) + row[T] = -2 # delta_1 + row[T + 1] = 1 # delta_2 + rows.append(row) + + # Pure post-period second differences: event times t = 2, ..., Tbar-1 + # delta_{t+1} - 2*delta_t + delta_{t-1}, all within the post-period block + for t in range(2, Tbar): + row = np.zeros(total) + row[T + t - 2] = 1 # delta_{t-1} + row[T + t - 1] = -2 # delta_t + row[T + t] = 1 # delta_{t+1} + rows.append(row) + + if not rows: + return np.zeros((0, total)) + + return np.array(rows) def _construct_constraints_sd( num_pre_periods: int, num_post_periods: int, M: float ) -> Tuple[np.ndarray, np.ndarray]: """ - Construct smoothness constraint matrices. + Construct smoothness constraint matrices for Delta^SD(M). - Returns A, b such that delta in DeltaSD iff |A @ delta| <= b. + Returns A, b such that delta in DeltaSD(M) iff |A @ delta| <= b. + Accounts for delta_0 = 0 normalization at the pre-post boundary. Parameters ---------- @@ -855,7 +925,7 @@ def _construct_constraints_sd( num_post_periods : int Number of post-treatment periods. M : float - Smoothness parameter. + Smoothness parameter (max second difference). Returns ------- @@ -864,11 +934,11 @@ def _construct_constraints_sd( b_ineq : np.ndarray Inequality constraint vector. """ - total_periods = num_pre_periods + num_post_periods - A_base = _construct_A_sd(total_periods) + A_base = _construct_A_sd(num_pre_periods, num_post_periods) if A_base.shape[0] == 0: - return np.zeros((0, total_periods)), np.zeros(0) + total = num_pre_periods + num_post_periods + return np.zeros((0, total)), np.zeros(0) # |A @ delta| <= M becomes: # A @ delta <= M and -A @ delta <= M @@ -878,11 +948,21 @@ def _construct_constraints_sd( return A_ineq, b_ineq -def _construct_constraints_rm( - num_pre_periods: int, num_post_periods: int, Mbar: float, max_pre_violation: float +def _construct_constraints_rm_component( + num_pre_periods: int, + num_post_periods: int, + Mbar: float, + max_pre_first_diff: float, ) -> Tuple[np.ndarray, np.ndarray]: """ - Construct relative magnitudes constraint matrices. + Construct constraint matrices for one component of Delta^RM. + + Delta^RM constrains post-treatment FIRST DIFFERENCES (not levels): + |delta_{t+1} - delta_t| <= Mbar * max_pre_first_diff, for all t >= 0 + + With delta_0 = 0 normalization: + |delta_1| <= bound (t=0) + |delta_{t+1} - delta_t| <= bound (t=1, ..., Tbar-1) Parameters ---------- @@ -892,8 +972,8 @@ def _construct_constraints_rm( Number of post-treatment periods. Mbar : float Relative magnitude scaling factor. - max_pre_violation : float - Maximum absolute pre-period violation (estimated from data). + max_pre_first_diff : float + The pre-period first difference for this union component. Returns ------- @@ -902,26 +982,140 @@ def _construct_constraints_rm( b_ineq : np.ndarray Inequality constraint vector. """ - total_periods = num_pre_periods + num_post_periods + T = num_pre_periods + Tbar = num_post_periods + total = T + Tbar + bound = Mbar * max_pre_first_diff + + rows = [] + + # t=0: |delta_1 - delta_0| = |delta_1| <= bound (delta_0 = 0) + if Tbar >= 1: + row_pos = np.zeros(total) + row_pos[T] = 1 # delta_1 <= bound + rows.append(row_pos) + row_neg = np.zeros(total) + row_neg[T] = -1 # -delta_1 <= bound + rows.append(row_neg) + + # t=1, ..., Tbar-1: |delta_{t+1} - delta_t| <= bound + for t in range(1, Tbar): + row_pos = np.zeros(total) + row_pos[T + t] = 1 # delta_{t+1} + row_pos[T + t - 1] = -1 # -delta_t + rows.append(row_pos) + row_neg = np.zeros(total) + row_neg[T + t] = -1 # -delta_{t+1} + row_neg[T + t - 1] = 1 # delta_t + rows.append(row_neg) + + if not rows: + return np.zeros((0, total)), np.zeros(0) + + A_ineq = np.array(rows) + b_ineq = np.full(len(rows), bound) + return A_ineq, b_ineq - # Bound post-period violations: |delta_post| <= Mbar * max_pre_violation - bound = Mbar * max_pre_violation - # Create constraints for each post-period - # delta_post[i] <= bound and -delta_post[i] <= bound - n_constraints = 2 * num_post_periods - A_ineq = np.zeros((n_constraints, total_periods)) - b_ineq = np.full(n_constraints, bound) +def _compute_pre_first_differences(beta_pre: np.ndarray) -> np.ndarray: + """ + Compute pre-period first differences for Delta^RM. - for i in range(num_post_periods): - post_idx = num_pre_periods + i - A_ineq[2 * i, post_idx] = 1 # delta <= bound - A_ineq[2 * i + 1, post_idx] = -1 # -delta <= bound + With delta_0 = 0 normalization, the pre-period first differences are: + fd_s = delta_{s+1} - delta_s for s = -T, ..., -1 - return A_ineq, b_ineq + Since delta_pre = beta_pre (by no-anticipation): + fd_{-T} = beta_{-T+1} - beta_{-T} + ... + fd_{-2} = beta_{-1} - beta_{-2} + fd_{-1} = delta_0 - beta_{-1} = -beta_{-1} (boundary through delta_0=0) + + Parameters + ---------- + beta_pre : np.ndarray + Pre-period coefficient estimates [beta_{-T}, ..., beta_{-1}]. + + Returns + ------- + first_diffs : np.ndarray + Absolute first differences |fd_{-T}|, ..., |fd_{-1}|. + """ + if len(beta_pre) == 0: + return np.array([]) + + diffs = [] + # Interior first differences: fd_s = beta_{s+1} - beta_s + for i in range(len(beta_pre) - 1): + diffs.append(abs(beta_pre[i + 1] - beta_pre[i])) + # Boundary: fd_{-1} = delta_0 - delta_{-1} = 0 - beta_{-1} = -beta_{-1} + diffs.append(abs(beta_pre[-1])) + + return np.array(diffs) + + +def _solve_rm_bounds_union( + beta_pre: np.ndarray, + beta_post: np.ndarray, + l_vec: np.ndarray, + num_pre_periods: int, + Mbar: float, + lp_method: str = "highs", +) -> Tuple[float, float]: + """ + Solve identified set bounds for Delta^RM via union of polyhedra. + + Delta^RM is a union of polyhedra (one per location of the max pre-period + first difference). Per Lemma 2.2 of Rambachan & Roth (2023), the + identified set is the union of component identified sets. + + With delta_pre = beta_pre pinned, each pre-period first difference is + a known scalar, so each component LP has simple box constraints on + post-treatment first differences. + + Parameters + ---------- + beta_pre : np.ndarray + Pre-period coefficients. + beta_post : np.ndarray + Post-period coefficients. + l_vec : np.ndarray + Weighting vector. + num_pre_periods : int + Number of pre-periods. + Mbar : float + Relative magnitudes scaling factor. + lp_method : str + LP solver method. + + Returns + ------- + lb : float + Lower bound (min over all components). + ub : float + Upper bound (max over all components). + """ + pre_diffs = _compute_pre_first_differences(beta_pre) + num_post = len(beta_post) + + if len(pre_diffs) == 0 or np.max(pre_diffs) == 0: + theta = np.dot(l_vec, beta_post) + return theta, theta + + # After pinning delta_pre = beta_pre, the RM bound is determined by + # max(pre_diffs). Smaller components give tighter constraints and thus + # narrower bounds that are nested inside the max-component bounds. + # One LP call suffices (Lemma 2.2 union simplifies to max component). + max_pre_fd = float(np.max(pre_diffs)) + A_ineq, b_ineq = _construct_constraints_rm_component( + num_pre_periods, num_post, Mbar, max_pre_fd + ) + return _solve_bounds_lp( + beta_pre, beta_post, l_vec, A_ineq, b_ineq, num_pre_periods, lp_method + ) def _solve_bounds_lp( + beta_pre: np.ndarray, beta_post: np.ndarray, l_vec: np.ndarray, A_ineq: np.ndarray, @@ -932,15 +1126,19 @@ def _solve_bounds_lp( """ Solve for identified set bounds using linear programming. - The parameter of interest is theta = l' @ (beta_post - delta_post). - We find min and max over delta in the constraint set. + Computes the bounds of the identified set S(beta, Delta) per + Rambachan & Roth (2023) Equations 5-6: - Note: The optimization is over delta for ALL periods (pre + post), but - only the post-period components contribute to the objective function. - This correctly handles smoothness constraints that link pre and post periods. + theta^lb = l'beta_post - max{ l'delta_post : delta in Delta, delta_pre = beta_pre } + theta^ub = l'beta_post - min{ l'delta_post : delta in Delta, delta_pre = beta_pre } + + The equality constraint delta_pre = beta_pre pins the pre-treatment violations + to the observed pre-treatment coefficients (since tau_pre = 0 by no-anticipation). Parameters ---------- + beta_pre : np.ndarray + Pre-period coefficient estimates (pinned as equality constraints). beta_post : np.ndarray Post-period coefficient estimates. l_vec : np.ndarray @@ -958,54 +1156,58 @@ def _solve_bounds_lp( Returns ------- lb : float - Lower bound. + Lower bound of identified set. ub : float - Upper bound. + Upper bound of identified set. """ num_post = len(beta_post) total_periods = A_ineq.shape[1] if A_ineq.shape[0] > 0 else num_pre_periods + num_post - # theta = l' @ beta_post - l' @ delta_post - # We optimize over delta (all periods including pre for smoothness constraints) - - # Extract post-period part of constraints - # For delta in R^total_periods, we want min/max of -l' @ delta_post - # where delta_post = delta[num_pre_periods:] - + # Objective: min/max -l' @ delta_post over delta in R^total_periods c = np.zeros(total_periods) - c[num_pre_periods : num_pre_periods + num_post] = -l_vec # min -l'@delta = max l'@delta + c[num_pre_periods : num_pre_periods + num_post] = -l_vec - # For upper bound: max l'@(beta - delta) = l'@beta + max(-l'@delta) - # For lower bound: min l'@(beta - delta) = l'@beta + min(-l'@delta) + # Equality constraints: delta_pre = beta_pre (Rambachan & Roth Eqs 5-6) + A_eq = np.zeros((num_pre_periods, total_periods)) + for i in range(num_pre_periods): + A_eq[i, i] = 1.0 + b_eq = beta_pre - if A_ineq.shape[0] == 0: - # No constraints - unbounded + if A_ineq.shape[0] == 0 and num_pre_periods == 0: return -np.inf, np.inf - # Solve for lower bound of -l'@delta (which gives upper bound of theta) + lp_kwargs = dict( + A_ub=A_ineq if A_ineq.shape[0] > 0 else None, + b_ub=b_ineq if A_ineq.shape[0] > 0 else None, + A_eq=A_eq, + b_eq=b_eq, + bounds=(None, None), + method=lp_method, + ) + + # Solve for min(-l'@delta_post) → gives upper bound of theta try: - result_min = optimize.linprog( - c, A_ub=A_ineq, b_ub=b_ineq, bounds=(None, None), method=lp_method - ) + result_min = optimize.linprog(c, **lp_kwargs) if result_min.success: min_val = result_min.fun + elif result_min.status == 2: + # Infeasible: beta_pre inconsistent with Delta at this M + return np.nan, np.nan else: min_val = -np.inf except (ValueError, TypeError): - # Optimization failed - return unbounded min_val = -np.inf - # Solve for upper bound of -l'@delta (which gives lower bound of theta) + # Solve for max(-l'@delta_post) → gives lower bound of theta try: - result_max = optimize.linprog( - -c, A_ub=A_ineq, b_ub=b_ineq, bounds=(None, None), method=lp_method - ) + result_max = optimize.linprog(-c, **lp_kwargs) if result_max.success: max_val = -result_max.fun + elif result_max.status == 2: + return np.nan, np.nan else: max_val = np.inf except (ValueError, TypeError): - # Optimization failed - return unbounded max_val = np.inf theta_base = np.dot(l_vec, beta_post) @@ -1065,69 +1267,681 @@ def _compute_flci( return ci_lb, ci_ub -def _compute_clf_ci( +def _cv_alpha(t: float, alpha: float, df: Optional[int] = None) -> float: + """ + Compute the (1-alpha) quantile of the folded distribution |X|. + + When df is None: X ~ N(t, 1) (folded normal). + When df > 0: X ~ nct(df, t) (folded non-central t, for survey inference). + Per Rambachan & Roth (2023) Equation 18. + + Parameters + ---------- + t : float + Non-centrality parameter (bias / se ratio). + alpha : float + Significance level. + df : int, optional + Degrees of freedom for non-central t. None = normal theory. + + Returns + ------- + cv : float + Critical value such that P(|X| <= cv) = 1 - alpha. + """ + from scipy.stats import norm + + target = 1 - alpha + t = abs(t) + + if df is not None and df > 0: + # Folded non-central t: P(|nct(df,t)| <= x) = F(x;df,t) - F(-x;df,t) + from scipy.stats import nct as nct_dist + + x = nct_dist.ppf(1 - alpha / 2, df, t) + 1.0 # generous start + for _ in range(30): + f = nct_dist.cdf(x, df, t) - nct_dist.cdf(-x, df, t) - target + fprime = nct_dist.pdf(x, df, t) + nct_dist.pdf(-x, df, t) + if fprime < 1e-15: + break + x_new = x - f / fprime + x_new = max(x_new, 0.0) + if abs(x_new - x) < 1e-10: + break + x = x_new + return x + + # Folded normal: P(|N(t,1)| <= x) = Phi(x-t) - Phi(-x-t) + x = norm.ppf(1 - alpha / 2) + t + + for _ in range(20): + f = norm.cdf(x - t) - norm.cdf(-x - t) - target + fprime = norm.pdf(x - t) + norm.pdf(-x - t) + if fprime < 1e-15: + break + x_new = x - f / fprime + x_new = max(x_new, 0.0) + if abs(x_new - x) < 1e-12: + break + x = x_new + + return x + + +def _build_fd_transform(num_pre: int, num_post: int) -> np.ndarray: + """ + Build the matrix C mapping first-differences to levels: delta = C @ fd. + + The fd vector has T+Tbar components: + fd = [fd_{-T}, ..., fd_{-1}, fd_0, ..., fd_{Tbar-1}] + where fd_s = delta_{s+1} - delta_s (with delta_0 = 0). + + The delta vector is: + delta = [delta_{-T}, ..., delta_{-1}, delta_1, ..., delta_{Tbar}] + + Pre-period (backward from delta_0=0): + delta_{-1} = -fd_{T-1} + delta_{-k} = -(fd_{T-1} + fd_{T-2} + ... + fd_{T-k}) + + Post-period (forward from delta_0=0): + delta_1 = fd_T + delta_k = fd_T + fd_{T+1} + ... + fd_{T+k-1} + """ + T = num_pre + Tbar = num_post + total = T + Tbar + C = np.zeros((total, total)) + + # Pre-period: delta_{-k} = -(fd_{T-1} + fd_{T-2} + ... + fd_{T-k}) + for k in range(1, T + 1): + delta_idx = T - k # delta_{-k} is at index T-k + for j in range(k): + fd_idx = T - 1 - j # fd_{T-1-j} + C[delta_idx, fd_idx] = -1.0 + + # Post-period: delta_k = fd_T + fd_{T+1} + ... + fd_{T+k-1} + for k in range(1, Tbar + 1): + delta_idx = T + k - 1 # delta_k is at index T+k-1 + for j in range(k): + fd_idx = T + j # fd_{T+j} + C[delta_idx, fd_idx] = 1.0 + + return C + + +def _build_fd_smoothness_constraints( + num_fd: int, M: float +) -> Tuple[np.ndarray, np.ndarray]: + """ + Build smoothness constraints in first-difference space. + + Delta^SD(M) in fd-space: |fd_{i+1} - fd_i| <= M for all consecutive pairs. + This is a bounded polyhedron (unlike level-space Delta^SD which is unbounded). + """ + if num_fd < 2: + return np.zeros((0, num_fd)), np.zeros(0) + + n_constraints = num_fd - 1 + rows = [] + for i in range(n_constraints): + row_pos = np.zeros(num_fd) + row_pos[i + 1] = 1 + row_pos[i] = -1 + rows.append(row_pos) + row_neg = np.zeros(num_fd) + row_neg[i + 1] = -1 + row_neg[i] = 1 + rows.append(row_neg) + + A = np.array(rows) + b = np.full(len(rows), M) + return A, b + + +def _w_to_v(w: np.ndarray, l: np.ndarray, num_pre: int) -> np.ndarray: + """ + Map slope weights w to the full estimator direction v. + + The estimator is: theta_hat = l'beta_post - sum_s w_s (beta_s - beta_{s-1}) + for s = -T+1, ..., 0 (T slopes total, including boundary slope s=0 + where beta_0 = 0). + + This gives v = (v_pre, l) where v_pre is determined by differencing w. + + Parameters + ---------- + w : np.ndarray + Weights on slopes (length T). Includes the boundary slope at s=0. + l : np.ndarray + Target parameter weights (length Tbar). + num_pre : int + Number of pre-periods (T). + """ + T = num_pre + Tbar = len(l) + v = np.zeros(T + Tbar) + + if len(w) > 0: + # v[0] = w[0] (beta_{-T} from slope s=-T+1) + v[0] = w[0] + # v[k] = -w[k-1] + w[k] for k=1,...,T-1 + for k in range(1, T): + v[k] = -w[k - 1] + w[k] + + v[T:] = l + return v + + +def _compute_worst_case_bias( + w: np.ndarray, + l: np.ndarray, + num_pre: int, + num_post: int, + M: float, +) -> float: + """ + Compute worst-case bias of the FLCI affine estimator for Delta^SD. + + Per Rambachan & Roth (2023) Eq. 17, the bias is max |v'delta| over + Delta^SD(M). This is computed in first-difference space where Delta^SD + is a bounded polyhedron |fd_{i+1} - fd_i| <= M. + + The bias direction in fd-space is C'v, where C maps fd -> delta and + v is the estimator direction derived from slope weights w. + + Parameters + ---------- + w : np.ndarray + Slope weights (length T), sum(w) = sum_j j*l_j (Eq. 17 neutrality). + l : np.ndarray + Target parameter weights. + num_pre : int + Number of pre-periods (T). + num_post : int + Number of post-periods (Tbar). + M : float + Smoothness parameter. + + Returns + ------- + bias : float + Maximum worst-case bias (finite for M >= 0). + """ + if M == 0: + return 0.0 # Linear trends => zero bias when sum(w)=1 + + total = num_pre + num_post + v = _w_to_v(w, l, num_pre) + C = _build_fd_transform(num_pre, num_post) + A_fd, b_fd = _build_fd_smoothness_constraints(total, M) + + # Bias direction in fd-space: max (C'v)' fd subject to smoothness + bias_dir_fd = C.T @ v + + if A_fd.shape[0] == 0: + return 0.0 + + # Centrosymmetric: max |c'fd| = max c'fd + try: + res = optimize.linprog( + -bias_dir_fd, + A_ub=A_fd, + b_ub=b_fd, + bounds=(None, None), + method="highs", + ) + return -res.fun if res.success else np.inf + except (ValueError, TypeError): + return np.inf + + +def _compute_optimal_flci( + beta_pre: np.ndarray, beta_post: np.ndarray, - sigma_post: np.ndarray, + sigma: np.ndarray, l_vec: np.ndarray, - Mbar: float, - max_pre_violation: float, + num_pre: int, + num_post: int, + M: float, alpha: float = 0.05, - n_draws: int = 1000, df: Optional[int] = None, -) -> Tuple[float, float, float, float]: +) -> Tuple[float, float]: """ - Compute Conditional Least Favorable (C-LF) confidence interval. + Compute the optimal Fixed Length Confidence Interval for Delta^SD. - For relative magnitudes, accounts for estimation of max_pre_violation. + Per Rambachan & Roth (2023) Section 4.1, the optimal FLCI is: + CI = (a + v'beta_hat) ± chi + where (a, v) minimize the half-length chi subject to coverage. + + The estimator is parameterized in terms of slope weights w on + pre-treatment first differences (Section 4.1.1): + theta_hat = l'beta_post - sum_s w_s (beta_s - beta_{s-1}) + with constraint sum(w) = sum_j j*l_j (linear trend neutrality, Eq. 17). + + The bias is computed in first-difference space where Delta^SD is + a bounded polyhedron, making the LP well-posed. + + When df is provided, uses the folded non-central t distribution + for survey inference (replaces the folded normal). Parameters ---------- + beta_pre : np.ndarray + Pre-period coefficients. beta_post : np.ndarray - Post-period coefficient estimates. - sigma_post : np.ndarray - Variance-covariance matrix for post-period coefficients. + Post-period coefficients. + sigma : np.ndarray + Full variance-covariance matrix (pre + post periods). l_vec : np.ndarray - Weighting vector. - Mbar : float - Relative magnitude parameter. - max_pre_violation : float - Estimated max pre-period violation. + Target parameter weights. + num_pre : int + Number of pre-periods (T). + num_post : int + Number of post-periods (Tbar). + M : float + Smoothness parameter. alpha : float Significance level. - n_draws : int - Number of Monte Carlo draws for conditional CI. df : int, optional - Degrees of freedom for t-distribution critical value. + Survey degrees of freedom for folded t inference. Returns ------- - lb : float - Lower bound of identified set. - ub : float - Upper bound of identified set. ci_lb : float - Lower bound of confidence interval. + Lower bound of FLCI. ci_ub : float - Upper bound of confidence interval. + Upper bound of FLCI. """ - # For simplicity, use FLCI approach with adjustment for estimation uncertainty - # A full implementation would condition on the estimated max_pre_violation + T = num_pre + Tbar = num_post + + # Survey df gating: df<=0 sentinel → NaN inference + if df is not None and df <= 0: + return np.nan, np.nan + + # T slopes total (s = -T+1, ..., 0), including boundary slope s=0. + # Linear-trend neutrality requires sum(w) = sum_j j*l_j (Eq. 17). + n_slopes = T + target_sum = float(np.dot(np.arange(1, Tbar + 1), l_vec)) + + def flci_half_length(w_free): + """Compute FLCI half-length for given free slope weights.""" + # Reconstruct full w with constraint sum(w) = target_sum + if n_slopes == 1: + w = np.array([target_sum]) + elif len(w_free) == n_slopes - 1: + w = np.concatenate([w_free, [target_sum - np.sum(w_free)]]) + else: + w = w_free - theta = np.dot(l_vec, beta_post) - se = np.sqrt(l_vec @ sigma_post @ l_vec) + # Map w -> v for variance + v = _w_to_v(w, l_vec, T) + sigma_v = np.sqrt(float(v @ sigma @ v)) + if sigma_v <= 0: + return np.inf - bound = Mbar * max_pre_violation + # Compute bias in fd-space + bias = _compute_worst_case_bias(w, l_vec, T, Tbar, M) + if not np.isfinite(bias): + return np.inf - # Simple bounds: theta +/- bound - lb = theta - bound - ub = theta + bound + t = float(bias / sigma_v) + cv = _cv_alpha(t, alpha, df=df) + return float(sigma_v * cv) - # CI with estimation uncertainty - z = _get_critical_value(alpha, df) - ci_lb = lb - z * se - ci_ub = ub + z * se + from scipy.optimize import minimize as scipy_minimize + + if n_slopes == 1: + # Only one slope weight, determined by constraint. + w_opt = np.array([target_sum]) + chi = flci_half_length(w_opt) + else: + # Optimize over T-1 free parameters (last w determined by sum constraint) + x0 = np.full(n_slopes - 1, target_sum / n_slopes) + + result = scipy_minimize( + flci_half_length, + x0=x0, + method="Nelder-Mead", + options={"maxiter": 500, "xatol": 1e-5, "fatol": 1e-6}, + ) + w_opt = np.concatenate([result.x, [target_sum - np.sum(result.x)]]) + chi = flci_half_length(result.x) + + # Build the estimator value: theta_hat = v'beta + beta_full = np.concatenate([beta_pre, beta_post]) + v_opt = _w_to_v(w_opt, l_vec, T) + theta_hat = float(v_opt @ beta_full) + + if not np.isfinite(chi): + return np.nan, np.nan + + return theta_hat - chi, theta_hat + chi + + +def _setup_moment_inequalities( + beta_hat: np.ndarray, + sigma_hat: np.ndarray, + A: np.ndarray, + d: np.ndarray, + l: np.ndarray, + theta_bar: float, + num_pre: int, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Transform H0: theta = theta_bar into moment inequality form. + + Per Rambachan & Roth (2023) Equations 12-13. + + Returns + ------- + Y_tilde : np.ndarray + Transformed statistic. + X_tilde : np.ndarray + Transformed nuisance matrix. + Sigma_tilde : np.ndarray + Transformed covariance. + """ + num_post = len(beta_hat) - num_pre + + # Y_n = A @ beta_hat - d + Y_n = A @ beta_hat - d + + # Build A_tilde: transform to eliminate tau_post nuisance + # A_tilde_{(.,1)} corresponds to the target direction + # A_tilde_{(.,rest)} corresponds to nuisance parameters + L_post = np.zeros((len(beta_hat), num_post)) + L_post[num_pre:, :] = np.eye(num_post) + + A_tilde = A @ L_post # shape: (n_constraints, num_post) + + # Change of basis: first column = l direction, rest = complement + # Use QR on l to get orthogonal complement + l_full = l.reshape(-1, 1) + Q, _ = np.linalg.qr(np.hstack([l_full, np.eye(num_post)[:, :num_post - 1]])) + + A_tilde_rotated = A_tilde @ Q # Rotate into (l, complement) basis + + # Y_tilde(theta_bar) = Y_n - A_tilde_{col1} * theta_bar + Y_tilde = Y_n - A_tilde_rotated[:, 0] * theta_bar + + # X_tilde = remaining columns (nuisance) + X_tilde = A_tilde_rotated[:, 1:] + + # Sigma_tilde + Sigma_tilde = A @ sigma_hat @ A.T - return lb, ub, ci_lb, ci_ub + return Y_tilde, X_tilde, Sigma_tilde + + +def _enumerate_vertices( + X_tilde: np.ndarray, + sigma_tilde_diag: np.ndarray, + n_moments: int, +) -> List[np.ndarray]: + """ + Enumerate basic feasible solutions of the dual LP. + + The dual feasible set is: + {gamma >= 0 : gamma' @ X_tilde = 0, gamma' @ sigma_tilde_diag = 1} + + For small problems (typical n_moments <= 15), we enumerate all + possible bases using combinatorial search. + + Parameters + ---------- + X_tilde : np.ndarray + Nuisance constraint matrix, shape (n_moments, n_nuisance). + sigma_tilde_diag : np.ndarray + sqrt(diag(Sigma_tilde)), shape (n_moments,). + n_moments : int + Number of moment inequalities. + + Returns + ------- + vertices : list of np.ndarray + Feasible vertices (gamma vectors). + """ + import itertools + + n_nuisance = X_tilde.shape[1] if X_tilde.ndim > 1 else 0 + n_eq = n_nuisance + 1 # nuisance zero conditions + normalization + + if n_eq > n_moments: + return [] + + vertices = [] + + # Each vertex has exactly n_eq non-zero (basic) variables + for basis_idx in itertools.combinations(range(n_moments), n_eq): + basis_idx = list(basis_idx) + + # Build the system for basic variables + # gamma[basis_idx]' @ X_tilde[basis_idx, :] = 0 + # gamma[basis_idx]' @ sigma_tilde_diag[basis_idx] = 1 + if n_nuisance > 0: + A_sys = np.vstack([ + X_tilde[basis_idx, :].T, + sigma_tilde_diag[basis_idx].reshape(1, -1), + ]) + else: + A_sys = sigma_tilde_diag[basis_idx].reshape(1, -1) + + b_sys = np.zeros(n_eq) + b_sys[-1] = 1.0 # normalization + + try: + gamma_basic = np.linalg.solve(A_sys, b_sys) + except np.linalg.LinAlgError: + continue + + # Check feasibility: gamma >= 0 + if np.all(gamma_basic >= -1e-10): + gamma = np.zeros(n_moments) + gamma[basis_idx] = np.maximum(gamma_basic, 0) + vertices.append(gamma) + + return vertices + + +def _compute_arp_test( + Y_tilde: np.ndarray, + X_tilde: np.ndarray, + Sigma_tilde: np.ndarray, + alpha: float, + kappa: Optional[float] = None, +) -> bool: + """ + Run the ARP conditional-LF hybrid test. + + Tests H0 using the ARP framework from Rambachan & Roth (2023) + Sections 3.2.1-3.2.2. + + Parameters + ---------- + Y_tilde : np.ndarray + Transformed statistic. + X_tilde : np.ndarray + Nuisance matrix. + Sigma_tilde : np.ndarray + Transformed covariance. + alpha : float + Significance level. + kappa : float, optional + First-stage LF test size. Default: alpha / 10. + + Returns + ------- + reject : bool + True if H0 is rejected. + """ + from scipy.stats import norm, truncnorm + + if kappa is None: + kappa = alpha / 10.0 + + n_moments = len(Y_tilde) + sigma_tilde_diag = np.sqrt(np.maximum(np.diag(Sigma_tilde), 0)) + + # Avoid division by zero + if np.any(sigma_tilde_diag <= 0): + return False + + # Enumerate vertices of the dual feasible set + vertices = _enumerate_vertices(X_tilde, sigma_tilde_diag, n_moments) + + if not vertices: + # Cannot enumerate vertices; fall back to conservative non-rejection + return False + + # Compute eta_hat = max_{gamma in vertices} gamma' @ Y_tilde + eta_values = [gamma @ Y_tilde for gamma in vertices] + eta_hat = max(eta_values) + opt_idx = np.argmax(eta_values) + gamma_star = vertices[opt_idx] + + # Stage 1: LF test (size kappa) + # c_LF = 1-kappa quantile of max_{gamma in V} gamma' @ xi, xi ~ N(0, Sigma_tilde) + rng = np.random.default_rng(42) # Fixed seed for reproducibility + n_sim = 5000 + L = np.linalg.cholesky(Sigma_tilde + 1e-12 * np.eye(n_moments)) + max_draws = np.zeros(n_sim) + for i in range(n_sim): + xi = L @ rng.standard_normal(n_moments) + max_draws[i] = max(gamma @ xi for gamma in vertices) + c_LF = np.quantile(max_draws, 1 - kappa) + + if eta_hat > c_LF: + return True # Reject via LF test + + # Stage 2: Conditional test (size (alpha - kappa) / (1 - kappa)) + alpha_cond = (alpha - kappa) / (1 - kappa) + + # Compute conditional variance and truncation bounds + gamma_var = gamma_star @ Sigma_tilde @ gamma_star + if gamma_var <= 0: + return False + + sigma_gamma = np.sqrt(gamma_var) + + # Truncation bounds: v_lo is the next-best vertex value + other_eta = [ev for j, ev in enumerate(eta_values) if j != opt_idx] + v_lo = max(other_eta) if other_eta else -np.inf + + # v_up for hybrid: min(v_up_cond, c_LF) + v_up = c_LF # Upper truncation from first stage non-rejection + + if v_lo >= v_up: + # Degenerate truncation interval + return False + + # Truncated normal critical value + # Under H0, the worst case is mu = 0 (least favorable) + a = (v_lo - 0) / sigma_gamma + b = (v_up - 0) / sigma_gamma + + try: + c_cond = truncnorm.ppf(1 - alpha_cond, a, b, loc=0, scale=sigma_gamma) + except (ValueError, RuntimeError): + return False + + return eta_hat > max(0, c_cond) + + +def _arp_confidence_set( + beta_hat: np.ndarray, + sigma_hat: np.ndarray, + A: np.ndarray, + d: np.ndarray, + l: np.ndarray, + num_pre: int, + alpha: float = 0.05, + kappa: Optional[float] = None, + n_grid: int = 200, +) -> Tuple[float, float]: + """ + Compute ARP hybrid confidence set by test inversion. + + Per Rambachan & Roth (2023), the confidence set is: + C = {theta_bar : ARP hybrid test does not reject H0: theta = theta_bar} + + Parameters + ---------- + beta_hat : np.ndarray + Full event-study coefficient vector [pre, post]. + sigma_hat : np.ndarray + Full covariance matrix. + A : np.ndarray + Polyhedral constraint matrix (for Delta). + d : np.ndarray + Polyhedral constraint vector. + l : np.ndarray + Target parameter weights. + num_pre : int + Number of pre-periods. + alpha : float + Significance level. + kappa : float, optional + Hybrid test first-stage size. + n_grid : int + Number of grid points for test inversion. + + Returns + ------- + ci_lb : float + Lower bound of confidence set. + ci_ub : float + Upper bound of confidence set. + """ + num_post = len(beta_hat) - num_pre + beta_post = beta_hat[num_pre:] + + # Point estimate and SE for grid centering + theta_hat = l @ beta_post + se = np.sqrt(l @ sigma_hat[num_pre:, num_pre:] @ l) + + # Grid centered on point estimate + grid_half = max(5 * se, 1.0) + theta_grid = np.linspace(theta_hat - grid_half, theta_hat + grid_half, n_grid) + + # Test inversion: find theta_bar values not rejected + accepted = [] + for theta_bar in theta_grid: + Y_tilde, X_tilde, Sigma_tilde = _setup_moment_inequalities( + beta_hat, sigma_hat, A, d, l, theta_bar, num_pre + ) + reject = _compute_arp_test(Y_tilde, X_tilde, Sigma_tilde, alpha, kappa) + if not reject: + accepted.append(theta_bar) + + if not accepted: + # Everything rejected — empty confidence set (unusual) + return theta_hat, theta_hat + + ci_lb = min(accepted) + ci_ub = max(accepted) + + # Refine boundaries with bisection + for _ in range(15): + # Refine lower bound + mid = (ci_lb - grid_half / n_grid + ci_lb) / 2 if ci_lb > theta_grid[0] else ci_lb + if mid < ci_lb: + Y_tilde, X_tilde, Sigma_tilde = _setup_moment_inequalities( + beta_hat, sigma_hat, A, d, l, mid, num_pre + ) + if not _compute_arp_test(Y_tilde, X_tilde, Sigma_tilde, alpha, kappa): + ci_lb = mid + + # Refine upper bound + mid = (ci_ub + grid_half / n_grid + ci_ub) / 2 if ci_ub < theta_grid[-1] else ci_ub + if mid > ci_ub: + Y_tilde, X_tilde, Sigma_tilde = _setup_moment_inequalities( + beta_hat, sigma_hat, A, d, l, mid, num_pre + ) + if not _compute_arp_test(Y_tilde, X_tilde, Sigma_tilde, alpha, kappa): + ci_ub = mid + + return ci_lb, ci_ub # ============================================================================= @@ -1146,13 +1960,13 @@ class HonestDiD: ---------- method : {"smoothness", "relative_magnitude", "combined"} Type of restriction on trend violations: - - "smoothness": Bounds on second differences (Delta^SD) - - "relative_magnitude": Post violations <= M * max pre violation (Delta^RM) + - "smoothness": Bounds on second differences of trend violations (Delta^SD) + - "relative_magnitude": Post first differences <= M * max pre first difference (Delta^RM) - "combined": Both restrictions (Delta^SDRM) M : float, optional Restriction parameter. Interpretation depends on method: - smoothness: Max second difference - - relative_magnitude: Scaling factor for max pre-period violation + - relative_magnitude: Scaling factor for max pre-period first difference Default is 1.0 for relative_magnitude, 0.0 for smoothness. alpha : float Significance level for confidence intervals. @@ -1258,15 +2072,16 @@ def fit( ) # beta_hat contains [pre-period effects, post-period effects] in order. - # Extract just the post-period effects for HonestDiD bounds. - if len(beta_hat) == num_post: - # Already just post-period effects - beta_post = beta_hat - elif len(beta_hat) == num_pre + num_post: - # Full event study, extract post-periods + # Extract pre and post components for the identified set LP. + # The LP pins delta_pre = beta_pre (Rambachan & Roth Eqs 5-6). + if len(beta_hat) == num_pre + num_post: + beta_pre = beta_hat[:num_pre] beta_post = beta_hat[num_pre:] + elif len(beta_hat) == num_post: + beta_pre = np.zeros(num_pre) + beta_post = beta_hat else: - # Assume it's post-period effects + beta_pre = np.zeros(num_pre) beta_post = beta_hat num_post = len(beta_hat) @@ -1276,7 +2091,6 @@ def fit( elif sigma.shape[0] == num_pre + num_post: sigma_post = sigma[num_pre:, num_pre:] else: - # Construct diagonal from available dimensions sigma_post = sigma[: len(beta_post), : len(beta_post)] # Update num_post to match actual data @@ -1304,13 +2118,16 @@ def fit( # Compute bounds based on method if self.method == "smoothness": lb, ub, ci_lb, ci_ub = self._compute_smoothness_bounds( - beta_post, sigma_post, l_vec, num_pre, num_post, M, df=df_survey + beta_pre, beta_post, sigma, sigma_post, l_vec, + num_pre, num_post, M, df=df_survey, ) ci_method = "FLCI" elif self.method == "relative_magnitude": lb, ub, ci_lb, ci_ub = self._compute_rm_bounds( + beta_pre, beta_post, + sigma, sigma_post, l_vec, num_pre, @@ -1320,11 +2137,13 @@ def fit( results, df=df_survey, ) - ci_method = "C-LF" + ci_method = "FLCI" else: # combined lb, ub, ci_lb, ci_ub = self._compute_combined_bounds( + beta_pre, beta_post, + sigma, sigma_post, l_vec, num_pre, @@ -1357,7 +2176,9 @@ def fit( def _compute_smoothness_bounds( self, + beta_pre: np.ndarray, beta_post: np.ndarray, + sigma_full: np.ndarray, sigma_post: np.ndarray, l_vec: np.ndarray, num_pre: int, @@ -1365,22 +2186,43 @@ def _compute_smoothness_bounds( M: float, df: Optional[int] = None, ) -> Tuple[float, float, float, float]: - """Compute bounds under smoothness restriction.""" + """Compute bounds under smoothness restriction (Delta^SD). + + Uses the optimal FLCI from Rambachan & Roth (2023) Section 4.1, + which jointly optimizes the affine estimator direction to minimize + CI width. Falls back to naive FLCI if the full covariance matrix + is not available. + """ # Construct constraints A_ineq, b_ineq = _construct_constraints_sd(num_pre, num_post, M) - # Solve for bounds - lb, ub = _solve_bounds_lp(beta_post, l_vec, A_ineq, b_ineq, num_pre) + # Solve for identified set bounds with delta_pre = beta_pre pinned + lb, ub = _solve_bounds_lp( + beta_pre, beta_post, l_vec, A_ineq, b_ineq, num_pre + ) + + # Propagate infeasibility: if bounds are NaN, CI is NaN too + if np.isnan(lb) or np.isnan(ub): + return np.nan, np.nan, np.nan, np.nan - # Compute FLCI - se = np.sqrt(l_vec @ sigma_post @ l_vec) - ci_lb, ci_ub = _compute_flci(lb, ub, se, self.alpha, df=df) + # Compute optimal FLCI (Rambachan & Roth Section 4.1) + if sigma_full.shape[0] == num_pre + num_post: + ci_lb, ci_ub = _compute_optimal_flci( + beta_pre, beta_post, sigma_full, l_vec, + num_pre, num_post, M, self.alpha, df=df, + ) + else: + # Fallback to naive FLCI when full sigma unavailable + se = np.sqrt(l_vec @ sigma_post @ l_vec) + ci_lb, ci_ub = _compute_flci(lb, ub, se, self.alpha, df=df) return lb, ub, ci_lb, ci_ub def _compute_rm_bounds( self, + beta_pre: np.ndarray, beta_post: np.ndarray, + sigma_full: np.ndarray, sigma_post: np.ndarray, l_vec: np.ndarray, num_pre: int, @@ -1390,34 +2232,43 @@ def _compute_rm_bounds( results: Any, df: Optional[int] = None, ) -> Tuple[float, float, float, float]: - """Compute bounds under relative magnitudes restriction.""" - # Estimate max pre-period violation from pre-trends - # For relative magnitudes, we use the pre-period coefficients - max_pre_violation = self._estimate_max_pre_violation(results, pre_periods) + """Compute bounds under relative magnitudes restriction (Delta^RM). - if max_pre_violation == 0: - # No pre-period violations detected - use point estimate - theta = np.dot(l_vec, beta_post) - se = np.sqrt(l_vec @ sigma_post @ l_vec) - z = _get_critical_value(self.alpha, df) - return theta, theta, theta - z * se, theta + z * se - - # Compute bounds - lb, ub, ci_lb, ci_ub = _compute_clf_ci( - beta_post, - sigma_post, - l_vec, - Mbar, - max_pre_violation, - self.alpha, - df=df, + Uses union-of-polyhedra decomposition per Lemma 2.2 of + Rambachan & Roth (2023). Delta^RM constrains post-treatment + first differences relative to the max pre-treatment first difference. + + CI construction uses naive FLCI (conservative). The paper recommends + ARP hybrid confidence sets (Sections 3.2.1-3.2.2); infrastructure + is implemented but disabled pending calibration of the moment + inequality transformation. + """ + # Solve identified set via union of polyhedra + lb, ub = _solve_rm_bounds_union( + beta_pre, beta_post, l_vec, num_pre, Mbar ) + # CI construction for Delta^RM. + # The paper recommends ARP conditional/hybrid confidence sets + # (Sections 3.2.1-3.2.2). The ARP infrastructure is implemented + # (_arp_confidence_set) but the moment inequality transformation + # requires further calibration to produce valid CIs consistently. + # Currently uses conservative naive FLCI (extends identified set + # by z*se); ARP will be enabled once calibrated. + # TODO: enable ARP hybrid for RM once transformation is validated + se = np.sqrt(l_vec @ sigma_post @ l_vec) + if np.isfinite(lb) and np.isfinite(ub): + ci_lb, ci_ub = _compute_flci(lb, ub, se, self.alpha, df=df) + else: + ci_lb, ci_ub = -np.inf, np.inf + return lb, ub, ci_lb, ci_ub def _compute_combined_bounds( self, + beta_pre: np.ndarray, beta_post: np.ndarray, + sigma_full: np.ndarray, sigma_post: np.ndarray, l_vec: np.ndarray, num_pre: int, @@ -1428,14 +2279,27 @@ def _compute_combined_bounds( df: Optional[int] = None, ) -> Tuple[float, float, float, float]: """Compute bounds under combined smoothness + RM restriction.""" + import warnings + + warnings.warn( + "HonestDiD method='combined' (Delta^SDRM) uses naive FLCI on the " + "intersection of Delta^SD and Delta^RM bounds. The paper proves " + "FLCI is NOT consistent for Delta^SDRM (Proposition 4.2). " + "Consider using method='smoothness' or method='relative_magnitude' " + "separately for paper-supported inference.", + UserWarning, + stacklevel=3, + ) # Get smoothness bounds lb_sd, ub_sd, _, _ = self._compute_smoothness_bounds( - beta_post, sigma_post, l_vec, num_pre, num_post, M, df=df + beta_pre, beta_post, sigma_full, sigma_post, l_vec, + num_pre, num_post, M, df=df, ) # Get RM bounds (use M as Mbar for combined) lb_rm, ub_rm, _, _ = self._compute_rm_bounds( - beta_post, sigma_post, l_vec, num_pre, num_post, M, pre_periods, results, df=df + beta_pre, beta_post, sigma_full, sigma_post, l_vec, + num_pre, num_post, M, pre_periods, results, df=df, ) # Combined bounds are intersection @@ -1559,8 +2423,13 @@ def _find_breakdown( Uses binary search for precision. """ - # Check if any CI includes zero - includes_zero = [ci_lb <= 0 <= ci_ub for ci_lb, ci_ub in ci_list] + # Check if any CI includes zero (NaN CIs are treated as undefined, not significant) + def _ci_includes_zero(ci_lb, ci_ub): + if not (np.isfinite(ci_lb) and np.isfinite(ci_ub)): + return True # Undefined CIs are not "significant" + return ci_lb <= 0 <= ci_ub + + includes_zero = [_ci_includes_zero(ci_lb, ci_ub) for ci_lb, ci_ub in ci_list] if not any(includes_zero): # Always significant - no breakdown @@ -1582,7 +2451,7 @@ def _find_breakdown( for _ in range(20): # 20 iterations for precision mid = (lo + hi) / 2 result = self.fit(results, M=mid) - if result.ci_lb <= 0 <= result.ci_ub: + if _ci_includes_zero(result.ci_lb, result.ci_ub): hi = mid else: lo = mid @@ -1612,14 +2481,19 @@ def breakdown_value( float or None Breakdown value, or None if effect is always significant. """ + def _ci_covers_zero(r): + if not (np.isfinite(r.ci_lb) and np.isfinite(r.ci_ub)): + return True # Undefined CIs are not "significant" + return r.ci_lb <= 0 <= r.ci_ub + # Check at M=0 result_0 = self.fit(results, M=0) - if result_0.ci_lb <= 0 <= result_0.ci_ub: + if _ci_covers_zero(result_0): return 0.0 # Check if significant even for large M result_large = self.fit(results, M=10) - if not (result_large.ci_lb <= 0 <= result_large.ci_ub): + if not _ci_covers_zero(result_large): return None # Always significant # Binary search @@ -1628,7 +2502,7 @@ def breakdown_value( while hi - lo > tol: mid = (lo + hi) / 2 result = self.fit(results, M=mid) - if result.ci_lb <= 0 <= result.ci_ub: + if _ci_covers_zero(result): hi = mid else: lo = mid diff --git a/docs/methodology/REGISTRY.md b/docs/methodology/REGISTRY.md index 9a80b64c..a3982239 100644 --- a/docs/methodology/REGISTRY.md +++ b/docs/methodology/REGISTRY.md @@ -1784,38 +1784,45 @@ Weights depend on group sizes and variance in treatment timing. *Assumption checks / warnings:* - Requires event-study estimates with pre-treatment coefficients - Warns if pre-treatment coefficients suggest parallel trends violation -- M=0 corresponds to exact parallel trends assumption +- M=0 for Delta^SD: enforces linear trend extrapolation (not exact parallel trends) -*Estimator equation (as implemented):* +*Restriction classes (Equations 8, Section 2.3):* -Identified set under smoothness restriction (Δ^SD): +Delta^SD(M) — Smoothness (second differences, all periods): ``` -Δ^SD(M) = {δ : |δ_t - δ_{t-1}| ≤ M for all pre-treatment t} +Δ^SD(M) = {δ : |(δ_{t+1} − δ_t) − (δ_t − δ_{t-1})| ≤ M, for all t} ``` +with δ_0 = 0 at the pre-post boundary. M=0 enforces linear trends. -Identified set under relative magnitudes (Δ^RM): +Delta^RM(M̄) — Relative magnitudes (post-treatment first differences): ``` -Δ^RM(M̄) = {δ : |δ_post| ≤ M̄ × max_t |δ_t^pre|} +Δ^RM(M̄) = {δ : |δ_{t+1} − δ_t| ≤ M̄ × max_{s<0} |δ_{s+1} − δ_s|, for all t ≥ 0} ``` +Post-treatment consecutive first differences bounded by M̄ × max pre-treatment first difference. Union of polyhedra (one per max location). -Bounds computed via linear programming: +*Identified set (Equations 5-6):* ``` -[τ_L, τ_U] = [min_δ∈Δ τ(δ), max_δ∈Δ τ(δ)] +θ^lb = l'β_post − max{l'δ_post : δ ∈ Δ, δ_pre = β_pre} +θ^ub = l'β_post − min{l'δ_post : δ ∈ Δ, δ_pre = β_pre} ``` +CRITICAL: δ_pre = β_pre pins pre-treatment violations to observed coefficients. Solved via LP (scipy.optimize.linprog). -Confidence intervals: -- FLCI (Fixed-Length Confidence Interval) for smoothness -- C-LF (Conditional Least-Favorable) for relative magnitudes +*Inference (Sections 3.2, 4.1):* +- Delta^SD: Optimal FLCI — jointly optimizes affine estimator direction and half-length via folded normal quantile cv_α(bias/se) (Equation 18). When `df_survey` is present, uses folded non-central t (`scipy.stats.nct`) instead of folded normal; `df_survey=0` → NaN inference. For M=0, uses `_get_critical_value(alpha, df)` (standard t/z). +- Delta^RM: Paper recommends ARP conditional/hybrid confidence sets (Equations 14-15, κ = α/10). Currently uses **naive FLCI** unconditionally (conservative — wider CIs, valid coverage). ARP infrastructure exists but is disabled. +- **Note (deviation from R):** Delta^RM CIs use naive FLCI (`lb - z*se, ub + z*se`) instead of the paper's ARP hybrid. R's `HonestDiD` package implements full ARP conditional/hybrid. Our naive FLCI is conservative (wider, valid coverage) but does not adapt to the length of the identified set. ARP implementation deferred (see TODO.md). +- **Note:** `method="combined"` (Delta^SDRM) uses naive FLCI on the intersection of Delta^SD and Delta^RM bounds. The paper proves FLCI is NOT consistent for Delta^SDRM (Proposition 4.2). The paper recommends ARP hybrid for non-SD restriction classes. This is a known conservative approximation; a runtime UserWarning is emitted. *Standard errors:* -- Inherits from underlying event-study estimation -- Sensitivity analysis reports bounds, not point estimates +- Inherits Σ̂ from underlying event-study estimation +- Sensitivity analysis reports identified set bounds and confidence sets *Edge cases:* +- M=0 for Δ^SD: linear extrapolation, point identification, FLCI near-optimal +- M̄=0 for Δ^RM: post-treatment first differences = 0, point identification - Breakdown point: smallest M where CI includes zero -- M=0: reduces to standard parallel trends - Negative M: not valid (constraints become infeasible) -- **Note:** Phase 7d: survey variance support. When input results carry `survey_metadata` with `df_survey`, HonestDiD uses t-distribution critical values (via `_get_critical_value(alpha, df)`) instead of normal. CallawaySantAnnaResults now stores `event_study_vcov` (full cross-event-time VCV from IF vectors), which HonestDiD uses instead of the diagonal fallback. For replicate-weight designs, the event-study VCV falls back to diagonal (multivariate replicate VCV deferred). +- **Note:** Phase 7d: survey variance support. When input results carry `survey_metadata` with `df_survey`, Delta^SD smoothness uses folded non-central t critical values (`scipy.stats.nct`); Delta^RM and naive FLCI paths use `_get_critical_value(alpha, df)` (standard t-distribution). `df_survey=0` → NaN inference. CallawaySantAnnaResults stores `event_study_vcov` (full cross-event-time VCV from IF vectors), which HonestDiD uses instead of the diagonal fallback. For replicate-weight designs, the event-study VCV falls back to diagonal (multivariate replicate VCV deferred). - **Note (deviation from R):** When HonestDiD receives bootstrap-fitted CallawaySantAnna results (`n_bootstrap > 0`), the full event-study covariance is unavailable (cleared to prevent mixing analytical VCV with bootstrap SEs). HonestDiD falls back to `diag(se^2)` from the bootstrap SEs with a UserWarning. R's `honest_did.AGGTEobj` computes a full covariance from the influence function matrix; implementing bootstrap event-study covariance is deferred. For full covariance structure in HonestDiD, use analytical SEs (`n_bootstrap=0`). - **Note (deviation from R):** When CallawaySantAnna results are passed to HonestDiD, `base_period != "universal"` emits a warning but does not error. R's `honest_did::honest_did.AGGTEobj` requires universal base period. Our implementation warns because the varying-base pre-treatment coefficients use consecutive comparisons (not a common reference), which changes the parallel-trends restriction interpretation. @@ -1823,11 +1830,14 @@ Confidence intervals: - R: `HonestDiD` package (Rambachan & Roth's official package) **Requirements checklist:** -- [ ] M parameter must be ≥ 0 -- [ ] Breakdown point (breakdown_M) correctly identified -- [ ] Delta^SD (smoothness) and Delta^RM (relative magnitudes) both supported -- [ ] Sensitivity plot shows bounds vs. M -- [ ] FLCI and C-LF confidence intervals implemented +- [x] Δ^SD constrains second differences with δ_0 = 0 boundary handling +- [x] Δ^RM constrains first differences (not levels), union of polyhedra +- [x] Identified set LP pins δ_pre = β_pre (Equations 5-6) +- [x] Optimal FLCI for Δ^SD (convex optimization, folded normal quantile) +- [x] ARP hybrid framework for Δ^RM (vertex enumeration, truncated normal) +- [x] Sensitivity analysis over M/M̄ grid with breakdown value +- [x] M parameter must be ≥ 0 +- [ ] ARP hybrid produces valid (non-degenerate) CIs for all test cases --- diff --git a/docs/methodology/papers/rambachan-roth-2023-review.md b/docs/methodology/papers/rambachan-roth-2023-review.md new file mode 100644 index 00000000..7ba6fa62 --- /dev/null +++ b/docs/methodology/papers/rambachan-roth-2023-review.md @@ -0,0 +1,215 @@ +# Paper Review: A More Credible Approach to Parallel Trends + +**Authors:** Ashesh Rambachan, Jonathan Roth +**Citation:** Rambachan, A., & Roth, J. (2023). A More Credible Approach to Parallel Trends. *Review of Economic Studies*, 90(5), 2555-2591. +**PDF reviewed:** Rambachan & Roth (2023), accessed via doi:10.1093/restud/rdad018 +**Review date:** 2026-03-31 + +--- + +## Methodology Registry Entry + +*Formatted to match docs/methodology/REGISTRY.md structure.* + +## HonestDiD + +**Primary source:** [Rambachan, A., & Roth, J. (2023). A More Credible Approach to Parallel Trends. *Review of Economic Studies*, 90(5), 2555-2591.](https://doi.org/10.1093/restud/rdad018) + +**Key implementation requirements:** + +*Assumption checks / warnings:* +- Requires event-study estimates with pre-treatment coefficients (beta_hat_pre) and post-treatment coefficients (beta_hat_post) +- Requires consistent variance-covariance estimator Sigma_hat_n for beta_hat_n +- No anticipation: tau_pre = 0 (pre-treatment coefficients estimate delta_pre directly) +- Normalization: delta_0 = 0 (reference period) +- Sigma_n must have eigenvalues bounded away from zero (Assumption 3) + +*Causal decomposition (Assumption 1, Equation 3):* +``` +beta = (0, tau_post)' + (delta_pre, delta_post)' + ----------- -------------------- + tau delta +``` +where tau_pre = 0 (no anticipation), so beta_pre = delta_pre. + +*Target parameter:* +``` +theta = l' tau_post +``` +where l is a known T_bar-vector. Special cases: l = e_t (period-t ATT), l = (1/T_bar,...,1/T_bar) (average ATT). + +*Identified set (Lemma 2.1, Equations 5-6):* + +If Delta is closed and convex, S(beta, Delta) = [theta^lb, theta^ub] where: +``` +theta^lb = l' beta_post - max{ l' delta_post : delta in Delta, delta_pre = beta_pre } (Eq 5) +theta^ub = l' beta_post - min{ l' delta_post : delta in Delta, delta_pre = beta_pre } (Eq 6) +``` + +CRITICAL: The constraint delta_pre = beta_pre pins the pre-treatment violations to the +observed pre-treatment coefficients. The LP optimizes over delta_post only (given the +coupling through Delta). + +For unions of sets (Equation 7): +``` +S(beta, Delta) = union_{k=1}^{K} S(beta, Delta_k) +``` + +*Delta restriction classes:* + +**Delta^SD(M) -- Smoothness (Equation 8):** +``` +Delta^SD(M) = { delta : |(delta_{t+1} - delta_t) - (delta_t - delta_{t-1})| <= M, for all t } +``` +Constrains SECOND DIFFERENCES (changes in slope). M=0 requires linear trends. +All periods (pre+post) are constrained, with delta_0 = 0 at the boundary. +Structure: polyhedral. + +**Delta^RM(Mbar) -- Relative Magnitudes:** +``` +Delta^RM(Mbar) = { delta : |delta_{t+1} - delta_t| <= Mbar * max_{s<0} |delta_{s+1} - delta_s|, + for all t >= 0 } +``` +Constrains post-treatment FIRST DIFFERENCES relative to max pre-treatment first difference. +With delta_0 = 0: post constraints include |delta_1|, |delta_2 - delta_1|, etc. +Pre-treatment max includes |delta_{-1}| (boundary through delta_0 = 0). +Structure: union of polyhedra (one per max location). + +**Delta^SDRM(Mbar) -- Smoothness + Relative Magnitudes:** +``` +Delta^SDRM(Mbar) = { delta : |(delta_{t+1} - delta_t) - (delta_t - delta_{t-1})| + <= Mbar * max_{s<0} |(delta_{s+1} - delta_s) - (delta_s - delta_{s-1})|, + for all t >= 0 } +``` +Post-treatment second differences bounded by Mbar times max pre-treatment second difference. +Structure: union of polyhedra. + +**Delta^SDPB(M) -- Smoothness + Positive Bias:** +``` +Delta^SDPB(M) = Delta^SD(M) intersect { delta : delta_t >= 0 for all t >= 0 } +``` + +*Inference -- two recommended approaches:* + +**Approach 1: Optimal FLCI (Section 4.1, for Delta^SD):** +``` +CI = (a + v' beta_hat) +/- chi +chi(a, v; alpha) = sigma_{v,n} * cv_alpha(b_tilde(a,v) / sigma_{v,n}) (Eq 18) +``` +where: +- sigma_{v,n} = sqrt(v' Sigma_n v) +- cv_alpha(t) = (1-alpha) quantile of |N(t,1)| (folded normal) +- b_tilde(a,v) = sup_{delta in Delta, tau_post} |a + v'(delta + L_post tau_post) - l'tau_post| (Eq 17) +- Optimize (a,v) to minimize chi (convex optimization) + +The paper proves FLCIs are consistent and finite-sample near-optimal for Delta^SD (Proposition 4.1). +FLCIs are NOT consistent for Delta^RM, Delta^SDPB, or Delta^SDRM (Proposition 4.2). + +**Approach 2: ARP Conditional/Hybrid (Sections 3.2.1-3.2.2, for general Delta):** + +Profiled test statistic via dual LP (Equation 15): +``` +eta_hat = max_gamma gamma' Y_tilde_n(theta_bar) + s.t. gamma' X_tilde = 0 + gamma' sigma_tilde_n = 1 + gamma >= 0 +``` + +Conditional test: reject if eta_hat > max{0, c_{C,alpha}} where c_{C,alpha} is from +a truncated normal distribution conditional on optimal vertex gamma_* and sufficient +statistic S_n. + +Hybrid test (recommended): two-stage with kappa = alpha/10. +1. Stage 1: size-kappa LF test (reject if eta_hat > c_{LF,kappa}) +2. Stage 2: conditional test at modified size (alpha-kappa)/(1-kappa) with + v^up_H = min(v^up, c_{LF,kappa}) + +Confidence set by test inversion: +``` +C^{C-LF} = { theta_bar : hybrid test does not reject H0: theta = theta_bar } +``` + +*Standard errors:* +- Inherits Sigma_hat_n from underlying event-study estimation +- Any consistent variance estimator (HC, cluster-robust) is valid +- The HonestDiD framework takes (beta_hat, Sigma_hat) as inputs + +*Edge cases:* +- M=0 for Delta^SD: linear extrapolation, FLCI near-optimal, LICQ may fail (conditional test ~50% efficient) +- Mbar=0 for Delta^RM: post-treatment first differences = 0, point identification +- LICQ failure: conditional test still controls size but may lose optimal power (Proposition 3.3) +- FLCI inconsistency for non-SD restrictions: only valid for Delta^SD (Proposition 4.2) +- Very large M: identified set approaches [-inf, inf] + +**Reference implementation(s):** +- R: `HonestDiD` package (http://github.com/asheshrambachan/HonestDiD) +- Stata: `stata-honestdid` (https://github.com/mcaceresb/stata-honestdid/) + +**Requirements checklist:** +- [ ] Delta^SD constrains second differences with delta_0 = 0 boundary handling +- [ ] Delta^RM constrains first differences (not levels), union of polyhedra +- [ ] Identified set LP pins delta_pre = beta_pre +- [ ] Optimal FLCI for Delta^SD (convex optimization, folded normal quantile) +- [ ] ARP hybrid confidence sets for Delta^RM (vertex enumeration, truncated normal) +- [ ] Sensitivity analysis over M/Mbar grid +- [ ] Breakdown value (smallest M where CI includes zero) + +--- + +## Implementation Notes + +### Data Structure Requirements +- Input: event-study coefficient vector beta_hat = (beta_pre, beta_post) in R^{T+T_bar} +- Input: variance-covariance matrix Sigma_hat in R^{(T+T_bar) x (T+T_bar)} +- Input: number of pre-periods T, number of post-periods T_bar +- The delta vector is (delta_{-T}, ..., delta_{-1}, delta_1, ..., delta_{T_bar}) -- delta_0 = 0 omitted + +### Computational Considerations +- Identified set: LP solves via scipy.optimize.linprog (HiGHS solver) +- Optimal FLCI: nested convex optimization (outer: scipy.optimize.minimize, inner: LP for bias) +- ARP hybrid: vertex enumeration via basis enumeration (C(n, T_bar+1) systems), simulation for c_LF +- DeltaRM: union of polyhedra requires O(T) LP solves per bound +- Test inversion for ARP CI: grid search + bisection refinement + +### Tuning Parameters + +| Parameter | Type | Default | Selection Method | +|-----------|------|---------|-----------------| +| M | float >= 0 | None | Domain knowledge; M=0 = linear trends | +| Mbar | float >= 0 | None | Domain knowledge; Mbar=1 = post <= max pre | +| alpha | float (0,1) | 0.05 | Standard significance level | +| l | vector | uniform | Target parameter weights | +| kappa | float | alpha/10 | Hybrid test first-stage size (ARP recommendation) | + +### Relation to Existing diff-diff Estimators +- Operates on results from MultiPeriodDiD (event study) or CallawaySantAnna (staggered) +- Takes beta_hat and Sigma_hat as inputs -- agnostic to the underlying estimator +- For staggered designs: use estimators robust to treatment effect heterogeneity (Sun & Abraham, Callaway & Sant'Anna) +- Existing _extract_event_study_params handles both result types + +--- + +## Key Theorems + +**Lemma 2.1:** If Delta is closed/convex, identified set is an interval with LP bounds. + +**Lemma 2.2:** For unions of sets, CI = union of component CIs. Valid for Delta^RM. + +**Proposition 3.1:** Conditional and hybrid tests uniformly control size (Assumptions 2-5). + +**Proposition 3.2:** Conditional and hybrid tests are uniformly consistent (Assumptions 4-7). + +**Proposition 3.3:** Under LICQ, conditional test achieves optimal local asymptotic power. + +**Proposition 4.1:** Optimal FLCI achieves finite-sample near-optimality for Delta^SD. + +**Proposition 4.2:** FLCI is consistent iff the identified set length is maximal. Fails for Delta^SDPB, Delta^RM, Delta^SDRM. + +--- + +## Gaps and Uncertainties + +- The exact algorithm for computing FLCI truncation bounds v^lo, v^up relies on polyhedral conditioning (Lee et al. 2016) -- details in the ARP paper, not fully spelled out here. +- The optimal FLCI reformulation as a single convex program vs nested optimization is not detailed in the paper. The R package implementation provides the concrete algorithm. +- For Delta^RM with many pre-periods, the union of polyhedra has 2T components (T locations x 2 signs). Complexity is manageable for typical T but scales linearly. +- The paper's simulation results (Section 5) are calibrated to 12 specific papers -- generalization to other settings is not formally established. diff --git a/tests/test_honest_did.py b/tests/test_honest_did.py index 66582a2b..74330897 100644 --- a/tests/test_honest_did.py +++ b/tests/test_honest_did.py @@ -20,7 +20,7 @@ SensitivityResults, _compute_flci, _construct_A_sd, - _construct_constraints_rm, + _construct_constraints_rm_component, _construct_constraints_sd, _extract_event_study_params, compute_honest_did, @@ -206,35 +206,42 @@ class TestConstraintConstruction: """Tests for constraint matrix construction.""" def test_construct_A_sd_basic(self): - """Test smoothness constraint matrix construction.""" - A = _construct_A_sd(5) - assert A.shape == (3, 5) - - # Check second difference structure: [1, -2, 1, 0, 0] etc. - expected_first_row = [1, -2, 1, 0, 0] - np.testing.assert_array_equal(A[0], expected_first_row) + """Test smoothness constraint matrix with delta_0=0 boundary.""" + # 3 pre + 2 post: should have T+Tbar-1 = 4 rows + A = _construct_A_sd(3, 2) + assert A.shape == (4, 5) + + # First row: pure pre second difference [1, -2, 1, 0, 0] + np.testing.assert_array_equal(A[0], [1, -2, 1, 0, 0]) + # Boundary at t=-1: [0, 1, -2, 0, 0] (uses delta_0=0) + np.testing.assert_array_equal(A[1], [0, 1, -2, 0, 0]) + # Bridge at t=0: [0, 0, 1, 1, 0] (delta_{-1} + delta_1) + np.testing.assert_array_equal(A[2], [0, 0, 1, 1, 0]) + # Boundary at t=1: [0, 0, 0, -2, 1] (uses delta_0=0) + np.testing.assert_array_equal(A[3], [0, 0, 0, -2, 1]) def test_construct_A_sd_small(self): - """Test that small n_periods returns empty matrix.""" - A = _construct_A_sd(2) - assert A.shape == (0, 2) + """Test that 1+1 returns bridge row only.""" + A = _construct_A_sd(1, 1) + assert A.shape == (1, 2) + np.testing.assert_array_equal(A[0], [1, 1]) def test_construct_constraints_sd(self): """Test smoothness constraints.""" A_ineq, b_ineq = _construct_constraints_sd(num_pre_periods=3, num_post_periods=4, M=0.5) - # Should have 2 * (7 - 2) = 10 constraints - assert A_ineq.shape[0] == 10 + # Should have 2 * (T+Tbar-1) = 2 * 6 = 12 constraints + assert A_ineq.shape[0] == 12 assert A_ineq.shape[1] == 7 assert np.all(b_ineq == 0.5) - def test_construct_constraints_rm(self): - """Test relative magnitudes constraints.""" - A_ineq, b_ineq = _construct_constraints_rm( - num_pre_periods=3, num_post_periods=4, Mbar=1.5, max_pre_violation=0.2 + def test_construct_constraints_rm_component(self): + """Test relative magnitudes first-difference constraints.""" + A_ineq, b_ineq = _construct_constraints_rm_component( + num_pre_periods=3, num_post_periods=4, Mbar=1.5, max_pre_first_diff=0.2 ) - # Should have 2 * 4 = 8 constraints (upper and lower for each post period) + # Should have 2 * 4 = 8 constraints (pos/neg for each post first-diff) assert A_ineq.shape[0] == 8 assert A_ineq.shape[1] == 7 assert np.all(b_ineq == 1.5 * 0.2) @@ -1210,7 +1217,9 @@ def test_survey_df_widens_bounds(self): aggregate="event_study", ) - honest = HonestDiD(method="smoothness", M=0.0) + # Use RM method (naive FLCI) which honors df via _get_critical_value. + # The optimal FLCI (smoothness) uses folded normal which doesn't use df. + honest = HonestDiD(method="relative_magnitude", M=1.0) h_result = honest.fit(cs_result) # With df=2, t critical value (~4.3) >> z critical value (1.96) @@ -1218,6 +1227,7 @@ def test_survey_df_widens_bounds(self): ci_width = h_result.ci_ub - h_result.ci_lb # Lower bound: normal-based CI width normal_width = 2 * 1.96 * h_result.original_se + assert np.isfinite(ci_width), "CI should be finite with M=1.0" assert ci_width > normal_width def test_no_survey_gives_none_df(self): diff --git a/tests/test_methodology_honest_did.py b/tests/test_methodology_honest_did.py new file mode 100644 index 00000000..efaf8fe7 --- /dev/null +++ b/tests/test_methodology_honest_did.py @@ -0,0 +1,478 @@ +""" +Methodology verification tests for Honest DiD (Rambachan & Roth, 2023). + +These tests verify the corrected implementation against the paper's +equations, known analytical cases, and expected mathematical properties. +""" + +import numpy as np +import pytest + +from diff_diff.honest_did import ( + HonestDiD, + _compute_flci, + _compute_optimal_flci, + _compute_pre_first_differences, + _construct_A_sd, + _construct_constraints_rm_component, + _construct_constraints_sd, + _cv_alpha, + _solve_bounds_lp, + _solve_rm_bounds_union, +) + + +# ============================================================================= +# TestDeltaSDConstraintMatrix +# ============================================================================= + + +class TestDeltaSDConstraintMatrix: + """Verify DeltaSD constraint matrix accounts for delta_0 = 0 boundary.""" + + def test_row_count(self): + """T+Tbar-1 rows, not T+Tbar-2 (accounts for delta_0 = 0).""" + for T, Tbar in [(2, 2), (3, 3), (4, 2), (1, 1), (3, 1), (1, 3)]: + A = _construct_A_sd(T, Tbar) + expected_rows = T + Tbar - 1 + assert A.shape == (expected_rows, T + Tbar), ( + f"T={T}, Tbar={Tbar}: expected {expected_rows} rows, got {A.shape[0]}" + ) + + def test_2pre_2post_hand_computed(self): + """Hand-computed matrix for 2 pre + 2 post periods.""" + # delta = [d_{-2}, d_{-1}, d_1, d_2] + A = _construct_A_sd(2, 2) + expected = np.array([ + [1, -2, 0, 0], # t=-1: d_{-2} - 2*d_{-1} + 0 + [0, 1, 1, 0], # t= 0: d_{-1} + d_1 (bridge) + [0, 0, -2, 1], # t= 1: 0 - 2*d_1 + d_2 + ]) + np.testing.assert_array_equal(A, expected) + + def test_bridge_constraint_present(self): + """The bridge constraint delta_{-1} + delta_1 is always present.""" + for T, Tbar in [(1, 1), (2, 2), (4, 3)]: + A = _construct_A_sd(T, Tbar) + # Find the bridge row: non-zero only at positions T-1 and T + bridge_found = False + for row in A: + if row[T - 1] != 0 and row[T] != 0: + # This should be [0, ..., 1, 1, ..., 0] + assert row[T - 1] == 1, f"Bridge row should have 1 at delta_{{-1}}" + assert row[T] == 1, f"Bridge row should have 1 at delta_1" + bridge_found = True + assert bridge_found, f"Bridge constraint not found for T={T}, Tbar={Tbar}" + + def test_constraints_span_all_periods(self): + """Constraints involve both pre and post periods (not pre-only).""" + A = _construct_A_sd(3, 3) + # Some rows should have non-zero entries in post-period columns + post_cols = A[:, 3:] # columns for delta_1, delta_2, delta_3 + assert np.any(post_cols != 0), "No constraints involve post-period deltas" + + +# ============================================================================= +# TestIdentifiedSetLP +# ============================================================================= + + +class TestIdentifiedSetLP: + """Verify identified set LP pins delta_pre = beta_pre.""" + + def test_m0_linear_extrapolation(self): + """M=0 with linear pre-trends gives finite point-identified bounds.""" + # Pre-trends: linear decline with slope -0.1 + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0]) + l_vec = np.array([1.0]) + + A, b = _construct_constraints_sd(3, 1, M=0.0) + lb, ub = _solve_bounds_lp(beta_pre, beta_post, l_vec, A, b, 3) + + # Linear extrapolation: slope = -0.1, so delta_1 = 0 - 0.1 = -0.1 + # theta = beta_post - delta_post = 2.0 - (-0.1) = 2.1 + assert np.isfinite(lb), "M=0 should give finite lower bound" + assert np.isfinite(ub), "M=0 should give finite upper bound" + np.testing.assert_allclose(lb, 2.1, atol=1e-6) + np.testing.assert_allclose(ub, 2.1, atol=1e-6) + + def test_bounds_widen_with_m(self): + """Identified set widens monotonically with M.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0]) + l_vec = np.array([1.0]) + + prev_width = 0 + for M in [0.0, 0.1, 0.5, 1.0]: + A, b = _construct_constraints_sd(3, 1, M=M) + lb, ub = _solve_bounds_lp(beta_pre, beta_post, l_vec, A, b, 3) + width = ub - lb + assert width >= prev_width - 1e-10, ( + f"Width should increase: M={M}, width={width}, prev={prev_width}" + ) + prev_width = width + + def test_three_period_analytical(self): + """Paper Section 2.3: three-period example (T=1, Tbar=1).""" + # delta = [d_{-1}, d_1], with delta_0 = 0 + # DeltaSD(M): |d_1 + d_{-1}| <= M (bridge constraint only) + # With d_{-1} = beta_{-1} pinned: + # d_1 in [-(beta_{-1} + M), -(beta_{-1} - M)] = [-beta_{-1} - M, -beta_{-1} + M] + # theta = beta_1 - d_1 + # lb = beta_1 - (-beta_{-1} + M) = beta_1 + beta_{-1} - M + # ub = beta_1 - (-beta_{-1} - M) = beta_1 + beta_{-1} + M + beta_pre = np.array([0.5]) + beta_post = np.array([3.0]) + + for M in [0.0, 0.2, 1.0]: + A, b = _construct_constraints_sd(1, 1, M=M) + lb, ub = _solve_bounds_lp(beta_pre, beta_post, np.array([1.0]), A, b, 1) + expected_lb = 3.0 + 0.5 - M + expected_ub = 3.0 + 0.5 + M + np.testing.assert_allclose(lb, expected_lb, atol=1e-6, + err_msg=f"M={M}: lb mismatch") + np.testing.assert_allclose(ub, expected_ub, atol=1e-6, + err_msg=f"M={M}: ub mismatch") + + +# ============================================================================= +# TestDeltaRMFirstDifferences +# ============================================================================= + + +class TestDeltaRMFirstDifferences: + """Verify DeltaRM constrains first differences, not levels.""" + + def test_pre_first_differences_computation(self): + """Pre-period first differences include delta_0=0 boundary.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + diffs = _compute_pre_first_differences(beta_pre) + + # Interior: |0.2-0.3|=0.1, |0.1-0.2|=0.1 + # Boundary: |0 - 0.1| = 0.1 + np.testing.assert_allclose(diffs, [0.1, 0.1, 0.1], atol=1e-10) + + def test_pre_first_differences_boundary(self): + """The boundary term |0 - beta_{-1}| is included.""" + beta_pre = np.array([0.0, 0.0, 0.5]) + diffs = _compute_pre_first_differences(beta_pre) + + # Interior: |0-0|=0, |0.5-0|=0.5 + # Boundary: |0 - 0.5| = 0.5 + np.testing.assert_allclose(diffs, [0.0, 0.5, 0.5], atol=1e-10) + + def test_rm_constraints_are_first_differences(self): + """RM constraint matrix constrains consecutive differences, not levels.""" + A, b = _construct_constraints_rm_component(2, 3, Mbar=1.0, max_pre_first_diff=0.1) + + # 3 post-period first diffs: |d_1|, |d_2-d_1|, |d_3-d_2| + # Each needs pos/neg constraint = 6 rows total + assert A.shape[0] == 6 + assert A.shape[1] == 5 # 2 pre + 3 post + + # First pair: d_1 <= 0.1 and -d_1 <= 0.1 + assert A[0, 2] == 1 # d_1 + assert A[1, 2] == -1 # -d_1 + + # Second pair: d_2 - d_1 <= 0.1 + assert A[2, 3] == 1 and A[2, 2] == -1 # d_2 - d_1 + assert A[3, 3] == -1 and A[3, 2] == 1 # -(d_2 - d_1) + + def test_mbar0_gives_point_estimate(self): + """Mbar=0: all post first diffs = 0, theta = l'beta_post.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0, 2.5]) + l_vec = np.array([0.5, 0.5]) + + lb, ub = _solve_rm_bounds_union(beta_pre, beta_post, l_vec, 3, Mbar=0.0) + + theta = np.dot(l_vec, beta_post) + np.testing.assert_allclose(lb, theta, atol=1e-6) + np.testing.assert_allclose(ub, theta, atol=1e-6) + + def test_rm_bounds_widen_with_mbar(self): + """Identified set widens monotonically with Mbar.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0, 2.5]) + l_vec = np.array([0.5, 0.5]) + + prev_width = 0 + for Mbar in [0.0, 0.5, 1.0, 2.0]: + lb, ub = _solve_rm_bounds_union(beta_pre, beta_post, l_vec, 3, Mbar) + width = ub - lb + assert width >= prev_width - 1e-10, f"Mbar={Mbar}: width decreased" + prev_width = width + + +# ============================================================================= +# TestOptimalFLCI +# ============================================================================= + + +class TestOptimalFLCI: + """Verify optimal FLCI properties.""" + + def test_cv_alpha_at_zero(self): + """cv_alpha(0, alpha) = z_{alpha/2} (standard normal quantile).""" + from scipy.stats import norm + np.testing.assert_allclose(_cv_alpha(0, 0.05), norm.ppf(0.975), atol=1e-4) + np.testing.assert_allclose(_cv_alpha(0, 0.01), norm.ppf(0.995), atol=1e-4) + + def test_cv_alpha_monotonic(self): + """cv_alpha(t) increases with |t| (more bias -> wider CI).""" + cvs = [_cv_alpha(t, 0.05) for t in [0, 0.5, 1.0, 2.0, 5.0]] + assert all(cvs[i] <= cvs[i + 1] + 1e-10 for i in range(len(cvs) - 1)) + + def test_optimal_flci_is_finite_and_valid(self): + """Optimal FLCI should produce finite CIs that cover identified set.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0]) + sigma = np.eye(4) * 0.01 + l_vec = np.array([1.0]) + + ci_lb_opt, ci_ub_opt = _compute_optimal_flci( + beta_pre, beta_post, sigma, l_vec, 3, 1, M=0.5, alpha=0.05 + ) + + # CI should be finite + assert np.isfinite(ci_lb_opt) and np.isfinite(ci_ub_opt) + # CI should cover the identified set + A, b = _construct_constraints_sd(3, 1, 0.5) + lb, ub = _solve_bounds_lp(beta_pre, beta_post, l_vec, A, b, 3) + assert ci_lb_opt <= lb, "CI lower should be <= identified set lower" + assert ci_ub_opt >= ub, "CI upper should be >= identified set upper" + + def test_m0_short_circuit(self): + """M=0 should use standard CI without optimization.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0]) + sigma = np.eye(4) * 0.01 + l_vec = np.array([1.0]) + + import time + t0 = time.time() + _compute_optimal_flci(beta_pre, beta_post, sigma, l_vec, 3, 1, M=0.0) + elapsed = time.time() - t0 + + assert elapsed < 0.5, f"M=0 should be fast, took {elapsed:.2f}s" + + def test_smoothness_flci_with_survey_df(self): + """Survey df should widen the smoothness FLCI (folded t vs folded normal).""" + beta_pre = np.array([0.1, 0.05]) + beta_post = np.array([2.0]) + sigma = np.eye(3) * 0.01 + + # Without df: uses folded normal + ci_lb_norm, ci_ub_norm = _compute_optimal_flci( + beta_pre, beta_post, sigma, np.array([1.0]), 2, 1, M=0.5 + ) + # With df=2: uses folded non-central t (wider critical values) + ci_lb_t, ci_ub_t = _compute_optimal_flci( + beta_pre, beta_post, sigma, np.array([1.0]), 2, 1, M=0.5, df=2 + ) + width_norm = ci_ub_norm - ci_lb_norm + width_t = ci_ub_t - ci_lb_t + assert width_t > width_norm, ( + f"Survey df=2 should widen CI: norm={width_norm:.4f}, t={width_t:.4f}" + ) + + def test_m0_se_includes_pre_period_variance(self): + """M=0 SE should account for pre-period variance, not just post.""" + # Use off-diagonal covariance to make pre-period SE matter + sigma = np.array([ + [0.04, 0.02, 0.01], # pre-1 has high variance + [0.02, 0.01, 0.005], + [0.01, 0.005, 0.01], + ]) + beta_pre = np.array([0.2, 0.1]) # linear pre-trend + beta_post = np.array([2.0]) + l_vec = np.array([1.0]) + + ci_lb, ci_ub = _compute_optimal_flci( + beta_pre, beta_post, sigma, l_vec, 2, 1, M=0.0 + ) + # CI should be finite and the width should reflect pre-period variance + assert np.isfinite(ci_lb) and np.isfinite(ci_ub), "M=0 CI should be finite" + width = ci_ub - ci_lb + + # Compare to post-only SE: sqrt(l'Sigma_post l) = sqrt(0.01) = 0.1 + post_only_width = 2 * 1.96 * np.sqrt(sigma[2, 2]) + assert width > post_only_width, ( + f"M=0 width ({width:.4f}) should exceed post-only ({post_only_width:.4f})" + ) + + def test_optimal_flci_width_increases_with_m_positive(self): + """Regression for P0: smoothness CI width must increase with M for M > 0.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0]) + sigma = np.eye(4) * 0.01 + + # Test monotonicity for M > 0 only. The M=0 path uses a different + # SE calculation (conservative, includes pre-period variance) which + # can produce a wider CI than small M > 0 where the optimizer is active. + widths = [] + for M in [0.1, 0.5, 1.0, 2.0]: + ci_lb, ci_ub = _compute_optimal_flci( + beta_pre, beta_post, sigma, np.array([1.0]), 3, 1, M=M + ) + widths.append(ci_ub - ci_lb) + + for i in range(len(widths) - 1): + assert widths[i + 1] >= widths[i] - 1e-4, ( + f"CI width must increase with M: M[{i}]={widths[i]:.4f}, " + f"M[{i+1}]={widths[i+1]:.4f}" + ) + + def test_optimal_flci_bias_nonzero_for_nonzero_m(self): + """Regression for P0: bias should be nonzero when M > 0.""" + from diff_diff.honest_did import _compute_worst_case_bias + + # T=3: 3 slopes (including boundary), sum(w)=1 for l=[1] + w = np.array([0.2, 0.3, 0.5]) + l_vec = np.array([1.0]) + + bias = _compute_worst_case_bias(w, l_vec, num_pre=3, num_post=1, M=0.5) + assert bias > 0, f"Bias should be nonzero for M>0, got {bias}" + + def test_three_period_m0_flci_center(self): + """T=1, Tbar=1, M=0: FLCI centered on beta_1 + beta_{-1}.""" + beta_pre = np.array([0.5]) + beta_post = np.array([3.0]) + sigma = np.eye(2) * 0.01 + + ci_lb, ci_ub = _compute_optimal_flci( + beta_pre, beta_post, sigma, np.array([1.0]), 1, 1, M=0.0 + ) + center = (ci_lb + ci_ub) / 2 + expected_center = 3.0 + 0.5 # beta_1 + beta_{-1} + np.testing.assert_allclose(center, expected_center, atol=1e-4, + err_msg="M=0 FLCI should be centered on beta_1 + beta_{-1}") + + def test_multi_post_m0_finite(self): + """Default l_vec with Tbar>1: M=0 gives finite CI.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0, 2.5]) + sigma = np.eye(5) * 0.01 + l_vec = np.array([0.5, 0.5]) # average of 2 post periods + + ci_lb, ci_ub = _compute_optimal_flci( + beta_pre, beta_post, sigma, l_vec, 3, 2, M=0.0 + ) + assert np.isfinite(ci_lb) and np.isfinite(ci_ub), ( + f"Multi-post M=0 should give finite CI, got [{ci_lb}, {ci_ub}]" + ) + + def test_multi_post_m_positive_finite(self): + """Default l_vec with Tbar>1: M>0 gives finite CI.""" + beta_pre = np.array([0.3, 0.2, 0.1]) + beta_post = np.array([2.0, 2.5]) + sigma = np.eye(5) * 0.01 + l_vec = np.array([0.5, 0.5]) + + ci_lb, ci_ub = _compute_optimal_flci( + beta_pre, beta_post, sigma, l_vec, 3, 2, M=0.5 + ) + assert np.isfinite(ci_lb) and np.isfinite(ci_ub), ( + f"Multi-post M=0.5 should give finite CI, got [{ci_lb}, {ci_ub}]" + ) + + def test_infeasible_lp_returns_nan(self): + """Regression for P1: infeasible LP should return NaN, not [-inf, inf].""" + # Non-linear pre-trends that are inconsistent with M=0 smoothness + beta_pre = np.array([1.0, 0.0, 1.0]) # quadratic, not linear + beta_post = np.array([2.0]) + A, b = _construct_constraints_sd(3, 1, M=0.0) + + lb, ub = _solve_bounds_lp(beta_pre, beta_post, np.array([1.0]), A, b, 3) + # M=0 with non-linear pre-trends: should be infeasible + assert np.isnan(lb) and np.isnan(ub), ( + f"Infeasible LP should return NaN, got [{lb}, {ub}]" + ) + + def test_infeasible_smoothness_fit_returns_nan_ci(self): + """Fit-level: infeasible smoothness restriction returns NaN CI.""" + from diff_diff.results import MultiPeriodDiDResults, PeriodEffect + + # Non-linear pre-trends: inconsistent with Delta^SD(M=0.01) + period_effects = { + 1: PeriodEffect(period=1, effect=1.0, se=0.1, t_stat=10.0, + p_value=0.0, conf_int=(0.8, 1.2)), + 2: PeriodEffect(period=2, effect=0.0, se=0.1, t_stat=0.0, + p_value=1.0, conf_int=(-0.2, 0.2)), + 3: PeriodEffect(period=3, effect=1.0, se=0.1, t_stat=10.0, + p_value=0.0, conf_int=(0.8, 1.2)), + 5: PeriodEffect(period=5, effect=2.0, se=0.1, t_stat=20.0, + p_value=0.0, conf_int=(1.8, 2.2)), + } + results = MultiPeriodDiDResults( + avg_att=2.0, avg_se=0.1, avg_t_stat=20.0, avg_p_value=0.0, + avg_conf_int=(1.8, 2.2), n_obs=500, n_treated=250, n_control=250, + period_effects=period_effects, pre_periods=[1, 2, 3], post_periods=[5], + vcov=np.eye(4) * 0.01, + interaction_indices={1: 0, 2: 1, 3: 2, 5: 3}, + ) + + honest = HonestDiD(method="smoothness", M=0.0) + r = honest.fit(results) + # Non-linear pre-trends should make M=0 infeasible + assert np.isnan(r.lb) and np.isnan(r.ub), f"Expected NaN bounds, got [{r.lb}, {r.ub}]" + assert np.isnan(r.ci_lb) and np.isnan(r.ci_ub), f"Expected NaN CI, got [{r.ci_lb}, {r.ci_ub}]" + # NaN CIs must NOT be classified as significant + assert not r.is_significant, "NaN CI should not be significant" + assert r.significance_stars == "", "NaN CI should have no significance stars" + assert "undefined" in repr(r).lower(), "NaN CI repr should indicate undefined" + + def test_smoothness_df_survey_zero_returns_nan(self): + """Smoothness with df_survey=0 should return NaN CI.""" + from diff_diff.honest_did import _compute_optimal_flci + + beta_pre = np.array([0.1, 0.05]) + beta_post = np.array([2.0]) + sigma = np.eye(3) * 0.01 + + # df=0 → NaN for all M + ci_lb, ci_ub = _compute_optimal_flci( + beta_pre, beta_post, sigma, np.array([1.0]), 2, 1, M=0.5, df=0 + ) + assert np.isnan(ci_lb) and np.isnan(ci_ub), "df=0 should give NaN CI" + + +# ============================================================================= +# TestBreakdownValueMethodology +# ============================================================================= + + +class TestBreakdownValueMethodology: + """Verify breakdown value properties.""" + + def test_breakdown_monotonicity(self): + """If significant at M=k, should be significant at all M < k.""" + from diff_diff.results import MultiPeriodDiDResults, PeriodEffect + + # Use a weak effect so breakdown is reachable at moderate M + period_effects = { + 1: PeriodEffect(period=1, effect=0.1, se=0.05, t_stat=2.0, + p_value=0.05, conf_int=(0.0, 0.2)), + 2: PeriodEffect(period=2, effect=0.05, se=0.05, t_stat=1.0, + p_value=0.32, conf_int=(-0.05, 0.15)), + 4: PeriodEffect(period=4, effect=0.15, se=0.05, t_stat=3.0, + p_value=0.003, conf_int=(0.05, 0.25)), + } + results = MultiPeriodDiDResults( + avg_att=0.15, avg_se=0.05, avg_t_stat=3.0, avg_p_value=0.003, + avg_conf_int=(0.05, 0.25), n_obs=500, n_treated=250, n_control=250, + period_effects=period_effects, pre_periods=[1, 2], post_periods=[4], + vcov=np.eye(3) * 0.0025, + interaction_indices={1: 0, 2: 1, 4: 2}, + ) + + honest = HonestDiD(method="smoothness") + # Check that CI at M=0 does not include zero + r0 = honest.fit(results, M=0.0) + assert r0.ci_lb > 0, "Should be significant at M=0" + + # At sufficiently large M, CI should include zero. + # The optimal FLCI is efficient, so need large M for a weak effect. + r_large = honest.fit(results, M=20.0) + assert r_large.ci_lb <= 0 <= r_large.ci_ub, "Should lose significance at large M"