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.
103,983 characters
Robust Signal Maximization in Spillover Experiments
\title{Robust Signal Maximization in Spillover Experiments}
\author{\vspace{1.5cm}
}
\author{Kirill Borusyak\\
UC Berkeley\and Peter Hull\\
Brown\and Evan Munro\\
Chicago Booth\thanks{Contact: [email removed], peter\[email removed], and [email removed]. We thank David Atkin, Michael Best, and Eric Verhoogen for encouraging this project; we thank Gabriel Kreindler, Shuangning Li, and Davide Viviano for useful comments. OpenAI's ChatGPT contributed valuable insights. }}
\date{\vspace{0.25cm}
July 2026}
\maketitle
\vspace{0.25cm}
\begin{abstract}
\begin{singlespace} \noindent\begin{adjustwidth*}{1cm}{1cm} {\normalsize We
study the optimal design and analysis of experiments for estimating
spillover effects. Assuming a known (e.g., linear) exposure mapping,
we characterize the treatment-assignment distribution and regression-based
estimator that minimize worst-case asymptotic variance against a broad
class of distributions of unobservables. The design problem yields
an intuitive solution in which the planner trades off spillover signal
strength against diffusion of spillover variation. The analysis problem
yields a simple recentered instrumental variable estimator to best
leverage this variation. This framework produces natural solutions
in several benchmark cases---such as clustered exposure---and suggests
computationally tractable approximations for general networks, including
bipartite settings. We illustrate these new tools in semi-synthetic
experiments based on two applications from development economics.
Our approach yields large standard error reductions in both experiments,
increasing effective sample sizes by 50--100\% or more.}\end{adjustwidth*} \end{singlespace} \vfill{}
\thispagestyle{empty}
\end{abstract}
\newpage
\global\long\def\expec#1#2{\mathbb{E}_{#1}\left[#2\right]}
\global\long\def\Pr#1#2{\mathrm{Pr}_{#1}\left[#2\right]}
\global\long\def\var#1#2{\mathrm{Var}_{#1}\left[#2\right]}
\global\long\def\cov#1#2{\mathrm{Cov}_{#1}\left[#2\right]}
\global\long\def\corr#1#2{\mathrm{Corr}_{#1}\left[#2\right]}
\global\long\def\V#1#2{\mathcal{V}_{#1}\left[#2\right]}
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\def\op#1{\left\Vert #1\right\Vert _{\infty}}
\global\long\global\long\global\long\def\Frob#1{\left\Vert #1\right\Vert _{2}}
\global\long\global\long\global\long\global\long\global\long
\global\long\global\long\def\interior#1{\mathrm{int}\left(#1\right)}
\section{Introduction}
Researchers increasingly estimate spillovers in experiments. Sometimes,
an intervention is randomized to one set of units while the spillover
effects are measured in another set. In such “bipartite” experiments
the intervention and outcome units are connected via a network: e.g.,
of friendship across individuals \citep{cai_social_2015} or of supplying
relationships across firms \citep{best_spillover_2025}. In other
cases, the intervention and outcome units coincide and direct and
spillover effects are jointly estimated in a single experimental sample
\citep{Miguel2004,chaurey_social_2025}. Power is a common issue in
such studies: spillover effect estimates can be noisy, while increasing
the sample size is costly. It is therefore important to make the best
use of the experimental variation, as well as to design the randomization
to target the spillover parameter of interest.
This paper develops new tools for experimental design and estimation
when spillover effects are of primary interest. Our baseline assumption
is that the researcher adopts a linear model of spillovers, in which
the key treatment variable is the fraction or number of “friends”
(i.e., network connections) who have been randomly selected for the
intervention---perhaps augmented with weights reflecting the importance
of each connection.\footnote{Specifically, we assume the spillover treatment can be written $x_{i}=\sum_{k}w_{ik}g_{k}$
where $g_{k}$ is an indicator of unit $k$ assigned to the treatment
group. The $w_{ik}$ are fully unrestricted, and the intervention
units $k=1,\dots,K$ can differ from the outcome units $i=1,\dots,N$.
We also discuss how our approach can be used for nonlinear formula
treatments $x_{i}=X_{i}(g_{1},\dots,g_{K})$ for known $X_{i}(\cdot)$
that can implicitly depend on some $w_{ik}$ \citep{borusyak_design-based_2025}.} Such linear specifications are widespread in economics, including
in experiments (e.g., \citealp{Miguel2004,cai_social_2015,atkin_reducing_2024}).
We assume the researcher has access to the network measure at the
experimental design stage. We further assume the researcher wants
to identify the spillover effect using the experimental variation;
hence, they use some recentered estimator \citep{BH1}. It is increasingly
recognized that recentering is necessary to guarantee a causal interpretation
of spillover effect estimates without assuming that the network is
exogenous, with appropriate regression or instrumental variable (IV)
estimators employed either intuitively \citep{Miguel2004} or formally
\citep{chaurey_social_2025,atkin_reducing_2024}.
We then ask: which randomization design and which recentered estimator
are likely to yield good power? Our answer reflects two countervailing
intuitions. On the one hand, when an outcome unit $i$ is exposed
to the shocks of multiple intervention units $k$, it is desirable
to positively correlate the treatment assignments of those intervention
units. Doing so increases the variance of $i$'s spillover treatment
(i.e., the spillover “signal”). On the other hand, when multiple
outcome units are exposed to the same intervention unit, their exposures
are correlated and this can increase estimation noise if their unobservables
are also correlated; correlating assignments of multiple intervention
units exacerbates this clustering problem. This makes it more desirable
to increase the degree of independence across different exposures,
spreading experimental variation out across the sample (“isotropy”).
A likely-powerful design is thus one which maximizes spillover signals
while ensuring spillover variation is sufficiently diffuse. The estimator
choice can also help with isotropy, reweighting the data to spread
out the induced spillover variation, at the cost of a lower signal.
We formalize these intuitions in a minimax approach: we look for the
experimental design and recentered estimator that minimize the worst-case
approximate variance over a large class of error distributions. We
parameterize this class parsimoniously, by a scalar that captures
the extent to which the researcher is willing to rule out mutual dependencies,
dependence on the network, and heteroskedasticity of the errors. This
approach has three advantages. First, it reflects the reality that
researchers tend to have at best a rough sense of how errors are clustered,
especially at the experimental design stage. Second, it allows for
general network structures and yields a tractable closed-form characterization
of worst-case variance, both as a function of the estimator (holding
the design fixed) and as a function of the design (provided the optimal
estimator will be used). Although implementing our solution generally
involves numerical optimization, it is amenable to practical algorithms.
Finally, the minimax approach yields non-degenerate and intuitive
solutions in contrast to some corner solutions in the literature (e.g.,
\citet{kiefer_general_1974}). For example, independent randomization
and cluster-level randomization arise as special cases when, respectively,
there are no spillovers and when the exposure is the cluster-level
treatment saturation. When exposure is a leave-one-out average of
cluster shocks, the solution is to optimally mix between these two
designs.
We characterize the minimax-optimal design and estimator in two steps.
First, we derive the optimal recentered IV estimator for any experimental
design. Its construction entails taking the spillover treatment itself,
recentering it by its expectation over the experimental design, and
partially “whitening” it to spread out the remaining exogenous
variation. While the whitening step weakens the first stage, the induced
isotropy improves robustness against non-spherical errors. Second,
we characterize the best designs given the optimal recentered instrument
will be used. We show this optimal design satisfies natural properties,
such as a 50/50 marginal distribution for each shock and separability
across independent network blocks. The optimal IV is straightforward
to construct for any design, and we propose a computationally efficient
approximation to the optimal design; the entire procedure is collected
in Algorithm \ref{alg:full}, below. In the parameterization with
the highest robustness to adversarial errors, this procedure yields
an especially intuitive closed-form solution in which the correlation
of shock assignments is determined by the cosine similarity of exposure
to them. We show how our estimators are asymptotically normal with
a sparse exposure network---even when the errors have strong mutual
correlations---and give a simple inference procedure.
We extend these solutions in several directions, including by characterizing
minimax-optimal designs and estimators when spillover effects are
estimated alongside direct effects or other linear exposure measures.
Here the optimal design balances signal and isotropy of the residual
variation in the focal spillover treatment after accounting for the
other exposures. Other extensions include estimators with predetermined
covariates, nonlinear spillover treatments, and optimal design when
a researcher faces budget constraints in allocating shocks.
We then illustrate the power gains from using our approach, relative
to independent randomization, in semi-synthetic experiments calibrated
to the data from the bipartite experiment in \citet{cai_social_2015}
and the joint estimation of direct and spillover effects in \citet{Miguel2004}.
We find large declines in standard errors in both settings, equivalent
to increasing effective sample sizes by around 50--100\% or more
relative to the independent randomization in the original papers and
no whitening. These gains arise from both the optimized experimental
design and the use of optimal recentered instruments.
This paper contributes to several related literatures. Most directly,
we add to a growing literature on optimal experimental design under
interference. Many papers in this literature focus on specific settings,
a restricted class of designs, and a pre-specified estimator (with
different models of interference and estimands). For example, \citet{Baird2018}
and \citet{Cruces} focus on partial interference, where spillovers
are confined within clusters, and derive optimal saturation designs
for estimating direct and spillover effects, including by regression.
\citet{pouget-abadie_variance_2019} and \citet{harshaw_design_2023}
instead focus on bipartite experiments and derive optimal clustered
assignments under a linear exposure-response model with heterogeneous
treatment effects. \citet{thiyageswaran_optimal_2026} study optimal
worst-case estimation of the global average treatment effect (GATE)
over a general class of designs using the standard Horvitz-Thompson
estimator, which does not account for spillovers and is generally
biased.\footnote{Other work with general networks includes \citet{viviano_experimental_2025}
who studies two-wave experiments for estimands including average treatment
and spillover effects, using a pilot wave to estimate error variances,
and \citet{viviano_causal_2025} who choose treatment clusterings
to trade off worst-case asymptotic bias and variance when estimating
GATEs.}
Relative to this literature, our approach differs in four ways. First,
we provide a unified and tractable characterization of optimal experimental
designs for general network structures, nesting partial interference
and bipartite graphs as special cases. Second, we do not restrict
the class of designs \emph{a priori}. Third, we jointly optimize the
design and the estimator to achieve further power improvement. Like
\citet{harshaw_design_2023}, we employ recentering to get unbiased
estimators without assuming network exogeneity. Finally, we target
the parameters of a linear spillover model, rather than GATE-style
estimands. While more restrictive, such constant-effects models are
widely used in practice with evidence for linearity in different contexts
(e.g., in \citet{Miguel2004} and \citet{egger2022general}).\footnote{In the presence of heterogeneous treatment effects, our baseline estimators
identify their convex averages. This is shown in Appendix \ref{appx:HetFX},
in line with earlier results on design-based estimators \citep{borusyak_negative_2024}.} When both direct and indirect effects are included, our approach
focuses on isolating them rather than obtaining the total in a GATE-style
analysis. Our problem can also be reformulated to target the GATE;
by incorporating additional linear-exposure structure, it is expected
to deliver more power than more non-parametric GATE-targeting designs.
We also add to a literature using minimax arguments to justify different
experimental designs. In settings without interference, \citet{kallus_optimality_2021}
and \citet{Bai2021} show that complete or blockwise randomization
is minimax-optimal absent strong assumptions on error dependence.
We extend this logic to spillover estimation with recentered IV, where
we show how the minimax-optimal design depends on the allowed adversarial
clustering of errors. As in \citet{thiyageswaran_optimal_2026}, we
impose a Schatten $p$-norm constraint on the second-moment matrix
of model residuals to define our minimax problem; choosing a smaller
$p$ leads to greater robustness to cross-unit dependence and pushes
the optimal design towards more diffuse randomization. The ex-ante
choice of $p$ by the researcher is similar to choices made in the
broader literature on minimax-optimal estimators.\footnote{In non-parametric regression, for example, researchers choose the
smoothness class of the data-generating process---determining the
tradeoff between robustness and confidence interval length \citep{armstrong2018optimal}.}
The minimax approach allows us to avoid degenerate solutions that
often arise in other analyses. For example, a classic literature on
optimal regression design \citep[e.g., ][]{kiefer_general_1974} chooses
where to place observations in covariate space to minimize a functional
of the ordinary least squares (OLS) covariance matrix (e.g. D-optimality,
which minimizes its determinant; \citet{kiefer1959optimum}). This
literature assumes independent and homoskedastic errors, and optimal
designs tend to be concentrated on a few extreme points in the design
space.\footnote{Our focus on non-spherical errors echoes \citet{bickel_robustness_1979}
who recover non-degenerate designs when errors are serially correlated.} In the interference setting, the optimal design in \citet{pouget-abadie_variance_2019}
can similarly produce very coarse randomization: e.g., grouping all
units into two clusters only and assigning treatment at the cluster
level. \citet{harshaw_design_2023} employ a heuristic correction
to the optimization criterion to reduce the issue, while our minimax
approach resolves it naturally---yielding a guarantee of asymptotic
normality with a sparse exposure network and mild regularity conditions.
Finally, we add to a literature on robust and powerful spillover estimation
with recentered IV. \citet{BH1} show how recentering via knowledge
of the (quasi-)experimental design can identify spillover effects
under arbitrary endogeneity of the network, while \citet{borusyak_optimal_2026}
and \citet{borusyak_estimating_2025} show how asymptotically optimal
recentered instruments can be constructed given a particular design.
Here we optimize both the design and instrument, with our applications
showing that both levers can meaningfully improve estimator precision.
Moreover, while the implementation of the \citet{borusyak_optimal_2026}
optimal estimator is hindered by its dependence on the (likely unknown)
dependence of error terms, our minimax solution provides a parsimonious
solution robust to a wide range of error structures.
The remainder of this paper is organized as follows. In Section \ref{sec:optimal_designs}
we develop our theoretical framework and results. In Section \ref{sec:computation_clt},
we propose tractable methods of computing the optimal design and provide
a central limit theorem for the recentered IV estimator under the
feasible design. In Section \ref{sec:Applications} we illustrate
these tools in two semi-synthetic experiments. Section \ref{sec:Conclusion}
concludes. The proofs of main results are given in Appendix \ref{appx:Proofs-of-Main}.
Additional results are given in Appendix \ref{appx:Additional-Results},
with proofs in Appendix \ref{appx:Proofs-additional}.
\section{Theory\label{sec:optimal_designs}}
\subsection{Setting}
We consider estimation of the spillover effects of a set of treatment
shocks $g_{k}\in\{0,1\}$, $k=1,\dots,K$, on a set of outcomes $y_{i}$,
$i=1,\dots,N$. This might be in a bipartite setting, in which the
$K$ intervention units are distinct from the $N$ outcome units,
or in a setting like \citet{Miguel2004} in which we estimate spillovers
within a single experimental sample (with $i=k$).
We parameterize spillovers with a linear exposure model, which we
assume is correctly specified:
\begin{align}
y_{i} & =\beta x_{i}+\varepsilon_{i},\qquad x_{i}=w^{\prime}_{i}g,\label{eq:model_baseline}
\end{align}
where the spillover treatment $x_{i}$ combines the shocks $g=(g_{1},\dots,g_{K})^{\prime}$
with exposure weights $w_{i}\in\mathbb{R}^{K}$, which may or may
not add to one, and $\varepsilon_{i}$ is an unobserved error.\footnote{Here correct specification includes an implicit exclusion restriction:
that the shocks only affect outcomes through $x_{i}$. While we focus
on a linear model with constant effects, Appendix \ref{appx:HetFX}
discusses the interpretation of our proposed estimands under heterogeneous
treatment effects.} For example, \citet{cai_social_2015} estimate a spillover effect
$\beta$ from information sessions randomized to a set of $K$ farmers
in China on their peers' decisions to adopt weather insurance $y_{i}$;
here $x_{i}=w^{\prime}_{i}g$ gives the share of farmer $i$'s neighbors
who were assigned to the information session, with $w_{ik}$ being
entries of the row-normalized adjacency matrix of the peer graph.
In Section \ref{subsec:Direct-Effects} we extend the model to include
direct effects of treatment, for settings like \citet{Miguel2004}.
Below we discuss how our main results extend to spillover treatments
that are nonlinear functions of the shocks.
In most of our analysis, we condition on the exposure matrix $W=(w^{\prime}_{1},\dots,w^{\prime}_{N})^{\prime}$;
we make this conditioning implicit to simplify notation. We allow
the error vector $\varepsilon=(\varepsilon_{1},\dots,\varepsilon_{N})^{\prime}$
to be either fixed (as in a standard design-based setup) or drawn
from some distribution that can implicitly depend on $W$.\footnote{All of our results hold in both cases but the restriction on the error
distribution introduced below is more natural when $\varepsilon$
are viewed as stochastic. Nevertheless, our approach is still design-based
in the sense that the moment conditions we leverage derive their validity
from the randomization of $g$ alone.} Without loss of generality, we assume that no column of $W$ is entirely
zero.
We consider the problem of choosing an experimental design, i.e.,
a probability distribution for $g$: $\delta\in\mathcal{D}$ where
$\mathcal{D}=\Delta(\left\{ 0,1\right\} ^{K})$ and $\Delta(\cdot)$
denotes the simplex. By virtue of randomization, we have $g\protect\mathpalette{\protect\independenT}{\perp}\varepsilon$.
We then estimate $\beta$ by recentered IV \citep{BH1}: i.e., we
choose a set of instrument functions $z_{i}(\cdot):\{0,1\}^{K}\rightarrow\mathbb{R}$
where $\expec{\delta}{z_{i}(g)}=0$ for all $i$, and estimate:
\begin{align*}
\hat{\beta}\left[z\right] & =\frac{\sum^{N}_{i=1}z_{i}(g)y_{i}}{\sum^{N}_{i=1}z_{i}(g)x_{i}}.
\end{align*}
We let $\mathcal{Z}_{\delta}$ be the set of all such $z=(z_{i})^{N}_{i=1}$
for a given $\delta$.
Recentering (i.e., the mean-zero property of $z_{i}(g)$, implicitly
conditional on $W$) ensures we identify $\beta$ just using the experimental
variation in $g$: the network $W$ can be arbitrarily correlated
with the unobserved $\varepsilon$. Indeed, the moment condition $\expec{\delta}{\sum^{N}_{i=1}z_{i}\varepsilon_{i}}=0$
is satisfied for any such $z$. Recentered $z\in\mathcal{Z}_{\delta}$
can be obtained from any fixed function, $\tilde{z}_{i}$, of both
$g$ and $W$ by subtracting its (known) expectation under a given
design: i.e. $z_{i}=\tilde{z}_{i}-\expec{\delta}{\tilde{z}_{i}}$.
The set of recentered IV estimators $\hat{\beta}\left[z\right]$,
$z\in\mathcal{Z}_{\delta}$, is large, including estimators that weight
by functions of $W$ (by rescaling $z$), control for functions of
$W$ (by the Frisch--Waugh--Lovell theorem), as well as non-instrumented
estimators (e.g., OLS of $y_{i}$ on $\tilde{z}_{i}=x_{i}$ controlling
for $\expec{\delta}{\tilde{z}_{i}}$, again by the Frisch--Waugh--Lovell
theorem).\footnote{Appendix \ref{appx:non-recentered} considers estimation without recentering.
While non-recentered IVs can in principle improve on mean-squared
error grounds by introducing a small amount of bias to estimation,
we show the corresponding optimal designs can have undesirable properties
like unusually high treatment rates. We also show that when $W\mathbf{1}=1$
and an intercept is included in estimation (or, more generally, when
$W\mathbf{1}$ is linearly spanned by the included controls), there
is no loss in only considering recentered instruments since any bias
in non-recentered instruments is absorbed.}
With this setup, we ask: what design $\delta$ and instrument $z$
should we choose in order to get the most precise estimate of the
spillover effect $\beta$?
\subsection{Optimal Design and Instrument}
We follow \citet{borusyak_optimal_2026} in studying the finite-sample
approximate variance of the recentered IV estimator:
\begin{align*}
\V{\delta,\mathcal{E}}z & =\frac{\var{\delta,\mathcal{E}}{\sum^{N}_{i=1}z_{i}(g)\varepsilon_{i}}}{\expec{\delta}{\sum^{N}_{i=1}z_{i}(g)x_{i}}^{2}},
\end{align*}
where here we let $\mathcal{E}$ denote the (unknown) distribution
of $\varepsilon$, with $\V{\delta,\mathcal{E}}z=\infty$ whenever
the denominator is zero. Proposition 2 of \citet{borusyak_optimal_2026}
shows that $\V{\delta,\mathcal{E}}z$ gives a good approximation to
the appropriately scaled asymptotic variance of $\hat{\beta}[z]$
so long as the estimator is well-behaved in a particular sense, which
primarily means it converges to $\beta$ at some rate with a non-vanishing
first stage.\footnote{\label{fn:well-defined-estimator}As a just-identified IV estimator,
$\hat{\beta}\left[z\right]$ may have no finite-sample moments. More
precisely, in our setting with discrete $g$ the problems may arise
from $\hat{\beta}\left[z\right]$ not well-defined in a small-probability
event that $\sum^{N}_{i=1}z_{i}(g)x_{i}=0$.}
We look for the experimental design and instrument that minimize the
worst-case $\V{\delta,\mathcal{E}}z$ under a weak restriction on
the error distribution. Specifically, we suppose the researcher believes
$\mathcal{E}$ belongs to some class of distributions $\mathcal{F}_{p}(\sigma)$
and solve the minimax problem:
\begin{align}
\delta^{*},z^{*} & \in\arg\min_{\delta\in\mathcal{D},z\in\mathcal{Z}_{\delta}}\max_{\mathcal{E}\in\mathcal{F}_{p}(\sigma)}\V{\delta,\mathcal{E}}z.\label{eq:minimax}
\end{align}
The class of error distributions, parameterized by $p\in[1,\infty]$
and $\sigma>0$, is given by:
\begin{align}
\mathcal{F}_{p}(\sigma) & =\left\{ \mathcal{E}\colon\left\Vert \expec{\mathcal{E}}{\varepsilon\varepsilon'}\right\Vert _{p}\le N^{1/p}\sigma^{2}\right\} ,\label{eq:F_p}
\end{align}
where $\left\Vert \cdot\right\Vert _{p}$ is the Schatten-$p$ matrix
norm defined for positive semi-definite matrices as\footnote{For general matrices, it is defined as $\left(\operatorname{tr}\left(\sqrt{A'A}\right)^{p}\right)^{1/p}$.
Schatten norms should not be confused with other matrix norms that
can be denoted in the same way.}
\[
\left\Vert A\right\Vert _{p}=\left(\operatorname{tr}\left(A^{p}\right)\right)^{1/p}.
\]
Three important special cases are the nuclear norm (a.k.a., the trace
norm) with $p=1$, the Frobenius norm with $p=2$, and the operator
norm with $p=\infty$. Denoting the eigenvalues of $\expec{\mathcal{E}}{\varepsilon\varepsilon'}$
by $\lambda_{1},\dots,\lambda_{N}\ge0$, we can rewrite (\ref{eq:F_p})
as an upper bound on their power mean:
\[
\mathcal{F}_{p}(\sigma)=\left\{ \mathcal{E}\colon\left(\frac{1}{N}\sum^{N}_{i=1}\lambda^{p}_{i}\right)^{1/p}\le\sigma^{2}\right\} ,
\]
with the limit $\mathcal{F}_{p}(\sigma)=\left\{ \mathcal{E}\colon\max\left\{ \lambda_{1},\dots,\lambda_{N}\right\} \le\sigma^{2}\right\} $
for $p=\infty$.
\begin{figure}[t]
\caption{Boundary of $\mathcal{F}_{p}(\sigma)$ for $\sigma^{2}=1$, $N=2$,
for Different Choices of $p$\label{fig:Fp_illustration}}
\begin{centering}
\includegraphics[width=0.3\paperwidth]{Illustrations/power_countours_p}\smallskip{}
\par\end{centering}
{\small\emph{Notes}}{\small : This figure shows boundaries of $\mathcal{F}_{p}(\sigma)$
for different Schatten-$p$ norms, with $N=2$ observations and $\sigma=1$,
in terms of the two eigenvalues of the $\expec{\mathcal{E}}{\varepsilon\varepsilon'}$
matrix. The $(1,1)$ point corresponds to spherical errors.}{\small\par}
\end{figure}
To interpret $\mathcal{F}_{p}(\sigma)$, note first that spherical
(i.e., homoskedastic and mutually uncorrelated) errors with mean zero
(implicitly conditional on $W$) and variance $\sigma^{2}$, i.e.
$\expec{\mathcal{E}}{\varepsilon\varepsilon'}=\sigma^{2}I_{N}$, are
at the boundary of this set for any $p$. Different choices of $p$
then place different weights on the correlations among the errors.
When $p=1$, only the average second moment of $\varepsilon_{i}$
is constrained, allowing fully unrestricted correlations across observations.
At the other extreme, when $p=\infty$, correlations are strongly
penalized such that a worst-case scenario in (\ref{eq:minimax}) is
always with spherical errors. Intermediate values of $p$ correspond
to smaller penalties on correlations. Figure \ref{fig:Fp_illustration}
provides a visualization for $N=2$ by plotting the boundary of $\mathcal{F}_{p}(\sigma)$
in terms of the two eigenvalues of $\expec{\mathcal{E}}{\varepsilon\varepsilon'}$;
increasing the correlation across errors makes the eigenvalues more
asymmetric without changing their sum; see also Appendix \ref{appx:Equicorrelated}
for analytical results with equicorrelated $\varepsilon_{i}$. Overall,
$p$ captures the extent to which the researcher is willing to rule
out strong correlations across the errors when designing the experiment.
Similarly, deviations from homoskedasticity and from network exogeneity
(i.e., the zero conditional mean of the errors) are penalized more
strongly when $p$ is larger.
\paragraph{Optimal Instrument.}
Solving the problem (\ref{eq:minimax}) backwards, we first characterize
the optimal instrument for a given design $\delta$:
\begin{thm}
\label{prop:OptimalIV}Fix $\delta$ and let
\begin{align*}
\tilde{x}_{\delta} & =W\tilde{g}_{\delta},\qquad\tilde{g}_{\delta}=g-\expec{\delta}g,\\
S_{\delta} & =\var{\delta}x=W\var{\delta}gW^{\prime}.
\end{align*}
Suppose $S_{\delta}\neq0$. Then for any $p\in[1,\infty]$ and $\sigma>0$,
the optimal instrument is:
\begin{align*}
z^{*}_{\delta}(g) & \propto\left(S^{\dag}_{\delta}\right)^{\frac{1}{p+1}}\tilde{x}_{\delta}\in\arg\min_{z\in\mathcal{Z}_{\delta}}\max_{\mathcal{E}\in\mathcal{F}_{p}(\sigma)}\V{\delta,\mathcal{E}}z,
\end{align*}
where $\dag$ denotes the pseudoinverse. The worst-case approximate
variance is:
\[
\max_{\mathcal{E}\in\mathcal{F}_{p}(\sigma)}\V{\delta,\mathcal{E}}{z^{*}_{\delta}}=N^{1/p}\sigma^{2}\left(\operatorname{tr}(S^{p/(p+1)}_{\delta})\right)^{-(p+1)/p}.
\]
\end{thm}
\noindent The Appendix \ref{appx:Proof-Thm1} proof provides an explicit
characterization of the worst-case error distribution. In this proposition
and for the rest of the paper, we use the convention that $1/\infty=0$,
$p/(p+1)=1$ for $p=\infty$, and $A^{0}=A^{\dag}A$ for a matrix
$A$.
Intuitively, $z^{*}_{\delta}$ takes the spillover treatment $x$,
recenters it by $\expec{\delta}x$, and then reweights the resulting
$\tilde{x}_{\delta}$ inversely by a fractional power of $\var{\delta}x$,
$S^{1/(p+1)}_{\delta}$. The goal in reweighting is to spread out
the useful variation in $\tilde{x}_{\delta}$ across different directions
that could be adversarially targeted by clustering in $\varepsilon$.
When $p=1$ and $\var{\delta}x$ is an invertible matrix, this fully
“whitens” $\tilde{x}_{\delta}$ (i.e., rotates it to make $\var{\delta}{z^{*}_{\delta}}=I$)
in order to guard against the fully-unconstrained clustering of $\varepsilon$,\footnote{When $\var{\delta}x$ is not invertible, e.g. because $\delta$ perfectly
correlates $x_{i}$ for several $i$, $z^{\ast}$ preserves that dependence.} while as $p$ increases we do less whitening as errors with strong
mutual correlations are excluded. As a result, at $p=\infty$, the
optimal IV is the unwhitened $\tilde{x}_{\delta}$.\footnote{Theorem \ref{prop:OptimalIV} technically only shows this for invertible
$\var{\delta}x$; the general optimality of $\tilde{x}_{\delta}$
for $p=\infty$, i.e. when $\mathcal{F}_{p}(\sigma)$ bounds the operator
norm of the error second-moment matrix, was previously shown in \citet[Lemma 2]{borusyak_optimal_2026}.}
\paragraph{Optimal Design.}
We next characterize the optimal design, assuming the optimal instrument
will be used for estimation:
\begin{prop}
\label{prop:OptimalDesign}Define $\tilde{x}_{\delta}$ and $S_{\delta}$
as in Theorem \ref{prop:OptimalIV} and suppose there is $\delta\in\mathcal{D}$
such that $S_{\delta}\ne0$. Then
\[
\delta^{*},z^{*}\in\arg\min_{\delta\in\mathcal{\mathcal{D}},z\in\mathcal{Z}_{\delta}}\max_{\mathcal{E}\in\mathcal{F}_{p}(\sigma)}\V{\delta,\mathcal{E}}z.
\]
is solved by
\begin{align*}
\delta^{*} & \in\arg\max_{\delta\in\mathcal{D}}\operatorname{tr}\left(S^{\frac{p}{p+1}}_{\delta}\right)
\end{align*}
together with the optimal instrument $z^{\ast}=z^{\ast}_{\delta^{\ast}}$
from Theorem \ref{prop:OptimalIV}.
\end{prop}
\noindent This result follows directly from Theorem \ref{prop:OptimalIV}.
For $p=1$, the optimal design maximizes $\operatorname{tr}\left(\var{\delta}x^{1/2}\right)$.
At $p=\infty$, the sum of variances, $\operatorname{tr}\left(\var{\delta}x\right)$,
is maximized instead.
In what follows we avoid $p=\infty$ as the optimal design problem
has a degenerate rank-1 solution. Specifically, Appendix \ref{appx:Rank-1-Designs}
shows that there exists $\delta^{*}$ that puts all the mass on just
one pair of complementary assignments $g^{*}\in\left\{ 0,1\right\} ^{K}$
and $\mathbf{1}-g^{*}$. Moreover, if $W$ is non-negative, one solution
corresponds to $g^{\ast}=\mathbf{1}$, i.e. treating all intervention units
with probability 0.5 and treating none otherwise. While this is not
problematic under spherical errors, it is highly undesirable when
errors can be correlated since the estimator may not be consistent.
Choosing $p<\infty$ guards against this possibility and helps with
consistency, as we show in Section \ref{subsec:clt-inference}.\footnote{Another problem with degenerate designs is that the asymptotic approximation
that makes the approximate variance a useful concept may fail. Moreover,
the estimator itself may not be well-defined with a probability that
is not asymptotically small, as $\sum_{i}z_{i}(g)x_{i}=0$ whenever
$g=\boldsymbol{0}$; see footnote \ref{fn:well-defined-estimator}.}
A\textbf{ }further result helps characterize optimal designs:
\begin{prop}
\label{prop:DesignProperties}There is a $\delta^{*}$ from Proposition
\ref{prop:OptimalDesign} with $\expec{\delta^{\ast}}{g_{k}}=0.5$
for all $k$. Further, if $W$ is block-diagonal after a joint partition
of rows and columns, there is such a $\delta^{*}$ with independent
subvectors of $g$ across blocks and the optimal design can be solved
separately block-by-block.
\end{prop}
\noindent Intuitively, any design with \emph{$\expec{\delta^{\ast}}{g_{k}}\ne0.5$}
can be symmetrized by averaging it with its complement and this can
never hurt the criterion. Similarly, any across-block dependence in
$g_{k}$ can get removed without hurting the criterion, since $S_{\delta^{*}}=W\var{\delta^{\ast}}gW^{\prime}$.\footnote{While $\mathcal{F}_{p}(\sigma)$ allows for error dependence even
across blocks, Proposition \ref{prop:DesignProperties} implies that
the same design would be optimal even if one were to \emph{a priori
}rule out cross-block dependence.} This result is generalized in Appendix \ref{appx:Symmetry-Normalization}:
whenever $W$ is invariant under a group of relabelings over rows
and columns, there exists a $\delta^{*}$ which is invariant under
the same group. In particular, if $W$ is invariant to permutations
of rows and columns in a cluster, one may restrict attention to exchangeable
designs within that cluster.
\subsection{Special Cases\label{subsec:Special-Cases}}
We now illustrate the results of Theorem \ref{prop:OptimalIV} and
Proposition \ref{prop:OptimalDesign} in several insightful special
cases.
\paragraph{No Spillovers.}
A useful benchmark is the standard experimental case with $N=K$ and
$W=I$. Since $W$ is diagonal, Proposition \ref{prop:DesignProperties}
immediately implies that independent random assignment, $g_{i}\stackrel{iid}{\sim}Bernoulli(0.5)$,
is optimal in this case, and the optimal instrument is just $z^{*}_{i}(g)\propto g_{i}-0.5$.\footnote{Note that complete randomization, with exactly half of the units treated,
is not optimal for $p<\infty$ as it induces a slight negative correlation
across assignments which can be exploited by the adversary choosing
the errors.}
\paragraph{Group-Specific Shocks.}
This and remaining cases are visualized in Figure \ref{fig:SpecialCases}.
Consider a simple bipartite graph in which observations are partitioned
into $J$ groups, $C_{1},\dots,C_{J}$, with shocks assigned at the
group level: $K=J$ and $x_{i}=g_{j(i)}$ corresponding to $w_{i}=e_{j(i)}$,
where $e_{j}$ is the $j$th standard basis in $\mathbb{R}^{K}$ and
$j(i)\in\left\{ 1,\dots,K\right\} $ gives the group of observation
$i$. Again Proposition \ref{prop:DesignProperties} implies independent
$g_{j}\stackrel{iid}{\sim}Bernoulli(0.5)$ is optimal. The optimal
instrument is now proportional to $N^{-1/(1+p)}_{j(i)}\cdot(g_{j(i)}-0.5)$
where $N_{j}=\left|C_{j}\right|$. Thus when $p=1$ the estimator
down-weights large groups by a factor of $1/\sqrt{N_{j(i)}}$ to ensure
the random variation is equally spread out; not doing so would make
the estimator noisy if the errors are very correlated within clusters.
With higher $p$ less of this reweighting is needed, with no reweighting
in the limit $p=\infty$.
\begin{figure}
\caption{Graphs in Special Cases\label{fig:SpecialCases}}
\begin{centering}
\begin{tabular}{ccc}
(a) Group-specific shocks\medskip{}
& & (b) Group-average shocks\tabularnewline
\begin{tikzpicture}[
>=Stealth,
obs/.style={circle, fill=black, inner sep=0pt, minimum size=4.5pt},
grp/.style={circle, draw, line width=0.8pt, inner sep=1.2pt, minimum size=18pt, font=\small},
edge/.style={thin}
]
\node[obs] (i1) at (0,3.0) {};
\node[left=3pt of i1] {\small $i$};
\node[obs] (i2) at (0,2.4) {};
\node[obs] (i3) at (0,1.8) {};
\node[obs] (i4) at (0,1.2) {};
\node[obs] (i5) at (0,0.3) {};
\node[obs] (i6) at (0,-0.3) {};
\node[obs] (i7) at (0,-1.2) {};
\node[obs] (i8) at (0,-1.8) {};
\node[grp] (j1) at (4,2.1) {$C_1$};
\node[grp] (j2) at (4,0.3) {$C_2$};
\node[grp] (j3) at (4,-1.5) {$C_3$};
\draw[edge] (i1) -- (j1);
\draw[edge] (i2) -- (j1);
\draw[edge] (i3) -- (j1);
\draw[edge] (i4) -- (j1);
\draw[edge] (i5) -- (j2);
\draw[edge] (i6) -- (j2);
\draw[edge] (i7) -- (j3);
\draw[edge] (i8) -- (j3);
\end{tikzpicture}\vspace{1cm}
& & \begin{tikzpicture}[
>=Stealth,
obs/.style={circle, fill=black, inner sep=0pt, minimum size=4.5pt},
edge/.style={thin}
]
\node[obs] (l1) at (0,3.0) {};
\node[obs] (l2) at (0,2.4) {};
\node[obs] (l3) at (0,1.8) {};
\node[obs] (l4) at (0,1.2) {};
\node[obs] (l5) at (0,0.3) {};
\node[obs] (l6) at (0,-0.3) {};
\node[obs] (l7) at (0,-1.2) {};
\node[obs] (l8) at (0,-1.8) {};
\node[left=4pt of l1] {\small $i$};
\node[obs] (r1) at (4,3.0) {};
\node[obs] (r2) at (4,2.4) {};
\node[obs] (r3) at (4,1.8) {};
\node[obs] (r4) at (4,1.2) {};
\node[obs] (r5) at (4,0.3) {};
\node[obs] (r6) at (4,-0.3) {};
\node[obs] (r7) at (4,-1.2) {};
\node[obs] (r8) at (4,-1.8) {};
\node[right=4pt of r1] {\small $k$};
\foreach \a in {1,2,3,4}
\foreach \b in {1,2,3,4}
\draw[edge] (l\a) -- (r\b);
\foreach \a in {5,6}
\foreach \b in {5,6}
\draw[edge] (l\a) -- (r\b);
\foreach \a in {7,8}
\foreach \b in {7,8}
\draw[edge] (l\a) -- (r\b);
\end{tikzpicture}\tabularnewline
(c) Leave-out averages\medskip{}
& & (d) Spatial spillovers\tabularnewline
\begin{tikzpicture}[
>=Stealth,
obs/.style={circle, fill=black, inner sep=0pt, minimum size=4.5pt},
edge/.style={thin}
]
\node[obs] (l1) at (0,3.0) {};
\node[obs] (l2) at (0,2.4) {};
\node[obs] (l3) at (0,1.8) {};
\node[obs] (l4) at (0,1.2) {};
\node[obs] (l5) at (0,0.3) {};
\node[obs] (l6) at (0,-0.3) {};
\node[obs] (l7) at (0,-1.2) {};
\node[obs] (l8) at (0,-1.8) {};
\node[left=4pt of l1] {\small $i$};
\node[obs] (r1) at (4,3.0) {};
\node[obs] (r2) at (4,2.4) {};
\node[obs] (r3) at (4,1.8) {};
\node[obs] (r4) at (4,1.2) {};
\node[obs] (r5) at (4,0.3) {};
\node[obs] (r6) at (4,-0.3) {};
\node[obs] (r7) at (4,-1.2) {};
\node[obs] (r8) at (4,-1.8) {};
\node[right=4pt of r1] {\small $k$};
\draw[edge] (l1) -- (r2);
\draw[edge] (l1) -- (r3);
\draw[edge] (l1) -- (r4);
\draw[edge] (l2) -- (r1);
\draw[edge] (l2) -- (r3);
\draw[edge] (l2) -- (r4);
\draw[edge] (l3) -- (r1);
\draw[edge] (l3) -- (r2);
\draw[edge] (l3) -- (r4);
\draw[edge] (l4) -- (r1);
\draw[edge] (l4) -- (r2);
\draw[edge] (l4) -- (r3);
\draw[edge] (l5) -- (r6);
\draw[edge] (l6) -- (r5);
\draw[edge] (l7) -- (r8);
\draw[edge] (l8) -- (r7);
\end{tikzpicture} & & \raisebox{-1.5cm}[0pt][0pt]{
\begin{tikzpicture}[
x=1cm,y=1cm,
line cap=round,
line join=round,
scale=0.90,
transform shape
]
\draw[line width=0.35pt,draw=black!60,dash pattern=on 2pt off 1.5pt]
(0,0) circle (2.0);
\foreach \j in {0,...,9} {
\coordinate (v\j) at ({90-360*\j/10}:2.0);
}
\foreach \j in {1,...,9} {
\pgfmathtruncatemacro{\jmone}{mod(\j-1+10,10)}
\pgfmathtruncatemacro{\jmtwo}{mod(\j-2+10,10)}
\pgfmathtruncatemacro{\jpone}{mod(\j+1,10)}
\pgfmathtruncatemacro{\jptwo}{mod(\j+2,10)}
\pgfmathsetmacro{\ang}{90-360*(\j)/10}
\begin{scope}[draw=black!20,line width=0.45pt]
\draw (v\j) .. controls ($ (0,0)!2.08!(v\j) $) and ($ (0,0)!1.62!(v\jmtwo) $) .. (v\jmtwo);
\draw (v\j) .. controls ($ (0,0)!1.72!(v\j) $) and ($ (0,0)!1.40!(v\jmone) $) .. (v\jmone);
\draw (v\j) .. controls ($ (0,0)!1.72!(v\j) $) and ($ (0,0)!1.40!(v\jpone) $) .. (v\jpone);
\draw (v\j) .. controls ($ (0,0)!2.08!(v\j) $) and ($ (0,0)!1.62!(v\jptwo) $) .. (v\jptwo);
\coordinate (loopcenter) at ($ (0,0)!1.75!(v\j) $);
\draw (v\j)
.. controls ($ (loopcenter)+(\ang+90:0.60) $)
and ($ (loopcenter)+(\ang-90:0.60) $)
.. (v\j);
\end{scope}
}
\pgfmathtruncatemacro{\jmone}{mod(0-1+10,10)}
\pgfmathtruncatemacro{\jmtwo}{mod(0-2+10,10)}
\pgfmathtruncatemacro{\jpone}{mod(0+1,10)}
\pgfmathtruncatemacro{\jptwo}{mod(0+2,10)}
\pgfmathsetmacro{\ang}{90-360*(0)/10}
\begin{scope}[draw=black,line width=0.45pt]
\draw (v0) .. controls ($ (0,0)!2.08!(v0) $) and ($ (0,0)!1.62!(v\jmtwo) $) .. (v\jmtwo);
\draw (v0) .. controls ($ (0,0)!1.72!(v0) $) and ($ (0,0)!1.40!(v\jmone) $) .. (v\jmone);
\draw (v0) .. controls ($ (0,0)!1.72!(v0) $) and ($ (0,0)!1.40!(v\jpone) $) .. (v\jpone);
\draw (v0) .. controls ($ (0,0)!2.08!(v0) $) and ($ (0,0)!1.62!(v\jptwo) $) .. (v\jptwo);
\coordinate (loopcenter) at ($ (0,0)!1.75!(v0) $);
\draw (v0)
.. controls ($ (loopcenter)+(\ang+90:0.60) $)
and ($ (loopcenter)+(\ang-90:0.60) $)
.. (v0);
\end{scope}
\foreach \j in {0,...,9} {
\fill (v\j) circle (1.9pt);
}
\fill (v0) circle (2.5pt);
\end{tikzpicture}
}\tabularnewline
\end{tabular}
\par\end{centering}
\begin{centering}
\bigskip{}
\bigskip{}
\par\end{centering}
{\small\emph{Notes}}{\small : This figure shows network graphs in the
special cases considered in Section \ref{subsec:Special-Cases}.}{\small\par}
\end{figure}
\paragraph{Group-Average Shocks.}
Now consider a case in which observations are partitioned into $J$
groups $C_{1},\dots,C_{J}$, but the shocks $g_{i}$ are assigned
to individuals ($K=N$) and fully propagate through the cluster-average:
$x_{i}=\bar{g}_{j(i)}$ for $\bar{g}_{j}=\frac{1}{N_{j}}\sum_{i\in C_{j}}g_{i}$,
which corresponds to $W=\operatorname{diag}(P_{N_{1}},\dots,P_{N_{J}})$ for $P_{n}=\mathbf{1}_{n}\mathbf{1}^{\prime}_{n}/n$.
This is isomorphic to the previous case, with shocks that are perfectly
correlated within groups and independent across groups being optimal:
$g_{i}=\check{g}_{j(i)}$ with $\check{g}_{j}\stackrel{iid}{\sim}Bernoulli(0.5)$.
The optimal instrument is again proportional to $N^{-1/(1+p)}_{j(i)}(\check{g}_{j(i)}-0.5)$.
Intuitively, there is no need to vary the shocks within groups since
treatment is at the group level; perfect within-group shock correlation
ensures the maximal variability of the spillover treatment while we
still reweight by $N^{-1/(1+p)}_{j(i)}$ to more evenly spread out
this variation.
\paragraph{Leave-Out Averages.}
Suppose in the previous case we instead have $N_{j}>1$ for all $j$
and $x_{i}=\bar{g}_{-i}=\frac{1}{N_{j(i)}-1}\sum_{k\in C_{j(i)},k\neq i}g_{k}$.
Appendix \ref{appx:Leave-Out-Averages} shows that an optimal design
mixes between group-specific iid $Bernoulli(0.5)$ shocks with probability
$\psi_{j}=(1-r_{j})/(1+(N_{j}-1)r_{j})$ for $r_{j}=(N_{j}-1)^{-2p}$,
and independent individual-level $Bernoulli(0.5)$ shocks with probability
$1-\psi_{j}$, independently across groups.\footnote{Note that this is just one way to induce the optimal $\var{\delta^{*}}x$.
Other designs may deliver the same variance matrix and are therefore
also optimal.} This design is close to a more conventional saturation design (studied,
e.g., in \citet{Baird2018}) with saturation rates of $0$, $0.5$,
and $1$, except with independent (rather than complete) randomization
within groups. The optimal instrument is now proportional to $\left(1+(N_{j(i)}-1)\psi_{j(i)}\right){}^{-1/(1+p)}\cdot\left((\bar{g}_{j(i)}-\frac{1}{2})+(N_{j(i)}-1)(\bar{g}_{j(i)}-g_{i})\right)$.
It is easy to verify that it is optimal to correlate the shocks more
in larger groups: $\psi_{j}$ increases in $N_{j}$, from $\psi_{j}=0$
when $N_{j}=2$ to $\psi_{j}\to1$ when $N_{j}\to\infty$. Intuitively,
in the smallest group of size $N_{j}=2$, we are effectively back
in the no-spillover case, while with $N_{j}\to\infty$ the leave-out
average is equivalent to the group average from the previous case.
For intermediate sizes, within-group correlation increases in $p$:
as before, it is higher when less clustering is allowed in $\varepsilon$.
\paragraph{Spatial Spillovers.}
Finally, suppose the $N=K$ units are equal-spaced points on a circle,
with treatment being an average of shocks to units within the distance
of $L$: $x_{i}=\frac{1}{2L+1}\sum^{i+L}_{k=i-L}g_{k}$, where index
summation is understood modulo $N$. Appendix \ref{appx:Spatial-Spillovers}
shows that, for $p=1$, the optimal design correlates shock assignments
at the same spatial scale as the exposure mapping. Specifically, the
correlation between assignments at circular distance $d$ is $\max\left\{ 1-d/(2L+1),0\right\} $,
linearly decaying in $d$. The design can be implemented by splitting
the circle into blocks of $2L+1$ consecutive units, randomly shifting
this partition, and assigning shocks at the block level. The optimal
instrument is $z^{\ast}=W^{\dag}\tilde{g}$. Under the shifted-block
implementation of the design, it takes values $(2L+1)\tilde{g}_{i}-2L\bar{\tilde{g}}$
for units at the center of a block and $\bar{\tilde{g}}$ for all
other units, where $\bar{\tilde{g}}$ is the average of all recentered
shocks. Intuitively, the IV undoes the smoothing done by spatial averaging
to mitigate the correlations that could otherwise be exploited by
the adversary. These results can be extended to a two-dimensional
grid, with $x_{i}$ averaging the shocks within a certain distance
in each coordinate.
\subsection{\label{subsec:Direct-Effects}Multiple Exposures}
We now consider a model in which the shocks $g$ affect the outcomes
through multiple linear exposure measures. For brevity, we consider
two exposures but all results generalize immediately to three or more.
For a known fixed $N\times K$ matrix $U=(u_{i}')^{N}_{i=1}$, we
now have:
\begin{align}
y_{i} & =\beta x_{i}+\tau h_{i}+\varepsilon_{i},\label{eq:multiple_exposure}\\
x_{i} & =w^{\prime}_{i}g,\qquad h_{i}=u^{\prime}_{i}g.\nonumber
\end{align}
Our goal is still to minimize the worst-case approximate estimation
variance of the effect of $x_{i}$.
While we consider a general formulation, this setting is especially
relevant when the intervention and outcome units are the same and
$\tau$ corresponds to the direct effect of the intervention:
\begin{equation}
y_{i}=\beta x_{i}+\tau g_{i}+\varepsilon_{i}\label{eq:model_direct}
\end{equation}
which corresponds to $U=I_{N}$.\footnote{Another class of applications is to multiplex networks (cf. \citet{zenou2025peer}).}
In this context, there are three reasons why one would focus on the
estimation power of $\beta$. First, the researcher may be most interested
in the spillover effect, but still wants to acknowledge the presence
of the direct effect. Second, even if direct and spillover effects
are of equal interest, the spillover effect is usually more challenging
to estimate. A design that yields sufficient power for it would usually
deliver precise estimates of the direct effect, too. Finally, if the
researcher is interested in the Global Average Treatment Effect (GATE),
defined as the effect of switching from $g=\boldsymbol{0}_{K}$ to
$g=\boldsymbol{1}_{K}$ via both direct and indirect effects, that
problem can be reformulated as the one we solve here.\footnote{Specifically, the GATE equals $\gamma=\bar{w}\beta+\bar{u}\tau$ for
known constants $\bar{w}=\frac{1}{N}\sum_{i,k}w_{ik}$ and $\bar{u}=\frac{1}{N}\sum_{i,k}u_{ik}$.
One can then rewrite (\ref{eq:multiple_exposure}) as $y_{i}=\gamma\frac{x_{i}}{\bar{w}}+\tau\left(h_{i}-\frac{\bar{u}}{\bar{w}}x_{i}\right)+\varepsilon_{i}$
and apply our results to this reformulation, optimizing estimation
power for $\gamma$ while also including $h_{i}-\frac{\bar{u}}{\bar{w}}x_{i}$
in the model.}
To see how Theorem \ref{prop:OptimalIV} and Propositions \ref{prop:OptimalDesign}--\ref{prop:DesignProperties}
generalize, collect the parameters and explanatory variables in $\theta=(\beta,\tau)^{\prime}$
and $\boldsymbol{x}_{i}=(x_{i},h_{i})^{\prime}$. Absent restrictions
on how the error term can correlate with $W$ and $U$, we again use
recentered IV. For some recentered instrument vector $\boldsymbol{z}(\cdot)=(\boldsymbol{z}_{i}(\cdot))^{N}_{i=1}$,
where $\boldsymbol{z}_{i}(\cdot):\{0,1\}^{K}\rightarrow\mathbb{R}^{2}$
and $\expec{\delta}{\boldsymbol{z}_{i}(g)}=0$, consider the estimator:
\begin{align*}
\hat{\theta}\left[\boldsymbol{z}\right] & =\left(\sum^{N}_{i=1}\boldsymbol{z}_{i}(g)\boldsymbol{x}^{\prime}_{i}\right)^{-1}\sum^{N}_{i=1}\boldsymbol{z}_{i}(g)y_{i}.
\end{align*}
Define the approximate variance of $\hat{\beta}$ as:
\begin{align*}
\boldsymbol{\mathcal{V}}_{\delta,\mathcal{E}}\left[\boldsymbol{z}\right] & =e^{\prime}_{1}\expec{\delta}{\sum^{N}_{i=1}\boldsymbol{z}_{i}(g)\boldsymbol{x}^{\prime}_{i}}^{-1}\var{\delta,\mathcal{E}}{\sum^{N}_{i=1}\boldsymbol{z}_{i}(g)\varepsilon_{i}}\expec{\delta}{\sum^{N}_{i=1}\boldsymbol{z}_{i}(g)\boldsymbol{x}^{\prime}_{i}}^{-1\prime}e_{1},
\end{align*}
where $e_{1}=(1,0)^{\prime}$, with $\boldsymbol{\mathcal{V}}_{\delta,\mathcal{E}}\left[\boldsymbol{z}\right]=\infty$
whenever $\expec{\delta}{\sum^{N}_{i=1}\boldsymbol{z}_{i}(g)\boldsymbol{x}^{\prime}_{i}}$
is singular. Also, as before, define
\begin{align*}
\mathcal{V}_{\delta,\mathcal{E}}[z] & =\frac{\var{\delta,\mathcal{E}}{\sum^{N}_{i=1}z_{i}(g)\varepsilon_{i}}}{\expec{\delta}{\sum^{N}_{i=1}z_{i}(g)x_{i}}^{2}}.\hspace{0pt}
\end{align*}
An initial result shows that the problem of searching over recentered
IV vectors $\boldsymbol{z}$ is equivalent to a constrained version
of the original problem of searching over scalar $z$ that requires
$z$ to be orthogonal to the “nuisance” exposure $h$:
\begin{lem}
\label{lemma:directFX}Fix $\delta$ and $\mathcal{E}$. Let $\mathcal{Z}^{(2)}_{\delta}$
be the set of recentered $\boldsymbol{z}$ and suppose $\operatorname{tr}(\var{\delta}h)>0$
for $h=(h_{i})^{N}_{i=1}$. Then:
\begin{align*}
\left\{ \boldsymbol{\mathcal{V}}_{\delta,\mathcal{E}}\left[\boldsymbol{z}\right]:z\in\mathcal{Z}^{(2)}_{\delta}\right\} & =\left\{ \V{\delta,\mathcal{E}}z:z\in\mathcal{Z}_{\delta},\expec{\delta}{z^{\prime}h}=0\right\} .
\end{align*}
That is, the set of attainable $\boldsymbol{\mathcal{V}}_{\delta,\mathcal{E}}\left[\boldsymbol{z}\right]$
equals the set of attainable $\V{\delta,\mathcal{E}}z$ under the
constraint $\expec{\delta}{z^{\prime}h}=0$. Moreover, this equivalence
is constructive in both directions: any $\boldsymbol{z}\in\mathcal{Z}^{(2)}_{\delta}$
can be transformed into an instrument vector whose first component
is some scalar $z\in\mathcal{Z}_{\delta}$ with $\expec{\delta}{z'h}=0$,
without changing the value of the objective, while any such scalar
$z$ can be implemented by the vector instrument $\boldsymbol{z}_{i}(g)=\bigl(z_{i}(g),\ h_{i}-\expec{\delta}{h_{i}}\bigr)$.
\end{lem}
Imposing the orthogonality constraint yields the analog of Theorem
\ref{prop:OptimalIV}:
\begin{thm}
\label{prop:OptimalIV_direct}Fix $\delta$ and let
\begin{align}
c_{\delta} & \in\arg\min_{c\in\mathbb{R}}\operatorname{tr}\text{\ensuremath{\left(\left((W-cU)\Sigma_{\delta}(W-cU)^{\prime}\right)^{p/(p+1)}\right)},}\label{eq:c_delta}\\
\tilde{x}^{\perp}_{\delta} & =(W-c_{\delta}U)(g-\expec{\delta}g),\qquad S^{\perp}_{\delta}=(W-c_{\delta}U)\Sigma_{\delta}(W-c_{\delta}U)^{\prime}.\nonumber
\end{align}
Suppose $S^{\perp}_{\delta}\neq0$ and $\operatorname{tr}(U\Sigma_{\delta}\hspace{0pt} U')>0$.
Then, for any $p\in[1,\infty]$ and $\sigma>0$:
(a) The worst-case approximate variance with the optimal instrument
is:
\[
\min\limits_{\bm{z}\in\mathcal{Z}^{(2)}_{\delta}}\max_{\mathcal{E}\in\mathcal{F}_{p}(\sigma)}\mathcal{\boldsymbol{\mathcal{V}}}_{\delta,\mathcal{E}}\left[\bm{z}\right]=N^{1/p}\sigma^{2}\operatorname{tr}\left(\left(S^{\perp}_{\delta}\right)^{p/(p+1)}\right)^{-(p+1)/p}.
\]
(b) Assume that either $p>1$ or $p=1$ and $\operatorname{rank}{\!}\left((W-c_{\delta}U)\Sigma^{1/2}_{\delta}\right)=\min\left\{ N,\operatorname{rank}(\Sigma_{\delta})\right\} $.
Then the instrument $\boldsymbol{z}^{\ast}_{i}(g)=(z^{*}_{\delta}(g)_{i},h_{i}-\expec{\delta}{h_{i}})$
minimizes $\max_{\mathcal{E}\in\mathcal{F}_{p}(\sigma)}\boldsymbol{\mathcal{V}}_{\delta,\mathcal{E}}\left[\boldsymbol{z}\right]$
over $\mathcal{Z}^{(2)}_{\delta}$, where:
\begin{align*}
z^{*}_{\delta}(g) & =\left((S^{\perp}_{\delta})^{\dag}\right)^{1/(p+1)}\tilde{x}^{\perp}_{\delta}.
\end{align*}
\end{thm}
The key difference from the optimal instrument in the single-exposure
case is that the matrix $W$ is replaced with an appropriate residualization
$W-c_{\delta}U$ where $c_{\delta}$ (which plays the role of a Lagrange
multiplier on the $\expec{\delta}{z^{\prime}h}=0$ constraint) can
be solved for by the scalar optimization (\ref{eq:c_delta}). The
additional technical assumption for $p=1$ is imposed because the
objective (\ref{eq:c_delta}) can be non-differentiable, which may
lead to a different analytical form of $z^{*}_{\delta}(g)$. To avoid
this issue, we require that the residualized exposure matrix retains
as much assignment-induced variation as its dimensions allow. Given
orthogonality between the optimal IV for $x_{i}$ and the exposure
$h_{i}$, the IV for $h_{i}$ is simply $h_{i}-\expec{\delta}{h_{i}}$.
Analogs of Propositions \ref{prop:OptimalDesign}--\ref{prop:DesignProperties}
also follow:
\begin{prop}
\label{prop:OptimalDesign_direct}In the Theorem \ref{prop:OptimalIV_direct}
setting, suppose there is $\delta\in\mathcal{D}$ such that $\operatorname{tr}(U\Sigma_{\delta}\hspace{0pt} U')>0$
and $\min_{c\in\mathbb{R}}\operatorname{tr}\left(((W-cU)\Sigma_{\delta}(W-cU)^{\prime})^{p/(p+1)}\right)>0$.
Then
\[
\delta^{*},\boldsymbol{z}^{*}\in\arg\min_{\delta\in\mathcal{D},\boldsymbol{z}\in\mathcal{Z}^{(2)}_{\delta}}\max_{\mathcal{E}\in\mathcal{F}_{p}(\sigma)}\mathcal{\boldsymbol{\mathcal{V}}}_{\delta,\mathcal{E}}\left[\boldsymbol{z}\right]
\]
is solved by
\begin{align*}
\delta^{*} & \in\arg\max_{\delta\in\mathcal{D}}\operatorname{tr}\left((S^{\perp}_{\delta})^{p/(p+1)}\right).
\end{align*}
The optimal instrument is $\boldsymbol{z}^{\ast}=\boldsymbol{z}^{\ast}_{\delta^{\ast}}$
from Theorem \ref{prop:OptimalIV_direct}, as long as, for $p=1$,
\textup{$\operatorname{rank}{\!}\left((W-c_{\delta}U)\Sigma^{1/2}_{\delta}\right)=\min\left\{ N,\operatorname{rank}(\Sigma_{\delta})\right\} $
}.
\end{prop}
\begin{prop}
\label{prop:DesignProperties_direct}There is a $\delta^{*}$ from
Proposition \ref{prop:OptimalDesign_direct} with $\expec{\delta^{\ast}}{g_{k}}=0.5$
for all $k$. Further, if $W$ and $U$ are jointly block diagonal
under the same row-column partition, there is such a $\delta^{*}$
with independent subvectors of $g$ across blocks.
\end{prop}
\setcounter{prop}{2} \setcounter{thm}{1}
We note that, unlike with Proposition \ref{prop:DesignProperties},
the optimal design cannot generally be solved block-by-block in the
block-diagonal case because $c_{\delta}$ is common across all blocks.
We illustrate these results in a simple but nontrivial special case:
when $x_{i}$ is the leave-out average of shocks in the cluster to
which $i$ belongs and $h_{i}=g_{i}$ captures the direct effect.
Appendix \ref{appx:Leave-Out-Averages} shows that, at least when
$N_{j}$ is even, it is optimal to mix between cluster-specific $Bernoulli(0.5)$
shocks and complete randomization of the shocks within each cluster,
independently across clusters. In the special case where all clusters
are of the same size, $N_{j}=n$, one can obtain a solution more similar
to the one from Section \ref{subsec:Special-Cases}: mix between cluster-specific
shocks with probability $\psi^{\perp}=(1-r^{\perp})/(1+(n-1)r^{\perp})$
for $r^{\perp}=(n-1)^{-2p/(2p-1)}$ and independent $Bernoulli(0.5)$
randomization with probability $1-\psi^{\perp}$. The optimal spillover
instrument is proportional to $(\bar{g}_{j(i)}-0.5)+(n-1)^{1/(2p-1)}(\bar{g}_{j(i)}-g_{i})$.
Relative to the case without direct effects, there is generally more
randomization within clusters ($\psi^{\perp}\le\psi$, with strict
inequality for $n>2$ and $p>1$), which helps isolate spillover effects
from the direct effects. The exception is when $p=1$: with both homogeneous
and heterogeneous $N_{j}$, the design and instrument are then the
same as in the original problem, because the objective already creates
enough variation to identify direct effects.
\subsection{Other Extensions}
We derive three further extensions to the baseline results. First,
Appendix \ref{appx:predetermined_covs} shows how the results extend
when the researcher includes an intercept and possibly other predetermined
covariates $r_{i}$ in estimation. We formalize this idea by enriching
the class of error distributions to $\varepsilon_{i}=r^{\prime}_{i}\gamma+\tilde{\varepsilon}_{i}$
where $\tilde{\varepsilon}=\left(\tilde{\varepsilon}_{i}\right)$
satisfies $\expec{\mathcal{E}}{\tilde{\varepsilon}\tilde{\varepsilon}'}\in\mathcal{F}_{p}(\sigma)$
while $\gamma$ is unrestricted. The results in Theorem \ref{prop:OptimalIV}
and Propositions \ref{prop:OptimalDesign}--\ref{prop:DesignProperties}
continue to hold, with the exposure matrix $W$ residualized columnwise
on the covariates.
Second, Appendix \ref{appx:nonlinear} discusses how Theorem \ref{prop:OptimalIV}
and Proposition \ref{prop:OptimalDesign} extend to the general formula
treatment setting where the $x_{i}$ are nonlinear (but still known)
$i$-specific formulas combining $W$ and $g$. This allows, for example,
$x_{i}$ to be an indicator for $i$ having at least one treated friend
in their social network or the more elaborate formulas considered
in \citet{BH1}. The proofs to Theorem \ref{prop:OptimalIV} and Proposition
\ref{prop:OptimalDesign} turn out to extend verbatim after defining
the general $\tilde{x}_{\delta}=x-\expec{\delta}x$ and $S_{\delta}=\var{\delta}x$,
making them a robust characterization of optimal design with formula
treatments. Proposition \ref{prop:DesignProperties}, however, does
not extend cleanly: it is no longer without loss of generality to
consider designs with $0.5$ marginals. The following discussion of
feasible design computation and inference is also reliant on linearity
of the spillover treatment formula.
Finally, Appendix \ref{appx:budgets} considers settings where the
researcher faces a constraint of a marginal treatment probability
$q\neq0.5$ (common across units), e.g. due to a fixed budget or government
target. Theorem \ref{prop:OptimalIV} applies to any design and is
unaffected by the constraints on feasible designs. We therefore show
that Proposition \ref{prop:OptimalDesign} extends naturally with
such constraints, as well as the second claim of Proposition \ref{prop:DesignProperties}
and the results on computation, below.
\section{Feasible Designs and Inference \label{sec:computation_clt}}
We now introduce a relaxed version of the optimal design problem.
This relaxation admits a closed-form solution in special cases and
a computationally efficient approximation in general. Moreover, this
approximation can yield $\sqrt{K}$-consistent recentered IV estimators
with $p<\infty$, as well as a simple asymptotic inference procedure.
\subsection{The Relaxed Problem\label{subsec:Computation}}
Solving the Proposition \ref{prop:OptimalDesign} problem is computationally
intractable for even moderate $K$. While its objective depends on
the experimental design $\delta$ only through the $K\times K$ covariance
matrix $\var{\delta}g$, the class of all $K\times K$ covariance
matrices achievable with binary treatments, $\mathcal{C}_{K}=\left\{ \var{\delta}g\colon\delta\in\mathcal{D}\right\} $,
is high-dimensional. The computation time required to find an optimal
$\var{\delta}g$, by searching over convex combinations of the $2^{K}$
possible treatment vectors in $\left\{ 0,1\right\} ^{K}$, is generally
exponential in $K$.
To make progress, we follow \citet{goemans_improved_1995} and \citet{thiyageswaran_optimal_2026}
in considering a relaxed problem with a polynomial-time solution in
$K$. Specifically, we replace the optimization over $\mathcal{C}_{K}$
with optimization over $\mathcal{Q}_{K}=\left\{ \Sigma\in\mathbb{S}^{K}:\Sigma\succeq0,\Sigma_{kk}=1/4\text{ for all }k\right\} $,
where $\mathbb{S}^{K}$ is the set of $K\times K$ real and symmetric
matrices. Thus $\mathcal{Q}_{K}$ is the set of $K\times K$ covariance
matrices with $1/4$ on the diagonal, containing $\mathcal{C}_{K}$
when we restrict attention to designs with $\expec{\delta}{g_{k}}=0.5$
(without loss of generality, by Proposition \ref{prop:DesignProperties}).
The relaxed Proposition \ref{prop:OptimalDesign} problem is then:
\begin{equation}
\max_{\Sigma\in\mathcal{Q}_{K}}\operatorname{tr}\left((W\Sigma W')^{\frac{p}{p+1}}\right).\label{eq:relax2}
\end{equation}
Searching over $\mathcal{Q}_{K}$ is feasible for moderate $K$ using
an interior point solver for conic optimization problems. Two additional
results further improve computational efficiency without any additional
approximation cost. First, if $W$ can be split into independent blocks
$W_{(1)},\dots,W_{(B)}$ with $K(b)$ shocks in each, we can use Proposition
\ref{prop:DesignProperties} to solve the optimal design problem separately
by block. Second, as the following result shows, we can characterize
the solution explicitly in terms of a lower-dimensional convex optimization
problem over the $K-1$-dimensional simplex $\Delta_{K}=\left\{ \ell\in\mathbb{R}^{K}\colon\ell_{k}\ge0,\sum^{K}_{k=1}\ell_{k}=1\right\} $
rather than $\mathcal{Q}_{K}$; this problem is amenable to gradient-based
optimization approaches that scale more efficiently with $K$.
\begin{prop}
\label{prop:Optimize_l}For $p<\infty$, let
\[
\Sigma^{\ast}=\frac{1}{4}D^{-1/2}_{\ell^{\ast}}\frac{M^{p}_{\ell^{\ast}}}{\operatorname{tr}\left(M^{p}_{\ell^{\ast}}\right)}D^{-1/2}_{\ell^{\ast}},
\]
where, for $\ell>0$,
\[
D_{\ell}=\operatorname{diag}\left(\ell_{1},\dots,\ell_{K}\right),\qquad M_{\ell}=D^{-1/2}_{\ell}W'WD^{-1/2}_{\ell},
\]
and
\begin{equation}
\ell^{\ast}\in\arg\min_{\ell\in\Delta_{K}}\operatorname{tr}\left(M^{p}_{\ell}\right),\label{eq:lstar}
\end{equation}
with $\ell^{*}>0$. Then $\Sigma^{\ast}$ solves the relaxed problem
(\ref{eq:relax2}).
\emph{This characterization is especially helpful in two special cases
where it yields closed-form solutions of the relaxed problem:}
\end{prop}
\begin{cor}
\label{cor:p1}For $p=1$, the relaxed problem (\ref{eq:relax2})
is solved by:
\[
\Sigma^{\ast}_{kl}=\frac{1}{4}\frac{w_{\cdot k}'w_{\cdot l}}{\left\Vert w_{\cdot k}\right\Vert _{2}\left\Vert w_{\cdot l}\right\Vert _{2}},\qquad k,l\in\left\{ 1,\dots,K\right\} ,
\]
where $w_{\cdot k}$ is the $k$th column of $W$. That is, in the
relaxed problem's solution, the correlation of assignments of any
two shocks equals the cosine similarity of the vectors of exposures
to those shocks.
\end{cor}
\begin{cor}
\label{cor:sym_diagonal}For $p<\infty$, if all diagonal elements
of $\left(W'W\right)^{p}$ are equal, the relaxed problem (\ref{eq:relax2})
is solved by:
\[
\Sigma^{\ast}=\frac{1}{4}\frac{\left(W'W\right)^{p}}{\operatorname{tr}\left(\left(W'W\right)^{p}\right)/K}.
\]
The Theorem \ref{prop:OptimalIV} optimal IV corresponding to $\Sigma^{\ast}$
is $z^{\ast}=a\cdot\left(W'\right)^{\dag}\tilde{g}$ which depends
on $p$ only through $a=\left(4\operatorname{tr}\left(\left(W'W\right)^{p}\right)/K\right)^{1/(p+1)}$.
Moreover, when $W$ consists of multiple independent blocks with the
diagonal elements of $\left(W'W\right)^{p}$ equal within blocks,
these expressions apply block by block, with $a$ possibly varying
across blocks.
\end{cor}
Corollary \ref{cor:p1} captures a simple intuition: shocks to two
intervention units should be correlated to the extent that they are
connected, in the sense that the same outcome units are exposed to
them. This turns out to be the optimal solution to the relaxed problem
when $p=1$. Corollary \ref{cor:sym_diagonal} extends this result
to arbitrary $p$, as long as each block of $W$ is sufficiently symmetric---a
restrictive condition that nevertheless includes all of the special
cases in Section \ref{subsec:Special-Cases}. The corresponding solution
is particularly intuitive when $W$ is nonnegative and $p$ is an
integer. In that case, the matrix $(W'W)^{p}$ captures $p$th-degree
connections between intervention units $k$. When $p=2$, for example,
the solution suggests correlating two shocks when they are either
connected directly or connected to the same third shock. The worst-case
loss from such longer-range shock correlations is lower with higher
$p$, as error correlations are more restricted. As $p\to\infty$,
this solution correlates all shocks in the same connected component
of the bipartite network defined by $W$.
Corollary \ref{cor:sym_diagonal} further characterizes the IV that
is optimal if the relaxed problem's solution can be implemented (an
issue we discuss below). This optimal IV satisfies $W'z^{\ast}\propto\tilde{g}$,
meaning that partial whitening entails the inverse operation to the
exposure mapping, $\tilde{x}=W\tilde{g}$. For instance, when $x$
involves some averaging of the (correlated) shocks, $z^{\ast}$ performs
their deconvolution.\footnote{Interestingly, the optimal IV is invariant to $p$, up to the proportionality
constant $a$. With larger $p$, the shocks are more strongly correlated;
thus, the same IV that fully whitens the shocks when $p=1$ does only
partial whitening for larger $p$. We also note that when $W$ consists
of independent blocks the constant $a$ varies across blocks; thus,
$p$ affects the reweighting of blocks but not the optimal IV within
blocks.}
All of these results generalize immediately to settings with an intercept
and possibly other predetermined covariates in estimation, by first
residualizing $W$ on those covariates column by column. For instance,
for $p=1$ and when estimation includes an intercept as the only covariate,
Corollary \ref{cor:p1} applied to the demeaned $w_{\cdot k}$ implies
that shock correlations under the relaxed problem's solution are equal
to the correlations of the $N\times1$ exposures to those shocks (rather
than cosine similarities without the intercept). This solution positively
correlates shocks with similar exposure while introducing slight negative
correlations between shocks with non-overlapping exposure.\footnote{That is, if $w_{\cdot k},w_{\cdot l}\ge0$ and $w_{\cdot k}'w_{\cdot l}=0$
(no exposure overlap), $\Sigma^{\ast}_{kl}$ is proportional to the
sample correlation of $w_{\cdot k}$ and $w_{\cdot l}$, which is
negative.}
We can also generalize the relaxed problem to multiple exposures,
where it remains computationally efficient. By the same arguments
as in the Proposition \ref{prop:OptimalDesign_direct} proof, the
relaxed problem is:
\begin{equation}
\min\limits_{c\in\mathbb{R}}\max_{\Sigma\in\mathcal{Q}_{K}}\operatorname{tr}\Big(\left((W-cU)\Sigma(W-cU)')\right)^{\frac{p}{p+1}}\Big).\label{eq:relaxp}
\end{equation}
Here the inner problem has the same structure as before with a single
exposure, with $W$ replaced by $W-cU$; the computationally efficient
form of Proposition \ref{prop:Optimize_l} and its closed-form special
cases thus apply.\footnote{The Proposition \ref{prop:Optimize_l} problem becomes $\min_{c\in\mathbb{R}}\min_{\ell\in\interior{\Delta_{K}}}\operatorname{tr}\left\{ \left(D^{-1/2}_{\ell}\left(W-cU\right)'\left(W-cU\right)D^{-1/2}_{\ell}\right)^{p}\right\} $.}
Additionally, when $W$ and $U$ are jointly block-diagonal, the inner
problem can be split block-by-block for any $c$. The outer problem
is then a scalar optimization, amenable to a golden section search
over a compact subset of $\mathbb{R}$.\footnote{The minimizer is guaranteed to exist in a compact subset of $\mathbb{R}$
as long as $\max\limits_{\Sigma\in\mathcal{Q_{K}}}\operatorname{tr}\Big(\left((W-cU)\Sigma(W-cU)')\right)^{\frac{p}{p+1}}\Big)\ge\operatorname{tr}\Big(\left((W-cU)I/4(W-cU)')\right)^{\frac{p}{p+1}}\Big)\to\infty$
when $\left|c\right|\to\infty$, which requires only that $U\neq0$.}
\subsection{Feasible Implementation\label{subsec:Feasible-Implementation}}
Two challenges remain once a solution $\Sigma^{\ast}$ to the relaxed
problems in Equation (\ref{eq:relax2}) or (\ref{eq:relaxp}) is found.
First, this $\Sigma^{\ast}$ may not correspond to an element of $\mathcal{C}_{K}$:
i.e., it may not be implementable via binary shock vectors. Second,
even if it is theoretically implementable, there is no known computationally
feasible procedure for sampling binary vectors from a generic covariance
matrix. These problems can be overcome in simpler special cases: the
designs in Section \ref{subsec:Special-Cases} are indeed implementations
of the relaxed problem solution.
Outside special cases, we address both challenges by using a simple
Gaussian rounding procedure to sample binary assignments with the
variance matrix approximating $\Sigma^{\ast}$, following \citet{goemans_improved_1995}.
Specifically, we draw $\xi\sim\mathcal{N}(0,\Sigma^{\ast})$ and set
$g^{GR}_{k}=\mathbf{1}\left[\xi_{k}\ge0\right]$.\footnote{This procedure simplifies further in the special case of Corollary
\ref{cor:p1}: draw $\eta\sim\mathcal{N}(0,I_{N})$ and set $g^{GR}_{k}=\mathbf{1}\left[(W'\eta)_{k}\ge0\right]$.
Indeed, we can write $g^{GR}_{k}=\mathbf{1}\left[\xi_{k}\ge0\right]$ for
$\xi_{k}=\sum_{i}w_{ik}\eta_{i}/2\left\Vert w_{\cdot k}\right\Vert _{2}$,
where $(\xi_{1},\dots,\xi_{K})\sim\mathcal{N}(0,\Sigma^{\ast})$ for
$\Sigma^{\ast}$ from Corollary \ref{cor:p1}.} This procedure distorts correlations in a known way: $\corr{}{g^{GR}_{k},g^{GR}_{l}}=\arcsin\left[\corr{}{\xi_{k},\xi_{l}}\right]/\frac{\pi}{2}\ne\corr{}{\xi_{k},\xi_{l}}$;
see \citet{goemans_improved_1995} and the illustration in Appendix
Figure \ref{fig:arcsin}. Given this Gaussian rounding design, we
find the optimal IV based on $S_{GR}=W\var{}{g^{GR}}W^{\prime}$ for
\begin{equation}
\Sigma^{GR}=\var{}{g^{GR}}=\frac{1}{2\pi}\arcsin\left[4\Sigma^{\ast}\right],\label{eq:V}
\end{equation}
with $\arcsin\left[\cdot\right]$ applied entrywise. Algorithm \ref{alg:full}
summarizes the entire procedure.
\begin{algorithm}
\caption{\label{alg:full}Feasible Optimal Experimental Design and IV}
\begin{enumerate}
\item[0.] \textbf{Input:}
\begin{itemize}
\item $N\times K$ matrix of the indirect exposure of interest, $W$
\item \emph{Optional:} $N\times K$ matrix of the other included exposure,
$U$; e.g., $U=I$ for direct exposure
\item \emph{Optional:} $N\times L$ matrix of predetermined covariates,
$R$
\item Tuning parameter $p\in[1,\infty)$. Set lower $p$ for higher robustness
to non-spherical errors
\end{itemize}
\item[1.] If $R\ne\emptyset$, replace $W$ and $U$ with their columnwise
residuals after projecting on $R$
\item[2.] Split $(W,U)$ jointly into independent blocks of rows and columns,
$(W_{(b)},U_{(b)})$, if any
\begin{itemize}
\item \emph{Optional: }to ease computation, at the cost of some approximation
error, also split approximately independent blocks (e.g., if $R=\mathbf{1}$,
the zero elements in the original $(W,U)$ for rows and columns from
different blocks will be replaced in step 1 with $-1/N\approx0$)
\end{itemize}
\item[3.] Fix $c=0$. Solve for the relaxed-optimal shock covariance matrix
$\Sigma^{\ast}$ block-by-block:
\begin{itemize}
\item If $p=1$, use the Corollary \ref{cor:p1} closed-form solution
\item If all diagonal elements of $(W'_{(b)}W_{(b)})^{p}$ are equal, use
the Corollary \ref{cor:sym_diagonal} closed-form solution
\item Else, use the Proposition \ref{prop:Optimize_l} numerical optimization
\end{itemize}
\item[4.] If $U\ne\emptyset$, choose $c_{\delta}$ that solves the optimization
problem in (\ref{eq:relaxp}) by repeating Step 3 for $c\ne0$ and
$W-cU$ replacing $W$ (on a grid or using gradient-based methods)
\item[5.] Randomize shocks using the Gaussian rounding design $\delta$ from
Section \ref{subsec:Feasible-Implementation}. \textbf{Output $g$}
\item[6.] Set the recentered instrument as in Theorem \ref{prop:OptimalIV}
with $\var{\delta}g=\Sigma^{GR}$ from (\ref{eq:V}). \textbf{Output
$z^{\ast}_{\delta}(g)$}
\begin{itemize}
\item If $U\ne\emptyset$, add the simple IV from Theorem \ref{prop:OptimalIV_direct}
for the other exposure. \textbf{Output $\boldsymbol{z}^{\ast}_{\delta}(g)$}
\end{itemize}
\end{enumerate}
\end{algorithm}
Beyond simplicity in implementation, this procedure offers two advantages.
First, Appendix \ref{appx:approx_bound} shows the cost of using this
relaxed and rounded version of the problem is bounded: the worst-case
approximate variance of the resulting estimator is at most $\pi/2$
times that of the original problem.\footnote{We find much smaller approximation errors in applications: for all
designs in Section \ref{sec:Applications}, the worst-case approximate
variance exceeds the variance of the original problem by no more than
6\%. This is computed by comparing the worst-case variance arising
from the Gaussian rounded design to that of the unrounded solution
to the relaxed problem (\ref{eq:relax2}), which is a lower bound
for the variance of the optimal solution to the original problem.} Second, as we next show, the fact that our shocks are a simple transformation
of a correlated Gaussian vector is helpful for establishing a central
limit theorem for the recentered IV estimator and asymptotically valid
inference, even when all errors are mutually correlated.
\subsection{\label{subsec:clt-inference}Asymptotic Normality and Inference}
We now consider the asymptotic behavior of the estimator $\hat{\beta}=\hat{\beta}\left[z^{\ast}\right]$
corresponding to the Gaussian rounding implementation of the Section
\ref{subsec:Computation} relaxed-optimal design, with the Theorem
\ref{prop:OptimalIV} optimal IV $z^{\ast}$, as in Algorithm \ref{alg:full}.
We show that choosing a small integer $p$ can yield sparsity of the
implied shock covariance matrix, which in turn can yield a central
limit theorem for $\hat{\beta}$ and a simple inference procedure.
Let $d(\Sigma)=\max_{k}\sum^{K}_{l=1}\mathbf{1}\left\{ \Sigma_{kl}\neq0\right\} $
measure the sparsity of a covariance matrix $\Sigma$ by the maximum
number of non-zero covariances across all rows. We consider a sequence
of data-generating processes indexed by $K$. For $K\to\infty$ (which
under Assumption \ref{ass:asymptotic} will require $N\to\infty$
as well) we assume:
\begin{assumption}
\label{ass:asymptotic}There are constants $\underline{\lambda}>0$,
\textup{$\bar{\lambda}\in(0,\infty)$,} $\bar{d}<\infty$, $h^{*}\in(0,\infty)$,
and $v^{*}\in(0,\infty)$ such that:
\begin{enumerate}
\item[(a)] \textup{$d(\Sigma^{\ast})\le\bar{d}$ and }$\lambda_{\max}(W'W)\le\bar{\lambda}$
\textup{;}
\item[(b)] \textup{}$\lambda^{+}_{\min}(S_{GR})>\underline{\lambda}$ with $\lambda^{+}_{\min}(\cdot)$
denoting the smallest positive eigenvalue;
\item[(c)] \textup{}$h_{K}=K^{-1}\operatorname{tr}(S^{p/(p+1)}_{GR})\xrightarrow{p} h^{*}$;
\item[(d)] \textup{}$v_{K}=K^{-1}\varepsilon'S^{(p-1)/(p+1)}_{GR}\varepsilon\stackrel{p}{\rightarrow}v^{*}$;
\item[(e)] \textup{}$\frac{1}{K^{3/2}}\sum^{K}_{k=1}|b_{k}|^{3}=o_{p}(1)$ for
$b=W'(S^{\dag}_{GR})^{1/(p+1)}\varepsilon$,
\end{enumerate}
where statements (a) and (b) hold with probability $1-o(1)$ and all
probability statements are with respect to the joint distribution
of $(W,\varepsilon,\xi)$.
\end{assumption}
The first part of Assumption \ref{ass:asymptotic}(a) is the key sparsity
condition imposed on the solution to the relaxed problem in Equation
(\ref{eq:relax2}). While the sparsity of $\Sigma^{*}$ is formally
an asymptotic condition, the researcher can heuristically check whether
$d(\Sigma^{\ast})$ is small compared to $K$; at the cost of more
complex primitive conditions, we can accommodate $d(\Sigma^{\ast})$
growing slowly in $K$.\footnote{Theorem \ref{thm:main} can be extended to cases where $\Sigma^{*}$
is approximately sparse, in that it contains many small entries but
has row and column sums that are bounded or grow very slowly. This
extension is particularly useful when predetermined covariates are
included, since then the effective $W$ and corresponding $\Sigma^{*}$
will not generally be sparse. In such cases we could allow the operator
norm of $\Sigma^{*}$ to grow slowly by a smoothing argument and the
second order Poincaré inequality of \citet{chatterjee2009fluctuations},
as applied in \citet{lei2018asymptotics}. We thank Lihua Lei for
pointing this out.} Since $\frac{1}{N}\left\Vert x\right\Vert ^{2}_{2}=\frac{1}{N}\left\Vert Wg\right\Vert ^{2}_{2}\le\lambda_{\max}(W'W)\cdot\frac{1}{K}\left\Vert g\right\Vert ^{2}\cdot\frac{K}{N}$,
the second part precludes $x$ from diverging provided $K\asymp N$.
Appendix \ref{app:Theorem1} provides sufficient conditions for both
parts of this assumption. It shows that $d(\Sigma^{\ast})$ is upper-bounded
by $\left(d(W'W)\right)^{p}$ under the assumptions of Proposition
\ref{prop:Optimize_l}. With $p=\infty$, no sparsity of $\Sigma^{\ast}$
is guaranteed within any connected component of the graph, and so
designs that are optimal for $p=\infty$ may not yield consistent
$\hat{\beta}$ when $W$ does not consist of many independent blocks
and when errors are also strongly dependent. But for smaller integer
$p$, the optimal design only correlates units that are connected
by paths in $W'W$ of lengths $p$ or less. Thus, for finite integer
$p$, both parts of Assumption \ref{ass:asymptotic}(a) are guaranteed
when maximum row and column degree of $W$ as well as the maximum
absolute value of its elements are all bounded. This set of conditions
in turn requires $K\asymp N$ in non-degenerate cases, as also shown
in Appendix \ref{app:Theorem1}.
In addition to sparsity, Assumption \ref{ass:asymptotic} imposes
regularity conditions to ensure the joint distribution of spillover
treatments $x_{i}$ and errors $\varepsilon_{i}$ is well-behaved.
Assumption \ref{ass:asymptotic}(b) precludes near-collinearity of
exposures to ensure that the whitening matrix remains well-conditioned.
Yet, this condition allows $S_{GR}$ to be degenerate, which happens,
e.g., when some observations have the same exposure to all shocks.
Noting that $h_{K}=\expec{\delta}{\frac{1}{K}x'z^{\ast}\mid W}$,
Assumption \ref{ass:asymptotic}(c) implies that $W$ and the experimental
design are such that there is sufficient variation in the spillover
treatments $x_{i}$ so that the first stage is not degenerate asymptotically.
Noting further that $v_{K}=\frac{1}{K}\var{\delta}{\varepsilon'z^{\ast}\mid W,\varepsilon}$,
Assumption \ref{ass:asymptotic}(d) ensures that the finite-sample
variance of $\varepsilon'z^{\ast}$ conditional on $W$ and $\varepsilon$
(scaled appropriately) converges to a deterministic constant. This
is a law of large numbers for the quadratic form $v_{K}$ and it imposes
regularity on the joint distribution of $W$ and $\varepsilon$. For
example, conditional on a sequence of $W$, if $\varepsilon$ is drawn
from a mixture of two deterministic sequences that lead to different
limits for $v_{K}$, then this assumption is violated; when dependence
in $\varepsilon$ is sufficiently weak or local, then this assumption
is satisfied. Assumption \ref{ass:asymptotic}(e) is an anti-concentration
condition for the implicit weights of the numerator $\frac{1}{K}\varepsilon'z^{\ast}$;
combined with the sparsity condition, it ensures that the numerator
can be asymptotically approximated by sums of independent components.
Under these conditions, $\widehat{\beta}$ is $\sqrt{K}$-consistent:
\begin{thm}
\label{thm:main} Under Assumption~ \ref{ass:asymptotic} and the
Gaussian rounding design based on the relaxed problem's solution from
Proposition \ref{prop:Optimize_l}, the Theorem \ref{prop:OptimalIV}
optimal recentered IV estimator $\hat{\beta}$ satisfies
\[
\sqrt{K}(\widehat{\beta}-\beta)\Rightarrow\mathcal{N}\!\left(0,\frac{v^{*}}{h^{*2}}\right).
\]
\end{thm}
\noindent Asymptotically valid inference is straightforward from
this result using the normal approximation and a simple plug-in estimator
for the standard error, $\sqrt{\hat{v}/(Kh^{2}_{K})}$. Here $S_{GR}$
and therefore $h_{K}$ are known. Moreover, the plug-in estimator
for $v^{\ast}$ is straightforward to compute: $\hat{v}=\frac{1}{K}\hat{\varepsilon}\,'S^{(p-1)/(p+1)}_{GR}\hat{\varepsilon}$
for $\hat{\varepsilon}=y-\hat{\beta}x$. The convergence in probability
of $\hat{v}$ to $v^{*}$ follows from consistency of $\widehat{\beta}$
and Assumption \ref{ass:asymptotic} (see Appendix \ref{app:Theorem1}).
\section{Applications\label{sec:Applications}}
\subsection{Bipartite Experiment: Cai et al. (2015)}
\paragraph{Setup.}
Our first application is the bipartite experiment of \citet{cai_social_2015}
who study how social networks affect the adoption of weather insurance
by randomly assigning rice farmers in villages in China to simple
vs. intensive information sessions over two rounds. They collect information
on demographics, friendship, and insurance take-up for each farmer.
We calibrate a simulation to a simplified version of the social network
effect specification in Table 2, Column 2 of their paper. For the
sample of farmers assigned to either the simple or intensive information
session in the second round, it examines how insurance take-up varies
with the fraction of friends that are treated in the first round.
Specifically, our outcome units $i$ are those randomized in the second
round ($N=1,274$),\footnote{More precisely, outcome units are those used to estimate spillover
effects on insurance take-up (from “Type 1” villages) rather
than price effects.} and our intervention units $k$ are those randomized in the first
round with at least one friend in the second round ($K=995$). As
in \citet{cai_social_2015}, we define $x_{i}=\sum_{k}w_{ik}g_{k}$
where $w_{ik}$ equals an indicator that $i$ and $k$ are friends
divided by the total number of friends $i$ has (including friends
outside the first-round sample). Thus, $x_{i}$ measures the percentage
of the second round unit's friends who were treated in the first round.
Each $i$ and $k$ is assigned to one of 44 administrative villages,
which are administrative units that each combine multiple local communities;
throughout this analysis we drop 10 of 6,251 friendships that are
across administrative villages so $W$ has block-diagonal structure
that eases computational burden. The top row of Figure \ref{fig:Applications}(a)
shows the exposure matrix $W$ for one example village, along with
the $W'W$ matrix that captures similarity in exposure weights across
intervention units.
We estimate the specification $y_{i}=\alpha+\beta x_{i}+\upsilon_{i}$
using recentered IV without whitening. Our simple specification has
a smaller estimated coefficient of $\beta^{\ast}=0.1408$ compared
to the specification in \citet{cai_social_2015}, which uses OLS as
an estimator and includes an additional set of controls and village
fixed effects, rather than recentering.
For $S=1000$ simulations, we generate the outcome as $y_{i}=\beta^{*}x_{i}+\varepsilon_{i}$,
where $g_{k}$ varies by the experimental design, $W$ is fixed, and
we draw $\varepsilon=(\varepsilon_{i})$ from six different data-generating
processes (DGPs), i.e., joint distributions that vary the structure
of the $N\times N$ matrix $\expec{\mathcal{E}}{\varepsilon\varepsilon^{\prime}}$.\footnote{While the original outcomes in \citet{cai_social_2015} are binary,
our generated outcomes are continuous.} For the first five DGPs, the distribution of the errors is calibrated
so that the standard error of the unwhitened recentered IV estimator
of $\beta$ under a 50\% Bernoulli RCT matches the standard error
of the estimated coefficient in the data. The error distributions
for the simulation are as follows:
\begin{itemize}
\item Homoskedastic: $\varepsilon_{i}\stackrel{iid}{\sim}\mathcal{N}(0,\nu^{2})$,
with $\nu=0.56$;
\item Heteroskedastic: $\varepsilon_{i}\sim\mathcal{N}(0,\nu^{2}(1+\gamma d_{i}))$,
independent across units $i$, where $d_{i}$ is the number of $i$'s
first-round friends, $\nu=0.34$ and $\gamma=0.82$, which implies
50\% of the variance is driven by the heteroskedastic component;
\item Degree-in-Mean: $\varepsilon_{i}=\gamma d_{i}+\nu\zeta_{i}$, where
$\zeta_{i}\stackrel{iid}{\sim}\mathcal{N}(0,1)$, $\gamma=0.29$,
and $\nu=0.31$, with $\gamma$ calibrated such that 50\% of the cross-sectional
variation of $\varepsilon_{i}$ arises from the mean;
\item Village-Correlated: $\varepsilon_{i}=\nu\cdot(\zeta_{i}+u_{v(i)})$,
where $\nu=0.33$, $v(i)$ is the village of individual $i$, and
$u_{v}\stackrel{iid}{\sim}\mathcal{N}(0,1)$, $\zeta_{i}\stackrel{iid}{\sim}\mathcal{N}(0,1)$;
\item Network-Correlated: $\varepsilon=\nu\zeta+\gamma W_{\text{raw}}\eta$,
where $\zeta_{i}\stackrel{iid}{\sim}\mathcal{N}(0,1)$, $\eta_{k}\stackrel{iid}{\sim}\mathcal{N}(0,1)$,
and $W_{\text{raw}}$ is the un-normalized friendship matrix, so $W_{\text{raw},i,k}=1$
if $i$ and $k$ are friends and is otherwise 0. The parameters $\gamma=0.29$
and $\nu=0.18$ are calibrated so that 50\% of the variance comes
from the network component;
\item Estimated Residuals: $\varepsilon_{i}=\hat{\upsilon}_{i}$, the residuals
from the regression on the real data.
\end{itemize}
\begin{figure}
\caption{Spillover Structure and Optimal Designs in Applications\label{fig:Applications}}
\begin{centering}
(a) An example village in \citet{cai_social_2015}, $K=34$, $N=44$
\par\end{centering}
\begin{centering}
\includegraphics[width=0.95\textwidth]{cai_designs}
\par\end{centering}
\begin{centering}
\medskip{}
(b) \citet{Miguel2004}, $K=N=49$
\par\end{centering}
\begin{centering}
\includegraphics[width=0.95\textwidth]{mg_designs}\smallskip{}
\par\end{centering}
{\small\emph{Notes}}{\small : Panel (a) shows an example of one village
(Xinlian) in the \citet{cai_social_2015} data. The first row shows
the exposure matrix $W$ and the matrix $W'W$ with the $kl$ element
measuring the extent to which $k$ and $l$ are friends of the same
farmers $i$. The second row shows the correlation matrices of the
Algorithm \ref{alg:full} feasible optimal designs corresponding to
$p=2$ and $p=1$, along with the analogous solutions for $p=\infty$.
Panel (b) reports the same objects for the full set of schools in
\citet{Miguel2004}, except $W$ and $W'W$ are normalized to be between
0 and 1 by dividing each matrix element-wise by its maximum value.}{\small\par}
\end{figure}
\paragraph{Results.}
Following Algorithm \ref{alg:full}, we solve for the feasible optimal
design and recentered IV estimator of the indirect effect for several
choices of $p$. The bottom row of Figure \ref{fig:Applications}(a)
plots a correlation heatmap for these designs. As our theory suggests,
$p=\infty$ involves perfectly correlating all treatments within clusters
of connected farmers (with independent randomization across villages),
while lower $p$ implies more diffuse correlation structures.\footnote{Appendix Figure \ref{fig:frontier} shows how the optimal designs
in both applications trade off isotropy, as measured by $\operatorname{tr}(S_{\delta^{*}})/\lambda_{\max}(S_{\delta^{*}})$,
and signal $\operatorname{tr}(S_{\delta^{*}})$. As expected, the frontier moves
down and to the right as $p\rightarrow\infty$. Notably, the baseline
independent Bernoulli design (indicated by the red RCT dot) has higher
isotropy than any optimal design but also much lower signal.}
\begin{table}[tp]
\caption{RMSE and Effective Sample Size Gains for the \citet{cai_social_2015}
Application\label{tab:cai}}
\begin{centering}
\begin{tabular}{l*{6}{c}}
\toprule
& \shortstack{Homo- \\ skedastic} & \shortstack{Hetero- \\ skedastic} & \shortstack{Deg.-in- \\ Mean} & \shortstack{Village- \\ Corr.} & \shortstack{Network- \\ Corr.} & \shortstack{Estimated \\ Residuals} \\
Design and Estimator: & (1) & (2) & (3) & (4) & (5) & (6) \\
\midrule
\shortstack[l]{RCT Benchmark \\ \strut} & \shortstack{0.150 \\ {[+0\%]}} & \shortstack{0.139 \\ {[+0\%]}} & \shortstack{0.142 \\ {[+0\%]}} & \shortstack{0.138 \\ {[+0\%]}} & \shortstack{0.141 \\ {[+0\%]}} & \shortstack{0.129 \\ {[+0\%]}} \\
\shortstack[l]{Optimal, $p=1$ \\ \strut} & \shortstack{0.131 \\ {[+32\%]}} & \shortstack{0.123 \\ {[+27\%]}} & \shortstack{0.121 \\ {[+37\%]}} & \shortstack{0.115 \\ {[+44\%]}} & \shortstack{0.107 \\ {[+73\%]}} & \shortstack{0.111 \\ {[+35\%]}} \\
\shortstack[l]{Optimal, $p=2$ \\ \strut} & \shortstack{0.115 \\ {[+70\%]}} & \shortstack{0.111 \\ {[+57\%]}} & \shortstack{0.119 \\ {[+41\%]}} & \shortstack{0.105 \\ {[+72\%]}} & \shortstack{0.106 \\ {[+75\%]}} & \shortstack{0.103 \\ {[+55\%]}} \\
\shortstack[l]{Optimal, $p=\infty$ \\ \strut} & \shortstack{0.095 \\ {[+151\%]}} & \shortstack{0.099 \\ {[+96\%]}} & \shortstack{0.133 \\ {[+13\%]}} & \shortstack{0.117 \\ {[+39\%]}} & \shortstack{0.107 \\ {[+73\%]}} & \shortstack{0.099 \\ {[+71\%]}} \\
\bottomrule
\end{tabular}
\smallskip{}
\par\end{centering}
{\small\emph{Notes}}{\small : For the \citet{cai_social_2015} application,
this table reports the root mean-squared error of different design-estimator
pairs (rows) under different data-generating processes for the errors
(columns). The first row is the baseline of independent Bernoulli
assignment and unwhitened recentered IV. The other rows are the Algorithm
\ref{alg:full} feasible optimal design and optimal IV for different
finite values of $p$ (and the analogous solutions for $p=\infty$),
with $R=\mathbf{1}$. The text describes the six error DGPs. Increases in
effective sample sizes, given in brackets, are computed as the squared
ratio of inverse RMSE under the focal scenario relative to the RCT
baseline with the same error distributions, minus one.}{\small\par}
\end{table}
Table \ref{tab:cai} reports root mean-squared error (RMSE) for the
optimal recentered IV estimator from Theorem \ref{prop:OptimalIV};
RMSE represents the standard deviation of the recentered IV estimator
since the bias is negligible. All estimators include an intercept.\footnote{Correspondingly, Algorithm \ref{alg:full} is applied with $R=\mathbf{1}$.
We use the optional part of step 2 to block-diagonalize the residualized
$W$. This allows us to solve for the optimal design village by village,
with independent randomization across villages.} The columns correspond to the error structures described above, while
the rows correspond to the optimal designs, along with a benchmark
of a simple Bernoulli experiment with a 50\% treatment probability
and a recentered IV estimator without whitening. In brackets, we also
include the gain in effective sample size, defined as the squared
RMSE ratio relative to the benchmark.
All optimal designs beat a standard RCT, which is not optimized to
estimate spillover effects. Under homoskedastic errors, the design
with $p=\infty$ (i.e., perfectly correlating treatments within blocks)
works best. But under more complex and realistic error structures,
the advantage of extreme shock correlations deteriorates: the $p=2$
design is best for the DGP in columns 3--5 with network-dependent
or correlated errors. This choice of $p$ also shows robust power
gains across all columns, equivalent to increasing sample size of
41--75\% in all columns.
To isolate the gain from the optimal design vs. the optimal estimator,
Appendix Table \ref{tab:cai-extra} reports RMSE of the recentered
IV estimator without whitening for the $p$-optimal designs. On one
hand, the $p=\infty$ design involves no whitening yet achieves substantial
improvements across most error distributions---suggesting a key role
of design. Moreover, for simple (e.g., homoskedastic) errors whitening
only reduces precision gains. On the other hand, column 5 shows that
with more challenging error distributions extra power gains due to
the $p=1$ and $p=2$ designs are only realized if the instrument
is also whitened. Appendix Table \ref{tab:Coverage-Cai} also includes
an analysis of coverage for confidence intervals built using the normal
approximation from Theorem \ref{thm:main}. Since the network of \citet{cai_social_2015}
is made up of many independent villages, all designs and simulations
have close to nominal coverage, even when shocks are perfectly correlated
within villages under $p=\infty$.
\subsection{Direct and Indirect Effects: Miguel and Kremer (2004)}
\paragraph{Setup.}
Our second application is in the setting of \citet{Miguel2004} who
estimate the direct and spillover effects of a randomized school-level
deworming treatment on helminth infection rates among children in
Kenya. The intervention shocks $g_{i}$ are at the school level; slightly
simplifying the original analysis, we define the outcome $y_{i}$
as the average infection rate in school $i$. Thus, the outcome and
intervention units are the same set of $N=K=49$ schools (specifically,
Groups 1 and 2 of the experiment for which health outcomes are observed).
Following \citet{Miguel2004}, the spillover treatment $x_{i}$ is
the number (rather than the share) of students in nearby schools that
were treated, rescaled by 1,000. This is formalized by the $N\times N$
adjacency matrix $W$ with entry $w_{ik}=E_{k}/1000$ if school $i$
is within 3km of school $k\ne i$, where $E_{k}$ is the number of
pupils eligible for the treatment in school $k$, and 0 otherwise
(including $w_{ii}=0$ for all $i$). The top row of Figure \ref{fig:Applications}(b)
shows the exposure matrix $W$ and the matrix of similarity in exposure
weights $W'W$.
We calibrate our simulation to a simplified version of column (1)
in Table VII of \citet{Miguel2004}. We estimate a school-level regression
$y_{i}=\alpha+\beta x_{i}+\tau g_{i}+\upsilon_{i}$ using recentered
IV without whitening. This yields coefficients $\beta^{\ast}=-0.226$
and $\tau^{\ast}=-0.258$.\footnote{Our specification differs from \citet{Miguel2004} in several ways.
We use a linear model rather than probit, we recenter rather than
including additional controls, and we drop a third treatment variable
based on exposure to schools between 3 and 6 km away. Despite these
simplifications, our estimates are very similar to the corresponding
average marginal effects in the original paper, of -0.26 for the indirect
effect and -0.25 for the direct effect.}
For $S=1,000$ simulations, we generate the outcome as
\begin{equation}
y_{i}=\beta^{*}x_{i}+\tau^{*}g_{i}+\varepsilon_{i},\label{eq:MK_specification}
\end{equation}
where again the distribution of $g$ depends on the experimental design,
$W$ is fixed, and we vary $\expec{\mathcal{E}}{\varepsilon\varepsilon'}$
according to the same six DGPs described above, with a few minor changes.
We replace the village-correlated DGP with a “cluster-correlated”
one by creating five clusters of schools using spectral clustering
on $W$ and use random cluster shocks. For the network-correlated
DGP, we set $W_{\text{raw}}=W$. We calibrate the simulations like
in the first application: to match the estimated standard error of
the unwhitened recentered IV estimator and targeting the same shares
of variation explained by different error components.\footnote{The corresponding parameters are{\small{} as follows. Homoskedastic:
$\nu=0.17$; heteroskedastic: $\gamma=0.4$, $\nu=0.12$; degree-in-mean:
$\gamma=0.056$, $\nu=0.085$, cluster-correlated: $\nu=0.10$; network-correlated:
$\gamma=0.12,\nu=0.05$.}}
\begin{table}[tp]
\caption{RMSE and Effective Sample Size Gains for the \citet{Miguel2004} Application\label{tab:RMSE-MK}}
\begin{centering}
(a) Spillover Effect\smallskip{}
\par\end{centering}
\begin{centering}
\begin{tabular}{l*{6}{c}}
\toprule
& \shortstack{Homo- \\ skedastic} & \shortstack{Hetero- \\ skedastic} & \shortstack{Deg.-in- \\ Mean} & \shortstack{Cluster- \\ Corr.} & \shortstack{Network- \\ Corr.} & \shortstack{Estimated \\ Residuals} \\
Design and Estimator: & (1) & (2) & (3) & (4) & (5) & (6) \\
\midrule
\shortstack[l]{RCT Benchmark \\ \strut} & \shortstack{0.101 \\ {[+0\%]}} & \shortstack{0.099 \\ {[+0\%]}} & \shortstack{0.128 \\ {[+0\%]}} & \shortstack{0.106 \\ {[+0\%]}} & \shortstack{0.109 \\ {[+0\%]}} & \shortstack{0.179 \\ {[+0\%]}} \\
\shortstack[l]{Optimal, $p=1$ \\ \strut} & \shortstack{0.085 \\ {[+43\%]}} & \shortstack{0.087 \\ {[+30\%]}} & \shortstack{0.087 \\ {[+114\%]}} & \shortstack{0.075 \\ {[+98\%]}} & \shortstack{0.064 \\ {[+191\%]}} & \shortstack{0.120 \\ {[+124\%]}} \\
\shortstack[l]{Optimal, $p=2$ \\ \strut} & \shortstack{0.070 \\ {[+110\%]}} & \shortstack{0.075 \\ {[+76\%]}} & \shortstack{0.069 \\ {[+241\%]}} & \shortstack{0.070 \\ {[+128\%]}} & \shortstack{0.064 \\ {[+190\%]}} & \shortstack{0.114 \\ {[+148\%]}} \\
\shortstack[l]{Optimal, $p=\infty$ \\ \strut} & \shortstack{0.052 \\ {[+277\%]}} & \shortstack{0.056 \\ {[+215\%]}} & \shortstack{0.049 \\ {[+570\%]}} & \shortstack{0.067 \\ {[+147\%]}} & \shortstack{0.057 \\ {[+267\%]}} & \shortstack{0.111 \\ {[+160\%]}} \\
\bottomrule
\end{tabular}\bigskip{}
\par\end{centering}
\begin{centering}
(b) Direct Effect\smallskip{}
\par\end{centering}
\begin{centering}
\begin{tabular}{l*{6}{c}}
\toprule
& \shortstack{Homo- \\ skedastic} & \shortstack{Hetero- \\ skedastic} & \shortstack{Deg.-in- \\ Mean} & \shortstack{Cluster- \\ Corr.} & \shortstack{Network- \\ Corr.} & \shortstack{Estimated \\ Residuals} \\
Design and Estimator: & (1) & (2) & (3) & (4) & (5) & (6) \\
\midrule
\shortstack[l]{RCT Benchmark \\ \strut} & \shortstack{0.052 \\ {[+0\%]}} & \shortstack{0.049 \\ {[+0\%]}} & \shortstack{0.040 \\ {[+0\%]}} & \shortstack{0.041 \\ {[+0\%]}} & \shortstack{0.034 \\ {[+0\%]}} & \shortstack{0.064 \\ {[+0\%]}} \\
\shortstack[l]{Optimal, $p=1$ \\ \strut} & \shortstack{0.054 \\ {[-7\%]}} & \shortstack{0.054 \\ {[-18\%]}} & \shortstack{0.042 \\ {[-7\%]}} & \shortstack{0.046 \\ {[-22\%]}} & \shortstack{0.036 \\ {[-14\%]}} & \shortstack{0.067 \\ {[-10\%]}} \\
\shortstack[l]{Optimal, $p=2$ \\ \strut} & \shortstack{0.054 \\ {[-7\%]}} & \shortstack{0.051 \\ {[-7\%]}} & \shortstack{0.039 \\ {[+6\%]}} & \shortstack{0.045 \\ {[-18\%]}} & \shortstack{0.040 \\ {[-27\%]}} & \shortstack{0.068 \\ {[-12\%]}} \\
\shortstack[l]{Optimal, $p=\infty$ \\ \strut} & \shortstack{0.051 \\ {[+6\%]}} & \shortstack{0.048 \\ {[+5\%]}} & \shortstack{0.040 \\ {[+1\%]}} & \shortstack{0.046 \\ {[-20\%]}} & \shortstack{0.038 \\ {[-23\%]}} & \shortstack{0.100 \\ {[-59\%]}} \\
\bottomrule
\end{tabular}\smallskip{}
\par\end{centering}
{\small\emph{Notes}}{\small : For the \citet{Miguel2004} application,
Panel (a) reports the root mean-squared error of different design-estimator
pairs (rows) for the spillover effect under different data-generating
processes for the errors (columns). Panel (b) reports the same for
the direct effect. The first row is the baseline of independent Bernoulli
assignment and unwhitened recentered IV. The other rows are the Algorithm
\ref{alg:full} feasible optimal design and optimal IV for different
finite values of $p$ (and the analogous solutions for $p=\infty$),
with $R=\mathbf{1}$. The text describes the six error DGPs. Increases in
effective sample sizes, given in brackets, are computed as the squared
ratio of inverse RMSE under the focal scenario relative to the RCT
baseline with the same error distributions, minus one.}{\small\par}
\end{table}
\paragraph{Results.}
We follow Algorithm \ref{alg:full} to solve for the feasible design
optimized for the \emph{indirect} effect in a linear model (\ref{eq:MK_specification})
that also includes a direct effect. As in the first application, we
include the intercept for all estimators.\footnote{We again use the algorithm with $R=\mathbf{1}$. Here, however, we do not
need the optional part of step 2.}
The bottom row of Figure \ref{fig:Applications}(b) presents the correlation
heatmap for the optimal designs used in the simulations. Now that
the direct effect is in the model, we no longer have a perfectly correlated
design when $p=\infty$. As $p$ decreases, the optimal design from
Proposition \ref{prop:OptimalDesign_direct} weakens both positive
and negative correlations of the assignments to protect against adversarially
correlated (again, positively or negatively) error distributions.
The RMSE results for the indirect and direct effects are in Table
\ref{tab:RMSE-MK}, including in brackets the effective sample size
gain compared to the baseline RCT design with an unwhitened recentered
IV estimator. For estimating the spillover effect, Panel (a) shows
that all optimal designs beat the simple RCT under all error structures.
Selecting $p=\infty$ gives the best performance in all columns. Error
correlations are not adversarial enough for $p=2$ to dominate the
other designs, as it did in the previous application,\footnote{A possible reason is that here $x_{i}$ is the number (rather than
share) of treated neighbors, so there is more scope for increasing
signal by correlating the assignments---which is what the $p=\infty$
design targets.} but it remains a good choice, especially with more complex error
structures. Overall, the power gains are even more substantial than
in the first application: with the $p=2$ choice, they range from
+76\% in column 2 to +241\% in column 3, while with $p=\infty$ they
reach +570\% in column 3.\footnote{To isolate the performance impact of the optimal design vs. the optimal
estimator, Appendix Table \ref{tab:MK-extra} reports RMSE of the
recentered IV estimator without whitening for the $p$-optimal designs.} Yet, Appendix Table \ref{tab:Coverage-MK} provides a reason to prefer
a conservative choice of $p$: under $p=\infty$, the resulting covariance
matrix of the shocks is too dense to meet the conditions of Theorem
\ref{thm:main}. For the simulations with correlated errors, confidence
intervals built using the normal approximation under the $p=\infty$
design severely undercover the target parameter, while the $p=1$
design has close to nominal coverage under all DGPs.
Since our proposed designs are optimized for the indirect effect,
they can lead to some deterioration of performance of the direct effect
estimator, as reported in Panel (b) of Table \ref{tab:RMSE-MK}. For
$p=\infty$ this deterioration can be substantial, equivalent to a
59\% reduction in sample size in column 6. However, for $p=1$ and
$p=2$ the deterioration is mild, of at most 22\% and 27\%, respectively.
\section{Conclusion\label{sec:Conclusion}}
We have shown how researchers can design experiments and choose recentered
estimators to more precisely estimate spillover effects, under weak
restrictions on the distribution of unobservables. The core insight
of this approach is that spillover signal can be increased, by correlating
the assignments of connected units, without overly concentrating variation
in directions where the unobservables may be adversarially clustered.
Balancing this tradeoff between signal and isotropy yields a tractable
minimax problem with intuitive special cases and practical approaches
to computation. In two semi-synthetic experiments, based on high-profile
applications, we find large decreases in spillover effect standard
errors relative to conventional randomization and estimation that
are possible without large deterioration in the standard errors for
the direct effects in the application with such effects.
Several open questions remain. First, while we have focused on linear
spillover treatments, and derived theoretical extensions to nonlinear
spillover treatment constructions, computation of optimal designs
in nonlinear settings remains a challenge. Second, though we have
characterized the optimal design and estimator for spillovers in the
presence of direct effects, we have not presented the more general
problem of minimizing a weighted average of their variances. Third,
we have not shown how the optimal design problem changes when a researcher
has access to some prior information on the distribution of unobservables,
such as from a first-wave pilot (\citealp{viviano_experimental_2025}).
Finally, we have not studied how the minimax approach changes when
a researcher is interested in estimating a more complicated parametric
model of spillovers, such as from a structural economic model (\citealp{borusyak_estimating_2025}).
We hope to address some of these extensions in future drafts.
\begin{singlespace}
\addcontentsline{toc}{section}{References}
\bibliographystyle{aer}
\bibliography{OptimalRCT}
\end{singlespace}