*! did_bootstrap.mata - Bootstrap variance estimation for DID estimators *! *! Implements block bootstrap for variance-covariance estimation of the *! standard DID and sequential DID estimators. The VCOV matrix serves as *! the optimal weight matrix in the GMM framework for the double DID estimator. *! *! Main functions: *! did_boot_std() - Block bootstrap for standard DID design *! sample_panel() - Unit-level bootstrap for staggered adoption design *! bootstrap_se() - Extract standard errors from bootstrap VCOV *! bootstrap_ci() - Compute percentile confidence intervals version 16.0 mata: mata set matastrict on // ============================================================================ // Block Bootstrap for Variance Estimation // ============================================================================ // The variance-covariance matrix of (tau_DID, tau_sDID) is computed via // cluster bootstrap with the unbiased (B-1) denominator (see // compute_vcov() in did_utils.mata): // // Var(tau_DID, tau_sDID) = (1/(B-1)) * sum_b (tau^{(b)} - tau_bar)(tau^{(b)} - tau_bar)' // // This VCOV matrix is used as the optimal weight matrix W* in GMM estimation // of the double DID estimator: // // tau_ddid = argmin (tau - tau_DID, tau - tau_sDID)' W* (tau - tau_DID, tau - tau_sDID) // // Sampling is performed at the cluster level to preserve within-cluster // correlation structure required for valid inference. // ============================================================================ // ---------------------------------------------------------------------------- // Bootstrap Configuration Constants // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * _BOOT_FAIL_WARN_PCT() - Warning threshold for bootstrap failure rate *---------------------------------------------------------------------------*/ real scalar _BOOT_FAIL_WARN_PCT() { return(0.10) } /*--------------------------------------------------------------------------- * _BOOT_MIN_SUCCESS() - Minimum successful iterations for reliable variance *---------------------------------------------------------------------------*/ real scalar _BOOT_MIN_SUCCESS() { return(10) } // ---------------------------------------------------------------------------- // Data Structures // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * struct boot_result - Bootstrap result container * * Stores bootstrap estimates and metadata for variance computation. *---------------------------------------------------------------------------*/ struct boot_result { real matrix estimates // Bootstrap estimates (n_successful x 2*n_lead) pointer vector vcov // VCOV matrix for each lead (pointer to 2x2 matrix) real scalar n_successful // Number of successful iterations real scalar n_failed // Number of failed iterations } /*--------------------------------------------------------------------------- * struct boot_result_k - K-dimensional bootstrap result container * * Stores bootstrap estimates for K component estimators per lead. *---------------------------------------------------------------------------*/ struct boot_result_k { real matrix estimates // Bootstrap estimates (n_successful x kmax*n_lead) pointer vector vcov // VCOV matrix for each lead (pointer to kmax x kmax matrix) real scalar n_successful // Number of successful iterations real scalar n_failed // Number of failed iterations real scalar kmax // Number of components per lead } // ---------------------------------------------------------------------------- // Cluster Sampling // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * _boot_sample_clusters() - Sample cluster IDs with replacement * * Bootstrap sampling is performed at the cluster level to preserve * within-cluster correlation structure required for valid inference. * * Arguments: * id_cluster_vec : real colvector - unique cluster identifiers * n_clusters : real scalar - number of clusters * * Returns: * real colvector - sampled cluster IDs (same length as input) *---------------------------------------------------------------------------*/ real colvector _boot_sample_clusters(real colvector id_cluster_vec, real scalar n_clusters) { real colvector random_idx, id_boot random_idx = safe_sample_idx(n_clusters, n_clusters) id_boot = id_cluster_vec[random_idx] return(id_boot) } // ---------------------------------------------------------------------------- // Bootstrap Dataset Construction // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * _boot_create_dataset() - Construct bootstrap dataset from sampled clusters * * A new dataset is created by selecting observations belonging to sampled * clusters and reassigning sequential unit IDs. For panel data without an * explicit cluster variable, clustering is performed by id_unit. * * Arguments: * data : struct did_data - original data * id_boot : real colvector - sampled cluster IDs * id_cluster_vec : real colvector - unique cluster identifiers * * Returns: * struct did_data - bootstrap sample with reassigned unit IDs *---------------------------------------------------------------------------*/ struct did_data scalar _boot_create_dataset(struct did_data scalar data, real colvector id_boot, real colvector id_cluster_vec) { struct did_data scalar dat_boot real scalar n_clusters, j, n_obs, start_idx, has_cluster, n_j, use_id_unit real scalar next_unit_id, k, h, n_units_boot real colvector idx, cluster_col, unit_ids_j, unit_source n_clusters = rows(id_boot) has_cluster = (rows(data.cluster_var) > 0 && cols(data.cluster_var) > 0) // Determine clustering column: explicit cluster, id_unit, or row-level use_id_unit = (!has_cluster && data.is_panel && rows(data.id_unit) > 0) if (has_cluster) { cluster_col = data.cluster_var } else if (use_id_unit) { cluster_col = data.id_unit } else if (!data.is_panel) { cluster_col = (1::data.N) } else { cluster_col = J(0, 1, .) } // First pass: count total observations n_obs = 0 for (j = 1; j <= n_clusters; j++) { if (rows(cluster_col) == 0) { n_obs = n_obs + 1 } else { n_obs = n_obs + sum(cluster_col :== id_boot[j]) } } if (n_obs == 0) { dat_boot.N = 0 return(dat_boot) } // Initialize bootstrap data structure dat_boot.outcome = J(n_obs, 1, .) dat_boot.outcome_delta = J(n_obs, 1, .) dat_boot.treatment = J(n_obs, 1, .) dat_boot.id_unit = J(n_obs, 1, .) dat_boot.id_time = J(n_obs, 1, .) dat_boot.id_time_std = J(n_obs, 1, .) dat_boot.Gi = J(n_obs, 1, .) dat_boot.It = J(n_obs, 1, .) if (cols(data.covariates) > 0) { dat_boot.covariates = J(n_obs, cols(data.covariates), .) } else { dat_boot.covariates = J(0, 0, .) } if (has_cluster) { dat_boot.cluster_var = J(n_obs, 1, .) } else { dat_boot.cluster_var = J(0, 1, .) } // Second pass: fill data start_idx = 1 next_unit_id = 1 for (j = 1; j <= n_clusters; j++) { if (rows(cluster_col) == 0) { dat_boot.outcome[start_idx] = data.outcome[id_boot[j]] dat_boot.treatment[start_idx] = data.treatment[id_boot[j]] dat_boot.id_time[start_idx] = data.id_time[id_boot[j]] dat_boot.id_unit[start_idx] = next_unit_id next_unit_id = next_unit_id + 1 if (cols(data.covariates) > 0) { dat_boot.covariates[start_idx, .] = data.covariates[id_boot[j], .] } start_idx = start_idx + 1 } else { idx = selectindex(cluster_col :== id_boot[j]) n_j = rows(idx) if (n_j == 0) continue if (n_j == 1) { dat_boot.outcome[start_idx] = data.outcome[idx[1]] dat_boot.treatment[start_idx] = data.treatment[idx[1]] dat_boot.id_time[start_idx] = data.id_time[idx[1]] dat_boot.id_unit[start_idx] = next_unit_id next_unit_id = next_unit_id + 1 if (has_cluster) { dat_boot.cluster_var[start_idx] = data.cluster_var[idx[1]] } if (cols(data.covariates) > 0) { dat_boot.covariates[start_idx, .] = data.covariates[idx[1], .] } } else { dat_boot.outcome[start_idx::(start_idx+n_j-1)] = data.outcome[idx] dat_boot.treatment[start_idx::(start_idx+n_j-1)] = data.treatment[idx] dat_boot.id_time[start_idx::(start_idx+n_j-1)] = data.id_time[idx] if (has_cluster) { dat_boot.cluster_var[start_idx::(start_idx+n_j-1)] = data.cluster_var[idx] } if (cols(data.covariates) > 0) { dat_boot.covariates[start_idx::(start_idx+n_j-1), .] = data.covariates[idx, .] } if (data.is_panel && has_cluster) { // Preserve unit-level panel structure within each sampled cluster. unit_source = data.id_unit[idx] unit_ids_j = uniqrows(unit_source) for (k = 1; k <= rows(unit_ids_j); k++) { for (h = 1; h <= n_j; h++) { if (unit_source[h] == unit_ids_j[k]) { dat_boot.id_unit[start_idx + h - 1] = next_unit_id } } next_unit_id = next_unit_id + 1 } } else { dat_boot.id_unit[start_idx::(start_idx+n_j-1)] = J(n_j, 1, j) } } start_idx = start_idx + n_j } } // Set metadata dat_boot.N = n_obs if (data.is_panel) { n_units_boot = next_unit_id - 1 } else { n_units_boot = n_clusters } dat_boot.n_units = n_units_boot dat_boot.n_periods = data.n_periods dat_boot.is_panel = data.is_panel dat_boot.treat_year = data.treat_year return(dat_boot) } // ---------------------------------------------------------------------------- // Outcome Delta Computation // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * _compute_outcome_delta() - Outcome transformation for sequential DID * * The outcome transformation required for the sequential DID estimator * is computed. This estimator is consistent under the parallel trends-in- * trends assumption: * * DeltaY_{it} = Y_{it} - E[Y_{i,t-1} | G_i] * * where E[Y_{i,t-1} | G_i] is the group mean outcome in the previous period. * * The sequential DID estimator subtracts the pre-treatment DID from the * standard DID to remove bias when trends are not parallel but the change * in trends is the same across groups (parallel trends-in-trends). * * Arguments: * data : struct did_data - data with outcome, Gi, id_time_std * * Returns: * real colvector - transformed outcome (Y - lagged group mean) *---------------------------------------------------------------------------*/ real colvector _compute_outcome_delta(struct did_data scalar data) { real scalar n, i, g, t real colvector outcome_delta, Ymean string scalar key, lag_key transmorphic scalar group_sum_map, group_count_map real scalar sum_y, count_y n = data.N outcome_delta = J(n, 1, .) Ymean = J(n, 1, .) // Hash tables for O(N) complexity group_sum_map = asarray_create("string", 1) group_count_map = asarray_create("string", 1) // First pass: accumulate sums and counts for each (Gi, id_time_std) group for (i = 1; i <= n; i++) { g = data.Gi[i] t = data.id_time_std[i] if (missing(g) || missing(t)) continue key = strofreal(g) + "_" + strofreal(t) if (missing(data.outcome[i])) { continue } if (asarray_contains(group_sum_map, key)) { asarray(group_sum_map, key, asarray(group_sum_map, key) + data.outcome[i]) asarray(group_count_map, key, asarray(group_count_map, key) + 1) } else { asarray(group_sum_map, key, data.outcome[i]) asarray(group_count_map, key, 1) } } // Second pass: assign lag group means using O(1) lookups for (i = 1; i <= n; i++) { g = data.Gi[i] t = data.id_time_std[i] if (missing(g) || missing(t)) continue lag_key = strofreal(g) + "_" + strofreal(t - 1) if (asarray_contains(group_sum_map, lag_key)) { sum_y = asarray(group_sum_map, lag_key) count_y = asarray(group_count_map, lag_key) if (count_y > 0) { Ymean[i] = sum_y / count_y } } } outcome_delta = data.outcome - Ymean return(outcome_delta) } // ---------------------------------------------------------------------------- // Bootstrap Data Preparation // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * _normalize_time() - Normalize time index to consecutive integers * * Time values are compressed to consecutive integers 1, 2, 3, ..., preserving * the temporal ordering. This is required for bootstrap samples that may have * non-consecutive time periods due to resampling. * * Arguments: * time_values : real colvector - original time values * * Returns: * real colvector - normalized time values (1, 2, 3, ...) *---------------------------------------------------------------------------*/ real colvector _normalize_time(real colvector time_values) { real colvector unique_times, result real scalar n, n_unique, i transmorphic scalar time_map n = rows(time_values) if (n == 0) return(J(0, 1, .)) unique_times = uniqrows(time_values) n_unique = rows(unique_times) result = J(n, 1, .) // Hash map for O(1) lookup time_map = asarray_create("real", 1) for (i = 1; i <= n_unique; i++) { asarray(time_map, unique_times[i], i) } for (i = 1; i <= n; i++) { if (!missing(time_values[i])) { result[i] = asarray(time_map, time_values[i]) } } return(result) } /*--------------------------------------------------------------------------- * _boot_panel_prep() - Re-run panel data preparation for bootstrap sample * * Derived variables (Gi, It, id_time_std, outcome_delta) are recomputed * for a bootstrap sample. Panel bootstrap draws must preserve the original * observed-sample treatment calendar; compressing time inside a draw can * fabricate valid-looking {-1,0} windows when the true t-1 period is absent. * * DID identification requires both treated and control units, and * observations in both pre- and post-treatment periods. An empty Gi * vector signals identification failure. * * Arguments: * dat_boot : struct did_data - bootstrap sample * * Returns: * struct did_data - bootstrap sample with recomputed derived variables *---------------------------------------------------------------------------*/ struct did_data scalar _boot_panel_prep(struct did_data scalar dat_boot) { real scalar n, i, max_treat real colvector unit_ids, unit_idx n = dat_boot.N if (n == 0) return(dat_boot) // Compute Gi = max(treatment) by id_unit unit_ids = uniqrows(dat_boot.id_unit) dat_boot.Gi = J(n, 1, .) for (i = 1; i <= rows(unit_ids); i++) { unit_idx = selectindex(dat_boot.id_unit :== unit_ids[i]) max_treat = max(dat_boot.treatment[unit_idx]) dat_boot.Gi[unit_idx] = J(rows(unit_idx), 1, max_treat) } // Preserve the original treatment calendar. Re-normalizing a panel // bootstrap draw can turn a missing t-1 period into a pseudo-valid {-1,0} // window and pollute valid-draw counts plus downstream VCOV/CI paths. if (missing(dat_boot.treat_year)) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } // Validate identification: require both treated and control units if (sum(dat_boot.Gi :== 0) == 0) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } if (sum(dat_boot.Gi :== 1) == 0) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } // Compute standardized time index dat_boot.id_time_std = dat_boot.id_time :- dat_boot.treat_year // Compute post-treatment indicator dat_boot.It = (dat_boot.id_time :>= dat_boot.treat_year) // Validate identification: require both pre and post periods if (sum(dat_boot.It :== 0) == 0 || sum(dat_boot.It :== 1) == 0) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } // Compute outcome_delta for sequential DID dat_boot.outcome_delta = _compute_outcome_delta(dat_boot) return(dat_boot) } /*--------------------------------------------------------------------------- * _boot_rcs_prep() - Re-run RCS data preparation for bootstrap sample * * For repeated cross-sectional data, Gi and It are input variables (not * computed from treatment timing). This function copies Gi and It from * original data, preserves the observed-sample calendar-time coding, and * recomputes id_time_std and outcome_delta for the sequential DID estimator. * * Arguments: * dat_boot : struct did_data - bootstrap sample * data_orig : struct did_data - original data (source of Gi, It) * id_boot : real colvector - sampled indices * * Returns: * struct did_data - bootstrap sample with recomputed time-relative variables *---------------------------------------------------------------------------*/ struct did_data scalar _boot_rcs_prep(struct did_data scalar dat_boot, struct did_data scalar data_orig, real colvector id_boot) { real scalar n, j, start_idx, n_j, has_cluster real scalar n_orig, boot_idx real colvector idx, cluster_source n = dat_boot.N if (n == 0) return(dat_boot) has_cluster = (rows(data_orig.cluster_var) > 0 && cols(data_orig.cluster_var) > 0) n_orig = rows(data_orig.Gi) if (has_cluster) { cluster_source = data_orig.cluster_var } else { cluster_source = (1::n_orig) } // Copy Gi and It from original data dat_boot.Gi = J(n, 1, .) dat_boot.It = J(n, 1, .) start_idx = 1 for (j = 1; j <= rows(id_boot); j++) { if (!has_cluster) { boot_idx = id_boot[j] if (boot_idx < 1 || boot_idx > n_orig) { continue } if (start_idx > n) { break } dat_boot.Gi[start_idx] = data_orig.Gi[boot_idx] dat_boot.It[start_idx] = data_orig.It[boot_idx] start_idx = start_idx + 1 } else { idx = selectindex(cluster_source :== id_boot[j]) n_j = rows(idx) if (n_j == 0) continue if (start_idx + n_j - 1 > n) { n_j = n - start_idx + 1 if (n_j <= 0) break } if (n_j == 1) { dat_boot.Gi[start_idx] = data_orig.Gi[idx[1]] dat_boot.It[start_idx] = data_orig.It[idx[1]] } else { dat_boot.Gi[start_idx::(start_idx+n_j-1)] = data_orig.Gi[idx] dat_boot.It[start_idx::(start_idx+n_j-1)] = data_orig.It[idx] } start_idx = start_idx + n_j } } dat_boot.Gi = normalize_binary01(dat_boot.Gi, 1e-6) dat_boot.It = normalize_binary01(dat_boot.It, 1e-6) // Verify index consistency if (start_idx - 1 != n) { errprintf("{txt}Warning: _boot_rcs_prep index mismatch: expected %g rows, filled %g\n", n, start_idx - 1) dat_boot.Gi = J(0, 1, .) return(dat_boot) } // Preserve the original calendar-time support. Compressing time inside a // bootstrap draw can fabricate valid-looking placebo / lead windows that // do not exist under the paper's original treatment calendar. if (missing(data_orig.treat_year)) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } dat_boot.treat_year = data_orig.treat_year // Validate identification if (sum(dat_boot.Gi :== 0) == 0) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } if (sum(dat_boot.Gi :== 1) == 0) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } if (sum(dat_boot.It :== 0) == 0) { dat_boot.Gi = J(0, 1, .) return(dat_boot) } // Recompute id_time_std and outcome_delta dat_boot.id_time_std = dat_boot.id_time :- data_orig.treat_year dat_boot.outcome_delta = _compute_outcome_delta(dat_boot) return(dat_boot) } // ---------------------------------------------------------------------------- // Main Bootstrap Function // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * did_boot_std() - Block bootstrap for standard DID design * * The variance-covariance matrix Sigma of (tau_DID, tau_sDID) is estimated * via cluster bootstrap. This VCOV matrix is used to construct the optimal * weight matrix in the GMM framework for the double DID estimator: * * tau_dDID = argmin (tau - tau_DID, tau - tau_sDID)' W (tau - tau_DID, tau - tau_sDID) * * where W = Sigma^{-1} is the precision matrix (optimal GMM weight). * The double DID estimator combines the standard DID and sequential DID * to achieve efficiency under parallel trends and robustness under the * weaker parallel trends-in-trends assumption. * * Sampling is performed at the cluster level with replacement to preserve * within-cluster correlation structure for valid inference. * * Arguments: * data : struct did_data - prepared data * lead : real rowvector - post-treatment period indices * n_boot : real scalar - number of bootstrap iterations * seed : real scalar - random seed (optional) * * Returns: * struct boot_result - bootstrap estimates and VCOV matrices *---------------------------------------------------------------------------*/ struct boot_result scalar did_boot_std(struct did_data scalar data, real rowvector lead, real scalar n_boot, | real scalar seed) { struct boot_result scalar result struct did_data scalar dat_boot real colvector id_cluster_vec, id_boot, valid_idx real scalar n_clusters, n_lead, b, l, n_successful, n_failed real matrix boot_est real rowvector est_b if (args() >= 4 && !missing(seed)) { rseed(seed) } n_lead = cols(lead) // Handle n_boot=0 case if (n_boot == 0) { printf("{txt}Note: n_boot=0, no bootstrap inference available\n") result.estimates = J(0, 2 * n_lead, .) result.vcov = J(n_lead, 1, NULL) result.n_successful = 0 result.n_failed = 0 return(result) } // Determine cluster structure for bootstrap sampling if (rows(data.cluster_var) > 0 && cols(data.cluster_var) > 0) { real colvector valid_cluster_mask, valid_clusters real scalar n_missing_clusters valid_cluster_mask = (data.cluster_var :< .) valid_clusters = select(data.cluster_var, valid_cluster_mask) n_missing_clusters = rows(data.cluster_var) - rows(valid_clusters) if (n_missing_clusters > 0) { printf("{txt}Warning: cluster variable contains %g missing values (excluded from bootstrap)\n", n_missing_clusters) } if (rows(valid_clusters) > 0) { id_cluster_vec = unique_in_order(valid_clusters) } else { errprintf("Error: cluster variable contains only missing values\n") result.estimates = J(0, 2 * n_lead, .) result.vcov = J(n_lead, 1, NULL) result.n_successful = 0 result.n_failed = n_boot return(result) } } else if (data.is_panel && rows(data.id_unit) > 0) { id_cluster_vec = unique_in_order(data.id_unit) } else { id_cluster_vec = (1::data.N) } n_clusters = rows(id_cluster_vec) // Initialize bootstrap estimates matrix boot_est = J(n_boot, 2 * n_lead, .) valid_idx = J(n_boot, 1, 0) // Bootstrap loop for (b = 1; b <= n_boot; b++) { // Sample clusters with replacement id_boot = _boot_sample_clusters(id_cluster_vec, n_clusters) // Construct bootstrap dataset dat_boot = _boot_create_dataset(data, id_boot, id_cluster_vec) if (dat_boot.N == 0) { continue } // Re-run data preparation if (data.is_panel) { dat_boot = _boot_panel_prep(dat_boot) } else { dat_boot = _boot_rcs_prep(dat_boot, data, id_boot) } if (rows(dat_boot.Gi) == 0) { continue } if (missing(dat_boot.Gi[1])) { continue } // Compute estimates for each lead est_b = J(1, 2 * n_lead, .) for (l = 1; l <= n_lead; l++) { est_b[1, (2*l-1)..(2*l)] = did_fit( dat_boot.outcome, dat_boot.outcome_delta, dat_boot.Gi, dat_boot.It, dat_boot.id_unit, dat_boot.covariates, dat_boot.id_time_std, lead[l], dat_boot.is_panel ) } // Check if at least one lead has valid estimates real scalar has_valid_est, ll_check has_valid_est = 0 for (ll_check = 1; ll_check <= n_lead; ll_check++) { if (!missing(est_b[1, 2*ll_check-1]) || !missing(est_b[1, 2*ll_check])) { has_valid_est = 1 break } } if (has_valid_est) { boot_est[b, .] = est_b valid_idx[b] = 1 } } // Summarize results n_successful = sum(valid_idx) n_failed = n_boot - n_successful if (n_successful > 0) { boot_est = select(boot_est, valid_idx) } else { boot_est = J(0, 2 * n_lead, .) } // Warn if many iterations failed if (n_failed > _BOOT_FAIL_WARN_PCT() * n_boot) { real scalar pct_failed pct_failed = 100 * n_failed / n_boot errprintf("Warning: %g / %g bootstrap iterations failed (%g percent)\n", n_failed, n_boot, pct_failed) } if (n_successful < _BOOT_MIN_SUCCESS() && n_successful > 0) { if (n_failed > 0) { errprintf("Warning: Only %g bootstrap iterations succeeded, results may be unreliable\n", n_successful) } else { errprintf("Warning: All %g bootstrap iterations succeeded, but nboot() is below the recommended minimum of %g\n", n_successful, _BOOT_MIN_SUCCESS()) } } // Compute VCOV for each lead result.estimates = boot_est result.vcov = J(n_lead, 1, NULL) result.n_successful = n_successful result.n_failed = n_failed if (n_successful >= 2) { for (l = 1; l <= n_lead; l++) { real matrix boot_est_lead real scalar n_valid_did, n_valid_sdid real matrix vcov_joint boot_est_lead = boot_est[., (2*l-1)..(2*l)] n_valid_did = rows(select(boot_est_lead[., 1], boot_est_lead[., 1] :< .)) n_valid_sdid = rows(select(boot_est_lead[., 2], boot_est_lead[., 2] :< .)) if (n_valid_did >= 2 && n_valid_sdid >= 2) { vcov_joint = compute_vcov_joint_valid(boot_est_lead) // Allocate a distinct matrix per lead so multi-lead runs do not // alias the last loop-local vcov_joint value across pointers. result.vcov[l] = &(J(2, 2, .)) *result.vcov[l] = vcov_joint } else { result.vcov[l] = &(J(2, 2, .)) } } } return(result) } // ---------------------------------------------------------------------------- // Standard Error and Confidence Interval Functions // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * bootstrap_se() - Extract standard errors from bootstrap VCOV * * Standard errors are computed as square roots of diagonal elements: * SE(tau_DID) = sqrt(Var[1,1]) * SE(tau_sDID) = sqrt(Var[2,2]) * * Arguments: * vcov : real matrix (2 x 2) - variance-covariance matrix * * Returns: * real rowvector (1 x 2) - (SE_DID, SE_sDID) *---------------------------------------------------------------------------*/ real rowvector bootstrap_se(real matrix vcov) { real rowvector se if (rows(vcov) != 2 || cols(vcov) != 2) { return((., .)) } if (missing(vcov[1,1]) || missing(vcov[2,2])) { return((., .)) } se = (sqrt(vcov[1,1]), sqrt(vcov[2,2])) return(se) } /*--------------------------------------------------------------------------- * bootstrap_ci() - Compute percentile bootstrap confidence intervals * * Confidence intervals are computed using the percentile method: * ci_low = quantile(boot_est, alpha/2) * ci_high = quantile(boot_est, 1 - alpha/2) * where alpha = 1 - level/100. * * Arguments: * boot_est : real colvector - bootstrap estimates for one parameter * level : real scalar - confidence level (e.g., 95) * * Returns: * real rowvector (1 x 2) - (ci_low, ci_high) *---------------------------------------------------------------------------*/ real rowvector bootstrap_ci(real colvector boot_est, real scalar level) { real scalar alpha, p_low, p_high real rowvector ci if (level <= 0 || level >= 100) { return((., .)) } alpha = 1 - level / 100 p_low = alpha / 2 p_high = 1 - alpha / 2 ci = (quantile_sorted(boot_est, p_low), quantile_sorted(boot_est, p_high)) return(ci) } // ============================================================================ // Panel Bootstrap for Staggered Adoption Design // ============================================================================ // Unit-level bootstrap is implemented for the staggered adoption (SA) design // where treatment timing varies across units. Unlike the standard DID bootstrap // which samples clusters, the SA design requires sampling entire units with // replacement to preserve the time series structure needed for: // // 1. Reconstructing the treatment timing matrix (Gmat) encoding adoption times // 2. Computing period-specific SA-ATT estimates: tau_DID(t), tau_sDID(t) // 3. Aggregating via time weights: tau_bar = sum_t pi_t * tau(t) // // The SA double DID extends the basic double DID framework by applying the // GMM combination of DID and sequential DID estimators at each treatment time, // then aggregating across time periods using appropriate weights. // ============================================================================ /*--------------------------------------------------------------------------- * sample_panel() - Unit-level bootstrap for staggered adoption design * * Entire units are sampled with replacement, preserving the time series * structure required for period-specific SA-ATT estimation. Each sampled * unit is assigned a new sequential ID to handle duplicate units. * * Derived fields (Gi, It, id_time_std, outcome_delta) are NOT populated here. * These must be recomputed downstream after the treatment timing matrix (Gmat) * is rebuilt from the bootstrap sample. * * Arguments: * data : struct did_data - panel data structure * * Returns: * struct did_data - bootstrap sample with reassigned unit IDs *---------------------------------------------------------------------------*/ struct did_data scalar sample_panel(struct did_data scalar data) { struct did_data scalar boot_data real colvector block_vec, block_boot, idx, new_id_unit real matrix boot_outcome, boot_treatment, boot_time, boot_covariates real colvector boot_cluster real scalar n_blocks, i, j, n_obs_total, n_obs_i, n_units_boot real scalar has_covariates, has_cluster real colvector unit_ids_block, unique_units_block, mapped_units_block real scalar u, r transmorphic scalar unit_map // Validate panel data structure if (!data.is_panel) { errprintf("sample_panel(): Requires panel data (is_panel=1)\n") errprintf(" For RCS data, use cluster bootstrap instead\n") boot_data.outcome = J(0, 1, .) boot_data.treatment = J(0, 1, .) boot_data.id_unit = J(0, 1, .) boot_data.id_time = J(0, 1, .) boot_data.covariates = J(0, 0, .) boot_data.cluster_var = J(0, 1, .) boot_data.Gi = J(0, 1, .) boot_data.It = J(0, 1, .) boot_data.id_time_std = J(0, 1, .) boot_data.outcome_delta = J(0, 1, .) boot_data.N = 0 boot_data.n_units = 0 boot_data.n_periods = 0 boot_data.treat_year = . boot_data.is_panel = 0 return(boot_data) } // Handle empty dataset if (data.N == 0 || data.n_units == 0) { boot_data.outcome = J(0, 1, .) boot_data.treatment = J(0, 1, .) boot_data.id_unit = J(0, 1, .) boot_data.id_time = J(0, 1, .) boot_data.covariates = J(0, 0, .) boot_data.cluster_var = J(0, 1, .) boot_data.Gi = J(0, 1, .) boot_data.It = J(0, 1, .) boot_data.id_time_std = J(0, 1, .) boot_data.outcome_delta = J(0, 1, .) boot_data.N = 0 boot_data.n_units = 0 boot_data.n_periods = data.n_periods boot_data.treat_year = . boot_data.is_panel = data.is_panel return(boot_data) } has_covariates = (cols(data.covariates) > 0 && rows(data.covariates) == data.N) has_cluster = (rows(data.cluster_var) == data.N && cols(data.cluster_var) > 0) // Sample cluster IDs with replacement when cluster_var is available. // Otherwise fall back to the unit-level bootstrap. if (has_cluster) { block_vec = unique_in_order(data.cluster_var) } else { block_vec = unique_in_order(data.id_unit) } n_blocks = rows(block_vec) block_boot = block_vec[safe_sample_idx(n_blocks, n_blocks)] // Cache selectindex() results for efficiency transmorphic idx_cache idx_cache = asarray_create("real", 1) asarray_notfound(idx_cache, J(0, 1, .)) // First pass: count total observations in the sampled blocks n_obs_total = 0 for (i = 1; i <= n_blocks; i++) { if (has_cluster) { idx = selectindex(data.cluster_var :== block_boot[i]) } else { idx = selectindex(data.id_unit :== block_boot[i]) } asarray(idx_cache, i, idx) n_obs_total = n_obs_total + rows(idx) } // Allocate result matrices boot_outcome = J(n_obs_total, 1, .) boot_treatment = J(n_obs_total, 1, .) boot_time = J(n_obs_total, 1, .) new_id_unit = J(n_obs_total, 1, .) if (has_covariates) { boot_covariates = J(n_obs_total, cols(data.covariates), .) } else { boot_covariates = J(0, 0, .) } if (has_cluster) { boot_cluster = J(n_obs_total, 1, .) } else { boot_cluster = J(0, 1, .) } // Second pass: fill data using cached indices j = 1 n_units_boot = 0 for (i = 1; i <= n_blocks; i++) { idx = asarray(idx_cache, i) n_obs_i = rows(idx) if (n_obs_i == 0) continue // Renormalize unit identifiers within each sampled block so that the // bootstrap sample remains a valid balanced panel on consecutive ids. unit_ids_block = data.id_unit[idx] unique_units_block = uniqrows(unit_ids_block) unit_map = asarray_create("real") mapped_units_block = J(n_obs_i, 1, .) for (u = 1; u <= rows(unique_units_block); u++) { n_units_boot++ asarray(unit_map, unique_units_block[u], n_units_boot) } for (r = 1; r <= n_obs_i; r++) { mapped_units_block[r] = asarray(unit_map, unit_ids_block[r]) } if (n_obs_i == 1) { boot_outcome[j] = data.outcome[idx[1]] boot_treatment[j] = data.treatment[idx[1]] boot_time[j] = data.id_time[idx[1]] new_id_unit[j] = mapped_units_block[1] if (has_covariates) { boot_covariates[j, .] = data.covariates[idx[1], .] } if (has_cluster) { boot_cluster[j] = data.cluster_var[idx[1]] } } else { boot_outcome[j::(j+n_obs_i-1)] = data.outcome[idx] boot_treatment[j::(j+n_obs_i-1)] = data.treatment[idx] boot_time[j::(j+n_obs_i-1)] = data.id_time[idx] new_id_unit[j::(j+n_obs_i-1)] = mapped_units_block if (has_covariates) { boot_covariates[j::(j+n_obs_i-1), .] = data.covariates[idx, .] } if (has_cluster) { boot_cluster[j::(j+n_obs_i-1)] = data.cluster_var[idx] } } j = j + n_obs_i } // Populate result structure boot_data.outcome = boot_outcome boot_data.treatment = boot_treatment boot_data.id_time = boot_time boot_data.id_unit = new_id_unit boot_data.covariates = boot_covariates boot_data.cluster_var = boot_cluster // Update metadata boot_data.N = n_obs_total boot_data.n_units = n_units_boot boot_data.n_periods = data.n_periods boot_data.is_panel = data.is_panel // Derived fields are NOT populated - must be computed downstream boot_data.Gi = J(0, 1, .) boot_data.It = J(0, 1, .) boot_data.id_time_std = J(0, 1, .) boot_data.outcome_delta = J(0, 1, .) boot_data.treat_year = . return(boot_data) } // ============================================================================ // K-DIMENSIONAL BOOTSTRAP FOR GENERALIZED K-DID // ============================================================================ /*--------------------------------------------------------------------------- * did_boot_std_k() - Block bootstrap for generalized K-DID * * K-dimensional extension of did_boot_std(). Each bootstrap replication * produces a K-dimensional component vector (τ̂_1, ..., τ̂_K) per lead, * stored in an n_success × (kmax * n_lead) matrix. * * Arguments: * data : struct did_data - prepared data * lead : real rowvector - post-treatment lead values * n_boot : real scalar - number of bootstrap iterations * kmax : real scalar - max number of K-DID components * seed : real scalar (optional) - random seed * * Returns: * struct boot_result_k - K-dimensional bootstrap estimates and VCOV *---------------------------------------------------------------------------*/ struct boot_result_k scalar did_boot_std_k(struct did_data scalar data, real rowvector lead, real scalar n_boot, real scalar kmax, | real scalar seed) { struct boot_result_k scalar result struct did_data scalar dat_boot real colvector id_cluster_vec, id_boot, valid_idx real scalar n_clusters, n_lead, b, l, k_comp, n_successful, n_failed real matrix boot_est real rowvector est_b, comp_k if (args() >= 5 && !missing(seed)) { rseed(seed) } n_lead = cols(lead) result.kmax = kmax // Handle n_boot=0 case if (n_boot == 0) { printf("{txt}Note: n_boot=0, no bootstrap inference available\n") result.estimates = J(0, kmax * n_lead, .) result.vcov = J(n_lead, 1, NULL) result.n_successful = 0 result.n_failed = 0 return(result) } // Determine cluster structure (same as did_boot_std) if (rows(data.cluster_var) > 0 && cols(data.cluster_var) > 0) { real colvector valid_cluster_mask, valid_clusters real scalar n_missing_clusters valid_cluster_mask = (data.cluster_var :< .) valid_clusters = select(data.cluster_var, valid_cluster_mask) n_missing_clusters = rows(data.cluster_var) - rows(valid_clusters) if (n_missing_clusters > 0) { printf("{txt}Warning: cluster variable contains %g missing values (excluded from bootstrap)\n", n_missing_clusters) } if (rows(valid_clusters) > 0) { id_cluster_vec = unique_in_order(valid_clusters) } else { errprintf("Error: cluster variable contains only missing values\n") result.estimates = J(0, kmax * n_lead, .) result.vcov = J(n_lead, 1, NULL) result.n_successful = 0 result.n_failed = n_boot return(result) } } else if (data.is_panel && rows(data.id_unit) > 0) { id_cluster_vec = unique_in_order(data.id_unit) } else { id_cluster_vec = (1::data.N) } n_clusters = rows(id_cluster_vec) // Initialize bootstrap estimates matrix: kmax columns per lead boot_est = J(n_boot, kmax * n_lead, .) valid_idx = J(n_boot, 1, 0) // Bootstrap loop for (b = 1; b <= n_boot; b++) { id_boot = _boot_sample_clusters(id_cluster_vec, n_clusters) dat_boot = _boot_create_dataset(data, id_boot, id_cluster_vec) if (dat_boot.N == 0) { continue } if (data.is_panel) { dat_boot = _boot_panel_prep(dat_boot) } else { dat_boot = _boot_rcs_prep(dat_boot, data, id_boot) } if (rows(dat_boot.Gi) == 0 || missing(dat_boot.Gi[1])) { continue } // Compute K-dimensional component estimates for each lead est_b = J(1, kmax * n_lead, .) for (l = 1; l <= n_lead; l++) { comp_k = did_fit_k( dat_boot.outcome, dat_boot.Gi, dat_boot.It, dat_boot.id_unit, dat_boot.covariates, dat_boot.id_time_std, lead[l], kmax, dat_boot.is_panel ) est_b[1, (kmax*(l-1)+1)..(kmax*l)] = comp_k } // Check if at least one lead has a valid estimate real scalar has_valid_est_k, ll_k, kk has_valid_est_k = 0 for (ll_k = 1; ll_k <= n_lead; ll_k++) { for (kk = 1; kk <= kmax; kk++) { if (!missing(est_b[1, kmax*(ll_k-1)+kk])) { has_valid_est_k = 1 break } } if (has_valid_est_k) break } if (has_valid_est_k) { boot_est[b, .] = est_b valid_idx[b] = 1 } } // Summarize results n_successful = sum(valid_idx) n_failed = n_boot - n_successful if (n_successful > 0) { boot_est = select(boot_est, valid_idx) } else { boot_est = J(0, kmax * n_lead, .) } if (n_failed > _BOOT_FAIL_WARN_PCT() * n_boot) { real scalar pct_failed_k pct_failed_k = 100 * n_failed / n_boot errprintf("Warning: %g / %g bootstrap iterations failed (%g percent)\n", n_failed, n_boot, pct_failed_k) } if (n_successful < _BOOT_MIN_SUCCESS() && n_successful > 0) { errprintf("Warning: Only %g bootstrap iterations succeeded, results may be unreliable\n", n_successful) } // Compute VCOV for each lead (K × K) result.estimates = boot_est result.vcov = J(n_lead, 1, NULL) result.n_successful = n_successful result.n_failed = n_failed if (n_successful >= 2) { for (l = 1; l <= n_lead; l++) { real matrix boot_est_lead_k, vcov_k real scalar col_start, col_end col_start = kmax * (l - 1) + 1 col_end = kmax * l boot_est_lead_k = boot_est[., col_start..col_end] // Compute K×K VCOV using joint-valid observations vcov_k = compute_vcov_joint_valid(boot_est_lead_k) result.vcov[l] = &(J(kmax, kmax, .)) *result.vcov[l] = vcov_k } } return(result) } // ---------------------------------------------------------------------------- // Parallel Bootstrap: Chunked Execution Functions // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * did_boot_std_chunk() - Chunked block bootstrap for standard DID * * Executes bootstrap iterations [b_start, b_end] using the same algorithm * as did_boot_std(). The optional seed_chunk is set before execution to * ensure reproducibility of this chunk across parallel runs. * * Arguments: * data : struct did_data - prepared data (same as did_boot_std) * lead : real rowvector - post-treatment lead values * b_start : real scalar - first iteration index of this chunk (>=1) * b_end : real scalar - last iteration index of this chunk (>=b_start) * seed_chunk : real scalar - (optional) RNG seed for this chunk; * if missing, no seed is set * * Returns: * struct boot_result - same format as did_boot_std() *---------------------------------------------------------------------------*/ struct boot_result scalar did_boot_std_chunk( struct did_data scalar data, real rowvector lead, real scalar b_start, real scalar b_end, | real scalar seed_chunk) { struct boot_result scalar result real scalar n_iter // Set RNG seed for this chunk if specified if (args() >= 5 && !missing(seed_chunk)) { rseed(seed_chunk) } // Compute number of iterations for this chunk n_iter = b_end - b_start + 1 if (n_iter <= 0) { result.estimates = J(0, 2 * cols(lead), .) result.vcov = J(cols(lead), 1, NULL) result.n_successful = 0 result.n_failed = 0 return(result) } // Delegate to did_boot_std() with n_iter iterations. // Seed is already set above; do NOT pass seed to did_boot_std() // to avoid overwriting the chunk seed. result = did_boot_std(data, lead, n_iter) return(result) } /*--------------------------------------------------------------------------- * did_boot_sa_chunk() - Chunked unit bootstrap for staggered adoption design * * Executes SA bootstrap iterations [b_start, b_end]. Each iteration resamples * entire units via sample_panel(), re-estimates the time-weighted SA-DID and * SA-sDID for all leads via sa_double_did(), and accumulates raw (DID[l], sDID[l]) * pairs. The resulting matrix preserves the bootstrap covariance structure * required by sa_to_ddid_matrix() for GMM-optimal weighting. * * Arguments: * data : struct did_data - panel data (is_panel must equal 1) * option : struct did_option - estimation options (lead, thres, level, quiet) * b_start : real scalar - first iteration index of this chunk (>=1) * b_end : real scalar - last iteration index of this chunk (>=b_start) * seed_chunk : real scalar - (optional) RNG seed for this chunk * * Returns: * struct boot_result * .estimates : n_successful x (2*n_lead) matrix; * columns 2l-1, 2l = (DID[l], sDID[l]) for lead index l * .n_successful : count of iterations producing at least one valid estimate * .n_failed : count of failed iterations * .vcov : J(n_lead, 1, NULL) -- filled by sa_to_ddid_matrix() *---------------------------------------------------------------------------*/ struct boot_result scalar did_boot_sa_chunk( struct did_data scalar data, struct did_option scalar option, real scalar b_start, real scalar b_end, | real scalar seed_chunk) { struct boot_result scalar result struct did_data scalar dat_boot struct sa_point scalar boot_pt real matrix boot_mat real rowvector row_b real scalar n_iter, n_lead, b, l, n_failed, boot_has_valid // Set RNG seed for this chunk if specified if (args() >= 5 && !missing(seed_chunk)) { rseed(seed_chunk) } n_iter = b_end - b_start + 1 n_lead = cols(option.lead) // Return empty result for degenerate chunk if (n_iter <= 0) { result.estimates = J(0, 2 * n_lead, .) result.vcov = J(n_lead, 1, NULL) result.n_successful = 0 result.n_failed = 0 return(result) } boot_mat = J(0, 2 * n_lead, .) n_failed = 0 for (b = 1; b <= n_iter; b++) { // Unit-level bootstrap: resample entire units with replacement dat_boot = sample_panel(data) if (dat_boot.N == 0 || dat_boot.n_units < 2) { n_failed++ continue } // Compute time-weighted SA-DID and SA-sDID for all leads jointly. // sa_double_did() rebuilds Gmat, selects valid periods, applies pi_t // weights, and returns aggregated (DID[l], sDID[l]) for each lead l. boot_pt = sa_double_did(dat_boot, option) // Accept iteration if at least one lead has a valid estimate pair boot_has_valid = 0 for (l = 1; l <= n_lead; l++) { if (!missing(boot_pt.DID[l]) || !missing(boot_pt.sDID[l])) { boot_has_valid = 1 break } } if (!boot_has_valid) { n_failed++ continue } // Store (DID[l], sDID[l]) for each lead as a 2*n_lead row vector. // Preserving raw pairs (not pre-combined dDID) is essential: the // bootstrap covariance Sigma = Cov(DID^b, sDID^b) determines the // GMM optimal weight W* = Sigma^{-1} (paper equation 13). row_b = J(1, 2 * n_lead, .) for (l = 1; l <= n_lead; l++) { row_b[2*l - 1] = boot_pt.DID[l] row_b[2*l] = boot_pt.sDID[l] } boot_mat = boot_mat \ row_b } result.estimates = boot_mat result.vcov = J(n_lead, 1, NULL) result.n_successful = rows(boot_mat) result.n_failed = n_failed return(result) } // ---------------------------------------------------------------------------- // Module Verification Function // ---------------------------------------------------------------------------- /*--------------------------------------------------------------------------- * _did_bootstrap_loaded() - Verify module is loaded *---------------------------------------------------------------------------*/ void _did_bootstrap_loaded() { printf("{txt}did_bootstrap.mata loaded successfully\n") } end