Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
142,145 characters · 75 sections · 92 citation commands
Big Wins, Small Net Gains: Direct and Spillover Effects of First Industry Entries in Puerto Rico
\thispagestyle{empty}
\pagestyle{plain}
\onehalfspacing
Large, persistent entries of new industries are central objects of local development policy. Governments devote substantial resources to attracting and retaining firms---through tax incentives, subsidies, and complementary business support programs---with the hope that major entries will generate jobs not only in the entering industry, but also in related sectors and the broader local economy TrinajsticEtAl2022_BusinessIncentives,SlatteryZidar2020_JEP. Evidence from regional firm subsidies, often targeted at manufacturing but with documented spillovers to other sectors and economically connected regions, shows that such programs can meaningfully raise local employment and reshape the composition of economic activity SieglochWehrhoeferEtzel2025_AEJPol,AobdiaKoesterPetacchi2024_JLE. Yet it remains unclear how these large entries reshape local labor markets over time, how far their effects travel across space, and whether they primarily benefit the entering industry or diffuse more broadly across the local economy.
Empirically, two challenges make these questions difficult to answer. First, much of the existing work on regional shocks and place-based policies relies on difference-in-differences designs with staggered treatment timing. A rapidly expanding econometric literature shows that commonly used implementations of these designs can perform poorly when treatment effects vary across units and over time, leading to misleading dynamic patterns and pre-trend diagnostics SunAbraham2021,RothSantAnnaBilinskiPoe2023,Baker2025. Evaluations of large industry entries that rely on staggered adoption must therefore use designs that remain valid under heterogeneous effects.
Second, standard evaluation frameworks typically assume no interference across units: one municipality's treatment is assumed not to affect another municipality's outcomes. This assumption is implausible in settings where firm entry in one location can reshape employment in nearby municipios or in related industries. Recent work develops approaches that explicitly allow for spillover effects and emphasize the need to separate own-treatment effects from spillovers on other units LiuHudgensSaulClemensAliEmch2019,Xu2023,FioriniLeePfeifer2024,ShahnZivichRenson2024. These developments are especially relevant when the policy motivation is explicitly about regional spillovers and cross-industry linkages.
This paper brings these insights to a new empirical setting: large, persistent industry entries in Puerto Rico, observed at the municipality--industry level using administrative employment records. I treat first sizable entry events---defined formally in Section (ref)---as shocks to local industry presence. The empirical framework separately tracks outcomes in the treated municipality--industry cell and in linked units connected through space and industry, with the central question of whether these “big wins” for the host municipio generate equally large net gains for the broader regional economy or mainly reallocate activity across space within the same industry. The analysis combines rich spatially disaggregated data with recent advances in policy evaluation designs that accommodate staggered treatment timing and spillovers across units.
The analysis is organized around the following pre-specified questions:
The remainder of the paper proceeds as follows. Section (ref) reviews the conceptual framework and related econometric literature on staggered DiD, interference-aware DiD, and spatially robust inference. Section (ref) describes the data, variables, spatial aggregation, and construction of heterogeneity strata. Section (ref) formalizes the treatment assignment and exposure mapping procedures. Section (ref) defines the estimands of interest, distinguishing direct effects from spillover channels. Section (ref) sets out the identifying assumptions, including conditional parallel trends under interference and exposure-positivity conditions. Section (ref) details the estimation procedures for direct and spillover effects and the auxiliary heterogeneity specification. Section (ref) presents the cluster- and spatially robust inference procedures. Section (ref) reports pre-trend tests, overlap and exposure-positivity diagnostics, spatial autocorrelation checks, and multiple-testing adjustments. Section (ref) presents the main results, including dynamic and cumulative effects for all exposure channels. Section (ref) explores heterogeneous effects across industry types and regions. Section (ref) discusses the economic interpretation and policy implications of the findings. Section (ref) reports sensitivity analyses for estimator choice, exposure-history definitions, and inference methods. Section (ref) concludes.
The empirical setting is a quarterly municipality--industry panel for Puerto Rico's 78 municipios, constructed from administrative employment and wage records (see Section (ref)). I treat first sizable entry events---defined formally in Section (ref)---as shocks to local industry presence.
Answering the research questions in Section (ref) requires methods that accommodate staggered adoption, heterogeneous treatment effects, interference across space and industries, and spatially correlated shocks. The remainder of this section situates the estimators and inference procedures used in the paper within three strands of the recent econometric literature: modern staggered DiD, interference-aware DiD and exposure mappings, and spatially robust inference. Formal estimands, identifying assumptions, and implementation details are developed in Sections (ref)--(ref).
My methodology is situated at the intersection of three recent advances in the econometrics literature: (1) heterogeneity-robust estimators for staggered difference-in-differences (DiD) designs; (2) DiD models that explicitly account for interference and spillovers; and (3) finite-sample-robust inference methods for spatial and panel data. This section provides a conceptual overview of these strands; the paper-specific implementation is developed in later sections.
The canonical two-way fixed effects (TWFE) regression, long the standard for DiD analysis, has been shown to be problematic in staggered adoption settings with heterogeneous treatment effects. When treatment timing varies, the TWFE estimator identifies a weighted average of individual treatment effects where some weights can be negative, because already-treated units are used as controls in “forbidden comparisons” RothSantAnnaBilinskiPoe2023,Baker2025. This can lead to biased estimates or even incorrectly signed coefficients.
Simulation studies confirm that while traditional TWFE performs well under homogeneous effects, its performance deteriorates in the presence of dynamic and heterogeneous effects WangHamadWhite2024. These findings motivate the use of “modern” DiD estimators designed for staggered adoption with heterogeneous treatment effects that avoid using already-treated units as controls.
For the primary analysis of direct effects, I use the estimator of CallawaySantAnna2021. This approach takes group-time average treatment effects $\mathrm{ATT}(g,t)$---the average effect for cohort $g$ at time $t$---as the central building block. By structuring the design around cohort-specific comparisons to clean control groups, the CS framework avoids the contamination issues that arise in conventional TWFE event studies.
As a cross-check and robustness exercise, I also employ the imputation-based estimator of BorusyakJaravelSpiess2024. This “BJS” method fits a TWFE model on clean (never- or not-yet-treated) observations and uses the estimated coefficients to impute counterfactual outcomes for treated observations. SunAbraham2021 emphasize that traditional event-study coefficients (TWFE with leads and lags) are contaminated by effects from other periods, so apparent pre-trends can arise solely from treatment effect heterogeneity. The CS and BJS estimators address this issue by relying only on clean control groups. Their implementation in my setting, including dynamic aggregation and event-time windows, is described in Section (ref).
A core contribution of the paper is to measure spillovers, so I build on the literature that relaxes the Stable Unit Treatment Value Assumption (SUTVA). The foundational concept is the exposure mapping, which formalizes interference by defining a unit's outcome as a function of the treatment status of its neighbors AronowEcklesSamiiZonszein2020. Spillover effects are defined as contrasts in expected outcomes across exposure levels rather than across treatment status alone Leung2024_ExposureContrasts.
In a DiD framework with interference, Xu2023 show that direct and spillover average treatment effects can be identified under a modified parallel trends assumption that holds conditional on exposure history. ShahnZivichRenson2024 connect this logic to the Structural Nested Mean Model (SNMM) literature and formalize how spillover effects can depend on both past and concurrent exposure, leading to a “network conditional parallel trends” assumption based on exposure history. I adopt this exposure-history perspective and specify a small number of interpretable exposure mappings that capture economically relevant channels of interference; these mappings are defined in Section (ref) and linked to estimands and identifying assumptions in Sections (ref)--(ref).
This approach follows the recommendation of Savje2023_MisspecExposure to separate the role of an exposure mapping---defining an interpretable estimand---from the impossible task of fully capturing the underlying causal structure. Related work includes FioriniLeePfeifer2024, who identify a direct effect under assumptions that limit spillovers, and PengYeZheng2025, who motivate neighborhood-structured interference via their Differences-in-Neighbors estimator. I draw on these contributions conceptually but employ a different implementation tailored to the Puerto Rico setting.
To conduct inference, I must account for spatial correlation in the municipality--industry panel. A standard approach is the Conley1999 spatial HAC estimator, which is a spatial analogue of time-series HAC methods ConleyMolinari2005. However, these and related methods can suffer from size distortions in non-stationary settings such as panel data and DiD designs Mueller_Watson_JBES_2023.
My primary inference strategy is the Conditional Spatial Correlation Principal Components (C-SCPC) method proposed by Mueller_Watson_JBES_2023, a robustified version of SCPC designed to have good size properties in empirically relevant panel settings. As an additional sensitivity check, I use the Thresholding Multiple Outcomes (TMO) procedure of DellaVignaImbensKimRitzwoller2025, which adjusts standard errors for cross-sectional dependence using auxiliary outcomes. The construction of spatial weight matrices, the implementation of C-SCPC, and the TMO adjustment are detailed in Section (ref).
This section describes the construction of the dataset used in the analysis, including the ingestion of raw employment files, the construction of a quarterly index, the harmonization of industry classifications, the creation of the municipios geometry layer, and the aggregation to a balanced municipality--industry--quarter panel.
The main data source is the Puerto Rico Department of Labor and Human Resources' Composición Industrial por Municipio portal DTRH_ComposicionIndustrial_XLSX. The tables are produced by the Programa Censo Trimestral de Empleo y Salarios (Puerto Rico's implementation of the U.S. Quarterly Census of Employment and Wages, QCEW; DTRH_QCEW_nd) and compile administrative records from the unemployment insurance system and quarterly employer tax reports, covering nearly all wage-and-salary employment in Puerto Rico.
Each quarterly release reports, for every municipality and NAICS industry, the number of establishments, total employment, total wages, and average wages. The dataset excludes self-employed and unpaid workers, but includes full- and part-time payroll employees, corporate officers, and federal personnel assigned to Puerto Rico. Establishments are assigned to municipios by physical location, with special codes identifying multi-municipality aggregates and unclassified sites. Detailed variable definitions, NAICS versioning, and coverage notes follow the official quarterly publication DTRH_ComposicionIndustrial_2025Q1 and the program's methodological report DTRH_QCEW_nd. All files are processed using a consistent schema, cleaned, and stacked into a single quarterly dataset.
I construct a complete quarterly time index spanning from 2014Q1 to 2025Q1. The index contains the calendar year, quarter number, a label (e.g., “2018Q3”), the start date of the quarter, and a sequential integer variable $t$. These identifiers provide a harmonized timeline used for event-time alignment and temporal joins across sources.
Industry codes appear in different formats and levels of aggregation across files. I harmonize all NAICS labels by normalizing the underlying text and standardizing numeric strings, and then merge the cleaned records against a manually curated NAICS mapping table. This table specifies both exact codes and two-digit ranges, with range entries expanded to cover the full set of corresponding sectors. When English titles are present for a mapped range, they are propagated to all sectors covered by that range. The resulting variable naics_harmonized combines observed and mapped identifiers into a single consistent classification used throughout the analysis.
I build a municipios geometry layer for Puerto Rico by filtering TIGER/Line county shapefiles to the state code STATEFP=72, ensuring one observation per municipality. I standardize identifiers (STATEFP, COUNTYFP, GEOID) and rename the municipality name field to NAME, trimming whitespace while preserving official casing. I de-duplicate on GEOID, sort by \texttt{GEOID} to stabilize joins, and assert exactly 78 municipios post-processing.
The canonical layer is saved as a GeoPackage (layer name municipios), with a companion table recording attribute metadata and the coordinate reference system used. To compute area and perimeter diagnostics and support distance calculations, I project geometries to the Puerto Rico StatePlane coordinate system (EPSG:32161).\footnote{epsg32161 defines EPSG:32161---NAD83 / Puerto Rico & Virgin Islands, units in meters.} This geometry layer is later used to construct spatial weights for the exposure mappings and inference procedures described in Sections (ref) and (ref). Geometry validation and repair procedures follow standard OGC and GEOS practices.\footnote{Administrative boundary semantics follow census_tiger_county, which defines STATEFP, COUNTYFP, and GEOID for county-level shapefiles.}
I aggregate the harmonized industry data to the municipality--industry--quarter level, denoted $(\text{GEOID}, \text{naics\_code}, \text{quarter})$, where GEOID follows the 5-digit FIPS format (72xxx). Before aggregation, I apply exclusion filters to remove multi-municipality aggregates and unclassified observations. Specifically, I drop municipality codes 995 and 999 (which represent “Multi Municipio” and “No Codificado” records) and municipality names matching these labels, and exclude observations with NAICS codes 99 or 10 (unclassified establishments). Municipality codes are standardized by extracting the last three digits from the original code field and prepending “72” to construct valid 5-digit FIPS identifiers, ensuring consistency with the canonical geometry layer.
For each $(\text{GEOID}, \text{naics\_code}, \text{quarter})$ triplet, I compute four core outcome variables by summing over all reporting units: establishments (count of establishments), covered_emp (total covered employment), total_wages (aggregate payroll), and avg_wage_level (defined as total wages divided by covered employment when the latter is positive, and missing otherwise).
To construct real wage measures, I load the monthly CPI series (Concatenación) from the Puerto Rico Department of Labor and Human Resources DTRH_IPC_Concatenacion_XLS. The CPI corresponds to the official Índice de Precios al Consumidor en Puerto Rico DTRH_IPC_Sep2025. I compute quarterly means by averaging the three months in each quarter and rebase the index so that the 2020 average equals 100. Nominal wages are then deflated to obtain total_wages_real_2020 and avg_wage_level_real_2020. All outcome variables (establishments, employment, nominal wages, and real wages) are log-transformed where positive, yielding log_establishments, log_covered_emp, \texttt{log_total_wages}, and \texttt{log_total_wages_real_2020}.
To ensure a balanced time structure, I construct a complete skeleton by crossing all observed $(\text{GEOID}, \text{naics\_code})$ pairs with the full quarterly index from Section (ref). This guarantees that every municipality--industry pair that appears in the data has a record for every quarter from 2014Q1 to 2025Q1, with missing observations explicitly represented as NA. A stable industry label (industry_name) is retained for each naics_code by selecting the most frequently observed English industry name, which is used in figures and tables. The resulting panel is the main analysis dataset.
The heterogeneity analysis relies on time-invariant strata defined at the municipality--industry level $(\text{GEOID}, \text{naics\_code})$. I construct three such strata: (i) tradable vs.\ non-tradable sectors, (ii) metro vs.\ non-metro municipios, and (iii) high- vs.\ low-wage sectors. These labels are computed once from the pre-period data and then merged onto the full panel.
\paragraph{Tradable vs.\ non-tradable sectors.}
Using the harmonized NAICS codes from Section (ref), I map each industry to its 2-digit NAICS sector and classify sectors as tradable or non-tradable. Sectors 11, 21, 31--33, 42, 48--49, 51, 52, 54, and 55 are treated as tradable, reflecting goods and services that can be supplied across regions. Sectors 22, 23, 44--45, 53, 56, 61, 62, 71, 72, 81, 92, and 99 are treated as non-tradable, reflecting local services and public administration. Any remaining 2-digit sectors not appearing in this list are conservatively assigned to the non-tradable group. The resulting tradable indicator is constant over time and across municipios for a given industry.
\paragraph{Metro vs.\ non-metro municipios.}
I classify municipios as metro or non-metro using Core Based Statistical Area delineations from OMB Bulletin No. 23-01 OMB_Bulletin_23_01. This table contains a binary metro indicator, which I merge onto the panel by municipality identifier GEOID. Municipios lacking a valid classification are treated as non-metro. The metro status is thus time-invariant and shared by all industries within a municipio.
\paragraph{High- vs.\ low-wage sectors.}
Finally, I construct wage-based strata using the real average wage level avg_wage_level_real_2020 defined in Section (ref). I restrict attention to the pre-program period 2014Q1--2019Q3 and, for each industry $k$, compute the median of avg_wage_level_real_2020 across all municipios and quarters in this window. Let $\tilde{w}_k$ denote this industry-level pre-period median and let $\tilde{w}$ denote the median of $\{\tilde{w}_k\}$ across all industries with sufficient data. Industries with $\tilde{w}_k \ge \tilde{w}$ are classified as high-wage and the remainder as low-wage. Industries with insufficient pre-period wage data are assigned to a residual “unknown” group, which I exclude from the high vs.\ low comparisons in the heterogeneity analysis.
For the Thresholding Multiple Outcomes (TMO) sensitivity analysis of DellaVignaImbensKimRitzwoller2025, I construct a separate municipal-level panel using tables from the ACS 5-year estimates for 2013--2023 Census_ACS_5yr. These tables provide socioeconomic metrics on demographics, employment, and housing. I build a set of 44 indicators, compute year-over-year changes, and upsample the resulting annual series to the quarterly frequency of the main analysis. Because the ACS 5-year series is currently available only through 2023, I carry the 2023 year-over-year changes in each auxiliary outcome forward to 2024--2025 so that the quarterly auxiliary panel spans the same time horizon as the main employment data. This auxiliary panel is used solely to estimate the cross-location correlation matrix that enters the TMO adjustment in the spatial inference procedures described in Section (ref).
Table (ref) summarizes the pre-program distribution of outcomes at the municipality--industry--quarter level over 2014Q1--2019Q3. Panel A reports means and standard deviations of the logarithms of covered employment, establishments, and the real wage bill. Municipality--industry cells are typically small but exhibit substantial dispersion, with a heavy upper tail of large sectors. Panel B shows that a large majority of municipality--industry cells belong to tradable sectors, most are located in metro municipios, and high- and low-wage industries are roughly evenly split, indicating that the balanced panel covers a broad cross-section of Puerto Rico's local production structure.
Figure (ref) maps the municipios geometry and highlights the metro versus non-metro split used in the heterogeneity analysis. Metro municipios, as defined by OMB Bulletin 23-01, cluster along the north coast and southern corridor, while non-metro municipios are concentrated in the central mountainous region and offshore islands.
Before defining exposure mappings, I first define the direct treatment event $G_{ik}$. The unit of analysis is a municipality--industry pair $(i,k)$, and “treatment” is a large entry event derived from establishment counts, employment, and real wage-bill data (total_wages_real_2020; Section (ref)). In the baseline, the event is the first quarter $t$ in which $(i,k)$ satisfies any of three establishment-based persistent trigger conditions (T1--T3).
When multiple rules are met in the same quarter, I apply a fixed priority order $\mathrm{T1}>\mathrm{T2}>\mathrm{T3}$. The cohort identifier $G_{ik}$ is defined as the earliest quarter that satisfies any of T1--T3. The persistence requirements distinguish durable entries from temporary fluctuations. I interpret these large, persistent jumps as first sizable entries that proxy for major expansions of industry presence in a municipality. The cohort variable $G_{ik}$ underpins the staggered-adoption design used throughout the paper; formal estimands based on $G_{ik}$ are defined in Section (ref), and event-time windows are discussed in Section (ref).
The spillover analysis is based on explicit exposure mappings that translate neighbors' treatment into exposure variables for each unit. I adopt the exposure-mapping framework of AronowEcklesSamiiZonszein2020, in which a unit's exposure is a function of the treatment status of its neighbors and a specified network structure, and the exposure-contrast perspective of Leung2024_ExposureContrasts. The conceptual discussion and identification logic for interference-aware DiD with exposure histories follow Xu2023 and ShahnZivichRenson2024 and are summarized in Section (ref). Here I define the specific exposure channels used in the analysis.
I specify three distinct exposure channels for municipality--industry cells $(i,k)$:
These mappings are designed to distinguish distinct economic channels. By defining separate exposure histories for same-industry versus cross-industry neighbors, I allow for the possibility that competitive effects within the same industry SlatteryZidar2020_JEP and agglomeration or demand spillovers across industries SieglochWehrhoeferEtzel2025_AEJPol operate simultaneously but with different signs and magnitudes. Formal estimands based on contrasts under each channel are introduced in Section (ref).
The spatial weight matrix $W$ is constructed from the municipios geometry layer for Puerto Rico (Section (ref)). I start from a Queen contiguity representation based on the municipios polygons and, for municipios that are isolates under pure contiguity, add $k=3$ centroid-based nearest-neighbor (KNN) links. These additional edges are added only for isolates, following a “Queen + KNN plug-in” rule.
After adding these links, the neighbor graph is made symmetric (if $i$ is a neighbor of $j$, then $j$ is a neighbor of $i$), and the resulting weight matrix is then row-standardized. Row-standardization implies that each row of $W$ sums to one, so exposure measures formed as weighted averages of neighbors' treatment can be interpreted as mean treatment among a unit's neighbors. The underlying neighbor relation is undirected, even though the row-standardized matrix need not be numerically symmetric. These weights serve as the basis for the same-industry and all-industries neighbor exposure measures described above.
In addition to contemporaneous exposure, I summarize each channel's exposure history using current exposure shares and short lagged histories. For the direct own treatment and each spillover channel (same-industry neighbor exposure, within-municipality cross-industry exposure, and neighbor all-industries exposure), I construct two such summaries at each quarter $t$:
For each channel, I also include up to four lags of these exposure shares as additional history controls in the nuisance functions.
These summaries are defined analogously across channels, with the underlying exposure at each date determined by the relevant mapping ($\mathrm{Expo}$, $\mathrm{XExpo}$, or $\mathrm{NExpo}$). The resulting history variables are used both to define estimands that depend on exposure duration and to condition on exposure history in the identification and estimation procedures. Their role in the identifying assumptions and doubly-robust estimation strategy is detailed in Sections (ref) and (ref).
The treatment rule in Section (ref) yields 1{,}871 municipality--industry cells with a first sizable entry over 2014Q2--2024Q4. Figure (ref) plots the distribution of cohort dates $G_{ik}$ by calendar quarter. Entries are heavily front-loaded: the bulk of first sizable entries occur between 2014 and 2016, after which the number of new treated municipality--industry cells declines steadily and falls to relatively low levels by the early 2020s. There are no first sizable entries in the initial quarter 2014Q1 or after 2024Q4.
Figure (ref) maps the cumulative number of first sizable entries by municipio. Every municipio experiences at least ten entry events, but the intensity of entry is unevenly distributed across space: municipios receive between 10 and 40 first sizable entries, with a median of 24 and an interquartile range of 20 to 27 entries. The map bins municipios into five quintile groups based on their entry counts (10--18, 19--22, 23--26, 27--28, and 29--40 entries), highlighting the concentration of entering industries in coastal and metropolitan municipios and the relatively lower intensity in many interior municipios.
Tables (ref) and (ref) decompose treated municipality--industry cells by trigger type. About three quarters of first sizable entries arise from large, top-decile establishment jumps (T3), while the remainder are split between “new” entries (T1) from previously empty cells and “small-to-large” expansions (T2) that move from at most one establishment to at least two. Roughly one third of treated cells under each trigger contribute to the balanced $[-8,+16]$ event-time window used in the main event-study analysis, and the median pre-period employment levels differ sharply across trigger types.
Figure (ref) provides a simple check on the treatment definition by plotting average establishment counts around the entry event for each trigger type, using only municipality--industry cells that belong to the balanced window. T1 (new entries) and T2 (small-to-large growth) show discrete increases in establishments at $t=0$ from very low pre-period levels, while T3 (top-decile jumps) exhibits a large, persistent increase from roughly 26 establishments pre-event to around 80 post-event, consistent with the interpretation of these triggers as sizable, durable entry events.
I next formalize the treatment effect parameters that correspond to the research questions in Section (ref). For each municipality--industry cell $(i,k)$, let $G_{ik}$ denote the first quarter of a sizable entry (Section (ref)) and let $\ell = t - G_{ik}$ be time since entry. The primary direct-effect target is the dynamic path $\mathrm{ATT}^{\mathrm{own}}(\ell)$, which measures the average effect of a first sizable entry on the three main log outcomes---$\log\_\mathrm{establishments}$, $\log\_\mathrm{covered\_emp}$, and the real wage bill $\log\_\mathrm{total\_wages\_real\_2020}$ (2020 = 100)---in the treated municipality--industry cell at horizon $\ell$ (Q1). I also define slice-level direct and spillover effects that summarize the average impact over post-treatment windows such as 0--4, 5--8, 9--16, and 0--16 quarters, corresponding to the “policy slices” used to summarize cumulative gains (Q5).
Following CallawaySantAnna2021, I use group-time average treatment effects $\mathrm{ATT}(g,t)$ as the basic building blocks and aggregate these effects across cohorts sharing the same relative event time $\ell = t - g$ to obtain $\mathrm{ATT}^{\mathrm{own}}(\ell)$. I then extend this framework to spillover channels by defining exposure-based slice estimands that hold own treatment fixed while varying exposure to neighbors' entries along three pre-specified channels: same-industry neighbors, within-municipality cross-industry, and neighbor all-industries (Q2--Q3b). Unless stated otherwise, dynamic and slice-level parameters are defined with unit weights over municipality--industry cells, so the target population is the average treated municipality--industry cell rather than the average worker.
In the context of Q1, for a given outcome $Y$ (log covered employment, log establishments, or log real total wages), $\mathrm{ATT}^{\mathrm{own}}(\ell)$ can be read as the average percentage effect on $Y$ in a municipality--industry cell that experienced a first sizable entry $\ell$ quarters ago, relative to what would have happened in the absence of such an entry. The cumulative quantity \[ \mathrm{ATT}^{\mathrm{own}}_{\mathrm{cumu}}(0,16) = \sum_{\ell=0}^{16}\mathrm{ATT}^{\mathrm{own}}(\ell) \] summarizes the total impact over the first 16 post-entry quarters and is the main direct-effect object used in the budgeting and return-on-investment discussion in Q5.
For notational simplicity, I suppress the industry index $k$ and write $i$ for the municipality--industry cell $(i,k)$. Let $Y_{i,t}(1)$ and $Y_{i,t}(0)$ denote the potential outcomes for unit $i$ at time $t$ with and without a first sizable entry at $G_i=g$. The primary estimand of interest is the direct average treatment effect on the treated (ATT) at event-time horizon $\ell$, denoted $\mathrm{ATT}^{\mathrm{own}}(\ell)$, defined as \[ \mathrm{ATT}^{\mathrm{own}}(\ell) \equiv \mathbb{E}\!\left[ Y_{i,g+\ell}(1) - Y_{i,g+\ell}(0) ~\middle|~ G_i=g \right], \] with aggregation across cohorts sharing the same $\ell$ implemented via the group-time ATT structure of CallawaySantAnna2021. The cumulative $\mathrm{ATT}^{\mathrm{own}}_{\mathrm{cumu}}(0,16)$ can be viewed as a summary of these horizon-specific effects over the post-entry window of interest. Practical implementation choices, including reference periods and limited-anticipation windows, are described in Sections (ref) and (ref).
The spillover estimands answer Q2--Q3b. For a given slice $[a,b]$, $\mathrm{DATT}[a,b]$ captures the average direct effect of own entry over horizons $\ell\in[a,b]$ on the treated municipality--industry cell, while $\mathrm{SATT}_{\mathrm{same}}[a,b]$, $\mathrm{SATT}_{\mathrm{cross}}[a,b]$, and $\mathrm{SATT}_{\mathrm{nall}}[a,b]$ capture the average effects of increasing exposure along the same-industry neighbor, within-municipality cross-industry, and neighbor all-industries channels, respectively, holding own treatment fixed. Thus:
The aggregate spillover $\mathrm{SATT}[a,b]$ and total effect $\mathrm{TATT}[a,b]$ then decompose the post-entry impact into direct and spillover components that feed into the cumulative-value calculations in Q5.
I implement these objects using the exposure-mapping framework described in Sections (ref) and (ref), summarizing neighbors' assignments by low-dimensional exposure variables and their histories and defining spillover estimands as exposure contrasts that hold a unit's own treatment fixed.
\paragraph{Slices and notation.}
Let the event-time slices be $[a,b]\in\{[0,4],[5,8],[9,16],[0,16]\}$ and denote the slice length by $L_s = b-a+1$. For each slice, I define an own-treatment indicator \[ D_{\mathrm{slice}}=\mathbf{1}\{t-g\in[a,b],\,A_{i,t}=1\}, \] where $A_{i,t}$ is the indicator of own treatment for unit $i$ at time $t$. For the spillover channels, I summarize exposure over the slice by summing the underlying per-period exposure processes. For example, with $E^{\mathrm{same}}_{i,t}(\ell)$ denoting the same-industry neighbor exposure at event time $\ell$, \[ S_{\mathrm{same},\mathrm{slice}} = \sum_{\ell\in[a,b]} E^{\mathrm{same}}_{i,t}(\ell), \] and analogous definitions apply to $S_{\mathrm{cross},\mathrm{slice}}$ and $S_{\mathrm{nall},\mathrm{slice}}$, pooling across source industries where appropriate (consistent with the exposure construction in Section (ref)). Because each $E^{\cdot}_{i,t}(\ell)$ is a share in $[0,1]$, the slice sums $S_{\cdot,\mathrm{slice}}$ lie in $[0,L_s]$. A one-unit increase in, for example, $S_{\mathrm{same},\mathrm{slice}}$ corresponds to an additional quarter of full exposure; moving from zero exposure to 100% exposure in every quarter of the slice corresponds to an $L_s$-unit increase, so the effect of that full-slice contrast is $L_s$ times the corresponding slice coefficient. For overlap diagnostics and the $\epsilon$-trimming rule in Section (ref), I also work with slice-average shares $\bar S_{\cdot,\mathrm{slice}} = S_{\cdot,\mathrm{slice}}/L_s \in [0,1]$. In the DR--DiD specification in Section (ref), the slice sums $S_{\cdot,\mathrm{slice}}$ serve as the main spillover regressors, while the averages $\bar S_{\cdot,\mathrm{slice}}$ are used only in overlap diagnostics.
\paragraph{Slice-level spillover estimands.}
I define direct and channel-specific spillover targets as average partial effects within a slice, holding the other channels fixed. Let \[ Y_{i,t}(d,s_{\mathrm{same}},s_{\mathrm{cross}},s_{\mathrm{nall}}) \] denote the potential outcome for unit $i$ at time $t$ under own-treatment status $d$ and exposure slices $(s_{\mathrm{same}},s_{\mathrm{cross}},s_{\mathrm{nall}})$. Then
I summarize total spillovers by \[ \mathrm{SATT}[a,b]=\mathrm{SATT}_{\mathrm{same}}[a,b]+\mathrm{SATT}_{\mathrm{cross}}[a,b]+\mathrm{SATT}_{\mathrm{nall}}[a,b], \] and the total effect by \[ \mathrm{TATT}[a,b]=\mathrm{DATT}[a,b]+\mathrm{SATT}[a,b]. \]
The assumptions under which these slice-level DATT and SATT objects admit a causal interpretation are discussed in Section (ref), and the estimation strategy that delivers consistent estimators of these parameters is presented in Section (ref).
I follow Xu2023 in adopting a parallel trends under interference, conditional on exposure histories assumption as the core identifying restriction. Intuitively, conditional on past exposure histories and observed controls, potential outcome trends would have evolved in parallel across exposure groups in the absence of treatment, generalizing the classical DiD assumption to settings with interference (see also ShahnZivichRenson2024).
This formulation nests the standard staggered DiD setting of CallawaySantAnna2021 as a special case with no interference, in which conditioning is on pre-treatment covariates rather than on exposure paths. When the exposure mapping is misspecified, I follow Xu2023 and Savje2023_MisspecExposure in interpreting the direct path as an expected direct ATT (EDATT) and the spillover targets as expected exposure contrasts with respect to the chosen mappings. In practice, the exposure histories that index this assumption enter the estimator via the low-dimensional history summaries described in Section (ref); implementation using doubly robust DiD and cross-fitting is discussed in Section (ref).
I assume a limited anticipation window of $\delta$ quarters (baseline: $\delta=2$, robustness: $\delta=1$). Economically, this allows for short-run adjustments in the last few pre-treatment quarters but rules out systematic anticipation far in advance of a first sizable entry. Following CallawaySantAnna2021, I operationalize this by using a cohort-specific pre-treatment period, separated from the treatment date by a buffer of $\delta$ quarters, as the reference outcome.
For a unit first treated in quarter $g$, event-time contrasts at horizon $\ell$ compare outcomes $\ell$ periods after adoption to outcomes in this cohort-specific reference period. Anticipatory leads $\ell \in \{-\delta,\ldots,-1\}$ are treated as diagnostics for pre-trends rather than as identifying moments for the dynamic effect path. The choice of reference period and treatment of leads in estimation and inference are detailed in Sections (ref) and (ref).
I assume that interference operates only through the three pre-specified exposure channels defined in Section (ref): same-industry neighbors, within-municipality cross-industry, and neighbor all-industries. Exposure is constructed on the queen-contiguity plus KNN graph with a single-hop neighborhood ($h=1$) derived from the spatial weights matrix $W$ in Section (ref).
This corresponds to a partial interference structure in the sense of LiuHudgensSaulClemensAliEmch2019, in which spillovers are confined to predefined networks rather than allowed to propagate arbitrarily across the economy. Multi-hop exposure measures are used for diagnostics but do not define the primary estimands.
Identification additionally requires exposure positivity: for each exposure contrast of interest, there must be sufficient variation in exposure levels among otherwise comparable units so that counterfactual outcomes are well defined and estimable. I follow Xu2023 and ShahnZivichRenson2024 in treating overlap in exposure histories as necessary for DiD with interference, paralleling the support requirements in AronowEcklesSamiiZonszein2020 and LiuHudgensSaulClemensAliEmch2019.
In practice, I enforce exposure positivity by trimming observations in regions of very low or very high exposure where empirical support is weak. Specifically, I bound the slice-average exposure shares $\bar S_{\cdot,\mathrm{slice}}$ defined in Section (ref) away from zero and one, retaining only units whose exposure lies in an interior interval $[\epsilon,1-\epsilon]$ for each channel and slice, with a small $\epsilon>0$. This trimming rule is motivated by the overlap conditions in crump2009 and the exposure-positivity discussion in Savje2023_MisspecExposure and Xu2023. The exact choice of $\epsilon$ and the associated minimum-sample thresholds, as well as graphical diagnostics for exposure overlap, are described in Section (ref).
The assumptions above provide a baseline causal interpretation for the direct and spillover estimands but may be violated in practice. To assess robustness to deviations from parallel trends, I use the Honest DiD framework of RambachanRoth2023, which constructs robust confidence sets under explicit restrictions on how much the path of treatment effects can deviate from parallel trends in pre-treatment periods.
I focus on the smoothness restriction $\mathcal{S}\mathcal{D}(M)$, which bounds the curvature (second differences) of the effect path and is well suited to settings with potentially accelerating dynamic responses. This choice aligns with recommendations in recent DiD syntheses RothSantAnnaBilinskiPoe2023 and applied evaluations WangHamadWhite2024. Following Xu2023, I interpret these sensitivity analyses in the presence of interference as robustness checks for both the parallel trends and exposure-mapping assumptions. Implementation details and the presentation of robust bounds are discussed in Section (ref) and revisited in the robustness analysis in Section (ref).
For clarity, I summarize the main assumptions underlying the identification of direct and spillover effects; each corresponds to the subsections above.
Throughout, I treat log_covered_emp as the primary outcome. Complementary outcomes are log_establishments and the log real wage bill, log_total_wages_real_2020. Nominal wage measures, including log_total_wages, are used only for internal checks and are not reported as main outcomes. All specifications include municipality--industry and quarter fixed effects, implemented via within-demeaning.
I distinguish between balanced and unbalanced event-time horizons following CallawaySantAnna2021. A balanced horizon restricts attention to event times that are observed for all treated cohorts, so each relative event time $\ell$ is estimated from the same set of treated units. BorusyakJaravelSpiess2024 emphasize that such balanced aggregation improves efficiency and avoids composition changes when treatment timing is staggered, while DubeGirardiJordaTaylor2025 show that unbalanced horizons can distort dynamic patterns. I therefore adopt the balanced horizon as the primary window and use unbalanced paths as diagnostic overlays, following the recommendation of Baker2025 to guard against imbalance in event time when aggregating treatment effects.
In practice, treated municipality--industry cells must have sufficient pre- and post-treatment observations to contribute to the full event-time window, while all available control observations are retained. Implementation details are documented in Section (ref).
For the primary analysis of direct effects, I estimate the dynamic path $\mathrm{ATT}^{\mathrm{own}}(\ell)$ defined in Section (ref) using the group--time ATT estimator of CallawaySantAnna2021. The comparison set consists of municipality--industry cells that are untreated at time $t$, combining not-yet-treated and never-treated units within the same industry (an NY design). I use the unconditional version of the estimator (no additional covariates) and impose a limited-anticipation window by differencing outcomes relative to a cohort-specific pre-event baseline quarter $g-(\delta+1)$, as described in Section (ref).
The balanced event-time horizon from the previous subsection is used as the primary window, with unbalanced paths reported in figures for comparison and composition diagnostics. Uncertainty is summarized using uniform (simultaneous) confidence bands based on multiplier bootstrap procedures; inference details are provided in Section (ref).
As a complementary check, I implement the imputation-based estimator of BorusyakJaravelSpiess2024. I first fit a two-way fixed-effects model (unit and time effects) to untreated observations only (within industry), use the fitted model to impute untreated potential outcomes $\widehat{Y}_{it}(0)$ for all cells, and define estimated treatment effects on treated cells as $Y_{it} - \widehat{Y}_{it}(0)$. Dynamic event-time effects $\mathrm{ATT}^{\mathrm{own}}(\ell)$ are then constructed by averaging these imputed treatment effects across treated municipality--industry cells at each event time $\ell$, with equal weight on each treated unit.
I estimate these effects on the same balanced event-time horizon as in the CS specification and report unbalanced paths as robustness overlays. Standard errors and simultaneous confidence bands for the BJS-based event-time path are obtained using GEOID-cluster bootstrap procedures (Section (ref)). Cumulative effects over $\ell \in [0,16]$ are computed by summing the dynamic coefficients, and the same bootstrap draws are used to obtain uncertainty for these cumulative statistics.
Spillover effects are estimated using a doubly robust (DR) DiD estimator with cross-fitting, following Xu2023 and ShahnZivichRenson2024. The analysis is organized in terms of event-time slices $[a,b] \in \{[0,4],[5,8],[9,16],[0,16]\}$, corresponding to the policy windows and slice-level estimands in Section (ref). For each slice, I use the own-treatment indicator $D_{\text{slice}}$ and the three exposure-slice variables $S_{\text{same},\text{slice}}$, $S_{\text{cross},\text{slice}}$, and $S_{\text{nall},\text{slice}}$ as defined in Section (ref). The interpretation and scaling of these exposure-slice variables, and their link to $\mathrm{DATT}[a,b]$ and the channel-specific $\mathrm{SATT}$ parameters, follow Section (ref).
For each outcome and slice, I estimate a linear model in which both the outcome and the regressors $(D_{\text{slice}}, S_{\text{same},\text{slice}}, S_{\text{cross},\text{slice}}, S_{\text{nall},\text{slice}})$ are first residualized using cross-fitted nuisance models (Section (ref)), and then regressed in a single OLS step with GEOID-clustered inference. The resulting coefficients provide estimates of $\mathrm{DATT}[a,b]$ and the channel-specific spillover parameters.
Balanced horizons are used as the primary specification: treated or exposed units must have sufficient support over the full post-treatment window, while control observations are retained. Before estimation, I enforce the overlap conditions described in Section (ref) by trimming observations with extreme exposure shares and imposing minimum post-trim sample sizes per slice; the trimming rule and diagnostics are detailed in Section (ref). Unbalanced versions of the slice estimates are reported for comparison and as additional diagnostics.
Cross-fitting is used to obtain doubly robust estimators while respecting the clustering structure of the data. I implement cross-fitting with group-based folds defined at the municipality level (GroupKFold by GEOID), so that all observations from a given municipality are assigned to the same fold. This strategy is supported by Cao2024_ClusterCV, who show that cluster-respecting cross-validation is important for dependent data and that buffer zones are not required when using asymptotically linear learners.
I employ simple parametric learners as nuisance models. The outcome regression is estimated via ridge regression, the own-treatment propensity via logistic regression, and the exposure-slice nuisance functions via ridge regression. Exposure-history controls, constructed as described in Section (ref), enter the nuisance models to implement the “parallel trends conditional on exposure histories” assumption. Positivity and availability flags are enforced as sample filters but are not used as features in the nuisance models. After cross-fitted residualization of the outcome and regressors, I apply a two-way within transformation at the municipality--industry and quarter levels and estimate the final slice regression with GEOID-clustered standard errors.
The heterogeneity analysis is based on an auxiliary two-way fixed-effects specification that uses the strata defined in Section (ref) but does not re-estimate the full DR--DiD model for each subgroup. For each outcome $Y_{ikt}$ and event-time slice $[a,b]\in\{[0,4],[5,8],[9,16],[0,16]\}$, I construct the same direct-treatment and same-industry spillover regressors as in the main DR--DiD specification, namely $D_{\text{slice}}$ and a same-industry exposure slice $S_{\text{slice}}$ based on the same-industry neighbor shares.
I then estimate a linear model with municipality--industry and quarter fixed effects using within-demeaning. For each stratum variable $Z\in\{\text{tradable},\text{metro},\text{wage\_stratum}\}$, I interact both $D_{\text{slice}}$ and $S_{\text{slice}}$ with indicator functions for the strata in $Z$ and read off stratum-specific direct and spillover coefficients from these interactions. The heterogeneity analysis focuses on the direct and same-industry neighbor spillover channels; the within-municipality cross-industry and neighbor all-industries spillovers are not re-estimated by stratum.
The estimation sample for the heterogeneity specification is restricted to the balanced event-time window so that each slice is estimated on a common set of units. There is no additional trimming or cross-fitting in this auxiliary model. When available, exposure-history controls enter additively to support the same “parallel trends conditional on exposure histories” restriction, but these controls are not interacted with the strata. Standard errors are computed using cluster-robust inference at the municipality level, and I form stratum-specific DATT, SATT, and TATT estimates for each outcome--slice--stratum combination, as well as aggregate summaries that average across strata using the distribution of municipality--industry units as weights. Adjustments for multiple testing across heterogeneous effects are discussed in Section (ref).
All regression specifications include municipality--industry and quarter fixed effects, as described in Section (ref). Inference combines cluster-robust standard errors at the level of treatment assignment with event-study bootstrap procedures and, when warranted, spatially robust covariance estimators.
For the direct-effect event studies, I construct uniform (simultaneous) confidence bands for the dynamic path $\mathrm{ATT}^{\mathrm{own}}(\ell)$ under both the Callaway--Sant'Anna (CS) and Borusyak--Jaravel--Spiess (BJS) estimators.
For the CS estimator, I follow CallawaySantAnna2021 and use a Rademacher multiplier bootstrap applied to the influence-function-based aggregation over group--time cells. Each bootstrap draw perturbs the group--time contributions with independent Rademacher multipliers drawn at the municipality level, recomputes the aggregated $\mathrm{ATT}^{\mathrm{own}}(\ell)$, and forms a max-$|t|$ statistic over all event times in the estimation window, excluding the omitted reference bin $\ell=-1$ and the anticipatory leads $\ell\in\{-\delta,\ldots,-1\}$ implied by the limited-anticipation window. The $(1-\alpha)$ quantile of this max-$|t|$ distribution is used as a critical value to construct uniform bands for all reported event times.
For the BJS estimator, I obtain uncertainty from a GEOID-cluster bootstrap applied to the municipality-level effect series. Each bootstrap sample resamples municipalities with replacement, recomputes the cluster-mean $\mathrm{ATT}^{\mathrm{own}}(\ell)$ path, and forms a max-$|t|$ statistic over all event times in the estimation window, again excluding the omitted reference bin and anticipatory leads. Its $(1-\alpha)$ quantile is then used to construct uniform bands for the BJS-based event-time profile.
Cumulative effects over $\ell\in[0,16]$ are computed by summing the dynamic coefficients, \[ \widehat{\mathrm{ATT}}^{\mathrm{own}}_{\mathrm{cumu}}(0,16) = \sum_{\ell=0}^{16}\widehat{\mathrm{ATT}}^{\mathrm{own}}(\ell). \] For the BJS estimator, their uncertainty is obtained from the same GEOID-cluster bootstrap draws used for the dynamic bands: for each draw, I sum the corresponding event-time coefficients and report percentile intervals from the bootstrap distribution of the cumulative effect. For the CS estimator, I focus inference on the dynamic path and its uniform bands, treating the corresponding cumulative sums as descriptive summaries.
For the slice-level DR--DiD spillover estimands (Section (ref)), I use GEOID-cluster-robust standard errors as the baseline and report spatially robust alternatives described in the subsections below.
I report standard errors clustered at the municipality level---the level of independent treatment assignment---following design-based reasoning in RothSantAnnaBilinskiPoe2023. They recommend clustering at the level where treatment is independently assigned, a consideration that is particularly important when the number of clusters is not large. This clustering choice serves as the baseline for uncertainty quantification for the direct and spillover effects in Q1--Q4, as well as the cumulative-value objects in Q5.
To determine whether spatially robust inference is warranted, I apply a residual spatial autocorrelation diagnostic based on the global Moran's $I$ statistic. As emphasized by Chen2016, the Durbin--Watson test is inappropriate for spatial cross-sectional data because its value depends on the ordering of observations, whereas Moran's $I$ provides a permutation-invariant measure of spatial autocorrelation.
I compute Moran's $I$ on the regression residuals using the spatial-weights matrix $W$ from Section (ref) and obtain $p$-values from a permutation reference distribution. I define the trigger for spatial inference as a significant Moran's $I$ result ($p < 0.05$) in this permutation test.
For multivariate settings (e.g., stacked residuals across multiple outcomes), I use the multivariate extension of Moran's $I$ proposed by Yamada2024, which summarizes spatial autocorrelation across multiple residual series as a single statistic. Whenever this diagnostic indicates significant spatial dependence, I switch from conventional GEOID-clustered standard errors to the spatially robust procedures described in the following subsections.
When the Moran's-$I$ trigger indicates spatial autocorrelation, I conduct inference using Conditional Spatial Correlation Projection Control (C--SCPC), following Mueller_Watson_JBES_2023. C--SCPC controls size conditional on regressors and spatial locations and is valid in panel or difference-in-differences designs with fixed effects. It extends the spatial-correlation-robust method developed by MuellerWatson2022_Ecta, who establish finite-sample size control and asymptotic validity under weak spatial dependence.
I follow their guidance in selecting the number of principal components $q$ according to the expected-length rule and calibrate the worst-case correlation bound using plausible average spatial correlation. In practice, all specifications for which the Moran's-$I$ trigger is activated are reported in the main tables with C--SCPC standard errors and associated confidence intervals; conventional cluster-robust and SHAC standard errors are relegated to robustness tables and the appendix.
As an appendix robustness check, I compute Conley/Spatial HAC (SHAC) standard errors using distances between municipality centroids measured in the Puerto Rico StatePlane coordinate system (EPSG:32161) and converted to kilometers. The primary specification uses a 75 km cutoff radius for the spatial Bartlett kernel, calibrated to the geographic scale of Puerto Rico. This radius is large enough to capture economically relevant cross-municipality dependence on an island of roughly 170 km by 60 km while remaining small enough to avoid treating all locations as effectively equidistant. I implement this estimator following the methodology outlined in the conleyreg package Duben2025, specifying a Bartlett kernel and lag cutoffs appropriate for the panel structure.
While asymptotic SHAC estimators are standard in applied work, recent contributions highlight that they can suffer from size distortions in finite samples, motivating the use of procedures like C--SCPC as the primary spatial inference tool. To evaluate the sensitivity of spatial SE adjustments, I also include results adjusted by the Thresholding Multiple Outcomes (TMO) method of DellaVignaImbensKimRitzwoller2025, which uses auxiliary outcomes to identify cross-location dependence and typically increases standard errors. In the robustness tables, I report SHAC, SHAC+$\pm$TMO, and C--SCPC standard errors side by side.
As a further robustness check, I estimate two-way cluster-robust standard errors. The canonical approach of CameronGelbachMiller2009_MultiWayCluster constructs two-way clustered variances as the sum of one-way cluster variances minus their intersection term, and it applies readily to non-nested structures such as municipality $\times$ time. To address potential serial correlation in time effects, I also compute the variant proposed by ChiangHansenSasaki2023, which augments the two-way covariance estimator with kernel-weighted time autocovariances and selects the lag truncation using the Andrews rule. I treat these two-way clustered estimates as a sensitivity check relative to the baseline GEOID-clustered and spatially robust procedures described above.
I implement a series of diagnostic and falsification checks designed to assess identification credibility and robustness along several dimensions: pre-trends, sensitivity to violations of parallel trends, overlap and exposure positivity, spatial dependence, and multiple-testing adjustments.
Following SunAbraham2021, I estimate the interaction-weighted (IW) event-study model to diagnose potential violations of parallel trends. This approach corrects the contamination of naive two-way fixed effects (TWFE) estimates when treatment effects are heterogeneous across cohorts. The IW estimator isolates cohort-specific dynamics and recovers unbiased lead coefficients that serve as valid pre-trend diagnostics, avoiding the negative-weight and “forbidden comparison” problems emphasized by RothSantAnnaBilinskiPoe2023. These IW lead coefficients are reported alongside the CS and BJS paths and used purely for diagnostic purposes.
To evaluate robustness to deviations from parallel trends, I apply the HonestDiD sensitivity approach of RambachanRoth2023, focusing on the smoothness (SD) restriction that constrains second differences of the event-study path, $|\Delta^2 \theta_\ell| \le M$, over a grid of curvature budgets $M$. This framework constructs bounds on the ATT that remain valid under bounded violations of parallel trends, calibrated to observed pre-trends.
In the presence of interference, I follow Xu2023 and interpret these bounds as sensitivity checks for both the parallel-trends and exposure-mapping assumptions discussed in Section (ref). The choice of the $M$ grid and the visualization of robust bounds are reported with the main event-study figures and robustness tables.
For settings with interference, the validity of causal identification relies on the overlap, or positivity, of exposure histories. Xu2023 formalize exposure-positivity for DiD with interference, paralleling the support conditions in AronowEcklesSamiiZonszein2020, ShahnZivichRenson2024, and LiuHudgensSaulClemensAliEmch2019. This requirement underlies Assumption A5 in Section (ref).
Building on Savje2023_MisspecExposure, I diagnose overlap using the distributions of the slice-average exposure shares $\bar S_{\cdot,\mathrm{slice}}$ defined in Section (ref) and apply a two-sided $\epsilon$-trim that retains observations only if their slice-average exposure lies in $[\epsilon, 1-\epsilon]$ (default $\epsilon=0.02$), following the trimming logic of crump2009. In the DR--DiD specification in Section (ref), the corresponding slice-sum exposures $S_{\cdot,\mathrm{slice}}$ serve as the main spillover regressors, while trimming is implemented on the slice-average shares $\bar S_{\cdot,\mathrm{slice}}$. I report coverage impacts by channel and slice, alongside histogram and ECDF diagnostics comparing exposure distributions before and after trimming and balancing, and I enforce a minimum post-trim sample-size threshold per slice in the main DR--DiD estimations.
Residual spatial dependence is assessed using the Moran's $I$ residual diagnostic described in Section (ref). For each main specification, I compute Moran's $I$ (or its multivariate extension for stacked residuals) on the model residuals using the spatial-weights matrix $W$ and a permutation reference distribution. Whenever this diagnostic rejects the null of no spatial autocorrelation at the 5% level, I report C--SCPC and SHAC spatially robust standard errors in the main and robustness tables, respectively, instead of relying solely on conventional GEOID-clustered inference.
Finally, I correct for multiple inference separately for the primary dynamic event-study paths and for the heterogeneous-effects summaries.
For the main dynamic event-study paths of the CS and BJS estimators (in particular, the full set of post-treatment coefficients for the primary outcome), I treat the post-treatment horizons as a single family and use the uniform (simultaneous) confidence bands constructed in Section (ref) based on studentized max-$|t|$ (“sup-$t$”) statistics. These bands provide joint coverage over all post-treatment horizons and act as a family-wise multiple-testing adjustment for the event-study path. For cumulative effects such as the 0--16 quarter ATT, I treat the cumulative estimand as a linear combination of the event-time coefficients and construct bootstrap confidence intervals from the same resampling scheme applied to this linear combination.
For the heterogeneous-effects analysis, which uses the auxiliary TWFE specification described in Section (ref), I instead control the false discovery rate (FDR). I treat all heterogeneous DATT, SATT, and TATT coefficients within a given combination of outcome and event-time slice as a single testing family and compute Benjamini--Hochberg (BH) and Benjamini--Yekutieli (BY) $q$-values BenjaminiHochberg1995,BenjaminiYekutieli2001. The BH procedure is reported as the default FDR adjustment, while BY provides a robustness check under arbitrary dependence across coefficients. These adjustments are implemented on the heterogeneity results table and reported alongside the corresponding point estimates and standard errors in the main and appendix tables.
Figure (ref) reports the Sun--Abraham interaction-weighted event-study coefficients for log covered employment around the first sizable entry. Pre-treatment coefficients for $\ell \leq -2$ are close to zero and display wide confidence intervals, with no clear systematic trend before treatment. This pattern is consistent with the conditional parallel trends assumption discussed in Section (ref) and provides the baseline calibration for the HonestDiD sensitivity analysis.
Figure (ref) examines overlap for the main spillover channel. The left panel shows the distribution of slice-average same-industry neighbor exposure shares $\bar S_{\mathrm{same},0\text{--}16}$ before trimming, while the right panel reports the corresponding empirical CDF. Most observations exhibit very low exposure shares, and the two-sided $\epsilon$-rule at $0.02$ and $0.98$ trims a large mass of near-zero exposure cells. Table (ref) quantifies this trimming by event-time slice: between 60% and 82% of raw observations are removed, but all $78$ municipios remain in the trimmed samples, ensuring geographic coverage while enforcing exposure positivity.
Table (ref) summarizes the Moran's $I$ residual diagnostics. For the direct exposure-augmented event-study and for all DR--DiD spillover slices, Moran's $I$ is positive and statistically significant at conventional levels using the Queen+KNN spatial-weights matrix. Examining the time-series of quarter-specific local Moran statistics, between 4% and 18% of quarters show statistically significant spatial clustering across the five specifications reported here (detailed quarter-by-quarter diagnostics are in Appendix Table D1). These findings activate the spatial-dependence “trigger” described in Section (ref) and imply that, for both the direct exposure-augmented event study and the DR--DiD spillover slices, the main reported inference is based on C--SCPC spatially robust standard errors, with cluster-robust and SHAC standard errors reported only as complementary robustness checks.
This section presents the main estimates for the direct and spillover effects of first sizable entries. I begin with dynamic Callaway--Sant'Anna (CS) event-study estimates for the three primary outcomes, then summarize cumulative 0--16 quarter effects. I then turn to the doubly robust DiD (DR--DiD) slice estimands that decompose the total impact on employment into direct and spillover components across the three exposure channels.
Figure (ref) reports the CS event-study for $\log\_\mathrm{covered\_emp}$. All coefficients are normalized to the cohort-specific baseline quarter $g-(\delta+1)$; in the baseline specification with $\delta=2$ this corresponds to event time $\ell=-3$, three quarters before entry, while the last two pre-event quarters $\ell\in\{-2,-1\}$ are treated as a limited-anticipation window used for pre-trend diagnostics. Pre-treatment leads are near zero and well inside the simultaneous confidence band, consistent with the pre-trend diagnostics in Section (ref). At the time of entry, covered employment jumps by roughly 13 log points and continues to rise over the first year, stabilizing between 14 and 17 log points above the pre-entry level during the subsequent three to four years. The 95% uniform band excludes zero for all post-entry horizons in the balanced window.
Figure (ref) shows analogous event studies for log establishments and the log real wage bill. The response of establishments is large and immediate: the number of establishments in the treated municipality--industry cell increases by about 25--30 log points on impact and remains at least 25 log points above the pre-entry level through 16 quarters. The real wage bill also rises quickly, with a pattern similar to covered employment but somewhat more volatile, reflecting both changes in headcounts and in average wages. Across all three outcomes, the dynamic paths display a sharp break at entry, persistent positive effects over four years, and no evidence of systematic pre-trends.
Taken together, these dynamic estimates answer the first part of the research question on direct effects: first sizable entries are followed by large, abrupt, and persistent increases in covered employment, establishments, and wage bills in the treated municipality--industry cell, with no evidence of pre-existing trends that could account for the post-entry jump.
Table (ref) summarizes the CS estimates as cumulative 0--4, 5--8, 9--16, and 0--16 quarter effects, obtained by summing the event-time coefficients $\mathrm{ATT}^{\mathrm{own}}(\ell)$ over all quarters in each slice. For each outcome and slice, I report the sum of the dynamic ATT coefficients, together with bootstrap standard errors and percentile-based 95% confidence intervals.
For covered employment, the cumulative effect is 0.74 log points over quarters 0--4 and 0.55 log points over quarters 5--8, with both slices precisely estimated. The 9--16 quarter slice adds another 1.03 log points, yielding a total 0--16 cumulative effect of 2.32 log points. I return to the economic interpretation of this cumulative effect in Section (ref), where I translate it into approximate job counts and wage-bill changes.
The establishment counts show even larger cumulative responses, with a 0--16 quarter effect of 4.67 log points, while the real wage bill exhibits a 0--16 cumulative effect of 2.33 log points that closely mirrors the employment path. The close alignment between covered employment and the real wage bill suggests that the main margin of adjustment is quantities rather than wages.
Taken together, the cumulative CS estimates indicate that a first sizable entry is associated with a multi-year expansion in activity in the treated municipality--industry cell: total employment-quarters and wage bills increase by roughly an order of magnitude relative to the counterfactual, and the number of establishments rises even more sharply, providing a clear quantitative answer on the magnitude and persistence of the direct effects.
Table (ref) reports the DR--DiD slice-average estimands for $\log\_\mathrm{covered\_emp}$ along the four slices: each entry is a per-period average effect over the corresponding event-time window rather than a cumulative sum. The first column shows the direct effect $\mathrm{DATT}[a,b]$, while the next three columns report the same-industry neighbor, within-municipality cross-industry, and neighbor all-industries spillover components. The final two columns summarize total spillovers $\mathrm{SATT}[a,b]$ and the total effect $\mathrm{TATT}[a,b]=\mathrm{DATT}[a,b]+\mathrm{SATT}[a,b]$. Each estimate is shown with its C--SCPC spatially robust standard error; conventional municipality-clustered standard errors are reported in Appendix (ref) as a comparison, and 95% confidence intervals (using a normal approximation) are based on the C--SCPC standard errors.
Over the first four post-entry quarters, the point estimate of the direct effect is positive, $\widehat{\mathrm{DATT}}[0,4]=0.156$, but under C--SCPC the associated standard error is sizable (SE $\approx 0.13$), so the 95% confidence interval includes zero. The same-industry neighbor spillover remains large in magnitude, $\widehat{\mathrm{SATT}}\_{\mathrm{same}}[0,4]\approx 0.84$, and the implied total 0--4 quarter effect $\widehat{\mathrm{TATT}}[0,4]\approx 0.99$ is similarly positive, but C--SCPC standard errors are much larger than their clustered counterparts and the corresponding confidence intervals no longer exclude zero. Thus, early spillover gains are economically meaningful but not statistically precise once spatial dependence is accounted for.
In the second year after entry (slice 5--8), both the direct and spillover components are small and imprecisely estimated, and the total effect is close to zero. In contrast, the 9--16 quarter slice continues to show a negative same-industry neighbor spillover, $\widehat{\mathrm{SATT}}\_{\mathrm{same}}[9,16]\approx -0.33$. However, the C--SCPC standard error for this effect is large (SE $\approx 0.43$; see Appendix (ref)), so the resulting 95% confidence interval, roughly $[-1.18, 0.51]$, includes zero. The implied total effect $\widehat{\mathrm{TATT}}[9,16]\approx -0.33$ is therefore not statistically distinguishable from zero under our primary spatial inference procedure, even though the point estimate remains negative. These patterns suggest that early gains in neighboring same-industry employment may be partially reversed over longer horizons, consistent with competitive or displacement effects, but the medium-run spillover evidence is statistically weak once spatial dependence is fully accounted for.
Aggregating over the full 0--16 quarter window, the direct effect remains positive in magnitude, $\widehat{\mathrm{DATT}}[0,16]=0.214$, but the C--SCPC standard error (SE $\approx 0.16$) implies a wide confidence interval that includes zero. Net spillovers over the same horizon are negative but very imprecisely estimated, and the total effect $\widehat{\mathrm{TATT}}[0,16]\approx 0.10$ (about 11%) is also far from statistically significant once spatial dependence is taken into account. Appendix (ref) compares these C--SCPC standard errors to conventional cluster-robust and SHAC estimators and shows that clustered inference tends to understate uncertainty in the presence of spatial correlation. Figure (ref) visualizes this decomposition, showing that most of the 0--16 quarter total effect is accounted for by the direct component, with sizable positive short-run spillovers that are offset by longer-run negative same-industry effects.
Finally, it is worth noting that the statistical insignificance of the aggregate cross-industry spillover ($\mathrm{SATT}_{\mathrm{cross}}$) masks heterogeneous linkages between specific sectors. An exploratory pairwise analysis of the underlying source--target variation reveals that entries in manufacturing industries---specifically Apparel (NAICS 315), Printing (323), and Miscellaneous Manufacturing (339)---are associated with statistically significant increases in Educational Services (NAICS 61) establishments in the same municipality (mean post-entry effects ranging from 0.12 to 0.18 log points). This pattern is consistent with induced demand for vocational training or workforce development services following the arrival of new manufacturing employers, suggesting that while broad cross-sector spillovers are weak, specific labor-market complementarities do operate in the local economy.
In terms of the research questions on regional and spillover effects, these DR--DiD estimates indicate that first sizable entries initially generate large positive spillovers to same-industry neighbors, but that over a four-year horizon most of the net change in covered employment comes from the direct effect in the treated municipality--industry cell, with net regional gains that are modest in size and estimated imprecisely once spillovers are taken into account.
This section summarizes how the 0--16 quarter effects on covered employment vary across the pre-defined strata described in Section (ref). I focus on the auxiliary TWFE specification in Section (ref), which interacts the direct-treatment and same-industry neighbor exposure slices with time-invariant indicators for metro status, tradability, and pre-period wage stratum. Figures (ref)--(ref) plot the resulting stratum-specific direct and spillover effects for $\log\_\mathrm{covered\_emp}$ by tradable status, metro status, and wage stratum, respectively, and Table (ref) reports the corresponding estimates, standard errors, confidence intervals, and FDR-adjusted $q$-values.
Across all strata, the 0--16 quarter direct effect of a first sizable entry on covered employment is large and precisely estimated. Point estimates range from about 1.15 to 2.08 log points, implying roughly 3.2--8.0-fold increases in cumulated employment-quarters relative to the counterfactual.\footnote{For example, $\exp(1.146)\approx 3.1$ and $\exp(2.085)\approx 8.0$.} Direct effects are notably larger in tradable industries (2.08 vs 1.15 log points for non-tradable, with barely overlapping confidence intervals) and somewhat larger in high-wage industries (1.95 vs 1.30 log points). Differences between metro and non-metro municipios are more modest relative to their uncertainty bands. In contrast, same-industry neighbor spillovers are uniformly negative: in every stratum, the 0--16 quarter same-industry spillover effect lies between $-1.10$ and $-1.35$ log points, with Benjamini--Hochberg $q$-values below 0.10. These estimates suggest that while first sizable entries generate substantial own-industry gains in both metro and non-metro municipios and in both tradable and non-tradable sectors, they are accompanied by sizable reductions in neighboring same-industry employment, consistent with competitive or business-stealing forces. Population-weighted decompositions in the appendix confirm that these qualitative patterns are robust when weighting municipality--industry cells by their pre-period size.
\footnotesizeNotes: Figures (ref)--(ref) report stratum-specific 0--16 quarter effects on log covered employment from the auxiliary TWFE specification in Section (ref). Blue markers show direct effects (DATT$[0,16]$), and orange markers show same-industry neighbor spillovers (SATT\textsubscript{same}$[0,16]$); vertical lines are 95% confidence intervals based on municipality-clustered standard errors. Filled black stars indicate effects with Benjamini--Hochberg $q$-values below 0.10 within the outcome--slice family. Numerical values are reported in Table (ref).
This section interprets the magnitudes of the estimated direct and spillover effects and discusses their implications for local development policy in Puerto Rico. I focus on covered employment as the primary outcome, using the CS dynamic and cumulative estimates (Section (ref)) together with the DR--DiD decomposition and the heterogeneity patterns in Section (ref).
The CS event-study shows that a first sizable entry has a large and persistent impact on covered employment in the treated municipality--industry cell. Summing the dynamic effects over the 0--16 quarter window (Table (ref)), the cumulative CS estimate for covered employment is \[ \widehat{\mathrm{ATT}}^{\mathrm{own}}_{\mathrm{cumu}}(0,16) ~\approx~ 2.32 \text{ log points}, \] with a narrow confidence interval that lies well above zero. Interpreting this as the cumulative log change in the sum of covered employment over the first four post-entry years, the point estimate corresponds to roughly a $914\%$ increase in employment-quarters: \[ \exp\!\bigl(2.32\bigr) - 1 \;\approx\; 9.1, \] so that the cumulative number of worker-quarters in the treated municipality--industry cell is about ten times larger than in the counterfactual without a sizable entry.
To give one concrete benchmark, suppose a municipality--industry cell would otherwise generate 400 employment-quarters over four years (for example, 100 workers per quarter). The CS estimates imply that, after a first sizable entry, this cumulative total would instead be on the order of \[ (1 + 9.1)\times 400 \;\approx\; 4{,}000 \] employment-quarters over the same horizon. That corresponds to roughly 3{,}600 additional worker-quarters, or about 900 additional job-years.\footnote{These back-of-the-envelope calculations are illustrative and use round numbers for clarity; all inference is based on the log-point estimates and confidence intervals in Table (ref).} The cumulative effect on the real wage bill is almost identical in magnitude (2.33 log points), suggesting that the dominant margin of adjustment is the number of jobs rather than average wages.
The DR--DiD slice estimates provide a complementary perspective. Aggregating over the 0--16 quarter window, the direct slice estimand is \[ \widehat{\mathrm{DATT}}[0,16] \approx 0.214 \quad (\text{SE } 0.161), \] which corresponds to roughly a $24\%$ increase in covered employment when moving from zero to full own-treatment intensity, $\exp(0.214)-1\approx 24\%$. While the point estimate is similar across inference methods, the spatially robust standard error (0.161) is substantially larger than the conventional cluster-robust SE (0.036), reflecting considerably greater uncertainty once spatial correlation is accounted for. This object is not directly comparable to the CS cumulative ATT, because it is defined as a slice-level contrast in the DR--DiD framework rather than as the sum of horizon-specific ATTs. Nonetheless, both sets of estimates indicate that first sizable entries generate large own-industry gains within the treated municipality--industry cell.
The DR--DiD decomposition helps assess whether these direct gains are amplified or offset by spillovers to neighboring and related industries. Over the first four post-entry quarters, the same-industry neighbor channel exhibits marked positive spillovers: \[ \widehat{\mathrm{SATT}}_{\mathrm{same}}[0,4] \approx 0.84, \] while the within-municipality cross-industry and neighbor all-industries channels are close to zero. In this short-run window, the total spillover effect is therefore strongly positive and comparable in magnitude to the direct effect, so that the combined 0--4 quarter total effect, \[ \widehat{\mathrm{TATT}}[0,4] \approx 0.99, \] reflects both own- and neighbor-industry gains (Table (ref)).
In contrast, the 9--16 quarter slice yields negative total spillover effects, \[ \widehat{\mathrm{SATT}}[9,16] \approx -0.34, \] driven almost entirely by the same-industry neighbor channel ($\widehat{\mathrm{SATT}}_{\mathrm{same}}[9,16] \approx -0.33$). Under GEOID-clustered standard errors, this effect is tightly estimated, but it becomes much less precise once I apply the spatially robust C--SCPC adjustments, with a 95% confidence interval that widens to approximately $[-1.18,0.51]$ and includes zero. I therefore interpret these medium-run spillovers as suggestive of crowding-out or business-stealing dynamics---where early neighbor gains are partially reversed as new entrants expand and draw demand or workers away from nearby incumbents---rather than as a definitive, statistically sharp displacement effect.
Aggregating over the full 0--16 quarter window, the net spillover effect is negative but imprecise, \[ \widehat{\mathrm{SATT}}[0,16] \approx -0.11 \quad (\text{C--SCPC SE } 0.908), \] so that the total effect \[ \widehat{\mathrm{TATT}}[0,16] \approx 0.10 \] corresponds to about an $11\%$ increase in covered employment, $\exp(0.10)-1\approx 11\%$, with a wide confidence interval that includes zero. Figure (ref) illustrates that, in this decomposition, most of the 0--16 quarter total effect is accounted for by the direct component, while short-run positive spillovers are offset by longer-run negative same-industry effects.
Taken together, these results imply that first sizable entries substantially increase employment and wage bills in the treated municipality--industry cell, generate meaningful short-run gains for same-industry neighbors, and are consistent with longer-run crowding out within the same industry across space, although the medium-run same-industry spillovers are imprecisely estimated once spatial correlation is accounted for. From a regional perspective, the net 0--16 quarter effect on covered employment is positive but modest, and not sharply estimated, once these offsetting spillover forces are taken into account.
These interpretations should be read in light of the diagnostic and robustness exercises reported earlier in the paper. Pre-trend and overlap checks are broadly consistent with the identifying assumptions of the CS and DR--DiD estimators, and spatial diagnostics motivate the use of spatially robust inference. At the same time, the HonestDiD and alternative-variance exercises underscore that moderate violations of parallel trends or additional spatial dependence could attenuate the estimated effects. The discussion here therefore emphasizes qualitative patterns and the sign and relative magnitude of components, rather than precise point estimates.
The auxiliary heterogeneity analysis in Section (ref) shows that these patterns are broadly similar across metro and non-metro municipios, tradable and non-tradable sectors, and high- and low-wage industries. Stratum-specific direct effects from the TWFE specification are uniformly large and precisely estimated, with 0--16 quarter DATT estimates ranging from about 1.15 to 2.09 log points, implying that cumulative employment-quarters are roughly 3.2--8.0 times higher than the counterfactual. Direct gains tend to be somewhat larger in tradable and high-wage industries, but the differences across strata are small relative to the uncertainty bands.
Spillover effects are negative in every stratum, with 0--16 quarter SATT estimates between about $-1.10$ and $-1.35$ log points and Benjamini--Hochberg $q$-values below 0.10. These total spillover effects are dominated by the same-industry neighbor channel, as shown in the main decomposition. This pattern is consistent with the DR--DiD results: large own-industry gains in treated municipio--industry cells coexist with substantial contractions in neighboring same-industry employment, regardless of whether the sector is tradable or non-tradable, metro or non-metro, high- or low-wage. Population-weighted decompositions (reported in the appendix) confirm that these qualitative patterns are robust when weighting municipality--industry cells by their pre-period size rather than counting each cell equally.
The findings have three main implications for the design and evaluation of local industrial policy in Puerto Rico.
First, at the level of the treated municipality--industry cell, first sizable entries are associated with very large and persistent employment and wage-bill gains. From a narrow cost--benefit perspective that focuses on jobs and payroll in the entering industry and location, the returns to attracting such entries can be substantial. This is the perspective most often emphasized in program evaluations that track jobs created in subsidized plants or zones.
Second, once spillovers across municipios and industries are incorporated, the net gains in covered employment are much smaller and less precisely estimated over the 0--16 quarter horizon. Short-run positive spillovers to same-industry neighbors are offset by longer-run negative effects, and cross-industry and neighbor-all channels are quantitatively modest. For budget and return-on-investment discussions (Q5), this implies that counting only the direct jobs in the entering municipality--industry cell is likely to overstate the net employment gains for the broader regional economy.
Third, the similarity of patterns across metro and non-metro municipios and across tradable and non-tradable sectors suggests that these dynamics are not confined to a particular part of Puerto Rico's production structure. Rather, they appear to reflect general features of how large entries reshape local competition and labor demand across space. This has implications for how incentives are targeted: for example, coordinated policies that internalize cross-municipio competition within the same industry may yield higher net gains than fragmented, municipality-by-municipality bidding for the same firms.
Overall, the results point to a nuanced interpretation of large, persistent entries as development tools. They create substantial own-industry gains where they land, generate some short-run spillovers, but also reallocate activity away from neighboring same-industry locations over time. Evaluating such policies therefore requires moving beyond plant-level job counts to consider the full spatial and cross-industry incidence of both direct and spillover effects.
This section examines the sensitivity of my main findings to alternative estimators, inference procedures, and specification choices. I focus on three key robustness checks: (1) comparing the Callaway--Sant'Anna (CS) and Borusyak--Jaravel--Spiess (BJS) estimators for direct effects, (2) assessing sensitivity to violations of parallel trends using the HonestDiD framework of RambachanRoth2023, and (3) validating the DR--DiD spillover results across multiple inference methods and exposure-history specifications.
Figure (ref) overlays the balanced-horizon event-study estimates of the direct effect of first sizable entries on log covered employment using the CS and BJS estimators. The CS series (solid line) implements the aggregation-based approach with cohort-size weighting and uniform 95% max-$|t|$ confidence bands as described in Section (ref). The BJS series (dashed line) uses outcome-model imputation with cluster-robust pointwise bands. Both estimators produce qualitatively similar dynamic patterns: negligible pre-treatment effects (consistent with the parallel-trends diagnostics in Section (ref)), a sharp positive response beginning at $\ell=0$, and persistent gains through $\ell=16$. The two series track closely in magnitude, with the CS estimates slightly higher in the early post-treatment quarters and the BJS estimates slightly higher in later quarters. The overlap in confidence bands indicates that the choice of estimator does not materially alter my conclusions about the direct effects of entry. This consistency across aggregation-based and imputation-based approaches strengthens confidence in the main results.
To assess whether my findings are robust to violations of the strict parallel-trends assumption, I apply the HonestDiD framework of RambachanRoth2023, which constructs identification-robust confidence intervals under smoothness restrictions on the counterfactual trend. Figure (ref) reports bounds for the cumulative 0--16 quarter ATT on log covered employment under the SD($M$) restriction, which limits the curvature of the trend to at most $M$ in standard-deviation units. The horizontal axis varies $M$ from 0.00 (strict parallel trends) to 0.10 (modest nonlinearity), while the solid line shows the baseline CS estimate of $\text{ATT}^{\mathrm{own}}_{\mathrm{cumu}}(0,16) = 2.317$ log points (from Table (ref)). The shaded region gives the corresponding 95% confidence band. At $M=0.00$, the lower bound coincides with the baseline estimate, and the upper bound extends to approximately 2.55 log points, reflecting sampling uncertainty alone. As $M$ increases to 0.05, the bounds widen slightly but remain strictly positive, with a lower bound near 1.80 log points. Even at $M=0.10$---which allows for substantial curvature---the lower bound remains above 1.40 log points, and the upper bound reaches approximately 3.20 log points. Importantly, zero is excluded from the identification-robust confidence band at all curvature levels I consider. This result indicates that the positive direct effect of first sizable entries on covered employment is robust to modest violations of parallel trends, lending further credibility to the main findings.
Table (ref) compares 0--16 quarter slice effects on log covered employment from the DR--DiD specification across four inference procedures: municipality clustering, Conditional Spatial Correlation Projection Control (C--SCPC), Conley/SHAC spatial standard errors, and the Thresholding Multiple Outcomes (TMO) adjustment. The direct effect (DATT[0,16] = 0.214 log points) is precisely estimated and statistically significant under standard cluster-robust inference ($p < 0.001$) and the TMO adjustment ($p < 0.001$), with standard errors of 0.036 for both methods. However, when accounting for spatial dependence using C--SCPC (SE = 0.161, $p = 0.183$) or SHAC spatial standard errors (SE = 0.296), the direct effect is no longer statistically significant at conventional levels, though the point estimate remains positive and economically meaningful. The total spillover effect (SATT[0,16] = $-0.110$ log points) and the total effect (TATT[0,16] = 0.104 log points) are not statistically significant under any method, as evidenced by standard errors that exceed the point estimates in magnitude. The C--SCPC standard error for SATT[0,16] is 0.908, reflecting the spatial dependence documented in Section (ref). The TMO-adjusted standard errors (final column) are nearly identical to the cluster-robust standard errors, indicating that the multiple-testing adjustment does not substantively change inference for these key parameters. Overall, while the direct effect remains positive and economically meaningful across all inference methods, its statistical significance depends critically on whether spatial dependence is accounted for. The spillover and total effects remain imprecise under all methods. These findings suggest that the main gains from first sizable entries accrue locally within the treated municipality--industry cells, but the strength of this conclusion is sensitive to assumptions about spatial correlation in the error structure.
Table (ref) examines the sensitivity of the 0--16 quarter DR--DiD effects to alternative exposure-history summaries. The baseline specification uses an “any-since-adoption” history, which flags exposure as active if the neighbor (or within-municipality partner industry) has ever experienced a first sizable entry prior to the current quarter. An alternative specification uses a “last four quarters” history, which only counts exposure from entries that occurred in the four most recent quarters. The table reports direct (DATT[0,16]) and total (TATT[0,16]) effects under both definitions. Estimates and standard errors (both cluster-robust and C--SCPC) are virtually identical across the two specifications, with DATT[0,16] = 0.214 log points (cluster SE = 0.036, C--SCPC SE = 0.161) and TATT[0,16] = 0.104 log points (cluster SE = 0.149, C--SCPC SE = 0.965) in both cases. This stability suggests that the spillover dynamics are not driven by the choice of exposure-history window: whether I summarize exposure using the full adoption history or only the recent four quarters, the estimated effects remain unchanged. This result further validates the robustness of the DR--DiD spillover estimates and underscores the limited role of indirect effects in this context.
This paper has examined the employment and wage-bill consequences of first sizable industry entries in Puerto Rico. Rather than focusing on a specific program or subsidy, I use administrative data on establishments, covered employment, and wage bills to identify municipio--industry cells that experience a large and persistent jump in economic activity. I then combine modern staggered-adoption difference-in-differences methods with a doubly robust decomposition that explicitly allows for and separates spatial spillovers along three channels. The analysis addresses five questions about the dynamic direct effects of first entries, their spillovers across space and sectors, heterogeneity across pre-defined strata, and the implications for budgeting and return-on-investment discussions.
The first conclusion is that first sizable entries generate very large and persistent direct gains in the treated municipio--industry cell. Callaway--Sant'Anna event studies show sharp breaks at the time of entry, with covered employment, the number of establishments, and the real wage bill all jumping on impact and remaining well above pre-entry levels for at least four years. Cumulative 0--16 quarter estimates imply roughly a nine-fold increase in employment-quarters and a similarly sized increase in the wage bill, indicating that the dominant margin of adjustment is the number of jobs rather than average wages. From the perspective of the host municipio and industry, these events are therefore associated with substantial local development.
The second conclusion is that these large local gains coexist with non-trivial spillovers to neighboring municipios in the same industry. The DR--DiD decomposition indicates that same-industry neighbors experience sizeable short-run gains in the first post-entry year, but that these effects reverse over the 9--16 quarter window, with point estimates indicating declines in same-industry neighbor employment. Spillovers to other industries in the same municipio and to all industries in neighboring municipios are small and imprecisely estimated. Aggregating over all channels, net spillovers over the 0--16 quarter horizon are negative but noisy, and the implied total regional effect is modest in magnitude and not sharply estimated, especially once spatially robust inference is applied.
Third, auxiliary heterogeneity analyses suggest that these patterns are not confined to a particular part of Puerto Rico's production structure or geography. Stratum-specific estimates indicate large and positive 0--16 quarter direct effects across metro and non-metro municipios, tradable and non-tradable sectors, and high- and low-wage industries, with somewhat larger gains in tradable and high-wage strata. Same-industry neighbor spillovers are negative in every stratum and remain so after controlling the false discovery rate. Population-weighted decompositions confirm that these qualitative findings are robust when municipio--industry cells are weighted by their initial size. Overall, the data point to a general pattern in which first sizable entries strongly expand activity in the treated municipio--industry, while reallocating same-industry employment away from neighboring locations.
These findings have direct implications for the evaluation and design of local industrial policy in Puerto Rico. From a narrow standpoint that counts jobs and payroll in the entering municipio--industry cell, first sizable entries appear highly attractive: they deliver large and persistent increases in covered employment and the wage bill. From a broader regional perspective, however, the picture is more muted. Once losses in neighboring same-industry municipios and the limited cross-industry and neighbor-all spillovers are taken into account, the net 0--16 quarter gain in covered employment across affected areas is much smaller and estimated with substantial uncertainty. This suggests that plant-level job counts, as commonly reported in program evaluations, are likely to overstate the net employment gains relevant for budgeting and return-on-investment discussions at the island level. It also underscores the potential value of more coordinated incentive policies that internalize cross-municipio competition within the same industry, rather than municipality-by-municipality bidding for the same firms.
The analysis has several limitations. The outcomes are limited to covered private-sector employment and wage bills and do not capture informal work, self-employment, or broader measures of welfare such as household income or productivity. The identification strategy relies on conditional parallel trends and exposure-positivity assumptions that, while supported by pre-trend, overlap, and spatial-dependence diagnostics, cannot be verified directly. HonestDiD sensitivity analyses and spatially robust variance estimators show that moderate deviations from these assumptions could attenuate the estimated effects and widen the range of plausible net regional impacts. In addition, the study focuses on the medium-run 0--16 quarter horizon; it does not speak to longer-run adjustments in firm dynamics, human capital, or migration.
These limitations suggest several directions for future work. Linking the administrative data used here to household surveys or tax records would allow a richer assessment of how first sizable entries affect earnings, poverty, and distributional outcomes. Extending the framework to incorporate public finance outcomes—such as tax revenues and infrastructure costs—would speak more directly to the fiscal returns on industrial incentives. Finally, applying similar methods to other settings with different institutional environments and spatial structures would help assess the external validity of the Puerto Rican experience. Despite these caveats, the evidence presented here points to a nuanced view of large entries as development tools: they create substantial and persistent gains where they land but generate only modest and uncertain net employment gains for the broader regional economy, highlighting the importance of accounting for both direct and spillover effects when designing and evaluating place-based policies.
During the preparation of this work the author used Claude (Anthropic) to generate Python code for data processing, statistical analysis, and visualization. The author reviewed and edited the code and takes full responsibility for the content.