Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
100,436 characters · 9 sections · 87 citation commands
Quasi-Score Matching Estimation for Spatial Autoregressive Model with Random Weights Matrix and Regressors
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 {\spacingset{1.5} } \fi
\if10 {
} \fi
{\it Keywords:} Computational scalability, random weights matrix, random-X regression, score matching, spatial autoregressive model.
\spacingset{1.8}
Spatial econometrics focuses on modeling spatial cross-sectional dependence, with the spatial autoregressive (SAR) model, introduced by cliff1973spatial, being one of the most popular and widely used in the field (see, e.g., anselin1988spatial and lesage2009introduction). Its clear structure and meaningful interpretation have made it a key tool in both research and real-world applications.
The SAR model is defined for spatial locations $i=1,\cdots,n$, where $y_i$ is the response variable and $\boldsymbol{x}_i$ is the associated $p$-dimensional covariate vector. Let $\boldsymbol{y}=(y_1,\cdots,y_n)^\top$ be the response vector and $\boldsymbol{X}=(\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n)^\top$ be the $n\times p$ covariate matrix. The model is given by
where $\boldsymbol{\varepsilon}=(\varepsilon_1,\cdots,\varepsilon_n)^\top$ is the error vector with mean $\boldsymbol{0}_n=(0,\cdots,0)^\top\in\mathbb{R}^n$ and covariance matrix $\sigma^2 \boldsymbol{I}_n$, $\boldsymbol{I}_n$ is the $n$-dimensional identity matrix, $\boldsymbol{W}=(w_{ij})_{n\times n}$ is the $n\times n$ spatial weights matrix, $\lambda\in\mathbb{R}$ measures spatial dependence, and $\boldsymbol{\beta}=(\beta_1,\cdots,\beta_p)^\top\in\mathbb{R}^p$ is the coefficient vector. The parameter of interest is $\boldsymbol{\theta}=(\lambda, \boldsymbol{\beta}^\top,\sigma^2)^\top$.
To estimate $\boldsymbol{\theta}$, the quasi-maximum likelihood estimation (QMLE) is widely used in spatial econometrics lee2004asymptotic, lesage2009introduction. With growing data sizes, SAR models are increasingly applied to large $n$-dimensional response vector $\boldsymbol{y}$ and weights matrix $\boldsymbol{W}$ huang2019least, huang2020two, zhu2020multivariate, where $n$ can exceed thousands. However, computing the QMLE becomes challenging due to the $O(n^3)$ complexity of calculating the determinant $\det(\boldsymbol{I}_n-\lambda \boldsymbol{W})$ in the quasi-log-likelihood; see the discussion in Section (ref). {\color{black} To mitigate this, empirical researchers often consider alternative methods such as the two-stage least squares (2SLS, kelejian1998generalized) estimation.} Another approach introduces the matrix exponential SAR model lesage2007matrix in contrast to model ((ref)), to bypass determinant computation.
This article introduces a novel quasi-score matching estimation for the SAR model ((ref)), inspired by score matching in machine learning (see, e.g., hyvarinen2005estimation). Unlike existing score matching methods for independent and identically distributed data, to our understanding, our approach is the first to address spatially dependent data. It estimates $\boldsymbol{\theta}$ without computing the determinant, thereby largely reducing computational cost for large $n$. {\color{black} It is worth noting that although the quasi-score matching is established based on the likelihood (see Section (ref)), it does not maximize the likelihood function. Accordingly, the quasi-score matching reduces computational cost while trading off some estimation efficiency compared to the QMLE; see Corollary (ref) in Section (ref) of the supplementary material. However, based on the numerical estimation efficiencies demonstrated in various simulations in Section (ref), we find that the quasi-score matching achieves efficiency close to the QMLE, especially for large $n$. Additionally, as shown in Table (ref) of our case study, its standard errors closely match those of the QMLE.} Apart from the efficiency comparisons, our simulation results in Section (ref) highlight its significant computational advantage, requiring only about 1/16,000 to 1/26,000 of the QMLE’s time for $n=$ 10,000. {\color{black} Additional simulation studies in Section (ref) of the supplementary material compare our proposed quasi-score matching method with the generalized method of moments (GMM) estimation by lee2007gmm, the least squares estimation (LSE) by huang2019least, the 2SLS estimation by kelejian1998generalized, the best 2SLS estimation by lee2003best and the one-step estimation by gupta2023efficient, and shows that our method is computationally faster than all of these methods, especially for large $n$.}
In addition to the computational gains of the proposed quasi-score matching estimation for model ((ref)), this paper also pioneers the study of its asymptotic properties when both the spatial weights matrix $\boldsymbol{W}$ and regressors $\boldsymbol{X}$ are random. Traditionally, $\boldsymbol{W}=(w_{ij})_{n\times n}$ is assumed nonstochastic, with weights $w_{ij}$ based on fixed geographical distances $d_{ij}$ or neighborhood structures for $i,j=1,\cdots,n$ (see, e.g., kelejian1998generalized, lee2004asymptotic, lesage2009introduction). However, recent applications define the distance $d_{ij}$ using economic variables, such as GDP rankings (yu2016strategic), or set $\boldsymbol{W}$ as the normalized network adjacency matrix in social networks huang2019least,huang2020two,zhu2020multivariate. Since economic variables and network adjacency matrices are observations or realizations of random variables, assuming $\boldsymbol{W}$ to be stochastic is more realistic, aligning with stochastic network modeling frameworks (see, e.g., bickel2009nonparametric and schweinberger2020exponential).
Aside from the weights matrix $\boldsymbol{W}$, the models and numerical studies in the above literature also consider that the regressors $\boldsymbol{x}_i$ in $\boldsymbol{X}=(\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n)^\top$ are observed or realized random vectors. However, for simplicity, some studies of the above literature assume that once observed, $\boldsymbol{W}$ and $\boldsymbol{X}$ are nonstochastic, with asymptotic results relying on this assumption (see, e.g., kelejian1998generalized, kelejian2001asymptotic, lee2004asymptotic, huang2019least, huang2020two, zhu2020multivariate). This implicitly assumes $\boldsymbol{x}_i$s are independent across locations $i=1,\cdots, n$, as nonstochastic constants are mutually independent, and that $\boldsymbol{X}$ and $\boldsymbol{W}$ are also independent. However, in real applications, $\boldsymbol{x}_i$ is often spatially dependent in $i$ alongside $y_i$, making $\boldsymbol{X}$ potentially dependent on $\boldsymbol{W}$. As noted in Chapter 1.4 of wooldridge2002econometric, assuming nonstochastic $\boldsymbol{X}$ may introduce limitations in econometric practice.
Building on the above motivations, we propose a new theoretical framework to study the asymptotic properties of the proposed estimators under random $\boldsymbol{W}$ and $\boldsymbol{X}$, thereby largely broadening SAR model applications while making the classical nonstochastic setting a special case. To this end, we develop a new law of large numbers (LLN) and central limit theorem (CLT) for a general linear-quadratic form $\boldsymbol{s}_n=(\boldsymbol{\varepsilon}^\top \boldsymbol{A}_1\boldsymbol{\varepsilon},\cdots, \boldsymbol{\varepsilon}^\top \boldsymbol{A}_d\boldsymbol{\varepsilon})^\top +\boldsymbol{B}^\top \boldsymbol{\varepsilon}$, where without loss of generality, the random matrices $\boldsymbol{A}_j\in\mathbb{R}^{n\times n}$ and $\boldsymbol{B}\in\mathbb{R}^{n\times d}$ for $j=1,\cdots,d$ are assumed to be measurable functions of $\boldsymbol{W}$ and $\boldsymbol{X}$. Unlike kelejian2001asymptotic, which considers only \( d = 1 \) with nonstochastic \( \boldsymbol{A}_1 \) and \( \boldsymbol{B} \), {\color{black} as well as xu2015maximum and wu2022applications, which allow random \( \boldsymbol{X} \) but assume nonstochastic \( \boldsymbol{W} \)}, our framework accommodates randomness in both \( \boldsymbol{W} \) and \( \boldsymbol{X} \). It turns out that this new framework requires non-trivial conditions on $\boldsymbol{\varepsilon}$ beyond standard exogeneity, and the random $\boldsymbol{X}$ setting can relax the condition of uniform boundedness in absolute value for the elements of $\boldsymbol{X}$ imposed in classical spatial econometrics (see, e.g., kelejian1998generalized, kelejian2001asymptotic and lee2004asymptotic). More details are provided in Section (ref). Notably, both the quasi-score matching estimation and our theoretical framework can be applied in various spatial econometrics models when their $\boldsymbol{W}$ and $\boldsymbol{X}$ are random (see, e.g., lee2010estimation, zou2021network and wu2021inward).
The rest of this paper is organized as follows. Section (ref) introduces the quasi-score matching estimation for the SAR model ((ref)). Section (ref) establishes the new LLN and CLT to study the consistency and asymptotic normality of the quasi-score matching estimators. Sections (ref) and (ref) provide simulation studies and a case study, respectively. Section (ref) gives concluding remarks. To save space, additional simulation and case study results are presented in Sections (ref) and (ref) of the supplementary material, respectively. All theoretical proofs, technical lemmas, and additional theoretical results including corollaries, are also relegated to the supplementary material.
In this section, we denote $\boldsymbol{\theta}_0=(\lambda_0,\boldsymbol{\beta}_0^\top,\sigma_0^2)^\top$ as the true parameter vector in model ((ref)). Subsequently, the true SAR model becomes
where the error vector $\boldsymbol{\varepsilon}$ has mean $\boldsymbol{0}_n$ and covariance matrix $\sigma_0^2 \boldsymbol{I}_n$.
In order to estimate $\boldsymbol{\theta}_0$, the QMLE is commonly considered. Specifically, the estimator maximizes the log-likelihood function based on a conditional Gaussian density of $\boldsymbol{y}$ given $\boldsymbol{X}$ and $\boldsymbol{W}$, which is
where $\boldsymbol{S}(\lambda)=\boldsymbol{I}_n-\lambda \boldsymbol{W}$, $Z(\boldsymbol{\theta}|\boldsymbol{W})=(2\pi\sigma^2)^{n/2} \left| \det \{\boldsymbol{S}(\lambda)\}\right|^{-1}$ is a normalizing constant such that $\int p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W}){\rm d} \boldsymbol{y} =1$, and $\det \{\boldsymbol{S}(\lambda)\}$ is the determinant of the $n\times n$ matrix $\boldsymbol{S}(\lambda)$. The QMLE can be obtained by maximizing the log-likelihood function
Since the data $\boldsymbol{y}$ is not required to be Gaussian in the SAR model ((ref)), $\ell_n(\boldsymbol{\theta})$ is called a quasi-log-likelihood. In real practice, however, the concentrated likelihood approach in lee2004asymptotic is adopted instead of directly optimizing the quasi-log-likelihood function ((ref)). In particular, given $\lambda$, we first maximize $\ell_n(\lambda,\boldsymbol{\beta},\sigma^2)$ with respect to $\boldsymbol{\beta}$ and $\sigma^2$, which yields the optimal $\boldsymbol{\beta}$ and $\sigma^2$ given below:
where $\boldsymbol{M}_{\boldsymbol{X}}=\boldsymbol{I}_n-\boldsymbol{X}(\boldsymbol{X}^\top \boldsymbol{X})^{-1}\boldsymbol{X}$. Then we substitute $\boldsymbol{\beta}$ and $\sigma^2$ in $\ell_n(\lambda,\boldsymbol{\beta},\sigma^2)$ with $\tilde\boldsymbol{\beta}(\lambda)$ and $ \tilde\sigma^2(\lambda)$ to attain the concentrated log-likelihood function
where we omit those constants that are irrelevant to $\lambda$ in $\ell_n\big( \lambda, \tilde\boldsymbol{\beta}(\lambda), \tilde\sigma^2(\lambda)\big)$. Next, we maximize $\ell_n^c(\lambda)$ with respect to $\lambda$, substitute $\lambda$ in ((ref)) and ((ref)) with the optimal $\lambda$, and obtain the QMLE $\tilde\boldsymbol{\theta}=(\tilde\lambda,\tilde\boldsymbol{\beta}^\top,\tilde\sigma^2)^\top$, where
According to the above approach, one can see that the calculation of the QMLE $\tilde\boldsymbol{\theta}=(\tilde\lambda,\tilde\boldsymbol{\beta}^\top,\tilde\sigma^2)^\top$ involves computing $\det\{\boldsymbol{S}(\lambda)\}$, the determinant of the $n\times n$ matrix $\boldsymbol{S}(\lambda)$, which has a computational complexity of $O(n^3)$. As a consequence, the QMLE is not computationally scalable for a large $n$. This motivates us to propose an alternative approach to reduce the computational complexity of estimation. It is worth noting that the determinant $\det\{\boldsymbol{S}(\lambda)\}$ comes from the normalizing constant $Z(\boldsymbol{\theta}|\boldsymbol{W})$ in the conditional density ((ref)). If we can avoid computing the normalizing constant $Z(\boldsymbol{\theta}|\boldsymbol{W})$ in the estimation procedure, the computational complexity of the QMLE will be reduced. To this end, we adapt the score matching approach in the machine learning literature (see, e.g., hyvarinen2005estimation), which is originally proposed to remove the normalizing constant in estimation.
In what follows, we demonstrate three important differences between our proposed score matching approach in this paper and the original score matching in hyvarinen2005estimation. (I) Based on our understanding, this article is the first to propose a score matching approach for spatially dependent data $y_1,\cdots,y_n$ which follow the SAR model ((ref)), and the asymptotic results are established under $n\to\infty$; see Theorems (ref) -- (ref) below. In contrast, the original score matching is proposed for $T$ independent and identical replications of $\boldsymbol{y}=(y_1,\cdots,y_n)^\top$, say $\boldsymbol{y}^{(1)},\cdots,\boldsymbol{y}^{(T)}$, where the asymptotic result of hyvarinen2005estimation is based on fixed $n$ and $T\to\infty$. (II) The original score matching in hyvarinen2005estimation does not consider the regression setting. Therefore, their approach is formulated on the basis of the density of $\boldsymbol{y}$, say $p_{\boldsymbol{\theta}_0}(\boldsymbol{y})$. On the contrary, the SAR model ((ref)) is under a regression setting and hence we build up our score matching approach via the conditional density $p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})$ in ((ref)). (III) The original formulation of score matching assumes that $\boldsymbol{y}^{(t)}$ for $t=1,\cdots,T$ have a true density of $p_{\boldsymbol{\theta}_0}(\boldsymbol{y})$. However, following the idea of the QMLE, we do not require data $\boldsymbol{y}$ to be Gaussian, i.e., $\boldsymbol{y}$ only needs to satisfy the SAR model ((ref)) but does not necessarily follow the Gaussian density $p_{\boldsymbol{\theta}_0}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})$. Hence, we name our proposed approach quasi-score matching estimation. {\color{black} In fact, the proposed quasi-score matching is also developed based on the density $p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})$ in ((ref)), or equivalent, the likelihood, with details provided below.}
Let
and then we have $ p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})=q_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})/{Z(\boldsymbol{\theta}|\boldsymbol{W})}. $ To propose the quasi-score matching estimation, we define a squared distance by
which measures the distance between {\color{black} $\partial \log p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})/{\partial \boldsymbol{y}}=(\partial \log p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})/{\partial y_1},\cdots,$ $\partial \log p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})/{\partial y_n})^\top$} and $\partial \log p_{\boldsymbol{\theta}_0}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})/{\partial \boldsymbol{y}}$, where $\|\cdot\|_2$ is the Euclidean norm. It is worth noting that the machine learning literature refers $\partial \log p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})/{\partial \boldsymbol{y}}$ as a “score” function (see, e.g., hyvarinen2005estimation), which is different from the classical score function $\partial \log p_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})/{\partial \boldsymbol{\theta}}$ considered in econometrics. In spite of this, we still follow the machine literature and use the terminology “score matching" in this paper.
Motivated by the score matching literature, there are two motivations for considering the squared distance $J(\boldsymbol{\theta})$. First, the normalizing constant $Z(\boldsymbol{\theta}|\boldsymbol{W})$ in ((ref)) disappears within the norm $\|\cdot\|_2$; see ((ref)). This will eliminate the need to compute $Z(\boldsymbol{\theta}|\boldsymbol{W})$ in our quasi-score matching estimation introduced later. Second, {\color{black} $J(\boldsymbol{\theta})$ in ((ref)) achieves its minimum, $J(\boldsymbol{\theta}_0) = 0$, at $\boldsymbol{\theta} = \boldsymbol{\theta}_0$. Consequently, minimizing $J(\boldsymbol{\theta})$ yields the true parameter vector $\boldsymbol{\theta}_0$.}
Next, {\color{black} Lemma (ref) in the supplementary material demonstrates that ((ref)) can be transformed as follows:}
where the constant in ((ref)) is irrelevant to $\boldsymbol{\theta}$. {\color{black} This result is obtained via integration by parts; see the proof of Lemma (ref) in the supplementary material.} Since $\boldsymbol{\theta}_0$ minimizes $J(\boldsymbol{\theta})$, if we ignore the constant in ((ref)), we obtain that $\boldsymbol{\theta}_0$ should also minimize \[ D(\boldsymbol{\theta}) =\int \sum_{i=1}^n \left[\frac{\partial^2}{\partial y_i^2}\log q_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})+\frac{1}{2}\left\{ \frac{\partial}{\partial y_i}\log q_{\boldsymbol{\theta}}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})\right\}^2\right] p_{\boldsymbol{\theta}_0}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W}){\rm d} \boldsymbol{y}. \] Let
We subsequently have ${\rm E}_{\boldsymbol{y}\sim p_{\boldsymbol{\theta}_0}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})}\{D_n(\boldsymbol{\theta})\}=D(\boldsymbol{\theta})$, where ${\rm E}_{\boldsymbol{y}\sim p_{\boldsymbol{\theta}_0}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})}\{\cdot\}$ indicates the expectation taken over $\boldsymbol{y}$ with $\boldsymbol{y}$ following the conditional density $p_{\boldsymbol{\theta}_0}(\boldsymbol{y}|\boldsymbol{X}, \boldsymbol{W})$. As a consequence, minimizing the objective function $D_n(\boldsymbol{\theta})$ leads to an $M$-estimator (van1998asymptotic). This estimator is then the quasi-score matching estimator proposed in this paper. Since the normalizing constant $Z(\boldsymbol{\theta}|\boldsymbol{W})=(2\pi\sigma^2)^{n/2} \left|\det\{\boldsymbol{S}(\lambda)\}\right|^{-1}$ also disappears in the quasi-score matching objective function $D_n(\boldsymbol{\theta})$ in ((ref)), obtaining the quasi-score matching estimator does not require the computation of the $n\times n$ matrix $\boldsymbol{S}(\lambda)$'s determinant, and hence reduces the computational complexity.
To further demonstrate the computational complexity of the quasi-score matching estimation, we derive the closed-form of $D_n(\boldsymbol{\theta})$ based on ((ref)) and ((ref)), which is given below:
Similar to the QMLE, we adopt the concentrated approach when minimizing the objective function $D_n(\boldsymbol{\theta})$, and obtain the quasi-score matching estimator. In particular, {\color{black} we obtain the first-order condition (FOC) of minimizing $D_n(\boldsymbol{\theta})$, which is ${\partial D_n(\boldsymbol{\theta})}/{\partial \boldsymbol{\theta}}=\big( {\partial D_n(\boldsymbol{\theta})}/{\partial \lambda},{\partial D_n(\boldsymbol{\theta})}/{\partial \boldsymbol{\beta}^\top},{\partial D_n(\boldsymbol{\theta})}/{\partial \sigma^2}\big)^\top=\boldsymbol{0}_{p+2}$, where the closed form of ${\partial D_n(\boldsymbol{\theta})}/{\partial \boldsymbol{\theta}}$ is given in Lemma (ref) of the supplementary material.} Solving $\big( {\partial D_n(\boldsymbol{\theta})}/{\partial \boldsymbol{\beta}^\top},{\partial D_n(\boldsymbol{\theta})}/{\partial \sigma^2}\big)^\top=\boldsymbol{0}_{p+1}$ for a given $\lambda$, yields the optimal $\boldsymbol{\beta}$ and $\sigma^2$:
where $\boldsymbol{Q}_{\boldsymbol{X}}(\lambda)=\boldsymbol{I}_n-\boldsymbol{S}(\lambda)^\top \boldsymbol{X}\left\{ \boldsymbol{X}^\top \boldsymbol{S}(\lambda)\boldsymbol{S}(\lambda)^\top \boldsymbol{X}\right\}^{-1}\boldsymbol{X}^\top \boldsymbol{S}(\lambda)$. We subsequently substitute $\boldsymbol{\beta}$ and $\sigma^2$ in $(\ref{eq:obj})$ with $\hat\boldsymbol{\beta}(\lambda)$ and $ \hat\sigma^2(\lambda)$, to attain the concentrated $D_n(\lambda,\boldsymbol{\beta},\sigma^2)$, which is
Then the quasi-score matching estimators of $\lambda$, $\boldsymbol{\beta}$ and $\sigma^2$ are obtained via
{\color{black} In Remark (ref) of Section (ref) in the supplementary material, we demonstrate the numerical visualization of the function $D^c_n(\lambda)$ based on the simulation settings considered in Section (ref), showing that minimizing $D_n^c(\lambda)$ yields a unique solution numerically.}
{\color{black} Based on the above concentrated procedure, for a given $\lambda$, the quasi-score matching estimator $\hat\boldsymbol{\beta}(\lambda)$ in ((ref)) is a weighted least squares estimator, whereas the QMLE $\tilde\boldsymbol{\beta}(\lambda)$ in ((ref)) is an ordinary least squares estimator. Introducing the unnecessary weight $\boldsymbol{S}(\lambda)\boldsymbol{S}(\lambda)^\top$ in $\hat\boldsymbol{\beta}(\lambda)$ naturally reduces the asymptotic efficiency of the quasi-score matching estimator; also see Section (ref) and Corollary (ref) of the supplementary material, which theoretically show that the QMLE $\tilde\boldsymbol{\theta}$ is asymptotically more efficient than the quasi-score matching estimator $\hat\boldsymbol{\theta}=(\hat\lambda,\hat\boldsymbol{\beta}^\top,\hat\sigma^2)^\top$. To improve the asymptotic efficiency, we can substitute $\lambda$ in QMLE ((ref)) and ((ref)) with the quasi-score matching estimator $\hat\lambda$, yielding an improved estimator introduced in Section (ref). }
Apart from trading off some estimation efficiency compared to the QMLE, one can see that the computational complexity of obtaining the quasi-score matching estimator ((ref)) is of order $O(n^2)$ for a large $n$; this is because the computation of $\det\{\boldsymbol{S}(\lambda)\}$ is not needed, and the computational complexity is mainly determined by $n$-dimensional matrices' multiplications. If we consider that the weights matrix $\boldsymbol{W}$ has at most $c$ non-zeros entries per row and per column (see, e.g., some $\boldsymbol{W}$ considered in kelejian1998generalized and yu2016strategic given each location having at most $c$ neighbors, and some normalized network adjacency matrices $\boldsymbol{W}$ considered in huang2019least, huang2020two and zhu2020multivariate), the computational complexity can be further reduced to $O(c n)$, which is linear to $n$ for a finite $c$. Accordingly, the quasi-score matching estimation is more computationally scalable than the QMLE as $n$ is large. It is worth noting that our proposed quasi-score matching estimator is also superior to the GMM estimator by lee2007gmm and the one-step estimator by gupta2023efficient in terms of computation. Notably, both methods require computing $\boldsymbol{S}^{-1}(\lambda)$, which has a computational complexity of $O(n^3)$. Numerical comparisons with these and other existing methods are provided in Section (ref) of the supplementary material.
In addition to analyzing the computational complexity, we next study the asymptotic properties of the quasi-score matching estimator $\hat\boldsymbol{\theta}=(\hat\lambda,\hat\boldsymbol{\beta}^\top,\hat\sigma^2)^\top$ for the SAR model ((ref)) with both random weights matrix $\boldsymbol{W}$ and regressors $\boldsymbol{X}$.
To develop the asymptotic properties of the quasi-score matching estimator $\hat\boldsymbol{\theta}$, the main steps include establishing the CLT for $n^{-1/2}\partial D_n(\boldsymbol{\theta}_0)/\partial \boldsymbol{\theta}$ and the LLN for $n^{-1}\partial^2 D_n(\boldsymbol{\theta}_0)/$ $(\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^\top)$. Denote $[\boldsymbol{G}]_s=(\boldsymbol{G}+\boldsymbol{G}^\top)/2$ as the symmetrization of any generic matrix $\boldsymbol{G}$, and $\boldsymbol{0}_{n_1\times n_2}$ as the $n_1\times n_2$ matrix of zeros. By Lemma (ref) of the supplementary material, we obtain that
where $ \boldsymbol{A}_1(\boldsymbol{\theta}_0)=-\left\{ \left[ \boldsymbol{S}(\lambda_0) \boldsymbol{S}(\lambda_0)^\top \boldsymbol{W} \boldsymbol{S}^{-1}(\lambda_0)\right]_s + \left[ \boldsymbol{S}(\lambda_0)\boldsymbol{W}^\top \right]_s\right\}/\sigma_0^4$, $ \boldsymbol{A}_2(\boldsymbol{\theta}_0)=\cdots=\boldsymbol{A}_{p+1}$ $(\boldsymbol{\theta}_0)=\boldsymbol{0}_{n\times n}$, $\boldsymbol{A}_{p+2}(\boldsymbol{\theta}_0)=-\boldsymbol{S}(\lambda_0)\boldsymbol{S}(\lambda_0)^\top/\sigma_0^6$, and
Note that the first two terms on the right hand side of ((ref)) constitute a linear-quadratic form of $\boldsymbol{\varepsilon}$ with matrices $\boldsymbol{A}_j(\boldsymbol{\theta}_0)$ and $\boldsymbol{B}(\boldsymbol{\theta}_0)$ for $j=1,\cdots,(p+2)$, where $\boldsymbol{A}_j(\boldsymbol{\theta}_0)$ and $\boldsymbol{B}(\boldsymbol{\theta}_0)$ are measurable functions of $\boldsymbol{W}$ and $\boldsymbol{X}$. Meanwhile, $\partial^2 D_n(\boldsymbol{\theta}_0)/(\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^\top)$ also contains a linear-quadratic form of $\boldsymbol{\varepsilon}$ with matrices that are functions of $\boldsymbol{W}$ and $\boldsymbol{X}$; see the closed form of $\partial^2 D_n(\boldsymbol{\theta}_0)/(\partial \boldsymbol{\theta}\partial\boldsymbol{\theta}^\top)$ given in ((ref)) of the supplementary material. Therefore, the key to obtain the asymptotic properties of the quasi-score matching estimator is to establish the LLN and CLT for a general linear-quadratic form
where matrices $\boldsymbol{A}_j$ and $\boldsymbol{B}$ are measurable functions of $\boldsymbol{W}$ and $\boldsymbol{X}$ for $j=1,\cdots, d$ and $1\leq d<\infty$. Without loss of generality, we only need to consider $\boldsymbol{A}_j$ being symmetric in the linear-quadratic form of ((ref)); this is because that if $\boldsymbol{A}_j$ is not symmetric, the quadratic form $\boldsymbol{\varepsilon}^\top\boldsymbol{A}_j\boldsymbol{\varepsilon}\equiv \boldsymbol{\varepsilon}^\top[\boldsymbol{A}_j]_s\boldsymbol{\varepsilon}$, where the symmetrization $[\boldsymbol{A}_j]_s$ is always a symmetric matrix.
As discussed in Introduction, in real data analysis it is often the case that $\boldsymbol{W}$ and $\boldsymbol{X}$ are stochastic, which leads to matrices $\boldsymbol{A}_j$ and $\boldsymbol{B}$ in ((ref)) being random. In Section (ref), we develop a new LLN and a new CLT, which are suitable for the linear-quadratic form $\boldsymbol{s}_n$ with such random matrices. It is worthing noting that the CLT developed in kelejian2001asymptotic only considers $d=1$ and matrices $\boldsymbol{A}_1$ and $\boldsymbol{B}$ being nonstochastic. Hence, our theoretical results are more general, and can be applied to the models with random $\boldsymbol{W}$ and $\boldsymbol{X}$. In addition, it turns out that a non-trivial condition is needed for the error vector $\boldsymbol{\varepsilon}$ compared to the nonstochastic case, which is also discussed in Section (ref) below and Section (ref) of the supplementary material. Next, we apply the new LLN and CLT, and obtain the asymptotic properties of the quasi-score matching estimator in Section (ref). Lastly, we make use of the QMLE ((ref)) and ((ref)), and the quasi-score matching estimator of $\lambda$, to obtain new estimators of $\boldsymbol{\beta}$ and $\sigma^2$ in Section (ref). We further demonstrate that doing this will remain the same computational complexity, but can improve the asymptotic efficiency of $\boldsymbol{\beta}$ and $\sigma^2$ compared to the quasi-score matching estimation.
To establish the LLN and CLT of the linear-quadratic form $\boldsymbol{s}_n$ in ((ref)), note that $\boldsymbol{A}_j$ and $\boldsymbol{B}$ are measurable functions of random $\boldsymbol{W}$ and $\boldsymbol{X}$ for $j=1,\cdots,d$ mentioned below ((ref)). Hence, $\boldsymbol{A}_j$ and $\boldsymbol{B}$ are random and $\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}$-measurable, where $\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}=\sigma\langle \boldsymbol{X},\boldsymbol{W}\rangle$ is the $\sigma$-algebra generated by random matrices $\boldsymbol{X}$ and $\boldsymbol{W}$. In econometrics, it is common to impose the exogenous condition for random regressors $\boldsymbol{X}$ with respect to random errors $\varepsilon_i$. Adapting it to spatial econometrics with both random $\boldsymbol{X}$ and $\boldsymbol{W}$, we consider the exogenous condition given below:
(C1) ${\rm E}[\varepsilon_i|\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}]={\rm E}[\varepsilon_i|\boldsymbol{X},\boldsymbol{W}]=0$ for $i=1,\cdots,n$.
Next, we introduce a non-trivial condition for random errors $\varepsilon_i$.
(C2) Assume that $\varepsilon_i$ are conditional independent in $i$ given $\sigma$-algebra $\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}$ for $i=1,\cdots,n$.
It is worth noting that the classical spatial econometrics (see, e.g., lee2004asymptotic and lesage2009introduction) often considers the more commonly assumed condition given below for random errors $\varepsilon_i$ under the nonstochastic setting of $\boldsymbol{W}$ and $\boldsymbol{X}$:
(C2$'$) Assume that $\varepsilon_i$ are independent in $i$ for $i=1,\cdots,n$.
However, independence condition (C2$'$) does not imply conditional independence condition (C2), and the opposite implication is also not valid (see stoyanov2013counterexamples for counterexamples). Further discussion on the necessity of the non-trivial Condition (C2) for establishing the CLT of the linear-quadratic form with random matrices is presented in Section (ref) of the supplementary material.
After discussing the conditions for random errors $\varepsilon_i$ in Conditions (C1) -- (C2) above, we impose conditions for random matrices $\boldsymbol{A}_j$ and $\boldsymbol{B}$ in ((ref)) as follows.
(C3) Suppose that $\boldsymbol{A}_j=(a_{i_1i_2,j})_{n\times n}$ are $n\times n$ random symmetric matrices and $\boldsymbol{B}=(\boldsymbol{b}_1,\cdots,\boldsymbol{b}_d)=(b_{ij})_{n\times d}$ is an $n\times d$ random matrix for $i,i_1,i_2\in\{1,\cdots,n\}$ and $j=1,\cdots, d$, where $\boldsymbol{b}_j$ is the $j$-th column of $\boldsymbol{B}$. Assume that $\boldsymbol{A}_j$, $j=1,\cdots, d$, and $\boldsymbol{B}$ are all $\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}$-measurable. In addition, assume that (i) there exists $\eta_1>0$ such that ${\rm E}|\varepsilon_i|^{4+\eta_1}<\infty$; (ii) there exists $\eta_2\geq 8/\eta_1$ such that ${\rm E}|a_{i_1i_2,j}|^{2+\eta_2}<\infty$; and (iii) there exists $\eta_3\geq 4/(2+\eta_1)$ such that ${\rm E}|b_{ij}|^{2+\eta_3}<\infty$.
We next introduce a stronger condition for random matrices $\boldsymbol{A}_j$ and $\boldsymbol{B}$ below, which is used to obtain the new LLN for the linear-quadratic form ((ref)). Before introducing this condition, let $\|\cdot\|_\eta$ be the vector $\eta$-norm or the matrix $\eta$-norm for $1\leq \eta\leq \infty$. Specifically, for any generic vector $\boldsymbol{\alpha}=(\alpha_1,\cdots,\alpha_q)\in\mathbb{R}^q$, $\|\boldsymbol{\alpha}\|_\eta=(\sum_{k=1}^q|\alpha_k|^\eta)^{1/\eta}$, and, for any generic matrix $\boldsymbol{G}\in\mathbb{R}^{m\times q}$,
If matrix $\boldsymbol{G}=(g_{k_1k_2})_{m\times q}$ is random, then, for $1\leq \eta\leq\infty$, denote $\|g_{k_1k_2}\|_{L^\eta}=({\rm E}|g_{k_1k_2}|^\eta)^{1/\eta}$ as the $L^\eta$-norm of the random variable $g_{k_1k_2}$. Subsequently, let
and it is an $m\times q$ matrix consisting of the $L^\eta$-norm of $g_{k_1k_2}$. Moreover, denote
for $\eta_1,\eta_2\in[1,\infty]$, where $\eta_1$ and $\eta_2$ correspond to the $L^{\eta_1}$-norm ((ref)) and the matrix $\eta_2$-norm ((ref)), respectively. Lastly, let ${\rm vec}(\boldsymbol{G})$ be the vectorization for any generic matrix $\boldsymbol{G}$.
(C4) For $\eta_1$, $\eta_2$ and $\eta_3$ in (C3), assume that (i) $\sup_i{\rm E}|\varepsilon_i|^{4+\eta_1}<\infty$; (ii) $\sup_{n\geq 1}{\vert\kern-0.25ex\vert\kern-0.25ex\vert \boldsymbol{A}_j \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{L^{2+\eta_2},1}$ $<\infty$ for $j=1,\cdots,d$; and (iii) $\sup_{n\geq 1}n^{-1}{\rm E}\|{\rm vec}(\boldsymbol{B})\|_{2+\eta_3}^{2+\eta_3}<\infty$.
Conditions (C3) and (C4) only impose some moment conditions on random matrices $\boldsymbol{A}_j$ and $\boldsymbol{B}$. Compared to the literature, Assumption 2 in kelejian2001asymptotic is the nonstochastic version of these conditions for the linear-quadratic form $\boldsymbol{s}_n$ in ((ref)) with fixed $\boldsymbol{A}_j$ and $\boldsymbol{B}$, and $d=1$. However, since we consider matrices $\boldsymbol{A}_j$ and $\boldsymbol{B}$ being random, their $(2+\eta_2)$-th and $(2+\eta_3)$-th order moments, respectively, need to be coordinated with $\varepsilon_i$'s $(4+\eta_1)$-th order moments in a way of those demonstrated in Conditions (C3) -- (C4).
Let $\boldsymbol{\Upsilon}^{(2)}={\rm diag}\{\mu_1^{(2)},\cdots,\mu_n^{(2)}\}$, $\boldsymbol{\Upsilon}^{(3)}={\rm diag}\{\mu_1^{(3)},\cdots,$ $\mu_n^{(3)}\}$ and $\boldsymbol{\Upsilon}^{(4)}={\rm diag}\{\mu_1^{(4)}-3(\mu_1^{(2)})^2,\cdots,\mu_n^{(4)}-3(\mu_n^{(2)})^2\}$ with $\mu_i^{(s)}={\rm E}[\varepsilon_i^s|\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}]$ for $s=2,3,4$, where ${\rm diag}\{d_{1},\cdots,d_{n}\}$ denotes an $n\times n$ diagonal matrix with diagonals being $d_{1},\cdots,d_{n}$. In addition, let $(g_{k_1k_2})_{m\times q}$ denote a generic $m\times q$ matrix with the $(k_1,k_2)$-th element being $g_{k_1k_2}$ for $k_1=1,\cdots, m$ and $k_2=1,\cdots,q$, and $\circ$ represent the Hadamard product of any two matrices with the same dimensions. The above notation and Conditions (C1) -- (C4) yield the following theorem for the linear-quadratic form $\boldsymbol{s}_n$ of ((ref)).
Theorem (ref) (i) provides the closed forms of conditional mean and conditional covariance matrix of the linear-quadratic form $\boldsymbol{s}_n$. This is further utilized to show the LLN in Theorem (ref) (ii), and to obtain the closed form of ${\rm Cov}[\boldsymbol{s}_n|\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}]$, the latter of which is required in the CLT for the linear-quadratic form $\boldsymbol{s}_n$ in Theorem (ref) given below. It is worth noting that $\eta$ in Theorem (ref) (ii) is a general positive constant. When we apply Theorem (ref) (ii) in developing the theorems in the following subsections, we mainly use the result $n^{-1}\left( \boldsymbol{s}_n-{\rm E}\left[ \boldsymbol{s}_n|\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}\right]\right)\stackrel{{\rm P}}\longrightarrow \boldsymbol{0}_d$ with $\eta=1/2$.
In addition to Conditions (C1) -- (C4) introduced earlier, the CLT of Theorem (ref) requires an extra condition:
(C5) Assume that (C4) holds for $\eta_2>8/\eta_1$ and $\eta_3>4/(2+\eta_1)$.
The above condition, together with Conditions (C1) -- (C4), leads to the CLT for $\boldsymbol{s}_n$.
As mentioned in the beginning of this section, the CLT of Theorem (ref) is developed to show the asymptotic normality of $n^{-1/2}\partial D_n(\boldsymbol{\theta}_0)/\partial \boldsymbol{\theta}$, and the LLN of Theorem (ref) (ii) is established to prove that $n^{-1}\partial^2 D_n(\boldsymbol{\theta}_0)/(\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^\top)$ converges in probability. Both of these results are given in Lemma (ref) of the supplementary material, and they further lead to the asymptotic normality of the quasi-score matching estimator $\hat\boldsymbol{\theta}$ in Theorem (ref) of the next subsection. It is worth noting that for the generality reason, both Theorems (ref) and (ref) are developed under the heterogeneity of random errors $\varepsilon_i$. The classical SAR model ((ref)), on the other hand, assumes the homoskedasticity of random errors $\varepsilon_i$, which is a special case. In addition, Theorems (ref) and (ref) can be applied in many other spatial econometrics models with random $\boldsymbol{W}$ and $\boldsymbol{X}$, as long as the linear-quadratic form is involved in deriving the asymptotic properties of estimators (see, e.g., lee2010estimation, zou2021network and wu2021inward).
In this subsection, we apply Theorems (ref) and (ref) to establish the asymptotic normality of the proposed quasi-score matching estimator $\hat\boldsymbol{\theta}$ for the SAR model ((ref)) with both random weights matrix $\boldsymbol{W}$ and regressors $\boldsymbol{X}$ in Theorem (ref) below. We first lay out the technical assumptions required in this theorem.
The assumptions for the true parameter value $\lambda_0$ and $\boldsymbol{S}(\lambda)$ are the same as those in the classical spatial econometrics literature (see, e.g., lee2004asymptotic and lesage2009introduction), {\color{black} but the assumption for $\boldsymbol{S}(\lambda)$ needs to hold almost surely for random $\boldsymbol{S}(\lambda)$.} In addition, we impose the same conditions of Theorems (ref) and (ref) for the random errors $\varepsilon_i$ such that we can apply those theorems in developing the asymptotic normality of $\hat\boldsymbol{\theta}$ later. We further assume that the random errors have the same conditional variance. This assumption can lead to the homoskedasticity of the classical SAR model ((ref)) by the tower property of conditional expectation.
Before we introduce the next assumption, we introduce some notation below. For a generic matrix $\boldsymbol{G}=( g_{k_1k_2})_{m\times q}$, if $g_{k_1k_2}$ is a function of $\lambda\in \Lambda$, namely $\boldsymbol{G}=\boldsymbol{G}(\lambda)=\big( g_{k_1k_2}(\lambda)\big)_{m\times q}$, then we denote $\sup_{\lambda\in \Lambda}|\boldsymbol{G}(\lambda)|=\big( \sup_{\lambda\in \Lambda}|g_{k_1k_2}(\lambda)|\big)_{m\times q}$. Based on this notation, we introduce the assumption associated to the random weights matrix $\boldsymbol{W}$.
If $\boldsymbol{W}$ is nonstochastic, Assumption (ref) reduces to the uniform boundedness of $\|\boldsymbol{W}\|_1$, $\|\boldsymbol{W}\|_\infty$, $\|\boldsymbol{S}^{-1}(\lambda)\|_1$ and $\|\boldsymbol{S}^{-1}(\lambda)\|_\infty$, which is often considered in spatial econometrics (e.g., see Assumption 5 of lee2004asymptotic). {\color{black} More interpretations of this assumption, including the relationships between ${\vert\kern-0.25ex\vert\kern-0.25ex\vert \boldsymbol{W} \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{L^{\delta_0(2+\eta_2)},1}$ (or ${\vert\kern-0.25ex\vert\kern-0.25ex\vert \boldsymbol{W} \vert\kern-0.25ex\vert\kern-0.25ex\vert}_{L^{\delta_0(2+\eta_2)},\infty}$) and $\|\boldsymbol{W}\|_1$ (or $\|\boldsymbol{W}\|_\infty$), are provided in Section (ref) of the supplementary material.} {\color{black} Based on these interpretations, the examination of Assumption (ref) is discussed in Section (ref) of the supplementary material.} {\color{black} Since $\boldsymbol{W}$ is random, some moment conditions for random matrices $\boldsymbol{W}$ and $\boldsymbol{S}^{-1}(\lambda)$ via the definition of ${\vert\kern-0.25ex\vert\kern-0.25ex\vert \cdot \vert\kern-0.25ex\vert\kern-0.25ex\vert}$ in ((ref)), are required in order to establish the asymptotic normality of $\hat\boldsymbol{\theta}$. These moment conditions also need to be coordinated with $\varepsilon_i$'s $(4+\eta_1)$-th order moments in a way of those demonstrated in Assumption (ref).} Similarly, we next introduce the moment condition for the random regressors $\boldsymbol{X}$ below.
Assumption (ref) does not impose any independence conditions on $\boldsymbol{x}_i$ for $i=1,\cdots, n$, or $x_{ij}$ for $j=1,\cdots, p$. {\color{black} Instead, we only assume some moment condition on $x_{ij}$.} This allows the dependence of $x_{ij}$ in both $i$ and $j$, and is more general than the nonstochastic $\boldsymbol{X}$ setting, since any nonstochastic constants $x_{ij}$ are mutually independent for $i=1,\cdots,n $ and $j=1,\cdots,p$. In addition, this assumption allows the dependence between $\boldsymbol{X}$ and $\boldsymbol{W}$. {\color{black} Moreover, the classical spatial econometrics literature (see, e.g., kelejian2001asymptotic and lee2004asymptotic) requires the condition of uniform boundedness in absolute value for the elements of $\boldsymbol{X}$ under the nonstochastic $\boldsymbol{X}$ setting, i.e., $\max_{1\leq i\leq n;1\leq j\leq p}|x_{ij}|\leq C_x$ for some finite positive constant $C_x$. However, if $x_{ij}$ are observations or realizations of random variables considered in some of the models and numerical studies of the above literature, this condition may not be valid. For example, if $x_{ij}\stackrel{iid}\sim N(0,\sigma_x^2)$, we have ${\rm E}(\max_{1\leq i\leq n;1\leq j\leq p}|x_{ij}|)/\sqrt{2\sigma_x^2\log(np)}\to 1$ as $n\to\infty$ (wainwright2019high). Therefore, there exists a realization of the random variables $x_{ij}$ such that $\max_{1\leq i\leq n;1\leq j\leq p}|x_{ij}|/\sqrt{2\sigma_x^2\log(np)}\to 1$ (see, e.g., the similar technique in vershynin2018high), which indicates that this realization of $\max_{1\leq i\leq n;1\leq j\leq p}|x_{ij}|$ is not uniformly bounded as $n\to\infty$. On the other hand, under such a random setting of $x_{ij}$, the moment condition in Assumption (ref) still holds, and hence Assumption (ref) releases the strong uniform boundedness condition for the elements of $\boldsymbol{X}$ in the literature.}
Under Assumptions (ref) -- (ref) and using Theorem (ref) (i), we subsequently obtain that \[ \frac{1}{n}{\rm E}\left[ -\frac{\partial D_n(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}}\biggl{|}\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}\right]=\boldsymbol{0}_{p+2}, \frac{1}{n}{\rm Cov}\left[ -\frac{\partial D_n(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}}\biggl{|}\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}\right]=\boldsymbol{V}_{{\rm S},n}+\boldsymbol{\Omega}_{{\rm S},n}, \textrm{ and } \] \[ \frac{1}{n}{\rm E}\left[ \frac{\partial^2 {\color{black} D_n}(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}\partial\boldsymbol{\theta}^\top}\biggl{|}\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}\right]=\boldsymbol{U}_{{\rm S},n}, \] where the closed forms of $(p+2)\times (p+2)$ matrices $\boldsymbol{V}_{{\rm S},n}$, $\boldsymbol{\Omega}_{{\rm S},n}$ and $\boldsymbol{U}_{{\rm S},n}$ are given by ((ref)), ((ref)) and ((ref)), respectively, in Section (ref) of the supplementary material. Based on $\boldsymbol{V}_{{\rm S}, n}$, $\boldsymbol{\Omega}_{{\rm S},n}$ and $\boldsymbol{U}_{{\rm S},n}$, we introduce the next assumption.
Assumption (ref) is similar to Assumptions 6 and 8 of lee2004asymptotic but under the random $\boldsymbol{W}$ and $\boldsymbol{X}$ setting, where in this case $\boldsymbol{V}_{{\rm S}, n}$, $\boldsymbol{\Omega}_{{\rm S},n}$ and $\boldsymbol{U}_{{\rm S},n}$ obtained in Section (ref) are all random but not nonstochastic. In addition, this assumption, together with Assumptions (ref) -- (ref), leads to the asymptotic covariance matrix of $-n^{-1/2}\partial D_n(\boldsymbol{\theta}_0)/\partial \boldsymbol{\theta}$, and the limit of $n^{-1}\partial^2 D_n(\boldsymbol{\theta}_0)/(\partial \boldsymbol{\theta}\partial \boldsymbol{\theta}^\top)$ in probability, both of which are given in Lemma (ref) of the supplementary material. {\color{black} Finally, we introduce the assumption for identification of $\boldsymbol{\theta}_0$.
Assumption (ref), together with Assumptions (ref) -- (ref), lead to the identification of $\boldsymbol{\theta}_0$ and the consistency of the quasi-score matching estimator $\hat\boldsymbol{\theta}$ given in Lemma (ref) of the supplementary material. Based on Lemma (ref), we demonstrate that when the QMLE identifies $\boldsymbol{\theta}_0$ or is consistent, the quasi-score matching does so as well; see Remark (ref) in Section (ref) of the supplementary material. The intuitive reasons for the consistency of $\hat\boldsymbol{\theta}$ are provided in Remarks (ref) and (ref) in Section (ref) of the supplementary material. Under heteroskedasticity, i.e., when ${\rm E}[\varepsilon_i^2|\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}]=\sigma_0^2$ in Assumption (ref) is violated, $\hat\boldsymbol{\theta}$ is inconsistent, whereas the 2SLS estimator kelejian1998generalized remains consistent and robust; however, the 2SLS estimation does not work in the case that ${\rm E}[\boldsymbol{W} \boldsymbol{y}|\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}]$ and $\boldsymbol{X}$ are linearly dependent.} {\color{black} On the other hand, we show that the quasi-score matching estimation works in this case; see Section (ref) of the supplementary material.}
Lemma (ref), in conjunction with Assumption (ref), yields the asymptotic normality of $\hat\boldsymbol{\theta}$ given in the theorem below.
Theorem (ref) can lead to the asymptotic variance of the quasi-score matching estimator $\hat\boldsymbol{\theta}$, and it is $n^{-1}(\boldsymbol{U}_{{\rm S}}^{-1}\boldsymbol{V}_{{\rm S}}\boldsymbol{U}_{{\rm S}}^{-1}+\boldsymbol{U}_{{\rm S}}^{-1}\boldsymbol{\Omega}_{{\rm S}}\boldsymbol{U}_{{\rm S}}^{-1})$. {\color{black} Although the quasi-score matching approach can alleviate the computational complexity of the QMLE when $n$ is large, we utilize the asymptotic variance obtained from Theorem (ref) to show that it trades off some estimation efficiency compared to the QMLE; see Corollary (ref) in Section (ref) of the supplementary material. Nevertheless, based on the numerical estimation efficiencies demonstrated in various simulations in Section (ref), we find that the quasi-score matching achieves efficiency close to the QMLE, especially for large $n$.} In the following subsection, we propose an approach to improve the asymptotic efficiency of the quasi-score matching estimators $\hat\boldsymbol{\beta}$ and $\hat\sigma^2$ while retaining the computational complexity as $n$ is large.
Recall that the concentrated likelihood approach for the QMLE leads to two estimating equations ((ref)) and ((ref)) that are used to estimate $\boldsymbol{\beta}$ and $\sigma^2$ given $\lambda$. As a consequence, if we substitute $\lambda$ in ((ref)) and ((ref)) with the quasi-score matching estimator $\hat\lambda$, we can have two new estimators of $\boldsymbol{\beta}$ and $\sigma^2$:
One can see that the computational complexity of obtaining these estimators is not increased compared to the calculations of the quasi-score matching estimators $\hat\boldsymbol{\beta}$ and $\hat\sigma^2$ in ((ref)). However, since the information of the QMLE's estimating equations is involved in these new estimators, one can expect an improvement of the efficiency for the estimation of $\boldsymbol{\beta}$ and $\sigma^2$.
In what follows, we study the asymptotic properties of the efficiency improved estimators $\hat{\tilde\boldsymbol{\beta}}=\tilde\boldsymbol{\beta}(\hat\lambda)$ and $\hat{\tilde\sigma}^2=\tilde\sigma^2(\hat\lambda)$ under the random $\boldsymbol{W}$ and $\boldsymbol{X}$ setting. In particular, we show a more general result for the estimator $\hat{\tilde\boldsymbol{\theta}}=(\hat\lambda, \hat{\tilde\boldsymbol{\beta}}^\top,\hat{\tilde\sigma}^2)^\top$, where $( \hat{\tilde\boldsymbol{\beta}}^\top,\hat{\tilde\sigma}^2)^\top$ is just a sub-vector of $\hat{\tilde\boldsymbol{\theta}}$. Since the efficiency improved estimator $\hat{\tilde\boldsymbol{\theta}}=(\hat\lambda, \hat{\tilde\boldsymbol{\beta}}^\top,\hat{\tilde\sigma}^2)^\top$ includes the information of both the quasi-score matching estimator $\hat\lambda$ and the QMLE's estimating equations ((ref)) and ((ref)), the asymptotic distribution of $\hat{\tilde\boldsymbol{\theta}}$ involves the moments of both ${\partial D_n(\boldsymbol{\theta}_0)}/{\partial \boldsymbol{\theta}}$ and ${\partial \ell_n(\boldsymbol{\theta}_0)}/{\partial \boldsymbol{\theta}}$, where $D_n(\boldsymbol{\theta})$ is in ((ref)) from the quasi-score matching, and $\ell_n(\boldsymbol{\theta})$ is in ((ref)) from the QMLE. Specifically, we use Theorem (ref) (i) and Assumptions (ref) -- (ref) to obtain that \[ \frac{1}{n}{\rm Cov}\left[
\Biggl{|}\mathcal{F}_{\boldsymbol{X},\boldsymbol{W}}\right]=\boldsymbol{V}_{n}+\boldsymbol{\Omega}_{n}, \] where the closed forms of $2(p+2)\times 2(p+2)$ matrices $\boldsymbol{V}_{n}$ and $\boldsymbol{\Omega}_{n}$ are given by ((ref)) in Section (ref) of the supplementary material. Similar to Theorem (ref), to show the asymptotic normality of $\hat{\tilde\boldsymbol{\theta}}$, we require an extra assumption on $\boldsymbol{V}_n$ and $\boldsymbol{\Omega}_n$, which is given below.
\begin{as:4'} As $n\to\infty$, assume that \[\boldsymbol{V}_{n}\stackrel{{\rm P}}\longrightarrow \boldsymbol{V}=
, and \boldsymbol{\Omega}_{n}\stackrel{{\rm P}}\longrightarrow \boldsymbol{\Omega}=
,\]where $\boldsymbol{V}$ and $\boldsymbol{\Omega}$ are finite $2(p+2)\times 2(p+2)$ matrices such that $\boldsymbol{V}+\boldsymbol{\Omega}$ and $\boldsymbol{V}$ are positive definite, $\boldsymbol{V}_{{\rm S}}$ and $\boldsymbol{\Omega}_{{\rm S}}$ are $(p+2)\times (p+2)$ matrices defined in Assumption (ref), and $\boldsymbol{V}_{{\rm M}}$, $\boldsymbol{\Omega}_{{\rm M}}$, $\boldsymbol{V}_{{\rm S}{\rm M}}$ and $\boldsymbol{\Omega}_{{\rm S}{\rm M}}$ are finite $(p+2)\times (p+2)$ matrices. \end{as:4'}
Based on this assumption and $\boldsymbol{U}_{{\rm S}}$ in Assumption (ref), let
for $V_{{\rm M},\lambda,\lambda}\in\mathbb{R}$, $\boldsymbol{V}_{{\rm M},\lambda,-\lambda}\in\mathbb{R}^{1\times (p+1)}$, $\boldsymbol{V}_{{\rm M},-\lambda,\lambda}\in\mathbb{R}^{(p+1)\times 1}$ and $\boldsymbol{V}_{{\rm M},-\lambda,-\lambda}\in\mathbb{R}^{(p+1)\times (p+1)}$, and denote
The notation, together with Assumptions (ref) -- (ref) and (ref)$'$, leads to the asymptotic normality of $\hat{\tilde\boldsymbol{\theta}}$ given below.
{\color{black} Based on this theorem, Theorem (ref), and Lemma (ref) in Section (ref) of the supplementary material, we show that the QMLE $\tilde\boldsymbol{\theta}$ is asymptotically more efficient than either the quasi-score matching estimator $\hat\boldsymbol{\theta}=(\hat\lambda,\hat\boldsymbol{\beta}^\top,\hat\sigma^2)^\top$ or the efficiency improved estimator $\hat{\tilde\boldsymbol{\theta}}=(\hat\lambda, \hat{\tilde\boldsymbol{\beta}}^\top,\hat{\tilde\sigma}^2)^\top$; see Section (ref) and Corollary (ref) of the supplementary material. In addition, Section (ref) provides numerical justification that the asymptotic efficiency of $\hat{\tilde\boldsymbol{\theta}}$ is better than that of $\hat\boldsymbol{\theta}$, and our simulation studies in Section (ref) further confirm this result.} {\color{black} As suggested by an anonymous reviewer, we also prove that the efficiency improved estimator is asymptotically more efficient than the best 2SLS estimator lee2003best under some conditions; see Section (ref) and Corollary (ref) of the supplementary material.}
Note that both the calculations of $\hat{\tilde\boldsymbol{\theta}}$ and $\hat\boldsymbol{\theta}$ are simpler and faster than that of $\tilde\boldsymbol{\theta}$ as $n$ is large, while the asymptotic variances of the three estimators are all of order $O(n^{-1})$. We next provide the simulation studies to compare the performance of $\hat\boldsymbol{\theta}$, $\hat{\tilde\boldsymbol{\theta}}$ and $\tilde\boldsymbol{\theta}$ in Section (ref) below.
In this section, we conduct various simulation studies to evaluate the finite sample performance of the quasi-score matching estimator $\hat\boldsymbol{\theta}=(\hat\lambda,\hat\boldsymbol{\beta}^\top,\hat\sigma^2)^\top$ and the efficiency improved estimator $\hat{\tilde\boldsymbol{\theta}}=(\hat\lambda, \hat{\tilde\boldsymbol{\beta}}^\top,\hat{\tilde\sigma}^2)^\top$. Subsequently, we compare them to the QMLE $\tilde\boldsymbol{\theta}=(\tilde\lambda,\tilde\boldsymbol{\beta}^\top,\tilde\sigma^2)^\top$. It is worth noting that the efficiency improved estimator $\hat{\tilde\boldsymbol{\theta}}$ is a variant of the quasi-score matching estimator, whose calculation relies on obtaining the estimate $\hat\lambda$ in $\hat\boldsymbol{\theta}$. Without causing confusion, we call both $\hat\boldsymbol{\theta}$ and $\hat{\tilde\boldsymbol{\theta}}$ quasi-score matching estimators in the following studies. All the tables in this section, except Table (ref), are provided in Section (ref) of the supplementary material to save space.
In particular, the simulated data are generated from model ((ref)), where $\lambda_0=0.3$, $\boldsymbol{\beta}_0=(2,1)^\top$, $\boldsymbol{X}=(x_{ij})_{n\times 2}$ with $x_{i1}\equiv 1$ and $x_{i2}$ being independently generated from the standard normal distribution for $i=1,\cdots, n$, and the random errors $\varepsilon_i$ are independently generated from the standard normal distribution and the mixture normal distribution $0.9N(0,5/9)+0.1N(0, 5)$, respectively. Hence, $\sigma_0^2=1$. Moreover, we follow the recent literature of network data analysis (see, e.g., zhu2020multivariate) to construct the spatial weights matrix $\boldsymbol{W}=(w_{ij})_{n\times n}$ by the row-normalization of network adjacency matrix $\boldsymbol{A}=(a_{ij})_{n\times n}$, i.e., $w_{ij}=a_{ij}/\sum_{j'=1}^na_{ij'}$, where $a_{ij}$ is either one or zero for $i,j=1,\cdots, n$, and the diagonals $a_{ii}$ are all zeros. In this study, we consider two types of $\boldsymbol{A}$ which are generated from
For the convenience of referring to these two settings of simulation, we name the weights matrices constructed by the simple network $\boldsymbol{A}$ and the SBM $\boldsymbol{A}$, respectively, as the Bernoulli-type and SBM-type weights matrices. {\color{black} It is worth noting that both the Bernoulli-type and SBM-type weights matrices are asymmetric under the above settings.} Lastly, we consider four sample sizes $n=500$, 1,000, 5,000, and 10,000 in the simulation studies.
For each of the above simulation settings, we generate $r=1,\cdots,1000$ replications of simulated data. In the $r$-th replication, we estimate the parameter $ \boldsymbol{\theta}_0=(\lambda_0,\boldsymbol{\beta}_0^\top, \sigma_0^2)^\top=(\theta_{1,0},\cdots,\theta_{4,0})^\top$ by $\hat\boldsymbol{\theta}=(\hat\lambda,\hat\boldsymbol{\beta}^\top, \hat\sigma^2)^\top$, $\hat{\tilde\boldsymbol{\theta}}=(\hat\lambda,\hat{\tilde\boldsymbol{\beta}}^\top, \hat{\tilde\sigma}^2)^\top$, and $\tilde\boldsymbol{\theta}=(\tilde\lambda,\tilde\boldsymbol{\beta}^\top, \tilde\sigma^2)^\top$. For convenience, we denote the estimates in the $r$-th replication by $\hat\boldsymbol{\theta}^{(r)}=(\hat\theta_1^{(r)},\cdots,\hat\theta_4^{(r)})^\top$, $\hat{\tilde\boldsymbol{\theta}}^{(r)}=(\hat{\tilde \theta}_1^{(r)},\cdots,\hat{\tilde \theta}_4^{(r)})^\top$, and $\tilde\boldsymbol{\theta}^{(r)}=(\tilde\theta_1^{(r)},\cdots,\tilde\theta_4^{(r)})^\top$, respectively. To evaluate the performance of these estimates, we consider three measurements: the empirical bias (BIAS) $1000^{-1}\sum_{r=1}^{1000}(\hat\theta^{(r)}_j$ $- \theta_{j,0})$, the empirical standard deviation (SD) $\{1000^{-1}\sum_{r=1}^{1000}(\hat\theta_j^{(r)}- 1000^{-1}\sum_{r=1}^{1000}\hat\theta_j^{(r)})^2\}^{1/2}$, and the root mean squared error (RMSE) $(\textrm{BIAS}^2$ $+\textrm{SD}^2)^{1/2}$, for $\hat\boldsymbol{\theta}$ and $j=1,\cdots, 4$; and similarly for $\hat{\tilde\boldsymbol{\theta}}$ and $\tilde\boldsymbol{\theta}$. To assess the computational complexities of obtaining these estimates, we calculate the average computation times $t_{{\rm S}}=1000^{-1}\sum_{r=1}^{1000}t_{{\rm S}}^{(r)}$ and $t_{{\rm M}}=1000^{-1}\sum_{r=1}^{1000}t_{{\rm M}}^{(r)}$, where for the $r$-th replication, $t_{{\rm S}}^{(r)}$ is the computation time of obtaining the quasi-score matching estimates $\hat\boldsymbol{\theta}$ and $\hat{\tilde\boldsymbol{\theta}}$, and $t_{{\rm M}}^{(r)}$ is the computation time of obtaining the QMLE $\tilde\boldsymbol{\theta}$.
Under the standard normal random errors, Tables (ref) and (ref) present the simulation results for the settings of Bernoulli-type weights matrix and SBM-type weights matrix, respectively, and reveal four important findings. (I) In general, the magnitudes of BIAS and SD become smaller for the two quasi-score matching estimates and the QMLE as $n$ gets larger. It is not surprising that the RMSE shows the same pattern, which indicates the consistency of all three estimators. (II) The difference of SD and RMSE between the quasi-score matching estimate $\hat\boldsymbol{\theta}$ and the QMLE $\tilde\boldsymbol{\theta}$ is small, and gets even smaller as $n$ increases. In Table (ref), for example, the RMSE difference between $\hat\lambda$ and $\tilde\lambda$ is $0.61\times 10^{-2}$ for $n=500$, and is $0.11\times 10^{-2}$ for $n=10,000$. Similarly, the RMSE difference between $\hat\beta_2$ and $\tilde\beta_2$ reduces from $0.14\times 10^{-2}$ to $0.03\times 10^{-2}$ as $n$ increases from 500 to 10,000. (III) {\color{black} The SD and RMSE of the efficiency improved estimates $\hat {\tilde\beta}_1$, $\hat {\tilde\beta}_2$ and $\hat {\tilde\sigma}^2$ are indeed smaller than those of the quasi-score matching estimates $\hat {\beta}_1$, $\hat {\beta}_2$ and $\hat {\sigma}^2$, which confirms the remarks after Theorem (ref).} In addition, the difference of SD and RMSE between the efficiency improved estimate $\hat {\tilde\beta}_2$ and the QMLE $\tilde \beta_2$ is not visible in Tables (ref) and (ref)'s decimal places ($0.01\times 10^{-2}$) under the smallest sample size $n=500$, not to mention under other larger settings of $n$. Similar findings can be found for $\hat {\tilde\sigma}^2$ and $\tilde \sigma^2$. (IV) In terms of the computation time, calculating the quasi-score matching estimates $\hat\boldsymbol{\theta}$ and $\hat{\tilde\boldsymbol{\theta}}$ takes much less than calculating the QMLE $\tilde\boldsymbol{\theta}$. For instance, Table (ref) shows that it takes 1051.27 seconds on average to obtain the QMLE for $n=$ 10,000, while the computation of the quasi-score matching estimates only requires 0.04 seconds (1/26119 of 1051.27). In Table (ref), the average computation time becomes 985.61 seconds for the QMLE, but is only 0.059 seconds (1/16691 of 985.61) for the quasi-score matching estimates. In summary, Tables (ref) and (ref) indicate the huge computation gain of the quasi-score matching estimation for a large $n$, while its efficiency sacrifice is minor, and is negligible for non-intercept regression coefficient $\beta_2$ and variance of random errors $\sigma^2$.
After analyzing the estimation performance under $\lambda_0=0.3$ and the standard normal random errors, we examine the performance of the parameter estimates and computation times under other settings. To save space, these simulation results are relegated to Section (ref) of the supplementary material. Specifically, Tables (ref) -- (ref) in the supplementary material report the results when the random errors are generated from the mixture normal distribution, but $\lambda_0$ remains to be 0.3. In addition, Tables (ref) -- (ref) and Tables (ref) -- (ref) in the supplementary material present the results under settings of $\lambda_0=0.5$ and 0.7, respectively, and the rest of the simulation settings are the same as those in Tables (ref) and (ref). These tables show qualitatively similar findings to those in Tables (ref) and (ref).
In addition to these results, additional simulation studies according to anonymous referees’ suggestions are presented in Sections (ref) -- (ref) of the supplementary material. {\color{black} Specifically, Section (ref) generates the QQ plots for our proposed estimates in simulations and confirms the asymptotic normality in Theorems (ref) -- (ref).} {\color{black} Section (ref) examines the convergence rate $n^{-1}$ in the asymptotic variances of estimators, and demonstrates that the convergence rate of the RMSEs under the Bernoulli-type and SBM-type weights matrices above is approximately $1/\sqrt{n}$. This section also introduces an alternative SBM setting that generates a denser network adjacency matrix $\boldsymbol{A}$, whose resulting weights matrix violates Assumptions (ref) and (ref)$'$, leading to a slower RMSE convergence rate.} {\color{black} Section (ref) shows that the SDs of the quasi-score matching estimates $\hat\boldsymbol{\theta}$ and the efficiency-improved estimates $\hat{\tilde\boldsymbol{\theta}}$, are close to their corresponding theoretical asymptotic standard deviations obtained from Theorems (ref) -- (ref).} {\color{black} Section (ref) demonstrates the computational advantage of the quasi-score matching estimation over the QMLE under a setting of the symmetric weights matrix, where the determinant in the QMLE log-likelihood function is computed using the eigenvalue decomposition algorithm.} {\color{black} Section (ref) compares our proposed quasi-score matching method with the GMM estimation by lee2007gmm, the LSE by huang2019least, the 2SLS estimation by kelejian1998generalized, the best 2SLS estimation by lee2003best and the one-step estimation by gupta2023efficient, and shows that our method is computationally faster than all of these methods, especially for large $n$.}
In this study, we consider the social network experiment conducted in paluck2016changing. The purpose of the experiment was to test the theories that individuals attend to the behavior of certain people in their community, and to understand what is socially normative to adjust individuals' own behavior in their responses. This experiment was conducted over a school year from September 2012 to June 2013, and 28 middle schools were randomly selected to receive an anti-conflict intervention (treatment) via organizing an education workshop about conflicts from 56 middle schools in New Jersey. More details can be found in paluck2016changing.
To measure students' social connections, a survey was conducted at the start of the school year. In this survey, students were asked to nominate up to ten peers with whom they spent most time in the last few weeks. Based on this survey, if student $i$ nominates student $j$, for $i,j=1,\cdots,n$, there is an edge from student $i$ to student $j$ in the social network of $n$ students. Accordingly, $a_{ij}$ in the associated network adjacency matrix $\boldsymbol{A}=(a_{ij})_{n\times n}$ is set to be 1; otherwise, $a_{ij}$ is set to be 0. For completeness, $a_{ii}$ are set to be zeros for $i=1,\cdots,n$. It is worth noting that the nominations are not symmetric, and hence the constructed adjacency matrix $\boldsymbol{A}$ is asymmetric. These students' nominations in the survey result in a realization of the row-normalized adjacency matrix $\boldsymbol{W}=(w_{ij})_{n\times n}$ as the weights matrix in the SAR model ((ref)), where $w_{ij}=a_{ij}/\sum_{j'=1}^na_{ij'}$ for $i,j=1,\cdots,n$. Alongside the edge information between any two students $i,j\in\{1,\cdots, n\}$, there were various questions in a second survey after the experiment about each student $i$'s perception of their school's conflict situations, as well as their demographics. Aside from the surveys, schools' administrative records collected each student $i$'s conflict-related disciplinary events across the school year.
To utilize the SAR model ((ref)) to study the impacts of the anti-conflict intervention and students' social network, on their perception of the conflict situations at school, we construct the response variable $y_i$, $i=1,\cdots, n$, as follow. First, there were 13 questions in the second survey for students' view on their school's conflict situations. Eight of them were positively worded questions. For example, “How many students at this school think it's good to be friendly and nice with all students at this school, no matter who?” The remaining five were negatively worded questions. For instance, “How many students at this school think it's funny, or not a big deal, to post something mean about someone online?" For each of these questions, students were asked to give a score from 0 (“almost nobody”) to 5 (“almost everyone”). Second, we reverse the scores of those negatively worded questions. By doing so, a high score indicates a more friendly and positive school environment consistently for all the 13 questions. Lastly, the response variable $y_i$ is the average score of these 13 questions, which measures the overall feeling of student $i$ towards their school's conflict atmosphere. Preliminary results show that the distribution of $y_i$s is roughly Gaussian.
In addition to the response variable, we construct the regressors $\boldsymbol{x}_i=(1, x_{i2},\cdots, x_{ip})^\top$ for students $i=1,\cdots,n$ given below. We first use a 0/1 binary variable to indicate the treatment, i.e., if student $i$ was from the 28 treated schools, then $x_{i2}=1$; otherwise, $x_{i2}=0$. The other regressors $(x_{i3},\cdots, x_{ip})^\top$ consist of the demographic variables: gender (girl or boy), grade, race (White, Black, Hispanic/Latino, Asian-American, South Asian or Other), whether student $i$ lived with both parents, whether student $i$ was a return student from the previous year, and the number of times that student $i$ was disciplined for a conflict during the school year in which the intervention was conducted.
We then fit the SAR model ((ref)) with the weights matrix $\boldsymbol{W}$ constructed based on the survey conducted at the start of the school year, while we remove the observations with missing values or no network connections in the model fitting, which results in $n=$ 14,732 students in total. We follow Section (ref) to obtain the QMLE $\tilde\boldsymbol{\theta}=(\tilde\lambda,\tilde\boldsymbol{\beta}^\top, \tilde\sigma^2)^\top$, the quasi-score matching estimate $\hat\boldsymbol{\theta}=(\hat\lambda,\hat\boldsymbol{\beta}^\top, \hat\sigma^2)^\top$ and the efficiency improved estimate $\hat{\tilde\boldsymbol{\theta}}=(\hat\lambda,\hat{\tilde\boldsymbol{\beta}}^\top, \hat{\tilde\sigma}^2)^\top$. The estimation results are summarized in Table (ref). Overall, the three estimates and their standard errors are similar. In particular, the efficiency improved estimate $\hat{\tilde \boldsymbol{\theta}}$ is numerically closer to the QMLE $\tilde\boldsymbol{\theta}$ compared to the quasi-score matching estimate $\hat\boldsymbol{\theta}$. Moreover, the standard errors of $\hat{\tilde\boldsymbol{\beta}}$ and $ \hat{\tilde\sigma}^2$ are almost identical to those of $\tilde\boldsymbol{\beta}$ and $\tilde\sigma^2$, which is consistent with the findings in the simulation studies.
In terms of testing the significance of parameters, the three estimation approaches lead to the same conclusions for parameters $\lambda$, $\sigma^2$ and $\boldsymbol{\beta}$ at the significance level 0.10, but with two exceptions for $\boldsymbol{\beta}$ at the significance level 0.05: the coefficients of gender and Asian-American. In these two exceptions, the QMLE and the efficiency improved estimate both imply the same conclusions, i.e., the coefficient of gender is not significant and the coefficient of Asian-American is significant. On the other hand, the quasi-score matching estimate gives different conclusions based on the $p$-values. This finding is sensible since the efficiency improved estimator has a smaller estimation error for $\boldsymbol{\beta}$ compared to the quasi-score matching estimator, which is demonstrated in the simulation studies of Section (ref). It is worth noting that the efficiency improved estimator is designed to alleviate the computation for a large $n$, though in this study we also calculate the QMLE as a benchmark. Based on the above comparisons, we recommend using the efficiency improved estimate in practice for fitting the SAR model with a large $n$.
We subsequently focus on interpreting the results of the efficiency improved estimate $\hat{\tilde \boldsymbol{\theta}}$ in Table (ref), which show seven interesting findings at the significance level 0.05. First, the positive significance of the spatial autoregressive coefficient $\lambda$ indicates that the social network effects exist and the students' perception of their school's conflict situations highly depends on the view of their connected peers. Second, the treatment is significantly negative. This is sensible because students became more aware of their school's conflicts than what they used to be due to the organized education workshop about conflicts. This result is consistent with the findings in paluck2016changing, in which they find that the introduction of an anti-conflict intervention increases the community's attention to conflicts. Third, the significant negative effect of students' grade shows that senior students tend to perceive more conflicts than junior students. Fourth, it can be seen that White and Asian-American students have more significant positive ratings on their school's atmosphere, while Black and Hispanic/Latino students have more negative views, especially for Black students based on the strongly negative coefficient of this group. Such a finding may be due to the potential racial discrimination at school. Fifth, students living with both parents tend to have higher ratings on their school's atmosphere; this is probably because of the good impacts from their family relationships. Sixth, it turns out that the more time students are disciplined, the more negative perception they have toward their school's atmosphere. Lastly, the other variables including gender and whether the student was a return student from the previous year are not significant.
For comparison, we also fit the SAR model ((ref)) with $\boldsymbol{W}$ constructed based on another survey conducted at the end of the school year. To save space, the fitting results and their associated findings are provided in Table (ref) and Section (ref) of the supplementary material. {\color{black} As recommended by an anonymous referee, we also present the results of using the 2SLS estimation by kelejian1998generalized and the best 2SLS estimation by lee2003best in Table (ref) of the supplementary material. Compared to the efficiency improved estimate $\hat{\tilde \boldsymbol{\theta}}$ in Tables (ref) and (ref), we find that both 2SLS estimates yield consistent results in terms of signs and significance. Moreover, the efficiency improved estimates exhibit smaller standard errors and are notably closer to those of the QMLE.}
This paper proposes a new quasi-score matching estimation approach for the SAR model. Based on our understanding, this article is the first to introduce score matching for spatially dependent data, unlike its prior use for independent and identically distributed data in machine learning (see, e.g., hyvarinen2005estimation). In addition, we propose an efficiency improved estimator that further approaches the efficiency of the QMLE. Both theoretical and numerical studies show that as $n$ is large, the efficiency differences among our two proposed estimators and the QMLE are small, while our estimators offer significantly simpler and faster computations. To analyze the asymptotic properties of these estimators under random weights matrix $\boldsymbol{W}$ and regressors $\boldsymbol{X}$, we develop a pioneering theoretical framework, including a new LLN and CLT. Both the quasi-score matching estimation approach and the proposed theoretical framework can broad applications in many other spatial econometrics models (see, e.g., lee2010estimation, zou2021network and wu2021inward) for future research, in particularly when both $\boldsymbol{W}$ and $\boldsymbol{X}$ are random.
{\color{black} In addition to the asymptotic theory developed in this paper for random $\boldsymbol{W}$ and $\boldsymbol{X}$, as suggested by an anonymous reviewer, a promising direction for future research is to extend the theory of xu2015maximum and wu2022applications to allow for random $\boldsymbol{W}$. This would involve establishing near-epoch or spatial functional dependence properties for $y_i$ and related variables under model ((ref)), thereby generalizing, for example, Proposition 1 in xu2015maximum or Proposition 4.3 in wu2022applications to the random $\boldsymbol{W}$ setting.} We believe these efforts would strengthen the usefulness of the theoretical framework for random weights matrix $\boldsymbol{W}$ and regressors $\boldsymbol{X}$.
\spacingset{0.97} {4pt plus 0.3ex}