The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
74,667 characters
10mmInference for High-Dimensional Network Data
\maketitle
\begin{abstract}
\begin{spacing}{1.15}
We develop a novel method of inference for network-dependent high-dimensional random vectors. Dependence is characterized via a functional dependence measure based on graph distance, allowing the approximation theory to capture the interaction between the decay of dependence and the growth of network neighborhoods. We establish Gaussian approximation results for the maximum norm under finite-moment and sub-Weibull conditions, providing explicit conditions under which the dimension may increase with the network size. We also propose a high-dimensional network HAC covariance estimator and establish its convergence properties, yielding a feasible procedure for simultaneous inference. Simulation studies demonstrate favorable finite-sample performance of the proposed method. We apply the procedure to study how spillover effects vary with an index of network homophily by constructing confidence bands for the conditional spillover-effect function. The application reveals heterogeneity and local significance that would be obscured by conventional low-dimensional inference.
\smallskip\noindent
{\bf JEL Codes:} C12, C21
\smallskip\noindent
{\bf Keywords:} Gaussian approximation, high dimensional data, inference, network
\end{spacing}
\end{abstract}
\section{Introduction}
Network-dependent data arise in a wide range of empirical settings in which individuals, firms, regions, markets, or institutions are connected. Prominent examples include social interactions and peer effects, production and supply-chain networks, financial networks, trade networks, and interference in experiments conducted on social or geographic networks. In these settings, statistical dependence is governed by graph distance. Units that are close in the network may exhibit strong dependence, while the influence between units may decay as their graph distance increases. An important feature of network dependence is that its aggregate strength is determined jointly by the rate at which dependence decays with graph distance and the rate at which network neighborhoods grow.
While a growing literature has developed asymptotic theory for network-dependent data, much of the existing theory concerns low-dimensional statistics. At the same time, many empirical problems require simultaneous inference for high-dimensional parameters. Examples include inference for functional treatment-effect parameters (e.g., conditional average treatment effects or CATEs, continuous treatment effects, etc.) and inference for structural parameters that are (partially) identified by conditional moment (in)equalities. High-dimensional Gaussian approximation provides a powerful foundation for such inference, but existing results are primarily designed for independent observations, time series, spatial random fields, local dependence structures, or exchangeable arrays. These results do not directly accommodate dependence that propagates over a general observed network, where the accumulation of dependence is governed by both graph distance and network topology.
This paper develops high-dimensional Gaussian approximation and a method of simultaneous inference for sums of network-dependent random vectors. We consider observations indexed by the nodes of an observed network and characterize their dependence through a functional dependence measure based on perturbations of primitive shocks. The effect of a shock on an observation decays with their graph distance, while the network itself determines how many potentially dependent observations occur at each distance. This formulation provides a direct way to combine the strength and spatial propagation of dependence with the topology of the network.
Our main contribution is to establish Gaussian approximation for the maximum of a high-dimensional sum of network-dependent observations. We allow the dimension of the random vectors to increase with the network size and provide results under both finite-moment and sub-Weibull conditions. The approximation theory explicitly accommodates the interaction among the dimension of the statistic, moment conditions, decay of functional dependence, and growth of network neighborhoods. The proofs combine functional-dependence coupling and localization with a partition of the network into approximately independent interior clusters separated by buffer regions. This construction reduces the network-dependent problem to a high-dimensional Gaussian approximation for independent cluster sums while controlling the errors introduced by localization, buffers, and cross-cluster dependence.
Our second contribution is to make the Gaussian approximation feasible for simultaneous inference. We develop a high-dimensional network HAC covariance estimator in which observations are weighted according to graph distance and establish its consistency in high dimensions. Combining the covariance estimator with the Gaussian approximation yields feasible critical values for the maximum statistic, and hence simultaneous confidence intervals and hypothesis testing for high-dimensional network statistics.
Our framework is related to several strands of the literature. A large literature studies cross-sectional dependence through spatial random fields. When dependence is driven by geographic distance or a similar spatial metric, random fields indexed by Euclidean locations provide a natural framework. Important contributions include \cite{conley1999gmm}, \cite{KelejianPrucha2007}, \cite{kim2011spatial}, \cite{BesterConleyHansen2011}, \cite{JenishPrucha2009}, and \cite{jenish2012spatial}. Our setting instead takes the observed network and its graph distance as the primitive structure governing statistical dependence.
More closely related is the recent literature on network dependence. \cite{KMS2021} introduce a general framework based on $\psi$-dependence and graph distance, and establish LLNs, CLTs, and network HAC variance estimation. Building on this framework, \cite{Leung2022} studies causal inference under interference, while \cite{Leung2023} characterizes the validity of cluster-robust inference. See also \cite{kojevnikov2021bootstrap}, \cite{johnsson2021estimation}, and \cite{LeeSong2019} for related inference methods. This literature, however, has not investigated Gaussian approximation or inference for high-dimensional vectors. Our paper complements this literature
by developing such methods under network dependence.
Recently, \cite{gao2026coupling} study graph-dependent data and establish novel maximal inequalities for empirical processes. Our paper provides a complementary contribution with a different objective. While \cite{gao2026coupling} develop maximal inequalities for uniformly controlling empirical processes, our focus is on high-dimensional Gaussian approximation for processes that need not be stochastically equicontinuous.
Our work is also related to limit theory for cross-sectionally dependent data based on martingale and mixingale methods. \cite{KuersteinerPrucha2013} establish a CLT for martingale difference sequences. \cite{kuersteiner2020dynamic} incorporates network dependence into a martingale framework, while \cite{kuersteiner2019limit} develops conditional spatial mixingale limit theory accommodating stochastic network formation. Related results based on dependency graphs and local dependence include \cite{stein1972bound}, \cite{janson1988normal}, \cite{BaldiRinott1989}, \cite{rinott1996multivariate}, and \cite{chen2004normal}; applications to network data include \cite{AronowSamii2017}, \cite{Leung2020}, and \cite{song2018measuring}.
Our theoretical results build on the literature on high-dimensional Gaussian approximation. For independent random vectors and empirical processes, important contributions include \cite{CCK2013,CCK2014,chernozhukov2015comparison,CCK2016,CCK2017}, \cite{CCK2019}, \cite{deng2020beyond}, \cite{KuchibhotlaMukherjeeBanerjee2021}, and \cite{chernozhukov2023nearly}. High-dimensional Gaussian approximation has subsequently been extended to several dependent-data settings. For time series, \cite{ZW2017} establish Gaussian approximation under the physical dependence measure \citep{wu2005nonlinear}, while \cite{ZhangCheng2018} derive related results under functional dependence and develop inference based on the non-overlapping block bootstrap. For local dependence, \cite{fang2021high} establish high-dimensional CLTs under dependency graph structures. \cite{chiang2023inference} establish a high-dimensional Gaussian approximation for exchangeable arrays. For spatial dependence, \cite{kurisu2024gaussian} extend high-dimensional Gaussian approximation to spatial random fields. Our results complement this literature by allowing dependence to propagate over an observed network and by making explicit how network neighborhood growth interacts with dependence decay in determining the validity of high-dimensional Gaussian approximation.
We investigate the finite-sample performance of the proposed methods through Monte Carlo simulations. The empirical application examines heterogeneous spillover effects and demonstrates how the proposed procedure can be used to construct simultaneous confidence bands for network-based treatment- and spillover-effect parameters.
\medskip\noindent
{\bf Organization:}
Section \ref{sec:setup} introduces the setting, the statistic of interest, and analysis tools. Section \ref{sec:GA} presents the main Gaussian approximation results. Section \ref{sec:estimation-cov} develops a method of inference. Sections \ref{sec:simGA} and \ref{sec:HAC} study the finite sample performance through Monte Carlo simulations. Section \ref{sec:application} applies the method to heterogeneous spillover analysis in a network experiment. Proofs and additional technical details are collected in the appendix.
\medskip\noindent
{\bf Notation:}
For a random variable $X$ and $q>0$,
we write $X \in \mathcal{L}^q$ if
$
\|X\|_q := \big( \mathbb{E}|X|^q \big)^{1/q} < \infty,
$
and for a vector $v = (v_1,\dots,v_p)^\top$, let the norm-$s$ length be
$|v|_s = \left( \sum_{j=1}^p |v_j|^s \right)^{1/s}$ for $s \ge 1.$
Write the $p \times p$ identity matrix as $\mathrm{Id}_p$.
For two sequences of positive numbers $(a_n)$ and $(b_n)$,
we write $a_n \asymp b_n$ (resp., $a_n \lesssim b_n$ or
$a_n \ll b_n$) if there exists some constant $C>0$ such that
$C^{-1} \le \frac{a_n}{b_n} \le C
\quad
(\text{resp., } \frac{a_n}{b_n} \le C
\text{ or } \frac{a_n}{b_n} \to 0)$
for all large $n$.
We use $C, C_1, C_2, \dots$ to denote positive constants whose values may differ from place to place.
A constant with a symbolic subscript is used to emphasize the dependence of the value on the subscript.
We allow $p = p_n$ to increase with $n$,
and all asymptotic statements are taken as $n\to\infty$.
\section{Network Dependence Framework}{\label{sec:setup}}
\subsection{Network Topology}
Suppose that a researcher observes a set of cross-sectional units, denoted by $V_n=\{1,\dots,n\}$.
Modeling statistical dependence across these units typically requires a notion of distance.
In time-series applications, dependence is indexed by temporal distance, whereas in spatial applications, it is typically indexed by the Euclidean distance.
In the current paper, we consider dependence induced by the topology of an observed network.
Specifically, suppose that a researcher observes an undirected network $\mathcal{G}_n=(V_n,E_n)$ on $V_n$, where
$
E_n\subseteq \{\{i,u\}: i,u\in V_n,\ i\neq u\}
$
denotes the set of links.
For any $i,u\in V_n$, let $d_n(i,u)$ denote the distance between units $i$ and $u$ in $\mathcal{G}_n$, defined as the length of the shortest path connecting them.
The resulting function $d_n$ defines a metric on $V_n$ and serves as the fundamental notion of proximity throughout the paper.
Define the volume neighborhood $\mathcal N_n(i,r)$ as the set of nodes within distance $r$ of node $i$, and the shell neighborhood $\mathcal N_n^\partial(i,r)$ as the set of nodes at exactly distance $r$ from node $i$:
\[
\mathcal N_n(i,r)
=
\{u\in V_n:d_n(i,u)\le r\}
\qquad\text{and}\qquad
\mathcal N_n^\partial(i,r)
=
\{u\in V_n:d_n(i,u)=r\}.
\]
Their cardinality will be denoted by
\begin{equation}\label{eq:neighbor size}
N_n(i,r)
=
|\mathcal N_n(i,r)|
\qquad\text{and}\qquad
N_n^\partial(i,r)
=
|\mathcal N_n^\partial(i,r)|.
\end{equation}
The sizes of these local neighborhoods play a central role in our analysis.
Because dependence propagates through the network, the cumulative dependence surrounding a node depends not only on the strength of dependence but also on the rate at which neighborhoods expand.
\subsection{ Functional Dependence}
Network dependence can be characterized through a variety of weak dependence notions, including mixing conditions, dependency graphs, and functional dependence.\footnote{These alternative approaches are also encompassed by the $\psi$-dependence framework of \citet{kojevnikov2021bootstrap}.}
Among these alternatives, we focus on the functional dependence approach, originally introduced by \citet{wu2005nonlinear}, as it provides a convenient framework for deriving localization bounds, nonasymptotic probability inequalities, and high-dimensional Gaussian approximations required for our analysis.
Related functional-dependence ideas have been widely used in the high-dimensional time-series literature. See \citet{ZW2017} and \citet{ZhangCheng2018} for example.
The functional dependence framework represents each observation as a function of primitive shocks and measures dependence through the effects of shock perturbations.
In the network setting, this formulation is particularly natural, as the impact of a shock can be allowed to decay with its graph distance from the observation.
Specifically, let
\[
\varepsilon_n:=\{\varepsilon_{n,u}:u\in V_n\}
\]
be a collection of mutually independent primitive shocks indexed by the nodes of the network.
Each observation is generated as a measurable function of the entire shock field,
\[
X_{n,i}=H_{n,i}(\varepsilon_n)\in\mathbb R^p,
\qquad
\mathbb E(X_{n,i})=0
\quad
\text{for all } i\in V_n,
\]
where $H_{n,i}$ is allowed to vary across units, thereby accommodating heterogeneity.
Write
\[
X_{n,i}
=
(X_{n,i1},\ldots,X_{n,ip})^\top
\]
for the coordinate representation of the $p$-dimensional vector $X_{n,i}$.
This framework encompasses a broad class of linear and nonlinear network processes.
Throughout the paper, the network \(\mathcal{G}_n=(V_n,E_n)\) is treated as observed, and all probabilistic statements are understood conditionally on the realized network structure.
For any node \(u\in V_n\), let \(\varepsilon_n^{*(u)}\) denote the coupled shock field obtained by replacing \(\varepsilon_{n,u}\) with an independent copy \(\varepsilon'_{n,u}\) while leaving all remaining shocks unchanged.
Define
\[
X_{n,i}^{*(u)}
:=
H_{n,i}(\varepsilon_n^{*(u)}).
\]
Following the functional dependence approach of \citet{ZW2017}, we measure the influence of the shock at node \(u\) on the \(j\)-th coordinate of observation \(i\) by the \(L_q\) distance
\[
\delta_{n,i,u,q,j}
:=
\left\|
X_{n,ij}
-
X^{*(u)}_{n,ij}
\right\|_q
\]
between \(X_{n,ij}\) and its coupled counterpart $X_{n,ij}^{\ast(u)}$.
This quantity \(\delta_{n,i,u,q,j}\) measures the sensitivity of \(X_{n,ij}\) to a perturbation of the shock at node \(u\).
Large values indicate that the shock at \(u\) has a substantial effect on observation \(i\).
Although the dependence measure may vary across units and shocks, network dependence is fundamentally governed by network distance. We therefore summarize the maximal influence of shocks from distance $s$ through the distance-wise functional dependence measure
\[
\delta_{n,s,q,j}
:=
\sup_{i\in V_n}
\sup_{u:d_n(i,u)=s}
\delta_{n,i,u,q,j}.
\]
To quantify the cumulative influence of distant shocks,
we define the cumulative tail dependence envelop
\(\Delta_{n,m,q,j}\) and its high-dimensional envelope
\(\Psi_{n,q}(m)\) by
\begin{equation}\label{eq:psi}
\Delta_{n,m,q,j}
=
\sup_{i\in V_n}\sum_{s\ge m}
N_n^\partial(i,s)\,
\delta_{n,s,q,j},
\qquad
\Psi_{n,q}(m)
=
\max_{1\le j\le p}
\Delta_{n,m,q,j},
\end{equation}
The special case \(m=0\) plays a particularly important role.
Define
\[
\Theta_{n,q,j}
:=
\Delta_{n,0,q,j}.
\]
Following \citet{ZW2017}, we refer to \(\Theta_{n,q,j}\) as the dependence-adjusted norm of the \(j\)-th coordinate, as it combines both moment magnitude and cumulative network dependence.
The quantities \(\delta_{n,s,q,j}\), \(\Delta_{n,m,q,j}\), and \(\Psi_{n,q}(m)\) will serve as the primary measures of network dependence throughout the paper.
\section{Gaussian Approximations}{\label{sec:GA}}
This section presents the main assumptions and results on Gaussian approximations.
To facilitate this objective, we introduce the notation
$$
T_X=\sum_{i \in V_n}X_{n,i},
\qquad\text{and}\qquad
\Sigma_n=\mathrm {Var}(T_X/\sqrt{n})=n^{-1}\sum_{i,i'\in V_n}\mathrm {Cov}(X_{n,i},X_{n,i'})
$$
for the network sum and scaled variance, respectively.
Let $\Sigma_0=\mathrm{diag}(\Sigma_n)=\mathrm{diag}(\sigma_{11},\cdots,\sigma_{pp})$ denote the diagonal matrix of $\Sigma_n$, and let $D_0=\Sigma_0^{1/2}$ be its square root.
The location is normalized to $\mathbb E(X_{n,i})=0$ for all $i\in V_n$ throughout.
We aim to establish the Gaussian approximation:
\begin{equation}{\label{eq:GA}}
\rho_n
:= \sup_{u \ge 0}
\left|
\mathbb{P}\!\left(
\,\bigl| D_0^{-1} T_X/\sqrt{n} \bigr|_\infty \le u
\right)
-
\mathbb{P}\!\left(
\bigl| D_0^{-1} Z \bigr|_\infty \le u
\right)
\right|
\to 0
\end{equation}
where $Z \sim N(0,\Sigma_n)$.
\subsection{Cluster Partition}
The key step in establishing the Gaussian approximation \eqref{eq:GA} is to reduce the network-dependent statistic
to a sum of approximately independent cluster-level contributions.
Once such a representation is obtained, the Gaussian approximation theory for maxima of sums of independent random vectors \citep{CCK2013} becomes applicable.
To achieve this goal, we partition the network into many balanced clusters and show that dependence across clusters are asymptotically negligible after localization.
Specifically, we suppose that there exists a squence of cluster sets satisfying the following assumption.
\begin{assumption}[Regular Cluster Partition]
\label{ass:clusters}
The network $V_n$ can be partitioned\footnote{
That is, $\{\mathcal C_g\}_{g=1}^{G_n}$ satisfies $\bigcup_{g=1}^{G_n}\mathcal C_g=V_n$ and $\mathcal C_{g_1}\cap\mathcal C_{g_2}=\varnothing$ for $g_1\neq g_2.
$}
into clusters
$
\{\mathcal C_g\}_{g=1}^{G_n}
$
satisfying the following properties:
\begin{itemize}
\item[(i)] $G_n\to\infty$ and $G_n=o(n)$.
\item[(ii)] There exist a sequence $l_n\to\infty$ and constants $0<\underline c<\bar c<\infty$
such that $\underline c\,l_n \le l_g \le \bar c\,l_n$ for each $g=1,\ldots,G_n$, where $l_g := |\mathcal C_g|$.
\end{itemize}
\end{assumption}
Assumption \ref{ass:clusters} paves the way for a blocking structure that will be used for our analysis.
The required property $G_n\to\infty$ in part (i) of the assumption guarantees that the statistic of interest can be represented as a sum of many cluster-level contributions, while the balance condition (ii) prevents any single cluster from dominating the statistic.
Given $\{\mathcal C_g\}_{g=1}^{G_n}$, we define the boundary
\[
\partial\mathcal C_g
:=
\Bigl\{
i\in\mathcal C_g:
\exists\,k\notin\mathcal C_g
\text{ such that }
(i,k)\in E_n
\Bigr\}
\]
of a cluster $\mathcal C_g$.
This object will be used to quantify interactions across clusters.
The aggregate boundary size is denoted by
\[
\eta_n(G_n)
:=
\sum_{g=1}^{G_n}
|\partial\mathcal C_g|.
\]
The following examples illustrate how balanced cluster partitions can be constructed for commonly used network models.
\begin{example}[Ring Network]
\label{ex:ring}
Consider the ring graph
$\mathcal G_n=(V_n,E_n)$
where each node is connected to its two nearest index neighbors.
That is, node $1$ is connected to nodes $n$ and $2$, node $i$ is connected to nodes $i-1$ and $i+1$ for $i \in \{2,...,n-1\}$, and node $n$ is connected to nodes $n-1$ and $1$.
The graph distance can be written as
$d_n(i,i') = \min\{|i-i'|,\,n-|i-i'|\}.$
Partition the nodes into contiguous segments of equal
length $a_n$, i.e.,
\[
\mathcal C_g
=
\{(g-1)a_n+1,\ldots,ga_n\}
\qquad\text{and}\qquad
G_n={n}/{a_n}.
\]
By construction, all clusters have the same size
$
l_g=a_n.
$
Therefore, setting
$a_n\to\infty$
and
$a_n=o(n)$
satisfies Assumption~\ref{ass:clusters}.
Since each cluster has exactly two boundary nodes, we have
$|\partial \mathcal C_g|=2$
and
$\eta_n(G_n)=2G_n.$
\end{example}
\begin{example}[Regular Random Geometric Graph]
\label{ex:rgg}
Consider a random geometric graph
\(\mathcal G_n=(V_n,E_n)\),
where each node \(i\in V_n\) is associated with a random spatial location
$
U_{n,i}\in[0,1]^2.
$
Two nodes are connected whenever their Euclidean distance is no larger than a connection radius
$
r_n
=
c_r n^{-1/2}
$
for some fixed constant \(c_r>0\), that is,
\[
(i,i')\in E_n
\quad\Longleftrightarrow\quad
\lVert U_{n,i}-U_{n,i'}\rVert
\le r_n.
\]
The choice \(r_n\asymp n^{-1/2}\) makes the connection radius of the same order as the typical spacing between \(n\) nodes in a two-dimensional region.
Suppose further that there exist constants
\(c_{\mathrm{sep}},C_{\mathrm{Cov}}>0\), independent of \(n\), such that
\[
\min_{i\neq i'}
\lVert U_{n,i}-U_{n,i'}\rVert
\ge
c_{\mathrm{sep}}r_n
\qquad\text{and}\qquad
\sup_{x\in[0,1]^2}
\min_{i\in V_n}
\lVert x-U_{n,i}\rVert
\le
C_{\mathrm{Cov}}r_n.
\]
The first condition prevents excessive local crowding, while the second rules out spatial holes larger than the connection scale. Together, these two conditions allow irregular random configurations while ensuring that the node locations remain sufficiently regular at scale \(r_n\).
Partition the unit square into disjoint square regions
\(\{\mathcal R_g\}_{g=1}^{G_n}\)
with side length \(a_n\), and define
$
\mathcal C_g
=
\left\{
i\in V_n:
U_{n,i}\in\mathcal R_g
\right\}
$
and
$
G_n=a_n^{-2}.
$
Let
$
a_n\to0
$
and
$
\frac{a_n}{r_n}\to\infty.
$
Under the separation and covering conditions, the number of nodes contained in each square is of the same order as the square's area relative to the local spacing scale:
\[
l_g
=
|\mathcal C_g|
\asymp
\left(
\frac{a_n}{r_n}
\right)^2
\asymp
na_n^2
=
\frac{n}{G_n}
\]
uniformly over \(g\), with probability approaching one. Hence the resulting partition satisfies the balanced-cluster requirement in Assumption~\ref{ass:clusters}.
A node belongs to the cluster boundary if it lies within graph distance one of another cluster. Since every edge has Euclidean length at most \(r_n\), such nodes are contained in a spatial strip of width of order \(r_n\) along the boundary of the corresponding square region. The area of this strip is of order
$
a_n r_n.
$
By the same local-density regularity,
\[
|\partial\mathcal C_g|
\asymp
\frac{a_n r_n}{r_n^2}
=
\frac{a_n}{r_n}
\asymp
\left(
\frac{n}{G_n}
\right)^{1/2}
\]
uniformly over \(g\), with probability approaching one. Summing over all clusters gives
\[
\eta_n(G_n)
=
\sum_{g=1}^{G_n}
|\partial\mathcal C_g|
\asymp
n^{1/2}G_n^{1/2}.
\]
\end{example}
\subsection{ Localization and Approximate Independence}
Even if the network is partitioned into many clusters, observations belonging to different clusters generally remain dependent since shocks may propagate through arbitrarily long network paths.
In this section, we are going to localize such dependence inside the clusters and construct independent sums by imposing some conditions.
The first condition requires that the influence of a shock decays sufficiently fast with network distance, as formally stated below.
\begin{assumption}[Exponential Decay]
\label{ass:decay rate}
For some $q>4$, there exists $0<\rho<1$ such that
\[
\delta_{n,s,q,j}
\le
C_q\rho^s
\]
for all $s \ge 0$ and $j \in \{1,...,p\}$ for all $n$.
\end{assumption}
The second condition restricts the growth of network neighborhoods.
\begin{assumption}[Polynomial Shell and Volume Growth]
\label{ass:volume}
There exist constants $C_V>0$ and $d\ge1$, such that
\[
N_n(r)
=
\sup_{i\in V_n}
N_n(i,r)
\le
C_V(1+r)^{d}
\]
for all $r \ge 0$ and $n$.
\end{assumption}
Assumptions \ref{ass:decay rate} and \ref{ass:volume} are motivated by the fact that cumulative dependence depends on both the decay of shock propagator and the number of nodes affected, respectively.
Assumption \ref{ass:volume} is satisfied by bounded-growth networks such as rings and regular random geometric graphs.
However, it excludes graph sequences whose neighborhoods expand too rapidly relative to the dependence decay.
\begin{remark}\label{rem:local}
Under Assumption \ref{ass:decay rate} and \ref{ass:volume},
\begin{align*}
&\Theta_{n,q,j}
=
\Delta_{n,0,q,j}
\le
C_VC_q
\sum_{s\ge0}
(1+s)^d\rho^s
=
O(1)
\qquad\text{for all $j$,}\\
&\Psi_{n,q}(m)
:=
\max_{1\le j\le p}
\Delta_{n,m,q,j}
\lesssim
(1+m)^d\rho^m.
\end{align*}
\end{remark}
This bound motivates us to introduce the localized approximation
\begin{equation}{\label{eq:localize}}
X_{n,i}^{(m)}
=
\mathbb E
\bigl(
X_{n,i}
\mid
\mathcal F_{i,m}
\bigr),
\end{equation}
where
\[
\mathcal F_{i,m}
=
\sigma
\{
\varepsilon_{n,u}
:
d_n(i,u)\le m
\}.
\]
The localized approximation $X_{n,i}^{(m)}$ depends only on the shocks contained in the $m$-neighborhood of node $i$.
When $m$ is chosen of logarithmic order in $np$, the localized approximation error
\(
X_{n,i}-X_{n,i}^{(m)}
\)
is asymptotically negligible.
A formal statement is provided in Appendix \ref{sec:appendix}.
After this localization procedure, the dependence is effectively confined to local neighborhoods.
Consequently, cluster sums constructed from observations whose $m$-neighborhoods are entirely contained within different clusters become independent.
This observation forms the main bases for our blocking argument used in our proofs of Theorems \ref{th:GA_q} and \ref{th:GA_sub} to be presented in the following subsection.
\subsection{Gaussian Approximation Results}\label{sec:gaussian_approximation}
The assumptions in previous section ensure that localized observations can be organized into approximately independent cluster sums.
With this, we are now ready to apply the Gaussian approximation theory for independent observations \citep{CCK2013}, under the nondegeneracy condition on both the overall statistic and the cluster-level contributions.
\begin{assumption}[Variance Regularity]
\label{ass:variance lb}
There exist constants $c>0$ and $c_1>0$ such that
\begin{align*}
&\min_{1\le j\le p}\sigma_{jj}\ge c
\qquad\text{and}
\\
&\inf_{1\le g\le G_n}\inf_{1\le j\le p}
\mathrm{Var}\left(
l_g^{-1/2}\sum_{i\in \mathcal{C}_g}X_{n,ij}
\right)
\ge c_1.
\end{align*}
\end{assumption}
\begin{theorem}[Finite-Moment Gaussian Approximation]\label{th:GA_q}
Suppose that Assumption \ref{ass:clusters}, \ref{ass:decay rate}, \ref{ass:volume}, and \ref{ass:variance lb} hold.
If the partition satisfies
\begin{align*}
\mathrm{(i)} \ & \
G_n
\gg
\max\{
{\log(np)}^{7},
p^{2/(q-2)}
{\log(np)}^{3q/(q-2)}\},
\\
\mathrm{(ii)} \ & \ \max_{1\le g\le G_n}
\frac{|\partial \mathcal C_g|
(\log(np))^{d}
G_n}{n}
\to0,
\qquad\text{and}\\
\mathrm{(iii)} \ & \ \eta_n(G_n)
(\log(np))^{d}
\ll
n p^{-2/q}(\log p)^{-1},
\end{align*}
then the Gaussian approximation \eqref{eq:GA} holds.
\end{theorem}
Restriction (i) in the statement of this theorem requires that the partition contains sufficiently many clusters, whereas
restrictions (ii) and (iii) require that the contribution of boundary observations is asymptotically negligible both at the cluster level and in aggregate, respectively.
These general conditions are stated at a high-level.
Thus, let us revisit Examples~\ref{ex:ring} and~\ref{ex:rgg}, and discuss sufficient conditions in the contexts of these examples.
\bigskip\noindent
{\bf Example~\ref{ex:ring}} (Ring Network){\bf, Revisited.}
Consider the ring network introduced in Example~\ref{ex:ring}. Since
$|\partial\mathcal C_g|=2$ and $\eta_n(G_n)=2G_n$,
the lower and upper bounds on \(G_n\) given by conditions (i) and (iii) are
\[
G_n
\gg
\max\left\{
(\log np)^7,\,
p^{2/(q-2)}(\log np)^{3q/(q-2)}
\right\},
\]
and
\[
G_n
\ll
np^{-2/q}(\log np)^{-(d+1)}.
\]
For $p\ge n$, a sufficient
rate condition ensuring that a feasible choice of \(G_n\) exists is
\[
p^{\frac{2}{q-2}+\frac{2}{q}}
(\log p)^{d+1+3q/(q-2)}
=o(n).
\]
\bigskip\noindent
{\bf Example~\ref{ex:rgg}} (Regular Random Geometric Graph){\bf, Revisited.}
Consider the regular random geometric graph introduced in Example~\ref{ex:rgg}. Since
$
|\partial\mathcal C_g|
\asymp
\left({n}/{G_n}\right)^{1/2}
$
and
$
\eta_n(G_n)
\asymp
n^{1/2}G_n^{1/2},
$
the lower and upper bounds on \(G_n\) given by conditions (i) and (iii) are
\[
G_n
\gg
\max\left\{
(\log np)^7,\,
p^{2/(q-2)}(\log np)^{3q/(q-2)}
\right\},
\]
and
\[
G_n
\ll
np^{-4/q}(\log np)^{-(2d+2)}.
\]
For $p\ge n$, a sufficient
rate condition ensuring that a feasible choice of \(G_n\) exists is
\[
p^{\frac{2}{q-2}+\frac{4}{q}}
(\log p)^{2d+2+3q/(q-2)}
=o(n).
\]
\bigskip
Theorem \ref{th:GA_q}, and thus the discussions of the two examples above, focus on the case where the dependence-adjusted norm $\Theta_{n,q,j}$ exists for some $q>4$. This framework allows the dimension $p$ to grow with $n$ only at a polynomial rate. To have ultra high dimensionality for $p$, we can consider a stronger moment condition under which the dependence-adjusted norm $\Theta_{n,q,j}$ exists for all orders.
We extend Assumption \ref{ass:decay rate} to hold for all orders $q \ge 4$.
\begin{assumption}[Sub-Weibull Dependence Decay]
\label{ass:subweibull}
Assumption \ref{ass:decay rate} holds for all \(q\ge4\).
Moreover, there exist constants \(C>0\), \(\nu>0\), and \(\rho \in (0,1)\) such that
\[
\delta_{n,s,q,j}
\le
C q^\nu \rho^s
\]
holds for all $s\ge0$, $q\ge4$, and $j \in \{1,...,p\}$ for all $n$.
\end{assumption}
The parameter \(\nu\) governs the growth rate of higher-order moments and characterizes the tail behavior.
The next theorem states a counterpart of Theorem \ref{th:GA_q} for the case of sub-Weibull tails.
\begin{theorem}[Sub-Weibull Gaussian Approximation]{\label{th:GA_sub}}
Suppose that Assumption \ref{ass:clusters}, \ref{ass:volume}, \ref{ass:variance lb}, and \ref{ass:subweibull} hold.
If the partition satisfies
\begin{align*}
\mathrm{(i)} \ & \ G_n\gg (\log(np))^{\max\{7,\,4+2\nu\}},
\\
\mathrm{(ii)} \ & \ \max_{1\le g\le G_n}
\frac{|\partial \mathcal C_g|
(\log(np))^{d}
G_n}{n}
\to0,
\qquad\text{and}\\
\mathrm{(iii)} \ & \ \eta_n(G_n)\left(\log(np)\right)^{d}
\ll
n(\log p)^{-2-2\nu}.
\end{align*}
then the Gaussian Approximation in \eqref{eq:GA} holds.
\end{theorem}
\bigskip\noindent
{\bf Example~\ref{ex:ring}} (Ring Network){\bf, Revisited.}
Consider again the ring network introduced in Example~\ref{ex:ring}.
Now, suppose the conditions of Theorem~\ref{th:GA_sub} requiring sub-Weibull tails.
For $p\ge n$, all conditions reduce to
\[
(\log p)^{\max\{9+d+2\nu,\;6+d+4\nu\}}
=o(n).
\]
Equivalently, the admitted dimension can be ultra-high-dimensional:
\begin{align*}
\log p=o(n^c),
\qquad\text{where}\qquad
c=
\begin{cases}
1/(9+d+2\nu)
& \text{if }\frac{3}{2}\ge \nu\ge 0,\\[6pt]
1/(6+d+4\nu)
& \text{if }\nu\ge \frac{3}{2}.
\end{cases}
\end{align*}
\bigskip\noindent
{\bf Example~\ref{ex:rgg}} (Regular Random Geometric Graph){\bf, Revisited.}
Consider again the regular random geometric graph in Example~\ref{ex:rgg}.
Now, suppose the conditions of Theorem~\ref{th:GA_sub} requiring sub-Weibull tails.
For $p\ge n$, all conditions reduce to
\[
(\log p)^{\max\{11+2d+4\nu,\;8+2d+6\nu\}}
=o(n).
\]
Equivalently, the admitted dimension can be ultra-high-dimensional:
\begin{align*}
\log p=o(n^c),
\qquad\text{where}\qquad
c=
\begin{cases}
1/(11+2d+4\nu)
& \text{if }\dfrac{3}{2}\ge \nu\ge 0,\\[6pt]
1/(8+2d+6\nu)
& \text{if }\nu\ge \dfrac{3}{2}.
\end{cases}
\end{align*}
\section{Inference}\label{sec:estimation-cov}
To implement statistical inference based on the Gaussian approximation developed in Section \ref{sec:GA}, we require a consistent estimator of the network covariance matrix \(\Sigma_n\). Unlike the low-dimensional setting, however, consistency under conventional matrix norms is generally insufficient for high-dimensional inference.
Instead, we need uniform consistency under the matrix infinity norm.
Several covariance estimation procedures have been proposed for network-dependent data, including the HAC estimator \citep{KMS2021}, cluster-robust methods \citep{Leung2023}, and bootstrap procedures \citep{kojevnikov2021bootstrap}.
We focus on the network HAC estimation approach, and establish its convergence rate under the matrix infinity norm.
\subsection{Network HAC Estimator}
Recall the network covariance
\[
\Sigma_n = \mathrm{Var}(T_X/\sqrt{n})=\frac1n\sum_{i,i'\in V_n}\mathrm{Cov}(X_{n,i},X_{n,i'})\in \mathbb{R}^{p\times p}.
\]
For simplicity of exposition, we first assume
\(\mathbb E(X_{n,i})=0\) for all \(i\in V_n\).
Under Assumption \ref{ass:decay rate}, covariance contributions from pairs of nodes decay with network distance. Hence, nearby pairs carry the dominant contribution to \(\Sigma_n\), while distant pairs contribute primarily to a small tail component. This motivates a network HAC estimator that keeps local covariance terms, while down-weighting and eventually truncating covariance terms at larger network distances.
For \(s\ge0\), define the distance-specific covariance matrix
\[
\Gamma_n(s)=
\frac{1}{n}
\sum_{i\in V_n}
\sum_{i'\in \mathcal{N}_n^\partial(i,s)}
\mathbb{E}[X_{n,i}X_{n,i'}^\top] ,
\]
By construction,\[
\Sigma_n=\sum_{s\ge0}
\Gamma_n(s).
\]
This representation is analogous to the autocovariance decomposition in time-series analysis, where the graph distance \(s\) plays the role of a lag.
To implement the down-weighting of distant covariance terms, we use a kernel function $w$ satisfying the following conditions.
\begin{assumption}[Kernel]{\label{ass:kernel}}
A kernel function \(w:\mathbb R^+\to[-1,1]\) satisfies \(w(0)=1\) and \(w(x)=0\) for \(x>1\). Moreover, there exist
constants \(C_w<\infty\) and \(\kappa>0\) such that
\[
|w(x)-1|\le C_w x^\kappa,
\qquad 0\le x\le1.
\]
\end{assumption}
Let $b_n$ denote the bandwidth and write
\[
w_n(s)
=
w(s/b_n).
\]
With this notation, the network HAC estimator can be defined by
\begin{align}
\label{eq:hac}
\widehat\Sigma_n(b_n)
=&
\sum_{s\ge0}
w_n(s)\widehat\Gamma_n(s),
\qquad\text{where}
\\
\widehat\Gamma_n(s)
=&
\frac1n
\sum_{i\in V_n}
\sum_{i'\in\mathcal N_n^\partial(i,s)}
X_{n,i}X_{n,i'}^\top.
\notag
\end{align}
The bandwidth \(b_n\) controls the range of network distances incorporated into the estimator, while the kernel function $w$ controls how covariance contributions are attenuated as the distance increases.
The next theorem characterizes the deterministic approximation error induced by kernel weighting and bandwidth truncation.
\begin{theorem}{\label{lem:bias}}
If Assumptions \ref{ass:decay rate}, \ref{ass:volume}, \ref{ass:kernel} hold, then
\[
\left|
\mathbb{E}\widehat\Sigma_n(b_n)-\Sigma_n
\right|_\infty
=
O\left(b_n^{-\kappa}+b_n^{2d+1}\rho^{b_n}\right).
\]
In particular, if \(b_n\to\infty\), then
\[
\left|
\mathbb{E}\widehat\Sigma_n(b_n)-\Sigma_n
\right|_\infty=o(1).
\]
\end{theorem}
To establish consistency of the HAC estimator, it remains to control its stochastic fluctuations around its expectation. The following theorem establishes nonasymptotic concentration inequalities for the HAC estimator, both entrywise and uniformly over all matrix elements, under finite-moment conditions.
\begin{theorem}\label{lem:hac-doob-variance_q}
Suppose that Assumptions \ref{ass:decay rate}, \ref{ass:volume}, and \ref{ass:kernel} hold.
For each $(a,b)\in [p]^2$, there exists a constant $c>0$ such that for all $x>0$,
\begin{equation}\label{eq:hac-entry-tail}
\mathbb P\!\left(
n\bigl|\,\widehat\Sigma_{n,ab}(b_n)-\mathbb E\widehat\Sigma_{n,ab}(b_n)\,\bigr|\ge x
\right)
\lesssim
\frac{
n^{q/4} (N_n(b_n))^{q/2} \Psi_{n,q}^q(0)
}{x^{q/2}}.
\end{equation}
Consequently,
\begin{equation}\label{eq:hac-max-tail}
\mathbb P\!\left(
n\bigl|\widehat\Sigma_n(b_n)-\mathbb E\widehat\Sigma_n(b_n)\bigr|_{\infty}
\ge x
\right)
\lesssim\,\frac{p^2\,n^{q/4}(N_n(b_n))^{q/2}\Psi_{n,q}^q(0)}{x^{q/2}}.
\end{equation}
\end{theorem}
Under the stronger sub-Weibull moment condition, the concentration inequality in Theorem \ref{lem:hac-doob-variance_q} can be strengthened to an exponential tail bound, as formally stated below.
For convenience, we define the dependence-adjusted sub-Weibull norm
\[
\Theta_{n,\psi_\nu,j}
:=
\sup_{q\ge2}
\frac{\Theta_{n,q,j}}{q^\nu}
\qquad\text{and}\qquad
\Phi_{n,\psi_\nu}
:=
\max_{1\le j\le p}
\Theta_{n,\psi_\nu,j}.
\]
\begin{theorem}\label{lem:hac-doob-variance_subwei}
Suppose that Assumptions \ref{ass:volume}, \ref{ass:subweibull}, and \ref{ass:kernel} hold.
For each $(a,b)\in [p]^2$, there exists a constant $c>0$ such that for all $x>0$,
\begin{equation}\label{eq:hac-entry-tail}
\mathbb P\!\left(
n\bigl|\,\widehat\Sigma_{n,ab}(b_n)-\mathbb E\widehat\Sigma_{n,ab}(b_n)\,\bigr|\ge x
\right)
\lesssim
\exp\!\left(
- c \,
\frac{x^{\gamma}}
{ \big(\sqrt{n}\,N_n(b_n)\,\Phi_{n,\psi_\nu}^2\big)^{\gamma}}
\right).
\end{equation}
Consequently,
\begin{equation}\label{eq:hac-max-tail}
\mathbb P\!\left(
n\bigl|\widehat\Sigma_n(b_n)-\mathbb E\widehat\Sigma_n(b_n)\bigr|_{\infty}
\ge x
\right)
\lesssim\,
p^2 \exp\!\left(
- c \,
\frac{x^{\gamma}}
{ \big(\sqrt{n}\,N_n(b_n)\,\Phi_{n,\psi_\nu}^2\big)^{\gamma}}
\right),
\end{equation}
where $\gamma=1/(1+2\nu)$ for $\nu$ given in Assumption \ref{ass:subweibull}.
\end{theorem}
Combining the bias characterization with the concentration inequalities yields uniform consistency rates for the HAC estimator under the matrix infinity norm, as formally stated in the following corollary.
\begin{corollary}[Consistency Rate of the HAC Estimator]{\label{cor:sigmadis}}${}$\\
(i) Under the conditions of Theorems \ref{lem:bias} and \ref{lem:hac-doob-variance_q}, $|\widehat\Sigma_n(b_n)-\Sigma_n|_\infty=O_p(R_n)$ holds, where
\[R_n=p^{4/q}n^{-1/2}N_n(b_n)\Psi^2_{n,q}(0)+b_n^{-\kappa}+b_n^{2d+1}\rho^{b_n}.\]
(ii) Under the conditions of Theorems \ref{lem:bias} and \ref{lem:hac-doob-variance_subwei},
$|\widehat\Sigma_n(b_n)-\Sigma_n|_\infty=O_p(R^*_n)$ holds, where
\[R^*_n=n^{-1/2}\,N_n(b_n)\,\Phi_{n,\psi_\nu}^2(\log{p})^{1/\gamma}+b_n^{-\kappa}+b_n^{2d+1}\rho^{b_n}.\]
\end{corollary}
Thus far, we have focused on the HAC estimator \(\widehat\Sigma_n(b_n)\), which assumes that $\mathbb E(X_{n,i})=0$ for all $i\in V_n$.
Now, consider the more general case where $\mathbb E (X_{n,i}) = \mu_0$ (not necessarily 0) for all $i \in V_n$, where $\mu_0$ is unknown to a researcher.
Let the sample mean be denoted by $\bar X_n= T_X/n$.
The feasible HAC estimator in this general case is
\begin{align*}
\widetilde\Sigma_n(b_n)=&\sum_{s\ge 0}w_n(s)\widetilde\Gamma_n(s),
\qquad\text{where}
\\
\widetilde\Gamma_n(s)
=&
\frac{1}{n}
\sum_{i\in V_n}
\sum_{i'\in \mathcal{N}_n^\partial(i,s)}
(X_{n,i}-\bar X_n)(X_{n,i'}-\bar X_n)^\top.
\end{align*}
Under the conditions of Theorem \ref{lem:hac-doob-variance_q}, we have
\[|\widetilde \Sigma_n(b_n)-\widehat \Sigma_n(b_n)|_\infty= O_p\left(
N_n(b_n)\Psi^2_{n,q}(0) p^{2/q}n^{-1}
\right).\]
Likewise, under the conditions of Theorem \ref{lem:hac-doob-variance_subwei}, we have
\[
|\widetilde\Sigma_n(b_n)-\widehat\Sigma_n(b_n)|_\infty
=
O_p\left(
N_n(b_n)\Phi_{n,\psi_\nu}^2 (\log p)^{1/\gamma}n^{-1}
\right).
\]
Therefore, the residual errors $R_n$ or $R^*_n$ in Corollary \ref{cor:sigmadis} are of the same order between $|\widetilde\Sigma_n(b_n)-\Sigma_n|_\infty$ and $|\widehat\Sigma_n(b_n)-\Sigma_n|_\infty$.
We summarize the implications of these results formally as the following corollary.
\begin{corollary}[Consistency Rate of the Feasible HAC Estimator]\label{cor:sigmafea}${}$\\
(i) Under the conditions of Theorems \ref{lem:bias} and \ref{lem:hac-doob-variance_q}, $|\widetilde\Sigma_n(b_n)-\Sigma_n|_\infty=O_p(R_n)$ holds, where
\[R_n=p^{4/q}n^{-1/2}N_n(b_n)\Psi^2_{n,q}(0)+b_n^{-\kappa}+b_n^{2d+1}\rho^{b_n}.\]
(ii) Under the conditions of Theorems \ref{lem:bias} and \ref{lem:hac-doob-variance_subwei},
$|\widetilde\Sigma_n(b_n)-\Sigma_n|_\infty=O_p(R^*_n)$ holds, where
\[R^*_n=n^{-1/2}\,N_n(b_n)\,\Phi_{n,\psi_\nu}^2(\log{p})^{1/\gamma}+b_n^{-\kappa}+b_n^{2d+1}\rho^{b_n}.\]
\end{corollary}
\begin{remark}[Positive Semidefiniteness]
The network HAC estimator $\widetilde\Sigma_n(b_n)$ is not necessarily positive semidefinite in finite samples for a general network-distance kernel. This feature is shared by network HAC estimators more generally, including that of \citet{KMS2021}. When positive semidefiniteness is required for implementation, such as for Gaussian simulation, one may replace $\widetilde\Sigma_n(b_n)$ by a positive-semidefinite regularization, for example by replacing its negative eigenvalues with zero.
\end{remark}
\subsection{Simultaneous Inference}{\label{sec:inference}}
To conduct hypothesis testing and construct simultaneous confidence intervals, we need to approximate the critical value of the limiting Gaussian distribution in Theorems \ref{th:GA_q} and \ref{th:GA_sub}. Let \(\chi_\theta\) denote the \(\theta\)-quantile of
\[
|D_0^{-1}Z|_\infty,\qquad
Z\sim N(0,\Sigma_n),
\]
where \(0<\theta<1\).
If \(\Sigma_n\) were known, \(\chi_\theta\) could be computed directly by simulation.
In practice, \(\Sigma_n\) is unknown and must be replaced by the network HAC estimator
\(\widetilde\Sigma_n(b_n)\). Let $\widetilde D_0=[\mathrm{diag}(\widetilde\Sigma_n(b_n))]^{1/2}$. We estimate $\theta$-quantile $\chi_\theta$ by the conditional $\theta$-quantile $\widetilde{\chi}_\theta$ of
$
\left| \widetilde{D}_0^{-1} (\widetilde{\Sigma}_n(b_n))^{1/2} \eta \right|_\infty
$ given $(X_{n,i})_{i\in V_n}$, where $\eta \sim N(0,Id_p)$ is independent of $(X_{n,i})_{i\in V_n}$. Note that $\widetilde{\chi}_\theta$ can be computed by simulations where we use a Gaussian multiplier resampling method with the estimated network covariance matrices.
Given a significance level \(\alpha\in(0,1)\), we reject the null hypothesis
\(
H_0:\mu=\mu_0
\)
whenever
\[
\sqrt n
\left|
\widetilde D_0^{-1}
(\bar X_n-\mu_0)
\right|_\infty
>
\widetilde\chi_{1-\alpha}.
\]
The corresponding simultaneous \((1-\alpha)\)-confidence intervals for
\(\mu=(\mu_1,\ldots,\mu_p)^\top\) are
\[
\bar X_{n,j}
\pm
\frac{\widetilde\chi_{1-\alpha}}
{\sqrt n}
\widetilde\sigma_{jj}^{1/2},
\qquad
j=1,\ldots,p,
\]
where
\(
\widetilde\sigma_{jj}
=
(\widetilde\Sigma_n(b_n))_{jj}.
\)
The following corollary establishes the asymptotic validity of this inference procedure.
\begin{corollary}[Validity of Inference]{\label{cor:test}}
(i) Let conditions in Theorem \ref{th:GA_q} and Corollary \ref{cor:sigmafea}(i) hold. In addition, assume $R_n\log^2 p \to 0$ with $R_n$ in Corollary \ref{cor:sigmafea}(i). Then,
\begin{equation}\label{eq:test}
\sup_{\theta \in (0,1)}
\left|
\mathbb{P}\!\left(
\sqrt{n}\,\big|\widetilde D_0^{-1}(\bar X_n-\mu_0)\big|_\infty
\ge \widetilde\chi_{1-\theta}
\right)
- \theta
\right|
\;\to\; 0.
\end{equation}
(ii) Let conditions in Theorem \ref{th:GA_sub} and Corollary \ref{cor:sigmafea}(ii) hold. In addition, assume $R^*_n\log^2 p \to 0$ with $R^*_n$ in Corollary \ref{cor:sigmafea}(ii). Then,
we have validity of the test as \eqref{eq:test}.
\end{corollary}
\section{Simulation Studies I: Gaussian Approximation}
\label{sec:simGA}
In this section, we examine the finite-sample performance of the high-dimensional Gaussian approximation property using Monte Carlo simulations.
\subsection{Design Overview}
\label{subsec:sim_overview}
In each Monte Carlo repetition, we first generate a network $\mathcal{G}_n$, and then simulate node-level observations $\{X_{n,i}\}_{i\in V_n}$ from a network-dependent data generating process (DGP).
We compare the distribution of the standardized maximum statistic with its Gaussian benchmark. Performance of the approximation is assessed by using QQ plots.
\subsection{Data Generating Processes}
\label{subsec:sim_dgp}
We take the ring network and the regular random geometric graph introduced in Examples~\ref{ex:ring} and \ref{ex:rgg}, respectively. For the random geometric graph, we use the standard random geometric graph where node locations are generated independently from the uniform distribution on $[0,1]^2$, and choose the connection radius
\(
r_n
=
\left(\frac{\bar k}{\pi n}\right)^{1/2},
\)
where $\bar k$ is the target expected degree.
This construction preserves the scaling $r_n\asymp n^{-1/2}$ used in Example~\ref{ex:rgg}.
Conditional on each realized finite network, the network regularity conditions are satisfied for suitable realization-specific constants. In particular, neighborhood sizes satisfy a polynomial growth bound in graph distance, and the network admits a partition into approximately balanced clusters with small boundaries.
Conditional on the network realization, we generate $p$-dimensional observations $\{X_{n,i}\}$ according to
\begin{equation}
X_{n,i}
=
\sum_{s\ge 0}
\frac{\gamma^s}{N_n^{\partial}(i,s)} M_{n,s}
\left(
\sum_{i'\in N_n^{\partial}(i,s)}
\varepsilon_{n,i'}
\right),
\qquad i\in V_n,
\label{eq:sim_dgp_kms}
\end{equation}
where $N_n^{\partial}(i,s)$ is the neighborhood size defined in \eqref{eq:neighbor size}. The shocks $\{\varepsilon_{n,i}\}$ are independently generated from a standardized
Student-$t$ distribution, i.e.,
\(
\varepsilon_{n,i}\sim T_t/\sqrt{t/(t-2)},
\) where $T_t$ denotes a Student-$t$ random variable with $t$ degrees of freedom.
The matrices $\{M_{n,s}\}_{s\ge0}$ introduce cross-coordinate dependence in the
$p$-dimensional observations. They are generated once and then kept fixed across
Monte Carlo repetitions, with entries independently drawn from $N(0,1/p)$. Parameter $\gamma\in\{0.6,0.8\}$
controls the strength of network dependence.
This DGP is a linear network functional of independent node-level shocks. The
weight attached to shocks at graph distance $s$ decays exponentially through
$\gamma^s$, while the normalization by $N_n^{\partial}(i,s)$ prevents distant
shells with many nodes from mechanically dominating the variance. Hence the
design directly captures the interaction between dependence decay and network
distance growth.
Recall the notation $T_X = \sum_{i\in V_n} X_{n,i}$, and the covariance matrix of
$T_X/\sqrt n$ is given by
\begin{equation}
\Sigma_n
=
\mathrm {Var}(T_X/\sqrt n)
=
\frac{1}{n}\sum_{i,i'\in V_n}
\mathrm {Cov}(X_{n,i},X_{n,i'}).
\label{eq:Sigma-def}
\end{equation}
Equivalently, exchanging the order of summations yields the closed-form
representation
\begin{equation}
\Sigma_n
=
\frac{1}{n}
\sum_{i\in V_n}
v_{i} v_{i}^\top,
\qquad\text{where}\qquad
v_i
:=
\sum_{i'\in V_n}
\gamma^{d(i,i')} M_{n,d(i,i')}\mathbf 1\{d_n(i,i')<\infty\}.
\label{eq:Sigma-closed}
\end{equation}
The statistic of interest is the standardized maximum statistic
\begin{equation}
T_n
=
\sqrt{n}\left|
D_0^{-1} \sum_{i\in V_n} X_{n,i}
\right|_\infty,
\label{eq:sim_Tn}
\end{equation}
where $D_0=\left(\mathrm{diag}(\Sigma_n)\right)^{1/2}$. Let $Z\sim N(0,\Sigma_n)$ denote the Gaussian
benchmark, and let $T_Z=|D_0^{-1}Z|_\infty$ be its counterpart of the statistic $T_n$.
\subsection{Simulation results}
\label{subsec:sim_takeaways}
For each $(n,p,\gamma,t)$ configuration, we plot the empirical distributions of $T_n$ against the empirical distributions of $T_Z$ in 2000 Monte Carlo realizations. Simulation results are illustrated in Figures ~\ref{fig:ring_QQ_nu4_nu8} and ~\ref{fig:RGG_QQ_nu4_nu8} for the ring network and the random geometric graph, respectively. Figure \ref{fig:rgg_1x3} illustrates the random geometric graph networks.
\begin{figure}[p]
\centering
\begin{minipage}[t]{0.48\textwidth}
\centering
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n100_p200_rho0p600_nu8p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n100_p200_rho0p800_nu8p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n200_p200_rho0p600_nu8p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n200_p200_rho0p800_nu8p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n400_p200_rho0p600_nu8p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n400_p200_rho0p800_nu8p00_R2000ppng.png}
\end{subfigure}
\caption*{\footnotesize $t = 8$}
\end{minipage}
\hfill
\begin{minipage}[t]{0.48\textwidth}
\centering
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n100_p200_rho0p600_nu12p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n100_p200_rho0p800_nu12p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n200_p200_rho0p600_nu12p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n200_p200_rho0p800_nu12p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n400_p200_rho0p600_nu12p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/QQ_ring_n400_p200_rho0p800_nu12p00_R2000ppng.png}
\end{subfigure}
\caption*{\footnotesize $t = 12$}
\end{minipage}
\caption{
QQ plots for the ring network.
Left panels correspond to $t=8$ and right panels correspond to $t=12$.
Rows correspond to $n=\{100,200,400\}$ and columns correspond to
$\rho=\{0.6,0.8\}$ with $p=200$ and $R=2000$.
The horizontal axis reports the empirical distributions of $T_n$,
and the vertical axis reports the corresponding Gaussian distributions of $T_z$.
}
\label{fig:ring_QQ_nu4_nu8}
\end{figure}
\begin{figure}[p]
\centering
\begin{minipage}[t]{0.48\textwidth}
\centering
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n100_p200_rho0p600_nu8p00_dbar6p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n100_p200_rho0p800_nu8p00_dbar6p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n200_p200_rho0p600_nu8p00_dbar6p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n200_p200_rho0p800_nu8p00_dbar6p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n400_p200_rho0p600_nu8p00_dbar6p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n400_p200_rho0p800_nu8p00_dbar6p00_R2000ppng.png}
\end{subfigure}
\caption*{\footnotesize $t = 8$}
\end{minipage}
\hfill
\begin{minipage}[t]{0.48\textwidth}
\centering
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n100_p200_rho0p600_nu12p00_dbar6p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n100_p200_rho0p800_nu12p00_dbar6p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n200_p200_rho0p600_nu12p00_dbar6p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n200_p200_rho0p800_nu12p00_dbar6p00_R2000ppng.png}
\end{subfigure}
\vspace{2mm}
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n400_p200_rho0p600_nu12p00_dbar6p00_R2000ppng.png}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.48\linewidth}
\includegraphics[width=\linewidth]{figure/figure_rgg/QQ_rgg_n400_p200_rho0p800_nu12p00_dbar6p00_R2000ppng.png}
\end{subfigure}
\caption*{\footnotesize $t = 12$}
\end{minipage}
\caption{
QQ plots for the random geometric graph network.
Left panels correspond to $t=8$ and right panels correspond to $t=12$.
Rows correspond to $n=\{100,200,400\}$ and columns correspond to
$\rho=\{0.6,0.8\}$ with $p=200$ and $R=2000$.
The horizontal axis reports the empirical distributions of $T_n$,
and the vertical axis reports the corresponding Gaussian distributions of $T_z$.
}
\label{fig:RGG_QQ_nu4_nu8}
\end{figure}
\begin{figure}[htbp]
\centering
\begin{subfigure}[t]{0.32\linewidth}
\centering
\includegraphics[width=\linewidth]{figure/figure_rgg/RGG_n100_dbar6p0.png}
\caption{n=100}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.32\linewidth}
\centering
\includegraphics[width=\linewidth]{figure/figure_rgg/RGG_n200_dbar6p0.png}
\caption{n=200}
\end{subfigure}\hfill
\begin{subfigure}[t]{0.32\linewidth}
\centering
\includegraphics[width=\linewidth]{figure/figure_rgg/RGG_n400_dbar6p0.png}
\caption{n=400}
\end{subfigure}
\caption{
Random geometric graphs with $\bar k=6$.
Panels (a)–(c) correspond to increasing number of points.
}
\label{fig:rgg_1x3}
\end{figure}
The results are broadly consistent with our theoretical predictions.
First, the Gaussian approximation becomes less accurate when network dependence gets stronger, as reflected by larger values of $\gamma$.
Second, network topology matters.
The ring network generally produces a closer Gaussian approximation than the random geometric graph because the ring has a more homogeneous and sparser local structure.
In contrast, the random geometric graph introduces degree heterogeneity and more irregular shell sizes, which slow down the rate of finite-sample approximation precision.
The improvement from increasing $n$ is visible, but not dramatic, especially for
the random geometric graph.
This is consistent with the fact that approximation quality depends not only on sample size, but also on the interaction between dependence decay and network expansion.
When the graph contains more local neighbors or more irregular neighborhoods, the effective dependence accumulated over network shells can remain non-negligible in moderate samples.
Tail behavior also plays an important role.
Designs with $t=8$ display larger deviations from the Gaussian benchmark than designs with $t=12$, reflecting the effect of heavier-tailed shocks on high-dimensional maxima.
Overall, the approximation performs better under weaker dependence, sparser and more regular
network structure, and lighter-tailed innovations in accordance with our theoretical predictions.
\section{Simulation Studies II: Statistical Inference}\label{sec:HAC}
In this section, we investigate the finite-sample performance of the proposed method of inference based on the network HAC covariance estimator.
We evaluate the coverage accuracy of the HAC-based inference under different levels of network dependence and network density.
We focus on the more complicated random geometric graph network, and generate the data according to the DGP described in Section \ref{subsec:sim_dgp} with dimension $p=400$ and degrees of freedom $\nu=12$.
For covariance estimation, we employ the feasible network HAC estimator
\(\widetilde{\Sigma}_n(b_n)\), which does not rely on knowledge of the expected value $\mu_0$.
To compute the HAC estimator, we follow the implementation strategy of
\cite{KMS2021}.
Specifically, we employ the Parzen kernel
\[
w(x) =
\begin{cases}
1 - 6x^2 + 6x^3, & 0 \le x \le \tfrac{1}{2}, \\
2(1-x)^3, & \tfrac{1}{2} < x \le 1,\\
0, & \text{otherwise},
\end{cases}
\]
and choose the bandwidth according to
\[
b_n
=
c\frac{\log n}{\log(\mathrm{Avg.Deg.})}.
\]
Following the bandwidth specification suggested by \cite{KMS2021}, we allow the
constant \(c\) to vary across a range of values.
We conduct sets of Monte Carlo simulations across a variety of dependence and network configurations, with \(\rho\in\{0.2,0.3,0.4,0.5\}\) and \(\bar k\in\{2,3,4,5\}\).
For the bandwidth constant, we consider \(c\in\{1.2,1.3,1.4,1.5,1.6\}\).
Among the values examined, \(c=1.5\) delivers the most stable coverage performance and is therefore used in the results reported in Table \ref{tab:rgg_cov_c15}.
\begin{table}[tbp]
\centering
\caption{Average network statistics and simulated coverage probabilities of the
95\% HAC-based confidence bands for bandwidth constant \(c=1.5\).}
\label{tab:rgg_cov_c15}
\begin{tabular}{cccccccc}
\toprule
\(\bar k\) & \(n\) & Diam. & Avg.Deg. &
\multicolumn{4}{c}{Simulated coverage} \\
\cmidrule(lr){5-8}
& & & & \(\rho=0.2\) & \(\rho=0.3\) & \(\rho=0.4\) & \(\rho=0.5\) \\
\midrule
2.0 & 200 & 7 & 1.75 & 0.975 & 0.973 & 0.978 & 0.975 \\
2.0 & 400 & 16 & 1.84 & 0.965 & 0.960 & 0.961 & 0.963 \\
2.0 & 800 & 10 & 1.9275 & 0.961 & 0.960 & 0.953 & 0.952 \\
\addlinespace
3.0 & 200 & 17 & 2.77 & 0.977 & 0.976 & 0.976 & 0.971 \\
3.0 & 400 & 18 & 2.90 & 0.966 & 0.966 & 0.962 & 0.956 \\
3.0 & 800 & 20 & 2.85 & 0.960 & 0.955 & 0.956 & 0.948 \\
\addlinespace
4.0 & 200 & 14 & 3.64 & 0.977 & 0.976 & 0.974 & 0.968 \\
4.0 & 400 & 42 & 3.845 & 0.966 & 0.964 & 0.954 & 0.944 \\
4.0 & 800 & 39 & 3.8075 & 0.957 & 0.954 & 0.950 & 0.942 \\
\addlinespace
5.0 & 200 & 13 & 4.48 & 0.979 & 0.980 & 0.973 & 0.964 \\
5.0 & 400 & 40 & 4.715 & 0.965 & 0.962 & 0.953 & 0.940 \\
5.0 & 800 & 65 & 4.84 & 0.964 & 0.956 & 0.946 & 0.932 \\
\bottomrule
\end{tabular}
\end{table}
This table reports the simulated coverage probabilities of
the proposed HAC-based confidence band for a range of dependence strengths
and network densities. Overall, the procedure delivers accurate finite-sample
inference, with coverage probabilities concentrated around the nominal
\(95\%\) level across all configurations considered.
Several patterns emerge from the table. First, the empirical coverage remains
remarkably stable over different network densities. Even when the expected
degree increases from \(\bar k=2\) to \(\bar k=5\), the coverage distortion is
generally small.
Second, stronger dependence leads to a mild deterioration in coverage for most fixed network designs. For most combinations of
\(\bar k\)
and \(n\), the coverage probability tends to
decrease as \(\rho\) increases from $0.2$ to $0.5$, reflecting the greater
difficulty of estimating the long-run covariance matrix under stronger network
dependence. Nevertheless, the resulting coverage remains close to the nominal
level in most cases.
Third, the finite-sample performance improves with network size. For
$n=800$, coverage probabilities are typically very close to \(95\%\), whereas
for $n=200$ the procedure exhibits slight over-coverage in some designs. This
pattern is consistent with the asymptotic theory and suggests that the proposed
HAC estimator provides a reliable approximation even in moderately sized
networks.
\section{Application: Heterogeneous Spillover Effects}\label{sec:application}
In this section, we revisit the network experiments studied by \cite{paluck2016changing}, \cite{AronowSamii2017}, and \cite{Leung2022}, where they examine the effects of an anti-conflict intervention on adolescent social norms related to antagonistic behaviors such as harassment, rumor spreading, social exclusion, and bullying.
Focusing on average effects, \citet{Leung2022} found statistically insignificant spillover effects.
Since average effects aggregate potentially heterogeneous spillover effects across individuals, this finding may mask meaningful spillover effects present for particular subpopulations.
Motivated by this possibility, we investigate whether treatment spillovers exhibit systematic heterogeneity and whether significant conditional spillover effects emerge after conditioning on measures of network position and the social environment.
\subsection{Heterogeneous Spillover Effects: Estimation and Inference}
A researcher observes $(Y_i,D_i,W_i)_{i=1}^n$, where $Y_i$ is an outcome, $D_i$ is a binary treatment indicator, and $W_i$ is a covariate.
The outcome $Y_i$ indicates whether student $i$ reports wearing a wristband disseminated as
part of the program as a reward to students observed engaging in conflict-mitigating behavior, $D_i$ indicates whether the student $i$ is offered treatment, and $W_i$ captures measures of network position or the social environment.
In network settings, the outcome of unit $i$ is generally determined not only by its own treatment status but also by the treatment assignments of other units, that is, by the entire treatment vector $D=(D_i)_{i=1}^n$.
Following \citet{Leung2022}, we adopt the exposure mapping framework to provide a parsimonious characterization of both direct treatment and spillover exposures.
Let
\[
T_i=T(i,D,A)
\]
denote the exposure status of unit \(i\), taking values in the exposure space \(\mathcal{T}\subseteq\mathbb{R}\), where \(A\) denotes the observed network.
We can then index the potential outcomes by the exposure state, i.e., \(Y_i(t)\) for \(t \in \mathcal{T}\).
Our estimand of interest is the conditional average exposure effect
\begin{align*}
\tau(w;t,t')
=&
\mu(w,t)-\mu(w,t'),
\qquad\text{where}\\
\mu(w,t)
:=&
\mathbb E[Y_i(t)\mid W_i=w].
\end{align*}
\paragraph{Estimation.}
As in \citet{Leung2022}, we employ the inverse propensity score weighting (IPW) approach.
However, since our estimand $\tau(\,\cdot\,;t,t')$ is a nonparametric functional parameter rather than a scalar, we augment the conventional IPW estimator with kernel smoothing:
\begin{align*}
\hat\tau(w;t,t') =& \frac{ \sum_{i=1}^{n} K_{h}(W_i-w) Z_i(t,t') }{ \sum_{i=1}^{n} K_{h}(W_i-w) },
\qquad\text{where}
\\
Z_i(t,t') =& Y_i \left( \frac{\mathbf 1_i(t)}{\pi_i(t)} - \frac{\mathbf 1_i(t')}{\pi_i(t')} \right)
\qquad\text{and}\qquad
K_{h}(u)=h^{-1}K(u/h)
\end{align*}
with $K$ denoting the Gaussian kernel and $h$ denoting the bandwidth parameter.
\begin{comment}
\todo[inline,color=orange!50]{\textbf{Yuya: }The score process should also account for the fact that the denominator $\hat f_W(w)$ is estimated, like delta method. Did you account for that? Because of a first-order approximation, such a score process will have $\hat\tau(t,t') = \frac{1}{n}\sum_{i=1}^n \psi_i(t,t') + o_p(\sqrt{n})$, where $o_p$ here is uniform over the coordinates.}
\textcolor{red}{\paragraph{Note on the score process.}
For each evaluation point \(w\), write the kernel-IPW estimator as
\[
\hat\tau(w;t,t')
=
\frac{\hat g(w;t,t')}{\hat f_W(w)},
\]
where
\[
\hat g(w;t,t')
=
\frac1n\sum_{i=1}^n
K_{h}(W_i-w)Z_i(t,t')
\]
is the kernel-weighted sample moment of the IPW pseudo-outcome,
\[
\hat f_W(w)
=
\frac1n\sum_{i=1}^n
K_{h}(W_i-w)
\]
is the kernel density estimator of \(W_i\) at \(w\), and
\[
Z_i(t,t')
=
Y_i
\left\{
\frac{\mathbf 1_i(t)}{\pi_i(t)}
-
\frac{\mathbf 1_i(t')}{\pi_i(t')}
\right\}
\]
is the IPW pseudo-outcome associated with the exposure contrast
\((t,t')\).
Define the corresponding population quantities by
\[
g_{h}(w;t,t')
=
\mathbb E\!\left[
K_{h}(W_i-w)Z_i(t,t')
\right]
\]
and
\[
f_{h}(w)
=
\mathbb E\!\left[
K_{h}(W_i-w)
\right].
\]
Their ratio
\[
\tau_{h}(w;t,t')
=
\frac{
g_{h}(w;t,t')
}{
f_{h}(w)
}
\]
is the kernel-smoothed population target.
Since the ratio map
\[
h(g,f)=\frac{g}{f}
\]
has gradient
\[
\nabla h(g,f)
=
\begin{pmatrix}
1/f\\
-g/f^2
\end{pmatrix},
\]
the delta method gives
\[
\begin{aligned}
\hat\tau(w;t,t')
-
\tau_{h}(w;t,t')
\approx{}&
\frac{1}{f_{h}(w)}
\left\{
\hat g(w;t,t')
-
g_{h}(w;t,t')
\right\}
\\
&-
\frac{
g_{h}(w;t,t')
}{
f_{h}(w)^2
}
\left\{
\hat f_W(w)
-
f_{h}(w)
\right\}.
\end{aligned}
\]
The first term captures estimation error in the numerator, whereas the
second term captures estimation error in the kernel-density denominator.
Moreover,
\[
\hat g(w;t,t')
-
g_{h}(w;t,t')
=
\frac1n\sum_{i=1}^n
\left\{
K_{h}(W_i-w)Z_i(t,t')
-
g_{h}(w;t,t')
\right\},
\]
and
\[
\hat f_W(w)
-
f_{h}(w)
=
\frac1n\sum_{i=1}^n
\left\{
K_{h}(W_i-w)
-
f_{h}(w)
\right\}.
\]
Substituting these sample-mean representations into the delta-method
expansion and collecting the first-order contribution of observation
\(i\) yields the population influence score
\[
\begin{aligned}
\psi_i^0(w;t,t')
={}&
\frac{1}{f_{h}(w)}
\left\{
K_{h}(W_i-w)Z_i(t,t')
-
g_{h}(w;t,t')
\right\}
\\
&-
\frac{
g_{h}(w;t,t')
}{
f_{h}(w)^2
}
\left\{
K_{h}(W_i-w)
-
f_{h}(w)
\right\}.
\end{aligned}
\]
Using
\[
\tau_{h}(w;t,t')
=
\frac{
g_{h}(w;t,t')
}{
f_{h}(w)
},
\]
the population influence score simplifies to
\[
\psi_i^0(w;t,t')
=
\frac{
K_{h}(W_i-w)
}{
f_{h}(w)
}
\left\{
Z_i(t,t')
-
\tau_{h}(w;t,t')
\right\}.
\]
For feasible covariance estimation, the unknown population quantities
are replaced by their sample analogs, giving
\[
\hat\psi_i(w;t,t')
=
\frac{
K_{h}(W_i-w)
}{
\hat f_W(w)
}
\left\{
Z_i(t,t')
-
\hat\tau(w;t,t')
\right\}.
\]
Thus, placing \(\hat\tau(w;t,t')\) inside the kernel-weighted residual
is precisely the first-order correction induced by estimation of the
denominator \(\hat f_W(w)\).}
\end{comment}
To conduct statistical inference accounting for the first-order estimation effects of both the numerator and denominator in $\hat\tau(w,t,t')$, define the feasible score process
\[
\hat\psi_i(w;t,t')
=
\frac{
K_{h}(W_i-w)
}{
\hat f_W(w)
}
\left\{
Z_i(t,t')
-
\hat\tau(w;t,t')
\right\},
\]
where
$
\hat f_W(w)
=
\frac1n
\sum_{i=1}^{n}
K_{h}(W_i-w).
$
The centering term \(\hat\tau(w;t,t')\) accounts for the
first-order contribution of estimating the denominator of the
kernel ratio estimator.
We employ undersmoothing $h$ so that the smoothing bias is asymptotically
negligible relative to the stochastic error uniformly over the evaluation
points.
Let
\begin{align*}
\hat\psi_i(t,t')
&=
\Bigl(
\hat\psi_i(w_1;t,t'),
\dots,
\hat\psi_i(w_p;t,t')
\Bigr)^\top,
\\
\hat\tau(t,t')
&=
\Bigl(
\hat\tau(w_1;t,t'),
\dots,
\hat\tau(w_p;t,t')
\Bigr)^\top.
\end{align*}
The vector \(\hat\tau(t,t')\) summarizes the estimated conditional
spillover effects over the support of \(W_i\), while
\(\hat\psi_i(t,t')\) is used below for network HAC covariance
estimation and simultaneous inference.
\paragraph{Inference.}
Following the network HAC inference procedure developed in
Section~\ref{sec:inference}, we conduct simultaneous inference for the
vector of conditional spillover effects
\[
\hat\tau(t,t')
=
\Bigl(
\hat\tau(w_1;t,t'),
\dots,
\hat\tau(w_p;t,t')
\Bigr)^\top.
\]
Let
\[
\hat\psi_i(t,t')
=
\Bigl(
\hat\psi_i(w_1;t,t'),
\dots,
\hat\psi_i(w_p;t,t')
\Bigr)^\top
\]
denote the feasible score vector defined above, where each coordinate
accounts for the first-order effect of estimating the denominator of the
kernel ratio estimator. The covariance matrix of the first-order
stochastic component of
\[
\sqrt n
\Bigl(
\hat\tau(t,t')
-
\tau_{h}(t,t')
\Bigr)
\]
is estimated by the network HAC estimator
\[
\widetilde\Sigma_n(t,t')
=
\frac1n
\sum_{i=1}^{n}
\sum_{j=1}^{n}
\hat\psi_i(t,t')
\hat\psi_j(t,t')^\top
1\{\ell_A(i,j)\le b_n\},
\]
where
\[
\tau_{h}(t,t')
=
\Bigl(
\tau_{h}(w_1;t,t'),
\dots,
\tau_{h}(w_p;t,t')
\Bigr)^\top
\]
denotes the corresponding kernel-smoothed population target. We employ
the truncated kernel and the bandwidth \(b_n=2\), as recommended by
\citet{Leung2022}.
For each $\ell=1,\ldots,p$, let
$
\tilde\sigma_{\ell\ell}^2(t,t')
=
\widetilde\Sigma_{n,\ell\ell}(t,t')
$
denote the estimated variance of the \(\ell\)-th coordinate of
$
\sqrt n
\Bigl(
\hat\tau(t,t')
-
\tau(t,t')
\Bigr).
$
To obtain the simultaneous critical value, we simulate
\[
Z^*
\sim
N\!\left(
0,
\widetilde\Sigma_n(t,t')
\right)
\]
conditional on the data and compute
\[
c_{1-\alpha}
=
q_{1-\alpha}
\left(
\max_{1\le \ell\le p}
\left|
\frac{Z_\ell^*}
{\tilde\sigma_{\ell\ell}(t,t')}
\right|
\right),
\]
where \(q_{1-\alpha}(\cdot)\) denotes the conditional \((1-\alpha)\)-quantile.
The resulting \(1-\alpha\) confidence band for the conditional spillover function is
\[
\hat\tau(w_\ell;t,t')
\pm
\frac{
c_{1-\alpha}
\tilde\sigma_{\ell\ell}(t,t')
}{
\sqrt n
},
\qquad
\ell=1,\ldots,p.
\]
Under the linear-in-means model and complex contagion model considered in \cite{Leung2022}, network dependence decays exponentially with graph distance. Consequently, the score process underlying the kernel IPW estimator satisfies the dependence conditions required by the Gaussian approximation and network HAC inference results developed in Section~\ref{sec:inference}. We therefore apply the proposed simultaneous inference procedure to the conditional spillover estimator.
\subsection{Spillover Effects conditional on Socio-Demographic Homophily}
We explore heterogeneity in spillover effects along socio-demographic dimensions. For each student, we construct homophily measures based on the similarity between the student and his network neighbors in terms of race, gender, and grade.
For race, students may belong to multiple racial categories and some observations contain missing or non-applicable entries.
We define racial homophily as the proportion of observed neighbors whose race profile exactly matches that of student \(i\). Analogous measures are constructed for gender and grade homophily after excluding observations with missing values.
Figure~\ref{fig:homophily_dist} reports the distributions of the three homophily measures. Gender and grade homophily exhibit substantial mass points at one, indicating that many students are surrounded entirely by peers of the same gender or grade. Race homophily displays greater variation across individuals.
Because the three homophily measures are highly concentrated and partially overlapping, we summarize their common variation using the first principal component \(W_i\).
\begin{figure}[htbp]
\centering
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]{figure/race.png}
\caption{Race homophily}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]{figure/gender.png}
\caption{Gender homophily}
\end{subfigure}
\hfill
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width=\linewidth]{figure/grade.png}
\caption{Grade homophily}
\end{subfigure}
\caption{Distribution of homophily measures.}
\label{fig:homophily_dist}
\end{figure}
To avoid unstable kernel estimates in regions with sparse support, we restrict the evaluation range of the homophily index to \(w \in [-1,1]\). All reported spillover curves and confidence bands are constructed over this range. We then conduct inference on the heterogeneous spillover effect conditional on the social homophily index \(W_i=w\). Figure~\ref{fig:cate_homo} reports the estimated spillover curve together with the 95\% confidence band.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{socio.png}
\caption{Spillover effects conditional on social homophily score.}
\label{fig:cate_homo}
\end{figure}
The estimated spillover effect varies across the support of the homophily index, suggesting the presence of spillover effect heterogeneity. The curve is positive over most of the
reported range and reaches a local peak around \(W_i \approx -0.7\), but it does not
display a monotonic relationship with the homophily index. Although the point estimates
are positive for most evaluation points, the confidence band contains zero
over much of the support.
The confidence band excludes zero over a narrow region centered around \(W_i\approx -0.6\), providing evidence of positive spillover effects for students with these socio-demographic characteristics. Outside this region, the data do not provide sufficient evidence to distinguish the spillover effect from zero at the 95\% confidence level.
Overall, the results suggest that spillover effects are not constant across students and may depend on the socio-demographic composition of their friendship networks. However, the evidence for positive spillovers is concentrated in a localized region of the homophily index rather than being present uniformly across the population.
\section{Summary and Discussions}\label{sec:summary}
In this paper, we developed a Gaussian approximation theory for sums of high-dimensional network-dependent random vectors. Under suitable conditions on the decay of network dependence and the underlying network topology, we established Gaussian approximation results for the maximum norm and characterized admissible dimensionality under both finite-moment and sub-Weibull regimes. We also proposed a high-dimensional network HAC estimator and established its asymptotic convergence properties in high dimensions. Together, these theoretical results provide a foundation for asymptotically valid simultaneous inference for high-dimensional random vectors in the presence of network dependence.
Our simulation studies indicate that the proposed Gaussian approximation and inference procedure perform well in finite samples. We further illustrated the practical usefulness of the methodology through an empirical application examining how spillover effects vary with a continuous measure of network homophily, where the proposed procedure was used to construct confidence bands for the heterogeneous conditional spillover-effect function.
Gaussian approximation for high-dimensional vectors has proven useful in a
wide range of empirical applications. Beyond inference on the conditional
average treatment effect (CATE) function, as illustrated in our empirical
application, the proposed method is applicable to a broad class of functional
treatment-effect parameters, including continuous treatment effect functions, among others.
Beyond treatment-effect settings, high-dimensional inference methods have also
proven useful for inference on structural parameters that are partially
identified by conditional moment equalities or inequalities. Our results are
therefore expected to facilitate applications to inference in structural
network models as well.
\bibliographystyle{apalike}
\bibliography{reference}
\clearpage