# Inference for Dependent Data with Learned Clusters

Jianfei<EMAIL_ADDRESS>[ Chicago Booth  Northeastern
University
Department of Economics
360 Huntington Ave., 301 Lake Hall
Boston, MA 02115, USA

Damian<EMAIL_ADDRESS>[ University of Zürich
University of Zürich
Department of Economics
Schönberggasse 1, 8001 Zürich, Switzerland

Christian<EMAIL_ADDRESS>[ Chicago Booth  The
University of Chicago
Booth School of Business
5807 S. Woodlawn, Chicago, IL 60637, USA

Lucciano Villacortalabel=e4]<EMAIL_ADDRESS>[ Central Bank of Chile
Central Bank of Chile
Agustinas 1180, Santiago, Chile

###### Abstract

This paper presents and analyzes an approach to cluster-based inference for
dependent data. The primary setting considered here is with spatially indexed
data in which the dependence structure of observed random variables is
characterized by a known, observed dissimilarity measure over spatial indices.
Observations are partitioned into clusters with the use of an unsupervised
clustering algorithm applied to the dissimilarity measure. Once the partition
into clusters is learned, a cluster-based inference procedure is applied to a
statistical hypothesis testing procedure. The procedure proposed in the paper
allows the number of clusters to depend on the data, which gives researchers a
principled method for choosing an appropriate clustering level. The paper
gives conditions under which the proposed procedure asymptotically attains
correct size. A simulation study shows that the proposed procedure attains
near nominal size in finite samples in a variety of statistical testing
problems with dependent data.

Key Words: Unsupervised Learning, Cluster-based Inference, HAR inference.

JEL Codes: C1.

t1First version: October 2018. This version is of . Christian Hansen would
like to thank the National Science Foundation as well as The University of
Chicago Booth School of Business for financial support of this research.
Damian Kozbur would like to thank The University of Zürich for financial
support of this research.

##  Introduction

Conducting accurate statistical inference with data featuring serial or cross-
sectional dependence requires carefully accounting for the underlying
dependence structure. A variety of methods are available for researchers
analyzing dependent data. An important class of inferential methods is the
class of cluster-based methods. Cluster-based methods work with a partition of
observations into clusters. Inference proceeds by treating as negligible any
covariance between observations that fall in different clusters, followed by
performing appropriate tests for a statistical hypothesis of interest.
Cluster-based methods are known to deliver asymptotically valid inference, in
the sense of controlling size of tests and coverage of interval estimates, in
a variety of dependent-data settings; see, for example, Ibragimov and Müller,
(2010), Bester et al., (2011), Canay et al., (2017), and Hansen and Lee,
(2019).

Though the term can apply more generally, in empirical economics settings
“clustering” is most commonly used in the special case of an $n\times T$ panel
data setting, in which individual level observations $i=1,...,n$ are clustered
together across the $T$ dimension. This setting has been analyzed extensively,
and one comprehensive reference is the textbook Wooldridge, (2010). In the
case of a panel data setting, there is natural, “known” way to partition data
into clusters, which is justified from an assumption on the data, namely that
observations $i$ and $j$ are independent whenever $i\neq j$.

However, there is not an immediate “natural” set of partitions available in
many settings. Rather, researchers often use data where observations can be
thought of as being located in some space and where a measure of distance or
dissimilarity between observations is available. For example, researchers
routinely use data where observations are indexed geographically over an
irregular geographic region and spatial correlation is a concern. As another
example, researchers may be working with data where observations are equipped
with a notion of economic distance. Examples of such applications include
Conley and Dupor, (2003), which constructs a notion of economic distance based
on input-output tables, and Conley and Topa, (2002), which constructs a notion
of economic distance based on demographic composition of zip codes. In such
settings, researchers may wish to use cluster-based inference as a robust
method to account for dependence in the data but desire some guidance about
exactly which observations to cluster together and how many clusters to use.

In this paper, we consider settings in which natural partitions may not exist
but a distance or dissimilarity measure that is informative about the
underlying correlation structure is available. We assume that dependence
between observable quantities is suitably small when the distance between them
is suitably large. A main goal is then to verify that cluster-based methods
remain asymptotically valid when clusters are not pre-specified but are
constructed in a data-dependent fashion using the underlying notion of
dissimilarity and ‘unsupervised learning’ methods. There are many such tools
available for constructing groups from measures of dissimilarity including
$k$-medoids, $k$-means, and hierarchical clustering; see Hastie et al., (2009)
for a general review.

At a high level, the data-driven cluster-based inference approach we propose
for testing a hypothesis $H_{0}$ proceeds in three main steps. In the first
step, unsupervised learning methods are applied using the observed
dissimilarity between observations to produce a collection of partitions of
the observations into clusters,
$\mathscr{C}=\\{\mathcal{C}^{(k)}\\}_{k=2}^{k_{\max}}$ for some integer
$k_{\max}$. Each element $\mathcal{C}^{(k)}$ provides a candidate structure
for use with a cluster-based method and may be based on a different number of
clusters. Cluster-based methods for testing $H_{0}$ will involve decision
thresholds for determining whether to reject $H_{0}$. For example, we may
decide to reject $H_{0}$ when an empirical p-value is smaller than some
threshold. In the second step, we select a partition $\widehat{\mathcal{C}}$
from $\mathscr{C}$ and a rejection rule for $H_{0}$ chosen to maximize exact
weighted average power against a collection of alternatives subject to
maintaining exact size control within completely specified data generating
processes (DGPs) for the observed data. These DGPs may be data dependent; for
example, they may involve parameters whose values are estimated in the data.
While $\widehat{\mathcal{C}}$ and the rejection rule are selected based on
size and power calculations in pre-specified DGPs, the selection is done in
such a way that asymptotic size control is formally maintained regardless of
whether these DGPs are correctly specified. In the final step, inference in
the actual data is performed using the rejection rule and
$\widehat{\mathcal{C}}$ determined in the second step. As $\mathscr{C}$ will
typically contain candidate partitions with different numbers of groups, the
procedure outlined above encompasses a data-driven method for choosing the
appropriate level at which to cluster the data as well as a data-driven way to
determine to which cluster each observation should be assigned.

The main theoretical results in this paper provide conditions under which
inference based on the general procedure outlined in either Ibragimov and
Müller, (2010) or Canay et al., (2017) remains valid when data-dependent
clusters produced by applying $k$-medoids are used in lieu of relying on a
fixed number of pre-specified clusters. The regularity conditions involve
moment and mixing restrictions, weak homogeneity assumptions on second moments
of regressors and unobservables across groups, and restrictions on group
boundaries. These moment and mixing conditions are implied by routine
assumptions necessary for use of central limit approximations and the required
homogeneity is less restrictive than covariance stationarity. The assumptions
allow for relatively general heterogeneous dependent processes and are akin to
conditions routinely imposed in other heteroskedasticity and autocorrelation
robust (HAR) inference approaches; see, e.g. Andrews, (1991) and Kelejian and
Prucha, (2007).

In the course of providing the main inferential results, we also develop novel
insights about the finite-sample behavior of the $k$-medoids clustering
algorithm that are important for cluster-based inference. Specifically, we
verify that the partitions produced by $k$-medoids satisfy suitably defined
balance and small boundary conditions without requiring a notion of
consistency to a “true” partition. Establishing these properties for the
groups produced by $k$-medoids is important for verifying their use as an
input into cluster-based inference as existing results for such procedures in
the setting where clusters are pre-specified rely on conditions that impose
that clusters have small boundaries and have sizes that do not diverge
relative to each other. We also note that these results on the behavior of
$k$-medoids clustering are, to our knowledge, new to the unsupervised learning
literature and may be of some independent interest.

As cluster-based methods provide a general approach to performing robust
inference in the presence of dependent and heteroskedastic data, they belong
to the general class of HAR inference procedures; see, e.g., Lazarus et al.,
(2018) for a recent review centered on time series analysis and Conley, (1999)
and Kelejian and Prucha, (2007) for seminal references in the spatial context.
Our results cover scenarios where the number of clusters used for inference is
bounded as in Bester et al., (2011), Ibragimov and Müller, (2010), and Canay
et al., (2017) and thus are in the spirit of the “fixed-b” approach - e.g.
Kiefer et al., (2000), Kiefer and Vogelsang (2002, 2005), Bester et al.,
(2016) \- and the projection approach with fixed number of projections - e.g.
Phillips, (2005), Müller, (2007), Sun, (2013) \- both of which are well-known
to lead to improved size control relative to more traditional approaches under
weak dependence. Within the general context of HAR inference, the choice of
cluster structure is analogous to the choice of smoothing parameters in other
HAR procedures, and our work thus complements the broad literature on data-
driven tuning parameter choice in HAR inference. Our proposal in this paper is
closely related to Sun and Kim, (2015), Lazarus et al., (2018), Lazarus et
al., (2021), and especially Müller and Watson, (2020). Sun and Kim, (2015)
proposes a method to select the number of projections to use for HAR inference
in a spatially dependent context by minimizing asymptotic mean-square error
(MSE) of the variance estimator where terms that enter the expression for the
asymptotic MSE are estimated using a parametric model. Lazarus et al., (2021)
provide the size-power frontier based on the “fixed-b” approach in a time
series context, and Lazarus et al., (2018) explicitly considers tuning
parameter choice by trading off size and power along this frontier. Müller and
Watson, (2020) considers a novel standard error estimator for use with
spatially dependent data based on spatial principal components and proposes a
method to choose both the number of components to use in the method and the
critical value to use in hypothesis testing based on minimizing confidence
interval length subject to exactly controlling size within benchmark
parametric models. Their procedure exactly controls size in the baseline model
and can be constructed to asymptotically control size under general
conditions. Relative to these approaches, we focus on data-dependent
construction of clusters for use with cluster-based inference methods. Similar
to Müller and Watson, (2020), we choose both a partition and rejection rule
for a hypothesis test to maximize weighted average power subject to exact size
control within a parametric model in such a way that asymptotic size control
is maintained. In a related paper, Bai et al., (2020) propose a new standard
error estimator for panel data models robust to serial and spatial dependence
when clusters are unknown, and dependence is weak. In their approach the time
series dependence is controlled using a Newey-West HAC estimator, while the
spatial dependence is controlled by assuming sparsity in the cross-section
where only a small subgroup of $\left(i,j\right)$ pairs have spatial
dependence different from zero. The $\left(i,j\right)$ pairs of this subgroup
are unknown and are estimated by applying the thresholding approach in Bickel
and Levina, (2008). A similar procedure is proposed by Cai, (2021).

The remainder of this paper is organized as follows. Section 2 presents a
high-level overview of our proposed inferential procedure. We illustrate a
concrete use of the method by applying it to the study of effect of insurgent
attacks on voter turnout from Condra et al., (2018) in Section 3. Section 4
then illustrates the performance of our proposal in simulated data. Formal
results are provided in Sections 5. Finally, the Appendix contains a detailed
discussion of the implementation we use in the empirical and simulation
examples. All proofs are contained in additional supplemental material.

##  Methodology: Inference with Unsupervised Cluster Learning

Consider data given by $\mathscr{D}=\\{\zeta_{i}\\}_{i\in\mathsf{X}}$. Here,
$\zeta_{i}$ are observable random variables or vectors and $\mathsf{X}$ is a
(spatial) indexing set of cardinality $n$. This paper assumes that
$\mathsf{X}$ is equipped with a known dissimilarity measure $\mathrm{d}$,
which is an $n\times n$ array of nonnegative real dissimilarities. When added
emphasis is helpful, $\mathsf{X}$ is written $(\mathsf{X},\mathrm{d})$. The
data $\mathscr{D}$ is distributed according to an unknown joint probability
distribution $\mathscr{D}\sim{\mathrm{P}}_{0}$. The object $\mathsf{X}$ will
be the main object used to characterize any dependence in the data
$\mathscr{D}$ over $i$. Finally, ${\mathrm{P}}_{0}$ is an element of a larger
known class of distributions ${\mathrm{P}}_{0}\in\textbf{P}$.

Consider testing a null hypothesis, $H_{0}$, at level $\alpha\in(0,1)$. As a
concrete example, we will focus on the problem of testing a hypothesis about a
coefficient in a linear regression model (i.e. testing $H_{0}:\theta_{0}=0$
where $\theta_{0}$ is a parameter in a linear regression). In the case of this
example, $\mathbf{P}$ consists of all possible data generating processes
(DGPs) for the regression data, which typically include DGPs in which
observations are correlated with each other. Failure to account for dependence
in $\mathscr{D}$ across $i\in\mathsf{X}$ may lead to substantial size
distortion when testing $H_{0}$. In the following, we outline an approach for
obtaining valid inference within this setting using cluster-based inferential
procedures with clusters generated by an unsupervised clustering algorithm.

Let $\mathcal{C}=\\{\mathsf{C}_{1},...,\mathsf{C}_{k}\\}$ be a partition of
$\mathsf{X}$ of cardinality $k\geqslant 2$. The elements
$\mathsf{C}_{1},...,\mathsf{C}_{k}$ are referred to as clusters. A cluster-
based inferential procedure for testing $H_{0}$ is a (possibly random)
assignment

$\displaystyle\mathtt{Test}:(\mathscr{D},\mathcal{C})\mapsto
T\in\big{\\{}\text{Fail to Reject, Reject}\big{\\}}.$ (2.1)

Here, the decision rule itself is called $\mathtt{Test}$ and will generally
depend on the level $\alpha\in(0,1)$ of the test. The outcome of the test
$\mathtt{Test}$ is referred to as $T$, ie.
$T=\mathtt{Test}(\mathscr{D},\mathcal{C})$. The set containing the pair
$(\mathscr{D},\mathcal{C})$ remains unnamed to avoid additional notation.

We focus on the following three cluster-based inferential procedures for
testing a scalar hypothesis in this paper: the procedure of Ibragimov and
Müller, (2010) (IM), the procedure of Canay et al., (2017) (CRS), and
inference based on the cluster covariance estimator as described in Bester et
al., (2011) (CCE).111Extension of formal results for testing joint hypotheses
using CRS is straightforward. For clarity, consider ‘$t$-statistic’ based
testing of the hypothesis $H_{0}:\theta_{0}=\theta^{\circledast}$ where
$\theta_{0}$ is some scalar parameter of interest associated to the data
generating process ${\mathrm{P}}_{0}$ and $\theta^{\circledast}$ is a
hypothesized possible value for $\theta_{0}$. Let
$\mathsf{C}\subseteq\mathsf{X}$ and $\widehat{\theta}_{\mathsf{C}}$ be an
estimator of $\theta_{0}$ using only data corresponding to observations in
$\mathsf{C}$. Now define $S\in\mathbb{R}^{\mathcal{C}}$ and the $t$-statistic
function $t:\mathbb{R}^{\mathcal{C}}\rightarrow\mathbb{R}$ such that

$\displaystyle S=(S_{\mathsf{C}})_{\mathsf{C}\in\mathcal{C}},\ \
S_{\mathsf{C}}=(n/k)^{1/2}(\widehat{\theta}_{\mathsf{C}}-\theta^{\circledast}),\phantom{\Bigg{|}}$
(2.2) $\displaystyle
t(S)=\frac{k^{-1/2}\sum_{\mathsf{C}\in\mathcal{C}}S_{\mathsf{C}}}{\sqrt{(k-1)^{-1}\sum_{\mathsf{C}\in\mathcal{C}}\left(S_{\mathsf{C}}-k^{-1}\sum_{\mathsf{C}^{\prime}\in\mathcal{C}}S_{\mathsf{C}^{\prime}}\right)^{2}}}.$
(2.3)

For a specified level $a\in(0,1)$ (which may distinct from $\alpha$), there
are the IM, CRS, and CCE tests, denoted by $\mathtt{Test}_{\text{IM}(a)}$,
$\mathtt{Test}_{\text{CRS}(a)}$ , $\mathtt{Test}_{\text{CCE}(a)}$. These tests
are defined by their outcomes given data $\mathscr{D}$ and a partition
$\mathcal{C}$:

$\displaystyle T_{\text{IM}(a),\mathcal{C}}\ \ =\ \text{Reject}\ \ \
\text{if}\ \ \ \left|t(S)\right|>t_{1-a/2,k-1},$ (2.4) $\displaystyle
T_{\text{CRS}(a),\mathcal{C}}=\ \text{Reject}\ \ \text{ if }\ \
\left|t(S)\right|>\mathrm{quantile}_{1-a}(\\{|t(hS)|\\}_{h\in\mathcal{H}_{\mathcal{C}}}),$
(2.5) $\displaystyle T_{\text{CCE}(a),\mathcal{C}}=\ \text{Reject}\ \ \text{
if }\ \
\left|\frac{\widehat{\theta}_{\mathsf{X}}-\theta^{\circledast}}{\widehat{V}_{\text{CCE},\mathcal{C}}^{1/2}}\right|>\sqrt{\frac{k}{k-1}}\times
t_{1-a/2,k-1},$ (2.6)

where $t_{1-a/2,k-1}$ is the $(1-a/2)$-quantile of a $t$-distribution with
$k-1$ degrees of freedom; the set $\\{hS\\}_{h\in\mathcal{H}_{\mathcal{C}}}$
is the orbit of the action of $\\{\pm 1\\}^{\mathcal{C}}$ on $S$, so that for
each $h$, $hS\in\mathbb{R}^{\mathcal{C}}$ has $\mathsf{C}^{\text{th}}$
component $\pm(n/k)^{1/2}(\widehat{\theta}_{\mathsf{C}}-\theta^{\circledast})$
for some sign in $\\{\pm 1\\}$; and $\widehat{V}_{\text{CCE},\mathcal{C}}$ is
the standard cluster covariance matrix estimator. Similarly
$\mathtt{Test}_{\bullet(a)}$ and $T_{\bullet(a),\mathcal{C}}$ is used when the
choice of IM, CRS, or CCE is unspecified. Note that explicit dependence on
$\mathscr{D}$ is supressed in $T_{\bullet(a),\mathcal{C}}$ but that
$T_{\bullet(a),\mathcal{C}}$ may be formally related back to $\mathscr{D}$ by
$T_{\bullet(a),\mathcal{C}}=\mathtt{Test}_{\bullet(a)}(\mathscr{D},\mathcal{C})$.
With prespecified, non-data-dependent clusters, each of the IM, CRS, and CCE
procedures has the favorable property of asymptotic nominal size control under
respective regularity conditions whenever dependence in observations
$\zeta_{i},\zeta_{j}$ with $i,j$ in different clusters is suitably negligible.
The appropriate formal definition of negligible is method-specific.

The second important definition is that of an unsupervised clustering
algorithm, which is an assignment that returns, to every
$\mathsf{X}=(\mathsf{X},\mathrm{d})$, a partition of $\mathsf{X}$ given by the
mapping

$\displaystyle\mathtt{Cluster}:\mathsf{X}\mapsto\mathcal{C}.$ (2.7)

The idea behind using an unsupervised clustering algorithm is that if the
dissimilarity $\mathrm{d}$ appropriately reflects the dependence in
$\zeta_{i}$, then the resulting partition $\mathcal{C}$ may have the desired
property that observations belonging to different clusters exhibit negligible
dependence. In the formal analysis in Section 5, the imposed mixing conditions
imply that dependence between $\zeta_{i}$ and $\zeta_{j}$ vanishes as
$\mathrm{d}(i,j)$ becomes large. Then, if $\mathcal{C}$ places distant
observations (as defined by $\mathrm{d}$) in different clusters, favorable
properties of the test $T$ may be anticipated.

Though there are many commonly used unsupervised clustering algorithms and we
expect most to be usable as methods for forming data-dependent clusters, we
consider only $k$-medoids in this paper for technical reasons discussed in
Section 5.4. We provide a full description of the version of $k$-medoids we
use in Appendix A. Note that by composition of specific $\mathtt{Test}$ and
$\mathtt{Cluster}$ procedures, it is already possible to define an outcome $T$
for $H_{0}$ given $\mathscr{D},\mathsf{X}$ by constructing
$T=\mathtt{Test}(\mathscr{D},\mathtt{Cluster}(\mathsf{X})).$

The final layer to our proposed testing procedure is a simple method for data-
dependent choice of both the cluster-based inferential procedure by
considering a collection of candidate testing and clustering procedures of the
form $\mathtt{Test}$ and $\mathtt{Cluster}$. We propose making this choice on
the basis of simultaneously controlling Type-I and Type-II error rates. Let
$\text{Err}_{\text{Type-I}}(\mathtt{Test},\mathtt{Cluster})$ denote type-I
error for the testing outcome defined by
$\mathtt{Test}(\mathscr{D},\mathtt{Cluster}(\mathsf{X}))$. For a set of
alternatives $\mathbf{P}_{\text{alt}}$, let $\text{Err}_{\text{Type-
II},\mathbf{P}_{\text{alt}}}(\mathtt{Test},\mathtt{Cluster})$ denote a
weighted average type-II error. Because
$\text{Err}_{\text{Type-I}}(\mathtt{Test},\mathtt{Cluster})$ and
$\text{Err}_{\text{Type-
II},\mathbf{P}_{\text{alt}}}(\mathtt{Test},\mathtt{Cluster})$ will typically
not be known, we consider a setting in which estimates
$\widehat{\text{Err}}_{\text{Type-I}}(\mathtt{Test},\mathtt{Cluster})$ and
$\widehat{\text{Err}}_{\text{Type-
II},\mathbf{P}_{\text{alt}}}(\mathtt{Test},\mathtt{Cluster})$ are available.

To finish the final layer, let $\mathscr{T}$ be a collection of pairs of the
form $(\mathtt{Test},\mathtt{Cluster})$. In this paper, interest will
primarily be in $\mathscr{T}$ of the form

$\displaystyle\mathscr{T}_{\bullet(\alpha),k_{\max}}=\Big{\\{}(\mathtt{Test}_{\bullet(a)},k\text{-}\texttt{medoids}):a\in[0,\alpha],k\in\\{2,...,{k_{\mathrm{max}}}\\}\Big{\\}}$
(2.8)

where $\alpha$ is the nominal testing level for $H_{0}$, ${k_{\mathrm{max}}}$
is a researcher-chosen upper bound on the number of clusters, and
$\mathtt{Test}_{\bullet(a)}$ is defined above.222One could also include pairs
involving pre-specified partitions in $\mathscr{T}$. We ignore this
possibility for notational convenience. We then choose
$(\widehat{\mathtt{Test}},\widehat{\mathtt{Cluster}})\in\mathscr{T}$ by
solving

$\displaystyle(\widehat{\mathtt{Test}},\widehat{\mathtt{Cluster}})\in\arg\min\
\widehat{\text{Err}}_{\text{Type-
II},\mathbf{P}_{\text{alt}}}(\mathtt{Test},\mathtt{Cluster})$ (2.9)
$\displaystyle\text{s.t. }{(\mathtt{Test},\mathtt{Cluster})\in\mathscr{T}\ ,\
\widehat{\text{Err}}_{\text{Type-I}}(\mathtt{Test},\mathtt{Cluster})\leqslant\alpha}.$

Set the final testing outcome for $H_{0}$ to

$\displaystyle\widehat{T}=\widehat{\mathtt{Test}}(\mathscr{D},\widehat{\mathtt{Cluster}}(\mathsf{X})).$
(2.10)

In the case of interest, $\mathscr{T}=\mathscr{T}_{\bullet(\alpha),k_{\max}}$,
the “parameter space” in (2.8) depends on two independent parameters $a$ and
$k$. It follows that (2.9) is a two-dimensional optimization problem with a
single constraint. The resulting optimization problem is therefore non-
degenerate. The solution is then determined by two parameters
$\widehat{\alpha}$ and $\widehat{k}$. Furthermore, the testing outcome can be
expressed
$\widehat{T}=\mathtt{Test}_{\bullet(\widehat{\alpha})}(\mathscr{D},\widehat{\mathcal{C}}\hskip
1.13809pt),\ \text{with
}\widehat{\mathcal{C}}=\mathcal{C}^{(\widehat{k})}=\widehat{k}\text{-}\texttt{medoids}(\mathsf{X}).$
When tests are based on $\bullet$ being IM, CRS, or CCE, and when it is
helpful to make the overall level $\alpha$ explicit, write also
$\widehat{T}=\widehat{T}_{\bullet(\alpha)}$ for added emphasis.

$\widehat{k}$ provides a data-dependent answer to how many clusters to use.
Optimization over $a$ allows the parameter entering the data-dependent
decision rule, $\widehat{\alpha}$, to be smaller than the nominal level of the
test, $\alpha$. That is, the data-dependent decision rule may be more
conservative than would be implied by conventional, fixed rules that
asymptotically control size but may fail to do so in finite samples. Finally,
the constraint on $a$ in the definition of $\mathscr{T}$ in (2.8) guarantees
that inference based on $\widehat{T}$ will maintain asymptotic size control
under sensible conditions that do not require the estimator
$\widehat{\text{Err}}_{\text{Type-I}}(\mathtt{Test},\mathtt{Cluster})$ to
agree with the true Type-I error rate in finite-samples or asymptotically.

In practice, estimates
$\widehat{\text{Err}}_{\text{Type-I}}(\mathtt{Test},\mathtt{Cluster})$ and
$\widehat{\text{Err}}_{\text{Type-
II},\mathbf{P}_{\text{alt}}}(\mathtt{Test},\mathtt{Cluster})$ of Type-I and
Type-II error rates are needed. In all results reported in the following
sections, we obtain these estimates by forming preliminary estimates of the
dependence structure for the data based on Gaussian Quasi Maximum Likelihood
Estimation (QMLE) using a simple exponential covariance function. We then use
the estimated dependence structure within a Gaussian model to obtain
$\widehat{\text{Err}}_{\text{Type-I}}(\mathtt{Test},\mathtt{Cluster})$ and
$\widehat{\text{Err}}_{\text{Type-
II},\mathbf{P}_{\text{alt}}}(\mathtt{Test},\mathtt{Cluster})$. Further details
are provided in Appendix A.2. We note that in principle any baseline model
could be used in implementing (2.9) and one could consider uniform size
control over classes of models, $\mathcal{M}$, by replacing the size
constraint in (2.9) with
$\max_{m\in\mathcal{M}}\widehat{\text{Err}}_{\text{Type-I},m}(\mathtt{Test},\mathtt{Cluster})\leqslant\alpha$.
We wish to reemphasize that our theoretical results do not require the
consistency of the estimated dependence structure of ${\mathrm{P}}_{0}$ in
order to control size. As a result, misspecification in the model for
dependence asymptotically leads only to potential loss of power.

(2.9) differs from several recent methods in the literature for choosing data-
dependent tuning parameters for use in conducting inference with dependent
data. Much of the existing literature suggests choosing a single tuning
parameter to optimize a weighted combination of size distortion and power;
see, for instance, Lazarus et al., (2021), Sun and Kim, (2015), and references
therein. Instead, our proposal leverages the fact that most commonly used
inferential procedures for dependent data depend on two parameters - nominal
size and a smoothing parameter - and focuses on maximizing power within
procedures that control size. Our proposal is closely related to Müller and
Watson, (2020) who consider an inference approach for spatially dependent data
that makes use of a tuning parameter and a critical value which are chosen by
minimizing confidence interval length subject to exactly controlling size.
Relative to the existing literature, both Müller and Watson, (2020) and our
approach offer additional flexibility by explicitly considering two choice
variables and make use of criteria, minimizing interval length or maximizing
power subject to maintaining size control, that we believe will be appealing
to many researchers.

For convenience of reference in the following sections, the above described
procedure is stated under Algorithm 1 below. For concreteness, Algorithm 1
references the procedures specialized to IM, CRS, CCE in conjunction with
$k$-medoids (ie for $\mathscr{T}=\mathscr{T}_{\bullet(\alpha),k_{\max}}$). To
simplify notation at the cost of some abuse of notation, write
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathtt{Test}_{\mathrm{\bullet}(a)},k\text{-}\texttt{medoids})=\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\bullet(a),k)$
and $\widehat{\mathrm{Err}}_{\mathrm{Type-
II},\mathbf{P}_{\mathrm{alt}}}(\mathtt{Test}_{\mathrm{\bullet}(a)},k\text{-}\texttt{medoids})=\widehat{\mathrm{Err}}_{\mathrm{Type-
II},\mathbf{P}_{\mathrm{alt}}}(\bullet(a),k)$. The concrete description in
Algorithm 1 is the only version of the procedure described above implemented
in the empirical study and simulation study that follow. The more general
procedure (i.e. for arbitrary $\mathscr{T}$) could be stated analogously.

######  Algorithm 1.

(Inference with Cluster Learning with $k$-medoids and IM, CRS or CCE). Testing
$H_{0}$ at level $0<\alpha<1$.

Data: $\mathscr{D},\mathsf{X}$.

Inputs: $k_{\max}$; $\bullet=$ IM, CRS, or CCE; Estimates
$\widehat{\text{Err}}_{\text{Type-I}}(\bullet(a),k),\widehat{\text{Err}}_{\text{Type-
II},\mathbf{P}_{\text{alt}}}(\bullet(a),k)$

Procedure: Solve (2.9) to obtain $(\widehat{\alpha},\widehat{k})$.

Output: Set
$\widehat{T}_{\bullet(\alpha)}={\mathtt{Test}}_{\bullet(\widehat{\alpha})}(\mathscr{D},\widehat{k}\text{-}\texttt{medoids}({\mathsf{X}}))$.

##  Empirical Application: Insurgent Electoral Violence

In this section, we illustrate the inferential method proposed in Section 2 by
reexamining estimates of the effect of insurgent attacks on voter from Condra
et al., (2018). We use the same data as Condra et al., (2018) which consist of
district-level observations for 205 voting districts in Afghanistan in 2014.
Each observation contains information on direct morning attacks, voter
turnout, and district level control variables for two separate election
rounds.

###  3.1 Model and Spatial Dependence

We focus on column (4) of Table 2 in Condra et al., (2018) which reports the
estimated effect of nearby morning attacks on voter turnout obtained from
instrumental variables (IV) estimation of the linear model

$\displaystyle Y_{de}$
$\displaystyle=\alpha_{0}+\theta_{0}Attacks_{de}+W_{de}^{\prime}\gamma_{0}+U_{de}$
(3.1) $\displaystyle Attacks_{de}$
$\displaystyle=\mu_{0}+\pi_{0}Z_{de}+W_{de}^{\prime}\xi_{0}+V_{de}$ (3.2)

where (3.1) is the structural equation with parameter of interest $\theta_{0}$
and (3.2) is the first stage representation. Here, $Y_{de}$ is the turnout in
district $d$ and election round $e$, and $Attacks_{de}$ is the number of
morning attacks in district $d$ and election round $e$. $U_{de}$ and $V_{de}$
are unobservables where $U_{de}$ and $Attacks_{de}$ are potentially
correlated. Control variables $W_{de}$ include election round dummy variables,
voting hour wind conditions, population, measures of precipitation and ambient
temperature, and the average of predawn and morning wind conditions during the
pre-election period. The excluded instrument, $Z_{de}$, is a measure of early
morning wind conditions.

The analysis in Condra et al., (2018) reports inference for the parameter
$\theta_{0}$ based on CCE with clusters defined by the cross-sectional unit of
observation, the district. For the sake of reference, from this point, the
inference and corresponding method from Condra et al., (2018) will be named
“UNIT”. This approach will be asymptotically valid in the presence of
dependence within district and heteroskedasticity as long as associations
across districts are ignorable. To provide some evidence on the plausibility
of this assumption, we provide Moran, (1950) $\mathcal{I}$ tests; see also
Kelejian and Prucha, (2001). To construct the tests, we first estimate model
(3.1) to obtain estimates of $U_{de}$ which are then used to construct the
observation specific score for $\widehat{\theta}$, the IV estimator of
$\theta_{0}$ based on (3.1)-(3.2). The scores are then used as the data input
into the Moran $\mathcal{I}$ test. We consider two formulations for the weight
matrix used in the Moran $\mathcal{I}$ test. In the first, we create the
weight matrix by assigning a value of one for the weight corresponding to the
two nearest neighbors for each observation on the basis of Euclidean distance
in the latitude and longitude coordinates of the centroid of each district and
set all other weights to zero. In the second, we create the weight matrix by
assigning a value of one to the weight corresponding to observations from the
same district and set all other weights to zero. Note that the first setup
should have power against alternatives where dependence is related to
geographic distance and the second against alternatives where dependence is
intertemporal.

In Table 1, we report the value of the Moran $\mathcal{I}$ test statistic and
the associated p-value for testing the null-hypothesis of zero correlation
across observations based on the weight matrix formed from geographic distance
and on the weight matrix for intertemporal dependence. For spatial weights, we
consider each election round separately as well as report results pooling both
election rounds. In each case, the calculated test statistic is large, and the
p-value for testing the null hypothesis of zero correlation is small. The
results are highly suggestive of the presence of both cross-sectional
dependence associated with geographic distance and intertemporal dependence.
As there is no immediately obvious unique way to group districts on the basis
of geographic proximity, it is potentially interesting to construct clusters
in a data-dependent fashion.

Table 1: Moran, (1950) $\mathcal{I}$ Tests based on IV Scores

| Moran $\mathcal{I}$ | p-value
---|---|---
Spatial weights - First election round | 4.338 |  $<$ 0.001
Spatial weights - Second election round | 2.546 | 0.011
Spatial weights - Both election rounds | 5.387 |  $<$ 0.001
Intertemporal weights | 2.755 | 0.006

  * •

Notes: This table provides Moran, (1950) $\mathcal{I}$ test results for the
Condra et al., (2018) example. The test statistic value and associated p-value
are respectively reported in the columns labeled “Moran $\mathcal{I}$” and
“p-value.” The first three rows report tests with weights for the Moran test
constructed on the basis of geographic distance as described in the text. The
first row uses only observations from the first election round, the second row
uses only observations from the second election round, and the third row uses
all observations to construct the test statistic. The final row uses a weight
matrix constructed solely on the basis of intertemporal distance as described
in the text.

###  3.2 Data-Dependent Cluster-Based Inference

We now apply the approach outlined in Section 2 to construct data-dependent
clusters of observations and produce inferential statements robust to the
presence of spatial and temporal correlation and heteroskedasticity. A key
input into this process is a dissimilarity measure $\mathrm{d}$ capturing the
“distance” between each pair of distinct districts, $d,d^{\prime}$. For this
application, we use dissimilarity defined by Euclidean distance in geography

$\displaystyle\mathrm{d}(d,d^{\prime})=\|L_{d}-L_{d^{\prime}}\|_{2}=\sqrt{(\mathrm{lat}_{d}-\mathrm{lat}_{d^{\prime}})^{2}+(\mathrm{long}_{d}-\mathrm{long}_{d^{\prime}})^{2}}$
(3.3)

where $L_{d}=(\mathrm{lat}_{d},\mathrm{long}_{d^{\prime}})$ consists of the
latitude and longitude coordinates of the centroid of the district and
$\|\cdot\|_{2}$ denotes Euclidean distance.

Because of the nonuniform locations of district centers (i.e. districts are
not on a rectangular grid), it is unclear ex ante how to partition them or how
many clusters to use. As one potential resolution to this ambiguity, we apply
$k$-medoids to the data using dissimilarity $\mathrm{d}$. Note that (3.3)
assigns dissimilarity zero to observations from the same district which means
that these observations will always be assigned to the same cluster. We thus
allow general dependence across election rounds within district as well as
capture dependence across election rounds and across districts that belong to
the same cluster. Specifically, for each $k=2,...,k_{\mathrm{max}}$ where
$k_{\mathrm{max}}$ is a user-specified upper bound on the number of clusters
to be considered, we apply $k$-medoids, yielding a set of partitions
$\mathscr{C}=\\{\mathcal{C}^{(2)},...,\mathcal{C}^{(k_{\mathrm{max}})}\\}$. As
an illustration, Figure 1 displays a partition of the data resulting from
applying $k$-medoids with $k=6$ with districts rendered with common markers
belonging to common clusters $\mathsf{C}\in\mathcal{C}^{(6)}$.333We display
the result with $k=6$ as this is the solution obtained in (2.9) in this
example for CRS. In this empirical illustration, we set $k_{\mathrm{max}}=8.$

Figure 1: Display of partition of 205 districts in Afghanistan by $k$-medoids
using final number of clusters given by $k=6$. Distances are Euclidean
distances based on latitude and longitude coordinates recorded at district
centroids. Different marks correspond to different clusters in the partition.
Marks are plotted at district centroids. Districts without symbols have
insufficient data and were excluded in the analysis Condra et al., (2018) as
well as in the current analysis.

In practice, we recommend choosing a small value for $k_{\mathrm{max}}$.This
recommendation is based on two considerations. First, IM and CRS, our
preferred procedures, rely on within-cluster estimates and the number of
observations per cluster decreases quickly as the number of clusters
considered increases, which results in unstable performance when large numbers
of clusters are considered. Second, for dependence structures with non-trivial
dependence, size distortions typically become pronounced with larger numbers
of groups while improvements to power tend to rapidly diminish once a modest
number of groups are considered.

With the set of potential partitions, $\mathscr{C}$, in hand, the final inputs
we need to be able to select a partition and testing procedure using (2.9) are
feasible guesses for Type I and average Type II error rates,
$\widehat{\text{Err}}_{\text{Type-I}}(\bullet(a),k)$ and
$\widehat{\text{Err}}_{\text{Type-II},\mathbf{P}_{\text{alt}}}(\bullet(a),k)$.
Let $U$ and $V$ respectively denote the $410\times 1$ vectors produced by
stacking $U_{de}$ and $V_{de}$. To obtain
$\widehat{\text{Err}}_{\text{Type-I}}(\bullet(a),k)$ and
$\widehat{\text{Err}}_{\text{Type-II},\mathbf{P}_{\text{alt}}}(\bullet(a),k)$,
we use a parametric covariance model

$\displaystyle\begin{split}{\mathrm{E}}\begin{bmatrix}UU^{\prime}&UV^{\prime}\\\
VU^{\prime}&VV^{\prime}\end{bmatrix}&=\begin{bmatrix}\Sigma_{U}(\tau^{U})&\rho(\Sigma_{U}(\tau^{U})^{1/2})(\Sigma_{V}(\tau^{V})^{1/2})^{\prime}\\\
\rho(\Sigma_{V}(\tau^{V}))^{1/2}(\Sigma_{U}(\tau^{U})^{1/2})^{\prime}&\Sigma_{V}(\tau^{V})\end{bmatrix}\\\
&=\Sigma(\tau^{U},\tau^{V},\rho)\end{split}$ (3.4)

with parameters $(\tau^{U},\tau^{V},\rho)$ where $\Sigma_{U}(\tau)$ and
$\Sigma_{V}(\tau)$ are covariance matrices with entries given by

$\displaystyle\left[\Sigma_{\varepsilon}(\tau)\right]_{de,d^{\prime}e^{\prime}}=f_{\tau}((L_{d},e),(L_{d}^{\prime},e^{\prime}))=\exp(\tau_{1})\exp(-\tau_{2}^{-1}\|L_{d}-L_{d^{\prime}}\|_{2}-\tau_{3}^{-1}|e-e^{\prime}|)$
(3.5)

for $\varepsilon\in\\{U,V\\}$. We obtain calibrated values for these
parameters, $(\widehat{\tau}^{U},\widehat{\tau}^{V},\widehat{\rho})$ by
forming initial estimates $\widehat{U}_{de}$ and $\widehat{V}_{de}$ of
$U_{de}$ and $V_{de}$ from IV estimation of (3.1)-(3.2) and then using these
residuals to estimate the covariance model via Gaussian QMLE. We can then
approximate the exact Type I and average Type II error rates that would result
from applying each inferential procedure with each potential partition within
a Gaussian model with covariance matrix given by
$\Sigma(\widehat{\tau}^{U},\widehat{\tau}^{V},\widehat{\rho})$ by simulation.
Given these estimates of the Type I and Type II error rates, we then solve
(2.9) to obtain the final partition and level for defining the cutoff in the
decision rule used in the data. Further details are provided in Appendix A.2.

Table 2: Impact of Early Morning Attacks on Voter Turnout during the 2014
Election

Cluster-based Inference
---
|  $\phantom{\Big{|}}\widehat{\theta}_{0}$ |  $s.e.$ |  $t$-stat | C.I. | C.I. (usual CV) |  ${\widehat{k}}$
UNIT |  -0.145 |  0.061 |  2.385 |  ( -0.312, |  0.022 ) |  ( -0.265, |  -0.025 ) |  205
CCE |  -0.145 |  0.090 |  1.609 |  ( -0.404, |  0.114 ) |  ( -0.358, |  0.068 ) |  8
IM |  -0.122 |  0.432 |  0.282 |  ( -1.320, |  1.077 ) |  ( -1.320, |  1.077 ) |  5
CRS |  -0.242 |  |  |  ( -1.497, |  0.084 ) |  ( -1.497, |  0.084 ) |  6

  * •

Notes: This table presents inferential results based on selected clusters. Row
labels indicate which procedure is used. The column labeled
$\widehat{\theta}_{0}$ reports the IV estimate of $\theta_{0}$ for the full
sample in the rows labeled UNIT and CCE and the average of IV estimators of
$\theta_{0}$ of the five or six data-generated clusters in the rows labeled IM
and CRS respectively. Column $s.e.$ reports the estimated standard errors
obtained in each procedure. Note that CRS does not rely on an explicit
standard error estimate. Column $t$-stat reports the $t$-statistic for testing
the null hypothesis that $\theta_{0}=0$ for each of the procedures. Column C.I
reports confidence intervals for $\theta_{0}$ obtained using
$\widehat{\alpha}$ and $\widehat{k}$ from (2.9) for each procedure. Column
C.I. (usual CV) reports confidence intervals of the IV estimate of
$\theta_{0}$ using fixed asymptotic thresholds. Column ${\widehat{k}}$
indicates the number of clusters selected in each procedure.

Table 2 reports inferential results for clustering at the district level as in
Condra et al., (2018) (UNIT) and for the data-dependent CCE, IM, and CRS
cluster-based inferential procedures. We report results for all procedures for
completeness and illustration, but note that the theory and simulation results
lead us to recommend that one uses IM or CRS in practice. For each method, we
report the corresponding estimate of $\theta_{0}$, the standard error,
$t$-statistic, 95% confidence interval using $\widehat{k}$ and
$\widehat{\alpha}$ from (2.9), 95% confidence interval using $\alpha=.05$ and
$\widehat{k}$, and $\widehat{k}$.

In this illustration, all 95% level confidence intervals include 0 with the
exception of the interval produced by district level clustering with the
asymptotic cutoff. That is, all of the considered procedures which allow for
spatial correlation would lead to failure to reject the null hypothesis that
insurgent attacks have no effect on voter turnout at the 5% level within the
posited linear IV structure. There is substantive variation in the interval
estimates provided by the various procedures more generally though.

Both procedures based on clustered standard errors (UNIT and CCE), while
producing intervals that cover zero when spatial correlation is allowed for,
are relatively narrow when compared to IM and CRS. For both UNIT and CCE, we
can see a mild discrepancy between the intervals produced under the asymptotic
approximation and those produced by allowing a data-dependent decision
threshold based on the parametric model. Specifically, within the parametric
model underlying the simulation, controlling size at the 5% level would
correspond to rejecting the null only when the p-value is smaller than 0.006
or 0.022 for UNIT and CCE respectively (using the fixed number of groups
approximation of Bester et al., (2011).) These differences suggest that the
asymptotic approximation is potentially failing and suggest some caution
should be made in relying on these results; see also Ferman, (2019). We do
note that the interval provided in column C.I. will exactly control size under
the specified parametric model. However, we also find that both CCE and
individual level clustering are relatively sensitive to misspecification of
the parametric model underlying construction of the estimated Type I and Type
II error rates when compared to IM and CRS in the simulation results reported
in Section 4.

Based on the theoretical and simulation results reported in the next sections,
our preferred procedures are IM and CRS. These procedures are provably valid
asymptotically under relatively weak conditions and are also relatively robust
to misspecification of the parametric model underlying the size and power
calculations in our simulation experiments. Looking to IM and CRS, we see
similar performance, in the sense of relatively wide 95% interval estimates,
across both, though IM is substantially wider than CRS in this example. For
both of these procedures, the asymptotic cutoffs coincide with those produced
by solving (2.9) using the simulated Type I and Type II error rates, which is
reassuring though not necessary for validity. We note that there is a danger
of a researcher trying both IM and CRS (and perhaps CCE) and reporting the
most optimistic result. For this reason, our recommendation is to either
default to, and compute only, IM or CRSor to report results from all
procedures to provide complete information to readers.

##  Simulation

This section examines the finite sample performance of inference based on
data-dependent clusters in a series of simulation experiments roughly based
upon the empirical illustration of Section 3. We present results for inference
on a coefficient in a linear model with all exogenous variables and results
for inference on the coefficient on an endogenous variable in a linear IV
model. We refer to the first case as the “OLS simulation” and to the second
case as the “IV simulation.” We present the data generating process (DGP) for
each case and outline the inferential procedures we consider in Subection 4.1.
We then present results in Subection 4.2.

###  4.1 Data Generating Processes and Inferential Procedures

For both the OLS and IV simulation, we generate observations indexed by $de$
for $d=1,...,205$ and $e=1,2$ where $d$ represents the cross-sectional unit
and $e$ the time period.444Additional results with 820 cross-sectional units
and two time periods are provided in a supplementary appendix. Each
observation is thus associated with a spatial location given by
$L_{d}=(\mathrm{lat}_{d},\mathrm{long}_{d^{\prime}})$, the latitude and
longitude of the centroid of a district from the observed data in the
empirical example, and a temporal index $e\in\\{1,2\\}$.

OLS Simulation DGP: For the OLS simulation, we generate observed data,
$\mathscr{D}=\\{Y_{de},X_{de},W_{de}\\}_{d=1,...,205,e=1,2}$ according to

$\displaystyle Y_{de}=\theta_{0}X_{de}+W_{de}^{\prime}\gamma_{0}+U_{de},$
(4.1)

where $X_{de}$ is the variable of interest, $W_{de}$ is a $10\times 1$ vector
of control variables, and $U_{de}$ is unobservable. We set coefficients
$\theta_{0}=0$ and $\gamma_{0}=0_{10}=(0,\dots,0)^{\prime}$.

We condition on a single realization of
$\\{X_{de},W_{de}\\}_{d=1,\dots,205,e=1,2}$ generated as follows: $\forall\
l\in\\{1,\dots,10\\}$, $\forall\ m\neq l$, $\forall\
(d,e)\in\\{1,\dots,205\\}\times\\{1,2\\}$, $\forall\
(d^{\prime},e^{\prime})\neq(d,e)$,

$\displaystyle\begin{split}&X_{de}\sim\mathrm{N}(0,1);W_{del}\sim\mathrm{N}(0,1),\\\
&\mathrm{corr}(X_{de},X_{d^{\prime}e^{\prime}})=\mathrm{corr}(W_{del},W_{d^{\prime}e^{\prime}l})=f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime})),\\\
&\mathrm{corr}(X_{de},W_{del})=\mathrm{corr}(W_{del},W_{dem})=0.5,\\\
&\mathrm{corr}(X_{de},W_{d^{\prime}e^{\prime}l})=\mathrm{corr}(W_{del},W_{d^{\prime}e^{\prime}m})=0,\end{split}$
(4.2)

for $f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime}))$ defined in (3.5) with
$\tau=(0,3,1)^{\prime}$. Unobservables $U_{de}$ are drawn independently across
simulation replications. We consider two settings for the distribution of
$U_{de}$:

  * A.

Homogeneous exponential covariance (BASELINE).

$\displaystyle
U_{de}\sim\mathrm{N}(0,1);\quad\mathrm{corr}(U_{de},U_{d^{\prime}e^{\prime}})=f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime})),$
(4.3)

for $f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime}))$ defined in (3.5) with
$\tau=(0,3,1)^{\prime}$.

  * B.

Spatial auto-regression (SAR).

$\displaystyle\begin{split}&U_{de}=0.15\textstyle\sum_{d^{\prime}\neq
d}U_{d^{\prime}e}\textbf{1}_{\\{\|L_{d}-L_{d^{\prime}}\|_{2}<0.3\\}}+\varepsilon_{de},\\\
&\varepsilon_{de}\sim\mathrm{N}(0,1);\quad\mathrm{corr}(\varepsilon_{d1},\varepsilon_{d2})=\exp(-1),\
\mathrm{corr}(\varepsilon_{de},\varepsilon_{d^{\prime}e^{\prime}})=0\text{ for
}d\neq d^{\prime}.\end{split}$ (4.4)

IV Simulation DGP: For the IV simulation, observed data are
$\mathscr{D}=\\{Y_{de},X_{de},W_{de},Z_{de}\\}_{d=1,...,205,e=1,2}$ with
outcomes $Y_{de}$ and endogenous variable $X_{de}$ generated from the system
of equations

$\displaystyle\begin{split}Y_{de}&=\theta_{0}X_{de}+W_{de}^{\prime}\gamma_{0}+U_{de},\\\
X_{de}&=\pi_{0}Z_{de}+W_{de}^{\prime}\xi_{0}+V_{de}\end{split}$ (4.5)

where $X_{de}$ is the endogenous variable of interest, $W_{de}$ is a $10\times
1$ vector of control variables, $Z_{de}$ is a scalar instrumental variable,
and $U_{de}$ and $V_{de}$ are structural unobservables. We set coefficients
$\theta_{0}=0$, $\pi_{0}=2$, and
$\gamma_{0}=\xi_{0}=0_{10}=(0,\dots,0)^{\prime}$.

We condition on a single realization of the exogenous variables,
$\\{Z_{de},W_{de}\\}_{d=1,...,205,e=1,2}$, generated as follows: $\forall\
l\in\\{1,\dots,10\\}$, $\forall\ m\neq l$, $\forall\
(d,e)\in\\{1,\dots,205\\}\times\\{1,2\\}$, $\forall\
(d^{\prime},e^{\prime})\neq(d,e)$,

$\displaystyle\begin{split}&Z_{de}\sim\mathrm{N}(0,1);W_{del}\sim\mathrm{N}(0,1),\\\
&\mathrm{corr}(Z_{de},Z_{d^{\prime}e^{\prime}})=\mathrm{corr}(W_{del},W_{d^{\prime}e^{\prime}l})=f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime})),\\\
&\mathrm{corr}(Z_{de},W_{del})=\mathrm{corr}(W_{del},W_{dem})=0.5,\\\
&\mathrm{corr}(Z_{de},W_{d^{\prime}e^{\prime}l})=\mathrm{corr}(W_{del},W_{d^{\prime}e^{\prime}m})=0,\end{split}$
(4.6)

for $f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime}))$ defined in (3.5) with
$\tau=(0,3,1)^{\prime}$. We consider two scenarios for the unobservables
$(U_{de},V_{de})$ which are drawn independently across simulation
replications:

  * A.

Homogeneous exponential covariance (BASELINE).

$\displaystyle\begin{split}&U_{de}\sim\mathrm{N}(0,1),\
V_{de}\sim\mathrm{N}(0,1),\\\
&\mathrm{corr}(U_{de},U_{d^{\prime}e^{\prime}})=\mathrm{corr}(V_{de},V_{d^{\prime}e^{\prime}})=f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime})),\\\
&\mathrm{corr}(U_{de},V_{de})=0.8,\\\
&\mathrm{corr}(U_{de},V_{d^{\prime}e^{\prime}})=0\text{ for
}(d^{\prime},e^{\prime})\neq(d,e),\end{split}$ (4.7)

for $f_{\tau}((L_{d},e),(L_{d^{\prime}},e^{\prime}))$ defined in (3.5) with
$\tau=(0,3,1)^{\prime}$.

  * B.

Spatial auto-regression (SAR).

$\displaystyle\begin{split}&U_{de}=0.15\textstyle\sum_{d^{\prime}\neq
d}U_{d^{\prime}e}\textbf{1}_{\\{\|L_{d}-L_{d^{\prime}}\|_{2}<0.3\\}}+\varepsilon_{de},\\\
&V_{de}=0.15\textstyle\sum_{d^{\prime}\neq
d}V_{d^{\prime}e}\textbf{1}_{\\{\|L_{d}-L_{d^{\prime}}\|_{2}<0.3\\}}+\eta_{de},\\\
&\varepsilon_{de}\sim\mathrm{N}(0,1),\ \eta_{de}\sim\mathrm{N}(0,1),\\\
&\mathrm{corr}(\varepsilon_{d1},\varepsilon_{d2})=\mathrm{corr}(\eta_{d1},\eta_{d2})=\exp(-1),\\\
&\mathrm{corr}(\varepsilon_{de},\eta_{de})=0.8,\\\
&\mathrm{corr}(\varepsilon_{de},\varepsilon_{d^{\prime}e^{\prime}})=\mathrm{corr}(\eta_{de},\eta_{d^{\prime}e^{\prime}})=0\text{
for }d^{\prime}\neq d,\\\
&\mathrm{corr}(\varepsilon_{de},\eta_{d^{\prime}e^{\prime}})=0\text{ for
}(d^{\prime},e^{\prime})\neq(d,e).\end{split}$ (4.8)

Inferential Procedures: Within each of the four simulation designs defined
above, we generate 1000 simulation replications. We then report results for
point estimation and for inference, focusing on size and power of hypothesis
tests, about the parameter $\theta_{0}$ based on the following procedures:

  * 1.

SK. Inference based on the spatial HAC estimator from Sun and Kim, (2015) with
bandwidth selection adapted from the proposal of Lazarus et al., (2018).

  * 2.

UNIT-U. Inference based on the cluster covariance estimator with clusters
defined as the cross-sectional unit of observation with critical value from a
$t$-distribution with $k-1=204$ degrees of freedom.

  * 3.

UNIT. Inference based on the cluster covariance estimator with clusters
defined as the cross-sectional unit of observation and critical value obtained
by solving (2.9) with cluster structure given by unit-level clustering.

  * 4.

CCE. Inference based on the cluster covariance estimator with clusters and
rejection threshold obtained by solving (2.9) as described in Section 2.

  * 5.

IM. Inference based on IM with clusters and rejection threshold obtained by
solving (2.9) as described in Section 2.

  * 6.

CRS. Inference based on CRS with clusters and rejection threshold obtained by
solving (2.9) as described in Section 2.

For UNIT, CCE, IM, and CRS, we obtain preliminary estimates of unobserved
components $U_{de}$ in the OLS setting and $(U_{de},V_{de})$ in the IV
setting. We then apply Gaussian QMLE with the covariance structure specified
in BASELINE with free parameters using these preliminary estimates as data to
obtain the structure for use in simulating the Type I and Type II error rates
to allow feasible solution of (2.9). Thus, results for UNIT, CCE, IM, and CRS
in the BASELINE OLS and BASELINE IV setting illustrate performance when tuning
parameters result from solving (2.9) with a correctly specified model for the
covariance structure with feasible estimates of the covariance parameters. In
contrast, results from the SAR OLS and SAR IV setting illustrate performance
when tuning parameters are obtained by solving (2.9) with a misspecified
model. We provide detailed descriptions of the implementation of SK, UNIT,
CCE, IM, and CRS in Appendix A.2.

###  4.2 Simulation Results

We report size of 5% level tests as well as provide results on point
estimation quality in Table 3. In the OLS simulation, we obtain point
estimates by applying OLS to estimate the parameters of (4.1), and we obtain
point estimates in the IV simulation by applying IV to estimate the parameters
in (4.5). Recall that both IM and CRS rely on first obtaining within cluster
estimates and have natural point estimator defined by
$\frac{1}{\widehat{k}}\sum_{\mathsf{C}\in\widehat{\mathcal{C}}}\widehat{\theta}_{\mathsf{C}}$
where $\widehat{k}$ and $\widehat{\mathcal{C}}$ are the results from solving
(2.9) for each procedure and $\widehat{\theta}_{\mathsf{C}}$ is a point
estimator that uses only the observations in cluster $\mathsf{C}$. In terms of
point estimation properties, we report bias and root mean square error (RMSE)
for the OLS simulation. For the IV simulation, we report median bias and
median absolute deviation (MAD) due to the IV estimator’s lack of finite
sample moments.

Table 3: Simulation Results

| OLS |  | IV
---|---|---|---
Method |  Bias |  RMSE |  Size |  |  Median Bias |  MAD |  Size
A. BASELINE
SK |  0.015 |  0.337 |  0.381 |  |  0.002 |  0.114 |  0.362
UNIT-U |  |  |  0.577 |  |  |  |  0.568
UNIT |  |  |  0.047 |  |  |  |  0.074
CCE |  |  |  0.046 |  |  |  |  0.062
IM |  0.014 |  0.213 |  0.044 |  |  -0.059 |  0.091 |  0.046
CRS |  0.014 |  0.213 |  0.042 |  |  -0.061 |  0.095 |  0.040
B. SAR
SK |  -0.013 |  0.866 |  0.447 |  |  0.002 |  0.280 |  0.324
UNIT-U |  |  |  0.639 |  |  |  |  0.502
UNIT |  |  |  0.359 |  |  |  |  0.256
CCE |  |  |  0.047 |  |  |  |  0.083
IM |  -0.006 |  0.385 |  0.038 |  |  -0.046 |  0.158 |  0.025
CRS |  -0.003 |  0.354 |  0.049 |  |  -0.095 |  0.213 |  0.041

  * •

Notes: Results from the OLS Simulation (in columns labeled “OLS”) and IV
Simulation (in columns labeled “IV”) described in Section 4.1. Row labels
indicate inferential method. Panel A corresponds to the BASELINE design where
tuning parameters for UNIT, CCE, IM, and CRS are selected based on calculating
size and power from feasible estimates of a correctly specific parametric
model. Panel B corresponds to the SAR design where tuning parameters for UNIT,
CCE, IM, and CRS are selected based on calculating size and power from
feasible estimates of a misspecified parametric model. For the OLS designs, we
report the bias and RMSE of the point estimator associated with each procedure
along with size of 5% level tests. For the IV designs, we report the median
bias and MAD of the point estimator associated with each procedure along with
size of 5% level tests.

The main feature of the results presented in Table 3 related to the proposal
in our paper is that both IM and CRS with data dependent groups and rejection
threshold determined by (2.9) control size relatively well across the four
considered designs. CCE with data dependent groups and rejection threshold
determined by (2.9) also does a reasonable job controlling size across the
designs, though it does somewhat less well than IM and CRS in the IV setting
and is also not covered by our theoretical results. The performance of these
procedures is similar in both BASELINE settings, where tuning parameters
$\widehat{k}$ and $\widehat{\alpha}$ are obtained using a correctly specified
covariance model, and in the SAR setting, where the problem solved to obtain
tuning parameters makes use of a misspecified covariance model.

The robustness of IM, CRS, and, to a lesser extent, CCE is not exhibited by
the remaining procedures. The poor behavior of UNIT-U, which ignores spatial
dependence entirely, is unsurprising. More surprising is that SK, which
attempts to account for spatial dependence, also does poorly across all
designs considered. We suspect that adapting SK to choose tuning parameters on
the basis of solving a problem similar to (2.9) would perform much better, but
wished to compare to a benchmark from the existing literature. The performance
of UNIT is interesting. This procedure treats cross-sectional units as
spatially uncorrelated in constructing the standard error estimator but then
uses a parametric model that has spatial correlation to adjust the decision
threshold for rejecting a hypothesis according to (2.9). In this case, we see
that using a correctly specified parametric model for this adjustment restores
size control but that size is not controlled under the misspecified parametric
structure. This behavior is in line with our theoretical results which rely on
the use of a small number of clusters to maintain size control while allowing
for relatively general dependence structures and not requiring correct
specification of the parametric model used for tuning parameter choice. The
simulation clearly demonstrates that robustness to misspecification is not
maintained when large numbers of clusters are used.

It is also interesting to look at the point estimation results, though these
mostly mirror results already available in the literature; see, e.g.,
Ibragimov and Müller, (2010), Bester et al., (2011), and Conley et al.,
(2018). Specifically, we see that both IM and CRS dominate the full sample OLS
estimator in terms of both bias and RMSE in our simulation designs. That is,
there appears to be a gain, in terms of both point estimation properties and
size control, to using the IM or CRS procedure. The results are more muddy in
the IV simulation where the full sample estimator exhibits lower median bias
at the cost of larger MAD and poorer size control. Our preference is still for
IM or CRS as having both more accurate statistical inference and more precise
(as measured by MAD) estimation seems worth the cost of larger bias. We also
note again that these properties appear across simulation studies reported in
the literature but are predicated on using a relatively small number of groups
as the performance of IM and CRS become erratic when small numbers of
observations are used to form the within-group estimators. See Conley et al.,
(2018) for further discussion.

(a) BASELINE (b) SAR

Figure 2: OLS power curves.

(a) BASELINE (b) SAR

Figure 3: IV power curves.

We report power curves for 5% level tests of the hypothesis
$H_{0}:\theta_{0}=\theta^{\mathrm{alt}}$ for alternatives
$\theta^{\mathrm{alt}}$ produced by the different procedures across the
simulations in Figures 3 and 3. In these figures, the horizontal axis gives
the hypothesized value, $\theta^{\mathrm{alt}}$, so size of the test is
captured by the point $\theta^{\mathrm{alt}}=0$. Figure 3 presents the results
from the OLS simulation. Here, we see that the power curves are symmetric and
that the highest power among procedures that control size is obtained by IM
and CRS, both of which perform similarly. Looking to the IV results in Figure
3, we see that power curves are asymmetric and slightly shifted due to the
finite sample behavior of the IV estimator. We also see that there is no
longer a clear picture about which of the procedures that controls size
performs better in terms of power. Specifically, each of CCE, IM, and CRS
exhibits higher power over different sets of alternative values. Exploring
these tradeoffs more deeply may potentially be interesting but is beyond the
scope of this paper. Overall, these figures reinforce the takeaways from the
size and point estimation results presented in Table 2 which suggest that IM
and CRS provide good default procedures which are robust, both in simulation
and in the formal analysis, easy-to-compute, and not clearly dominated by
other commonly used approaches.

Table 4: Distribution of $\widehat{k}$

|  $\widehat{k}=2$ |  3 |  4 |  5 |  6 |  7 |  8
---|---|---|---|---|---|---|---
A. OLS - BASELINE
CCE |  0.000 |  0.000 |  0.001 |  0.098 |  0.311 |  0.260 |  0.330
IM |  0.000 |  0.000 |  0.000 |  0.007 |  0.018 |  0.113 |  0.862
CRS |  0.000 |  0.000 |  0.000 |  0.000 |  0.001 |  0.099 |  0.900
B. OLS - SAR
CCE |  0.000 |  0.001 |  0.022 |  0.051 |  0.113 |  0.268 |  0.545
IM |  0.000 |  0.003 |  0.069 |  0.554 |  0.040 |  0.176 |  0.158
CRS |  0.000 |  0.000 |  0.000 |  0.000 |  0.001 |  0.401 |  0.598
C. IV - BASELINE
CCE |  0.000 |  0.001 |  0.001 |  0.039 |  0.242 |  0.261 |  0.456
IM |  0.000 |  0.000 |  0.005 |  0.051 |  0.062 |  0.066 |  0.816
CRS |  0.000 |  0.000 |  0.000 |  0.000 |  0.000 |  0.157 |  0.843
D. IV - SAR
CCE |  0.006 |  0.027 |  0.048 |  0.047 |  0.071 |  0.257 |  0.544
IM |  0.005 |  0.078 |  0.240 |  0.568 |  0.030 |  0.038 |  0.041
CRS |  0.000 |  0.000 |  0.000 |  0.000 |  0.032 |  0.480 |  0.488

  * •

Notes: Simulation results for $\widehat{k}$, the data-dependent number of
clusters obtained in (2.9). Panels correspond to the different designs
described in Section 4.1. Each entry reports the simulation probability of
$\widehat{k}$ being equal to the value given in the column label for the
procedure indicated in the row label.

We report properties of the data dependent number of clusters obtained from
(2.9), $\widehat{k}$, across our simulation designs in Table 4. We first note
that 5% level tests based on CRS have trivial power (power equal size) when
based on fewer than six groups, so the number of groups selected for CRS is
always chosen to be six or greater. We also see that, with two exceptions,
$\widehat{k}=8$ is the most likely selection. The exceptions are for IM in the
OLS and IV settings with SAR covariance structure where five groups are
selected in 55.4% and 56.8% of the simulation replications, respectively. We
also note that, while $\widehat{k}=8$ is the most common solution, there is
non-trivial weight on other numbers of groups for all procedures. Finally, we
wish to reiterate that the grouping structure is not literally correct in any
of our designs but is being used as a feasible way to downweight small
covariances to allow for the construction of informative and robust
inferential statements. Thus, the goal here is not to find the “correct”
number of clusters.

Table 5: Distribution of $\widehat{\alpha}$

| Quantile
---|---
|  0.1 |  0.25 |  0.5 |  0.75 |  0.9
A. OLS - BASELINE
CCE |  0.004 |  0.006 |  0.008 |  0.010 |  0.012
IM |  0.039 |  0.042 |  0.046 |  0.050 |  0.050
CRS |  0.039 |  0.043 |  0.047 |  0.050 |  0.050
B. OLS - SAR
CCE |  0.021 |  0.024 |  0.027 |  0.031 |  0.034
IM |  0.047 |  0.049 |  0.050 |  0.050 |  0.050
CRS |  0.047 |  0.047 |  0.050 |  0.050 |  0.050
C. IV - BASELINE
CCE |  0.004 |  0.005 |  0.007 |  0.008 |  0.010
IM |  0.045 |  0.048 |  0.050 |  0.050 |  0.050
CRS |  0.039 |  0.039 |  0.047 |  0.050 |  0.050
D. IV - SAR
CCE |  0.014 |  0.020 |  0.025 |  0.029 |  0.033
IM |  0.032 |  0.050 |  0.050 |  0.050 |  0.050
CRS |  0.008 |  0.043 |  0.050 |  0.050 |  0.050

  * •

Notes: Simulation results for $\widehat{\alpha}$, the data-dependent p-value
threshold for rejecting a hypothesis at the 5% level obtained in (2.9). Panels
correspond to the different designs described in Section 4.1. For each design
and DGP we report the .1, .25, .5, .75, and .9 quantiles of the simulation
distribution of $\widehat{\alpha}$.

Finally, we summarize the simulation distribution of the data dependent
p-value threshold for rejecting a 5% level test obtained from (2.9),
$\widehat{\alpha}$, in Table 5. Recall that using the asymptotic approximation
underlying each procedure would correspond to rejecting with a threshold of
.05. Here we see that for CCE, the distribution of the threshold that would be
used to provide exact size control at the 5% level is substantially shifted
away from .05 with most of the mass of the distribution over values that are
much smaller than .05. We note that this suggests the asymptotic approximation
does not provide a particularly reliable guide to the performance of inference
based on CCE in our settings which is likely due to strong departures from
homogeneity assumptions that are used to establish the asymptotic behavior of
CCE with small numbers of groups. For both IM and CRS, we see substantially
smaller shifts of $\widehat{\alpha}$ away from .05. Indeed, for both SAR
cases, the distribution of $\widehat{\alpha}$ for IM and CRS has a substantial
mass point at .05. In the BASELINE cases, solving (2.9) for both IM and CRS
does lead to systematic use of p-value thresholds that are smaller than .05 to
achieve 5% level size control. While still noticeable, these departures are
significantly less extreme than when considering inference based on CCE.
Again, this behavior is consistent with the more stable and robust performance
of IM and CRS relative to CCE. It also highlights how having two tuning
parameters is useful for maintaining finite sample properties as size control
may not be achievable in finite samples if one may only choose the number of
groups for clustering and is unable to adjust the decision threshold.

Overall, the simulation results suggest that IM and CRS with data dependent
clusters constructed as outlined in Section 2 control size in the simulation
designs we consider while maintaining non-trivial power. In the remaining
sections, we establish theoretical results that demonstrate that these
procedures provide asymptotically valid inference under reasonable conditions.

##  Formal Analysis of Inference with Learned Clusters

This section presents a formal analysis of the procedure defined by Algorithm
1.

The analysis of Algorithm 1 requires preliminary theoretical groundwork. This
section begins by developing a notion of regularity for dissimilarity measures
in Subsection 5.1. Subsections 5.2–5.6 then define an asymptotic theory which
is helpful for the analysis of Algorithm 1. The regularity conditions in those
subsections are sufficient for establishing a mixing central limit theorem to
hold, as well as establishing relevant properties of $k$-medoids.

Subsections 5.1–5.6 are organized such that they build a ladder of
increasingly higher-level regularity conditions. The penultimate rung of this
ladder, Theorem 5, concludes with concrete properties for
$\widehat{\theta}_{\mathsf{C}}$ for a sequence of linear models in which
Algorithm 1 with $\mathsf{C}\in\mathscr{C}$ obtained by $k$-medoids has been
applied. The conclusions of Theorem 5 are then used as a final rung of high-
level conditions in Subections 5.8 and 5.9 to allow application of Theorems 6
and 7 which analyze size control for IM and CRS procedures with inputs more
general than those arising from Algorithm 1. Versions of Theorems 6 and 7 were
previously established in Ibragimov and Müller, (2010) and Canay et al.,
(2017); and thus the primary theoretical contribution in this paper is the
development of the theory in Sections 5.1–5.6. Theorems 6 and 7 do contain
refinements relative to previous results which allow simultaneous analysis
using multiple partitions. Sections 5.10 and 5.11 give discussion of
properties of inference based on CCE as well as issues surrounding uniformity
of the high-level conditions.

###  5.1 A notion of regularity for dissimilarity measures

This section defines formally a measure of regularity for a spatial indexing
set $\mathsf{X}$ with dissimilarity $\mathrm{d}$. The most basic condition on
$\mathsf{X}$ is that the triangle inequality holds, in other words, each
$\mathsf{X}$ is a genuine metric space. The next definition is used to
restrict the growth rate of balls centered around any point of $\mathsf{X}$.
Let $|\mathsf{X}|$ be the cardinality of $\mathsf{X}$. Define $\mathsf{X}$ to
be $(C,\delta)$-finite-Ahlfors regular if it is a metric space and
$\min(|\mathsf{X}|,C^{-1}r^{\delta})\leqslant|\textsf{B}_{\mathsf{X},r}(i)|\leqslant\max(Cr^{\delta},1)$
for any $r>0$, for any $i$ in any $\mathsf{X}$, where
$\textsf{B}_{\mathsf{X},r}(i)$ is the $r$-ball centered at $i$ in the space
$\mathsf{X}$. Whenever $\mathsf{X}$ is both a metric space and satisfies the
above definition with constants $C,r$, write
$\mathsf{X}\in\textbf{Ahlf}_{C,\delta}$.

$(C,\delta)$-finite-Ahlfors regularity is a modification of the notion of
Ahlfors regularity encountered in the theory of metric spaces equipped with a
Borel measure $\mu$, in which case $|\textsf{B}_{\mathsf{X},r}(i)|$ is
replaced by $\mu(\textsf{B}_{\mathsf{X},r}(i))$ and the conditions
$\min(|\mathsf{X}|,...)$ and $\max(...,1)$ are dropped. This notion has
advantages relative to assuming that $\mathsf{X}$ be a subset of a Euclidean
space. First, the definition refers only to intrinsic properties of the space.
Second, this simple condition is sufficient both for realizing mixing central
limit theorems and for analyzing clustering.555The theory considers
$\mathsf{X}$ such that $(\mathrm{d}(i,j)=0)\Rightarrow(i=j)$, as this is an
immediate implication of the condition that
$\mathsf{X}\in\mathbf{Ahlf}_{C,\delta}$. To rectify this convention with the
data structure of the empirical example and simulation study in the previous
sections, it is necessary either to consider observations
$(Y_{d1},X_{d1},W_{d1},Z_{d1}),(Y_{d2},X_{d2},W_{d2},Z_{d2})$ as a common
observation arising from a single index $d$, or to consider an expanded
definition of $\textbf{Ahlf}_{C,\delta}$ which allows finitely many spatial
indexes to exist with mutual distance 0 to each other. Both possibilities are
straightforward and can be carried without any essential changes to the
arguments presented in this paper. For brevity, the details are omitted.

Not every metric dissimilarity measure $\mathsf{X}$ satisfying
$(C,\delta)$-finite-Ahlfors regularity admits an isometric embedding into
$\mathbb{R}^{\nu}$ for some $\nu$. In fact, the condition is not even
sufficient to guarantee that $\mathsf{X}$ admits a bi-Lipschitz embedding
(defined so that the maximum distortion is bounded by a constant) into some
$\mathbb{R}^{\nu}$. However, the condition is strong enough to ensure that
$\mathsf{X}$ can be “regularized.” The new space
$(\mathsf{X},\mathrm{d}^{1-\varepsilon})$, in which the exponent
$1-\varepsilon$ is applied element-wise to $\mathrm{d}$, is a metric space for
all $\varepsilon\in(0,1)$ and is called the $\varepsilon$-snowflake of
$\mathsf{X}$. The exponent $\varepsilon$ serves to regularize the distance
$\mathrm{d}$ so that it can be embedded into $\mathbb{R}^{\nu}$ with bounded
distortion. This follows from the fact that the doubling constant of finite
Ahlfors spaces can be effectively controlled, after which Assoud’s Embedding
Theorem (Assoud,, 1977) applies. More details are given in the supplemental
appendix.

###  5.2 Asymptotics A: Frames defined by an increasing sequence of
dissimilarity measures

The analysis of Algorithm 1 in relies on asymptotic theory. In this paper, the
task of indexing of data is delegated to a spatial indexing set. It is
therefore convenient to build an asymptotic frame by first deliniating the
classes of sequences of spatial indexing sets considered. Given the definition
of $(C,\delta)$-finite-Ahlfors regularity in the previous subsection, a
natural class of sequences are of the form
$\mathbb{N}\rightarrow\mathbf{Ahlf}_{C,\delta}$.

A sequence of finite dissimilarity measures will be denoted with
${\mathsf{X}}_{\shortrightarrow}=(\mathsf{X}_{n})_{n\in\mathbb{N}}$. Note that
in this section, unlike in the earlier sections, $n$ appears as an explicit
subscript in $\mathsf{X}_{n}$ precisely because $\mathsf{X}_{n}$ are elements
of a newly introduced sequence. In the earlier sections, $n$ was not needed to
define quantities like $\widehat{T}$. From this point, any $\mathsf{X}_{n}$
appearing with a subscript always implicitly belongs to some sequence
$\mathsf{X}_{\shortrightarrow}$.

######  Condition  1.

(Ahlfors Regularity) The sequence of dissimilarity measures
${\mathsf{X}}_{\shortrightarrow}$ satisfies $|\mathsf{X}_{n}|=n$ and is a
sequence of the form $\mathbb{N}\rightarrow\mathbf{Ahlf}_{C,\delta}$ for some
$0<C<\infty,1\leqslant\delta\leqslant\infty$.

Condition 1 defines a spatial asymptotic frame. Sequences satisfying Condition
1 will also be refered to as uniformly Ahlfors sequences. The utility of this
notion is that it gives almost enough structure to allow analysis of
$k$-medoids techniques analytically as well as derive dependent central limit
theorems and laws of large numbers. The condition is also parsimonious and
simple to express. Examples of metric spaces which belong to
$\mathbf{Ahlf}_{C,\delta}$ for a common $C,\delta$ include the sequence with
elements $\mathbb{Z}^{m}\cap\mathsf{Sq}_{n}$ where $\mathsf{Sq}_{n}$ is the
$m$-dimensional cube of side length $2n$ centered at the origin and the
sequence of annuli $\mathbb{Z}^{m}\cap\mathsf{A}_{n}$ where $\mathsf{A}_{n}$
is the $m$-dimensional annulus of outer radius length $2n$ and inner radius
length $n$ centered at the origin. An metric spaces which belong to
$\mathbf{Ahlf}_{C,\delta}$ for a common $C,\delta$ for which no constituent
metric space can be embedded into Euclidean space is the product
$\\{1,2,3,...,n\\}\times\mathsf{W}$ where $(\mathsf{W},d_{\mathsf{W}})$ is any
fixed finite metric space which cannot be embedded into Euclidean space, and
the metric on the product is given by the sum of $d_{\mathbb{Z}}$ applied to
the first component and $d_{\mathsf{W}}$ applied to the second component.

###  5.3 Asymptotics B: Mixing conditions and central limit theory

This section develops a central limit theorem for arrays of random variables
on uniformly Ahlfors sequences of metric spaces.

Let $\mathscr{D}_{\shortrightarrow}=(\mathscr{D}_{n})_{n\in\mathbb{N}}$ and
$\mathscr{D}_{n}=\\{\zeta_{i}\\}_{i\in\mathsf{X}_{n}}$ be an array of real
scalar random variables on a common probability space
$(\Omega,\mathscr{F},{\mathrm{P}})$ taking values in some measurable space
$(V,\mathscr{V})$ with spatial indices given by $\mathsf{X}_{n}$. In this
paper, $V$ will always be a finite dimensional real vector space
$\mathbb{R}^{p}$ for some $p$. Let $\mathscr{A}$ and $\mathscr{B}$ be two
(sub-)$\sigma$-algebras of $\mathscr{F}$. Define the following mixing
coefficients.

$\displaystyle\alpha^{\mathrm{mix}}(\mathscr{A},\mathscr{B})=\sup\\{|{\mathrm{P}}(A\cap
B)-{\mathrm{P}}(A){\mathrm{P}}(B)|:A\in\mathscr{A},B\in\mathscr{B}\\}.$ (5.1)

For $\mathsf{U},\mathsf{V}\subseteq\mathsf{X}_{n}$ and integers $g,l,r,$
define additionally

$\displaystyle\begin{split}\alpha_{n}^{\mathrm{mix}}(\mathsf{U},\mathsf{V})&=\alpha^{\mathrm{mix}}(\sigma(\zeta_{i}:i\in\mathsf{U}),\sigma(\zeta_{i}:i\in\mathsf{V})),\\\
\alpha_{g,l,n}^{\mathrm{mix}}(r)&=\sup_{\mathsf{U},\mathsf{V}\subseteq\mathsf{X}_{n}}\left\\{\alpha^{\mathrm{mix}}_{n}(\mathsf{U},\mathsf{V}),|\mathsf{U}|\leqslant
g,|\mathsf{V}|\leqslant l,\mathrm{d}(\mathsf{U},\mathsf{V})\geqslant
r\right\\},\\\
\bar{\alpha}_{g,l}^{\mathrm{mix}}(r)&=\sup_{n\in\mathbb{N}}\alpha^{\mathrm{mix}}_{g,l,n}(r).\end{split}$
(5.2)

The above definitions are the same definition as used in Jenish and Prucha,
(2009). The next condition, Condition 2, relies on the above definition as
well as the following fact. Let
$(\mathsf{X}_{n},\mathrm{d}_{n})\in\mathbf{Ahlf}_{C,\delta}$ for every $n$.
Then each $(\mathsf{X}_{n},\mathrm{d}_{n}^{3/4})$ has an $L$-bi-Lipschitz map
into $\mathbb{R}^{\nu}$ where $\nu$ and the Lipschitz constant $L$ depend only
on $C,\delta$.

######  Condition  2 (Mixing for an Array of Real Scalar Random Variables).

The array $\mathscr{D}_{\shortrightarrow}$ is an array of real random vectors
$\zeta_{i}$ on a common probability space $(\Omega,\mathscr{F},{\mathrm{P}})$
taking values in a common real space $\mathbb{R}^{p}$, with spatial indices
given by $\mathsf{X}_{\shortrightarrow}$ satisfying Condition 1 with Ahlfors
constants $C<\infty,\delta\geqslant 1$. Let $\nu$ correspond to $C,\delta$ as
defined in the preceding paragraph above. For any
$\mathsf{C}\subseteq\mathsf{X}$, where $\mathsf{X}=\mathsf{X}_{n}$ for some
$n$, let
$\sigma^{2}(\mathsf{C})=\mathrm{var}(\sum_{i\in\mathsf{C}}\zeta_{{i}})$. There
exists an array of positive constants
$\\{\\{c_{i}\\}_{i\in\mathsf{X}_{n}}\\}_{n\in\mathbb{N}}$ and a positive
$\mu>0$ such that,

  1. (i)

${\mathrm{E}}[\zeta_{i}]=0$.

  2. (ii)

$\lim_{k\rightarrow\infty}\sup_{n\in\mathbb{N}}\sup_{i\in\mathsf{X}_{n}}{\mathrm{E}}[\|\zeta_{i}/c_{i}\|_{2}^{2+\mu}\textbf{\em
1}_{\|\zeta_{i}/c_{i}\|_{2}>k}]=0.{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\Big{|}}$

  3. (iii)

$\sum_{m=1}^{\infty}\bar{\alpha}^{\mathrm{mix}}_{1,1}(m)m^{\nu\times\frac{\mu+2}{\mu}-1}<\infty.{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\Big{|}}$

  4. (iv)

$\sum_{m=1}^{\infty}m^{\nu-1}\bar{\alpha}^{\mathrm{mix}}_{g,l}(m)<\infty\ for\
g+l\leqslant
4.{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\Big{|}}$

  5. (v)

$\bar{\alpha}^{\mathrm{mix}}_{1,\infty}(m)=O(m^{-\nu-\frac{4}{3}\mu}).{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\Big{|}}$

  6. (vi)

$\inf_{n\in\mathbb{N}}\inf_{\mathsf{C}\subseteq\mathsf{X}_{n}}|\mathsf{C}|^{-1}(\max_{i\in\mathsf{C}}c_{i}^{-2})\lambda_{\min}\sigma(\mathsf{C})^{2}>0.{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\Big{|}}$

The conditions are similar to those given for the mixing central limit theorem
in Jenish and Prucha, (2009) which requires that each $\mathsf{X}_{n}$ be a
possibly uneven lattice in a finite dimensional Euclidean space with a minimum
separation between all points. Note that Condition 1 automatically implies a
minimum separation between points. Note also that there are only two
differences between Condition 2 and the conditions for Corollary 1 of Jenish
and Prucha, (2009) in the case that $\mathsf{X}_{n}$ embeds isometrically into
$\mathbb{R}^{\nu}$ without the need to apply the snowflake regularization
construction. First, here
$\bar{\alpha}^{\mathrm{mix}}_{1,\infty}(m)=O(m^{-\nu-\frac{4}{3}\mu})$ rather
than the slightly weaker
$\bar{\alpha}^{\mathrm{mix}}_{1,\infty}(m)=O(m^{-\nu-\mu})$. Second, Condition
2(vi) now entails an infimum over $n$ and over all subsets
$\mathsf{C}\subseteq\mathsf{X}_{n}$ while Jenish and Prucha, (2009) only
requires the condition hold for an infimum over $n$ with the particular choice
$\mathsf{C}=\mathsf{X}_{n}$. In terms of notation, here, $\zeta_{i}$ is used
rather than the more explicit $\zeta_{i,n}$, because formally the indexes $i$
belong to distinct sets, $\mathsf{X}_{n}$, for each $n$. As a result,
Condition 2 is a condition on arrays of random variables.

In the next proposition, let $\Phi$ be the Gaussian cumulative distribution
function. Also, as in the statement of Condition 2, for
$\mathsf{C}\subseteq\mathsf{X}$ for any $\mathsf{X}=\mathsf{X}_{n}$ for some
$n$, let $\sigma(\mathsf{C})=\mathrm{var}(\sum_{i\in\mathsf{C}}\zeta_{i})$

######  Theorem 1.

Suppose that $\mathsf{X}_{\shortrightarrow}$ satisfies Condition 1 and
$\mathscr{D}_{\shortrightarrow}$ satisfies Condition 2 with scalar
$\zeta_{i}$. Then for every $z\in\mathbb{R}$,
$\lim_{n\rightarrow\infty}{\mathrm{P}}(\sigma(\mathsf{X}_{n})^{-1}\sum_{i\in\mathsf{X}_{n}}\zeta_{i}\leqslant
z)=\Phi(z)$.

Theorem 1 proves convergence in distribution to a Gaussian measure for
spatially indexed arrays
$\mathscr{D}_{n}=\\{\zeta_{i}\\}_{i\in\mathsf{X}_{n}}$ which satisfy
Conditions 1 and 2. The proof is carried out in detail in the appendix. The
proposition can be proven by regularizing $\mathsf{X}_{n}$ using the snowflake
construction described above with $\varepsilon=1/4$, on which the Jenish and
Prucha, (2009) central limit theorem can then be applied. The same argument
can be carried out more generally for other $\varepsilon\in(0,1/2)$, however
the single $\varepsilon=1/4$ snowflake is referenced for added concreteness.

###  5.4 Asymptotics C: Balance and small common boundary conditions

A key observation in Bester et al., (2011) is that if a fixed sequence of
partitions of $\mathsf{X}_{n}$ have balanced cluster sizes and small
boundaries, then under mixing and additional regularity conditions, the CCE
procedure as described in Bester et al., (2011) achieves asymptotically
correct size. We show that this phenomenon also holds in the settings
considered in this paper. This section formalizes notions of balanced cluster
sizes and small boundaries appropriate for sequences
$\mathsf{X}_{\shortrightarrow}$ satisfying Condition 1. This section then
gives a proposition which shows that under the appropriate regularity
conditions, the endpoint of the $k$-medoids algorithm satisfies these two
desired properties. Let
$\mathcal{C}_{\shortrightarrow}=(\mathcal{C}_{n})_{n\in\mathbb{N}}$ be a
sequence such that $\mathcal{C}_{n}$ is a partition of $\mathsf{X}_{n}$ for
each $n$.

######  Condition  3 (Balance and Small Boundaries).

The sequence of partitions $\mathcal{C}_{\shortrightarrow}$ satisfies
$|\mathcal{C}_{n}|\leqslant k_{\max}$ for some $k_{\max}$ which does not
depend on $n$, and is asymptotically balanced with small boundaries in the
sense of the following conditions.

  1. (i)

The cluster sizes satisfy $\underset{n\rightarrow\infty}{\lim\inf}\hskip
2.84526pt\frac{\hskip
2.84526pt\min_{\mathsf{C}\in\mathcal{C}_{n}}{|\mathsf{C}|}\hskip
2.84526pt}{{\max_{\mathsf{D}\in\mathcal{C}_{n}}{|\mathsf{D}|}}}>0,$

  2. (ii)

There is a sequence
$\bar{r}_{\shortrightarrow}=(\bar{r}_{n})_{n\in\mathbb{N}}\nearrow\infty$ with
$\lim_{n\rightarrow\infty}\frac{\max_{\mathsf{C}\in\mathcal{C}_{n}}|\\{i\in\mathsf{C}:\mathrm{d}(i,\mathsf{X}_{n}\setminus\mathsf{C})\leqslant\bar{r}_{n}\\}|}{\min_{\mathsf{D}\in\mathcal{C}_{n}}|\mathsf{D}|}=0.$

The definition of small boundary provided in Condition 3(ii) differs slightly
from the definition given in Bester et al., (2011). In particular, Bester et
al., (2011) leverage the fact their spatial domain is a subset of the integer
lattice to define neighbor orders for pairs of locations. Their definition of
small boundaries entails a bound on the number of first order neighbors from
$\mathsf{C}$ to $\mathsf{X}_{n}\setminus\mathsf{C}$. They assume that their
given spatial clusters are contiguous and use that fact to bound the number of
higher order neighbors from $\mathsf{C}$ to
$\mathsf{X}_{n}\setminus\mathsf{C}$. In this context, there is no available
definition of first order neighbor since $\mathsf{X}_{n}$ can be irregular
(even non-Euclidean). As a result, this paper works instead with an asymptotic
notion of boundary which entails a sequence $\bar{r}_{\shortrightarrow}$ which
allows boundaries to widen as $n\rightarrow\infty.$ An implication of having
asymptotically balanced with small boundaries clusterings is the following
proposition.

######  Theorem 2.

Suppose that $\mathsf{X}_{\shortrightarrow}$ satisfies Condition 1,
$\mathscr{D}_{\shortrightarrow}$ satisfies Condition 2, and
$\mathcal{C}_{\shortrightarrow}$ satisfies Condition 3. Then
$\lim_{n\rightarrow\infty}\sup_{\mathsf{C}\neq\mathsf{D}\in\mathcal{C}_{n}}\mathrm{cov}_{{\mathrm{P}}}\left(\sigma({\mathsf{C}})^{-1}\sum_{i\in\mathsf{C}}\zeta_{i},\
\sigma({\mathsf{D}})^{-1}\sum_{i\in\mathsf{D}}\zeta_{i}\right)=0.$

The argument establishing the above fact is related to but not identical to
arguments in Bester et al., (2011), which were previously also given in Jenish
and Prucha, (2009) and Bolthausen, (1982). Instead of counting points in
“shells” around the boundaries of clusters, the proof of Theorem 2 instead
relies on the doubling structure implied by the fact that $\mathsf{X}_{n}$ are
Ahlfors regular. Both arguments leverage a bound on covariances
$\text{cov}(\zeta_{i},\zeta_{j})$ for sufficiently distant spatial locations
as implied by the mixing conditions stated in Condition 2.

Condition 3 is a taken as a high-level assumption on clustering structures in
Theorem 2. Next, we verify that clusters produced by $k$-medoids satisfy
Condition 3 under an additional convexity assumption. $k$-medoids is a popular
clustering technique which is related to $k$-means, and both produce similar
clustering results in many settings. However, $k$-means centroids are not
necessarily defined in the space $\mathsf{X}_{n}$ when dealing with a
dissimilarity which does not necessarily arise from a Euclidean space.
$k$-medoids differs by requiring that clusters are defined around elements of
$\mathsf{X}_{n}$, called medoids, which must themselves also be elements of
$\mathsf{X}_{n}$. For a detailed description of $k$-medoids see Hastie et al.,
(2009).

The next proposition derives relevant asymptotic properties for a $k$-medoids
algorithm by studying the sequence of partitions

$\displaystyle\mathcal{C}_{\shortrightarrow}=(\mathcal{C}_{n})_{n\in\mathbb{N}}\text{
defined by }\mathcal{C}_{n}=k\text{-}\mathtt{medoids}(\mathsf{X}_{n}).$ (5.3)

######  Theorem 3.

Suppose that $\mathsf{X}_{\shortrightarrow}$ satisfies Condition 1. Assume the
following additional convexity condition holds. There is a constant $K$
independent of $n$ such that for each $n$, (1) $\mathsf{X}_{n}$ is
$K$-coursely isometric666
$f:(\mathsf{Y},\mathrm{d}_{\mathsf{Y}})\rightarrow(\mathsf{Z},\mathrm{d}_{\mathsf{Z}})$
is a $K$-course isometry if
$\mathrm{d}_{\mathsf{Z}}(f(i),f(j))-K\leqslant\mathrm{d}_{\mathsf{Y}}(i,j)\leqslant\mathrm{d}_{\mathsf{Z}}(f(i),f(j))+K$.
to a subset of a Euclidean space with dimension $\nu$ independent of $n$, and
(2) for any two point $i,j\in\mathsf{X}_{n}$ and any $a\in[0,1]$ there is an
interpolant $g\in\mathsf{X}_{n}$ such that
$|\mathrm{d}_{n}(i,g)-a\mathrm{d}_{n}(j,g)|\leqslant K$ and
$|\mathrm{d}_{n}(i,g)-(1-a)\mathrm{d}_{n}(j,g)|\leqslant K$. Let
$\mathcal{C}_{n}=k\text{-}\mathtt{medoids}(\mathsf{X}_{n})$ and let
$\mathcal{C}_{\shortrightarrow}=(\mathcal{C}_{n})_{n\in\mathbb{N}}$. Then
$\mathcal{C}_{\shortrightarrow}$ satisfies Condition 3.

Theorem 3 does not require that the sequence of partitions $\mathcal{C}_{n}$
of $\mathsf{X}_{n}$ resulting from applying $k$-medoids to converge in any
way. It also does not require that there be any notion of true clusters or a
true partition associated to any of the $\mathsf{X}_{n}$.

The convexity assumption is needed to exclude cases such as the following
would-be counterexample:

$\displaystyle\begin{split}\mathsf{Rothko}_{n}&=\\{(i^{1},i^{2})\in\mathbb{Z}^{2}:|i^{1}|\leqslant
n,|i^{2}|\leqslant 2n\\}\\\ &\ \ \ \ \ \ \ \
\setminus\\{(i^{1},i^{2})\in\mathbb{Z}^{2}:|i^{1}|\leqslant
n-2,1\leqslant|i^{2}+\frac{1}{2}|\leqslant 2n-2\\}.\end{split}$ (5.4)

The naming for the above example is based on the fact that when plotted,
$\mathsf{Rothko}_{n}$ resembles the background orange region in the oil
painting Rothko, (1956), titled Orange and Yellow. A possible sequence of
endpoints of the $2$-medoids algorithm for the above example is with medoids
at $(0,0)$ and $(0,-1)$ and boundary (with $\bar{r}=1$) at
$\\{0,1\\}\times\\{-n,-n+1,...,n\\}$. Note that these medoid endpoints are not
unique, and other medoid sequences, like $(-n,2n),(n,-2n)$ do lead to small
boundaries.

###  5.5 Asymptotics D: Almost sure representation

The next preliminary for the analysis of Algorithm 1 is an almost sure
representation theorem. Suppose that $\mathscr{C}_{\shortrightarrow}$ is a
sequence of collections of partitions of $\mathsf{X}_{n}$. That is
$\mathscr{C}_{\rightarrow}=(\mathscr{C}_{n})_{n\in\mathbb{N}}$ and each
$\mathscr{C}_{n}$ is a set containing elements $\mathcal{C}$ which are
partitions of $\mathsf{X}_{n}$. The following condition concerns corresponding
collections $S_{\mathcal{C}}$ of random vectors taking values in
$\mathbb{R}^{\mathcal{C}}$ where $\mathcal{C}\in\mathscr{C}_{n}$ for some $n$.
Let $\\{S_{\mathcal{C}}\\}$ be a collection of random vectors in
$\mathbb{R}^{\mathcal{C}}$ where $\mathcal{C}$ range in
$\cup_{n\in\mathbb{N}}\mathscr{C}_{n}$.

######  Condition  4 (Almost Sure Asymptotic Gaussian Representation).

Associated to all $S_{\mathcal{C}}\in\mathbb{R}^{\mathcal{C}}$ for all
$\mathcal{C}\in\cup_{n\in\mathbb{N}}\mathscr{C}_{n}$ is a common probability
space $(\tilde{\Omega},\tilde{\mathscr{F}},\tilde{\mathrm{P}})$ with random
variables
$\tilde{S}_{\mathcal{C}},\tilde{S}_{\mathcal{C}}^{*}\in\mathbb{R}^{\mathcal{C}}$.
Furthermore, $\tilde{S}_{\mathcal{C}}$ and ${S}_{\mathcal{C}}$ have the same
distribution, $\tilde{S}_{\mathcal{C}}^{*}$ is zero-mean Gaussian with
independent components,
$\liminf_{n\rightarrow\infty}\min_{\mathsf{C}\in\mathcal{C}\in\mathscr{C}_{n}}\mathrm{var}(\tilde{S}_{\mathcal{C},\mathsf{C}}^{*})>0$,
$\limsup_{n\rightarrow\infty}\max_{\mathsf{C}\in\mathcal{C}\in\mathscr{C}_{n}}\mathrm{var}(\tilde{S}_{\mathcal{C}_{n},\mathsf{C}}^{*})<\infty$
and
$\lim_{n\rightarrow\infty}\sup_{\mathcal{C}\in\mathscr{C}_{n}}\|\tilde{S}_{\mathcal{C}}-\tilde{S}_{\mathcal{C}}^{*}\|_{2}=0$
$\tilde{\mathrm{P}}$-almost surely.

If the above condition holds, then $S_{\mathcal{C}}$ are said to be uniformly
asymptotically almost surely representable by non degenerate independent
Gaussians. When needed, the representors relative to $S_{\mathcal{C}}$ will
always be referred to by $\tilde{S}_{\mathcal{C}}^{*}$ and
$\tilde{S}_{\mathcal{C}}$ just as in the definition. It is a well known result
that convergence in distribution of a sequence of random variables implies
almost sure convergence for an auxiliary sequence of random variables with the
same distribution as the original sequence. Specifically, Theorem 2.19 in van
der Vaart, (1998) states: “Suppose that the sequence of random vectors $X_{n}$
converges in distribution to a random vector $X_{0}$. Then there exists a
probability space $(\tilde{\Omega},\tilde{\mathcal{U}},\tilde{P})$ and random
vectors $\tilde{X}_{n}$ defined on it such that $\tilde{X}_{n}$ is equal in
distribution to $X_{n}$ for every $n\geqslant 0$ and
$\tilde{X}_{n}\rightarrow\tilde{X}_{0}$ almost surely.”

A subtle problem that occurs in the current setting is that the random vectors
$S_{\mathcal{C}}$ take values in $\mathbb{R}^{\mathcal{C}}$, which are
technically different spaces for different $\mathcal{C}$. One collection of
$S_{\mathcal{C}}$ of interest are constructed from
$\mathscr{D}_{\shortrightarrow}$ and a sequence of partitions
$\mathcal{C}_{\shortrightarrow}$ of $\mathsf{X}_{\shortrightarrow}$ by

$\displaystyle
S_{\mathcal{C},\mathsf{C}}=(n/k)^{-1/2}\sum_{i\in\mathsf{C}}\zeta_{i},\text{
for }\mathcal{C}\text{ partitioning }\mathsf{X}_{n}\text{ and
}|\mathcal{C}|=k.$ (5.5)

The issue of $S_{\mathcal{C}}$ belonging to different spaces does not yet
exist in Theorem 2 because the covariance is relative to the measure
${\mathrm{P}}$, which is common for all $n$. The issue is also circumvented in
Theorem 1 because there, the terms
$\sigma(\mathsf{X}_{n})^{-1}\sum_{i\in\mathsf{X}_{n}}\zeta_{i}$ may be taken
to belong to a common set, namely $\mathbb{R}$, as each $\zeta_{i}$ belongs to
$\mathbb{R}$. Nevertheless, the following holds.

######  Theorem 4.

Suppose that $\mathsf{X}_{\rightarrow}$ satisfies Condition 1,
$\mathscr{D}_{\rightarrow}$ satisfies Condition 2, and
$\mathcal{C}_{\rightarrow}$ satisfies Condition 3 with the additional
condition that $|\mathcal{C}_{n}|\leqslant{k_{\mathrm{max}}}$ with
${k_{\mathrm{max}}}$ independent of $n$. For $\mathcal{C}$ with
$|\mathcal{C}|=k$, let
$S_{\mathcal{C},\mathsf{C}}=(n/k)^{-1/2}\sum_{i\in\mathsf{C}}\zeta_{i}$. Let
$\mathscr{C}_{n}$ be the singleton $\mathscr{C}_{n}=\\{\mathcal{C}_{n}\\}$.
Then
$\\{S_{\mathcal{C}}\\}_{\mathcal{C}\in\cup_{n\in\mathbb{N}}\mathscr{C}_{n}}$
satisfies Condition 4.

###  5.6 Asymptotics E: Central Limit theory in the cases of OLS and IV

The previous theorems establish asymptotic normality and diminishing across-
group correlation for means. This section gives a similar limiting result in
the context of ordinary least squares (OLS) and instrumental variables (IV)
estimators. These are the precise central limit theorems which are relevant to
the subsequent analysis of cluster-based inferential procedures.

Consider data $\mathscr{D}_{n}=\\{\zeta_{i}\\}_{i\in\mathsf{X}_{n}}$ defined
on $(\Omega,\mathscr{F},{\mathrm{P}})$ with
$\zeta_{i}=(Y_{i},X_{i},W_{i},Z_{i})$ in which for each $i\in\mathsf{X}_{n}$
(and for each $n$) a linear model holds as follows:

$\displaystyle Y_{i}=\theta_{0}X_{i}+W_{i}^{\prime}\gamma_{0}+U_{i},$ (5.6)

where $Y_{i}$ is a scalar outcome variable, $X_{i}$ is a scalar variable of
interest, $U_{i}$ is a scalar idiosyncratic disturbance term, $W_{i}$ is a
$p$-dimensional vector of control variables, and $Z_{i}$ is a
$(p+1)$-dimensional vector of instruments. The following regularity condition
is helpful.

######  Condition  5.

$\mathscr{D}_{\shortrightarrow}$ consists of
$\mathscr{D}_{n}=\\{(Y_{i},X_{i},W_{i},Z_{i})\\}_{i\in\mathsf{X}_{n}}$
satisfying the linear model
$Y_{i}=\theta_{0}X_{i}+W_{i}^{\prime}\gamma_{0}+U_{i}$ for some $\theta_{0}$
and $\gamma_{0}$. ${\mathrm{E}}[Z_{i}U_{i}]=0$,
${\mathrm{E}}[Z_{i}X_{i}^{\prime}]$ exists and equals some $\mathrm{M}$ which
does not depend on $i$ and which has positive eigenvalues.
$(\\{Z_{i}U_{i},Z_{i}X_{i}^{\prime}-{\mathrm{E}}[Z_{i}X_{i}^{\prime}]\\}_{i\in\mathsf{X}_{n}})_{n\in\mathbb{N}}$
satisfies Condition 2.

Consider the IV estimator for $\theta_{0}$ using only observations
$\mathsf{C}\subseteq\mathsf{X}_{n}$, and, as earlier, denote this by
$\widehat{\theta}_{\mathsf{C}}$. Note, the OLS estimate is a special case
where $Z_{i}=X_{i}$. Let $\mathcal{C}$ be a partition of $\mathsf{X}_{n}$ for
some $n$ with $k$ clusters. Define
$S_{\mathcal{C}}\in\mathbb{R}^{\mathcal{C}}$ by

$\displaystyle S_{\mathcal{C},\mathsf{C}}=(n/k)^{1/2}\hskip
1.42262pt(\widehat{\theta}_{\mathsf{C}}-\theta_{0})\ \text{for}\
\mathsf{C}\in\mathcal{C}.$ (5.7)

Let $\mathscr{C}_{n}$ again be the singleton
$\mathscr{C}_{n}=\\{\mathcal{C}_{n}\\}$. The next theorem is the particular
asymptotic normality result about $\widehat{\theta}_{\mathsf{C}}$,
$\mathsf{C}\in\mathcal{C}_{n}$ which is useful for input into the analysis of
Algorithm 1 in the context of testing hypotheses about coefficients in a
linear model.

######  Theorem 5.

Suppose that $\mathsf{X}_{\shortrightarrow}$ satisfies Condition 1. Suppose
linear regression data $\mathscr{D}_{\shortrightarrow}$ satisfies Condition 5.
Suppose that $\mathcal{C}_{\shortrightarrow}$ satisfies Condition 3. Then for
$S_{\mathcal{C}}$ defined above in this subsection and for
$\mathscr{C}_{n}=\\{\mathcal{C}_{n}\\}$, the collection
$\\{S_{\mathcal{C}}\\}_{\mathcal{C}\in\cup_{n\in\mathbb{N}}\mathscr{C}_{n}}$
satisfies Condition 4.

###  5.7 Analysis of Cluster-Based Inference with Learned Clusters

Consider an array $\mathscr{D}_{\shortrightarrow}$ of spatially indexed
datasets $\mathscr{D}_{n}=\\{\zeta_{i}\\}_{i\in\mathsf{X}_{n}}$ all
distributed according to one common ${\mathrm{P}}_{0}$. Consider a sequence of
statistical hypotheses $H_{0\shortrightarrow}=\\{H_{0n}\\}_{n\in\mathbb{N}}$
where each $H_{0n}$ is a restriction on only the distribution of the $n$th
sample of $\mathscr{D}_{\shortrightarrow}$ (i.e., on $\mathscr{D}_{n}$). The
coming sections discuss formal properties of sequences of tests
$\widehat{T}_{\shortrightarrow}=(\widehat{T}_{n})_{n\in\mathbb{N}}$ with
primary interest in
$\widehat{T}_{\shortrightarrow}=\widehat{T}_{\bullet(\alpha)\shortrightarrow}$
arising from $\widehat{T}_{\bullet(\alpha)}$ defined by applying Algorithm 1
to data $\mathscr{D}_{n},\mathsf{X}_{n}$ with $\bullet=$ IM, CRS, or CCE as in
Section 2.

For any $n$, let $\mathscr{C}_{n}$ be its associated collection of partitions

$\displaystyle\mathscr{C}_{n}=\\{\mathcal{C}:\mathcal{C}=\mathtt{Cluster}(\mathsf{X}_{n})\text{
and }(\mathtt{Test},\mathtt{Cluster})\in\mathscr{T}_{n}\text{ for some
}\mathtt{Test}\\}.$ (5.8)

For each $n$, there is an element
$\widehat{\texttt{Cluster}}(\mathsf{X}_{n})\in\mathscr{C}_{n}$ which is
defined by the solution to the optimization problem in equation (2.9).

An important property used in the analysis of Algorithm 1 in Sections 5.8 and
5.9 is that the act of selecting partitions $\widehat{\mathcal{C}}_{n}$ has
negligible effect on $\widehat{T}_{n}$. This property is formalized in the
next condition.

######  Condition  6.

Corresponding to $\widehat{T}_{\bullet(\alpha)\shortrightarrow}$, there are
sets $\mathscr{C}_{1,n}\subseteq\mathscr{C}_{n}$ for all $n$ such that the
following are satisfied:

$\displaystyle(i)$ $\displaystyle\
\lim_{n\rightarrow\infty}{{\mathrm{P}}_{0}}(\widehat{\mathtt{Cluster}}(\mathsf{X}_{n})\notin\mathscr{C}_{1,n})=0,$
$\displaystyle(ii)$ $\displaystyle\
\lim_{n\rightarrow\infty}\sup_{\mathcal{C}\in\mathscr{C}_{1,n}}\Big{|}{{\mathrm{P}}_{0}}(T_{\bullet(\alpha),\mathcal{C}}=\mathrm{Reject}|\widehat{\mathtt{Cluster}}(\mathsf{X}_{n})=\mathcal{C})-{{\mathrm{P}}_{0}}(T_{\bullet(\alpha),\mathcal{C}}=\mathrm{Reject})\Big{|}=0.$

When specializing to the context of Algorithm 1 in Section 2, Condition 6 is
equivalently a high-level condition on the properties
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-}\mathrm{I}}$ and
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-}\mathrm{II},\mathbf{P}_{\mathrm{alt}}}$.
Recall that the method for selecting $\widehat{\mathcal{C}}$ described in
Section 2 is through a combination of QMLE, simulation, and constrained
minimization of a loss function (defined as estimated weighted average power).
Under sufficient regularity conditions, the QMLE procedure defined in the
appendix is consistent for a fixed object which may depend on $n$. If the loss
function is sufficiently regular, than $\mathscr{C}_{1,n}$ satisfying
Condition 6 can be taken to be a singleton. Alternatively, if a single fixed
$k$ is used in conjunction with $k$-medoids clustering, Condition 6 is
automatically satisfied.

###  5.8 General analysis for IM with learned clusters

The IM procedure at level $a$ was defined in Section 2 for a scalar hypothesis
about a regression coefficient in a linear model of the form
$H_{0}:\theta_{0}=\theta^{\circledast}$, depending on a partition
$\mathcal{C}$ of the data’s spatial indexing set. In particular, recall that
in Section 2, the IM test was defined for $a\in(0,1)$ and given by
$T_{\text{IM}(a),\mathcal{C}}=\text{Reject}\ \ \ \ \text{if}\ \ \
|t(S)|>t_{1-a/2,k-1},$ where $t_{1-a/2,k-1}$ is the $(1-a/2)$-quantile of
$t$-distribution with $k-1$ degrees of freedom and with
$S=(S_{\mathsf{C}})_{\mathsf{C}\in\mathcal{C}}\in\mathbb{R}^{\mathcal{C}}$,
$|\mathcal{C}|=k$,
$S_{\mathsf{C}}=(n/k)^{1/2}(\widehat{\theta}_{\mathsf{C}}-\theta^{\circledast})$.

######  Theorem 6.

Suppose that
$\\{S_{\mathcal{C}}\\}_{\mathcal{C}\in\cup_{n\in\mathbb{N}}\mathscr{C}_{n}}$
satisfy Condition 4 under ${\mathrm{P}}_{0}$. Suppose
$\widehat{T}_{\mathrm{IM}(\alpha)\shortrightarrow}$ satisfies Condition 6.
Suppose further that $\alpha\leq 2\Phi(-\sqrt{3})=0.083...$, where $\Phi$ the
standard Gaussian cumulative distribution function. Then
$\lim\sup_{n\rightarrow\infty}{{\mathrm{P}}_{0}}(\widehat{T}_{\mathrm{IM}(\alpha),n}=\text{\em
Reject})\leqslant\alpha.$

The condition that $\alpha\leq 2\Phi(-\sqrt{3})$ is also needed in Ibragimov
and Müller, (2010). Note that the conditions above were verified for the OLS
and IV models under regularity conditions given in the previous section with
partitions generated by $k$-medoids with fixed $k$. The properties of
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-I}}$ and
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-I},\mathbf{P}_{\mathrm{alt}}}$ do
not appear explicitly in the conditions of Theorem 6. However, their
properties are reflected through Condition 6. The proof of Theorem 6 builds
largely on the arguments in Ibragimov and Müller, (2010). Relative to
Ibragimov and Müller, (2010), Theorem 6 requires establishing uniform limit
results over tests corresponding to partitions in
$\mathscr{C}_{\shortrightarrow}.$ Additional details can be found in the
discussions provided within the proof in the supplemental appendix.

Theorem 6 is helpful in understanding the statistical performance of Algorithm
1 with IM in the context of the hypothesis tests considered in the empirical
application in Section 3 and simulation study in Section 4. In particular,
consider a sequence of regression datasets $\mathscr{D}_{\shortrightarrow}$
with $\mathscr{D}_{n}=\\{\zeta_{i}\\}_{i\in\mathsf{X}_{n}}$ with
$\zeta_{i}=(Y_{i},X_{i},W_{i},Z_{i})$ as in Section 5.6 satisfying
$Y_{i}=\theta_{0}X_{i}+W_{i}^{\prime}\gamma_{0}+U_{i}$.

######  Corollary 1.

($\mathrm{IM}$ and $k$-medoids in the linear model). Suppose that
$\mathsf{X}_{\shortrightarrow}$ satisfies Condition 1. Suppose that
$\mathscr{D}_{\shortrightarrow}$ satisfies Condition 5. Let
$H_{0\shortrightarrow}=(H_{0,n})_{n\in\mathbb{N}}$ be the sequence of
hypotheses $H_{0,n}:\theta_{0}=\theta^{\circledast}$. Let
$\widehat{T}_{\mathrm{IM}(\alpha)\shortrightarrow}=(\widehat{T}_{\mathrm{IM}(\alpha),n})_{n\in\mathbb{N}}$
be defined by Algorithm 1 applied to $(\mathscr{D}_{n},\mathsf{X}_{n})$ for
hypotheses $H_{0,n}$ for each $n$ with $k$-medoids and $k_{\max}$ fixed
independent of $n$. Suppose
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-}\mathrm{I}}$ and
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-}\mathrm{II},\mathbf{P}_{\mathrm{alt}}}$
are such that $\widehat{T}_{\mathrm{IM}(\alpha)\shortrightarrow}$ satisfies
Condition 6. Then
$\lim_{n\rightarrow\infty}{\mathrm{P}}_{0}(\widehat{T}_{\mathrm{IM}(\alpha),n}=\
\mathrm{Reject}\ )\leqslant\alpha$.

###  5.9 General analysis for CRS with learned clusters

The procedure defining $T_{\text{CRS}(\alpha)}$ can be stated in slightly more
generality. In this more general setup, CRS relative to a partition
$\mathcal{C}$ depends on a real-valued function
$w:\mathbb{R}^{\mathcal{C}}\rightarrow\mathbb{R}$, and a statistic
$S\in\mathbb{R}^{\mathcal{C}}$, which is chosen such that large values of $w$
provide evidence against $H_{0}$. As was described in Section 2, in all
implementations in this paper, $w$ is chosen to be $w(S)=|t({S})|$ and $S$ has
components
$S_{\mathsf{C}}=(n/k)^{1/2}(\widehat{\theta}_{\mathsf{C}}-\theta^{\circledast})$
with $k=|\mathcal{C}|$. Note, the CRS test was developed in Canay et al.,
(2017) for a single fixed partition and relies on a notion of approximate
symmetry. Specifically, given a partition $\mathcal{C}$ of $\mathsf{X}_{n}$,
the CRS procedure also depends on the availability of a finite group
$\mathcal{H}_{\mathcal{C}}$ of symmetries which act within component
$\mathsf{C}$ of $\mathcal{C}$. These symmetries are formalized by a group
action $\mathscr{D}_{n}\mapsto h\mathscr{D}_{n}$ for all
$h\in\mathcal{H}_{\mathcal{C}}$. The action on the data is defined in such a
way that

$\displaystyle\mathscr{D}_{n}\mapsto h\mathscr{D}_{n}\ \text{ induces }\
S_{\mathcal{C}}\mapsto hS_{\mathcal{C}}$ (5.9)

which respects the original action in the sense
$S_{\mathcal{C}}(h\mathscr{D}_{n})=hS_{\mathcal{C}}(\mathscr{D}_{n})$. With
slight abuse of notation, write $w(\mathscr{D}_{n})=w(S_{\mathcal{C}})$.

In the context of cluster-based inference in this paper, the only type of
group action of interest will be set of signs
$\mathcal{H}_{\mathcal{C}}=\\{-1,+1\\}^{\mathcal{C}}$ that operates within
clusters. Elements, $h$, which are tuples of signs
$h_{\mathsf{C}}\in\\{-1,+1\\}$ indexed by $\mathsf{C}\in\mathcal{C}$ can be
defined to act in a way such that $S_{\mathcal{C},\mathsf{C}}\mapsto
h_{\mathsf{C}}S_{\mathcal{C},\mathsf{C}}.$ In the context of the OLS and IV
models, $(\widehat{\theta}_{\mathsf{C}}-\theta^{\circledast})\mapsto
h_{\mathsf{C}}(\widehat{\theta}_{\mathsf{C}}-\theta^{\circledast})$

To complete the general description of CRS, let
$M=|\mathcal{H}_{\mathcal{C}}|$ and let $w^{1}(S_{\mathcal{C}})\leqslant
w^{2}(S_{\mathcal{C}})\leqslant\cdots\leqslant w^{M}(S_{\mathcal{C}})$ denote
the order statistics of the orbit
$\left\\{w\left(h\mathscr{D}_{n}\right):h\in\mathcal{H}\right\\}$. Let
$j_{a}=M(1-a)$, $M^{+}(X)=|\left\\{1\leqslant j\leqslant
M:w^{j}(S_{\mathcal{C}})>w^{{j_{a}}}S_{\mathcal{C}}\right\\}|$ and
$M^{0}(X)=|\left\\{1\leqslant j\leqslant
M:w^{j}(X)=w^{{j_{a}}}(S_{\mathcal{C}})\right\\}|$. Let
$\tilde{a}\left(S_{\mathcal{C}}\right)=\frac{Ma-M^{+}(S_{\mathcal{C}})}{M^{0}(S_{\mathcal{C}})}$.
Set $T_{\text{CRS}(a),\mathcal{C}}=\text{Reject}\ \ \ \ \text{if}\ \ \
w(S_{\mathcal{C}})>w^{j_{a}}(S_{\mathcal{C}})\ \ \ \text{or}\ \
\left(w(S_{\mathcal{C}})=w^{j_{a}}(S_{\mathcal{C}})\text{ and
}\tilde{a}(S_{\mathcal{C}})=1\right).$ $\widehat{T}_{\text{CRS}(\alpha)}$ is
then defined as in Section 2.

######  Condition  7 (Regularity Conditions for CRS).

  1. (i)

$\\{S_{\mathcal{C}}\\}_{\mathcal{C}\in\cup_{n\in\mathbb{N}}\mathscr{C}_{n}}$
satisfy Condition 4.

  2. (ii)

$h\tilde{S}^{*}_{\mathcal{C}}$ has the same distribution as
$\tilde{S}^{*}_{\mathcal{C}}$ for all $h\in\mathcal{H}_{\mathcal{C}}$, for all
$\mathcal{C}$ in all $\mathscr{C}_{n}$; for distinct $h,h^{\prime}$, either
$w\circ h=w\circ h^{\prime}$ or
$\tilde{\mathrm{P}}(w(h\tilde{S}^{*}_{\mathcal{C}})\neq
w(h^{\prime}\tilde{S}^{*}_{\mathcal{C}}))=1$; $w$ is continuous and the action
of $h$ is continuous for each $h$ in each $\mathcal{H}_{\mathcal{C}}.$

  3. (iii)

$|\mathscr{C}_{n}|$ and $\max_{\mathcal{C}\in\mathscr{C}_{n}}|\mathcal{C}|$
are each bounded by a constant independent of $n$.

######  Theorem 7.

Suppose that
$\\{S_{\mathcal{C}}\\}_{\mathcal{C}\in\cup_{n\in\mathbb{N}}\mathscr{C}_{n}}$
satisfies Condition 7 under ${\mathrm{P}}_{0}$. Suppose
$\widehat{T}_{\mathrm{CRS}(\alpha)\shortrightarrow}$ satisfies Condition 6.
Then
$\lim\sup_{n\rightarrow\infty}{{\mathrm{P}}_{0}}(\widehat{T}_{{\mathrm{CRS}(\alpha)},n}=\text{\em
Reject})\leqslant\alpha.$

In Canay et al., (2017), it was shown that the conclusion of Theorem 7 holds
under conditions which are subsumed by Condition 7, working in the case that
$\mathscr{C}_{n}$ are singletons, i.e.,
$\mathscr{C}_{n}=\\{\mathcal{C}_{n}\\}$ consist of a unique predetermined
partition of $\\{1,...,n\\}$ without any described dependence on any spatial
indexing $\mathsf{X}_{n}$.

Condition 7(iii) is stronger than necessary. Discussion about weakening
Condition 7(iii) is postponed until Subsection 5.11 about uniformity. In
particular, no requirement on $|\mathscr{C}_{n}|$ are needed when
$w(S)=|t(S)|$. The supplemental material proves Theorem 7 under a strictly
weaker condition, Condition 7(iii)’, which is stated formally there. The
argument accommodating the ability to generalize to this more general case is
new relative to Canay et al., (2017).

As was the case in previous section with the analysis of IM, Theorem 7 has an
important corollary which specializes to the specific case of Algorithm 1
(with $k$-medoids) for test $H_{0}:\theta_{0}=\theta^{\circledast}$, where
$\theta_{0}$ is an unknown coefficient in a linear regression model. Once
again, consider a sequence of regression datasets
$\mathscr{D}_{\shortrightarrow}$ with
$\mathscr{D}_{n}=\\{\zeta_{i}\\}_{i\in\mathsf{X}_{n}}$ with
$\zeta_{i}=(Y_{i},X_{i},W_{i},Z_{i})$ as in Subsection 5.6 above satisfying
$Y_{i}=\theta_{0}X_{i}+W_{i}^{\prime}\gamma_{0}+U_{i}$.

######  Corollary 2.

($\mathrm{CRS}$ and $k$-medoids in the linear model). Suppose that
$\mathsf{X}_{\shortrightarrow}$ satisfies Condition 1. Suppose that
$\mathscr{D}_{\shortrightarrow}$ satisfies Condition 5. Let
$H_{0\shortrightarrow}=(H_{0,n})_{n\in\mathbb{N}}$ be the sequence of
hypotheses $H_{0,n}:\theta_{0}=\theta^{\circledast}$. Let
$\widehat{T}_{\mathrm{CRS}(\alpha)\shortrightarrow}=(\widehat{T}_{\mathrm{CRS}(\alpha),n})_{n\in\mathbb{N}}$
be defined by Algorithm 1 applied to $(\mathscr{D}_{n},\mathsf{X}_{n})$ for
hypotheses $H_{0,n}$ for each $n$ with $k$-medoids and $k_{\max}$ fixed
independent of $n$. Suppose
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-}\mathrm{I}}$ and
$\widehat{\mathrm{Err}}_{\mathrm{Type}\text{-}\mathrm{II},\mathbf{P}_{\mathrm{alt}}}$
are such that $\widehat{T}_{\mathrm{CRS}(\alpha)\shortrightarrow}$ satisfies
Condition 6. Then
$\lim_{n\rightarrow\infty}{\mathrm{P}}_{0}(\widehat{T}_{\mathrm{IM}(\alpha),n}=\
\mathrm{Reject}\ )\leqslant\alpha$.

###  5.10 Discussion for BCH with learned clusters

Recall that the CCE test for $H_{0}:\theta_{0}=\theta^{\circledast}$ using a
partition $\mathcal{C}_{n}$ with $k$ clusters is defined by
$T_{\text{CCE}(a),n}=\ \text{Reject}\ \ \text{ if }\ \
\frac{\widehat{\theta}_{\mathsf{X}_{n}}-\theta^{\circledast}}{\widehat{V}_{\text{CCE},\mathcal{C}_{n}}^{1/2}}>\sqrt{\frac{k}{k-1}}\times
t_{1-a/2,k-1}$ where $\widehat{V}_{\text{CCE},\mathcal{C}_{n}}$ is the
standard cluster covariance estimator (without a degrees of freedom
correction). BCH analyzed the CCE estimator under an asymptotic frame with a
fixed, finite number of clusters $k$ (i.e. $k$ independent of $n$). The
resulting inference is based on calculating a $t$ statistic based on an
estimated standard error. Under regularity conditions, BCH show that the
resulting $t$ statistic is asymptotically pivotal, but distributed according
to $\sqrt{k/(k-1)}\times t_{k-1}$ where $t_{k-1}$ is the $t$ distribution with
$k-1$ degrees of freedom.

The regularity conditions required in Bester et al., (2011) are strong, and in
particular, require that the clusters have asymptotically equal numbers of
observations. The $k$-medoids algorithm does not generally return such a
partition. As a result, the CCE estimator is not anticipated to have
asymptotically correct size.

###  5.11 Discussion of uniformity

The high-level conditions used for proving Theorems 6 and 7 involve the
assumption of uniform central limit theorems for convergences
$\sup_{\mathcal{C}\in\mathscr{C}_{n}}\|\tilde{S}_{\mathcal{C}}-\tilde{S}^{*}_{\mathcal{C}}\|_{2}\rightarrow
0$ for Gaussian random variables $\tilde{S}_{\mathcal{C}}$. Uniform
convergence results may be explicitly derived in several cases of interest.
First, when each $\mathscr{C}_{n}$ contains only finitely many partitions
$\mathcal{C}$ all of which have cardinality bounded by some
${k_{\mathrm{max}}}$ independent of $n$, then uniform convergence follows
immediately from a pointwise convergence result like that presented in Section
5.6.

Another case of interest involves increasing sequences ${k_{\mathrm{max},n}}$.
In particular, if ${k_{\mathrm{max},n}}$ increases sufficiently slowly, then
uniform analogues of the results in Theorems 6 and 7 may be anticipated by
establishing Berry-Esseen-type bounds for dependent processes. Note, for
instance, that Jirak, (2016) establishes a unidimensional Berry-Esseen bound
for sums of weakly dependent random variables. Note also that standard
dimension-dependent Berry-Esseen bounds for independent data scale as
$O(k^{1/4}n^{-1/2})$ for total variation distance, implying bounds
$O(k^{3/4}n^{-1/2})$ when passing to Euclidean distance, i.e.,
$\|\tilde{S}_{\mathcal{C}}-\tilde{S}^{*}_{\mathcal{C}}\|_{2}$. Note also that
Theorem 3 can be extended to sequences ${k_{\mathrm{max},n}}$ sufficiently by
the same argument as given in the supplemental material with
$\check{n}:=\lfloor n/k\rfloor$ replacing $n$ in the appropriate places.

A final case to consider is to fix $k$ but allow $\mathscr{C}_{n}$ to contain
many different partitions. For example, this case could arise by considering
the output of $k$-medoids algorithm with many different initial starting
medoids. In this case, arguments from standard empirical process theory may be
used. For brevity, this paper does not carry out such a program further. Note,
however, that the set of $k$-tuples of points in a bounded subset of Euclidean
space has known metric entropy properties. These metric entropies may be used
to control the metric entropy of tuples in $\mathsf{X}_{n}$ after embedding
$\mathsf{X}_{n}^{1-\varepsilon}$ into a Euclidean space for suitable
$\varepsilon>0$.

The conditions imposed for Theorem 6 do not explicitly require any bounds on
$|\mathscr{C}_{n}|$ or $|\mathcal{C}|,\mathcal{C}\in\mathscr{C}_{n}$.
Condition 7(iii) used in Theorem 7 does require bounds on the
$|\mathscr{C}_{n}|$ or $|\mathcal{C}|,\mathcal{C}\in\mathscr{C}_{n}$ which are
independent of $n$, but this requirement can be dispensed if $w(S)$ is
sufficiently regular. This regularity condition is given in Condition 7(iii)’
in the supplemental material. There also is provided a proposition showing
that $w(S)=|t(S)|$ satisfies the required regularity condition.

##  Appendix A Implementation Details

###  A.1 Implementation of $k$-medoids clustering

This section states the $k$-medoids algorithm used in the main text. For
finite $(\mathsf{X},\mathrm{d})$, let $\mathsf{C}\subseteq\mathsf{X}$ and
define the cost of a cluster $\mathsf{C}$ with medoid $i\in\mathsf{X}$ to be
$\text{cost}(\mathsf{C},i)=\sum_{j\in\mathsf{C}}\mathrm{d}(i,j)^{2}.$ The
total cost for a partition $\mathcal{C}$ and set of medoids
$\\{i_{\mathsf{C}}\\}_{\mathsf{C}\in\mathcal{C}}$ is defined by summing over
clusters $\text{total
cost}(\mathcal{C},\\{i_{\mathsf{C}}\\}_{\mathsf{C}\in\mathcal{C}})=\sum_{\mathsf{C}\in\mathcal{C}}\text{cost}(\mathsf{C},i_{\mathsf{C}})=\sum_{\mathsf{C}\in\mathcal{C}}\sum_{j\in\mathsf{C}}\mathrm{d}(i,j)^{2}.$

$k$-medoids

Input. $(\mathsf{X},\mathrm{d})$, $k$.

Procedure.

Initialize cluster centroids $\\{i_{1},...,i_{k}\\}\subset\mathsf{X}_{n}$
arbitrarily.

While total cost decreases,

(a) For each $g\leqslant k$, for each $j\notin\\{i_{1},...,i_{k}\\}$ compute
the cost with new medoids

$\\{i_{1},...,i_{g-1},j,i_{g+1},...i_{k}\\}$;

(b) Update to new medoids if total cost decreases.

Output. $\mathcal{C}$ with $|\mathcal{C}|=k$ .

###  A.2 Implementation of cluster-based inference

This section gives complete implementation details for CCE, IM, and CRS in the
empirical example and simulation sections as well as SK in the simulation
section.

####  A.2.1 Cluster-based methods (CCE, IM, CRS)

To describe the implementation, introducing vector and matrix notation is
helpful. Let $N_{\mathrm{pan}}$ and $T_{\mathrm{pan}}$ be the cross-sectional
dimension and time dimensions (i.e., sample sizes) of the given panel dataset.
Write $n$ as the product $n=N_{\mathrm{pan}}T_{\mathrm{pan}}$. Let $Y$, $X$,
$W$, and $U$ be $n$-row matrices obtained by stacking $Y_{de}$, $X_{de}$,
$(W_{de}^{\prime},1)$, and $U_{de}$, respectively. In the IV model, let $Z$
and $V$ be $n$-row matrices obtained by stacking $Z_{de}$ and $V_{de}$,
respectively. Let $M_{A}=I_{n}-A(A^{\prime}A)^{-1}A^{\prime}$ for some matrix
$A$. For some $B\in\mathbb{R}^{n\times n}$ with rank $\ell$, let $P_{B}$ be
some bijective linear transformation from the subspace orthogonal to the one
spanned by the columns of $B$ (thus $(n-\ell)$-dimensional), to
$\mathbb{R}^{n-\ell}$.

  1. Step 1.

Form candidate clusters. For $k=2,...,k_{max}=\lceil n^{1/3}\rceil$ where $n$
is the total sample size, apply $k$-medoids with dissimilarity matrix from the
data to obtain $\mathcal{C}^{(k)}$, a collection of $k$ partitions of the
observations.

  2. Step 2.

Fit a parametric model to covariance structure. The implementation details in
this step differ slightly in the cases of OLS and IV estimation.

     * •

In the OLS model, let $\widehat{U}$ be the vector of residuals from the full-
sample least-squares estimation. Then, under $U\sim\mathrm{N}(0,\Sigma)$,

$P_{M_{[X,W]}}\widehat{U}\sim\mathrm{N}(0,\Sigma_{PMU}),$

where

$\Sigma_{PMU}=P_{M_{[X,W]}}M_{[X,W]}\Sigma M_{[X,W]}P_{M_{[X,W]}}^{\prime}.$

Note that the covariance matrix $\Sigma_{PMU}$ is made non-singular by
applying the matrix $P_{M_{[X,W]}}$ to $\widehat{U}$. $\Sigma$ is estimated by
QMLE using an exponential covariance model with parameter
$\tau=(\tau_{1},\tau_{2},\tau_{3})^{\prime}$,

$\Sigma_{de,d^{\prime}e^{\prime}}(\tau)=cov[U_{de},U_{d^{\prime}e^{\prime}};\tau]=\exp(\tau_{1})\exp(-\tau_{2}^{-1}\|L_{d}-L_{d^{\prime}}\|_{2}-\tau_{3}^{-1}|e-e^{\prime}|),$

which, in the BASELINE case, is the correct model. To implement, calculate

$\widehat{\tau}=\underset{\tau}{\arg\max}\left\\{\frac{1}{2}\log\det(\Sigma_{PMU}(\tau))+\frac{1}{2}\widehat{U}^{\prime}P_{M_{[X,W]}}^{\prime}(\Sigma_{PMU}(\tau))^{-1}P_{M_{[X,W]}}\widehat{U}\right\\},$

where

$\Sigma_{PMU}(\tau)=P_{M_{[X,W]}}M_{[X,W]}\Sigma(\tau)M_{[X,W]}P_{M_{[X,W]}}^{\prime}$

and
$\Sigma(\tau)=(\Sigma_{de,d^{\prime}e^{\prime}}(\tau))_{de,d^{\prime}e^{\prime}}$
is the implied covariance matrix of $U$ under $\tau$. The covariance matrix
estimator is thus $\Sigma(\widehat{\tau})$.

     * •

In the IV model, the covariance matrices for the structural and first-stage
equations are estimated separately. Let
$\widehat{U}=M_{W}Y-M_{W}X\widehat{\theta}$ and
$\widehat{V}=M_{W}X-M_{W}Z\widehat{\pi}$, where $\widehat{\theta}$ is the 2SLS
estimator for $\theta_{0}$ and $\widehat{\pi}$ is the least-square estimator
for $\pi$. Then, the covariance matrices for $U$ and $V$ are estimated by
solving

$\widehat{\tau}^{\varepsilon}=\underset{\tau}{\arg\max}\left\\{\frac{1}{2}\log\det(\Sigma_{PM\varepsilon}(\tau))+\frac{1}{2}\widehat{\varepsilon}^{\prime}P_{M_{W}}^{\prime}(\Sigma_{PM\varepsilon}(\tau))^{-1}P_{M_{W}}\widehat{\varepsilon}\right\\},$

where
$\Sigma_{PM\varepsilon}=P_{M_{W}}M_{W}\Sigma(\tau)M_{W}P_{M_{W}}^{\prime}$,
$\Sigma(\tau)=(\Sigma_{de,d^{\prime}e^{\prime}}(\tau))_{de,d^{\prime}e^{\prime}}$,
$\Sigma_{de,d^{\prime}e^{\prime}}(\tau)$ is as previously defined, and
$\varepsilon$ is either $U$ or $V$. Then, the covariance estimators for $U$
and $V$ are $\widehat{\Sigma}_{U}=\Sigma(\widehat{\tau}^{U})$ and
$\widehat{\Sigma}_{V}=\Sigma(\widehat{\tau}^{V})$, respectively. Finally,
estimate the correlation between first and second stage errors with
$\widehat{\rho}$, the empirical correlation between
$\widehat{\Sigma}_{U}^{-1/2}\widehat{U}$ and
$\widehat{\Sigma}_{V}^{-1/2}\widehat{V}$.

  3. Step 3.

Simulate data. This step simulates size and power for all candidate partitions
$\mathcal{C}\in\mathscr{C}=\\{\mathcal{C}^{(2)},...,\mathcal{C}^{(k_{\mathrm{max}})}\\}$.
Given the covariance estimator(s) from Step 2, simulate independent copies of
the observable data for each $b=1,...,B$ as follows for each
$\theta\in\\{-10/\sqrt{n},-9/\sqrt{n},\dots,-1/\sqrt{n},0,1/\sqrt{n},2/\sqrt{n},\dots,10/\sqrt{n}\\}$.
(We use $B=10000$ in the empirical example $B=1000$ in the simulation
experiments.)

     * •

In the OLS model, draw $U^{b}$ from the distribution
$N(0,\Sigma(\widehat{\tau}))$. Reproduce data by

$Y_{de}^{b}=\widehat{\alpha}+\theta
X_{de}+W_{de}^{\prime}\widehat{\gamma}+U_{de}^{b},$

where $\widehat{\alpha}$ and $\widehat{\gamma}$ are full-sample least-square
estimators, and $U_{de}^{b}$ is the $de$ element on $U^{b}$.

     * •

In the IV model, draw $(U^{b},V^{b})$ such that

$\begin{pmatrix}U^{b}\\\
V^{b}\end{pmatrix}\sim\mathrm{N}\left(0,\begin{bmatrix}\widehat{\Sigma}_{U}&\widehat{\rho}\widehat{\Sigma}_{U}^{1/2}(\widehat{\Sigma}_{V}^{1/2})^{\prime}\\\
\widehat{\rho}\widehat{\Sigma}_{V}^{1/2}(\widehat{\Sigma}_{U}^{1/2})^{\prime}&\widehat{\Sigma}_{V}\end{bmatrix}\right).$

Reproduce data by

$\begin{cases}Y_{de}^{b}=\widehat{\alpha}+\theta
X_{de}^{b}+W_{de}^{\prime}\widehat{\gamma}+U_{de}^{b}\\\
X_{de}^{b}=\widehat{\mu}+\widehat{\pi}Z_{de}+X_{de}^{\prime}\widehat{\xi}+V_{de}^{b}\end{cases},$

where $\widehat{\alpha}$ and $\widehat{\gamma}$ are full-sample 2SLS
estimators, $\widehat{\mu}$, $\widehat{\pi}$, $\widehat{\xi}$ are full-sample
least-square estimators for the first-stage equation, and $U_{de}^{b}$, and
$V_{de}^{b}$ are respectively the $de$ elements of $U^{b}$ and $V^{b}$.

  4. Step 4.

Calculate Type I and Type II error rates. For each partition size
$k=2,...,k_{\max}$ and for each $a\in[0,.05]$, compute simulated Type I error
rate $\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{IM}(a),k)$,
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{CRS}(a),k)$,
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{CCE}(a),k)$ by conducting
the corresponding test $H_{0}:\theta_{0}=0$ on each simulated dataset from
Step 3. Set $\widehat{\alpha}_{\text{IM},k}$,
$\widehat{\alpha}_{\text{IM},k}$, $\widehat{\alpha}_{\text{IM},k}$ to be the
largest value $a\in[0,\alpha]$ such that
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{IM}(a),k)\leqslant\alpha$,
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{CRS}(a),k)\leqslant\alpha$,
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{CCE}(a),k)\leqslant\alpha$.

For each partition $k=2,...,k_{\max}=8$, compute simulated average Type II
error rate
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{IM}(\widehat{a}_{\text{IM},k}),k)$,
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{CRS}(\widehat{a}_{\text{CRS},k}),k)$,
$\widehat{\mathrm{Err}}_{\mathrm{Type-I}}(\mathrm{CCE}(\widehat{a}_{\text{CCE},k}),k)$
using
$\widehat{a}_{\text{IM},\mathcal{C}},\widehat{a}_{\text{CRS},\mathcal{C}},\widehat{a}_{\text{CCE},\mathcal{C}}$
and by conducting the corresponding test
$H_{0}:\theta_{0}=\theta^{\mathrm{alt}}$ for
$\theta^{\mathrm{alt}}\in\\{-10/\sqrt{n},-9/\sqrt{n},\dots$
$\dots,-1/\sqrt{n},1/\sqrt{n},2/\sqrt{n},\dots,10/\sqrt{n}\\}$ on each
simulated dataset from Step 3 and averaging the Type 2 error obtained at each
hypothesized value.

  5. Step 5.

Solution to (2.9). Set $\widehat{k}_{\mathrm{CCE}}$,
$\widehat{k}_{\mathrm{IM}}$, $\widehat{k}_{\mathrm{CRS}}$ corresponding to the
highest simulated average power. Set
$\widehat{T}_{\text{IM}(\alpha)}=T_{\text{IM}(\widehat{a}_{\text{IM},k}),\mathcal{C}^{\widehat{k}_{\text{IM}}}},\widehat{T}_{\text{CRS}(\alpha)}=T_{\text{CRS}(\widehat{a}_{\text{CRS},\mathcal{C}}),\mathcal{C}^{\widehat{k}_{CRS}}},\widehat{T}_{\text{CCE}(\alpha)}=\widehat{T}_{\text{CCE}(\widehat{a}_{\text{CCE},\mathcal{C}}),\mathcal{C}^{\widehat{k}_{\text{CCE}}}}.$

####  A.2.2 Spatial-HAC (SK)

The method SK is based on the spatial-HAC estimator proposed by Sun and Kim,
(2015). We apply the methods in Sun and Kim, (2015) to transform an irregular
lattice to a regular integer lattice and to deal with locations that do not
form a rectangle. Sun and Kim, (2015) require the input of smoothing
parameters $(K_{1},K_{2})$. We apply the method recommended by Lazarus et al.,
(2018); i.e., we let
$K_{1}=K_{2}=\left[\sqrt{0.4N_{\mathrm{pan}}^{2/3}}\right]$ where
$N_{\mathrm{pan}}$ is the number of locations.

## References

  * Andrews, (1991) Andrews, D. W. K. (1991).  Asymptotic normality of series estimators for nonparametric and semiparametric regression models.  Econometrica, 59(2):307–345.
  * Assoud, (1977) Assoud, P. (1977).  Espaces Métriques, Plongements, Facteurs.  Doctoral Dissertation, Université de Paris XI, 91405 Orsay France.
  * Bai et al., (2020) Bai, J., Choi, S. H., and Liao, Y. (2020).  Standard errors for panel data models with unknown clusters.  Journal of Econometrics.
  * Bester et al., (2011) Bester, C. A., Conley, T. G., and Hansen, C. B. (2011).  Inference with dependent data using cluster covariance estimators.  Journal of Econometrics, 165(2):137 – 151.
  * Bester et al., (2016) Bester, C. A., Conley, T. G., Hansen, C. B., and Vogelsang, T. J. (2016).  Fixed-b asymptotics for spatially dependent robust nonparametric covariance matrix estimators.  Econometric Theory, 32(1):154–186.
  * Bickel and Levina, (2008) Bickel, P. J. and Levina, E. (2008).  Covariance regularization by thresholding.  The Annals of Statistics, 36(6):2577–2604.
  * Bolthausen, (1982) Bolthausen, E. (1982).  On the central limit theorem for stationary mixing random fields.  Ann. Probab., 10(4):1047–1050.
  * Cai, (2021) Cai, Y. (2021).  Panel Data with Unknown Clusters.  ArXiv e-prints.
  * Canay et al., (2017) Canay, I. A., Romano, J. P., and Shaikh, A. M. (2017).  Randomization tests under an approximate symmetry assumption.  Econometrica, 85(3):1013–1030.
  * Condra et al., (2018) Condra, L. N., Long, J. D., Shaver, A. C., and Wright, A. L. (2018).  The logic of insurgent electoral violence.  American Economic Review, 108(11):3199–3231.
  * Conley, (1999) Conley, T. G. (1999).  GMM estimation with cross sectional dependence.  Journal of Econometrics, 92:1–45.
  * Conley and Dupor, (2003) Conley, T. G. and Dupor, B. (2003).  A spatial analysis of sectoral complementarity.  Journal of Political Economy, 111(2):311–352.
  * Conley et al., (2018) Conley, T. G., Gonçalves, S., and Hansen, C. (2018).  Inference with dependent data in accounting and finance applications.  Journal of Accounting Research, 56:1139–1203.
  * Conley and Topa, (2002) Conley, T. G. and Topa, G. (2002).  Socio-economic distance and spatial patterns in unemployment.  Journal of Applied Econometrics, 17(4):303–327.
  * Ferman, (2019) Ferman, B. (2019).  A simple way to assess inference methods.  ArXiv e-prints.
  * Hansen and Lee, (2019) Hansen, B. E. and Lee, S. (2019).  Asymptotic theory for clustered samples.  Journal of Econometrics, 210:268 – 290.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009).  Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer, New York, NY.
  * Ibragimov and Müller, (2010) Ibragimov, R. and Müller, U. K. (2010).  $t$-statistic based correlation and heterogeneity robust inference.  Journal of Business & Economic Statistics, 28(4):453–468.
  * Jenish and Prucha, (2009) Jenish, N. and Prucha, I. R. (2009).  Central limit theorems and uniform laws of large numbers for arrays of random fields.  Journal of Econometrics, 150(1):86–98.
  * Jirak, (2016) Jirak, M. (2016).  Berry–Esseen theorems under weak dependence.  Ann. Probab., 44(3):2024–2063.
  * Kelejian and Prucha, (2001) Kelejian, H. H. and Prucha, I. (2001).  On the asymptotic distribution of the Moran I test statistic with applications.  Journal of Econometrics, 104:219–257.
  * Kelejian and Prucha, (2007) Kelejian, H. H. and Prucha, I. R. (2007).  HAC estimation in a spatial framework.  Journal of Econometrics, 140(1):131–154.