Title

trop -- Triply Robust Panel Estimator

Syntax

trop depvar treatvar [if] [in] [weight], panelvar(varname) timevar(varname)
[options]

pweights are allowed and must be constant within a panel unit; see
Survey weights.

trop, version displays the installed package version and exits without
requiring depvar / treatvar.

Required
panelvar(varname) panel unit identifier

timevar(varname) time period identifier

Estimation Method
method(twostep|joint|local|global) estimation method; default is twostep, the paper's main Algorithm 2 estimator. local is an alias for twostep and global for joint, matching the paper terminology.

Lambda Grid Settings
grid_style(default|fine|extended) preset lambda grid; default is default

lambda_time_grid(numlist) custom lambda_time grid values

lambda_unit_grid(numlist) custom lambda_unit grid values

lambda_nn_grid(numlist) custom lambda_nn grid values

LOOCV Control
fixedlambda(numlist) skip LOOCV; use fixed lambda values (3 values)

twostep_loocv(cycling|exhaustive) twostep-method LOOCV search strategy; default is cycling

joint_loocv(cycling|exhaustive) joint-method LOOCV search strategy; default is exhaustive

tol(#) convergence tolerance; default is 1e-6

maxiter(#) maximum iterations; default is 500

Bootstrap Inference
bootstrap(#) bootstrap replications (paper Alg 3); default is 200; set 0 to skip

bsalpha(#) deprecated significance level alias; overrides level() when supplied

bsvariance(sample|paper) bootstrap variance denominator; default is sample (1/(B-1)); use paper for the Alg 3 denominator (1/B)

cimethod(percentile|t|normal) primary confidence interval type; default is percentile (Alg 3 step 6) when bootstrap is enabled, otherwise t

Covariates
covariates(varlist) covariates for Eq. 14 adjustment (paper Section 6.2)

Survey Design
strata(varname) stratification variable for Rao-Wu bootstrap

psu(varname) primary sampling unit variable

fpc(varname) finite population correction variable

nest nest PSU within strata

singleunit(centered|skip) lonely PSU handling strategy; default is skip

Other
seed(#) random number seed; default is 42

level(#) confidence level; default is c(level)

vlevel(#) verbose output level (0-4); default is 0

verbose display detailed progress information (equivalent to vlevel(2))

notiming suppress timing estimate display

Description

trop implements the Triply Robust Panel (TROP) estimator proposed by
Athey, Imbens, Qu, and Viviano (2025) for estimating the average treatment
effect on the treated (ATT) in panel data settings. The estimator combines
unit weights, time weights, and a nuclear-norm regularized low-rank factor
model to predict counterfactual outcomes for treated unit-time pairs.

The triple robustness property (Theorem 5.1) ensures that the bias of the
estimator is bounded by the product of three components: unit imbalance,
time imbalance, and misspecification of the regression adjustment. The
estimator is consistent if any one of the three components is zero
(Corollary 1).

Two estimation methods are available:

twostep (default) implements Algorithm 2 of the paper, which estimates
individual treatment effects for each treated unit-time pair (Eq. 13 /
Eq. 2 at the cell level), allowing for heterogeneous effects across units
and time periods, and then aggregates to ATT via Eq. 1. This is the
primary estimator exposed by trop and the recommended choice unless
a homogeneous-effect restriction is substantively justified.

joint implements the aggregation strategy described in Remark 6.1,
which assumes homogeneous treatment effects and estimates a single scalar
tau via weighted least squares. Because Remark 6.1 does not specify
concrete time/unit kernels, this method uses a natural adaptation of
Eq. 3: time weights decay from the centre of the treated block, and unit
weights decay in the RMSE of each unit's pre-treatment trajectory against
the average treated trajectory (see Remarks below). This shared-weight
extension is most appropriate under simultaneous adoption, where the same
treated block is meaningful across all treated units.

Options

Required

panelvar(varname) specifies the variable that identifies the panel
units (e.g., individuals, firms, countries).

timevar(varname) specifies the variable that identifies the time
periods.

Estimation Method

method(twostep|joint|local|global) specifies the estimation method.
Paper-consistent aliases are accepted: local is identical to twostep
(per-observation weights, Algorithm 2) and global is identical to
joint (shared-weight estimator, Remark 6.1). Any other value is
rejected with r(198).

twostep (default) uses per-observation weights following Algorithm 2.
This allows for heterogeneous treatment effects and is the recommended
approach. Use this when treatment effects may vary across units, across
post-treatment periods, or both. twostep supports arbitrary 0/1
treatment matrices, including staggered adoption and switching treatments
(units turning treatment on and off).

joint uses global weights following Remark 6.1. This assumes
homogeneous treatment effects but provides higher precision when the
assumption holds. Because the same weight system is shared across treated
cells, joint requires a simultaneous-adoption, absorbing design and
rejects staggered adoption and switching treatments with r(459); use
twostep for those patterns.

Lambda Grid Settings

grid_style(default|fine|extended) specifies a preset lambda grid for
LOOCV hyperparameter selection.

default uses 180 combinations (6x6x5) and is recommended for
most applications. The lambda_nn dimension is a five-point
log-decade ladder {0, 0.01, 0.1, 1, 10}, which excludes the
DID/TWFE corner (lambda_nn = +infinity) for predictable LOOCV cost.

fine uses 343 combinations (7x7x7) and inserts half-decade steps
(0.0316, 0.316) into the critical 0.01-1 band of the
lambda_nn grid and a 0.3 point into the lambda_time and
lambda_unit grids. Recommended for small panels (e.g. the paper's
Basque and West Germany applications) where the LOOCV objective
Q(lambda) surface is non-convex and coarser grids can produce
platform-dependent lambda selections.

extended uses 4,256 combinations (14x16x19) covering all optimal
values from Table 2 of the paper. This provides a finer search but is
computationally more expensive.

lambda_time_grid(numlist) specifies custom values for lambda_time.
Overrides grid_style() for this parameter. Values must be finite and
non-negative: Stata missing (.) is rejected at parse time because the
corresponding degenerate kernel (all weight collapsed onto the target
period) has no counterpart in paper Eq. (3).

lambda_unit_grid(numlist) specifies custom values for lambda_unit.
Same finiteness requirement as lambda_time_grid(): Stata missing is
rejected at parse time.

lambda_nn_grid(numlist) specifies custom values for lambda_nn (nuclear
norm penalty). Overrides grid_style() for this parameter. Stata
missing (.) is accepted and interpreted as +infinity (paper Eq. 2
remark), which zeros the low-rank component L and yields the classical
DID / TWFE special case. The grid_style(default) preset does not
include this corner; opt in via grid_style(extended) or a custom
lambda_nn_grid() that contains . or a sufficiently large value.

Performance caveat (large panels): lambda_nn values in the open interval
(0, 0.1) - most notably the 0.01 point of the default grid - are by
far the most expensive to evaluate. For 0 < lambda_nn < 0.1 the
nuclear-norm subproblem is solved by FISTA with the inner-iteration cap raised
to 50, and every inner step costs one full T-by-N SVD; by contrast
lambda_nn = 0 and lambda_nn >= 0.1 take cheap closed-form paths.
On a large panel a single interior candidate can therefore dominate the whole
search: at roughly 8,300 control cells one such point was measured tens of
thousands of times slower per cell than lambda_nn = 0 (about 12,600
ms/cell versus about 0.2 ms/cell), i.e. tens of hours for that one candidate
alone. Large-panel users are advised to supply a custom grid that avoids the
open interval, e.g. lambda_nn_grid(0 1 10), which keeps only the cheap
regimes.

joint_loocv(cycling|exhaustive) selects the LOOCV search strategy for
the joint method. Only effective when method(joint) is combined with
LOOCV (that is, fixedlambda() is not specified).

exhaustive (default) evaluates every (lambda_time, lambda_unit,
lambda_nn) combination in the Cartesian product of the three grids in
parallel and returns the exact argmin of Q(lambda) over the grid.
Cost is O(|grid|^3) and may be expensive on dense grids.

cycling runs a two-stage coordinate-descent search adapted from
Footnote 2 of Athey et al. (2025): univariate sweeps with extreme fixed
values produce initial estimates, then each parameter is updated cyclically
while the others are held at their current optimum. Cost is
O(|grid| * max_cycles). Prefer cycling on dense grids (grid_style(extended)
or larger custom grids) where the O(|grid|^3) exhaustive cost becomes
prohibitive; it matches exhaustive to within the LOOCV tie tolerance
on all benchmark panels but may miss non-convex local minima on small panels.

twostep_loocv(cycling|exhaustive) selects the LOOCV search strategy
for the twostep method. Only effective when method(twostep) is
combined with LOOCV.

cycling (default) runs coordinate descent over the per-observation
LOOCV objective, matching the historical Stata behaviour. Cost is
O(|grid| * max_cycles) per evaluated control cell, so it scales well even
with grid_style(extended) and is the recommended default on
medium-to-large panels (N * T &gt;= 400).

exhaustive enumerates the full Cartesian product of the three
lambda grids and returns the global grid minimum. A deterministic
lex-order tie-breaker
(largest lambda_nn &gt; smallest lambda_time &gt; smallest
lambda_unit) resolves near-ties within 1e-12 of the best score so
that cross-BLAS or cross-platform runs return identical lambdas.
Recommended for small panels (e.g. the paper's Basque and West Germany
applications, N * T &lt; 200) where the Q(lambda) surface is non-convex and
coordinate descent may stall at platform-dependent local minima. Example:
trop y d, panelvar(id) timevar(t) grid_style(fine) twostep_loocv(exhaustive).

Infinity values: Use . (missing) or 1e10 to represent
infinity in custom grids. For lambda_time and lambda_unit, infinity means
uniform weights. For lambda_nn, infinity disables the low-rank factor model.

LOOCV Control

Computational cost: LOOCV evaluates the cross-validation criterion for each
control observation across all grid points. For panels with N*T > 500 control
cells, consider using fixedlambda() with theoretically motivated values,
or reducing the grid with grid_style(default).

fixedlambda(numlist) specifies fixed values for (lambda_time,
lambda_unit, lambda_nn), bypassing LOOCV hyperparameter selection entirely.
Exactly 3 non-negative values must be provided.

Example: fixedlambda(0.5 1.0 0.1) sets lambda_time=0.5,
lambda_unit=1.0, lambda_nn=0.1. When this option is specified,
grid_style() is ignored.
tol() and maxiter() still apply to the alternating minimization
in the estimation step.

tol(#) specifies the convergence tolerance for the alternating
minimization algorithm. Default is 1e-6.

maxiter(#) specifies the maximum number of iterations for the
alternating minimization algorithm. Default is 500.

For large panels (N > 30 or T > 25), the default maxiter(500) may be
insufficient for convergence. Consider increasing to maxiter(1000) or
maxiter(2000) for such cases. Convergence status is reported in
e(converged).

Bootstrap Inference

bootstrap(#) specifies the number of bootstrap replications for
variance estimation following Algorithm 3 of the paper. Default is 200.
Set to 0 to skip bootstrap inference; the point estimate is still
returned but standard errors, p-values, and confidence intervals are not
computed. The bootstrap constructs each replicate by sampling N_0 control
units with replacement and N_1 treated units with replacement separately.

Performance note: Each bootstrap replication requires a full model re-estimation.
For panels with N > 20 units, expect ~5-10 seconds per replication.
Use bootstrap(30) for quick diagnostics or bootstrap(200) for
publication-quality inference.

Confidence intervals. Three candidate intervals are always computed
and stored whenever bootstrap is enabled; the option cimethod()
selects which one is promoted to the primary pair
e(ci_lower)/e(ci_upper).

(a) t-based CI e(ci_lower_t)/e(ci_upper_t) computed as
att ± t_df_r(1-alpha/2) * se, with the
bootstrap standard error and a t(N_1 - 1) reference distribution, where
N_1 is the number of ever-treated units. The unit-level stratified
resampling in Algorithm 3 makes N_1 the cluster count that governs the
small-sample df. When N_1 < 2 the t-wrap collapses to the normal wrap.

(b) normal-wrap CI e(ci_lower_normal)/e(ci_upper_normal)
using the standard-normal 1-alpha/2 quantile. Defined whenever the SE
is defined; identical to the t-wrap in the large-sample limit.

(c) percentile CI
e(ci_lower_percentile)/e(ci_upper_percentile)
taken directly from the alpha/2 and 1-alpha/2 quantiles of the bootstrap
empirical distribution (step 6 of Algorithm 3). Distribution-free under
the bootstrap exchangeability assumption, and invariant under monotone
reparameterisations of tau.

All three intervals share the same confidence level level(). The
primary interval is echoed first in the results table with a
[<method>] tag; the two non-primary intervals are echoed directly
below for comparison. Large discrepancies between the parametric and
percentile intervals suggest asymmetry or heavy tails in the bootstrap
distribution and warrant investigation.

cimethod(percentile | t | normal) selects which interval populates
the primary e(ci_lower)/e(ci_upper) pair. The default is
percentile whenever bootstrap > 0 (Algorithm 3 step 6) and
t otherwise. When bootstrap = 0 an explicit
cimethod(percentile) downgrades to t and emits a note; in that
case the downgrade trace is recorded in e(cimethod) as
"percentile->t". The three candidate pairs listed above are always
written to e() independently of cimethod() so downstream code
can switch primary without re-estimating.

bsvariance(sample|paper) selects the denominator used to build the
bootstrap standard error e(se).

sample (default) divides the sum of squared centred replicate
estimates by B - 1, yielding the Bessel-corrected sample variance.
This is the unbiased estimator for the bootstrap variance and is the
conventional choice in bootstrap practice and in Stata's other bootstrap
commands.

paper divides by B, matching the denominator that appears
explicitly in Algorithm 3 of Athey, Imbens, Qu, and Viviano (2025). Select
this option to reproduce the paper's reported standard errors exactly.

The two choices differ by a factor of B / (B - 1) inside the
variance (equivalently, sqrt(B / (B - 1)) inside the SE), which is
below 0.5% when B = 200. Percentile confidence bounds are unaffected.
The chosen denominator is echoed in e(bsvariance).

bsalpha(#) is a deprecated alias for the significance level
alpha = 1 - level()/100. It is retained for backward compatibility.
When specified, bsalpha() overrides level() and a note is
printed at estimation time. New code should prefer level().

Covariates

covariates(varlist) specifies covariates to include in estimation
following the Equation 14 adjustment of Athey et al. (2025) Section 6.2.
Covariates enter the alternating minimisation as an additional WLS step that
updates a gamma coefficient vector. All specified variables must be numeric
and must not contain missing values in the estimation sample. The covariate
variables must not overlap with the outcome or treatment variables.

Stored results: e(gamma) (1 x p matrix of covariate coefficients),
e(n_covariates) (scalar), e(covariates) (macro listing the
covariate variable names).

Example: trop y d, panelvar(id) timevar(t) covariates(x1 x2 x3)

Survey Design

strata(varname) specifies the stratification variable for Rao-Wu
bootstrap variance estimation. Required when psu() or fpc() is
specified.

psu(varname) specifies the primary sampling unit variable. Required
when strata() is specified.

fpc(varname) specifies the finite population correction variable.
Optional; when supplied, the Rao-Wu bootstrap incorporates the sampling
fraction adjustment.

nest requests that PSUs be nested within strata (relevant when PSU
identifiers are not globally unique).

singleunit(centered|skip) specifies the strategy for handling strata
containing a single primary sampling unit (lonely PSU) during Rao-Wu
bootstrap variance estimation.

skip (default) omits the lonely-PSU stratum from the bootstrap
variance contribution entirely. This preserves backward compatibility with
pre-1.2 behaviour and is appropriate when the lonely stratum contributes
negligible variance.

centered uses the grand-mean centring correction: the single PSU's
contribution is computed relative to the overall weighted mean rather than
the within-stratum mean. This is the Stata svyset convention for
singleunit(centered) and avoids discarding information from the
stratum.

Other

seed(#) specifies the random number seed for reproducibility.
Default is 42.

level(#) specifies the confidence level for confidence intervals.
Default is c(level), typically 95.

vlevel(#) specifies the verbosity level for output control.
Five levels are available:

0 -- silent (no progress messages)

1 -- minimal (start/completion only)

2 -- detailed (LOOCV diagnostics; same as verbose)

3 -- debug (intermediate values)

4 -- trace (all internal steps)

Default is 0. When verbose is specified without vlevel(),
the level is set to 2.

verbose displays detailed progress information including LOOCV
diagnostics, convergence status, and timing information. Equivalent to
vlevel(2).

notiming suppresses the timing estimate display that trop
normally emits on large panels (N*T >= 50,000). Use this when the
progress estimate is not wanted (e.g., in automated scripts).

Remarks

Weight structures

The two estimation methods use fundamentally different weight structures,
following Algorithm 2 and Remark 6.1 of the paper respectively. The
twostep path is the paper's main estimator; joint is an
additional shared-weight restriction layered on top of the same TROP
ingredients:

method(twostep) — Per-observation weights (Eq. 3, Algorithm 2):

Each treated observation (i,t) receives its own set of weights. Time weights
measure distance from the specific target period:
theta_s = exp(-lambda_time * |t - s|). Unit weights measure the RMSE
between each control unit j and the specific target unit i over common
control periods, excluding the target period t:
omega_j = exp(-lambda_unit * dist_{-t}(j, i)). This permits heterogeneous
treatment effects across units and periods.

method(joint) — Global weights (Remark 6.1):

A single set of weights is shared across all treated observations. Time
weights measure distance to the center of the treated block:
delta_time[t] = exp(-lambda_time * |t - center|). Unit weights measure
the RMSE from each unit's pre-treatment trajectory to the average treated
trajectory: delta_unit[j] = exp(-lambda_unit * RMSE(j, avg_treated)).
This assumes homogeneous treatment effects, is computationally more
efficient, and is intended for panels where the treated units share the
same adoption window.

Note. Paper Remark 6.1 defines the homogeneous-tau aggregation but
does not prescribe a specific time/unit kernel under the shared-weight
setting; the post-block midpoint and the pre-period trajectory RMSE used
above are engineering choices adopted by trop, not formulas
restated from the paper. See "Methods and formulas" for the explicit
delta_time / delta_unit definitions.

When to use each method:

Use twostep (default) when treatment effects may vary across
units and time periods, or when the panel has complex treatment assignment
patterns.

Use joint when a single shared treatment effect and a common
treated block are substantively appropriate.

Interpretation of e(mu)

The content of e(mu) differs between estimation methods:

method(twostep): e(mu) is set to . (missing). The twostep
model is Y(0) = alpha_i + beta_t + L_ without an explicit global
intercept. Unit fixed effects are stored in e(alpha).

method(joint): e(mu) returns the global intercept mu (a scalar).
The joint model is Y(0) = mu + alpha_i + beta_t + L_ with
identification constraints alpha_1 = beta_1 = 0.

Both parameterizations yield mathematically equivalent counterfactual
predictions.

Lambda grid comparison

Grid Style lambda_time lambda_unit lambda_nn Total

default 6 values 6 values 5 values 180
fine 7 values 7 values 7 values 343
extended 14 values 16 values 19 values 4,256

The extended grid covers all optimal values from Table 2 of the paper:

Dataset lambda_unit lambda_time lambda_nn

CPS log-wage 0 0.1 0.9
CPS urate 1.6 0.35 0.011
PWT 0.3 0.4 0.006
Germany 1.2 0.2 0.011
Basque 0 0.35 0.006
Smoking 0.25 0.4 0.011
Boatlift 0.2 0.2 0.151

Survey weights

trop accepts pweight (probability / sampling weights) to support
survey-design-aware ATT aggregation. When [pweight=wgt] is
supplied, the per-cell treatment effects are aggregated with a weighted mean:

ATT_hat = sum_ W=1 w_i * tau_ / sum_ W=1 w_i

where w_i is the pweight attached to unit i. The pweight variable
must be strictly positive and constant within each panel unit; a cell-
varying weight triggers a hard error (r(459)).

Scope of the weighted path:

- ATT aggregation and the bootstrap empirical distribution use the weighted
mean. Each bootstrap draw inherits the original pweight of every resampled
unit.

- The nuisance-parameter fits (alpha, beta, L) and the LOOCV objective
Q(lambda) remain unweighted. This mirrors the design-based construction:
the cross-validation score is an in-sample prediction criterion that does
not depend on the target population.

- Only pweights are supported. aweight, fweight, and
iweight are rejected (r(101)). Bootstrap inference is pweight-
only — strata, PSU, and finite-population-correction Rao-Wu resampling are
not implemented.

- Running with [pweight=1] reproduces the unweighted output to
machine precision; this is a regression-test guarantee.

Stored metadata: e(weight_var) records the pweight variable name,
e(wtype) is set to pweight, and e(wexp) is = wgt.
trop_bootstrap re-reads these e() entries to dispatch the weighted
bootstrap path automatically.

Examples

Setup: Generate test panel data

. clear all

. set seed 12345

. set obs 1000

. gen id = ceil(_n/10)

. bysort id: gen t = _n

. gen y = rnormal() + 0.5*(id>80)*(t>7)

. gen d = (id > 80 & t > 7)

Quick start (fixed lambda, runs in seconds)

. trop y d, panelvar(id) timevar(t) fixedlambda(0.5 1.0 0.1) seed(42)

Use fixedlambda() to bypass LOOCV and get immediate results.
Recommended for first-time exploration and large panels.

Twostep method with LOOCV (default)

. trop y d, panelvar(id) timevar(t) seed(42)

Note: LOOCV hyperparameter search can be slow on large panels.
Consider fixedlambda() or use a custom lambda_nn_grid()
that avoids the (0, 0.1) region for faster convergence.

Joint method

. trop y d, panelvar(id) timevar(t) method(joint) seed(42)

Use this only when a single shared treatment effect and a common
treated block are substantively appropriate.

Joint method with exhaustive LOOCV (guaranteed grid argmin)

. trop y d, panelvar(id) timevar(t) method(joint) joint_loocv(exhaustive) seed(42)

With bootstrap inference

. trop y d, panelvar(id) timevar(t) bootstrap(200) seed(42)

Using extended grid for finer search

. trop y d, panelvar(id) timevar(t) grid_style(extended) seed(42)

Custom lambda grid

. trop y d, panelvar(id) timevar(t) lambda_nn_grid(0 0.1 1e10) seed(42)

Fixed lambda (skip LOOCV)

. trop y d, panelvar(id) timevar(t) fixedlambda(0.5 1.0 0.1) seed(42)

Verbose output

. trop y d, panelvar(id) timevar(t) seed(42) verbose

Row-labelled fixed effects

e(alpha) and e(beta) carry the sorted unique values of
panelvar / timevar as matrix row names, so the fixed effects
can be read directly against their original identifiers:

. trop y d, panelvar(id) timevar(t) seed(42)

. matrix list e(alpha) // rows: _1, _2, ... (numeric IDs get a leading underscore)

. matrix list e(beta) // rows: time labels (e.g. _1989, _1990, ...)

Postestimation diagnostics

After trop, the companion command estat
exposes thirteen diagnostics including hyperparameter sensitivity,
variance-covariance display, event-study analysis, pre-trend testing,
and a triple-robustness bias decomposition anchored to paper Theorem 5.1:

. trop y d, panelvar(id) timevar(t) bootstrap(200) seed(42)

. estat summarize // sample + treatment overview

. estat loocv // hyperparameter search trace

. estat weights // unit / time weight distribution

. estat factors // SVD of the factor matrix L

. estat bootstrap // bootstrap distribution summary

. estat sensitivity // LOOCV grid sensitivity summary

. estat vce // variance-covariance display

. estat triplerob // Theorem 5.1 bias decomposition

Per-cell treatment effects via e(tau_matrix)

e(tau_matrix) is a T x N matrix indexed by (time, panel) with
treatment effects in treated cells and Stata missing (.) elsewhere.
Use it to retrieve tau by coordinate or to aggregate effects by row/column:

. trop y d, panelvar(id) timevar(t) seed(42)

. matrix tauM = e(tau_matrix) // T x N panel-shaped tau

. di "tau at (t=8, i=85) = " tauM[8, 85]

Average treatment effect within each post-treatment period (event-time path):

. mata: tau = st_matrix("e(tau_matrix)")

. mata: present = (tau :< .) // 1 if treated cell

. mata: att_t = rowsum(editmissing(tau, 0)) :/ rowsum(present)

. mata: att_t // T x 1 (missing for pre-treatment rows)

Bundled Datasets

The trop package ships six panel datasets. After installing from
SSC, download them once into your working directory with net get:

. net get trop

This places the six .dta files in the current directory, ready to
use:

Dataset Source N T Treatment

cps_logwage.dta CPS log-wage (min wage) 50 40 real
cps_urate.dta CPS urate (min wage) 50 40 real
pwt_loggdp.dta PWT log-GDP (democracy) 111 48 real
germany_gdp.dta Abadie (2003) West Germany 17 44 d=0
basque_gdp.dta Abadie (2003) Basque Country 18 43 d=0
smoking_packs.dta California Prop 99 39 31 d=0

Important: germany_gdp, basque_gdp, and smoking_packs
have d = 0 throughout. They are raw outcome panels designed for
semi-synthetic simulation (as in Athey et al. 2025, Table 1). You must
assign a treatment indicator before running trop. See Example 4
below for a worked demonstration.

Complete Data Analysis Examples

The following four examples demonstrate end-to-end workflows using real
and simulated panel data. The bundled datasets are downloaded once with
net get trop (see above) and then loaded with use. As a
convenience, trop_data name downloads and loads a single
dataset in one step (used in Example 1 below).

Example 1: CPS wage data -- basic twostep (local) analysis

This example reproduces part of the paper's Table 1 analysis of minimum
wage effects on log wages across U.S. states. The CPS dataset contains
50 states observed from 1979 to 2018; 8 states are treated at t=2018.
The outcome variable y is the log average weekly wage.

Step 1: Download example datasets (one-time, after installation)

. net get trop

Load the CPS log-wage panel

. use cps_logwage.dta, clear

Alternatively, trop_data cps_logwage downloads and loads the
dataset in one step.

Step 2: Estimate ATT using the twostep (local) method with bootstrap

. trop y d, panelvar(id) timevar(t) method(twostep) bootstrap(200) seed(42)

Step 3: Examine stored results

. ereturn list

. display "ATT = " e(att)

. display "SE = " e(se)

. display "95%% CI: [" e(ci_lower) ", " e(ci_upper) "]"

. display "p-value = " e(pvalue)

Step 4: Inspect selected hyperparameters

. display "lambda_time = " e(lambda_time)

. display "lambda_unit = " e(lambda_unit)

. display "lambda_nn = " e(lambda_nn)

Step 5: Postestimation diagnostics

. estat summarize

. estat loocv

. estat bootstrap

Interpretation. The ATT estimates the average treatment effect on
the treated (paper Eq. 1): the mean counterfactual-adjusted impact of
the policy change on log wages for the 8 treated states in t=2018.
A positive ATT implies the treatment raised wages relative to the
predicted no-treatment trajectory. The three-way robust construction
(Theorem 5.1) ensures that the estimate remains valid even if the unit
weights, time weights, or the nuclear-norm model alone fail to fully
balance the treated and control groups.

Example 2: PWT GDP data -- joint (global) method with advanced options

The Penn World Table dataset measures log GDP per capita across 111
countries from 1960 to 2007; 29 countries receive treatment at t=2007
(a simultaneous-adoption design). The larger panel and common adoption
window make this an ideal setting for the joint estimator.

Step 1: Load PWT data

. use pwt_loggdp.dta, clear

. summarize

. display "N_units = " r(N) / 48 " | T = 48 | Treated = 29"

Step 2: Joint estimation with fixed lambda (paper Table 2 values)

. trop y d, panelvar(id) timevar(t) method(joint) fixedlambda(0.4 0.3 0.006) bootstrap(200) seed(42)

. display "ATT (joint, fixed lambda) = " e(att)

Step 3: Compare with LOOCV-selected lambda using fine grid

. trop y d, panelvar(id) timevar(t) method(joint) grid_style(fine) bootstrap(200) seed(42)

. display "ATT (joint, fine grid LOOCV) = " e(att)

. display "Selected: lt=" e(lambda_time) " lu=" e(lambda_unit) " lnn=" e(lambda_nn)

Step 4: Exhaustive vs. cycling LOOCV comparison

. trop y d, panelvar(id) timevar(t) method(joint) joint_loocv(exhaustive) seed(42)

. scalar att_exh = e(att)

. trop y d, panelvar(id) timevar(t) method(joint) joint_loocv(cycling) seed(42)

. scalar att_cyc = e(att)

. display "ATT difference (exhaustive - cycling) = " att_exh - att_cyc

Step 5: Weight diagnostics

. trop y d, panelvar(id) timevar(t) method(joint) bootstrap(200) seed(42)

. estat weights

. estat triplerob

Interpretation. The joint method assumes a homogeneous treatment
effect tau across all treated cells (Remark 6.1). Because 29 countries
adopt simultaneously, the shared-weight construction is substantively
appropriate. The estat triplerob decomposition shows whether the
triple-robustness property delivers bias reduction via unit balance,
time balance, or regression adjustment alone. Comparing fixed lambda
(from the paper's Table 2) with LOOCV-selected lambda verifies that
the cross-validation procedure finds values close to the paper's optimal.

Example 3: Simulated data -- verification with known true effect

This example constructs a panel with a known treatment effect (tau=2.0)
and demonstrates that trop recovers it accurately. The simulation
illustrates the three-way robustness by showing convergence to the
true effect even under various model misspecification scenarios.

Step 1: Generate a balanced panel with known DGP

. clear all

. set seed 2025

. local N_units = 30

. local T_periods = 20

. local N_treated = 5

. local T_treat = 16

. local true_tau = 2.0

. set obs `= `N_units' * `T_periods''

. gen id = ceil(_n / `T_periods')

. bysort id: gen t = _n

Step 2: Build Y(0) = alpha_i + beta_t + L_ + noise

. gen alpha_i = id * 0.3

. gen beta_t = t * 0.1

. gen L_it = sin(id * 0.5) * cos(t * 0.3)

. gen epsilon = rnormal(0, 0.2)

. gen d = (id > `N_units' - `N_treated') & (t >= `T_treat')

. gen y = alpha_i + beta_t + L_it + epsilon + `true_tau' * d

Step 3: Estimate with twostep method

. trop y d, panelvar(id) timevar(t) method(twostep) bootstrap(200) seed(42)

. display "True tau = `true_tau'"

. display "Estimated ATT = " e(att)

. display "Bias = " e(att) - `true_tau'

. display "SE = " e(se)

. display "Coverage: " cond(e(ci_lower) <= `true_tau' & `true_tau' <= e(ci_upper), "YES", "NO")

Step 4: Repeat with joint method for comparison

. trop y d, panelvar(id) timevar(t) method(joint) bootstrap(200) seed(42)

. display "Joint ATT = " e(att)

. display "Joint Bias = " e(att) - `true_tau'

Step 5: Verify triple robustness with the bias decomposition

. trop y d, panelvar(id) timevar(t) bootstrap(200) seed(42)

. estat triplerob

Interpretation. With a known true_tau = 2.0, the estimated ATT
should be close to 2.0 and the 95%% confidence interval should cover it.
The bias decomposition (estat triplerob) reveals the three
components of Theorem 5.1:

||Delta^u||_2 -- unit imbalance

||Delta^t||_2 -- time imbalance

||B||_* -- regression adjustment misspecification

The TROP estimator is consistent if any one of these three quantities is
zero (Corollary 1). In this simulated example, the true DGP satisfies
Y(0) = alpha_i + beta_t + L_, so the nuclear-norm regression
adjustment is correctly specified (||B||_* approx 0), guaranteeing
consistency regardless of finite-sample weight imbalance.

Example 4: Basque GDP data -- semi-synthetic simulation (d=0 panel)

The Basque Country GDP dataset (Abadie & Gardeazabal 2003) ships with
d = 0 throughout because it is a raw outcome panel. To replicate
the paper's Table 1 Monte Carlo exercise, you assign treatment randomly
or designate a specific unit as treated. This example demonstrates both
approaches.

Approach A: Treated-unit simulation (designate Basque Country)

. use basque_gdp.dta, clear

. tab d

. * All d=0 -- assign treatment to unit 17 (Basque Country) post-1970

. replace d = (id == 17) & (t >= 1970)

. tab d

. trop y d, panelvar(id) timevar(t) method(joint) fixedlambda(0.35 0 0.006) seed(42)

. display "ATT (Basque treated-unit) = " e(att)

Approach B: Random treatment assignment (Table 1 design)

. use basque_gdp.dta, clear

. * Randomly assign 3 units as treated from t >= 1975

. set seed 2025

. bysort id: gen tag = (_n == 1)

. gen u = runiform() if tag

. bysort id (u): replace u = u[1]

. sort u

. egen rank = group(u)

. replace d = (rank <= 3) & (t >= 1975)

. drop tag u rank

. trop y d, panelvar(id) timevar(t) method(twostep) fixedlambda(0.35 0 0.006) seed(42)

. display "ATT (random treatment) = " e(att)

Interpretation. Because the true treatment effect is zero (the
treatment was assigned artificially to a pre-existing observational
panel), the estimated ATT should be close to zero. Systematic deviation
from zero indicates finite-sample bias. The paper's Table 1 reports the
normalized RMSE across many such random draws to compare estimator
performance. The basque_gdp, germany_gdp, and
smoking_packs datasets are all used this way.

Stored results

trop stores the following in e():

Scalars

e(att) ATT point estimate

e(se) bootstrap standard error

e(t) t statistic (att/se)

e(ci_lower) primary confidence interval lower bound; matches one of the three candidates below according to e(cimethod)

e(ci_upper) primary confidence interval upper bound

e(pvalue) two-sided p-value for the primary interval (t-based p-value when e(cimethod)==percentile)

e(df_r) degrees of freedom for the t-based reference distribution = max(1, N_1 - 1) where N_1 is the ever-treated unit count; missing when N_1 < 2

e(ci_lower_t) t-wrap CI lower bound using SE and t(df_r)

e(ci_upper_t) t-wrap CI upper bound

e(pvalue_t) two-sided p-value from the t-wrap

e(ci_lower_normal) normal-wrap CI lower bound using SE and N(0,1)

e(ci_upper_normal) normal-wrap CI upper bound

e(pvalue_normal) two-sided p-value from the normal-wrap

e(ci_lower_percentile) bootstrap percentile CI lower bound (Algorithm 3 step 6)

e(ci_upper_percentile) bootstrap percentile CI upper bound

e(mu) global intercept (. for twostep, scalar for joint)

e(lambda_time) selected or fixed lambda_time

e(lambda_unit) selected or fixed lambda_unit

e(lambda_nn) selected or fixed lambda_nn

e(stage1_lambda_time) Stage-1 univariate argmin for lambda_time (paper Footnote 2; cycling LOOCV only, missing for exhaustive)

e(stage1_lambda_unit) Stage-1 univariate argmin for lambda_unit (paper Footnote 2; cycling LOOCV only)

e(stage1_lambda_nn) Stage-1 univariate argmin for lambda_nn (paper Footnote 2; cycling LOOCV only)

e(loocv_score) LOOCV objective Q(lambda_hat) = sum over (i,t) with W_=0 and finite Y_ of tau_hat_^loocv(lambda_hat)^2 (paper Eq. 5 restricted to observed control cells)

e(loocv_n_valid) number of successful LOOCV fits

e(loocv_n_attempted) total LOOCV fit attempts (= every D=0 cell, paper Eq. 5)

e(loocv_fail_rate) LOOCV failure rate (0 to 1)

e(loocv_rmse) LOOCV root mean squared error = sqrt(Q(lambda_hat) / n_valid)

e(loocv_used) 1 if LOOCV was performed, 0 if skipped

e(loocv_first_failed_t) time index of the first LOOCV fit that failed (0-based; -1 if none)

e(loocv_first_failed_i) unit index of the first LOOCV fit that failed (0-based; -1 if none)

e(n_lambda_time) number of lambda_time grid values

e(n_lambda_unit) number of lambda_unit grid values

e(n_lambda_nn) number of lambda_nn grid values

e(n_grid_combinations) total grid combinations

e(n_grid_per_cycle) grid points per coordinate descent cycle

e(effective_rank) effective rank of L (sum(s)/s[1])

e(n_iterations) number of iterations

e(converged) convergence indicator (1=yes, 0=no)

e(n_obs_estimated) successfully estimated observations (twostep)

e(n_obs_failed) failed observations (twostep)

e(N_units) number of panel units N

e(N_periods) number of time periods T

e(N_treated_units) number of ever-treated units N_1 (cluster count for Algorithm 3)

e(N_obs) total number of observations

e(N_treat) treated unit-period cell count = sum(W=1) across the estimation sample; legacy alias of e(N_treated_obs)

e(N_treated) length of e(tau) = successfully estimated treated cells; equals e(N_treated_obs) when every treated cell converges

e(N_treated_obs) treated unit-period cell count = sum(W=1) across the estimation sample. Each treated (i,t) cell contributes one tau; the ATT in Eq. 1 is a mean (or weighted mean under pweight) over these cells. Distinct from e(N_treated_units) which counts unique ever-treated units.

e(N_control) number of control observations (W=0)

e(N_control_units) number of never-treated units N_0

e(T_treat_periods) number of periods with treatment

e(bootstrap_reps) bootstrap replications requested

e(n_bootstrap_valid) successful bootstrap iterations

e(bootstrap_fail_rate) bootstrap failure rate (0 to 1); parallels e(loocv_fail_rate). An advisory note is printed when this exceeds 5%; exceeding 50% aborts with rc=504. Missing when bootstrap(0).

e(alpha_level) significance level for CI

e(level) confidence level (e.g., 95)

e(seed) random number seed used

e(balanced) 1 if panel is balanced, 0 otherwise

e(miss_rate) fraction of missing observations

e(min_pre_treated) minimum pre-treatment periods for treated units

e(min_valid_pairs) minimum common control periods across pairs

e(has_switching) 1 if treatment switching detected

e(max_switches) maximum treatment switches observed

e(data_validated) 1 if data validation passed

e(time_min) minimum time period value

e(time_max) maximum time period value

e(time_range) range of time periods

e(n_pre_periods) number of pre-treatment periods

e(n_post_periods) number of post-treatment periods

e(condition_number) condition number of the WLS design matrix; large values (>1e10) indicate ill-conditioning

e(n_covariates) number of covariates (0 if none)

e(deff_weights) Kish design effect of pweights = N * sum(w_i^2) / (sum(w_i))^2; measures efficiency loss from unequal weighting. Only stored when pweights are supplied.

e(max_fh) maximum finite-population fraction f_h across strata (Rao-Wu bootstrap); only stored when survey weights are used

e(n_high_fpc) number of strata where f_h > 0.5 (high FPC detection); only stored when survey weights are used

Matrices

e(b) coefficient vector (1x1, ATT)

e(V) variance-covariance matrix (1x1; bootstrap only)

e(alpha) unit fixed effects (N x 1). Row names are the sorted unique values of e(panelvar) on the estimation sample (sanitised to valid Stata matrix identifiers), so matrix list e(alpha) reports the fixed effects keyed by the original panel identifier. twostep: algorithmic mean across treated observations (Algorithm 2 step 2 fits a separate alpha^{i,t} for each treated cell; e(alpha) averages them); joint: single model estimate. See e(alpha_semantics).

e(beta) time fixed effects (T x 1). Row names are the sorted unique values of e(timevar) on the estimation sample (sanitised to valid Stata matrix identifiers). twostep: algorithmic mean across treated observations (same convention as e(alpha)); joint: single model estimate.

e(factor_matrix) low-rank factor matrix L (T x N). twostep: algorithmic mean across treated observations; joint: single model estimate.

e(tau) per-cell treatment effects (N_treated x 1); populated for both twostep and joint. For joint the vector carries the single scalar tau replicated, so mean(e(tau)) == e(att) to machine precision in either method.

e(tau_matrix) treatment effects arranged as a T x N panel-shaped matrix with Stata missing (.) in untreated cells; omitted when the underlying panel metadata is unavailable.

e(converged_by_obs) convergence flag per treated cell (N_treated x 1; 1 = converged, 0 = hit maxiter(), -1 = solver error); twostep only.

e(n_iters_by_obs) iteration count per treated cell (N_treated x 1); twostep only.

e(theta) time weights (T x 1; twostep only)

e(omega) unit weights (N x 1; twostep only)

e(delta_time) global time weights (T x 1; joint only)

e(delta_unit) global unit weights (N x 1; joint only)

e(bootstrap_estimates) bootstrap ATT estimates (B x 1)

e(lambda_time_grid) lambda_time grid values searched

e(lambda_unit_grid) lambda_unit grid values searched

e(lambda_nn_grid) lambda_nn grid values searched

e(gamma) covariate coefficients (1 x p; only when covariates() is specified)

e(lambda_grid) Cartesian product of lambda grids (K x 3)

e(cv_curve) LOOCV scores at grid points (K x 4)

Macros

e(cmd) "trop"

e(cmdline) full command line as typed

e(method) "twostep" or "joint"

e(depvar) dependent variable name

e(treatvar) treatment variable name

e(panelvar) panel variable name

e(timevar) time variable name

e(grid_style) grid style actually used: default, fine, extended, or custom when one or more lambda_*_grid() options override a preset dimension

e(twostep_loocv) twostep LOOCV strategy: cycling or exhaustive

e(joint_loocv) joint LOOCV strategy: cycling or exhaustive

e(alpha_semantics) interpretation tag for e(alpha)/e(beta)/e(factor_matrix): obs_average (twostep, averaged across treated cells) or single_model (joint)

e(treatment_pattern) treatment pattern detected

e(vcetype) "Bootstrap" (with bootstrap)

e(bsvariance) bootstrap variance denominator actually used: sample (default, 1/(B-1)) or paper (Alg 3, 1/B)

e(cimethod) primary confidence-interval method selected by cimethod(): percentile, t, or normal. When an explicit request was downgraded (e.g. cimethod(percentile) with bootstrap(0)), the trace is reported as "percentile->t".

e(estat_cmd) "trop_estat"

e(title) "TROP Estimator"

e(predict) "trop_p"

e(weight_var) pweight variable name (empty when no [pweight=] was supplied); see Survey weights

e(wtype) weight type tag; set to pweight when weights were supplied, otherwise empty

e(wexp) weight expression as typed (e.g. = wgt); empty without weights

e(covariates) space-separated list of covariate variable names (empty when none)

e(spec_string) specification string recording the estimation call parameters (method, lambda values, covariates) for reproducibility

e(strata_var) stratification variable name (survey design only)

e(psu_var) PSU variable name (survey design only)

e(fpc_var) FPC variable name (survey design only)

e(bootstrap_type) bootstrap type: standard or rao_wu

Functions

e(sample) marks estimation sample

Methods and formulas

The TROP estimator predicts counterfactual outcomes for treated unit-time
pairs using a working model Y_(0) = alpha_i + beta_t + L_, where
alpha_i are unit fixed effects, beta_t are time fixed effects, and L is a
low-rank factor component.

Unless otherwise noted, the formulas in this section describe the
twostep estimator, which is the paper's main Algorithm 2 path. The
joint option replaces the cell-specific tau_ fits with a single
shared tau under the homogeneous-effect restriction of Remark 6.1.

For each treated observation (i,t), the parameters are estimated by solving
a weighted nuclear-norm penalized regression (Eq. 2):

(alpha_hat, beta_hat, L_hat) = argmin_ sum_
theta_s^ * omega_j^ * (1-W_) *
(Y_ - alpha_j - beta_s - L_)^2 + lambda_nn * ||L||_*

where ||L||_* denotes the nuclear norm of L. The weights are defined as
exponential distance-decay functions (Eq. 3):

theta_s^ = exp(-lambda_time * |t - s|)

omega_j^ = exp(-lambda_unit * dist_{-t}(j, i))

The unit distance is the RMSE between the outcome trajectories of units j
and i over common control periods, excluding the target period t:

dist_{-t}(j, i) = sqrt( sum_u 1(1-W_)(1-W_)(Y_-Y_)^2
/ sum_u 1(1-W_)(1-W_) )

The treatment effect for each treated observation (i,t) is estimated as:

tau_hat_ = Y_ - alpha_hat_i - beta_hat_t - L_hat_

The ATT is the simple average of individual treatment effects over all
treated unit-time pairs (paper Eq. 1 / Algorithm 2 step 6):

tau_hat = (1 / sum_ W_) * sum_ W_ * tau_hat_

Joint method specifics (Remark 6.1). Under the homogeneous-treatment
restriction of paper Remark 6.1, joint replaces the cell-specific
fits with a single weighted least-squares problem over a shared weight
matrix delta = delta_time * delta_unit'. Paper Remark 6.1 sketches the
aggregation but does not prescribe concrete time / unit kernels, so
trop adopts the following engineering choice (pinned by the released
numerical baseline):

delta_time[s] = exp(-lambda_time * |s - center|),
center = T - T_post / 2

delta_unit[j] = exp(-lambda_unit * RMSE_pre(j, avg_treated))

where center is the midpoint of the treated block, T_post is the
shared post-period count, and RMSE_pre(j, avg_treated) is the RMSE
of unit j's pre-treatment outcomes against the period-wise average treated
trajectory. These distance definitions are a reasonable adaptation of
Eq. (3) to the shared-weight setting; the post-block midpoint and the
pre-period trajectory RMSE are specific engineering choices, not formulas
prescribed by the paper. After fitting (mu, alpha, beta, L) by
weighted nuclear-norm regression with weight matrix delta * (1 - W), the
scalar tau_hat is recovered post-hoc as the weighted mean residual on
treated cells (paper Eq. 10 identity for the homogeneous-tau case).

The triple robustness property (Theorem 5.1) states that the bias satisfies:

|E[tau_hat - tau | L]| <= ||Delta^u||_2 * ||Delta^t||_2 * ||B||_*

where Delta^u is unit imbalance (the discrepancy between the weighted
average of control unit loadings and the treated unit loading), Delta^t is
time imbalance (the analogous discrepancy for time factor loadings), and B
captures regression adjustment misspecification. The estimator is
consistent if any one of the three components is zero (Corollary 1):
(a) balance over unit loadings, (b) balance over factor loadings, or
(c) correct regression adjustment specification.

Tuning parameters (lambda_time, lambda_unit, lambda_nn) are selected via
leave-one-out cross-validation (LOOCV) minimizing (Eq. 5):

Q(lambda) = sum_ (1 - W_) * (tau_hat_(lambda))^2

The grid search uses coordinate descent (footnote 2 of the paper): each
parameter is optimized in turn while holding the other two at their current
optimal values. In trop, this paper-style cycling search remains
available for both methods, while exhaustive Cartesian search is also
exposed as an engineering option when users prefer the exact grid argmin.

Bootstrap variance estimation follows Algorithm 3: units are resampled
with replacement within the treated and control groups separately (paper
Alg 3 step 4), and the TROP estimator is recomputed on each bootstrap
sample holding the LOOCV-selected lambda_hat fixed (paper Alg 3 step 5;
strict re-LOOCV per draw is out of scope for the current release).

The paper writes the bootstrap variance as

V_hat_paper = (1/B) * sum_^ (tau_hat^{(b)} - tau_bar)^2

By default trop reports the Bessel-corrected sample variance
instead, which is the unbiased estimator of the bootstrap variance and
matches the default inferential convention used by the package for
bootstrap SEs:

V_hat_trop = (1/(B-1)) * sum_^ (tau_hat^{(b)} - tau_bar)^2

The denominator is user-selectable through bsvariance(sample|paper):
the default sample uses 1/(B-1) and matches V_hat_trop, while
paper uses 1/B and reproduces V_hat_paper exactly. The two
choices differ by a factor B/(B-1) inside the variance (about 0.5% at the
default B = 200) and coincide as B grows; e(bsvariance) records
the active choice. tau_hat^{(b)} is the ATT estimate from the b-th
bootstrap sample, tau_bar is the mean of the bootstrap estimates, and the
bootstrap standard error e(se) is the square root of the selected
variance. The percentile CI
e(ci_lower_percentile)/e(ci_upper_percentile) is taken from
the alpha/2 and 1-alpha/2 sample quantiles of the bootstrap distribution,
using linear interpolation between adjacent order statistics (the
(n-1)*p fractional-index convention) and is unaffected by the denominator
choice.

Implementation notes

The trop core is a precompiled Rust library invoked via a C plugin.
The numerical algorithm follows Athey, Imbens, Qu and Viviano (2025) Eq. 2
and Algorithms 1-3. Several implementation choices that a careful reader
might want to inspect are listed below; each one is pinned by a regression
test so future refactors cannot silently regress.

1. FISTA adaptive restart disabled (rust/src/estimation.rs). The
nuclear-norm proximal solver does not use the monotone gradient-restart
scheme of O'Donoghue & Candes (2015). Although the restart criterion can
eliminate momentum oscillations in theory, it fires too aggressively on
small panels and prevents convergence. The reference Python implementation
does not use restart either, so it is disabled to maintain numerical
consistency. Locked in by tests/test_fista_restart_stability.do.

2. LAPACK dgelsd for the weighted least-squares step
(rust/src/estimation.rs). The SVD-based minimum-norm solver
returns a Moore-Penrose pseudoinverse solution on rank-deficient
designs that arise when the weight vector zeroes out entire
rows/columns of the design matrix. The SVD truncation tolerance
rcond is max(eps * max(m, n), 1e-12) — the 1e-12 floor
stabilises alpha-hat / beta-hat on the smallest benchmark panels
(Basque N=17, West Germany N=16) without perturbing tau-hat.
Locked in by tests/test_dgelsd_rank_deficient_wls.do.

3. UnitDistanceCache (rust/src/distance.rs). The pairwise
Sum (Y_i - Y_j)^2 quantities are precomputed once; each leave-t-out
distance dist_{-t}(j, i) is then an O(1) subtraction instead
of an O(T) rescan. Locked in by
tests/test_unit_distance_cache_equivalence.do.

4. Deterministic LOOCV tie-breaker (rust/src/loocv.rs,
better_candidate). When two (lambda_time, lambda_unit, lambda_nn)
triples score within TIE_TOL = 1e-10 of each other, trop prefers
the largest lambda_nn, then smallest lambda_time, then smallest
lambda_unit. The 1e-10 threshold sits roughly 4 orders of magnitude
above the IEEE 754 double-precision epsilon (~2.2e-16) yet well below
typical finite LOOCV scores (Q ~ 1e-3 on the standard benchmarks), so
the tie-breaker activates only on genuine numerical ties. Without it
ULP-level BLAS differences can flip argmin Q(lambda) across
platforms. Locked in by tests/test_loocv_tiebreak_determinism.do.

5. Inference reference distribution (ado/trop.ado,
mata/trop_ereturn_store.mata). The bootstrap resamples units in
stratified fashion (Algorithm 3 step 3), so the cluster count that
governs the small-sample reference df is N_1, the number of
ever-treated units — not the number of treated cells. trop
therefore uses t(N_1 - 1) whenever N_1 &gt;= 2 and falls back to
the standard normal otherwise. The primary confidence interval
defaults to the paper-specified percentile interval whenever bootstrap
is enabled; the cimethod() option re-selects the primary pair
from the three candidates (percentile, t, normal). Locked in by
tests/test_inference_df_is_treated_units.do and
tests/test_cimethod_option.do.

Scope exclusions (paper-adjacent features that trop does not
implement): time-varying covariates X_(t)*beta(t) where the coefficient
vector varies by period. The current covariates() option implements
the time-invariant form from paper Section 6.2 Equation 14 where gamma is a
fixed p-vector. Survey-weighted bootstrap with strata()/psu()/
fpc() is available since v1.2.0; switching-treatment patterns under
method(joint) remain out of scope.

References

Athey, S., G. W. Imbens, Z. Qu, and D. Viviano. 2025.
Triply robust panel estimators.
arXiv preprint arXiv:2508.21536.

Author

Xuanyu Cai

City University of Macau

xuanyuCAI@outlook.com

Wenli Xu

City University of Macau

wlxu@cityu.edu.mo

Installation

From GitHub:

. net install trop, from("https://raw.githubusercontent.com/gorgeousfish/TROP/main") replace

Local installation:

. net install trop, from("/path/to/trop_stata") replace

Verify installation:

. trop, version

. trop_check

The package includes precompiled plugins for macOS ARM64 (Apple Silicon),
macOS Intel (x86-64), and Windows x64. No Rust toolchain is required for
these platforms. See the project repository for build instructions if you
need to compile from source.

Also see

Online: trop, trop_estat, trop_predict