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.
72,266 characters
Estimation and Inference for Peer Effects under Conditional Random Assignment
\title{Estimation and Inference for Peer Effects under Conditional Random Assignment\thanks{Companion software is available as an R package at \protect\url{https://github.com/zengying17/peergmm-r} and as a Stata command at \protect\url{https://github.com/zengying17/peergmm-stata}.}}
\author{Ying Zeng\thanks{School of Economics, Xiamen University, China. Email: \protect[email removed].}}
\maketitle
\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long\global\long
\begin{abstract}
Empirical studies of peer effects often exploit conditional random assignment to peer groups within urns. We develop a GMM framework for estimation and inference in this setting. The framework separately identifies endogenous and contextual peer effects and nests tests of random peer-group assignment as a special case. It permits unknown heteroskedasticity and corrects finite-urn bias in variance estimation. Its asymptotic theory allows the number of peer groups to grow through more urns, more groups within urns, or both. We establish the asymptotic validity of the procedures and evaluate their finite-sample performance through Monte Carlo simulations. We apply the method to study peer effects on personality among university students. For traits with positive reduced-form peer effects, the estimates indicate that positive contextual effects are partly offset by negative endogenous effects.
Keywords: Peer effects; Conditional random assignment; Randomization tests; Reflection problem; Spatial econometrics; Generalized method of moments.
\end{abstract}
\pagebreak{}
\section{Introduction}
To address the endogeneity of peer-group formation, empirical studies of peer effects often exploit conditional random assignment, whereby individuals are randomly assigned to peer groups within an urn (or selection pool). Examples include college roommates assigned within blocks defined by gender and roommate-matching characteristics \citep{sacerdote_peer_2001}, golfers assigned to competing groups within tournament-category cells \citep{guryan_peer_2009}, students assigned to classrooms within school-grade cells \citep{burke_classroom_2013,graham_identifying_2008}, and students assigned to study groups or sections within classes \citep{booij_ability_2017a,feld_understanding_2016,golsteyn_impact_2021}.\footnote{Recent peer-effects studies exploiting conditional random-assignment designs include \citet{barwick_digital_2026,bietenbeck_motivated_2025,chen_how_2024,hampole_peer_2026,huang_poverty_2025,radbruch_interview_2025,rivera_peers_2025,shan_peers_2025}.}
Despite the widespread use of these designs, two econometric gaps remain. First, existing tests of random peer-group assignment typically examine whether predetermined characteristics exhibit dependence among peers. Such tests can reveal departures from random assignment, but say little about the magnitude of peer dependence or the associated uncertainty. Second, methods for addressing the reflection problem under conditional random assignment remain limited. Beyond the i.i.d. error model of \citet{caeyers_exclusion_2024}, there are few methods for separately identifying the endogenous effects of peers' contemporaneous outcomes and the exogenous (contextual) effects of peers' characteristics. We develop a GMM framework that addresses both gaps.
Our first contribution is to treat the peer-dependence parameter $\lambda_{0}$ as an estimand rather than restricting attention to the null $H_{0}:\lambda_{0}=0$. Under random assignment, predetermined characteristics are uncorrelated among peers conditional on urn fixed effects, so random assignment corresponds to $\lambda_{0}=0$. Existing procedures differ in their implementation, assumptions on the error structure, requirements on urn-size variation, and asymptotic framework \citep{sacerdote_peer_2001,guryan_peer_2009,wang_peer_2010,stevenson_tests_2015,jochmans_testing_2023,caeyers_exclusion_2024}.\footnote{Section \ref{sec:Review of Tests} of the Online Appendix provides a detailed theoretical comparison. It also establishes asymptotic results for several existing procedures and shows that, under the null, the numerator of the homoskedastic test in \citet{jochmans_testing_2023} coincides up to normalization with the quadratic moment underlying our scalar GMM estimator.} Our GMM estimator measures the sign and magnitude of peer dependence and quantifies its uncertainty. It permits unknown heteroskedasticity and does not require variation in either peer-group size or urn size. Our asymptotic theory requires the total number of peer groups to diverge, with growth arising through more urns, more peer groups within urns, or both. This framework includes complete random assignment within a single urn.
Our second contribution is to develop estimation and inference for endogenous and exogenous (contextual) peer effects in the linear-in-means model of under conditional random assignment. Separating these effects is difficult because an individual's outcome and peers' contemporaneous outcomes are jointly determined, giving rise to the reflection problem \citep{manski_identification_1993}. Methods for separately estimating these effects in this setting are limited. \citet{caeyers_exclusion_2024} provide an approach under i.i.d. errors, while empirical studies in this setting often estimate reduced-form effects of peer characteristics instead \citep{burke_classroom_2013,feld_understanding_2016,guryan_peer_2009,shan_peers_2025,rivera_peers_2025}. Such reduced-form effects generally combine endogenous and contextual peer effects and therefore provide limited guidance for policy analysis \citep{fruehwirth_can_2013} and for assessing social multiplier effects \citep{glaeser_social_2003}. In our GMM framework, a quadratic moment implied by independence of the idiosyncratic errors identifies the endogenous peer effect, while linear moments implied by covariate exogeneity identify the contextual peer effects conditional on the endogenous effect. Related approaches to the reflection problem include the conditional QMLE of \citet{lee_identification_2007} with group fixed effects, and the GMM estimator of \citet{kuersteiner_efficient_2023} with random group effects. Relative to these approaches, our method does not require variation in group size and permits unknown individual-level heteroskedasticity.
Methodologically, our GMM framework builds on the linear and quadratic moment conditions introduced by \citet{kelejian_generalized_1998,kelejian_generalized_1999}. More directly, it builds on the GMM estimators for spatial autoregressive (SAR) models under unknown heteroskedasticity developed by \citet{kelejian_specification_2010} and \citet{lin_gmm_2010}.\footnote{\citet{kelejian_specification_2010} also allow for SAR disturbances.} Both papers maintain fixed-dimensional covariate frameworks and therefore do not directly accommodate a growing number of urn fixed effects. We instead eliminate the urn fixed effects by within-urn demeaning and adapt their GMM methods to the peer-group and urn structure. Exploiting this structure, we derive explicit identification conditions, including a covariate-free specification in which the endogenous peer-effect parameter is identified without external instruments. This specification is particularly useful for testing random assignment because the corresponding testing equation often contains no covariates. We also construct a heteroskedasticity-consistent variance estimator that corrects the finite-urn bias induced by within-urn demeaning. Our analysis thus extends spatial-econometric GMM methods to peer-group interactions with a growing number of urn fixed effects.
We complement these theoretical results with Monte Carlo simulations and an empirical application. The simulations examine the finite-sample performance under two growth patterns, an increasing number of urns and an increasing number of peer groups within urns. Bias and size distortions are small in the baseline designs at moderate sample sizes and decline as the total number of peer groups increases. We then apply the framework to the data of \citet{shan_peers_2025}, who study peer effects on personality in university study groups. Five of the six tests fail to reject random peer-group assignment within urns. For the traits with positive reduced-form peer effects, the structural estimates point to positive contextual effects that are partly offset by negative endogenous effects.
The remainder of the paper is organized as follows. Section \ref{sec:Motivating-Example} presents the model, estimator, and asymptotic theory. Section \ref{sec:Monte-Carlo-Simulations} reports Monte Carlo evidence. Section \ref{sec:Empirical-Application} presents an empirical application, and Section \ref{sec:Conclusion} concludes. The Appendix contains the main proofs, while the Online Appendix contains calculations, theoretical results for existing tests, additional simulations, and an additional empirical application.
\section{Model and Estimation}\label{sec:Motivating-Example}
\subsection{Model without Covariates}\label{subsec:Model-without-Covariates}
We first consider a model without covariates. This specification is particularly useful for testing conditional random assignment when the conditioning set consists only of urn fixed effects. In this case, the covariance structure of the innovations identifies $\lambda_{0}$ without instruments constructed from covariates. Identification of endogenous peer effects through the variance-covariance structure has also been studied by \citet{lee_identification_2007}, \citet{graham_identifying_2008}, \citet{kuersteiner_efficient_2023}, and \citet{caeyers_exclusion_2024}, although these approaches do not allow for unknown heteroskedasticity.
Suppose the sample includes $R$ urns, indexed by $r=1,\ldots,R$. Urn $r$ contains $G_{r}$ peer groups, indexed by $g=1,\ldots,G_{r}$. Peer group $g$ in urn $r$ has $m_{gr}$ members, and urn $r$ has $n_{r}=\sum^{G_{r}}_{g=1}m_{gr}$ members. The total sample size is $n=\sum^{R}_{r=1}n_{r}$, and the total number of peer groups is $G=\sum^{R}_{r=1}G_{r}$. Suppose that peers within the same group interact equally and that there is no interaction across peer groups. Then the peer-effects model is
\begin{equation}
Y_{igr}=\alpha_{r}+\lambda_{0}\bar{Y}_{(-i)gr}+\epsilon_{igr},\label{eq:scalar}
\end{equation}
Here $Y_{igr}$ is either an outcome or a characteristic of individual $i$ in peer group $g$ of urn $r$, and $\bar{Y}_{(-i)gr}=\sum_{i'\neq i}Y_{i'gr}/(m_{gr}-1)$ is the leave-one-out mean of $Y$ in the group, i.e., the average of $Y$ among $i$'s peers. The terms $\alpha_{r}$ and $\epsilon_{igr}$ denote the urn fixed effect and the idiosyncratic term. Throughout, the realized urn and peer-group structure and all matrices constructed from it are treated as nonstochastic arrays that may vary with $n$, with dependence on $n$ suppressed in the notation. When $Y$ is a predetermined characteristic, random assignment to peer groups within urns implies that $Y_{igr}$ is uncorrelated with $\bar{Y}_{(-i)gr}$ conditional on the urn fixed effects, and hence $\lambda_{0}=0$. Thus, testing $H_{0}:\lambda_{0}=0$ provides a test of conditional random assignment. When $Y$ is an outcome, the linear-in-means specification can arise as the Nash equilibrium of an interaction game with quadratic utility (see, e.g., \citealp{calvo-armengol_peer_2009,blume_identification_2010,pereda-fernandez_social_2017a}).
The matrix form of model (\ref{eq:scalar}) for group $g$ in urn $r$ is
\begin{eqnarray}
Y_{gr} & = & \alpha_{r}\mathbf{1}_{gr}+\lambda_{0}W_{gr}Y_{gr}+\epsilon_{gr},\label{eq:SAR_group}
\end{eqnarray}
where $Y_{gr}=(Y_{1gr},\ldots,Y_{m_{gr}gr})^{\prime}$, $\mathbf{1}_{gr}$ is an $m_{gr}\times1$ vector of ones, $W_{gr}=(J_{gr}-I_{gr})/(m_{gr}-1)$, $I_{gr}$ is the identity matrix of dimension $m_{gr}$, and $J_{gr}=\mathbf{1}_{gr}\mathbf{1}^{\prime}_{gr}$ is the $m_{gr}\times m_{gr}$ matrix of ones. In spatial econometrics, $W_{gr}$ is a spatial weight matrix and (\ref{eq:SAR_group}) is a Cliff--Ord type spatial model \citep{cliff_spatial_1973,cliff_spatial_1981a}. Stacking model (\ref{eq:SAR_group}) over peer groups gives the urn-level model
\begin{equation}
Y_{r}=\alpha_{r}\mathbf{1}_{r}+\lambda_{0}W_{r}Y_{r}+\epsilon_{r},\label{eq:SAR_urn}
\end{equation}
where $Y_{r}=(Y^{\prime}_{1r},\ldots,Y^{\prime}_{G_{r}r})^{\prime}$, $\epsilon_{r}=(\epsilon^{\prime}_{1r},\ldots,\epsilon^{\prime}_{G_{r}r})^{\prime}$, $W_{r}=\operatorname{diag}^{G_{r}}_{g=1}\{W_{gr}\}$, and $\mathbf{1}_{r}$ is the $n_{r}\times1$ vector of ones.
We next discuss estimation and inference for $\lambda_{0}$ in (\ref{eq:SAR_urn}), which naturally nests the test of $H_{0}:\lambda_{0}=0$. OLS is biased because of the reflection problem arising from the reciprocal relationship between $Y_{igr}$ and $\bar{Y}_{(-i)gr}$. We instead exploit the covariance structure of the innovations to identify and estimate $\lambda_{0}$. Let $\Lambda$ be the parameter space for $\lambda$.
\begin{assumption}
\label{assu:lambda}Suppose $\lambda_{0}$ is in the interior of $\Lambda$, and $\Lambda$ is a compact subset of $(-1,1)$.
\end{assumption}
Urns with only one group and groups with only one member do not contribute to estimation. We therefore restrict attention to the following case.
\begin{assumption}
\label{assu:size}For all $g,r$, $G_{r}\geqslant2$ and $m_{gr}\geqslant2$. The average peer-group size satisfies $n/G\leqslant C_{m}<\infty$ for some positive constant $C_{m}$ that does not depend on $n$.
\end{assumption}
Assumption \ref{assu:size} implies $2\leqslant n/G\leqslant C_{m}$, so the total sample size $n$ tends to infinity if and only if the total number of peer groups $G$ tends to infinity. We impose no further restrictions on the asymptotic behavior of $G_{r}$ or $n_{r}$. Each may remain bounded or tend to infinity. In particular, $G_{r}$ and $n_{r}$ may diverge for some urns while remaining bounded for others, so the framework accommodates settings with both large and small urns. For example, \citet{sacerdote_peer_2001} has twenty-five nonempty urns, but 99\% of the sample is in the sixteen largest urns. The number of urns $R$ may also remain bounded or tend to infinity. This flexibility contrasts with \citet{kelejian_specification_2010} and \citet{lin_gmm_2010}, whose frameworks keep the $R$ bounded, and with \citet{jochmans_testing_2023}, whose asymptotic theory requires $R$ to tend to infinity. The asymptotic results we establish in the Online Appendix for the procedures of \citet{guryan_peer_2009} and \citet{caeyers_exclusion_2024} likewise require $R$ to tend to infinity.
We assume that the innovations are independent while allowing for unknown heteroskedasticity.
\begin{assumption}
\label{assu:epsilon}Suppose the innovations $\epsilon_{igr}$ are totally independent, with $\operatorname{E}(\epsilon_{igr})=0$ and $\operatorname{E}(\epsilon^{2}_{igr})=\sigma^{2}_{igr}$, where $0<c_{\sigma}\leqslant\sigma^{2}_{igr}\leqslant C_{\sigma}<\infty$ for some constants $c_{\sigma}$ and $C_{\sigma}$. In addition, $\sup_{n\geqslant1}\sup_{i,g,r}\operatorname{E}|\epsilon_{igr}|^{4+c_{\epsilon}}<\infty$ for some $c_{\epsilon}>0$.
\end{assumption}
The independence assumption is common in the spatial and peer-effects literature and is maintained, for example, by \citet{kelejian_specification_2010}, \citet{lin_gmm_2010} and \citet{jochmans_testing_2023}. It is important for our covariance-based identification strategy. If dependence among the innovations is unrestricted, covariance generated by peer dependence cannot in general be distinguished from covariance already present in the innovations, so the covariance structure alone cannot identify $\lambda_{0}$.
A possible violation of the independence assumption is a common group-level shock. This concern is less relevant when testing random assignment using variables determined before group assignment, which cannot be affected by post-assignment group-level shocks. When $Y$ is an outcome, however, such shocks may be present. If valid instruments for $\bar{Y}_{(-i)gr}$ are available, however, $\lambda_{0}$ may still be identified from linear IV moments under appropriate exogeneity conditions, as discussed in the next section. When suitable instruments are unavailable, possible alternatives include extending the group-random-effects model of \citet{kuersteiner_efficient_2023} to accommodate urn fixed effects or using the group-fixed-effects approach of \citet{lee_identification_2007}, both of which rely on additional identifying restrictions.
Rather than estimating the urn fixed effects, we remove them by within-urn demeaning. Define $I^{\ast}_{r}=I_{r}-J_{r}/n_{r}$, where $J_{r}=\mathbf{1}_{r}\mathbf{1}^{\prime}_{r}$ is the $n_{r}\times n_{r}$ matrix of ones. Then $I^{\ast}_{r}\mathbf{1}_{r}=0$, and premultiplying a vector by $I^{\ast}_{r}$ subtracts its urn mean. For example, $I^{\ast}_{r}Y_{r}=Y_{r}-\bar{Y}_{r}\mathbf{1}_{r}$, where $\bar{Y}_{r}$ is the urn mean of $Y_{r}$. Premultiplying both sides of (\ref{eq:SAR_urn}) by $I^{\ast}_{r}$ gives the within-urn equation
\begin{equation}
I^{\ast}_{r}Y_{r}=\lambda_{0}I^{\ast}_{r}W_{r}Y_{r}+I^{\ast}_{r}\epsilon_{r}.\label{eq:SAR_urn_star}
\end{equation}
Let $\Omega_{gr}=\operatorname{E}\left(\epsilon_{gr}\epsilon^{\prime}_{gr}\right)$ and $\Omega_{r}=\operatorname{E}\left(\epsilon_{r}\epsilon^{\prime}_{r}\right)$ denote the variance-covariance matrices of group $g$ and urn $r$, respectively. Under Assumption \ref{assu:epsilon}, $\Omega_{r}=\operatorname{diag}^{G_{r}}_{g=1}\{\Omega_{gr}\}$ is diagonal. For any $n_{r}\times n_{r}$ matrix $A_{r}$ such that $A^{\ast}_{r}=I^{\ast}_{r}A_{r}I^{\ast}_{r}$ has zero diagonals, $\operatorname{E}\left[\left(I^{\ast}_{r}\epsilon_{r}\right)^{\prime}A_{r}\left(I^{\ast}_{r}\epsilon_{r}\right)\right]=\operatorname{tr}\left[\left(I^{\ast}_{r}A_{r}I^{\ast}_{r}\right)\Omega_{r}\right]=0$. This construction is closely related to the quadratic moment conditions in \citet{kelejian_specification_2010} and \citet{lin_gmm_2010}. A key distinction is that, after eliminating the urn fixed effects by within-urn demeaning, the zero-diagonal restriction is imposed on the transformed matrix $I^{\ast}_{r}A_{r}I^{\ast}_{r}$ rather than directly on $A_{r}$. Many matrices satisfy this condition. Following \citet{jochmans_testing_2023}, this paper uses
\[
A_{r}=W_{r}+\frac{1}{(n_{r}-1)}I_{r}.
\]
As shown in Section \ref{subsec:diag_Atilde} of the Online Appendix, $\operatorname{diag}(I^{\ast}_{r}A_{r}I^{\ast}_{r})=0$ for this choice of $A_{r}$. This special form of $A_{r}$ also allows us to construct a heteroskedasticity-consistent covariance estimator that corrects for the finite-urn bias induced by within-urn demeaning.
For any $\lambda\in\Lambda$, define $\epsilon^{\ast}_{r}(\lambda)=I^{\ast}_{r}(I_{r}-\lambda W_{r})Y_{r}$. At the true value, $\epsilon^{\ast}_{r}(\lambda_{0})=I^{\ast}_{r}\epsilon_{r}$, so the preceding zero-diagonal condition yields the quadratic moment condition $\operatorname{E} h^{q}_{r}(\lambda_{0})=0$, where $h^{q}_{r}(\lambda)=\epsilon^{\ast}_{r}(\lambda)^{\prime}A_{r}\epsilon^{\ast}_{r}(\lambda)$. Let $h^{q}_{n}(\lambda)=\frac{1}{n}\sum^{R}_{r=1}h^{q}_{r}(\lambda)=\frac{1}{n}\sum^{R}_{r=1}\epsilon^{\ast}_{r}(\lambda)^{\prime}A_{r}\epsilon^{\ast}_{r}(\lambda)$, and let $Q^{q}_{n}(\lambda)=h^{q}_{n}(\lambda)^{2}$ be the criterion function. The GMM estimator is $\hat{\lambda}=\operatorname{argmin}{}_{\lambda\in\Lambda}Q^{q}_{n}(\lambda)$.\footnote{\begin{singlespace}
Following \citet{kelejian_specification_2010} and \citet{lin_gmm_2010}, one could increase efficiency by adding multiple quadratic moments, for example, $h^{q}_{r}(\lambda)=(\epsilon^{\ast}_{r}(\lambda)^{\prime}A_{1r}\epsilon^{\ast}_{r}(\lambda),\ldots,\epsilon^{\ast}_{r}(\lambda)^{\prime}A_{Lr}\epsilon^{\ast}_{r}(\lambda))^{\prime}$ with $\operatorname{diag}(I^{\ast}_{r}A_{lr}I^{\ast}_{r})=0$ for $l=1,\ldots,L$. Such an extension would require verifying the identifying contribution of each moment and deriving a more complex variance-covariance estimator.
\end{singlespace}
}
\begin{thm}[Consistency and Asymptotic Normality]
\label{thm:Consistency_lambda}Suppose Assumptions \ref{assu:lambda} to \ref{assu:epsilon} hold. Suppose further that $\operatorname{E} h^{q}_{n}(\lambda)\rightarrow\bar{h}^{q}(\lambda)$ for every $\lambda\in\Lambda$. Then
(i) $\hat{\lambda}\xrightarrow{p}\lambda_{0}$ as $n\rightarrow\infty$.
(ii) Let $A^{\ast}_{r}=I^{\ast}_{r}A_{r}I^{\ast}_{r}$,
\begin{align*}
V_{\lambda,n} & =\frac{2}{n}\sum^{R}_{r=1}\operatorname{tr}\left(A^{\ast}_{r}\Omega_{r}A^{\ast}_{r}\Omega_{r}\right),\\
D_{\lambda,n} & =-\frac{2}{n}\sum^{R}_{r=1}\operatorname{tr}\left(A^{\ast}_{r}W_{r}(I_{r}-\lambda_{0}W_{r})^{-1}\Omega_{r}\right),
\end{align*}
and suppose $\lim_{n\rightarrow\infty}V_{\lambda,n}=\bar{V}_{\lambda}$, $\lim_{n\rightarrow\infty}D_{\lambda,n}=\bar{D}_{\lambda}$. Then $\bar{V}_{\lambda}>0$, $\bar{D}_{\lambda}<0$, and $\sqrt{n}\left(\hat{\lambda}-\lambda_{0}\right)\xrightarrow{d}N(0,\Sigma_{\lambda})$ as $n\rightarrow\infty$, where $\Sigma_{\lambda}=\bar{D}^{-2}_{\lambda}\bar{V}_{\lambda}$.
\end{thm}
The theorem follows by specializing the corresponding proofs of Theorems \ref{thm:consistency_theta} and \ref{thm:asym_theta} for the general model below to the model without covariates. The theorem provides the basis for testing random assignment through $H_{0}:\lambda_{0}=0$ when $Y$ is predetermined. Section \ref{subsec:VC-estimation} defines estimators $\hat{V}_{\lambda}$ and $\hat{D}_{\lambda}$ of $\bar{V}_{\lambda}$ and $\bar{D}_{\lambda}$, respectively, which are consistent under the conditions in Theorem \ref{thm:convergence_Sigma}. With $\hat{\Sigma}_{\lambda}=\hat{D}^{-2}_{\lambda}\hat{V}_{\lambda}$, $z_{\lambda}=\hat{\lambda}/\sqrt{\hat{\Sigma}_{\lambda}/n}\xrightarrow{d}N(0,1)$ under $H_{0}:\lambda_{0}=0$, yielding a feasible Wald test.
\subsection{General Model}
We now extend the model to include covariates, with the preceding model without covariates as a special case. In scalar form, the model is
\begin{equation}
Y_{igr}=\alpha_{r}+\lambda_{0}\bar{Y}_{(-i)gr}+X^{(1)\prime}_{igr}\beta_{1,0}+\bar{X}^{(1)\prime}_{(-i)gr}\beta_{2,0}+X^{(2)\prime}_{igr}\beta_{3,0}+\epsilon_{igr},\label{eq:general}
\end{equation}
where $X^{(1)}_{igr}$ is a $k_{1}\times1$ vector of exogenous nonstochastic covariates, which may include past outcomes, $\bar{X}^{(1)}_{(-i)gr}=\sum^{m_{gr}}_{i'\neq i}X^{(1)}_{i'gr}/(m_{gr}-1)$ is the corresponding peer average. The $k_{2}\times1$ vector $X^{(2)}_{igr}$ contains exogenous nonstochastic covariates whose peer averages are excluded and may include group-level characteristics. All other variables are defined as in (\ref{eq:scalar}). In the terminology of \citet{manski_identification_1993}, $\lambda_{0}$ represents endogenous peer effects, whereas $\beta_{2,0}$ represents exogenous, or contextual, peer effects.
Let $X_{igr}=(X^{(1)\prime}_{igr},\bar{X}^{(1)\prime}_{(-i)gr},X^{(2)\prime}_{igr})^{\prime}$ denote the $k_{x}\times1$ vector of exogenous covariates and let $\beta_{0}=(\beta^{\prime}_{1,0},\beta^{\prime}_{2,0},\beta^{\prime}_{3,0})^{\prime}$ collect the corresponding coefficients. The model can then be written compactly as
\[
Y_{igr}=\alpha_{r}+\lambda_{0}\bar{Y}_{(-i)gr}+X^{\prime}_{igr}\beta_{0}+\epsilon_{igr}.
\]
Let $X_{r}$ be the $n_{r}\times k_{x}$ matrix of all exogenous covariates. Using the notation in (\ref{eq:SAR_urn}), the matrix form of the general model for urn $r$ is
\begin{equation}
Y_{r}=\alpha_{r}\mathbf{1}_{r}+\lambda_{0}W_{r}Y_{r}+X_{r}\beta_{0}+\epsilon_{r}.\label{eq:before_demean}
\end{equation}
Premultiplying both sides of the equation by $I^{\ast}_{r}$ eliminates the urn fixed effect and gives
\begin{equation}
I^{\ast}_{r}Y_{r}=\lambda_{0}I^{\ast}_{r}W_{r}Y_{r}+I^{\ast}_{r}X_{r}\beta_{0}+I^{\ast}_{r}\epsilon_{r}.\label{eq:general_urn_demeaned}
\end{equation}
Let $\theta=(\lambda,\beta^{\prime})^{\prime}$ be the $\left(k_{x}+1\right)\times1$ parameter vector, and define
\begin{equation}
\epsilon^{\ast}_{r}(\theta)=I^{\ast}_{r}(I_{r}-\lambda W_{r})Y_{r}-I^{\ast}_{r}X_{r}\beta.\label{eq:epsstar_general}
\end{equation}
At the true parameter value, $\epsilon^{\ast}_{r}(\theta_{0})=I^{\ast}_{r}\epsilon_{r}$. The quadratic moment condition from Section \ref{subsec:Model-without-Covariates} therefore remains valid, $\operatorname{E}\left[\epsilon^{\ast}_{r}(\theta_{0})^{\prime}A_{r}\epsilon^{\ast}_{r}(\theta_{0})\right]=0$. In addition, exogeneity of the covariates implies the linear moment condition $\operatorname{E}\left[H^{\prime}_{r}\epsilon^{\ast}_{r}(\theta_{0})\right]=0$, where $H_{r}$ is a nonstochastic $n_{r}\times k_{h}$ IV matrix that contains $I^{\ast}_{r}X_{r}$ and may include additional columns constructed from $I^{\ast}_{r}W_{r}X_{r},I^{\ast}_{r}W^{2}_{r}X_{r},\ldots$, with redundant columns omitted after stacking.\footnote{When $X_{r}=[X^{(1)}_{r},W_{r}X^{(1)}_{r}]$, the distinct candidate blocks in $[I^{\ast}_{r}X_{r},I^{\ast}_{r}W_{r}X_{r}]$ are $I^{\ast}_{r}X^{(1)}_{r}$, $I^{\ast}_{r}W_{r}X^{(1)}_{r}$, and $I^{\ast}_{r}W^{2}_{r}X^{(1)}_{r}$.} We assume that $k_{h}$ and $k_{x}$ are fixed and allow either $k_{x}=k_{h}=0$ or $1\leqslant k_{x}\leqslant k_{h}<\infty$.
The moment function for urn $r$ is $h_{r}(\theta)=\left(\epsilon^{\ast}_{r}(\theta)^{\prime}A_{r}\epsilon^{\ast}_{r}(\theta),\begin{array}{c}
\epsilon^{\ast}_{r}(\theta)^{\prime}H_{r}\end{array}\right)^{\prime}$. The sum of $h_{r}(\theta)$ across urns, normalized by the total sample size, is
\begin{equation}
h_{n}(\theta)=\frac{1}{n}\sum^{R}_{r=1}h_{r}(\theta)=\frac{1}{n}\left(\begin{array}{c}
\sum^{R}_{r=1}\epsilon^{\ast}_{r}(\theta)^{\prime}A_{r}\epsilon^{\ast}_{r}(\theta)\\
\sum^{R}_{r=1}H^{\prime}_{r}\epsilon^{\ast}_{r}(\theta)
\end{array}\right).\label{eq:h_n_theta}
\end{equation}
When the model contains no covariates, $k_{x}=k_{h}=0$, $X_{r}$, $H_{r}$, and $\beta$ are omitted and $\theta=\lambda$. Hence $\epsilon^{\ast}_{r}(\theta)=\epsilon^{\ast}_{r}(\lambda)=I^{\ast}_{r}(I_{r}-\lambda W_{r})Y_{r}$, and the sample moment reduces to the scalar quadratic moment $h^{q}_{n}(\lambda)=\frac{1}{n}\sum^{R}_{r=1}\epsilon^{\ast}_{r}(\lambda)^{\prime}A_{r}\epsilon^{\ast}_{r}(\lambda)$, as in Section \ref{subsec:Model-without-Covariates}. The GMM estimator is defined as $\hat{\theta}=\operatorname{argmin}_{\theta\in\Theta}Q_{n}(\theta)$, where $Q_{n}(\theta)=h_{n}(\theta)^{\prime}\Xi_{n}h_{n}(\theta)$ and $\Xi_{n}$ is a symmetric positive-definite $(k_{h}+1)\times(k_{h}+1)$ weighting matrix. Conditions on $\Xi_{n}$ and specific choices of the weighting matrix are discussed below. When $k_{x}=k_{h}=0$, we set $\Xi_{n}=1$, so the criterion function reduces to $Q^{q}_{n}(\lambda)=h^{q}_{n}(\lambda)^{2}$, as in Section \ref{subsec:Model-without-Covariates}.
\subsection{Asymptotic Properties of the GMM Estimator}
We next establish conditions for consistency and asymptotic normality of $\hat{\theta}$. Let $\Theta$ denote the parameter space of $\theta=(\lambda,\beta^{\prime})^{\prime}$. Assumption \ref{assu:Theta} extends Assumption \ref{assu:lambda} to the model with covariates.
\begin{assumption}
\label{assu:Theta}Let $\Theta=\Lambda\times\mathcal{B}$, where $\Lambda$ is a compact subset of $(-1,1)$ and $\mathcal{B}$ is a compact subset of $\mathbb{R}^{k_{x}}$. The true value $\theta_{0}=(\lambda_{0},\beta_{0}')'$ lies in the interior of $\Theta$. When $k_{x}=0$, $\mathcal{B}=\mathbb{R}^{0}$ and $\Theta=\Lambda$.
\end{assumption}
Define the whole-sample matrices $Y=(Y^{\prime}_{1},\ldots,Y^{\prime}_{R})^{\prime}$, $X=(X^{\prime}_{1},\ldots,X^{\prime}_{R})^{\prime}$, $\epsilon=(\epsilon^{\prime}_{1},\ldots,\epsilon^{\prime}_{R})^{\prime}$, $I^{\ast}=\operatorname{diag}^{R}_{r=1}\{I^{\ast}_{r}\}$, $W=\operatorname{diag}^{R}_{r=1}\{W_{r}\}$, and $H=(H^{\prime}_{1},\ldots,H^{\prime}_{R})^{\prime}$. Stacking (\ref{eq:general_urn_demeaned}) across urns gives
\begin{equation}
I^{\ast}Y=\lambda_{0}I^{\ast}WY+I^{\ast}X\beta_{0}+I^{\ast}\epsilon.\label{eq:general_sample}
\end{equation}
Define $\tilde{X}=\left[W(I-\lambda_{0}W)^{-1}X\beta_{0},X\right]$. It follows that $I^{\ast}\tilde{X}=\operatorname{E}\left(I^{\ast}WY,I^{\ast}X\right)$, which is the expected matrix of right-hand-side variables in (\ref{eq:general_sample}).\footnote{Note that $I-\lambda_{0}W$ is invertible by (\ref{eq:I_lW_inv}), and $I^{\ast}$, $W$ and $(I-\lambda_{0}W)^{-1}$ commute by Remark \ref{rem:calculation}, (\ref{eq:general_sample}) implies $I^{\ast}Y=I^{\ast}(I-\lambda_{0}W)^{-1}\left(X\beta_{0}+\epsilon\right)$, and hence $I^{\ast}WY=I^{\ast}W(I-\lambda_{0}W)^{-1}\left(X\beta_{0}+\epsilon\right)$. Therefore, $\operatorname{E}\left(I^{\ast}WY\right)=I^{\ast}W(I-\lambda_{0}W)^{-1}X\beta_{0}$.} Identification depends on the rank of $I^{\ast}\tilde{X}$.
\begin{assumption}
\label{assu:identification}(i) If $k_{x}>0$, suppose $H$ contains $I^{\ast}X$ and possibly additional instruments. The elements of $I^{\ast}X$ and $H$ are uniformly bounded in absolute value, and the smallest eigenvalues of $X^{\prime}I^{\ast}X/n$ and $H^{\prime}I^{\ast}H/n$ are uniformly bounded below by a positive constant.
(ii) If $k_{x}>0$, $Q_{H\tilde{X}}=\lim_{n\rightarrow\infty}H^{\prime}I^{\ast}\tilde{X}/n$ has the same rank as $Q_{\tilde{X}\tilde{X}}=\lim_{n\rightarrow\infty}\left(\tilde{X}^{\prime}I^{\ast}\tilde{X}\right)/n$.
When $k_{x}=k_{h}=0$, parts (i) and (ii) are not imposed.
\end{assumption}
When $k_{x}>0$, part (i) requires $I^{\ast}X$ in (\ref{eq:general_sample}) to serve as its own instrument. Part (ii) requires $H$ to preserve the rank of $I^{\ast}\tilde{X}$. Since $\tilde{X}$ contains $X$ and part (i) implies $\operatorname{rank}(\lim_{n\rightarrow\infty}X^{\prime}I^{\ast}X/n)=k_{x}$, it follows that $k_{x}\leqslant\operatorname{rank}(Q_{\tilde{X}\tilde{X}})\leqslant k_{x}+1$.
If $\operatorname{rank}(Q_{\tilde{X}\tilde{X}})=k_{x}$, then $I^{\ast}W(I-\lambda_{0}W)^{-1}X\beta_{0}$ is asymptotically linearly dependent on $I^{\ast}X$. One such case is $\beta_{0}=0$, for which $\operatorname{E}(I^{\ast}WY)=0$. Another arises when every peer group has the same size $m$ and $X=[X^{(1)},WX^{(1)}]$, so $W(I-\lambda_{0}W)^{-1}X\beta_{0}$ lies in the span of $X$ and $\operatorname{rank}(Q_{\tilde{X}\tilde{X}})=k_{x}$. Part (ii) gives $\operatorname{rank}(Q_{H\tilde{X}})=k_{x}$ hence $H=I^{\ast}X$ suffices. The quadratic moment identifies $\lambda_{0}$, while the linear moments identify $\beta_{0}$ given $\lambda_{0}$. Variation in peer-group size is therefore unnecessary for identification.
If instead $\operatorname{rank}(Q_{\tilde{X}\tilde{X}})=k_{x}+1$, part (ii) requires $\operatorname{rank}(Q_{H\tilde{X}})=k_{x}+1$, so $H$ must contain at least one additional instrument for $I^{\ast}WY$ beyond $I^{\ast}X$. In this case, $\operatorname{E}\left(I^{\ast}WY\right)=I^{\ast}W(I-\lambda_{0}W)^{-1}X\beta_{0}$ is asymptotically linearly independent of $I^{\ast}X$. Since $I^{\ast}W(I-\lambda_{0}W)^{-1}X\beta_{0}=I^{\ast}WX\beta_{0}+\lambda_{0}I^{\ast}W^{2}X\beta_{0}+\cdots$, columns of $I^{\ast}WX$, $I^{\ast}W^{2}X$, and higher-order spatial lags provide natural candidates for the additional instruments needed to satisfy the rank condition. In this full-rank case, the proof of Theorem \ref{thm:consistency_theta} shows that the linear moment alone identifies $\theta$. Thus, with valid instruments for $I^{\ast}_{r}W_{r}Y_{r}$, the quadratic moment may be dropped and the independence assumption on the idiosyncratic terms in Assumption \ref{assu:epsilon} can be relaxed. We do not pursue this extension here and maintain Assumption \ref{assu:epsilon} throughout the paper.
\begin{assumption}
\label{assu:lim}Suppose $\Xi_{n}\xrightarrow{p}\bar{\Xi}$, where $\bar{\Xi}$ is finite and positive definite, and $\lim_{n\rightarrow\infty}\operatorname{E} h_{n}(\theta)=\bar{h}(\theta)$ for every $\theta\in\Theta$.
\end{assumption}
When $k_{x}\geqslant1$, Assumptions \ref{assu:identification} and \ref{assu:lim} ensure identification through the linear and quadratic moments. When $k_{x}=k_{h}=0$, identification follows from the quadratic moment and Assumption \ref{assu:lim}. As shown in the proof of Theorem \ref{thm:consistency_theta}, $\operatorname{E} h_{n}(\theta)$ is well defined and finite for every $n$ and $\theta\in\Theta$. Thus, with respect to the expected moment, Assumption \ref{assu:lim} only requires the existence of its pointwise limit $\bar{h}(\theta)$. The proof further establish that this convergence is uniform over $\Theta$, $\sup_{\theta\in\Theta}\left\Vert \operatorname{E} h_{n}(\theta)-\bar{h}(\theta)\right\Vert \rightarrow0$.
\begin{thm}[Consistency]
\label{thm:consistency_theta}Suppose Assumptions \ref{assu:size}--\ref{assu:lim} hold. Then $\hat{\theta}\xrightarrow{p}\theta_{0}$ as $n\rightarrow\infty$.
\end{thm}
Appendix \ref{subsec:Proof_Consistency_Theta} contains the proof.
For the asymptotic distribution of $\hat{\theta}$, let $H^{\ast}=I^{\ast}H$, $A^{\ast}=\operatorname{diag}^{R}_{r=1}\{A^{\ast}_{r}\}$, and $\Omega=\operatorname{diag}^{R}_{r=1}\{\Omega_{r}\}$, and define
\begin{align}
V_{n} & =\frac{1}{n}\left[\begin{array}{cc}
2\operatorname{tr}\left(A^{\ast}\Omega A^{\ast}\Omega\right) & 0\\
0 & H^{\ast\prime}\Omega H^{\ast}
\end{array}\right],\label{eq:Vn}\\
D_{n} & =-\frac{1}{n}\left(\begin{array}{c}
\begin{array}{cc}
2\operatorname{tr}\left(A^{\ast}W(I-\lambda_{0}W)^{-1}\Omega\right) & 0\\
H^{\prime}I^{\ast}W(I-\lambda_{0}W)^{-1}X\beta_{0} & H^{\prime}I^{\ast}X
\end{array}\end{array}\right).\label{eq:Dn}
\end{align}
When $k_{x}=0$, the zero-dimensional blocks in (\ref{eq:Vn}) and (\ref{eq:Dn}) are omitted, so $V_{n}=2\operatorname{tr}\left(A^{\ast}\Omega A^{\ast}\Omega\right)/n$ and $D_{n}=-2\operatorname{tr}\left(A^{\ast}W(I-\lambda_{0}W)^{-1}\Omega\right)/n$. As shown in the proof of Theorem \ref{thm:asym_theta}, $V_{n}=\operatorname{Var}\left(\sqrt{n}h_{n}(\theta_{0})\right)$ and $D_{n}=\operatorname{E}\left(\nabla_{\theta}h_{n}(\theta)\left|_{\theta=\theta_{0}}\right.\right)$.
\begin{thm}[Asymptotic Normality]
\label{thm:asym_theta}Suppose Assumptions \ref{assu:size}--\ref{assu:lim} hold and that $\lim_{n\rightarrow\infty}V_{n}=\bar{V}$, $\lim_{n\rightarrow\infty}D_{n}=\bar{D}$. Then $\bar{V}$ is positive definite, $\bar{D}$ has full column rank, and $\sqrt{n}(\hat{\theta}-\theta_{0})\xrightarrow{d}N(0,\Sigma_{\theta})$, where $\Sigma_{\theta}=(\bar{D}^{\prime}\bar{\Xi}\bar{D})^{-1}\bar{D}^{\prime}\bar{\Xi}\bar{V}\bar{\Xi}\bar{D}(\bar{D}^{\prime}\bar{\Xi}\bar{D})^{-1}$.
\end{thm}
The proof is provided in Appendix \ref{subsec:Proof_Normality_Theta}. When $k_{x}=0$, we set $\Xi_{n}=1$. When $k_{x}>0$, the optimal weighting matrix is $\bar{\Xi}=\bar{V}^{-1}$, yielding $\Sigma_{\theta}=(\bar{D}^{\prime}\bar{V}^{-1}\bar{D})^{-1}$. Following \citet{kuersteiner_dynamic_2020}, for the first step of the GMM estimator, we replace $\Omega$ in $V_{n}$ with $I$, define $V^{(1)}_{n}=\frac{1}{n}\operatorname{diag}\left[2\operatorname{tr}(A^{\ast2}),H^{\ast\prime}H^{\ast}\right]$, and use $\Xi^{(1)}_{n}=(V^{(1)}_{n})^{-1}$ as the weighting matrix.\footnote{Provided that $V^{(1)}_{n}\rightarrow\bar{V}^{(1)}$, the proof of Theorem \ref{thm:asym_theta} shows that $\bar{V}^{(1)}$ is positive definite. Hence $\Xi^{(1)}_{n}\xrightarrow{p}\bar{\Xi}^{(1)}=(\bar{V}^{(1)})^{-1}$, where $\bar{\Xi}^{(1)}$ is finite and positive definite, thereby satisfying the weighting-matrix condition in Assumption \ref{assu:lim}.} In the second step, $\Xi^{(2)}_{n}=\hat{V}^{-1}$, where $\hat{V}$ is a consistent estimator of $\bar{V}$ based on the first-step estimates, as discussed below.\footnote{In finite samples, $\hat{V}$ evaluated at the first-step estimate may not be positive definite. In that case, we retain the one-step GMM estimate as the final estimate. Because $\hat{V}\xrightarrow{p}\bar{V}$ and $\bar{V}$ is positive definite, this event has probability approaching zero and does not affect the asymptotic results.}
\subsection{Estimation of the Variance-Covariance Matrix}\label{subsec:VC-estimation}
Feasible estimation of the variance-covariance matrix $\Sigma_{\theta}$ in Theorem \ref{thm:asym_theta} requires estimators of $\bar{V}$ and $\bar{D}$, the limits of $V_{n}$ and $D_{n}$ in (\ref{eq:Vn}) and (\ref{eq:Dn}). The main difficulty is that within-urn demeaning makes the transformed innovations dependent, even though the original innovations are independent under Assumption \ref{assu:epsilon}. We correct for this dependence and obtain a consistent estimator of $\Sigma_{\theta}$ under the same asymptotic sequences considered above, where the total number of peer groups $G$ grows through more urns, more groups within urns, or both. Our approach differs from the urn-clustered estimator of \citet{jochmans_testing_2023}, whose validity requires $R\rightarrow\infty$ and bounded urn sizes. It also differs from the variance estimators of \citet{kelejian_specification_2010} and \citet{lin_gmm_2010}. When applied to our model with urn fixed effects, their fixed-dimensional regressor frameworks require $R$ to remain fixed and urn sizes to grow, so the dependence induced by within-urn demeaning vanishes asymptotically.
Let $\varsigma_{r}=(\sigma^{2}_{1r},\ldots,\sigma^{2}_{n_{r}r})^{\prime}$ so that $\Omega_{r}=\operatorname{diag}(\varsigma_{r})$, and let $\epsilon^{\ast}_{r}=I^{\ast}_{r}\epsilon_{r}=\epsilon_{r}-\bar{\epsilon}_{r}\mathbf{1}_{r}$. When $n_{r}$ is finite, $\operatorname{E}\left(\epsilon^{\ast}_{r}\odot\epsilon^{\ast}_{r}\right)\neq\varsigma_{r}$, where $\odot$ denotes element-wise multiplication. We therefore first define a bias-corrected estimator of $\varsigma_{r}$ by $\tilde{\varsigma}_{r}=P_{r}\left(\epsilon^{\ast}_{r}\odot\epsilon^{\ast}_{r}\right)$, where
\begin{equation}
P_{r}=\frac{n_{r}}{n_{r}-2}\left[I_{r}-\frac{1}{n_{r}(n_{r}-1)}J_{r}\right].\label{eq:P_r}
\end{equation}
The corresponding estimator of $\Omega_{r}$ is $\tilde{\Omega}_{r}=\operatorname{diag}(\tilde{\varsigma}_{r})$. Under Assumption \ref{assu:epsilon}, $\operatorname{E}(\tilde{\varsigma}_{r})=\varsigma_{r}$ and $\operatorname{E}(\tilde{\Omega}_{r})=\Omega_{r}$. To see this, note that the $j$th element of $\epsilon^{\ast}_{r}$ is $\epsilon^{\ast}_{jr}=\epsilon_{jr}-\bar{\epsilon}_{r}$, where $j=1,\ldots,n_{r}$ indexes individuals within urn $r$. By Assumption \ref{assu:epsilon},
\begin{align*}
\operatorname{E}(\epsilon^{\ast2}_{jr}) & =\operatorname{E}(\epsilon^{2}_{jr}-2\epsilon_{jr}\bar{\epsilon}_{r}+\bar{\epsilon}^{2}_{r})=(1-\frac{2}{n_{r}})\sigma^{2}_{jr}+\frac{1}{n^{2}_{r}}\sum^{n_{r}}_{j'=1}\sigma^{2}_{j'r},\\
\operatorname{E}(\epsilon^{\ast}_{r}\odot\epsilon^{\ast}_{r}) & =(1-\frac{2}{n_{r}})\varsigma_{r}+\frac{1}{n^{2}_{r}}\mathbf{1}_{r}\mathbf{1}^{\prime}_{r}\varsigma_{r}=\left(\frac{n_{r}-2}{n_{r}}I_{r}+\frac{1}{n^{2}_{r}}J_{r}\right)\varsigma_{r}.
\end{align*}
For $P_{r}$ in (\ref{eq:P_r}), it is readily verified that $P_{r}\left(\frac{n_{r}-2}{n_{r}}I_{r}+\frac{1}{n^{2}_{r}}J_{r}\right)=I_{r}$. Hence $\operatorname{E}(\tilde{\varsigma}_{r})=\operatorname{E}\left[P_{r}\left(\epsilon^{\ast}_{r}\odot\epsilon^{\ast}_{r}\right)\right]=\varsigma_{r}$ and $\operatorname{E}(\tilde{\Omega}_{r})=\Omega_{r}$. Thus, $P_{r}$ corrects the finite-urn bias induced by replacing $\epsilon_{jr}$ with its within-urn demeaned counterpart $\epsilon^{\ast}_{jr}$. This bias vanishes elementwise as $n_{r}\rightarrow\infty$, since $\operatorname{E}(\epsilon^{\ast}_{r}\odot\epsilon^{\ast}_{r})-\varsigma_{r}\rightarrow0$. Consequently, under the large-urn asymptotics of \citet{kelejian_specification_2010} and \citet{lin_gmm_2010}, $\epsilon^{\ast}_{r}\odot\epsilon^{\ast}_{r}$ can be used directly in place of $\varsigma_{r}$ in the variance estimators. By contrast, \citet{jochmans_testing_2023} avoids estimating $\varsigma_{r}$ by using urn-clustered standard errors, whose validity relies on $R\rightarrow\infty$.
Correcting the bias in $\epsilon^{\ast}_{r}\odot\epsilon^{\ast}_{r}$ is not sufficient for the $(1,1)$ element of $V_{n}$. In particular, $\operatorname{E}\left[\operatorname{tr}(A^{\ast}_{r}\tilde{\Omega}_{r}A^{\ast}_{r}\tilde{\Omega}_{r})\right]\neq\operatorname{tr}(A^{\ast}_{r}\Omega_{r}A^{\ast}_{r}\Omega_{r})$ because the elements of $\tilde{\varsigma}_{r}$ are correlated.\footnote{When $R$ is fixed and $n_{r}$ grows, the correlations among $\epsilon^{\ast2}_{jr}$ also vanish, so the variance-covariance estimators of \citet{kelejian_specification_2010} and \citet{lin_gmm_2010} remain valid.} We therefore apply an additional correction to the quadratic component of $V_{n}$. First, observe that for any $n_{r}\times n_{r}$ matrices $M_{1r}$ and $M_{2r}$,
\begin{equation}
\operatorname{tr}(M_{1r}\Omega_{r}M_{2r}\Omega_{r})=\sum^{n_{r}}_{j=1}\sum^{n_{r}}_{j'=1}M_{1r,jj'}\sigma^{2}_{j'r}M_{2r,j'j}\sigma^{2}_{jr}=\varsigma^{\prime}_{r}(M_{1r}\odot M^{\prime}_{2r})\varsigma_{r}.\label{eq:gamma_quadratic}
\end{equation}
Since $A^{\ast}_{r}=A^{\ast\prime}_{r}$,
\begin{equation}
\operatorname{tr}(A^{\ast}\Omega A^{\ast}\Omega)=\sum^{R}_{r=1}\operatorname{tr}(A^{\ast}_{r}\Omega_{r}A^{\ast}_{r}\Omega_{r})=\sum^{R}_{r=1}\varsigma^{\prime}_{r}\left(A^{\ast}_{r}\odot A^{\ast}_{r}\right)\varsigma_{r}.\label{eq:trAOAO}
\end{equation}
To estimate $\lim_{n\rightarrow\infty}\operatorname{tr}(A^{\ast}\Omega A^{\ast}\Omega)/n$, we use $\sum^{R}_{r=1}\tilde{\varsigma}^{\prime}_{r}A^{\dagger}_{r}\tilde{\varsigma}_{r}/n$, where
\begin{align}
A^{\dagger}_{r} & =\frac{(n_{r}-2)^{2}}{\left[4+(n_{r}-2)^{2}\right]}\left[A^{\ast}_{r}\odot A^{\ast}_{r}+\frac{4}{(n_{r}-2)^{2}(n_{r}-1)+4}T_{r}\right]+t_{r}(J_{r}-I_{r}),\label{eq:A_dagger}\\
T_{r} & =\left(\mathcal{D}_{r}-\frac{\operatorname{tr}(\mathcal{D}_{r})}{4\left(n_{r}-1\right)}I_{r}\right)(J_{r}-I_{r})+(J_{r}-I_{r})\left(\mathcal{D}_{r}-\frac{\operatorname{tr}(\mathcal{D}_{r})}{4\left(n_{r}-1\right)}I_{r}\right),\nonumber \\
t_{r} & =\frac{4(n_{r}-2)(3n_{r}-4)\operatorname{tr}(\mathcal{D}_{r})}{n_{r}(n_{r}-1)(n_{r}-3)\left[4+(n_{r}-2)^{2}\right]\left[(n_{r}-2)^{2}(n_{r}-1)+4\right]},\nonumber
\end{align}
and $\mathcal{D}_{r}=\operatorname{diag}^{G_{r}}_{g=1}\{\frac{(n_{r}-m_{gr})}{(m_{gr}-1)(n_{r}-1)}I_{gr}\}$. As shown in the proof of Theorem \ref{thm:convergence_Sigma}, this correction satisfies $\operatorname{E}\left(\tilde{\varsigma}^{\prime}_{r}A^{\dagger}_{r}\tilde{\varsigma}_{r}\right)=\varsigma^{\prime}_{r}\left(A^{\ast}_{r}\odot A^{\ast}_{r}\right)\varsigma_{r}$ given the special form of $A_{r}$ we choose.
To obtain feasible estimators, replace $\epsilon^{\ast}_{r}$ with the estimation residual $\hat{\epsilon}^{\ast}_{r}=\epsilon^{\ast}_{r}(\hat{\theta})$, where $\epsilon^{\ast}_{r}(\theta)$ is defined in (\ref{eq:epsstar_general}). Define $\hat{\varsigma}_{r}=P_{r}\left(\hat{\epsilon}^{\ast}_{r}\odot\hat{\epsilon}^{\ast}_{r}\right)$, with $P_{r}$ given in (\ref{eq:P_r}), and let $\hat{\Omega}_{r}=\operatorname{diag}(\hat{\varsigma}_{r})$ and $\hat{\Omega}=\operatorname{diag}^{R}_{r=1}\{\hat{\Omega}_{r}\}$. Finally, let $A^{\dagger}=\operatorname{diag}^{R}_{r=1}\{A^{\dagger}_{r}\}$, where $A^{\dagger}_{r}$ is defined in (\ref{eq:A_dagger}). Our estimators of $\bar{V}$ and $\bar{D}$ are
\begin{align}
\hat{V} & =\frac{1}{n}\left(\begin{array}{cc}
2\hat{\varsigma}^{\prime}A^{\dagger}\hat{\varsigma} & 0\\
0 & H^{\ast\prime}\hat{\Omega}H^{\ast}
\end{array}\right),\label{eq:Vhat}\\
\hat{D} & =-\frac{1}{n}\left(\begin{array}{c}
\begin{array}{cc}
2\operatorname{tr}\left(A^{\ast}W(I-\hat{\lambda}W)^{-1}\hat{\Omega}\right) & 0\\
H^{\prime}I^{\ast}W(I-\hat{\lambda}W)^{-1}X\hat{\beta} & H^{\prime}I^{\ast}X
\end{array}\end{array}\right).\label{eq:Dhat}
\end{align}
When $k_{x}=k_{h}=0$, $\hat{V}=\hat{V}_{\lambda}$ and $\hat{D}=\hat{D}_{\lambda}$ are their scalar $(1,1)$ entries. For the two-step GMM estimator described above, $\hat{V}$ evaluated at the first-step estimates provides the second-step weight $\Xi^{(2)}_{n}=\hat{V}^{-1}$. The variance-covariance matrix is then estimated using $\hat{V}$ and $\hat{D}$ evaluated at the final estimates. For a given weighting matrix $\Xi_{n}$, the resulting estimator of $\Sigma_{\theta}$ is $\hat{\Sigma}_{\theta}=(\hat{D}^{\prime}\Xi_{n}\hat{D})^{-1}\hat{D}^{\prime}\Xi_{n}\hat{V}\Xi_{n}\hat{D}(\hat{D}^{\prime}\Xi_{n}\hat{D})^{-1}$. Under optimal second-step weighting, $\hat{\Sigma}_{\theta}=(\hat{D}^{\prime}\hat{V}^{-1}\hat{D})^{-1}$.
\begin{thm}
\label{thm:convergence_Sigma}Suppose all conditions in Theorem \ref{thm:asym_theta} hold. In addition, suppose $\sup_{n\geqslant1}\sup_{i,g,r}\operatorname{E}\left|\epsilon_{igr}\right|^{8}<\infty$. Then $\hat{D}\xrightarrow{p}\bar{D}$, $\hat{V}\xrightarrow{p}\bar{V}$, and $\hat{\Sigma}_{\theta}\xrightarrow{p}\Sigma_{\theta}$ as $n\rightarrow\infty$.
\end{thm}
Appendix \ref{subsec:Proof_Consistency_Sigma} proves the theorem.
\section{Monte Carlo Simulations}\label{sec:Monte-Carlo-Simulations}
We conduct Monte Carlo simulations to examine the finite-sample properties of our estimators under two growth sequences covered by the asymptotic theory. In the first design, the number of urns is fixed at $R=2$, while the number of peer groups in each urn increases over $G_{r}\in\{10,20,50,100,200,500\}$. In the second design, each urn contains $G_{r}=2$ peer groups, while the number of urns increases over $R\in\{10,20,50,100,200,500\}$. In both designs, half of the peer groups in each urn have size $m_{gr}=2$ and half have size $m_{gr}=4$. Each urn therefore contains $n_{r}=3G_{r}$ observations, and the total sample size is $n=3G$. The growing-$R$ design holds urn size fixed at $n_{r}=6$, whereas the growing-$G_{r}$ design lets urn size diverge.
In the baseline design, data are generated according to
\[
Y_{igr}=\alpha_{r}+\lambda_{0}\bar{Y}_{(-i)gr}+X^{(1)\prime}_{igr}\beta_{1,0}+\bar{X}^{(1)}_{(-i)gr}\beta_{2,0}+\epsilon_{igr},
\]
where $\alpha_{r}$ and $X^{(1)}_{igr}$ are drawn independently from $N(0,1)$. The innovations are drawn independently from $N(0,m_{gr}+1)$, and are thus heteroskedastic in peer-group size. We set $\beta_{1,0}=\beta_{2,0}=1$ and consider $\lambda_{0}\in\{-0.4,0,0.4\}$. In each replication, the urn-level outcome vector is generated as $Y_{r}=(I_{r}-\lambda_{0}W_{r})^{-1}(\alpha_{r}\mathbf{1}_{r}+X^{(1)}_{r}\beta_{1,0}+W_{r}X^{(1)}_{r}\beta_{2,0}+\epsilon_{r})$. We use the two-step GMM procedure introduced after Theorem \ref{thm:asym_theta}. The IV matrix is $H_{r}=[I^{\ast}_{r}W^{2}_{r}X^{(1)}_{r},I^{\ast}_{r}W^{3}_{r}X^{(1)}_{r}]$. Each design uses 5,000 replications.
\begin{table}[!htbp]
\centering
\begin{threeparttable}
\caption{Monte Carlo results with $R=2$ and growing $G_r$}
\label{tab:mc-g}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{llccc@{\hspace{0.5em}}c@{\hspace{0.5em}}ccc@{\hspace{0.5em}}c@{\hspace{0.5em}}ccc}
\hline\hline
& & \multicolumn{3}{c}{$\lambda_0=-0.4$} & & \multicolumn{3}{c}{$\lambda_0=0$} & & \multicolumn{3}{c}{$\lambda_0=0.4$} \\
\cline{3-5}\cline{7-9}\cline{11-13}
$G$ & Statistic & $\lambda$ & $\beta_1$ & $\beta_2$ & & $\lambda$ & $\beta_1$ & $\beta_2$ & & $\lambda$ & $\beta_1$ & $\beta_2$ \\
\hline
20 & Bias & -0.030 & 0.004 & 0.017 & & -0.024 & 0.010 & 0.013 & & -0.019 & 0.014 & 0.027 \\
& MC SD & 0.139 & 0.287 & 0.380 & & 0.135 & 0.298 & 0.400 & & 0.103 & 0.304 & 0.419 \\
& Avg. SE & 0.121 & 0.259 & 0.339 & & 0.123 & 0.268 & 0.357 & & 0.093 & 0.277 & 0.375 \\
& Rej. freq. & 0.111 & 0.086 & 0.093 & & 0.089 & 0.086 & 0.088 & & 0.084 & 0.082 & 0.084 \\
\noalign{\smallskip}
40 & Bias & -0.011 & 0.007 & 0.008 & & -0.012 & 0.008 & 0.010 & & -0.009 & 0.011 & 0.006 \\
& MC SD & 0.093 & 0.198 & 0.259 & & 0.092 & 0.205 & 0.271 & & 0.069 & 0.213 & 0.288 \\
& Avg. SE & 0.087 & 0.187 & 0.245 & & 0.087 & 0.192 & 0.254 & & 0.065 & 0.198 & 0.265 \\
& Rej. freq. & 0.073 & 0.068 & 0.075 & & 0.070 & 0.067 & 0.068 & & 0.066 & 0.067 & 0.073 \\
\noalign{\smallskip}
100 & Bias & -0.005 & 0.003 & 0.002 & & -0.004 & 0.001 & 0.003 & & -0.003 & 0.006 & 0.003 \\
& MC SD & 0.056 & 0.123 & 0.157 & & 0.058 & 0.125 & 0.163 & & 0.042 & 0.131 & 0.172 \\
& Avg. SE & 0.056 & 0.119 & 0.156 & & 0.055 & 0.122 & 0.162 & & 0.041 & 0.126 & 0.169 \\
& Rej. freq. & 0.054 & 0.057 & 0.055 & & 0.061 & 0.053 & 0.048 & & 0.057 & 0.064 & 0.056 \\
\noalign{\smallskip}
200 & Bias & -0.002 & 0.001 & 0.001 & & -0.003 & 0.001 & 0.001 & & -0.002 & 0.001 & 0.003 \\
& MC SD & 0.040 & 0.085 & 0.109 & & 0.039 & 0.089 & 0.117 & & 0.029 & 0.089 & 0.122 \\
& Avg. SE & 0.039 & 0.085 & 0.110 & & 0.039 & 0.087 & 0.115 & & 0.029 & 0.090 & 0.119 \\
& Rej. freq. & 0.054 & 0.053 & 0.052 & & 0.053 & 0.058 & 0.056 & & 0.046 & 0.049 & 0.057 \\
\noalign{\smallskip}
400 & Bias & -0.001 & 0.000 & -0.000 & & -0.001 & 0.003 & 0.001 & & -0.001 & 0.000 & 0.001 \\
& MC SD & 0.028 & 0.060 & 0.078 & & 0.028 & 0.062 & 0.082 & & 0.021 & 0.063 & 0.087 \\
& Avg. SE & 0.028 & 0.060 & 0.078 & & 0.028 & 0.062 & 0.081 & & 0.021 & 0.063 & 0.085 \\
& Rej. freq. & 0.050 & 0.049 & 0.052 & & 0.052 & 0.054 & 0.050 & & 0.053 & 0.049 & 0.058 \\
\noalign{\smallskip}
1,000 & Bias & -0.000 & 0.000 & -0.000 & & -0.000 & 0.000 & -0.000 & & -0.000 & 0.000 & -0.000 \\
& MC SD & 0.018 & 0.038 & 0.050 & & 0.018 & 0.039 & 0.052 & & 0.013 & 0.040 & 0.054 \\
& Avg. SE & 0.018 & 0.038 & 0.050 & & 0.018 & 0.039 & 0.051 & & 0.013 & 0.040 & 0.053 \\
& Rej. freq. & 0.053 & 0.051 & 0.053 & & 0.053 & 0.053 & 0.055 & & 0.051 & 0.053 & 0.053 \\
\noalign{\hrule height 1pt}
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item[] \textit{Notes:} The first column reports the total number of peer groups, $G=\sum_r G_r$. For each value of $G$, the four rows report bias, Monte Carlo standard deviation, average estimated standard error, and the rejection frequency of the nominal 5 percent two-sided Wald test of the corresponding true parameter value. Within each urn, peer groups are split evenly between sizes 2 and 4. The design draws $\alpha_r$ and $X^{(1)}_{igr}$ independently from $N(0,1)$, sets $\beta_{1,0}=\beta_{2,0}=1$, and draws independent $\epsilon_{igr}\sim N(0,m_{gr}+1)$. Each design cell has 5,000 replications.
\end{tablenotes}
\end{threeparttable}
\end{table}
\begin{table}[!htbp]
\centering
\begin{threeparttable}
\caption{Monte Carlo results with $G_r=2$ and growing $R$}
\label{tab:mc-r}
\setlength{\tabcolsep}{3pt}
\begin{tabular}{llccc@{\hspace{0.5em}}c@{\hspace{0.5em}}ccc@{\hspace{0.5em}}c@{\hspace{0.5em}}ccc}
\hline\hline
& & \multicolumn{3}{c}{$\lambda_0=-0.4$} & & \multicolumn{3}{c}{$\lambda_0=0$} & & \multicolumn{3}{c}{$\lambda_0=0.4$} \\
\cline{3-5}\cline{7-9}\cline{11-13}
$G$ & Statistic & $\lambda$ & $\beta_1$ & $\beta_2$ & & $\lambda$ & $\beta_1$ & $\beta_2$ & & $\lambda$ & $\beta_1$ & $\beta_2$ \\
\hline
20 & Bias & -0.045 & 0.007 & 0.001 & & -0.048 & 0.019 & 0.025 & & -0.033 & 0.018 & 0.032 \\
& MC SD & 0.161 & 0.340 & 0.463 & & 0.168 & 0.343 & 0.484 & & 0.131 & 0.359 & 0.509 \\
& Avg. SE & 0.129 & 0.282 & 0.376 & & 0.136 & 0.292 & 0.393 & & 0.107 & 0.304 & 0.419 \\
& Rej. freq. & 0.147 & 0.115 & 0.129 & & 0.138 & 0.105 & 0.134 & & 0.121 & 0.102 & 0.122 \\
\noalign{\smallskip}
40 & Bias & -0.022 & 0.005 & 0.003 & & -0.020 & 0.013 & 0.010 & & -0.016 & 0.013 & 0.010 \\
& MC SD & 0.105 & 0.226 & 0.304 & & 0.111 & 0.233 & 0.323 & & 0.085 & 0.241 & 0.335 \\
& Avg. SE & 0.096 & 0.207 & 0.278 & & 0.101 & 0.213 & 0.293 & & 0.078 & 0.222 & 0.307 \\
& Rej. freq. & 0.094 & 0.085 & 0.088 & & 0.094 & 0.082 & 0.091 & & 0.079 & 0.079 & 0.085 \\
\noalign{\smallskip}
100 & Bias & -0.008 & 0.004 & 0.003 & & -0.008 & 0.005 & 0.001 & & -0.007 & 0.005 & 0.005 \\
& MC SD & 0.065 & 0.139 & 0.190 & & 0.068 & 0.141 & 0.202 & & 0.053 & 0.145 & 0.206 \\
& Avg. SE & 0.062 & 0.134 & 0.183 & & 0.066 & 0.137 & 0.190 & & 0.051 & 0.143 & 0.199 \\
& Rej. freq. & 0.069 & 0.061 & 0.068 & & 0.063 & 0.062 & 0.068 & & 0.065 & 0.059 & 0.059 \\
\noalign{\smallskip}
200 & Bias & -0.004 & 0.002 & 0.001 & & -0.004 & 0.002 & -0.000 & & -0.003 & 0.003 & 0.004 \\
& MC SD & 0.045 & 0.099 & 0.134 & & 0.048 & 0.099 & 0.140 & & 0.037 & 0.103 & 0.145 \\
& Avg. SE & 0.045 & 0.095 & 0.131 & & 0.047 & 0.098 & 0.136 & & 0.036 & 0.102 & 0.142 \\
& Rej. freq. & 0.060 & 0.062 & 0.061 & & 0.057 & 0.056 & 0.061 & & 0.059 & 0.054 & 0.056 \\
\noalign{\smallskip}
400 & Bias & -0.003 & 0.002 & 0.003 & & -0.002 & 0.001 & -0.000 & & -0.002 & 0.002 & 0.001 \\
& MC SD & 0.033 & 0.068 & 0.095 & & 0.033 & 0.071 & 0.100 & & 0.026 & 0.071 & 0.100 \\
& Avg. SE & 0.032 & 0.068 & 0.093 & & 0.033 & 0.070 & 0.097 & & 0.026 & 0.072 & 0.101 \\
& Rej. freq. & 0.063 & 0.048 & 0.058 & & 0.054 & 0.055 & 0.054 & & 0.054 & 0.050 & 0.046 \\
\noalign{\smallskip}
1,000 & Bias & -0.001 & -0.000 & 0.000 & & -0.001 & 0.000 & -0.002 & & -0.001 & 0.000 & 0.000 \\
& MC SD & 0.020 & 0.044 & 0.060 & & 0.021 & 0.045 & 0.062 & & 0.016 & 0.046 & 0.064 \\
& Avg. SE & 0.020 & 0.043 & 0.059 & & 0.021 & 0.044 & 0.061 & & 0.016 & 0.046 & 0.064 \\
& Rej. freq. & 0.051 & 0.052 & 0.055 & & 0.057 & 0.056 & 0.054 & & 0.049 & 0.052 & 0.049 \\
\noalign{\hrule height 1pt}
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item[] \textit{Notes:} The first column reports the total number of peer groups, $G=\sum_r G_r$. For each value of $G$, the four rows report bias, Monte Carlo standard deviation, average estimated standard error, and the rejection frequency of the nominal 5 percent two-sided Wald test of the corresponding true parameter value. Within each urn, peer groups are split evenly between sizes 2 and 4. The design draws $\alpha_r$ and $X^{(1)}_{igr}$ independently from $N(0,1)$, sets $\beta_{1,0}=\beta_{2,0}=1$, and draws independent $\epsilon_{igr}\sim N(0,m_{gr}+1)$. Each design cell has 5,000 replications.
\end{tablenotes}
\end{threeparttable}
\end{table}
Table \ref{tab:mc-g} reports results for the designs with $R=2$ and growing $G_{r}$, while Table \ref{tab:mc-r} reports results for the designs with $G_{r}=2$ and growing $R$. The first column of each table reports the total number of groups, $G=\sum_{r}G_{r}$. For each value of $G$, four successive rows report bias, Monte Carlo standard deviation, average estimated standard error, and the rejection frequency of the nominal 5 percent two-sided Wald test.
The estimators perform well under both growth sequences, with bias and Monte Carlo standard deviations generally declining as $G$ increases. At $G=100$, the absolute biases of the estimators of $\lambda$, $\beta_{1}$, and $\beta_{2}$ are below 0.01 in both designs. The bias of $\hat{\lambda}$ is somewhat smaller when $G_{r}$ grows than when $R$ grows at smaller values of $G$, while the differences for $\hat{\beta}_{1}$ and $\hat{\beta}_{2}$ are less systematic.
Inference is less accurate at smaller values of $G$. The average estimated standard errors understate the corresponding Monte Carlo standard deviations, particularly when $G_{r}=2$ and $R$ grows, leading to some over-rejection by the Wald tests. For the test of $\lambda=\lambda_{0}$, rejection frequencies across the two designs and the three values of $\lambda_{0}$ range from 8.4 to 14.7 percent at $G=20$. The estimated standard errors become close to the Monte Carlo standard deviations as $G$ increases, and the corresponding rejection frequencies range from 4.6 to 6.0 percent at $G=200$.
The Online Appendix reports additional simulation results for designs without covariates, with non-normal innovations, with constant peer-group sizes and homoskedastic innovations, and with larger peer groups. Across these designs, bias decreases and rejection frequencies approach the nominal level as $G$ increases. At smaller values of $G$, however, inference is more distorted under non-normal innovations, particularly when $G_{r}=2$ and $R$ increases.
\section{Empirical Application}\label{sec:Empirical-Application}
We apply our GMM estimator to the field experiment of \citet{shan_peers_2025}, which studies peer effects in personality development among university students.
The experiment covers six cohorts of students enrolled in an introductory economics course from the 2018/19 through 2023/24 academic years. Within each cohort, students were randomly assigned to four-person study groups within three broad study-program strata. For the 2020/21 cohort, assignment was additionally stratified by the last digit of each student's ID. We treat each cohort-by-stratum cell as an urn. The baseline data contain 1,776 students in 444 study groups and 45 urns. The 30 urns from 2020/21 contain between 2 and 34 students, whereas the 15 urns from the other five cohorts contain between 12 and 283 students. We exclude 23 groups whose members span more than one urn and eight urns with only one remaining group, since such urns do not contribute to estimation. The baseline estimation sample then contains 1,652 students in 413 groups and 35 urns. It combines small and large urns, as permitted by our asymptotic framework, while all peer groups have the same size, illustrating that identification does not require variation in group size. Before study-group assignment, the experiment measured six personality traits---competitiveness, openness, conscientiousness, extraversion, agreeableness, and neuroticism---together with gender, age, major, course-retaking status, and high-school characteristics. These baseline variables are observed for all 1,776 students. At the end of the semester (endline), 1,229 students, or 69 percent of the baseline sample, reported all six personality traits. \citet{shan_peers_2025} provide further details on the experimental design and data.
\begin{table}[!htbp]
\centering
\begin{threeparttable}
\setlength{\tabcolsep}{3pt}
\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}
\caption{Randomization test}
\label{tab:sz-randomization}
\begin{tabular}{@{}l*{6}{c}@{}}
\hline\hline
&\multicolumn{1}{c}{(1)}&\multicolumn{1}{c}{(2)}&\multicolumn{1}{c}{(3)}&\multicolumn{1}{c}{(4)}&\multicolumn{1}{c}{(5)}&\multicolumn{1}{c}{(6)}\\
&\multicolumn{1}{c}{\shortstack{Competi-\\tiveness}}&\multicolumn{1}{c}{\shortstack{Open-\\ness}}&\multicolumn{1}{c}{\shortstack{Conscien-\\tiousness}}&\multicolumn{1}{c}{\shortstack{Extra-\\version}}&\multicolumn{1}{c}{\shortstack{Agreeable-\\ness}}&\multicolumn{1}{c}{\shortstack{Neuroti-\\cism}}\\
\hline
$\hat{\lambda}$& 0.009 & -0.024 & -0.063\sym{**} & -0.006 & -0.027 & 0.033 \\
& (0.032) & (0.032) & (0.031) & (0.032) & (0.030) & (0.031) \\
\hline
Observations& 1,652 & 1,652 & 1,652 & 1,652 & 1,652 & 1,652 \\
Urns & 35 & 35 & 35 & 35 & 35 & 35 \\
Peer groups & 413 & 413 & 413 & 413 & 413 & 413 \\
\noalign{\hrule height 1pt}
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item[] \textit{Notes:} Each column reports the scalar GMM estimate $\hat{\lambda}$ for the indicated baseline personality trait. The sample excludes peer groups spanning multiple urns and urns containing fewer than two remaining peer groups. Standard errors are in parentheses below the estimates. \sym{*} \(p<0.10\), \sym{**} \(p<0.05\), \sym{***} \(p<0.01\).
\end{tablenotes}
\end{threeparttable}
\end{table}
Table 2 of \citet{shan_peers_2025} reports the randomization test proposed by \citet{guryan_peer_2009} for the six baseline personality traits. We revisit the same randomization implication using the scalar GMM estimator developed in Section \ref{subsec:Model-without-Covariates}. Table \ref{tab:sz-randomization} reports the estimates. For five of the six traits, $\hat{\lambda}$ ranges from $-0.027$ to $0.033$ and is not statistically different from zero at the 5 percent level. For conscientiousness, $\hat{\lambda}=-0.063$, and $H_{0}:\lambda_{0}=0$ is rejected at the 5 percent level. Because the table tests six related outcomes, however, one rejection at this level can arise by chance and does not by itself constitute strong evidence against random assignment.\footnote{\citet{jochmans_testing_2023} provides a multivariate test of random assignment. Our GMM framework could likewise be extended to multiple variables by stacking the variable-specific quadratic moments and estimating their full covariance matrix. We leave the formal development of this extension to future research.}
Table 3 of \citet{shan_peers_2025} estimates the effect of baseline peer personality on endline personality using a reduced-form specification that omits contemporaneous peer personality. As discussed above, the resulting coefficient on baseline peer personality generally combines endogenous and contextual peer effects. Specifically, when all peer groups have the same size $m$, as in this application, the structural model in (\ref{eq:general}) implies that the reduced-form coefficient on $\bar{X}^{(1)}_{(-i)gr}$ is\footnote{To see this, the reduced form corresponding to (\ref{eq:general}), with $X^{(2)}$ omitted, is
\begin{align*}
Y_{gr} & =(I_{gr}-\lambda_{0}W_{gr})^{-1}\left(\alpha_{r}\mathbf{1}_{gr}+X^{(1)}_{gr}\beta_{1,0}+W_{gr}X^{(1)}_{gr}\beta_{2,0}+\epsilon_{gr}\right)\\
& =\frac{\alpha_{r}}{1-\lambda_{0}}\mathbf{1}_{gr}+X^{(1)}_{gr}\widetilde{\beta}_{1,0}+W_{gr}X^{(1)}_{gr}\widetilde{\beta}_{2,0}+(I_{gr}-\lambda_{0}W_{gr})^{-1}\epsilon_{gr},
\end{align*}
where $\widetilde{\beta}_{1,0}=\frac{\left[(m-1)-(m-2)\lambda_{0}\right]\beta_{1,0}+\lambda_{0}\beta_{2,0}}{(1-\lambda_{0})(m-1+\lambda_{0})}$, and $\widetilde{\beta}_{2,0}=\frac{(m-1)(\beta_{2,0}+\lambda_{0}\beta_{1,0})}{(1-\lambda_{0})(m-1+\lambda_{0})}$ is the coefficients on the peer mean $W_{gr}X^{(1)}_{gr}$.}
\[
\widetilde{\beta}_{2,0}=\frac{(m-1)(\beta_{2,0}+\lambda_{0}\beta_{1,0})}{(1-\lambda_{0})(m-1+\lambda_{0})}.
\]
To separate the two peer effects, we estimate the structural model for each endline personality trait. Each specification includes the following variables and their peer averages: (i) the six baseline personality traits; (ii) high-school characteristics, including math and language grades, study hours, and an indicator for German as the language of instruction; and (iii) gender, course-retaking status, age fixed effects, and major fixed effects. This rich set of predetermined own and peer controls absorbs systematic variation associated with observed characteristics, making the mutual independence of the remaining idiosyncratic innovations in Assumption \ref{assu:epsilon} more plausible. For comparability with \citet{shan_peers_2025}, we standardize the baseline and endline personality traits.\footnote{Standardizing either the outcome or a covariate also standardizes its peer average. In both cases, $\lambda_{0}$ is unchanged and only the relevant $\beta$ coefficients are rescaled.} Constructing peer-average outcomes requires complete endline data for all members of a group. We therefore retain only groups for which endline personality is observed for every member. After further removing urns that contain only one remaining group, the estimation sample comprises 344 students in 86 groups and 13 urns.
\begin{table}[!htbp]
\centering
\begin{threeparttable}
\setlength{\tabcolsep}{2.5pt}
\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}
\caption{Peer effects on endline personality}
\label{tab:sz-personality}
\begin{tabular}{@{}l*{6}{c}@{}}
\hline\hline
&\multicolumn{1}{c}{(1)}&\multicolumn{1}{c}{(2)}&\multicolumn{1}{c}{(3)}&\multicolumn{1}{c}{(4)}&\multicolumn{1}{c}{(5)}&\multicolumn{1}{c}{(6)}\\
&\multicolumn{1}{c}{\shortstack{Competi-\\tiveness}}&\multicolumn{1}{c}{\shortstack{Open-\\ness}}&\multicolumn{1}{c}{\shortstack{Conscien-\\tiousness}}&\multicolumn{1}{c}{\shortstack{Extra-\\version}}&\multicolumn{1}{c}{\shortstack{Agreeable-\\ness}}&\multicolumn{1}{c}{\shortstack{Neuroti-\\cism}}\\
\hline
$\hat{\lambda}$& -0.160\sym{**} & -0.201\sym{**} & -0.137\sym{*} & -0.037 & -0.076 & -0.085 \\
& (0.072) & (0.084) & (0.072) & (0.075) & (0.078) & (0.072) \\
Peer competitiveness& \textbf{0.210\sym{**}} & 0.128 & 0.044 & -0.076 & -0.094 & -0.054 \\
& (0.097) & (0.085) & (0.070) & (0.063) & (0.086) & (0.078) \\
Peer openness& -0.149 & \textbf{0.236\sym{**}} & 0.104 & 0.073 & 0.167\sym{*} & -0.057 \\
& (0.092) & (0.104) & (0.083) & (0.083) & (0.089) & (0.091) \\
Peer conscientiousness& -0.112 & -0.221\sym{**} & \textbf{0.149} & 0.158\sym{*} & -0.020 & -0.093 \\
& (0.099) & (0.093) & (0.093) & (0.088) & (0.105) & (0.097) \\
Peer extraversion& 0.163\sym{**} & -0.005 & -0.027 & \textbf{0.018} & -0.040 & -0.005 \\
& (0.083) & (0.085) & (0.090) & (0.098) & (0.101) & (0.083) \\
Peer agreeableness& -0.068 & 0.065 & 0.005 & 0.020 & \textbf{0.047} & -0.037 \\
& (0.079) & (0.078) & (0.080) & (0.071) & (0.096) & (0.080) \\
Peer neuroticism& 0.129 & 0.172\sym{**} & 0.187\sym{**} & 0.042 & 0.026 & \textbf{-0.109} \\
& (0.081) & (0.085) & (0.083) & (0.074) & (0.096) & (0.099) \\
\hline
Implied reduced-form& 0.098 & 0.081 & 0.055 & -0.009 & -0.003 &-0.156\sym{*} \\
coefficient on peer mean& (0.076) & (0.071) & (0.071) & (0.078) & (0.077) & (0.086) \\
Observations& 344 & 344 & 344 & 344 & 344 & 344 \\
Urns & 13 & 13 & 13 & 13 & 13 & 13 \\
Peer groups & 86 & 86 & 86 & 86 & 86 & 86 \\
\noalign{\hrule height 1pt}
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item[] \textit{Notes:} Each column reports a separate specification for the indicated standardized endline personality trait. The $\hat{\lambda}$ row reports the estimated endogenous peer effect, while the peer-personality rows report estimated exogenous peer effects. Boldface marks the exogenous peer-effect estimate for the baseline trait corresponding to the column's endline outcome. The implied reduced-form rows report the coefficient on the peer mean of that baseline trait and its delta-method standard error. The sample excludes peer groups spanning multiple urns, retains only peer groups in which every member has all six endline personality traits observed, and excludes urns containing fewer than two remaining peer groups. Baseline and endline personality traits are standardized. All specifications include the six baseline personality traits, high-school characteristics, gender, course-retaking status, age fixed effects, and major fixed effects, together with their peer averages. Standard errors are in parentheses below the estimates. \sym{*} \(p<0.10\), \sym{**} \(p<0.05\), \sym{***} \(p<0.01\).
\end{tablenotes}
\end{threeparttable}
\end{table}
Table \ref{tab:sz-personality} reports the estimated endogenous peer effect and the six contextual effects of baseline peer personality. In each column, the outcome is a standardized endline personality trait. The bold entry is the coefficient on the peer average of the baseline trait corresponding to the column's endline outcome, and the bottom rows report the corresponding implied reduced-form coefficient $\widetilde{\beta}_{2,0}$ and its delta-method standard error.
\citet{shan_peers_2025} find positive and statistically significant reduced-form peer effects for competitiveness, openness, and conscientiousness, but not for the other three traits. Our structural estimates show how these reduced-form effects decompose into endogenous and contextual components. For competitiveness and openness, the endogenous-effect estimates are $-0.160$ and $-0.201$, while the corresponding same-trait contextual-effect estimates are $0.210$ and $0.236$. All four estimates are statistically significant at the 5 percent level. For conscientiousness, the endogenous-effect estimate is $-0.137$ and significant at the 10 percent level, whereas the contextual-effect estimate is $0.149$ and not statistically significant. The implied reduced-form coefficients for competitiveness, openness, and conscientiousness are $0.098$, $0.081$, and $0.055$, respectively. These point estimates are close to the corresponding fully controlled estimates of $0.077$, $0.067$, and $0.058$ reported by \citet{shan_peers_2025}, although they are less precisely estimated and not statistically significant in our smaller complete-group sample.
Our structural estimates reveal a pattern that is not apparent from the reduced-form evidence alone. For competitiveness, openness, and conscientiousness, positive same-trait contextual effects are partly offset by negative endogenous effects, yielding smaller positive implied reduced-form coefficients. Under an interaction-game interpretation of the linear-in-means model, this pattern is consistent with peer characteristics acting as complementary environmental inputs and contemporaneous peer outcomes acting as strategic substitutes within the group.
\section{Conclusion}\label{sec:Conclusion}
Existing work under conditional random assignment has largely focused on testing the assignment mechanism, while empirical applications often rely on reduced-form specifications that do not distinguish endogenous from contextual peer effects. We develop a GMM framework that separately identifies these two effects, with tests of random peer-group assignment arising as a special case of inference.
Our approach builds on spatial-econometric GMM methods. It combines a quadratic moment implied by the covariance structure of independent, heteroskedastic innovations with linear instrumental-variable moments. The asymptotic framework allows the total number of peer groups to diverge through an increasing number of urns, an increasing number of groups within urns, or both. We also construct a heteroskedasticity-consistent covariance estimator that corrects the finite-urn bias induced by within-urn demeaning. The Monte Carlo simulations show good finite-sample performance at moderate sample sizes, while an application to personality among university students illustrates how the proposed framework can uncover patterns that are obscured by reduced-form estimates.
The analysis focuses on linear-in-means interactions within peer groups. Extending the framework to more general forms of social-network interaction is a natural direction for future research.
\bibliographystyle{elsarticle-harv}
\bibliography{random_peer}
\newpage{}