The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
78,515 characters
Endogenous Selection and Spillovers: Bayesian Inference for Policy-Relevant Causal Effects
\maketitle
\thispagestyle{empty}
\begin{abstract}
\begin{spacing}{1.15}
This paper develops a new econometric framework to identify and estimate policy-relevant causal effects in contexts with endogenous selection into treatment and spillovers within single large networks or spatial settings. Conventional causal inference methods relying on either unconfoundedness or no-interference assumptions are generally inadequate in these scenarios. We introduce a Spillover Roy model that jointly models endogenous treatment selection and potential outcomes while allowing spillovers through a low-dimensional exposure mapping of neighbors' treatments. The model captures heterogeneous treatment responses across levels of latent resistance to treatment and neighborhood exposure. Within this framework, we define policy-relevant direct, spillover, and total effects under feasible policy changes and show that the total effect decomposes into a direct component from policy-induced participation and a spillover component from policy-induced changes in neighborhood treatment exposure. For estimation and inference, we develop a Bayesian data-augmentation algorithm with parameter expansion that enables efficient posterior computation and coherent uncertainty quantification for heterogeneous causal effects and policy counterfactuals. An application to the U.S. Opportunity Zones program finds positive direct effects on housing development but limited spillover benefits, while counterfactual policy analysis reveals diminishing returns from program expansion.
\smallskip
\noindent\textit{Keywords:} spillovers, network interference, spatial interference, endogenous selection, policy-relevant treatment effects, Bayesian inference, place-based policy.
\medskip
\noindent\textit{JEL classifications: C11, C31, C35, C36, R58.}
\end{spacing}
\end{abstract}
\newpage
\pagenumbering{arabic}
\section{Introduction}\label{introduction}
Spillovers, often referred to as interference, and endogenous selection into treatment are pervasive features of economic settings that complicate conventional causal inference approaches. Spillovers arise when a unit's outcome depends not only on its own treatment status but also on the treatments received by others within its social network or geographic area \citep[see, e.g.,][]{forastiere2021identification, giffin2022generalized}. For example, participants in after-school programs may influence the behavior or academic performance of non-participating peers, while place-based policies such as regional development incentives can affect neighboring communities through migration or business relocation. Simultaneously, participation in many programs is not randomly assigned \citep[see, e.g.,][]{abbring2007econometric}. Individuals or regions often enter treatment based on unobserved characteristics, such as motivation or growth potential, which are also related to potential outcomes, generating endogenous selection. In such environments, policy changes influence outcomes through two interconnected mechanisms: by altering who participates in the program and by reshaping the treatment exposure faced by others. Ignoring either mechanism can lead to biased estimates of program effectiveness and misleading conclusions regarding policy design and evaluation.
Nonetheless, the issue of policy-relevant causal inference in settings where both endogenous selection and spillover effects occur remains insufficiently explored due to the intertwined challenges. Existing approaches typically focus on identifying causal effects at fixed treatment or exposure states may not correspond directly to feasible policy changes. In practice, policymakers are interested in the consequences of modifying program rules, such as expanding eligibility, changing subsidies, or altering designation criteria, beyond the scope of existing treatment effects. As a result, conventional treatment-effect parameters do not fully capture multiple channels of policy interventions in such settings.
This paper develops a new econometric framework for policy-relevant causal effects that simultaneously handle endogenous selection and spillovers in large non-clustered networks or spatial settings. We extend the Generalised Roy model to accommodate network or spatial interactions by allowing potential outcomes to depend on both individual treatment status and an exposure measure summarizing neighbors' treatment. This structure preserves the economic interpretation of selection into treatment while incorporating spillovers in a tractable way. Our resulting Spillover Roy model captures heterogeneous treatment responses across the latent resistance distribution and across exposure levels. Within this framework, we define policy-relevant causal effects as contrasts between expected outcomes under alternative feasible policy regimes. A policy shift may affect outcomes through two distinct channels: induced changes in own treatment participation and induced changes in exposure to treatments experienced by neighbors. Our framework therefore decomposes the total policy impact into direct and spillover components, providing a transparent characterization of the mechanisms through which policies operate.
Our methodology contributes to three strands of research. \emph{First}, it relates to the literature on causal inference under interference. One approach imposes partial interference, whereby spillovers operate within exogenously defined clusters \citep[e.g.,][]{hudgens2008toward, sobel2006randomized, manski2013identification}, with recent work allowing for noncompliance or endogenous treatment take-up \citep{ditraglia2023identifying, vazquez2023causal}. Our setting is closer to work on general interference in a single network or spatial environment, where a low-dimensional exposure mapping summarizes the relevant treatment configuration. Existing methods typically rely on randomized treatment or unconfoundedness conditional on observed covariates \citep[e.g.,][]{aronow2017estimating, leung2020treatment, forastiere2021identification, forastiere2022estimating}. Recent studies relax treatment exogeneity: \citet{hoshino2024causal} use instrumental exposure mappings to identify local direct and indirect effects under noncompliance, while \citet{chen2025heterogeneous} model treatment choices as a network equilibrium and identify heterogeneous marginal exposure effects. Related work studies policy effects when interference is mediated through market-equilibrium variables under randomized treatment \citep{munro2025treatment}. Our setting instead combines endogenous treatment selection with network or spatial exposure and focuses on counterfactual changes in the rule governing treatment participation.
\emph{Second}, our policy estimands build on the literature concerning marginal and policy-relevant treatment effects under endogenous selection \citep[e.g.,][]{heckman2005structural, heckman2007econometric, carneiro2011estimating, mogstad2018using, sasaki2023estimation, opper2024late}. In the Generalized Roy framework, treatment is governed by a latent-index selection rule, and the treatment effects may vary with the unobserved determinants of treatment choice. The marginal treatment effect (MTE) characterizes treatment gains along this latent resistance margin and provides a building block for policy evaluation. A central insight of this literature is that policy evaluation requires specifying the policy-relevant target population: the individuals whose treatment choices would change under the counterfactual policy may not coincide with those whose choices are shifted by the available instrument. Among various measures or treatment effects, the PRTE has the advantage of directly evaluating alternative policy scenarios under consideration and can be represented as policy-specific weighted averages of the underlying MTE. PRTE evaluates a change from a baseline to a counterfactual treatment-selection regime and, when normalized by the change in participation, measures the average outcome gain per net participant induced by the policy. We extend this policy-evaluation approach to settings with interference.
A policy-induced change in treatment selection now affects outcomes through both own participation and the resulting change in neighbors' treatment exposure. We therefore define policy-relevant direct and spillover effects and show how the total policy impact decomposes into these two components, linking policy evaluation to heterogeneity along both the latent resistance and neighborhood-exposure margins.
\emph{Third}, our estimation strategy relates to Bayesian methods for endogenous selection and latent-index models, including inference on potential-outcome distributions in the presence of unidentified dependence parameters \citep{poirier2003predictive} and parameter-expanded data-augmentation techniques for handling covariance restrictions and identifying normalizations \citep[e.g.,][]{ding2014bayesian, dougan2018bayesian, zhang2026parameter}. Building on this literature, we develop a parameter-expanded Gibbs sampler tailored to the richer latent structure of the Spillover Roy model, which jointly estimates endogenous treatment selection and regime-specific potential outcomes and delivers posterior inference for heterogeneous treatment, spillover, and policy-relevant effects. Monte Carlo experiments show that our proposed proposed procedure yeilds valid inference, while naive approaches that ignore either endogenous selection or spillovers can exhibit substantial bias.
We apply the proposed framework to evaluate the causal effects of the Opportunity Zones (OZ) program on housing growth in U.S. census tracts. We model the designation of treated areas as an endogenous selection process influenced by local economic characteristics and political decisions. Our empirical results reveal substantial heterogeneity in treatment gains consistent with selection-on-gains behavior: areas with higher expected returns are more likely to receive the program. We find positive direct effects of designation on housing development but limited evidence of beneficial spillovers for neighboring non-designated areas. Furthermore, policy counterfactual analysis shows that expanding the program induces diminishing returns as marginal entrants generate smaller --- and eventually negative --- direct gains, while spillover benefits increase but remain insufficient to sustain positive net effects under large expansions.
The remainder of this paper is structured as follows. In Section \ref{BCIES-Section2}, we present Spillover Roy Model and define causal estimands with key identification assumptions. In Section \ref{BCIES-Section3}, we propose Bayesian data-augmentation approach to estimate the model and conduct inference. We then evaluate our method using simulations in Section \ref{BCIES-Section4} and investigate the causal impact of the U.S. Opportunity Zones (OZ) program on economic outcomes in Section \ref{BCIES-Section5}. Finally, Section \ref{BCIES-Section6} concludes with brief remarks and policy recommendations.
\section{The Spillover Roy Model}\label{BCIES-Section2}
\subsection{General Model Setup}\label{BCIES-Section2_1}
We consider a general setting for \(n\) agents (\(i = 1,\ldots,n\)) which involves treatment selection and outcome determination with spillovers.
\textbf{Treatment selection}
Let \(D_i\) be the observed binary treatment decision, which takes the value of \(1\) if the unit receives the treatment and \(0\) otherwise. This can be regarded as individual treatment and determined by a latent-index representation as follows
\begin{equation}
\begin{split}
D_i^* &= \mu(Z_i,X_i) - U_i,\\
D_i &= 1 \text{ if } D_i^* \geq 0, \quad D_i = 0 \text{ otherwise},
\end{split}
\end{equation}
where \(D_i^*\) denotes the net benefit, or latent utility, from receiving the treatment and \(U_i\) captures an unobserved component of the treatment choice. Vector \(X_i\) contains observed characteristics that may jointly influence treatment participation and potential outcomes. To identify causal effects under endogenous selection, we additionally observe a vector of excluded variables \(Z_i\), which shifts treatment participation without directly affecting potential outcomes. Throughout the paper, \(Z_i\) plays the role of an instrumental variable in the structural selection equation. Variation in \(Z_i\) provides exogenous changes in treatment propensity while satisfying the exclusion restriction imposed later in the identification analysis.
Assume that \(U_i\) is continuously distributed with a strictly increasing cumulative distribution function \(F_U\). Define \(V_i \coloneqq F_U(U_i)\), then it has uniformly distribution and indicates different quantile level of \(U_i\). Also define \(\nu(Z_i,X_i) \coloneqq F_U(\mu(Z_i,X_i))\), which is the mean scale utility function in discrete choice theory, we can thereby rewriting the treatment rule as:
\begin{equation}
D_i = \mathbbm{1} \{\mu(Z_i,X_i) \geq U_i\} = \mathbbm{1} \{\nu(Z_i,X_i) \geq V_i\}.
\end{equation}
The variable \(V_i\) indexes the unit's latent resistance to treatment: units with higher values of are less likely to participate for a given value of \(\nu(Z_i,X_i)\). This representation is standard in the Generalized Roy model \citep{heckman2005structural} and will be central for defining marginal treatment effects and policy-relevant effects. Because the latent resistance \(V_i\) may be statistically dependent on the potential-outcome disturbances, treatment selection mechanism is generally endogenous.
\textbf{Outcome determination with spillovers}
Under interference, potential outcomes of unit \(i\) may depend not only on its individual treatment but also on the treatment assignments of other units connected to \(i\) through a network or spatial interaction structure. Let \(\mathbf{D} = (D_1,\ldots,D_n)\) denote the population treatment vector, and let \(Y_i(D_i,\mathbf{D}_{-i})\) denote the potential outcome of unit \(i\) under the treatment assignment \(\mathbf{D} = (D_1,\ldots,D_n)\), where \(\mathbf{D}_{-i} = (D_1, \ldots, D_{i-1}, D_{i+1}, \ldots, D_n)\) collects the treatment assignments of all units other than \(i\). Without additional restrictions, unit \(i\) may have a distinct potential outcome for each possible realization of \(\mathbf{D}\). Because the number of potential outcomes grows exponentially with the number of units, such a framework is generally infeasible for both identification and estimation in large networks. Following the network causal inference literature \citep[see, e.g.,][]{aronow2017estimating, manski2013identification}, we impose an exposure-mapping restriction that summarizes the aspects of neighbors' treatment assignments relevant for unit \(i\)'s outcome.
\begin{assumption}{[Exposure mapping]}
\label{Assp1}
There exists a known measurable mapping $T:\{1,\ldots,n\}\times\{0,1\}^{n-1}\times\mathcal W\rightarrow\mathcal T$ such that, for each unit $i$,
$$
T_i=T(i,\mathbf{D}_{-i},\mathbf{W}),
$$
where $\mathbf W\in\mathcal W$ denotes the observed network or spatial structure and $\mathcal T$ is the exposure space. Depending on the application, $\mathbf W$ may be either a network adjacency matrix or a spatial weights matrix. For every $d\in\{0,1\}$ and every pair of treatment vectors $\mathbf d_{-i}$ and $\mathbf d'_{-i}$,
$$
T(i,\mathbf d_{-i},\mathbf W) = T(i,\mathbf d'_{-i},\mathbf W) \quad \Longrightarrow \quad Y_i(d,\mathbf d_{-i}) = Y_i(d,\mathbf d'_{-i}).
$$
\end{assumption}
Thus, conditional on unit \(i\)'s own treatment status, the potential outcome depends on the treatment assignments of all other units only through the exposure state \(T_i\), rather than through the full treatment vector \(\mathbf D_{-i}\). Consequently, there exists a function \(Y_i^{(d)}: \mathcal T \rightarrow \mathbbm{R}\) such that
\[
Y_i(d,\mathbf D_{-i}) = Y_i^{(d)}(T_i),\qquad d\in\{0,1\}.
\]
Throughout the paper, the interaction matrix \(\mathbf{W}\) is assumed to be known and predetermined with respect to the latent disturbances in the treatment and outcome equations. This matrix characterizes the pattern of either network interference or spatial interference. For the remainder of the paper, we focus on the weighted neighborhood treatment as the exposure mapping
\[
T_i = \bar D_{\mathcal N i} = \sum_{j\neq i}w_{ij}D_j,
\qquad w_{ii}=0,
\qquad \sum_{j\neq i} w_{ij} = 1,
\]
The scalar \(\bar D_{\mathcal N i}\) summarizes the treatment intensity in unit \(i\)'s neighborhood.\footnote{When \(w_{ij} = 1/N_i\) for all neighbors \(j\in \mathcal N (i)\), where \(N_i\) denotes the number of neighbors of unit \(i\), \(\bar D_{\mathcal N i}\) reduces to the proportion of treated neighbors. More generally, the weights may reflect heterogeneous interaction strengths, geographic proximity, or other measures of network influence, yielding a weighted neighborhood treatment intensity.} Other commonly used exposure mappings include binary exposure indicators (e.g., at least one treated neighbor), unweighted treated-neighbor counts, higher-order neighborhood summaries, and other nonlinear exposure measures.
Under the weighted neighborhood exposure mapping, we specify the treated and untreated potential outcomes as
\begin{equation}
\begin{split}
Y_i^{(1)}(\bar{D}_{\mathcal{N}i}) &= \mu_1\left(\bar{D}_{\mathcal{N}i}, X_i\right) + \varepsilon_i^{(1)}, \quad \text{ and }\\
Y_i^{(0)}(\bar{D}_{\mathcal{N}i}) &= \mu_0\left(\bar{D}_{\mathcal{N}i}, X_i\right) + \varepsilon_i^{(0)},
\end{split}
\end{equation}
where \(X_i\) denotes a vector of observed individual characteristics, \(\mu_d\left(\bar{D}_{\mathcal{N}i},X_i\right) = \mathbbm{E} \left[ Y_i^{(d)} \mid \bar{D}_{\mathcal{N}i}, X_i\right]\) is the regime-specific mean response function, and \(\varepsilon_i^{(d)}\) is the corresponding idiosyncratic disturbance for each \(d\in\{0,1\}\).
Let \(Y_i\) be the revealed outcome, which equals the treated potential outcome when \(i\) is treated (\(D_i = 1\)) and equals untreated potential outcome when \(i\) is untreated (\(D_i = 0\))
\begin{equation}
Y_i = D_i Y_i^{(1)} + (1-D_i) Y_i^{(0)}.
\end{equation}
Combining endogenous treatment selection with spillovers operating through the specified exposure mapping, the general framework is
\begin{equation}
\label{KeyModel1}
\begin{split}
D_i &= \mathbbm{1} \{\nu(Z_i,X_i) \geq V_i\},
\\
\bar{D}_{\mathcal{N}i} &= \sum_{j=1, j\neq i}^{n} w_{ij} D_j, \quad \sum_{j=1, j\neq i}^{n} w_{ij} = 1,
\\
Y_i^{(1)} &= \mu_1\left(\bar{D}_{\mathcal{N}i}, X_i\right) + \varepsilon_i^{(1)},
\\
Y_i^{(0)} &= \mu_0\left(\bar{D}_{\mathcal{N}i}, X_i\right) + \varepsilon_i^{(0)},
\\
Y_i &= D_i Y_i^{(1)} + (1-D_i) Y_i^{(0)}.
\end{split}
\end{equation}
\begin{example}
(Place-based policies with spatial spillovers) Consider the Opportunity Zones (OZ) program, where $D_i$ indicates OZ designation of census tract $i$ and $Y_i$ denotes a local economic outcome, such as housing development. Designation may be endogenous because unobserved local characteristics, such as growth potential or political support, can affect both designation and potential outcomes. Thus, latent resistance $V_i$ may be correlated with the unobserved determinants of potential outcomes, while an excluded variable $Z_i$ provides exogenous variation in designation. OZ designation may also affect nearby tracts through housing-market, migration, or business-location responses, generating spatial spillovers through the neighborhood exposure defined above.
\end{example}
\begin{example}
(Educational programs with network spillovers) Consider an after-school program (ASP), where $D_i$ indicates participation of student $i$ and $Y_i$ denotes an outcome such as academic performance or social-emotional development. Because participation is self-selected, latent resistance $V_i$ may be correlated with unobserved determinants of potential outcomes, generating endogenous treatment selection. An excluded cost shifter $Z_i$ provides exogenous variation in participation. Program participation may also affect students indirectly through interactions with participating peers, generating network spillovers through the neighborhood exposure defined above.
\end{example}
\subsection{Parametric Spillover Roy Model and Identification}\label{BCIES-Section2_2}
We now introduce the assumptions under which the structural objects in \eqref{KeyModel1} are identified.
\begin{assumption}{[Parametric Spillover Roy model]}
\label{Assp2}
Suppose the data are generated by
\begin{equation}
\label{KeyModel2}
\begin{split}
D_i &= \mathbbm{1}\{Z_i\alpha + X_i\beta^{(D)} + \varepsilon_i^{(D)} > 0\},\\
\bar{D}_{\mathcal{N}i} &= \sum_{j=1, j\neq i}^{n} w_{ij} D_j, \quad \sum_{j=1, j\neq i} w_{ij} = 1,\\
Y_i^{(1)} &= \delta^{(1)}\bar{D}_{\mathcal{N}i} + X_i\beta^{(1)} + \varepsilon_i^{(1)},\\
Y_i^{(0)} &= \delta^{(0)}\bar{D}_{\mathcal{N}i} + X_i\beta^{(0)} + \varepsilon_i^{(0)}, \\
Y_i &= D_i Y_i^{(1)} + (1 - D_i) Y_i^{(0)}.\\
\end{split}
\end{equation}
\end{assumption}
This specification imposes linearity in the mean response functions and the treatment index while accommodating endogenous selection through unrestricted dependence between the treatment-selection disturbance, \(\varepsilon_i^{(D)}\) and the potential-outcome disturbances, \(\bigl(\varepsilon_i^{(1)},\varepsilon_i^{(0)}\bigr)\). Substituting the potential outcomes into the switching equation yields
\[
Y_i = X_i \beta^{(0)} + D_i X_i (\beta^{(1)} - \beta^{(0)}) + \delta^{(0)} \bar{D}_{\mathcal{N}i} + D_i(\delta^{(1)} - \delta^{(0)})\bar{D}_{\mathcal{N}i} + (1-D_i)\varepsilon_i^{(0)} + D_i\varepsilon_i^{(1)},
\]
which makes clear that the model allows the individual outcome to depend on own treatment, neighborhood treatment, and their interaction.
\emph{Remark}. The specification in Assumption \ref{Assp2} extends the Generalized Roy
framework \citep{heckman2005structural} by allowing potential outcomes to depend on neighborhood treatment exposure, with potentially different spillover effects across treatment states, \(\delta^{(1)}\) and \(\delta^{(0)}\). When \(\delta^{(0)}=\delta^{(1)}=0\), the model reduces to the canonical Generalized Roy model.
\begin{assumption}{[Instrument validity]\\}
\label{Assp3}
\textnormal{(i) Instrument exogeneity.}
$$
\left(
\varepsilon_i^{(D)},
\varepsilon_i^{(1)},
\varepsilon_i^{(0)}
\right)
\perp\!\!\!\perp
(X_i,Z_i),
$$
where $\perp\!\!\!\perp$ denotes the statistical independence.
\textnormal{(ii) Instrument relevance.}
The excluded variable $Z_i$ generates nondegenerate variation in the treatment-selection index
$$
\nu_i = Z_i\alpha+X_i\beta^{(D)}.
$$
\end{assumption}
Assumption \ref{Assp3} requires the instrumental variable to satisfy both exogeneity and relevance. Part (i) states that the observed covariates and the instrumental variable are jointly independent of the latent disturbances governing treatment selection and potential outcomes. Together with the structural specification in Assumption \ref{Assp2}, this implies that the instrumental variable affects outcomes only through its effect on treatment selection. Part (ii) requires the instrument to generate sufficient variation in the latent treatment-selection index, thereby ensuring identification of the treatment-selection equation. These conditions are standard in structural latent-index models of endogenous treatment selection and heterogeneous treatment effects \citetext{\citealp[see, e.g.,][]{carneiro2011estimating}; \citealp{brinch2017beyond}; \citealp[and][]{cornelissen2018benefits}}.
\begin{assumption}{[Finite-mixture distribution and cross-unit independence]}
\label{Assp4}
Conditional on $\mathbf{W}$, the disturbance vectors are independently distributed across units and,
for $i = 1,\ldots,n$
\begin{equation}
\label{AS-finitemixture}
\begin{split}
\boldsymbol{\varepsilon}_i = \begin{bmatrix}
\varepsilon_i^{(D)} \\ \varepsilon_i^{(1)} \\ \varepsilon_i^{(0)}
\end{bmatrix} &\overset{ind}{\sim} \sum_{g=1}^G \pi_g N\left(0, \Sigma_g \right) \\
\end{split}
\end{equation}
$$
\begin{aligned}
\text{where} \quad \sum_{g=1}^G \pi_g = 1;
\quad \Sigma_g = \begin{bmatrix}
1 & \sigma_{1D,g} & \sigma_{0D,g} \\
& \sigma_{1,g}^2 & \sigma_{10,g} \\
& & \sigma_{0,g}^2
\end{bmatrix}
\end{aligned}
$$
\end{assumption}
Assumption \ref{Assp4} models the joint distribution of the latent disturbances as a finite mixture of multivariate normal distributions. The first diagonal element of each component covariance matrix is normalized to one, reflecting the standard scale normalization in binary latent-index models \citep[see, e.g.,][]{cameron2005microeconometrics, chan2019bayesian}. Without this normalization, the parameters in the selection equation are identified only up to scale. This finite-mixture specification flexibly approximates the joint distribution of the latent disturbances and allows for non-Gaussian heterogeneity while preserving tractability for estimation and inference. In addition, we note that the latent resistance \(V_i\) in \eqref{KeyModel1} can now be represented as
\[
V_i = \Phi(-\varepsilon_i^{(D)}),
\]
where \(\Phi(\cdot)\) denotes the standard normal cdf. Under Assumption \(\ref{Assp4}\) and the normalization \(\mathop{\mathrm{Var}}(\varepsilon_i^{(D)})=1\), it follows that \(V_i\sim U(0,1)\).
The conditional independence across units implies that, conditional on the predetermined network structure \(\mathbf{W}\), neighborhood treatment carries no additional information about unit \(i\)'s latent disturbances beyond that contained in its own treatment decision and observed covariates. Accordingly,
\[
\mathop{\mathrm{\mathbbm{E}}}[\varepsilon_i^{(d)} \mid D_i, \bar{D}_{\mathcal{N}i}, X_i, Z_i, \mathbf{W}] = \mathop{\mathrm{\mathbbm{E}}}[\varepsilon_i^{(d)} \mid D_i, X_i, Z_i, \mathbf{W}], \quad d\in\{0,1\}.
\]
\begin{theorem}[Identification of the Spillover Roy Model]
\label{Theorem1}
Suppose Assumptions \ref{Assp1}--\ref{Assp4} hold. Then the parameters of the Spillover Roy model in \eqref{KeyModel2},
$$
\alpha,\;\beta^{(D)},\;\beta^{(1)},\;\beta^{(0)},\;\delta^{(1)},\;\delta^{(0)},
$$
and the aggregate covariance parameters
$$
\sigma_{1D}
\coloneqq
\mathop{\mathrm{Cov}}\!\left(\varepsilon_i^{(1)},\varepsilon_i^{(D)}\right),
\qquad
\sigma_{0D}
\coloneqq
\mathop{\mathrm{Cov}}\!\left(\varepsilon_i^{(0)},\varepsilon_i^{(D)}\right),
$$
are identified.
\end{theorem}
\emph{Proof.} See Appendix \ref{BCIES-APDX_ProofProp1}.
\subsection{Causal Estimands}\label{causal-estimands}
Our targeted estimands include a hierarchy of causal parameters: marginal (structural) effects, average effects, and policy-relevant effects.
\textbf{Marginal Structural Objects}
We begin with the primitive objects that characterize heterogeneity in both selection into treatment and neighborhood exposure.
We define the \emph{Marginal Treatment Effect (MTE)} under interference as
\begin{equation}
\label{eq:MTE-bcies}
\mathop{\mathrm{MTE}}(\bar d_{\mathcal N},v,x)
\coloneqq
\mathop{\mathrm{\mathbbm{E}}}\!\left[
Y_i^{(1)}(\bar d_{\mathcal N})-Y_i^{(0)}(\bar d_{\mathcal N})
\mid V_i=v,\; X_i=x
\right].
\end{equation}
This generalizes the classical marginal treatment effect to settings with interference by allowing treatment effects to vary with both latent selection heterogeneity \(v\) and the neighborhood exposure level \(\bar d_{\mathcal N}\).
We define the \emph{Marginal Spillover Effect (MSE)} as
\begin{equation}
\label{eq:MSE-bcies}
\mathop{\mathrm{MSE}}(d,\bar d_{\mathcal N},v,x)
\coloneqq
\frac{\partial}{\partial \bar d_{\mathcal N}}
\mathop{\mathrm{\mathbbm{E}}}\!\left[
Y_i^{(d)}(\bar d_{\mathcal N})
\mid V_i=v,\; X_i=x
\right],
\qquad d\in\{0,1\},
\end{equation}
This measures the local causal response of potential outcomes to a marginal increase in neighborhood treatment exposure for individuals with treatment status \(d\) and latent resistance \(v\).
\begin{theorem}[Identification of Marginal Structural Objects]
\label{Theorem2}
Suppose Assumptions \ref{Assp1}--\ref{Assp4} hold. Then the marginal structural objects
$$
\mathop{\mathrm{MTE}}(\bar d_{\mathcal N},v,x)
\quad\text{and}\quad
\mathop{\mathrm{MSE}}(d,\bar d_{\mathcal N},v,x),\; d\in\{0,1\},
$$
are identified. In particular,
\begin{align}
\label{eq:MTE-identified}
\mathop{\mathrm{MTE}}(\bar d_{\mathcal N},v,x)
&=
(\delta^{(1)}-\delta^{(0)})\bar d_{\mathcal N}
+
x(\beta^{(1)}-\beta^{(0)})
+
\mathop{\mathrm{\mathbbm{E}}}\!\left[
\varepsilon_i^{(1)}-\varepsilon_i^{(0)}
\mid V_i=v
\right],
\\
\label{eq:MSE-identified}
\mathop{\mathrm{MSE}}(d,\bar d_{\mathcal N},v,x)
&=\delta^{(d)}, \qquad d\in\{0,1\}.
\end{align}
Moreover, under Assumption \ref{Assp4},
\begin{equation}
\label{eq:eps-diff-generic}
\mathop{\mathrm{\mathbbm{E}}}\!\left[
\varepsilon_i^{(1)}-\varepsilon_i^{(0)}
\mid V_i=v
\right] = - (\sigma_{1D}-\sigma_{0D})\Phi^{-1}(v),
\end{equation}
is identified from the finite-mixture distribution of
$\left(\varepsilon_i^{(D)},\varepsilon_i^{(1)},\varepsilon_i^{(0)}\right)$.
Hence both $\mathop{\mathrm{MTE}}(\bar d_{\mathcal N},v,x)$ and $\mathop{\mathrm{MSE}}(d,\bar d_{\mathcal N},v,x)$ are identified functions of the model primitives.
\end{theorem}
\emph{Proof.} See Appendix \ref{BCIES-APDX_ProofProp2}.
Evaluated at mean values of the covariates \(x\), MTE would exhibit heterogeneity in treatment effects due to \(\bar{d}_{\mathcal{N}}\) if \(\delta^{(1)} - \delta^{(0)} \neq 0\). Furthermore, \(\delta^{(1)} - \delta^{(0)}\) implies the patterns of interaction effects between individual treatment and neighborhood treatment: Positive interaction (\(\delta^{(1)} - \delta^{(0)} > 0\)) means the treatment is more valuable when more of neighbors are treated. In contrast, negative interaction (\(\delta^{(1)} - \delta^{(0)} < 0\)) means the treatment is more valuable when less of neighbors are treated.
\textbf{Policy-Relevant Effects}
Policy changes typically modify eligibility rules, subsidies, or program intensity, thereby shifting the probability of treatment participation. Such changes do not affect all individuals equally: they primarily induce participation among individuals who are marginal with respect to treatment choice. In contexts of social or spatial interactions, policy changes may also alter outcomes indirectly through changes in neighborhood treatment exposure. Our aim is therefore to evaluate \emph{policy-relevant average effects per induced participant}. In particular, we decompose the impact of a policy change into two components: a \emph{direct effect}, capturing the gain for individuals who are induced into treatment by the policy, and a \emph{spillover effect}, capturing the gain generated through the induced change in neighborhood treatment exposure. This framework extends the policy-rerelevant treatment effect (PRTE) concept of \citet{heckman2005structural} to settings with network or spatial interactions. The normalization by the induced participation share ensures that all policy-relevant effects are interpreted as average gains per additional participant generated by the policy change, which is the natural metric for policy evaluation.
Let \(a\) and \(a'\) denote two policy regimes, with \(a'\) representing the more generous policy. We assume that each policy modifies the treatment-selection rule through a known transformation of the identified selection equation, while the latent resistance \(V_i\) remains invariant across policy regimes Thus, for each policy regime \(r \in \{a,a'\}\),
\[
D_i^r = \mathbbm{1}\{P^r(Z_i,X_i)\ge V_i\},
\]
where \(P^r(Z_i,X_i)\) denotes the counterfactual treatment propensity under policy \(r\). The corresponding neighborhood treatment exposure is
\[
\bar D_{\mathcal N i}^r = \sum_{j\neq i} w_{ij} D^r_j,
\]
and the realized outcome is
\[
Y_i^r = Y_i^{(D_i^r)}(\bar D_{\mathcal N i}^r).
\]
We assume that the policy change is \emph{pointwise monotone}, namely,
\[
P^{a'}(z,x) \ge P^{a}(z,x), \text{ for almost every } (z,x).
\]
Thus, policy \(a'\) weakly expands participation relative to policy \(a\).
Define the conditional share of \emph{induced participants} by
\[
\Delta P(x;a,a')
\coloneqq
\mathop{\mathrm{\mathbbm{E}}}[P^{a'}(Z_i,X_i) - P^{a}(Z_i,X_i)\mid X_i=x],
\]
and assume throughout that \(\Delta P(x;a,a')>0\). This quantity represents the expected increase in treatment participation generated by the policy among individuals with covariates \(X_i=x\).
\emph{Policy-Relevant Direct Effect (PRDE)}
The Policy-Relevant Direct Effect is defined as
\[
\mathop{\mathrm{PRDE}}(x;a,a') \coloneqq \mathop{\mathrm{\mathbbm{E}}}\left[Y_i^{(1)}(\bar D_{\mathcal N i}^{a})- Y_i^{(0)}(\bar D_{\mathcal N i}^{a}) \mid P^{a}(Z_i,X_i)<V_i\le P^{a'}(Z_i,X_i), X_i=x \right].
\]
This estimand measures the average treatment gain for individuals induced into treatment by the policy, evaluated at the neighborhood treatment exposure that would prevail under the baseline policy \(a\).
\emph{Policy-Relevant Spillover Effect (PRSE)}
The Policy-Relevant Spillover Effect is defined as
\[
\mathop{\mathrm{PRSE}}(x;a,a') \coloneqq
\frac{\mathop{\mathrm{\mathbbm{E}}}\left[Y_i^{(D_i^{a'})}(\bar D_{\mathcal N i}^{a'}) - Y_i^{(D_i^{a'})}(\bar D_{\mathcal N i}^{a}) \mid X_i=x\right]}{\Delta P(x;a,a')}.
\]
This estimand isolates the contribution of the policy-induced change in neighborhood treatment exposure while holding each individual's treatment status fixed at its value under the new policy \(a'\).
\emph{Policy-Relevant Total Effect (PRTOT)}
The overall policy effect per induced participant is
\[
\mathop{\mathrm{PRTOT}}(x;a,a') \coloneqq \frac{\mathop{\mathrm{\mathbbm{E}}}[Y_i^{a'}-Y_i^{a}\mid X_i=x]}{\Delta P(x;a,a')}.
\]
Using the decomposition,
\[
Y_i^{a'}-Y_i^{a}
=
\underbrace{
\Big(
Y_i^{(D_i^{a'})}(\bar D_{\mathcal N i}^{a'})
-
Y_i^{(D_i^{a'})}(\bar D_{\mathcal N i}^{a})
\Big)
}_{\text{spillover component}}
+
\underbrace{
\Big(
Y_i^{(D_i^{a'})}(\bar D_{\mathcal N i}^{a})
-
Y_i^{(D_i^{a})}(\bar D_{\mathcal N i}^{a})
\Big)
}_{\text{direct component}},
\]
it follows that
\[
\mathop{\mathrm{PRTOT}}(x;a,a') = \mathop{\mathrm{PRDE}}(x;a,a') + \mathop{\mathrm{PRSE}}(x;a,a').
\]
To characterize the policy-relevant direct effect, define the mean neighborhood exposure under policy \(a\) among policy-induced participants with latent resistance \(v\) by
\[
\bar d_{\mathcal N}^{a,S}(v,x)
\coloneqq
\mathop{\mathrm{\mathbbm{E}}}\!\left[
\bar D_{\mathcal N i}^{a}
\;\middle|\;
P^a(Z_i,X_i)<v\le P^{a'}(Z_i,X_i),
V_i=v,
X_i=x
\right].
\]
\begin{theorem}[Identification of Policy-Relevant Direct, Spillover, and Total Effects]
\label{Theorem3}
Suppose Assumptions \ref{Assp1}--\ref{Assp4} and the policy-counterfactual conditions stated above hold. Then,
\[
\mathop{\mathrm{PRDE}}(x;a,a')
=
\int_0^1
\mathop{\mathrm{MTE}}\!\left(
\bar d_{\mathcal N}^{a,S}(v,x),
v,
x
\right)
h_{PR}(v\mid x;a,a')\,dv,
\]
where the policy weights are
\[
h_{PR}(v\mid x;a,a')
\coloneqq
\frac{
F_{P^a\mid X}(v\mid x)
-
F_{P^{a'}\mid X}(v\mid x)
}{
\Delta P(x;a,a')
}.
\]
Moreover,
$$
\mathop{\mathrm{PRSE}}(x;a,a')
=
\frac{
\mathop{\mathrm{\mathbbm{E}}}\!\left[
\delta^{(D_i^{a'})}
\left(
\bar D_{\mathcal N i}^{a'}
-
\bar D_{\mathcal N i}^{a}
\right)
\mid X_i=x
\right]
}{
\Delta P(x;a,a')
},
$$
where $\delta^{(D_i^{a'})} \coloneqq D_i^{a'}\delta^{(1)} + \left(1-D_i^{a'}\right)\delta^{(0)}.$
Consequently,
$$
\mathop{\mathrm{PRTOT}}(x;a,a') = \mathop{\mathrm{PRDE}}(x;a,a') + \mathop{\mathrm{PRSE}}(x;a,a'),
$$
and the policy-relevant direct, spillover, and total effects are identified.
\end{theorem}
\emph{Proof.} See Appendix \ref{BCIES-APDX_ProofProp3}.
\section{Bayesian Estimation and Inference}\label{BCIES-Section3}
\subsection{Bayesian data augmentation}\label{bayesian-data-augmentation}
We conduct Bayesian inference for the structural parameters and the causal estimands defined in Section \ref{BCIES-Section2_2}. Posterior computation involves two latent-data features. First, treatment status reveals only the sign of the latent treatment-selection index \(D_i^*\). Second, only one of the two potential outcomes is observed for each unit. We therefore employ data augmentation, treating the latent selection index and the missing potential outcome as auxiliary variables.
Define \(\mathbf{P}_i= \bigl[Z_i \ X_i \bigr]\), \(\mathbf{Q}_i=\bigl[ \bar D_{\mathcal N i} \ X_i \bigr]\), \(\boldsymbol{\gamma}=\bigl[\alpha \ \beta^{(D)}\bigr]\), \(\boldsymbol{\kappa}_1=\bigl[\delta^{(1)} \ \beta^{(1)}\bigr]\), \(\boldsymbol{\kappa}_0=\bigl[\delta^{(0)} \ \beta^{(0)}\bigr]\). The Spillover Roy model can then be written as
\begin{equation}
\label{eq:GGRM-mixture-main}
\begin{aligned}
D_i^*
&=\mathbf P_i^\top\boldsymbol{\gamma}
+\varepsilon_i^{(D)},\\
Y_i^{(1)}
&=\mathbf Q_i^\top\boldsymbol{\kappa}_1
+\varepsilon_i^{(1)},\\
Y_i^{(0)}
&=\mathbf Q_i^\top\boldsymbol{\kappa}_0
+\varepsilon_i^{(0)},\\
D_i&=\mathbbm{1}\{D_i^*>0\},\\
Y_i&=D_iY_i^{(1)}+(1-D_i)Y_i^{(0)}.
\end{aligned}
\end{equation}
As in Section \ref{BCIES-Section2_2}, the joint distribution of the unobservables is represented by a finite mixture of multivariate normal distributions. Let \(c_i\in\{1,\ldots,G\}\) denote the latent mixture component for unit \(i\), with
\[
\Pr(c_i=g\mid\boldsymbol{\pi})=\pi_g,\qquad
\sum_{g=1}^G\pi_g=1,
\]
and
\begin{equation}
\label{eq:mixture-main}
\boldsymbol{\varepsilon}_i\mid c_i=g
\sim
\mathcal N(\mathbf 0,\mathbf\Sigma_g),
\end{equation}
where
\begin{equation}
\label{eq:Sigma-main}
\boldsymbol{\varepsilon}_i
= \begin{bmatrix}
\varepsilon_i^{(D)}\\
\varepsilon_i^{(1)}\\
\varepsilon_i^{(0)}
\end{bmatrix},
\qquad
\mathbf\Sigma_g=
\begin{bmatrix}
1 & \sigma_{1D,g} & \sigma_{0D,g}\\
\sigma_{1D,g} & \sigma_{1,g}^2 & \sigma_{10,g}\\
\sigma_{0D,g} & \sigma_{10,g} & \sigma_{0,g}^2
\end{bmatrix}.
\end{equation}
The normalization \(\Sigma_{g,11}=1\) fixes the scale of the latent
treatment-selection equation.
Let \(Y_i^{\mathrm{mis}}\) denote the unobserved potential outcome and
define the augmented outcome vector
\begin{equation*}
\mathbf{L}_i^* \coloneqq
\begin{bmatrix}
D_i^*\\
Y_i^{(1)}\\
Y_i^{(0)}
\end{bmatrix}
=
\begin{bmatrix}
D_i^*\\
D_iY_i+(1-D_i)Y_i^{miss}\\
D_iY_i^{miss}+(1-D_i)Y_i
\end{bmatrix},
\end{equation*}
With the corresponding block-diagonal design matrix and parameter vector
\[
\mathbf R_i=
\begin{bmatrix}
\mathbf P_i^\top & 0 & 0\\
0 & \mathbf Q_i^\top & 0\\
0 & 0 & \mathbf Q_i^\top
\end{bmatrix}, \qquad
\boldsymbol{\theta}
=
\begin{bmatrix}
\boldsymbol{\gamma}\\
\boldsymbol{\kappa}_1\\
\boldsymbol{\kappa}_0
\end{bmatrix},
\]
the augmented model is
\begin{equation}
\label{eq:augmented-model-main}
\mathbf L_i^*
=
\mathbf R_i\boldsymbol{\theta}
+
\boldsymbol{\varepsilon}_i,
\qquad
\boldsymbol{\varepsilon}_i\mid c_i=g
\sim\mathcal N(\mathbf0,\mathbf\Sigma_g).
\end{equation}
Conditional on the augmented data and mixture allocations, \eqref{eq:augmented-model-main} has a Gaussian regression representation. This representation is the basis of the posterior sampler developed below. Details of the complete-data likelihood and the conditional distributions of the augmented variables are provided in Appendix \ref{BCIES-APDX_sampling}.
\subsection{Prior specification and parameter expansion}\label{prior-specification-and-parameter-expansion}
We complete the model by assigning priors to the regression parameters, mixture probabilities, and component-specific covariance matrices. Specifically,
\begin{align}
\boldsymbol{\theta}&\sim\mathcal N(\underline{\boldsymbol{\mu}}_{\theta},\underline{\mathbf V}_{\theta}),
\\
\boldsymbol{\pi}&\sim\mathcal{D}ir(\underline{\omega}_1,\ldots,\underline{\omega}_G).
\end{align}
A complication arises in posterior simulation of the covariance matrices \(\mathbf\Sigma_g\). Because treatment depends only on the sign of \(D_i^*\), the scale of the latent selection equation is not identified and its disturbance variance is normalized to one, \(\Sigma_{g,11}=1\) (see Assumption \ref{Assp4}). Directly sampling \(\mathbf\Sigma_g\) subject to this restriction complicates covariance updating.
We address this problem using parameter expansion. For each mixture component \(g\), introduce a positive expansion parameter \(\tau_g\) and define
\begin{equation}
\mathbf A_g
\coloneqq
\mathop{\mathrm{diag}}(\tau_g,1,1),
\qquad
\widetilde{\mathbf\Sigma}_g
\coloneqq
\mathbf A_g\mathbf\Sigma_g\mathbf A_g.
\end{equation}
Unlike \(\mathbf\Sigma_g\), the expanded covariance matrix \(\widetilde{\mathbf\Sigma}_g\) is unrestricted and is assigned an inverse-Wishart prior,
\begin{equation}
\label{eq:IW-expanded-main}
\widetilde{\mathbf\Sigma}_g \sim\mathcal W^{-1}(\mathbf{I}_3,\underline{\nu}).
\end{equation}
Posterior simulation proceeds on this expanded parameter space. At each covariance update, an auxiliary value of \(\tau_g^2\) is first drawn from its conditional prior implied by equation \eqref{eq:IW-expanded-main}. This auxiliary draw rescales the selection-equation residuals and enters the expanded residual cross-product matrix. Conditional on the transformed residuals, the unrestricted covariance matrix \(\widetilde{\mathbf\Sigma}_g\) is then drawn from its inverse-Wishart conditional posterior. After this update, the scale associated with the new expanded covariance draw is \(\tau_g^{2}=\widetilde{\Sigma}_{g,11},\)
and the covariance matrix in the identified parameterization is recovered as
\begin{equation}
\label{eq:recoverSigma-main}
\mathbf\Sigma_g = \mathbf A_g^{-1}\widetilde{\mathbf\Sigma}_g\mathbf A_g^{-1},\qquad \text{where }\mathbf A_g = \mathop{\mathrm{diag}}(\tau_g,1,1)
\end{equation}
which restores the identifying normalization \(\Sigma_{g,11}=1\). The parameter expansion therefore permits standard inverse-Wishart updating on an unrestricted covariance space while enforcing the normalization through deterministic rescaling. Parameter expansion can also improve the mixing of data-augmentation algorithms for latent-variable and sample-selection models \citep{ding2014bayesian, dougan2018bayesian}. The induced prior on \((\tau_g^2,\mathbf\Sigma_g)\) and the corresponding derivations are provided in Appendix \ref{BCIES-APDX_PX}.
\subsection{Posterior computation}\label{posterior-computation}
Let \(\Theta=\left(\boldsymbol{\theta},\boldsymbol{\pi},\{\mathbf\Sigma_g\}_{g=1}^G\right)\) denote the model parameters. Augmenting the observed data with \(\mathbf D^*\), \(\mathbf Y^{\mathrm{mis}}\), and the mixture allocations
\(\mathbf c\), the posterior distribution is proportional to
\begin{equation}
\label{eq:posterior-mixture-main}
p(\Theta,\mathbf D^*,\mathbf Y^{\mathrm{mis}},\mathbf c
\mid\mathbf Y,\mathbf D) \propto
p(\mathbf Y,\mathbf D,\mathbf D^*,\mathbf Y^{\mathrm{mis}},\mathbf c
\mid\Theta)
p(\boldsymbol{\theta})
p(\boldsymbol{\pi})
\prod_{g=1}^G p(\widetilde{\mathbf\Sigma}_g).
\end{equation}
We construct a parameter-expanded Gibbs sampler that alternates between latent-data augmentation and parameter updating. Given the current parameter values, the missing potential outcomes are sampled from Gaussian conditional distributions, the latent treatment indices from Gaussian distributions truncated according to observed treatment status, and the component allocations from multinomial distributions. Conditional on the resulting complete data, the regression parameters have a Gaussian posterior and the mixture probabilities have a Dirichlet posterior. Component-specific covariance matrices are updated on the expanded scale and subsequently transformed back to the identified parameterization using \eqref{eq:recoverSigma-main}. Algorithm \ref{alg:GGRM_mixture_short} summarizes the resulting sampler. Closed-form expressions for all conditional posterior distributions and further implementation details are provided in Appendix \ref{BCIES-APDX_MCMC}.
\begin{algorithm}[H]
\caption{Parameter-expanded Gibbs sampler}
\label{alg:GGRM_mixture_short}
Initialize
$\boldsymbol{\theta}$,
$\boldsymbol{\pi}$,
$\{\mathbf\Sigma_g\}_{g=1}^G$,
and $\mathbf c$.
\For{$s=1,\ldots,S$}{
\textbf{1. Latent-data augmentation}
\Indp
Sample missing potential outcomes
$\mathbf Y^{\mathrm{mis}}$;
Sample latent treatment indices
$\mathbf D^*$ subject to
$D_i=\mathbbm 1\{D_i^*>0\}$;
Sample mixture allocations $\mathbf c$.
\Indm
\textbf{2. Parameter updates}
\Indp
Sample regression parameters $\boldsymbol{\theta}$;
Sample mixture probabilities $\boldsymbol{\pi}$;
For each $g$, sample the auxiliary expansion scale $\tau_g^2$,
update $\widetilde{\boldsymbol{\Sigma}}_g$ on the expanded scale,
and normalize to recover $\boldsymbol{\Sigma}_g$.
\Indm
}
\Return retained posterior draws of
$\boldsymbol{\theta}$,
$\boldsymbol{\pi}$,
and $\{\mathbf\Sigma_g\}_{g=1}^G$.
\end{algorithm}
\subsection{Posterior inference for causal and policy-relevant effects}\label{posterior-inference-for-causal-and-policy-relevant-effects}
The structural and policy-relevant causal quantities introduced in Section \ref{BCIES-Section2_2} are functions of the model parameters and, for policy counterfactuals, of the treatment and exposure distributions induced by alternative policy regimes. Bayesian inference for these quantities follows directly from the retained posterior draws. For each posterior draw \(\Theta^{[s]}\), we evaluate the corresponding marginal treatment and
spillover effects,
\[
{\mathop{\mathrm{MTE}}}^{[s]}(\bar d_{\mathcal N},v,x),
\qquad
{\mathop{\mathrm{MSE}}}^{[s]}(d,\bar d_{\mathcal N},x),
\]
using the expressions derived in Section \ref{BCIES-Section2_2}. Posterior means are used as point estimates, and posterior quantiles provide credible intervals. Thus, posterior uncertainty about treatment selection, outcome responses, and their dependence is propagated jointly to the heterogeneous causal effects. For a counterfactual policy regime \(\mathcal P\), we additionally construct the treatment decisions and neighborhood exposures implied by the policy at each posterior draw. Comparing these quantities across the baseline and counterfactual regimes yields posterior draws of the policy-relevant direct and spillover effects,
\[
{\mathop{\mathrm{PRDE}}}^{[s]},
\qquad
{\mathop{\mathrm{PRSE}}}^{[s]},
\qquad
{\mathop{\mathrm{PRTOT}}}^{[s]}
=
{\mathop{\mathrm{PRDE}}}^{[s]}+{\mathop{\mathrm{PRSE}}}^{[s]}.
\]
Hence, uncertainty in the structural parameters is propagated through both endogenous treatment participation and the neighborhood exposure generated by counterfactual policies.
\section{Simulation Study}\label{BCIES-Section4}
We assess the finite-sample performance of the proposed framework through a series of Monte Carlo experiments. The baseline design incorporates both selection on unobservables and spillovers, the two features that motivate the Spillover Roy model. We evaluate recovery of the structural parameters and heterogeneous marginal treatment effects and compare the correctly specified \emph{Spillover Roy model (SRM)} with a misspecified \emph{Non-Spillover Roy model (NSRM)} that omits neighborhood exposure. Additional robustness designs are reported in Appendix \ref{BCIES-APDX_Simulations}.
\subsection{Data Generating Processes}\label{BCIES-Section4_1}
For each replication, we generate five exogenous variables \(\widetilde{\mathbf X}_k\), \(k=1,\ldots,5\) independently from the standard normal distribution and set \(\mathbf X=\left[\boldsymbol{\iota}_n^\top,\widetilde{\mathbf X}^\top\right]^\top.\). The instrumental variable \(\mathbf Z\) is also generated independently from the same distribution. We construct the interaction matrix \(\mathbf W\) following the interaction design described in \citet{liu2010gmm}. The matrix is block diagonal, with each block representing a group-specific interaction network. The sample is partitioned into \(G = 30\) groups. Group sizes \(m_g\) are allowed to vary around \(n/G\): for the first \(G-1\) groups, \(m_g\) is drawn from
\[
\left\{
\lfloor n/G\rfloor-2,\ldots,\lfloor n/G\rfloor+3
\right\},
\]
and the size of the final group is chosen so that \(\sum_{g=1}^G m_g=n\). Within group \(g\), the interaction matrix \(\mathbf W_g\) is generated as follows. For each row \(i=1,\ldots,m_g\), we draw \(\tau_{ig}\) uniformly from \(\{1,2,3,4\}\) and connect unit \(i\) to the subsequent \(\tau_{ig}\) units, wrapping around the group boundary when necessary. We then symmetrize the group-specific matrices and construct
\[
\mathbf{W} \coloneqq \mathop{\mathrm{diag}}(\bf{W}_1^\top+\bf{W}_1,\ldots,\bf{W}_g^\top+\bf{W}_g)
\]
After row normalization, neighborhood exposure is given by \(\bar{\mathbf{D}}_{\mathcal{N}} =\mathbf{W}\mathbf{D}\). Treatment selection and potential outcomes follow the Spillover Roy model described in \eqref{KeyModel2}. We set
\[
\beta^{(D)} = \left[0,0.5,-0.5,0.3,-0.3,-0.2\right]^\top,
\]
\[
\beta^{(1)} = \left[2,0.4,-0.4,0.3,-0.3,0.2\right]^\top, \quad \beta^{(0)} = \left[1,0.4,-0.4,0.3,-0.3,0.2\right]^\top,
\]
with instrument strength \(\alpha = 1.5\). Spillovers are present in both potential-outcome regimes, with \(\delta^{(1)} = 1.5\) and \(\delta^{(0)} = 0.5\). The disturbance vector
\(\epsilon_i = \left[\epsilon_i^{(D)} , \epsilon_i^{(1)} , \epsilon_i^{(0)}\right]^\top\) is independently distributed across units as
\[
\begin{aligned}
\epsilon_i\overset{iid}{\sim}\mathcal N(0,\Sigma),
\qquad
\Sigma=
\begin{bmatrix}
1 & 0.9 & 0.7\\
0.9 & 1 & 0.6\\
0.7 & 0.6 & 1
\end{bmatrix},
\end{aligned}
\]
Thus, \(\sigma_{D}^2 = \sigma_{1}^2 = \sigma_{0}^2 = 1\) and \((\rho_{1D},\rho_{0D},\rho_{10}) = (0.9,0.7,0.6)\). The nonzero correlations between the selection disturbance and the potential-outcome disturbances generate selection on unobservables, while \(\rho_{10}>0\) induces positive cross-regime dependence between \(\bigl(Y^{(1)},Y^{(0)}\bigr)\).
We consider sample sizes \(n\in\{500,1{,}000,2{,}000\}\). For each simulated sample, we estimate two specifications. The \emph{SRM} includes neighborhood exposure and corresponds to the correctly specified model. The \emph{NSRM} omits \(\bar D_{\mathcal Ni}\) and therefore provides a benchmark for assessing the consequences of ignoring spillovers. For each specification, the MCMC sampler is run for \(11{,}000\) iterations, with the first \(1{,}000\) draws discarded as burn-in. The prior hyperparameters are
\[\underline{\boldsymbol{\mu}}_{\theta} = \mathbf{0}_{21};\quad \underline{\mathbf V}_{\theta} = 10^2\mathbf{I}_{21};\quad \underline{\nu} = 4;\quad \underline{\omega}_1 = \ldots = \underline{\omega}_G = 1/G.\]
Across \(N_{\mathrm{sim}}=1{,}000\) Monte Carlo replications, posterior means are used as point estimates and 95\% posterior credible intervals are used for interval estimation. We report Monte Carlo bias, root mean squared error (RMSE), and empirical coverage of the 95\% credible intervals.
\subsection{Simulation Results}\label{BCIESsection4.2}
Table \ref{tab:tab-SMCresults-Mar26-DGP2-SI} reports finite-sample performance for the structural parameters. Under the correctly specified SRM, biases are generally small, RMSEs decline with the sample size, and empirical coverage is close to the nominal 95\% level for most parameters. The results indicate that the proposed Bayesian procedure accurately recovers the principal features of the Spillover Roy model in samples of the sizes considered. The NSRM produces a markedly different pattern. By construction, it sets the spillover coefficients \(\delta^{(1)}\) and \(\delta^{(0)}\) to zero even though both are nonzero in the data-generating process. More importantly, this omission contaminates estimation of other structural parameters. In particular, the outcome coefficients and several covariance parameters exhibit persistent bias and substantial coverage distortions. These discrepancies do not disappear as \(n\) increases, consistent with misspecification bias rather than finite-sample variability.
\begin{spacing}{1.5}
\begingroup
\setlength{\tabcolsep}{3.5pt}
\begin{table}[H]
\centering
\caption{\label{tab:tab-SMCresults-Mar26-DGP2-SI}Simulation Results for Model Parameters}
\centering
\resizebox{\ifdim\width>\linewidth\linewidth\else\width\fi}{!}{
\begin{threeparttable}
\begin{tabular}[t]{cccccccccccccccc}
\toprule
\multicolumn{3}{c}{ } & \multicolumn{4}{c}{Quantities of Interest} & \multicolumn{9}{c}{Other Parameters} \\
\cmidrule(l{3pt}r{3pt}){4-7} \cmidrule(l{3pt}r{3pt}){8-16}
Model & Metric & n & $\delta^{(1)}$ & $\delta^{(0)}$ & $\delta^{(1)}-\delta^{(0)}$ & $\sigma_{1D}-\sigma_{0D}$ & $\alpha$ & $\beta^{(D)}$ & $\beta_1^{(1)}$ & $\beta_1^{(0)}$ & $\sigma_1^2$ & $\sigma_0^2$ & $\rho_{1D}$ & $\rho_{0D}$ & $\rho_{10}$\\
\midrule
& & True Value & 1.500 & 0.500 & 1.000 & 0.200 & 1.500 & 0.000 & 2.000 & 1.000 & 1.000 & 1.000 & 0.900 & 0.700 & 0.600\\
\cmidrule{1-16}
& & 500 & 0.003 & -0.002 & 0.006 & -0.040 & 0.064 & -0.020 & 0.014 & -0.001 & 0.009 & 0.007 & -0.046 & -0.006 & 0.105\\
& & 1000 & 0.007 & 0.000 & 0.007 & -0.024 & 0.028 & -0.012 & 0.005 & 0.000 & 0.002 & 0.002 & -0.024 & -0.001 & 0.110\\
& \multirow{-3}{*}{\centering\arraybackslash Bias} & 2000 & -0.006 & -0.001 & -0.006 & -0.013 & 0.017 & -0.005 & 0.010 & 0.000 & 0.000 & 0.001 & -0.012 & 0.000 & 0.110\\
& & 500 & 0.220 & 0.240 & 0.325 & 0.119 & 0.142 & 0.076 & 0.130 & 0.148 & 0.105 & 0.098 & 0.061 & 0.082 & 0.126\\
& & 1000 & 0.152 & 0.176 & 0.233 & 0.087 & 0.085 & 0.054 & 0.089 & 0.108 & 0.067 & 0.070 & 0.035 & 0.058 & 0.122\\
& \multirow{-3}{*}{\centering\arraybackslash RMSE} & 2000 & 0.108 & 0.119 & 0.163 & 0.063 & 0.059 & 0.037 & 0.064 & 0.072 & 0.048 & 0.049 & 0.023 & 0.042 & 0.119\\
& & 500 & 0.968 & 0.961 & 0.954 & 0.956 & 0.900 & 0.928 & 0.964 & 0.960 & 0.966 & 0.954 & 0.929 & 0.961 & 0.972\\
& & 1000 & 0.948 & 0.937 & 0.953 & 0.950 & 0.920 & 0.941 & 0.954 & 0.941 & 0.976 & 0.953 & 0.943 & 0.961 & 0.967\\
\multirow{-9}{*}[1\dimexpr\aboverulesep+\belowrulesep+\cmidrulewidth]{\centering\arraybackslash SRM} & \multirow{-3}{*}{\centering\arraybackslash Coverage} & 2000 & 0.942 & 0.948 & 0.956 & 0.952 & 0.921 & 0.941 & 0.954 & 0.952 & 0.979 & 0.957 & 0.958 & 0.960 & 0.960\\
\cmidrule{1-16}
& & 500 & -1.500 & -0.500 & -1.000 & -0.029 & 0.060 & -0.013 & 0.757 & 0.249 & 0.138 & 0.021 & -0.089 & -0.014 & 0.114\\
& & 1000 & -1.500 & -0.500 & -1.000 & -0.015 & 0.027 & -0.007 & 0.754 & 0.249 & 0.131 & 0.016 & -0.070 & -0.008 & 0.124\\
& \multirow{-3}{*}{\centering\arraybackslash Bias} & 2000 & -1.500 & -0.500 & -1.000 & -0.011 & 0.016 & -0.002 & 0.755 & 0.251 & 0.124 & 0.016 & -0.061 & -0.005 & 0.130\\
& & 500 & 1.500 & 0.500 & 1.000 & 0.124 & 0.144 & 0.075 & 0.763 & 0.264 & 0.180 & 0.102 & 0.103 & 0.084 & 0.132\\
& & 1000 & 1.500 & 0.500 & 1.000 & 0.091 & 0.087 & 0.054 & 0.757 & 0.257 & 0.152 & 0.073 & 0.078 & 0.058 & 0.135\\
& \multirow{-3}{*}{\centering\arraybackslash RMSE} & 2000 & 1.500 & 0.500 & 1.000 & 0.067 & 0.061 & 0.037 & 0.756 & 0.254 & 0.136 & 0.053 & 0.066 & 0.042 & 0.138\\
& & 500 & 0.000 & 0.000 & 0.000 & 0.963 & 0.891 & 0.932 & 0.000 & 0.173 & 0.805 & 0.962 & 0.676 & 0.964 & 0.987\\
& & 1000 & 0.000 & 0.000 & 0.000 & 0.957 & 0.924 & 0.941 & 0.000 & 0.016 & 0.656 & 0.944 & 0.508 & 0.969 & 0.986\\
\multirow{-9}{*}[1\dimexpr\aboverulesep+\belowrulesep+\cmidrulewidth]{\centering\arraybackslash NSRM} & \multirow{-3}{*}{\centering\arraybackslash Coverage} & 2000 & 0.000 & 0.000 & 0.000 & 0.957 & 0.922 & 0.935 & 0.000 & 0.000 & 0.422 & 0.946 & 0.296 & 0.965 & 0.963\\
\bottomrule
\end{tabular}
\begin{tablenotes}[para]
\item \textit{Notes:} This table displays the average bias (Bias), the Root Mean Squared Error (RMSE), and the coverage rate (Coverage) across $R=100$ replicates; where $\text{Bias}=R^{-1}\sum_{r=1}^R (\hat{\alpha}_r-\alpha),\text{ RMSE}=\sqrt{R^{-1}\sum_{r=1}^R (\hat{\alpha}_r-\alpha)^2},$ and $\text{ Coverage}=R^{-1}\sum_{r=1}^R \mathbbm{1}\{\alpha\in \widehat{CI}_{0.95,r}\}$. The rows contain results for models with/without spatial interference and for various sample size $n$.
\end{tablenotes}
\end{threeparttable}}
\end{table}
\endgroup
\end{spacing}
We next evaluate estimation of the marginal treatment effect over both latent resistance to treatment and neighborhood exposure. To keep the main presentation concise, Table \ref{tab:tab-SMCresults-DGP2-MTE-combined} reports results at three representative resistance values, \(v\in\{0.1,0.5,0.9\}\), for low, medium, and high neighborhood exposure, \(\bar d_{\mathcal N}\in\{0.1,0.5,0.9\}\). Complete results over \(v\in\{0.1,\ldots,0.9\}\) are reported in Supplementary Tables \ref{tab:tab-SMCresults-Mar26-DGP2-MTE1s}---\ref{tab:tab-SMCresults-Mar26-DGP2-MTE3s}.
\begin{table}[!ht]
\centering
\caption{\label{tab:tab-SMCresults-DGP2-MTE-combined}Finite-sample performance of MTE estimation across neighborhood exposure regimes.}
\centering
\fontsize{9}{11}\selectfont
\begin{threeparttable}
\begin{tabular}[t]{>{\raggedright\arraybackslash}p{1.1cm}>{\centering\arraybackslash}p{0.7cm}rrrrrrrrr}
\toprule
\multicolumn{2}{c}{ } & \multicolumn{3}{c}{$n=500$} & \multicolumn{3}{c}{$n=1000$} & \multicolumn{3}{c}{$n=2000$} \\
\cmidrule(l{3pt}r{3pt}){3-5} \cmidrule(l{3pt}r{3pt}){6-8} \cmidrule(l{3pt}r{3pt}){9-11}
Model & $v$ & Bias & RMSE & Coverage & Bias & RMSE & Coverage & Bias & RMSE & Coverage\\
\midrule
\addlinespace[0.3em]
\hline
\multicolumn{11}{l}{\textbf{Panel A. Low exposure}}\\
& 0.1 & 0.372 & 0.429 & 0.603 & 0.386 & 0.416 & 0.267 & 0.390 & 0.406 & 0.044\\
& 0.5 & 0.409 & 0.428 & 0.103 & 0.405 & 0.414 & 0.006 & 0.404 & 0.408 & 0.000\\
\multirow{-3}{*}{\raggedright\arraybackslash \hspace{1em}NSRM} & 0.9 & 0.446 & 0.482 & 0.452 & 0.424 & 0.443 & 0.154 & 0.418 & 0.428 & 0.012\\
\cmidrule{1-11}
& 0.1 & -0.036 & 0.237 & 0.961 & -0.024 & 0.176 & 0.958 & -0.008 & 0.124 & 0.951\\
& 0.5 & 0.016 & 0.168 & 0.963 & 0.006 & 0.118 & 0.961 & 0.009 & 0.085 & 0.964\\
\multirow{-3}{*}{\raggedright\arraybackslash \hspace{1em}SRM} & 0.9 & 0.068 & 0.217 & 0.974 & 0.037 & 0.149 & 0.969 & 0.026 & 0.110 & 0.953\\
\cmidrule{1-11}
\addlinespace[0.3em]
\hline
\multicolumn{11}{l}{\textbf{Panel B. Medium exposure}}\\
& 0.1 & -0.028 & 0.217 & 0.958 & -0.014 & 0.155 & 0.948 & -0.010 & 0.113 & 0.947\\
& 0.5 & 0.009 & 0.126 & 0.946 & 0.005 & 0.084 & 0.955 & 0.004 & 0.061 & 0.952\\
\multirow{-3}{*}{\raggedright\arraybackslash \hspace{1em}NSRM} & 0.9 & 0.046 & 0.187 & 0.974 & 0.024 & 0.131 & 0.970 & 0.018 & 0.097 & 0.959\\
\cmidrule{1-11}
& 0.1 & -0.034 & 0.204 & 0.963 & -0.021 & 0.148 & 0.951 & -0.010 & 0.108 & 0.949\\
& 0.5 & 0.018 & 0.114 & 0.958 & 0.009 & 0.076 & 0.961 & 0.007 & 0.055 & 0.960\\
\multirow{-3}{*}{\raggedright\arraybackslash \hspace{1em}SRM} & 0.9 & 0.070 & 0.175 & 0.968 & 0.040 & 0.121 & 0.962 & 0.024 & 0.086 & 0.961\\
\cmidrule{1-11}
\addlinespace[0.3em]
\hline
\multicolumn{11}{l}{\textbf{Panel C. High exposure}}\\
& 0.1 & -0.428 & 0.479 & 0.514 & -0.414 & 0.442 & 0.276 & -0.410 & 0.425 & 0.053\\
& 0.5 & -0.391 & 0.410 & 0.117 & -0.395 & 0.404 & 0.004 & -0.396 & 0.401 & 0.000\\
\multirow{-3}{*}{\raggedright\arraybackslash \hspace{1em}NSRM} & 0.9 & -0.354 & 0.398 & 0.610 & -0.376 & 0.397 & 0.237 & -0.382 & 0.394 & 0.030\\
\cmidrule{1-11}
& 0.1 & -0.032 & 0.247 & 0.955 & -0.018 & 0.175 & 0.960 & -0.012 & 0.129 & 0.931\\
& 0.5 & 0.020 & 0.177 & 0.948 & 0.012 & 0.122 & 0.954 & 0.005 & 0.086 & 0.951\\
\multirow{-3}{*}{\raggedright\arraybackslash \hspace{1em}SRM} & 0.9 & 0.072 & 0.220 & 0.959 & 0.043 & 0.156 & 0.960 & 0.022 & 0.106 & 0.961\\
\bottomrule
\end{tabular}
\begin{tablenotes}[para]
\item \textit{Notes:} Panels A--C report posterior bias, RMSE, and 95\% credible interval coverage for representative values of the latent resistance to treatment ($v = 0.1, 0.5, 0.9$). Results are based on $R=1{,}000$ Monte Carlo replications.
\end{tablenotes}
\end{threeparttable}
\end{table}
\vspace{-0.5\baselineskip}
The SRM performs well across all three exposure levels. Bias is small, RMSE generally declines with \(n\), and empirical coverage remains close to 95\% throughout most of the resistance distribution. The NSRM behaves differently. At low and high exposure, its MTE estimates exhibit substantial and persistent bias, and coverage deteriorates sharply as the sample size increases. At medium exposure, the misspecification is less consequential because the omitted exposure component happens to generate much smaller distortion under this design. The contrast across exposure levels illustrates an important feature of the problem: a model that ignores spillovers may appear adequate at particular exposure values while failing severely elsewhere. Increasing the sample size therefore improves precision under the correctly specified model but does not eliminate the distortions generated by omitting neighborhood exposure. The resulting undercoverage is particularly pronounced in regions of the exposure space where the omitted spillover component is economically important. These results demonstrate that correctly modeling interference is necessary for reliable inference on heterogeneous treatment effects.
This finding is also relevant for the policy analysis in our proposed framework. Policy-relevant direct effects are constructed by averaging MTEs over individuals induced into treatment under a policy change and over the neighborhood exposures generated by the baseline policy. Reliable policy evaluation therefore requires accurate estimation of the MTE over both the resistance and exposure dimensions. The simulation evidence shows that the SRM provides such recovery in the baseline design, whereas an analysis that omits spillovers can substantially distort the MTE surface. Accordingly, the results support using the estimated SRM as an input to the policy counterfactual analysis considered later in the paper.
\section{Empirical Application}\label{BCIES-Section5}
\subsection{Institutional Context and Empirical Design}\label{BCIES-Section5_1}
To demonstrate the empirical relevance of the proposed framework, we investigate the effects of the Opportunity Zones (OZ) program, a major U.S. place-based tax incentive introduced by the Tax Cuts and Jobs Act of 2017. This program offers preferential tax treatment for investments in designated census tracts with the objective of stimulating local economic activity. The OZ setting is particularly well suited to our framework because designation was potentially endogenous and its effects may extend beyond designated tracts. Eligibility for OZ designation was determined primarily using pre-program socioeconomic conditions from the 2011--2015 American Community Survey (ACS). Census tracts generally qualified if their poverty rate exceeded \(20\%\) or their median family income was below \(80\%\) of the area median income. Approximately \(40\%\) of U.S. census tracts were eligible. Importantly, eligibility did not imply designation. State governors were given substantial discretion to nominate up to \(25\%\) of eligible tracts within their states, after which the nominations were certified by the U.S. Treasury. This two-stage process-rule-based eligibility followed by discretionary selection among eligible tracts creates scope for endogenous selection into OZ designation. Our outcome of interest is growth in the number of housing units at the census-tract level. Tax incentives may stimulate construction and other investment within designated tracts, but the resulting effects need not stop at tract boundaries. Designation may generate positive spillovers if investment in an OZ raises demand for development in nearby areas, or negative spillovers if investment is reallocated from neighboring tracts toward tax-advantaged locations. Existing empirical studies report mixed evidence on the effects of the OZ program \citep{corinth2024opportunity, freedman2023jue, chen2023jue, wheeler2022locally}. These features motivate an empirical specification that allows both endogenous selection into designation and spatial spillovers.
We apply the Spillover Roy model in Section \ref{BCIES-Section2} to the OZ setting by letting \(D_i = QOZ_i, Z_i = Political_i, X_i = Demographic_i,\) and \(\bar{D}_{\mathcal{N}i} = \overline{QOZ}_i\), with the disturbance vector following the finite-mixture specification in \eqref{eq:mixture-main}. The outcome \(Y_i\) is the housing-unit growth between 2017 and 2022. Individual treatment \(QOZ_i\) is an indicator variable equal to one if the tract was designated as a Qualified Opportunity Zone and zero if it was eligible but not designated, so the analysis focuses on the designation margin among OZ-eligible tracts. Neighborhood treatment \(\overline{QOZ} = \sum_{j \neq i}QOZ_i\) is the share of neighboring tracts designated as QOZs, based on a row-normalized spatial adjacency matrix \(\mathbf{W} = (w_{ij})\). OZ eligibility and designation are obtained from the Urban Institute, and tract boundaries used to construct \(\mathbf{W}\) are obtained from the U.S. Census Bureau's TIGER/Line Shapefiles. The demographic covariates include the poverty rate, median earnings, and employment rate constructed from the ACS 2013--2017 five-year estimates and enter both the selection and outcome equations to account for observed characteristics associated with OZ designation and housing development.
We use partisan alignment between a tract's state legislative representative and the governor as the excluded variable in the treatment-selection equation. Specifically, \(Political_i\) equals one if the representative of tract \(i\) in the state's lower legislative chamber and the governor belong to the same political party, and zero otherwise. Previous studies document that political alignment is associated with the likelihood of OZ designation \citep{alm2021land, frank2022determines, eldar2022does}, supporting instrument relevance. The identifying restriction is that, conditional on the included pre-treatment tract characteristics, partisan alignment affects the housing-unit growth between 2017 and 2022 only through OZ designation but does directly affects the potential outcomes. Appendix \ref{BCIES-APDX_EmpiricalApp} provides additional details on instrument construction and sensitivity analyses.
We focus on California, for which we can assemble comprehensive tract-level data on OZ designation, housing outcomes, demographic characteristics, political affiliation, and spatial linkages. The final sample contains \(3{,}699\) OZ-eligible census tracts, comprising \(727\) designated QOZs and \(2{,}972\) eligible but non-designated tracts (Non-QOZs). Supplementary Figure \ref{fig:fig-OZ-California} displays their spatial distribution illustrating the close geographic proximity of designated and non-designated tracts. Detailed variable definitions and data sources are provided in Supplementary Tables \ref{tab-definition}--\ref{tab-source}. Designated tracts are systematically more disadvantaged along several pre-treatment socioeconomic dimensions; detailed summary statistics are reported in Supplementary Table \ref{tab:tab-OZ-summarystats}. This reinforces the importance of accounting for nonrandom selection into OZ designation.
\subsection{Estimation Results}\label{BCIES-Section5_2}
\begin{table}[!ht]
\centering
\caption{\label{tab:tab-OZ-est-results}Posterior Estimates of Key Model Parameters}
\centering
\fontsize{9}{11.5}\selectfont
\begin{threeparttable}
\begin{tabular}[t]{>{\raggedright\arraybackslash}p{5.2cm}>{\centering\arraybackslash}p{2.6cm}>{\centering\arraybackslash}p{1.8cm}>{\centering\arraybackslash}p{3.4cm}}
\toprule
& Posterior Mean & SD & 90\% Credible Interval\\
\midrule
\addlinespace[0.3em]
\multicolumn{4}{l}{\textbf{Treatment Selection}}\\
\hspace{1em}Partisan alignment ($\alpha$) & 0.162 & 0.070 & {}[0.049, 0.277]\\
\addlinespace[0.3em]
\multicolumn{4}{l}{\textbf{Neighborhood Exposure}}\\
\hspace{1em}QOZs ($\delta^{(1)}$) & 0.032 & 0.014 & {}[0.009, 0.055]\\
\hspace{1em}Non-QOZs ($\delta^{(0)}$) & 0.009 & 0.009 & {}[-0.006, 0.024]\\
\hspace{1em}$\delta^{(1)}-\delta^{(0)}$ & 0.023 & 0.016 & {}[-0.004, 0.051]\\
\addlinespace[0.3em]
\multicolumn{4}{l}{\textbf{Endogenous Selection}}\\
\hspace{1em}$\rho_{1D}$ & 0.182 & 0.134 & {}[-0.063, 0.412]\\
\hspace{1em}$\rho_{0D}$ & -0.130 & 0.052 & {}[-0.218, -0.045]\\
\addlinespace[0.3em]
\multicolumn{4}{l}{\textbf{Selection on Unobserved Gains}}\\
\hspace{1em}$\sigma_{1D}-\sigma_{0D}$ & 0.052 & 0.022 & {}[0.018, 0.092]\\
Observations & 3,699 & & \\
\bottomrule
\end{tabular}
\begin{tablenotes}[para]
\item \textit{Notes:} Posterior means, standard deviations, and 90\% credible intervals are reported for the preferred specification with pre-treatment demographic controls. $\rho_{dD}$ denotes the correlation between the treatment-selection disturbance and the potential-outcome disturbance under treatment state $d$. $\sigma_{1D}-\sigma_{0D}$ measures selection on unobserved treatment gains. Full parameter estimates and alternative specifications are reported in Supplementary Table S.X.
\end{tablenotes}
\end{threeparttable}
\end{table}
Table \ref{tab:tab-OZ-est-results} reports posterior estimates of the key parameters from our preferred specification with pre-treatment demographic controls; full parameter estimates and alternative specifications are reported in Appendix \ref{BCIES-APDX_EmpiricalApp-estimates}. Partisan alignment is positively associated with OZ designation, with a 90\% credible interval excluding zero, supporting instrument relevance. The neighborhood-treatment coefficient is positive for QOZs but smaller and imprecisely estimated for non-QOZs. The estimated dependence between the treatment-selection and potential-outcome disturbances provides evidence of endogenous selection, while the positive estimate of \(\sigma_{1D}-\sigma_{0D}\) implies indicates selection on treatment gains. We therefore next examine how the MTE varies jointly with latent resistance and neighborhood exposure.
\begin{figure}[H]
{\centering \includegraphics[width=0.8\linewidth]{./figures-BCIES/OZ_MTE_combined_figure}
}
\caption[Marginal treatment effects by latent resistance and exposure.]{Marginal Treatment Effects by Latent Resistance and Exposure. The estimated MTE declines with latent resistance $v$ at all exposure levels, indicating negative selection on gains. Higher exposure shifts the MTE upward, consistent with positive spillover effects, although the shift is modest. Shaded bands show 90\% credible intervals.}\label{fig:fig-OZ-MTE-combined}
\end{figure}
\vspace{-0.75\baselineskip}
Figure \ref{fig:fig-OZ-MTE-combined} reveals substantial heterogeneity in the effect of OZ designation. The MTE declines with latent resistance \(v\), implying that tracts more likely to be designated experience larger gains, whereas effects become negative toward the upper end of the resistance distribution. Higher neighborhood OZ exposure shifts the MTE upward, although the magnitude of this exposure-related heterogeneity is modest. Thus, treatment gains vary systematically with both endogenous selection and the surrounding treatment environment. Detailed posterior estimates of the MTE are reported in Supplementary Table \ref{tab:tab-OZ-est-MTE}.
Under the realized OZ assignment, the average direct effect on treated tracts is approximately \(4.5\) percentage points, (\(90\% \text{CI}: [1.4, 7.8]\)), while neighborhood spillovers add approximately \(1.4\) percentage points, yielding a positive average total effect of \(5.9\) percentage points. By contrast, the average spillover effect on untreated tracts is small and statistically insignificant. Full estimates of average causal effects are reported in Supplementary Table \ref{tab:tab-OZ-est-CP}. Thus, the estimated gains are concentrated primarily among designated tracts, with comparatively limited spillover benefits to non-QOZs.
\subsection{Policy Counterfactual Analysis}\label{BCIES-Section5_3}
We next evaluate counterfactual expansions of OZ designation using the policy-relevant effects defined in Section \ref{BCIES-Section2}. We consider policies that increase the baseline treatment probability \(P_a\) according to
\[
P_a^{\tau} = P_a + \tau (1 - P_a), \quad \tau \in [0,1],
\]
where \(\tau\) closes a fraction of the remaining gap between the baseline treatment probability and one. Thus, larger values of \(\tau\) induce progressively broader expansions while preserving treatment probabilities within the unit interval. Because an expansion not only changes census tracts' own designation status but also neighborhood OZ exposure, its total effect reflects both direct gains for newly treated units and indirect gains arising from changes in the surrounding treatment environment.
Figure \ref{fig:fig-OZ-PRTE-shifta} reports the policy-relevant effects across counterfactual expansions, showing pronounced diminishing returns to OZ expansion. The PRDE declines steadily with \(\tau\), becoming negative under sufficiently large expansions. This pattern follows from the declining MTE profile: broader policies induce tracts farther along the latent-resistance margin, for which expected gains from designation are progressively smaller. By contrast, the PRSE increases with policy intensity as additional designations raise neighborhood OZ exposure. These spillover gains only partially offset the declining direct gains, however, so the PRTOT also falls and eventually becomes negative. Thus, the estimated benefits of expanding OZ designation depend importantly on the scale of expansion. Posterior estimates and the corresponding shares of induced tracts are reported in Supplementary Table \ref{tab:tab-OZ-PRTE-shifta}.
\begin{figure}[H]
{\centering \includegraphics[width=0.8\linewidth]{figures-BCIES/fig-OZ-PRTE-shifta.pdf}
}
\caption{Policy-Relevant Effects under OZ Expansion. The figure reports posterior means of the policy-relevant direct (PRDE), spillover (PRSE), and total (PRTOT) effects per induced tract across counterfactual policy expansions. Shaded regions denote 90\% credible intervals.}\label{fig:fig-OZ-PRTE-shifta}
\end{figure}
\vspace{-0.75\baselineskip}
Figure \ref{fig:fig-OZ-PRSE-shifta-bygroup} further decomposes the spillover channel by treatment response type. Spillover gains are largest for policy-induced tracts and also increase for always-treated tracts as expansion raises neighborhood exposure. By contrast, estimated spillover effects for never-treated tracts remain small and imprecisely, with credible intervals including zero across the expansions considered. The estimated spillover benefits therefore accrue primarily to induced and already-treated tracts rather than broadly to tracts that remain untreated. Corresponding estimates are reported in Supplementary Table \ref{tab:tab-OZ-PRSE-shifta-bygroup}.
\begin{figure}[H]
{\centering \includegraphics[width=0.8\linewidth]{figures-BCIES/fig-OZ-PRSE-shifta-bygroup.pdf}
}
\caption{Group-Specific Spillover Effects under OZ Expansion. The figure reports posterior mean spillover effects by policy shift for always-treated units, induced entrants, and never-treated units. Shaded regions denote 90\% credible intervals.}\label{fig:fig-OZ-PRSE-shifta-bygroup}
\end{figure}
\vspace{-0.75\baselineskip}
Taken together, the counterfactual results are informative for recurring OZ designation decisions. The recent permanent extension of the program introduces new rounds of tract designation beginning in 2027, requiring states to select among eligible low-income communities. Our counterfactuals do not evaluate the new designation rules directly, but they illustrate an important trade-off relevant to such decisions. As designation expands, additional neighborhood spillovers coexist with diminishing direct gains as the policy reaches tracts with greater latent resistance; under sufficiently large expansions, the spillover gains are insufficient to offset the declining direct returns. Thus, the consequences of expanding a place-based program depend not only on how many additional areas are designated, but also on which areas are induced into treatment and how those designations alter surrounding treatment exposure.
\section{Conclusion}\label{BCIES-Section6}
This paper develops a framework for policy-relevant causal inference when treatment is endogenously selected and outcomes are subject to spillovers in a large network or spatial setting. The proposed Spillover Roy model extends the Generalized Roy framework by allowing potential outcomes to depend on both own treatment and neighborhood treatment exposure. This structure accommodates heterogeneity along the latent resistance-to-treatment margin and across exposure levels. We characterize the consequences of feasible policy changes that jointly alter treatment participation and neighborhood exposure. The resulting total policy effect decomposes into a direct effect operating through induced participation and a spillover effect operating through policy-induced changes in neighborhood exposure.
We develop a Bayesian data-augmentation approach for estimation and inference, using parameter expansion to accommodate the normalization of the latent selection equation and facilitate posterior computation. Simulations demonstrate reliable recovery of structural and heterogeneous causal effects and show that ignoring spillovers can substantially distort inference. In the application to the U.S. Opportunity Zones program, we find positive direct effects of designation on housing growth and heterogeneous treatment gains consistent with selection on gains. Spillover benefits are concentrated among designated and policy-induced tracts, whereas we find little evidence of benefits for neighboring tracts that remain untreated. Counterfactual policy experiments further indicate diminishing direct returns to program expansion, with spillover gains insufficient to offset these declines under large expansions.
Several extensions merit further study. One is to allow treatment choices themselves to interact strategically, so that policy interventions propagate through equilibrium participation responses as well as outcome spillovers. A second direction is to relax the parametric structure and develop semiparametric or nonparametric identification and inference for policy-relevant effects under endogenous selection and interference, thereby broadening the robustness and applicability of the approach.
\clearpage
\bibliography{BCIES.bib}
\clearpage