Privacy-Protected Spatial Autoregressive Model

\csection{INTRODUCTION}

A network refers to a set of nodes and their observed relationships (i.e., edges), and network data refers to the information collected from a network. Network data from different nodes are likely to be dependent on each other because different nodes in the networks are connected to each other by edges. This is arguably the most important type of dependence induced by a network structure. For convenience, we refer to it as network dependence. To study network dependence, the spatial autoregressive (SAR, ord1975estimation) model has been considered an extremely useful tool. The key idea of the SAR model is to assume that the response collected from each node is linearly related to its connected neighbors. This simple mechanism introduces a sophisticated and elegant network-dependency relationship, which leads to the frequent use of the SAR model and its variants in real practice anselin2008spatial_typical,chen2013impact,zhu2017network,beenstock2019econometric.

In fact, the classical SAR models have inspired numerous follow-up studies. For example, lee2010estimation introduced the SAR panel model with individual effects and SAR disturbances. yang2016 and zhu2020multivariate developed a multivariate SAR model that models multivariate responses collected from each node. Various partially linear semiparametric spatial models were also developed to allow the SAR parameter to meaningfully vary across different nodes Su2012,Malikov2017. The spatial dynamic panel data models were proposed to model time and network dependence simultaneously yu2008quasi,lee2014efficient,li2017fixed. To compute the maximum likelihood estimator for various SAR models, the determinant of a large-scale matrix must be computed. This creates a significantly high level of computational complexity. To alleviate the computational cost, various efficient computing algorithms were developed pace2000method,chen2013impact,zhou2017. In particular, the least squares method of huang2019 and zhu2020multivariate has been demonstrated to be practically useful.

Notably, the aforementioned studies assume that network data can be observed accurately. This is true for many real-world applications where data privacy is not a major concern. However, in many situations, privacy protection has become increasingly important. Consider a data company that seeks collaboration with researchers for data analysis. While the company aims to achieve precise analysis results, it also wishes to protect the privacy of the data on its platform. For instance, a third-party payment platform may need researchers to analyze network effects among merchants on their platform. Network links, typically defined by merchants located in the same commercial district or in close proximity, can be directly observed and do not require privacy protection. However, the platform does not want to disclose specific business information of the merchants, such as turnover or average transaction value. In such cases, data sharing with researchers can only occur if there is a robust privacy protection mechanism in place. Similar scenarios are common in various fields, including medical research collaborations, transaction data analysis, and social media studies. To protect individual privacy, appropriate measures have been developed and implemented to safeguard raw data.

These privacy protection measures include, but are not limited to, data swapping Reiter2005, data imputation Raghunathan2003,Reiter2005, posterior distribution sampling Hu2018,Wilde2021, and noise addition Dwork2006,Wilde2021,Ito2021. Among these measures, noise addition has arguably gained great popularity. Its popularity may be due to two reasons. First, noise addition is a simple method and can be easily implemented in practice. Second, different levels of privacy protection can be provided by specifying different noise levels for the added random noise. In fact, the relationship between the added noise level and privacy-protection strength can be analytically studied using the theory of differential privacy Dwork2014. However, extending this useful idea to privacy-protected spatial data analysis becomes a challenge. When data companies share data, they aim to obtain accurate estimates while also protecting sensitive variables, whether the covariates \(X\) or the response \(Y\). This necessitates adding noise to both \(X\) and \(Y\). Therefore, we start with a classical SAR model and consider protecting the raw response and some explanatory variables that require privacy-protection by adding artificially generated random noise.

This issue is related to measurement errors in variables in traditional statistical research, where measurement errors in covariates have received considerable attention Fuller1987,Stefanski1987,Staudenmayer2005,Carroll2006,Buonaccorsi2010,Hausman2021. Tools to deal with measurement error include, but are not limited to, simulation-extrapolation cook1994simulation, Novick2002, bertrand2017inference, instrumental variables schennach2007instrumental,gustafson2007measurement,jiang2020measurement,Luo2022, and the corrected score method Stefanski1987,Nakamura1990,tsiatis2004locally,wang2012corrected. In the analysis of SAR models, scholars have also focused on addressing covariate measurement error. For example, Li2009 adopted structural model-based approaches to integrate out \(X\) when the distribution of \(X\)'s measurement error is known. Luo2022 proposed the 3SLS method which uses instrumental variables to deal with measurement error in \(X\). However, current literature on addressing measurement error in SAR models rarely considers scenarios where both \(X\) and \(Y\) contain noise, a common situation in privacy-preserving data analysis collaborations with data protection requirements. This might seem trivial in linear models because measurement errors in the continuous variable \(Y\) do not affect the unbiasedness of the parameter estimates Hausman2001. In SAR models, however, adding artificial noise to \(Y\) inevitably leads to biased parameter estimates. Furthermore, in data collaborations, the sample size is typically large. However, existing literature on SAR with measurement errors seldom considers the computational complexity. Therefore, developing accurate and computationally efficient parameter estimation methods for SAR models in large-scale datasets for both noisy \(X\) and \(Y\) is of great importance in real practice.

In this study, we propose a privacy-protected SAR (PSAR) model to meet the privacy-protection requirements. To obtain more accurate estimates, we assume that the data company has predetermined and disclosed the level of noise. Unlike previous SAR models, the unavailability of true values for both the response and some explanatory variables makes the classical quasi-maximum likelihood estimator (QMLE) inapplicable. To address this challenge, inspired by the corrected score method Stefanski1987,Nakamura1990, we first derive the likelihood function based on the true explanatory variables and observed responses. This approach allows for an explicit expression of the log-likelihood function and the derivation of its first- and second-order derivatives, which are, however, biased. Nonetheless, the bias can be analytically expressed under appropriate regularity conditions and subsequently corrected. Subsequently, a Newton-Raphson-type algorithm can be designed to obtain a corrected likelihood estimator (CLE). Nonetheless, computing the CLE becomes highly challenging for large-scale networks due to the computational complexity associated with high-dimensional matrix operations, such as determinant and inverse computations. To address this issue, we propose a corrected least squares estimator (CLS) for the PSAR model inspired by huang2019, zhu2020multivariate, and the concept of bias correction. Additionally, the statistical inference of the proposed method has been carefully discussed. In the theoretical analysis, we establish the identifiability, asymptotic normality, and numerical convergence properties for the estimators. Finally, we expand our analysis to encompass more complex scenarios. These include extending the classical SAR model to allow for multiple responses and multiple classes of node types in the networks, as well as incorporating perturbations to the network structure $W$ to enhance privacy protection.

The remainder of this paper is organized as follows. Section 2 describes the estimation methods in detail. Section 3 presents the theoretical properties from both statistical and numerical perspectives and provides a discussion of generalization for these methods. Section 4 presents the numerical results, including those of simulation studies and a real data example. Section 5 presents the concluding remarks. All technical details are relegated to the supplementary material.

\csection{PRIVACY-PROTECTED SPATIAL AUTOREGRESSIVE MODEL AND ESTIMATORS}

Consider a network with a total of $N$ nodes. To describe the network structure, we define an adjacency matrix $A=(a_{i_1i_2})\in \{0,1\}^{N\times N}$ ($1\leq i_1,i_2\leq N$). Each element $a_{i_1i_2}=1$ if there exists an edge from node $i_1$ to $i_2$ ($i_1\neq i_2$); otherwise, $a_{i_1i_2}=0$. We assume $a_{ii}=0$ for $1\leq i\leq N$. We define the weighting matrix $W =(w_{ij}) \in \mathbb{R}^{N \times N}$ with element $w_{ij} = a_{ij}/d_i$, where $d_i =\sum_{j=1}^N a_{ij} $ represents the nodal out-degree of node $i$. Throughout the remainder of this paper, we assume $d_i>0$ for every $1\leq i\leq N$. Otherwise, nodes with $d_i=0$ should not be included for analysis.

Next, let $Y_{i}$ be the response collected from the $i$-th node ($1\leq i\leq N$) and $\mathbf{x}_{i}\in\mathbb{R}^p$ be the associated covariate vector. Write $\mathbb{Y} = ({Y}_{1}, \cdots,{Y}_{N})^\top\in \mathbb{R}^{N}$ as the response vector, $\mathbb{X}=(\mathbf{x}_{1},\cdots,\mathbf{x}_{N})^\top\in \mathbb{R}^{N\times p}$ as the associated covariate matrix representing all explanatory variables, and $\mathbb{E}=(e_{1}, \cdots,e_{N})^\top\in \mathbb{R}^{N}$ as the error vector, where $e_{i}$s are independent and identically distributed with mean 0 and covariance $\sigma^2_0$. If all the response and explanatory variables can be faithfully observed, we consider the following SAR model ord1975estimation,anselin2008spatial_typical,

where $\rho_0 \in \mathbb{R}$ is the network autocorrelation parameter measuring the network effect, and $\beta_0 \in \mathbb{R}^p$ is the regression-coefficient vector reflecting the effect of nodal covariates. We define $S_0=I_N-\rho_0 W$, where $I_N\in\mathbb{R}^{N\times N}$ is the identity matrix. We then have $\mathbb{Y} = S_0^{-1}(\mathbb{X}\beta_0+\mathbb{E})$. By omitting the constant term, the log-likelihood function for the parameter of interest $\theta=(\rho,\beta^\top,\sigma^2)^\top\in\mathbb{R}^{p+2}$ can be obtained easily. Accordingly, a QMLE for $\theta$ can be obtained as $ \hat{\theta}=\arg\max_\theta \ell(\theta,\mathbb{Y}, \mathbb{X})$. Its asymptotic properties have been well studied. For example, refer to Lee:2004 and anselin2013spatial.

However, in a privacy-protection scenario, we assume that the response $\mathbb{Y}$ and part of the explanatory covariates cannot be faithfully observed. Specifically, for the covariates, we assume that $\mathbb{X}=(\mathbb{X}_1,\mathbb{X}_2)$ with $\mathbb{X}_1\in\mathbb{R}^{N\times p_1}$, $\mathbb{X}_2\in\mathbb{R}^{N\times p_2}$, and $p_1+p_2=p$. Here, we define $\mathbb{X}_1$ to be the matrix collecting the faithfully observed covariates, and $ \mathbb{X}_2^*=\mathbb{X}_2+\mathcal{E}_x\in\mathbb{R}^{N\times p_2}$ to be the matrix collecting explanatory variables with artificially created error matrix $ \mathcal{E}_x=(\varepsilon_{x,ij})\in\mathbb{R}^{N\times p_2}$. Further, we assume that each element $ \varepsilon_{x,ij}$ ($ 1\leq i\leq N,1 \leq j\leq p_2$) in $\mathcal{E}_x$ is independently and identically distributed with mean 0 and predefined variance $\lambda_x^2$. Moreover, for the response, a data user can only observe $\mathbb{Y}^*=\mathbb{Y}+\mathcal{E}$ with an artificially created noise vector $\mathcal{E}=(\varepsilon_{1},\cdots,\varepsilon_{N})^\top\in\mathbb{R}^N$. Here, we assume that different $\varepsilon_{i}$s are independent and identically distributed with mean 0 and known variance $\lambda^2$. We assume here that data users cannot observe the true response variables $\mathbb{Y}$ and some explanatory variables $\mathbb{X}_2$. However, to better utilize the data, they are aware of the privacy-protection parameters $\lambda^2$ and $\lambda_x^2$. This assumption ensures the privacy of the data while simultaneously facilitating better estimation results for data users. Notably, if \(\lambda^2\) is unknown, an identification issue arises for the parameter when \(\rho = 0\). However, since \(\lambda^2\) is assumed to be known, this identification problem is eliminated.

Based on the above notations, we can derive the expression for $\mathbb{Y}^*$ based on the observed $\mathbb{X}^*=(\mathbb{X}_1,\mathbb{X}^*_2)$ as,

where $\beta_{01}\in\mathbb{R}^{p_1}$ represents the coefficient corresponding to the truly observed covariates, $\beta_{02}\in\mathbb{R}^{p_2}$ is the coefficient corresponding to those with noise added, and we have $\beta_0=(\beta_{01}^\top,\beta_{02}^\top)^\top$. For convenience, we refer to this model in (ref) as the privacy-protected spatial autoregressive (PSAR) model. Notably, model (ref) only represents the relationship between the observable $\mathbb{Y}^*$ and $\mathbb{X}^*$. However, the challenge here is that $\mathbb{Y}^*$ is generated by true $\mathbb{X}$, which is not fully observable to data users.

To accomplish model estimation, the foremost method worth considering is QMLE. However, note that the likelihood function for (ref) is hard to be spelled out because $\mathbb{Y}^*$ is generated according to the true $\mathbb{X}_2$. Unfortunately, as mentioned above, $\mathbb{X}_2$ cannot be observed. Consequently, the classical QMLE method is not applicable in this case. To address this issue, we begin the analysis by assuming that $\mathbb{X}_2$ is known, and we consider the underlying generating process of the observed $\mathbb{Y}^*$. Therefore, we have,

By (ref), the negative log-likelihood function for $\theta$ based on $\mathbb{Y}^*$ and $\mathbb{X}$ is,

where $\Omega=\sigma^2 I_N+\lambda^2SS^{\top}$. Ideally, we can apply the classical Newton-Raphson algorithm to solve the optimization problem $\arg\min_\theta \mathcal L(\theta)$, and obtain estimates with desirable properties. Unfortunately, equation ((ref)) is computationally infeasible in real practice because we can only observe $\mathbb{X}_2^*$. Direct substitution of $\mathbb{X}_2$ with $\mathbb{X}_2^*$ in ((ref)) to obtain $\mathcal L^*(\theta)=\mathcal L(\theta,\mathbb{Y}^*,\mathbb{X}^*)$ for feasible computation will inevitably result in biased estimation. This is because the key statistical estimators derived from $\mathcal L^*(\theta)$, including the first- and second-order derivatives, are all biased.

We first analyze the difference between the first-order derivatives of $\mathcal L(\theta)$ and $\mathcal L^*(\theta)$. Define $ \dot\mathcal L(\theta)$ and $ \dot\mathcal L^*(\theta)$ to be the first-order derivative of $\mathcal L(\cdot)$ and $\mathcal L^*(\cdot)$ with respect to $\theta$, respectively. Then, we can calculate that $E\big\{ \dot{\mathcal L}(\theta_0) \big\}=\mathbf{0}_{p+2}$ and

where $\dot\mathcal L^{*\rho}(\theta)$, $\dot\mathcal L^{*\beta}(\theta)$, and $\dot\mathcal L^{*\sigma^2}(\theta)$ are the first-order derivatives of $\mathcal L(\theta)$ with respect to $\rho$, $\beta$, and $\sigma^2$, respectively. See Appendix A.1 for detailed verifications. We define this difference for the first-order derivative caused by $\mathbb{X}^*$ as $\Delta\mathbf{S}_{\rm CL}(\theta_0)=E\big\{\dot\mathcal L^*(\theta_0)\big\}-E\big\{\dot\mathcal L(\theta_0)\big\}$.

Similarly, define the second-order derivative for $\mathcal L(\theta)$ and $\mathcal L^*(\theta)$ as $\ddot\mathcal L(\theta)$ and $\ddot\mathcal L^*(\theta)$, respectively. Define $\mathbb{W}_S=WS^{\top}+SW^{\top}$. In this way, the difference between the expectations of second-order derivatives for $\mathcal L^*(\theta_0)$ and $\mathcal L(\theta_0)$ could be calculated and defined as,

It could be verified that $\Delta\mathbf{H}_{\rm CL}(\theta_0)$ is symmetric and

See Appendix A.1 for detailed verifications. Based on the explicitly expressed bias in the derivatives led by $\mathbb{X}_2^*$, a natural method to estimate parameter $\theta_0$ is to reconstruct an estimator using the bias-corrected first- and second-order derivatives. We refer to this estimator as {\it the corrected likelihood estimator} (CLE) , which is denoted as $\hat\theta_{\rm CLE}$. Accordingly, an iterative algorithm can be developed. At each iteration, we conduct three steps: (1) (Calculating) compute the first- and second-order derivatives based on the observed $\mathbb{X}^*$ and $\mathbb{Y}^*$; (2) (Debiasing) perform bias correction for the derivatives; and (3) (Updating) conduct a Newton-Raphson-type iterative formula based on the corrected derivatives.

To be more specific, let $\hat\theta_{\rm CLE}^{(0)}=(\hat\rho_{\rm CLE}^{(0)},\hat\beta_{\rm CLE}^{(0)\top},\hat\sigma^{2(0)}_{\rm CLE})^\top\in\mathbb{R}^{p+2}$ be the initial estimator, which could be, for example, the QMLE estimator assuming that $\mathbb{X}_2^*$ as the true $\mathbb{X}_2$. Let $\hat\theta_{\rm CLE}^{(t)}=(\hat\rho_{\rm CLE}^{(t)},\hat\beta_{\rm CLE}^{(t)\top},\hat\sigma^{2(t)}_{\rm CLE})^\top\in\mathbb{R}^{p+2}$ be the estimator obtained in the $t$-th iteration, and $\hat S^{(t)}$, $\hat\Omega^{(t)}$ be $S$, $\Omega$ with the plugged-in estimator $\hat\theta_{\rm CLE}^{(t)}$. We next consider how to update $\hat\theta_{\rm CLE}^{(t)}$ to $\hat\theta_{\rm CLE}^{(t+1)}$. We start from the parameter $\gamma_0=(\rho_0,\beta_0^\top)^\top\in\mathbb{R}^{p+1}$ and its update $\hat\gamma_{\rm CLE}^{(t+1)}=(\hat\rho_{\rm CLE}^{(t+1)},\hat\beta_{\rm CLE}^{(t+1)\top})^\top$. It could be calculated that, for the first-order derivative, $\Delta\textbf{S}_{\rm CL}^{(t)} =\Big(\lambda^2\lambda_x^2\hat\beta_{2}^{(t)\top}\hat\beta_{2}^{(t)}\mbox{tr}\Big\{(\hat\Omega^{(t)})^{-1}W\hat S^{(t)\top}(\hat\Omega^{(t)})^{-1}\Big\}, {\bf0}_{p_1}^\top, \lambda_x^2\mbox{tr}\big\{(\hat\Omega^{(t)})^{-1}\big\} \hat\beta^{(t)\top}_{2}\Big)\in\mathbb{R}^{p+1}$. Similarly, define $\Delta\textbf{H}_{\rm CL}^{(t)} =\Big(\mathbf{H}_{\rm CL}^{\rho\rho(t)} ,\mathbf{H}_{\rm CL}^{\rho\beta(t)\top};\mathbf{H}_{\rm CL}^{\rho\beta(t)},\mathbf{H}_{\rm CL}^{\beta\beta(t)} \Big)\in\mathbb{R}^{(p+1)\times(p+1)} $, which could be calculated with plugged in $\hat\theta_{\rm CLE}^{(t)}$ correspondingly. Then, $\hat\gamma_{\rm CLE}^{(t+1)}$ can be obtained using a corrected Newton-Raphson-type algorithm as,

where $\dot\mathcal L^{*\gamma}(\hat\theta_{\rm CLE}^{(t)})=\Big(\dot\mathcal L^{*\rho}(\hat\theta_{\rm CLE}^{(t)}),\dot\mathcal L^{*\beta\top}(\hat\theta_{\rm CLE}^{(t)})\Big)^\top$ and $\ddot\mathcal L^{*\gamma}(\hat\theta_{\rm CLE}^{(t)})=\Big(\ddot\mathcal L^{*\rho\rho}(\hat\theta_{\rm CLE}^{(t)}),\ddot\mathcal L^{*\rho\beta\top}(\hat\theta_{\rm CLE}^{(t)});\\ \ddot\mathcal L^{*\rho\beta}(\hat\theta_{\rm CLE}^{(t)}),\ddot\mathcal L^{*\beta\beta}(\hat\theta_{\rm CLE}^{(t)})\Big)$. Then, in the $(t+1)$-th iteration, the estimate of $\sigma_0^2$ can be obtained by solving (ref), which is $\hat\sigma_{\rm CLE}^{2(t+1)}$. See Algorithm (ref) for details. This leads to the final estimator $\hat\theta_{\rm CLE}^{(t)}$. We establish the numerical convergence of $\hat\theta_{\rm CLE}^{(t)}$ in the next section.

The proposed CLE offers a feasible parameter-estimation method for the PSAR model. However, the calculation of the CLE involves the determinants of $S=I_N-\rho W$ and $\Omega=\sigma^2I_N+\lambda^2 SS^{\top}$. This makes the estimation process computationally expensive for large-scale networks with complexity of order $O(N^3)$. Consequently, inspired by the least squares estimation method by huang2019 and zhu2020multivariate, we propose a corrected least squares estimator for the PSAR model as follows.

To illustrate the idea, we start with model (ref) based on the true response $\mathbb{Y}$ and covariates $\mathbb{X}_2$. Consider the conditional expectation of $Y_i$, given the responses of all other nodes. Define $\mathcal F_{(-i)}=\sigma\{(X_{i'},Y_{i'}):i'\neq i\}$ to be the $\sigma$-field generated by all other nodes. Then, under the assumption of a normally distributed ${\mathbb E}$, we have, $E\{Y_{i}|\mathcal F_{(-i)}\}=\mu_{i}+\sum_{j\neq i}^{N}\alpha_{ij}(Y_{j}-\mu_{j})$, where $\mu_{i}=E(Y_i)$ and

Define $d_\rho=\mbox{diag}^{-1}(S^\top S)$ and recall that $\gamma=(\rho,\beta^\top)^\top$. Accordingly, a least-squares-type objective function can be constructed as,

Then, a least squares estimator can be obtained as $\hat\gamma_{\rm LS}=\mbox{argmin} \mathcal L_{\rm LS} (\gamma)$.

We make two remarks about this method. First, from equation (ref), we see that, for node $i$, only directly connected nodes (i.e., $a_{ij}+a_{ji}\neq 0$) and certain types of indirectly connected nodes with $\sum_{k}a_{ki}a_{kj}\neq 0$ are considered for the calculation of the conditional mean. These nodes will then further be involved in the computation of the objective function (ref). Thus, as long as the network is sufficiently sparse, the computation of the objective function will be efficient. Second, the construction of the objective function here is inspired by the assumption that ${\mathbb E}$ follows a normal distribution. However, the formal statistical analysis does not rely on the normal assumption. Specific assumptions of the error term are provided in the next section.

However, the key challenge here is that we can only observe the noise-added responses $\mathbb{Y}^*$ and covariates $\mathbb{X}_2^*$. This leads to the noise-added loss function $\mathcal L_{\rm LS}^*(\gamma)=\mathcal L_{\rm LS}(\gamma;\mathbb{Y}^*,\mathbb{X}^*)$ instead of $\mathcal L_{\rm LS}(\gamma)$. Similar to the analysis of the CLE, we consider a corrected least squares method, while simultaneously correcting the biases introduced by the observational errors of $\mathbb{X}^*$ and $\mathbb{Y}^*$ on both the first- and second-order derivatives. This results in the {\it corrected least squares estimator} (CLS) $\hat\gamma_{\rm CLS}$.

Specifically, let $\dot \mathcal L^*_{\rm LS}(\gamma) = \big(\dot \mathcal L^{*\rho}_{\rm LS}(\gamma),\dot \mathcal L^{*\beta}_{\rm LS}(\gamma)^\top\big)^\top\in\mathbb{R}^{p+1}$ be the first-order derivative of $\mathcal L_{\rm LS}(\gamma)$ with respect to $\gamma$. Here, $\dot \mathcal L^{*\rho}_{\rm LS}(\gamma)$ and $\dot \mathcal L^{*\beta}_{\rm LS}(\gamma)$ represent the first-order derivatives of $\mathcal L_{\rm LS}(\gamma)$ with respect to $\rho$ and $\beta$, respectively. Define $\dot d_\rho$ and $\ddot d_\rho$ to be the first- and second-order derivatives of $d_\rho$ with respect to $\rho$, respectively. Define $ \mathbb{W}=W^{\top}S+S^{\top}W$ and $\mathbb{V}^*=S\mathbb{Y}^*-\mathbb{X}^*\beta$. It can be verified that $E\big\{\dot\mathcal L_{\rm LS}(\gamma_0)\big\}=\mathbf{0}_{p+1}$, and

The verification details of equations (ref) and (ref) are given in Appendix A.2. Define the right-hand side of (ref)-(ref) as $ \Delta \textbf{S}_{\rm LS}(\gamma_0)=E\big\{\dot \mathcal L_{\rm LS}^*(\gamma_0)\big\}$. Consequently, the estimator obtained by directly minimizing $\mathcal L_{\rm LS}^*(\gamma)$ cannot be consistent.

To propose a Newton-Raphson-type algorithm, we still need to discuss the second-order derivatives. Define $\ddot\mathcal L_{\rm LS}(\gamma)=\big(\ddot\mathcal L_{\rm LS}^{\rho\rho}(\gamma),\ddot\mathcal L_{\rm LS}^{\rho\beta}(\gamma)^\top;\ddot\mathcal L_{\rm LS}^{\rho\beta}(\gamma),\ddot\mathcal L_{\rm LS}^{\beta\beta}(\gamma)\big)$ to be the second-order derivative of $\mathcal L_{\rm LS}(\gamma)$. We can further define the difference between the expectations of second-order derivatives for $\mathcal L^*_{\rm LS}(\gamma)$ and $\mathcal L_{\rm LS}(\gamma)$ as,

As a result, we can verify that $\Delta\mathbf{H}_{\rm LS}(\gamma_0)$ is symmetric and

The verification details of (ref) and the expressions for it are given in Appendix A.2.

An iterative algorithm (i.e., Algorithm (ref)) can then be developed, which can also be described in the three steps (calculating, debiasing, and updating). Here, we only list the updating equation for simplicity. Let $\hat\gamma_{\rm CLS}^{(0)}=(\hat\rho_{\rm CLS}^{(0)},\hat\beta_{\rm CLS}^{(0)\top})^\top\in\mathbb{R}^{p+1}$ be the initial estimator and $\hat\gamma_{\rm CLS}^{(t)}=(\hat\rho_{\rm CLS}^{(t)},\hat\beta_{\rm CLS}^{(t)\top})^\top\in\mathbb{R}^{p+1}$ be the estimator obtained in the $t$-th iteration. The estimator in the ($t+1$)-th iteration can then be updated as follows,

where $\Delta\textbf{S}_{\rm LS}^{(t)}=\Delta\textbf{S}_{\rm LS}\Big(\hat\gamma^{(t)}_{\rm CLS}\Big)$ and $\Delta\textbf{H}_{\rm LS}^{(t)}=\Delta\textbf{H}_{\rm LS}\Big(\hat\gamma^{(t)}_{\rm CLS}\Big)$. By the time of convergence, we obtain the final estimator. It will be shown in the next section that $\hat\gamma_{\rm CLS}^{(t)}$ numerically converges to $\hat\gamma_{\rm CLS}$. This is the second estimator we have developed in this work.

In terms of the computational advantage, the CLS is observed to avoid the need for large-scale matrix-inversion calculations throughout the entire algorithm. This includes computing the first-order derivative, second-order derivative, and correction terms in equations (ref)--(ref). Consequently, this approach significantly reduces the computational time required compared to the CLE. The substantial computational advantages of $\hat\gamma_{\rm CLS}$ are further illustrated in the numerical analysis. \csection{THEORETICAL PROPERTIES}

In this section, we first theoretically analyze both the numerical and statistical properties of the CLE and CLS. Then, the generalizations of PSAR are carefully presented to extend the method to more model forms, enabling support for a wider range of data analysis applications. First, we introduce some theoretical assumptions.

Define $\|B\|_1=\max_j\sum_i|b_{ij}|$ to be $L_1$ norm and $\|B\|_\infty=\max_i\sum_j|b_{ij}|$ to be $L_\infty$ norm for an arbitrary matrix $B=(b_{ij})$. For theoretical analysis, the following regularity conditions are required.

(C1) (Covariates) Assume that $\lim_{N\to\infty} N^{-1} \mathbb{X}^\top \mathbb{X}$ exists and is non-singular.

(C2) (Uniform Boundedness) The elements of $\mathbb{X}$ are uniformly bounded for all $N$. Moreover, the elements $w_{ij}$ of $W$ are uniformly bounded with the uniform order $O(1/h_n)$, where $h_n$ can be bounded or divergent with $h_n/N \to 0$ as $N \to \infty$. Furthermore, $W$ and $S_0^{-1}$ have bounded $L_1$ and $L_\infty$ norms.

(C3) (Noise Term) Assume that all the ${e_i}$s, $\varepsilon_i$s, and $\varepsilon _{ij}^x$s ($1\leq i,j\leq N$) are independently and identically distributed with mean zero. For the variance, we assume that $\text{Var}(\varepsilon_i)=\lambda^2$ and $\text{Var}(\varepsilon_{ij}^x)=\lambda_x^2$. To simplify the asymptotic covariance form of the estimators, we assume $E(e_i^3)=E(\varepsilon_i^3)=E(\varepsilon_{ij}^{x3})=0$. For higher moment conditions, we assume that $E(e_i^4)=\mu_4^e$, $E(\varepsilon_i^4)=\mu_4^\varepsilon$, $E(\varepsilon_{ij}^{x4})=\mu_4^{\varepsilon_x}$, $E(e_i^2-\sigma_0^2)^4=c'_e$, $E(\varepsilon_i^2-\lambda^2)^4=c'_\varepsilon$, and $E(\varepsilon_{ij}^{x2}-\lambda_x^2)^4=c'_{\varepsilon_x}$, for positive constants $\mu_4^e$, $\mu_4^\varepsilon$, $\mu_4^{\varepsilon_x}$, $c'_e,$ $c'_\varepsilon$, and $c'_{\varepsilon_x}$.

Condition (C1) requires the covariance to exist and be non-singular. The same condition was adopted by Lee:2004, yang2016, and zhu2020multivariate. Condition (C2) requires the uniform boundedness of features $\mathbb{X}$, the weighting matrix $W$, and $S_0^{-1}$; this is a classical regularity condition in the SAR model Lee:2004,yang2016. Condition (C3) is a typical assumption for the noise term huang2019, zhu2020multivariate, huang2020two. It is remarkable that the artificially created noise is allowed to be non-Gaussian, as long as the moment conditions are satisfied. Furthermore, condition $E(e_i^3)=E(\varepsilon_i^3)=E(\varepsilon_{ij}^{x3})=0$ is assumed for convenience. The aim is to simplify the form of the asymptotic covariances for the estimators. Notably, the theory to be presented can be softly generalized under the assumption of non-zero third-order moments using the same theoretical proof techniques in this work. Given these conditions, we then establish the theoretical results in the subsequent subsections.

Before establishing the theoretical properties of the CLE, we first address the challenges in the theoretical analysis of the PSAR model compared to the traditional SAR models without privacy protection. Next, we discuss the identification issue of $\hat\theta_{\text{CLE}}$. The consistency and asymptotic normality of $\hat\theta_{\text{CLE}}$ are then carefully established. Finally, the numerical convergence of the proposed iterative algorithm is rigorously proved.

In the theoretical analysis of the traditional SAR estimators based on the likelihood function, the proof core involves handling of the matrix $S^{-1}=(I_N-\rho W)^{-1}$. The existing literature presents two primary approaches to address this issue. The first method considers bounding both the $L_1$ norm and $L_\infty$ norm of $S^{-1}$. This can be verified by assuming that the row and column sums of $W$ and $S^{-1}$ are uniformly bounded; for example, refer to Lee:2004 and yang2016. Meanwhile, the second approach considers the application of Taylor's expansion to $S^{-1}$, based on the assumption that the weighting matrix $W$ is a transition-probability matrix, and the stationary distribution for the network nodes exists; for example, refer to huang2019 and zhu2020multivariate.

However, to estimate the parameters for the PSAR model based on the likelihood function, we need to analyze $\Omega^{-1}=\big(\sigma^2I_N+\lambda^2 SS^{\top}\big)^{-1}$, which distinguishes the estimation procedure from the previous two approaches. On the one hand, even if $W$ and $S^{-1}$ satisfy the boundedness conditions of the $L_1$ and $L_\infty$ norms, it cannot be directly proven that $\Omega^{-1}$ is bounded in $L_1$ and $L_\infty$ norms. Consequently, the first method cannot be employed. On the other hand, notice that $S^{\top}S$ in $\Omega^{-1}$ cannot be analyzed as a transition-probability matrix; thus, the assumption of the existence of a stationary distribution cannot be made either. This makes the second method inapplicable. Therefore, we need to reexamine the theoretical properties of the CLE $\hat\theta_{\text{CLE}}$ based on the PSAR model. We begin with the discussion of the identification issue. Before presenting the theoretical results, the following identification condition is required.

(C4) (Identification) Assume that at least one of the following conditions holds: (a) $\lim_{N \to \infty} N^{-1} \beta_0^\top \mathbb{X}^\top (WS_0^{-1})^\top \Omega^{-1} H_{\Omega} WS_0^{-1} \mathbb{X} \beta_0$ is positive for any $(\rho,\sigma^2)$ in parameter space; or (b) $\lim_{N \to \infty} N^{-1} \beta_0^\top \mathbb{X}^\top (WS_0^{-1})^\top \Omega^{-1}H_{\Omega} WS_0^{-1} \mathbb{X} \beta_0 = 0$, sequence $\{h_n\}$ is bounded, and for any $(\rho,\sigma^2) \neq (\rho_0,\sigma_0^2)$, ${ \lim_{N\rightarrow \infty}}N^{-1} \mbox{tr} \big\{S^\top \Omega^{-1} S S_0^{-1} \Omega_0 (S_0^{-1})^\top \big\} \neq 1. $

Condition (C4) assumes that either (a) covariates $\mathbb{X}$ and $WS_0^{-1}\mathbb{X}\beta_0^\top$ do not exhibit asymptotic multicollinearity, or (b) the variance matrix of $\mathbb{Y}^*$ is unique. A similar type of condition was also assumed by Lee:2004. Then, we have the following theorem.

The proof of Theorem (ref) is given in Appendix C.1. We then examine the statistical property of the CLE, which results in the following theorem.

The proof of Theorem (ref) is given in Appendix C.2. The theorem guarantees the asymptotic normality of $\hat \theta_{\text{CLE}}$. By combining the previous discussion on identification issues, we can employ the proposed Newton-Raphson-type method to iteratively obtain the CLE by Algorithm (ref). In real practice, we obtain $\hat\theta_{\text{CLE}}^{(t)}$ instead of $\hat \theta_{\text{CLE}}$. Thus, it is necessary to discuss the relationship between $\hat \theta_{\text{CLE}}$ and $\hat\theta_{\text{CLE}}^{(t)}$ in the proposed algorithm. In this regard, the following theorem could be established.

The proof of Theorem (ref) is given in Appendix C.3. Theorem (ref) guarantees the numerical convergence of the proposed algorithm. This suggests that the proposed iterative algorithm can obtain $\hat \theta_{\rm CLE}$ with asymptotic normality.

Remarkably, for the inference of $\hat \theta_{\text{CLE}}$, the asymptotic covariance in (ref) cannot be immediately obtained when $\mathbb{X}$ is unobserved. Therefore, we rewrite $\Sigma_k^{\text{CL}}=\Sigma_k^{\text{CL}}(\mathbb{X},\theta)$ $(k=1,2)$ as a function of $\mathbb{X}$ and $\theta$. In this subsection, for simplicity, we use $\theta$ to represent the true parameter. Then, the difference between $\Sigma_k^{\text{CL}}(\mathbb{X}^*,\theta)$ and $\Sigma_k^{\text{CL}}(\mathbb{X},\theta)$ should be considered. This is because data users can only observe $\mathbb{X}^*$ instead of $\mathbb{X}$, and statistical inference can be made relying only on $\mathbb{X}^*$. Thus, define $\Delta_k(\theta)=E\{\Sigma_k^{\text{CL}}(\mathbb{X}^*,\theta)\}-\Sigma_k^{\text{CL}}(\mathbb{X},\theta)$. See Appendix A.1 for detailed expression of $\Delta_k(\theta)$. Therefore, based on the idea of bias correction, we can obtain a bias-corrected estimator for each element of the covariance matrix in (ref) based on the observed $\mathbb{X}^*$. The estimators can be denoted as $\widehat{\Sigma}_k^{\text{CL}}=\Sigma_k^{\text{CL}}(\mathbb{X}^*,\hat \theta_{\text{CLE}})-\Delta_k^{\text{CL}}(\hat \theta_{\text{CLE}})$ $(k=1,2)$ with $\hat \theta_{\text{CLE}}$ plugged in. Employing a technique similar to that used in the proof of Lemma 6 in Appendix B, the consistency of $\widehat{\Sigma}_k^{\text{CL}}$ can be established. Then, the estimated standard error of $\hat \theta_{\text{CLE}}$ could be calculated by the estimated asymptotic covariance matrix $(\widehat{\Sigma}_2^{\text{CL}})^{-1}\widehat{\Sigma}_1^{\text{CL}}(\widehat{\Sigma}_2^{\text{CL}})^{-1}$.

Then, we consider gaining a more intuitive understanding of the privacy-protection level and statistical efficiency of $\hat \theta_{\text{CLE}}$. Notably, the covariance structure of the estimator is quite complex here. For simplicity, we consider a special case of the pure SAR with $\beta=0$, known $\sigma^2$, and Gaussian error terms. We focus on estimating the network effect $\rho$ to intuitively express the impact of adding noise. Considering that $\rho$ is usually sufficiently small chen2013impact, we are then motivated to conduct Taylor's expansion on $\Sigma_k^{\text{CL}}$s ($k=1,2$) to obtain their leading terms in order to approximate the asymptotic covariance of the CLE. Then, it could be calculated that $\Delta_{\rho\rho}^{\rm CL}=o(1)$ and $\Sigma_2^{\rm CL}(\Sigma_1^{\rm CL})^{-1}\Sigma_2^{\rm CL}=$

Thus, we can clearly observe that as the privacy-protection level $\lambda^2$ increases, the asymptotic variance $(\Sigma_2^{\text{CL}})^{-1} \Sigma_1^{\text{CL}} (\Sigma_2^{\text{CL}})^{-1}$ gradually increases. Increasing the privacy-protection level $\lambda^2$ enhances data security, but inevitably decreases the statistical efficiency of the resulting estimator. Therefore, in practice, data providers need to set a reasonable level of privacy protection to balance between data security and statistical efficiency. Efforts should be made to ensure that while protecting data, the statistical efficiency of the estimates obtained by data users remains acceptable. Note that for more general cases, providing the specific form of the asymptotic variance is difficult as the privacy-protection level $\lambda^2$ varies. We show more evidence of this through numerical simulations.

In this subsection, we establish the theoretical properties of the CLS. The analysis procedure is similar to those in Section 3.2. Therefore, for brevity and to avoid repetition, we will omit the detailed discussion of the technique here. We start with the identification issue. To this end, define $\mathbb{X}_s=WS_0^{-1}\mathbb{X}\beta_0$ and $\widetilde \mathbb{X}=(\mathbb{X}_s,\mathbb{X})\in\mathbb{R}^{N\times(p+1)}$. The following assumption is necessary for the identification of the parameter.

(C4$^*$) (Identification) Assume that (a) $\lim_{N\to\infty}N^{-1} \widetilde \mathbb{X}^\top \widetilde \mathbb{X}$ exists and is non-singular, and $\delta>0$ exists such that $\min_{|\rho|\leq 1-\delta}\lambda_{\min}(S S^\top)\geq \tau$, where $\tau$ is a positive constant; or (b) for the SAR model with no exogenous covariates, assume $I_N$, $W$, $W^\top$, and $W^\top W$ are linearly independent.

Here, the identification issue is discussed in two cases: with covariates and without covariates. We then have the following theorem.

The proof of Theorem (ref) is given in Appendix C.4. Next, we establish the following numerical and statistical theoretical properties of $\hat \gamma_{\text{CLS}}$.

Through Theorem (ref), we can conclude that the CLS estimator $\hat \gamma_{\text{CLS}}$ is $\sqrt{N}$-consistent, which has the same convergence rate as the CLE. For statistical inference, similar to the discussion for $\hat \theta_{\text{CLE}}$, we can also provide consistent estimators for each element in the asymptotic covariance matrix (ref) based on the observed $\mathbb{X}^*$. Rewrite $\Sigma_k^{\text{LS}}=\Sigma_k^{\text{LS}}(\mathbb{X},\gamma)$ $(k=1,2)$ as functions of $\mathbb{X}$ and $\gamma$, which are expressed in (ref). Then, we could plug $\mathbb{X}^*,\hat \gamma_{\text{CLS}}$ in to obtain $\Sigma_k^{\text{LS}}(\mathbb{X}^*,\hat \gamma_{\text{CLS}})$ for $k=1,2$. As a result, $\Delta_k^{\text{LS}}(\gamma)=E\{\Sigma_k^{\text{LS}}(\mathbb{X}^*,\theta)\}-\Sigma_k^{\text{LS}}(\mathbb{X},\theta)$ could be calculated. See Appendix A.2 for a detailed expression of $\Delta_k^{\text{LS}}(\gamma)$. With $\hat \gamma_{\text{CLS}}$ plugged in, the consistent estimators for $\Sigma_k^{\text{LS}}$ could be obtained as $\widehat{\Sigma}_k^{\text{LS}}=\Sigma_k^{\text{LS}}(\mathbb{X}^*,\hat \gamma_{\text{CLS}})-\Delta_k^{\text{LS}}(\hat \gamma_{\text{CLS}})$ $(k=1,2)$. In this way, we could obtain the estimated standard error of $\hat \gamma_{\text{CLS}}$ by the estimated asymptotic covariance matrix $(\widehat{\Sigma}_2^{\text{LS}})^{-1}\widehat{\Sigma}_1^{\text{LS}}(\widehat{\Sigma}_2^{\text{LS}})^{-1}$. Using a technique similar to that used in the proof of Lemma 8 in Appendix B, the consistency of $\widehat{\Sigma}_k^{\text{LS}}$ can be established.

The previously discussed model (ref) represents the simplest form of the SAR model. In this subsection, we study several extensions. Specifically, (1) we consider multiple dependent variables \(Y\) instead of a single one, extending the SAR model to the multivariate spatial autoregressive model zhu2020multivariate; (2) we consider networks with different types of nodes, extending the network to a multi-mode network regression model that accommodates varying node types huang2020two; (3) we consider perturbations in the network structure as a preliminary attempt, exploring how they impact the robustness and reliability of the privacy-preserving techniques Lewbel2024.

Privacy-Protected Multivariate SAR. First, we consider extending the PSAR model to the privacy-protected MSAR model, which can be adapted to model the mutual influence of different dependent variables through the network structure. Suppose there are multivariate responses $\mathbb{Y} = (Y_{ij})\in \mathbb{R}^{N\times q}$ and error matrix $\tilde{{\mathbb E}}=(e_{ij}) \in \mathbb{R}^{N\times q}$. Denote the network parameters as $D = (d_{j'j}) \in \mathbb{R}^{q\times q}$ and the covariates parameters as $B = (b_{kj}) \in \mathbb{R}^{p\times q}$. The multivariate spatial autoregressive model is defined as $\mathbb{Y} = W\mathbb{Y} D + \mathbb{X} B + \tilde{\mathbb E}$. Here, $d_{j_1,j_2}$ ($1\leq j_1\neq j_2\leq q$) represents the {\it extra-activity effect}, measuring network effect between different responses, and $d_{jj}$ ($1\leq j\leq q$) represents the {\it inter-activity effect}, measuring network effect within the same response. Further, we define $\mathcal{Y} = {\rm vec}(\mathbb{Y}) \in\mathbb{R}^{Nq},\mathcal{X} = I_q \otimes\mathbb{X}\in\mathbb{R}^{Nq\times pq}$. Then the vector norm of the privacy-protected MSAR is,

where the observed data is given by $\mathcal{Y}^*= \mathcal{Y} + \mathcal{E}$, $\mathcal{X}^*=(\mathcal{X}_1,\mathcal{X}_2^*)=(\mathcal{X}_1,\mathcal{X}_2+\mathcal{E}_x)$, $\mathcal{E}$ and $\mathcal{E}_x$ are artificially added independent noises with zero mean and variance $\lambda^2$ and $\lambda_x^2$ correspondingly. In this way, define \(\gamma_M \in \mathbb{R}^{q^2 + pq}\) as the parameter vector that incorporates all the information from \(D\) and \(B\). We construct the objective function and obtain the CLS estimator $\hat{\gamma}_{M,\rm CLS}$ for $\gamma_M$ by correcting for the first- and second-order derivatives to deal with observational errors in $\mathcal{X}^*$ and $\mathcal{Y}^*$. Consequently, we can establish the following corollary in a manner similar to Theorem (ref). Detailed notations are provided in Appendix D, and the proof is omitted.

Privacy-Protected Multi-Mode SAR. Next, we employ the privacy-protected multi-mode SAR model to measure the network influence among different types of nodes. Consider a multi-mode network with $n_k$ nodes in the $k$th group ($1\leq k\leq K$), and $N=\sum_{k}n_k$. Define $a_{i_1i_2}=0$ if $i_1$ and $i_2$ are in the same group. Let $Y_k=(Y_{k,1}, \cdots,Y_{k,n_k})^\top\in\mathbb{R}^{n_k}$ $(1\leq k\leq K)$ be the response collected for the $k$th group of nodes. And define $X_k=(X_{k,1}^\top, \cdots,X_{k,n_k}^\top)^\top\in\mathbb{R}^{ n_k\times p_k}$ as the exogenous covariates for the $k$th group. Furthermore, define $\varepsilon_k=(\varepsilon_{k,1}, \cdots,\varepsilon_{k,n_k})^\top\in\mathbb{R}^{n_k}$ as the noise vector. The multi-mode autoregressive model can be defined as

where \(\rho_{k_1k_2}\) (\(1\leq k_1,k_2\leq K, k_1\neq k_2\)) are the cross-mode effects, \(\beta_{k_1}\) are the parameters corresponding to \(X_{k_1}\). Similarly, define \(\gamma_T\in\mathbb{R}^{K(K-1)+\sum_k p_k}\) as the parameter vector that encompasses all the information from \(\rho_{k_1k_2}\) and \(\beta_{k_1}\). We can construct the least-square type objective function and obtain the CLS estimator $\hat{\gamma}_{T,\rm CLS}$ for $\gamma_T$ based on corrected derivatives to deal with observational errors in $\mathbb{Y}^*$ and $\mathbb{X}^*$. Using this approach, the following corollary can be established for a two-mode network as an example. Specifically, we observe $\mathbb{Y}^*,\mathbb{X}^*$ with $\mathbb{Y}^*=(Y_1^\top,Y_2^\top)^\top+\mathcal{E}$, $\mathbb{X}^*=[\mathbb{X}^*_1,0_{n_1\times p_2};0_{n_2\times p_1},\mathbb{X}^*_2]^\top$, $\mathbb{X}^*_1=(\mathbb{X}_{11},\mathbb{X}^*_{12})=(\mathbb{X}_{11},\mathbb{X}_{12}+\mathcal{E}_{x_1})$, $\mathbb{X}^*_2=(\mathbb{X}_{21},\mathbb{X}^*_{22})=(\mathbb{X}_{21},\mathbb{X}_{22}+\mathcal{E}_{x_2})$, and $\mathcal{E}\in\mathbb{R}^{n_1+n_2}$, $\mathcal{E}_{x_1}\in\mathbb{R}^{n_1\times p_{12}}$, $\mathcal{E}_{x_2}\in\mathbb{R}^{n_2\times p_{22}}$ are artificially added noises with variance $\lambda^2, \lambda_x^2$ respectively, $p_{12}$ and $p_{22}$ represent the dimensions of the covariates requiring privacy protection for the two types of nodes. Detailed expressions are provided in Appendix D and the proof is omitted.

Privacy-Protected SAR with Perturbed Network. As a further discussion of the PSAR model, we consider randomly perturbing the edges in the network to protect the privacy of the network structure. The true network is represented by \(A\), and we define the noise-perturbed adjacency matrix as \(A^*\). In this way, we obtain the observed weighting matrix \(W^*\), which is normalized from \(A^*\). Consequently, we calculate the corrected least squares estimator \(\hat\gamma^*_{\rm CLS}\) based on \(W^*\) instead of \(\hat\gamma_{\rm CLS}\) based on \(W\). However, we cannot allow a large number of edges in the network to change, as excessive noise would prevent us from obtaining consistent parameter estimates. To address this issue, and inspired by Lewbel2024, we propose the following constraint condition (C5). Based on (C5), we establish the theoretical properties of the CLS under the condition of perturbed network edges. This corollary indicates that once condition (C5) is satisfied, we can still ensure the consistency and asymptotic normality of the CLS estimator.

(C5) (Perturbed Network) Assume \(\sum_i\sum_jE(|A^*_{ij}-A_{ij}|)=O(N^s)\) for some positive constant \(s < 1/2\), and assume that \(W^*\) is uniformly bounded in both row and column sums in probability.

\csection{NUMERICAL STUDIES}

To demonstrate the finite-sample performance of the proposed methods, we present three simulation examples based on the following generating mechanisms: the network structure $A$, data generation, noise distribution, and privacy-protection level. The classical QMLE (treating the observed variables as the true ones) and the proposed estimators (i.e., CLE and CLS) are compared. The network structures are given as follows.

Example 1. (Dyad Independence Network) Following holland1981exponential, we define a dyad as $\mathbf{A}_{ij} = (a_{ij}, a_{ji})$ ($1\le i < j\le N$) and assume that different $\mathbf{A}_{ij}$s are independent. To allow for network sparsity, we set $P(\mathbf{A}_{ij} = (1,1)) = 10N^{-1}$ and $P(\mathbf{A}_{ij} = (1,0)) = P(\mathbf{A}_{ij} = (0,1)) = 0.5N^{-0.8}$. Then, the probability of a null dyad is $P(\mathbf{A}_{ij} = (0,0)) = 1-10N^{-1}-N^{-0.8}$, which is close to 1 when $N$ is large.

Example 2. (Stochastic Block Network) The next network type that is considered is the stochastic block network wang1987stochastic,nowicki2001estimation. Following nowicki2001estimation, we randomly assign a block label $k$ ($1\leq k\leq K$) for each node with $K= 20$ as the total number of blocks. Define $P(a_{ij}=1) = 20N^{-1}$ if $i$ and $j$ belong to the same block, and $P(a_{ij}=1) = 2N^{-1}$ otherwise. Thus, the nodes in the same block are more likely to be connected.

Example 3. (Power-Law Distribution Network) It is commonly observed in network analysis that the majority of nodes have few links but a small proportion have a large number of edges barabasi1999emergence. Therefore, we simulate the adjacency matrix $A$ according to clauset2009power. The in-degree $m_i = \sum_j a_{ji}$ for node $i$ is generated by the discrete power-law distribution with $P(m_i = k) = ck^{-\alpha}$, with a normalizing constant $c$ and $\alpha = 3$. For the $i$th node, $m_i$ nodes are randomly selected to be its followers.

Data Generation. For each node, we generate the true exogenous covariates $\mathbf{x}_i = (x_{i1}, x_{i2})^\top \in \mathbb{R}^2$ from a multivariate normal distribution with mean $\mathbf{0}_2$ and $\Sigma_x= I_{2}$. The corresponding network autoregression coefficient is fixed to be $\rho_0 = 0.2$ and coefficient $\beta_0=(0.3,0.3)^\top$. The response $\mathbb{Y}^*$ is generated based on $\mathbb{Y}^*= (I_N-\rho_0 W)^{-1}(\mathbb{X}\beta_0+{\mathbb E})+{\mathcal E}$. Define $X_j=(x_{1j},\cdots,x_{nj})^\top$ for $1\leq j\leq 2$. To simulate privacy-protected covariates, we fix $\mathbb{X}_1=X_1\in \mathbb{R}^{N}$ and assume that $\mathbb{X}_2=X_2\in \mathbb{R}^{N}$ cannot be faithfully observed. Then, $\mathbb{X}_2^*$ can be generated by $\mathbb{X}_2^*=\mathbb{X}_2+{\mathcal E}$ with ${\mathcal E}=(\varepsilon_{x,i})$, which follows the setting of the noise distribution.

Noise Distribution. We consider that $\varepsilon_{x,i}$ follows a normal distribution with mean $0$ and $\lambda_x^2=0.5$. We consider two different distributions to generate $\varepsilon_i$ independently with mean $0$ and $\lambda^2=0.5$: (1) a normal distribution, and (2) a $t$-distribution with degree 6. It is notable that for the $t$-distribution, each element is divided by $\sqrt{3}$ to make $\mbox{var}(\varepsilon_i)=\lambda^2=0.5$. For the noise $e_i$, we consider the same cases for $N(0,1)$ and $t(6)$. Here, we set the sample size as $N=$500, 1,000, and 2,000 and consider all the network models.

Privacy-Protection Level. To better illustrate the effect of the privacy-protection level, we consider a fixed sample size $N=1,000$ and the dyad independent network-generation model as an example. Different $\varepsilon_i$s and $\varepsilon_{x,i}$s are generated from normal distributions with variances of $\lambda^2$ and $\lambda^2_x$, respectively. We consider the impact of the variations in $\lambda_x^2$ and $\lambda^2$ on the estimation results. Specifically, for fixed $\lambda^2=0.5$, we consider $\lambda^2_x=(0.2,0.5,0.8)$. For fixed $\lambda^2_x=0.5$, we consider $\lambda^2=(0.2,0.5,0.8).$

For a better comparison, we focus here on the estimation results of $\rho_0$ and $\beta_0$. To gauge the finite-sample performance, we use the following metrics. Define $\hat D^{(r)} = \{\hat d_{j}^{(r)}\}_{j=1}^q \in \mathbb{R}^q$ as the estimator from the $r$th replication with $q=p+1$. For any $1 \leq j \leq q$, the bias can be evaluated as ${\rm Bias}_j=|\bar{d_{j}}-d_j|$, where $\bar{d_{j}}=R^{-1}\sum_{r}\hat d_{j}^{(r)}$, and $d_j$ is the $j$th element of the true parameter. The standard error can be estimated using $\widehat{\mbox{SE}}_j=R^{-1}\sum_{r} \widehat{\mbox{SE}}_j^{(r)}$. Notably, $\widehat{\mbox{SE}}_{j}^{(r)}$ represents the $j$th diagonal element of the estimated asymptotic covariance matrix, which can be computed using ((ref)) or ((ref)) with $\hat D^{(r)}$ plugged in, followed by a correction step, as discussed in Section 3.3 for the CLE and Section 3.4 for the CLS. Define ${\rm SE}_j=\{R^{-1}\sum_{r}(\hat d_{j}^{(r)}-\bar{d_{j}})^2\}^{1/2}$ as the Monte Carlo standard deviation of $\hat d_{j}^{(r)}$, and the estimation efficiency of $\widehat {\rm SE}_j$ can be evaluated by comparing $\widehat{\mbox{SE}}_j$ and ${\rm SE}_j$. In addition, for each $\hat d_{j}^{(r)}$, a 95% confidence interval can be constructed as $\mbox{CI}_{j}^{(r)} = [\hat d_{j}^{(r)}-z_{0.975}N^{-1}\widehat{\mbox{SE}}_{j}^{(r)},\hat d_{j}^{(r)}+z_{0.975}N^{-1}\widehat{\mbox{SE}}_{j}^{(r)}]$, where $z_\alpha$ is the lower $\alpha$th quantile of the standard normal distribution. The empirical coverage probability is then evaluated as $\mbox{CP}_{j} = R^{-1}\sum_{r = 1}^R I(\hat d_{j}^{(r)}\in\mbox{CI}_{j}^{(r)} )$, where $I(\cdot)$ is the indicator function.

Each experiment is replicated 500 times ($R = 500$). All simulations are conducted on a Linux server with a 3.60 GHz Intel Core i7-9700K CPU and 16 GB RAM. Due to the similar estimation performances across different noise distributions, we only present the results for the normal distributions of both ${\mathbb E}$ and ${\mathcal E}$ in Tables (ref), and the other results are provided in Appendix E. Moreover, the estimation outcomes for different values of $\lambda^2$ and $\lambda_x^2$ are illustrated in Table (ref). We present the Bias, SE, and CP in Tables (ref)-(ref). Then, we show the averaged CPU time for deriving CLE and CLS in Figure (ref) to evaluate the computation efficiency.

Estimation Performance. From Table (ref), we can see that the classical QMLE is seriously biased as expected. Thus we focus on the performance of the proposed estimator CLE and CLS. We draw the following conclusions. First, the estimation bias is sufficiently small for all sample sizes and both the methods. Second, as the sample size $N$ increases, Bias, SE, and $\widehat{\mbox{SE}}$ all decrease, which shows the consistency of both methods. Third, all the CP values are approximately 95% for $\alpha=0.05$ with $\widehat{\mbox{SE}}$ and SE clearly being close to each other. This corroborates the theoretical conclusions in Theorem (ref) for the CLE and in Theorem (ref) for the CLS. Finally, CLE has a smaller SE than CLS.

Effect of Privacy-Protection Level. From Table (ref), we can see the bias of $\rho$ or $\beta$ of the classical QMLE increases as the noise variance $\lambda^2$ or $\lambda_x^2$ increases, respectively. Based on the results of CLE and CLS in Table (ref), we can conclude that an increase in the privacy-protection level leads to a larger $\widehat{\rm SE}$ for the estimators. Specifically, when $\lambda_x^2$ is fixed, increasing $\lambda^2$ leads to an increase in $\widehat{\rm SE}$ for all estimators. However, when $\lambda^2$ is fixed, as $\lambda_x^2$ increases, only $\widehat{\rm SE}$ for the estimator corresponding to the unobserved $\mathbb{X}_2$ increases, whereas those for the other estimators remain relatively stable.

Computation Time. We conduct the same experiment as previously, but with a normal noise distribution and fixed $\sigma_0^2=1,$ $\lambda_x^2=1,$ and $\lambda^2=0.5$. Moreover, we change $N$ from 500 to 5,000 for $R = 100$ (100 replicates) to further compare the computational efficiency of CLE and CLS. Figure (ref) shows the results from a computational perspective. Evidently, as the sample size $N$ increases, the computational time of CLE increases dramatically, whereas that of CLS increases much more slowly. In practice, a method should be chosen based on the balance between computational time and statistical accuracy.

We compare the CLE and CLS using a transaction dataset of small- and medium-sized restaurants. The dataset is provided by Shouqianba, which is a leading company in China's mobile payment market (http://www.shouqianbao.com.cn/). This dataset contains information from $N=2024$ fast-food restaurants specializing in snacks in Guangzhou, China. For each restaurant $i$ $(1\leq i\leq N)$, $Y_i$ is defined as the transaction volume of the restaurant from April 1st to June 1st, 2024. Subsequently, for each restaurant, we consider three covariates: (1) repeat customers $X_1$, which is defined as the percentage of consumers who dined at the restaurant for two or more times in the previous two months; (2) operating hours $X_2$, wherein $x_{2i}=1$ indicates that the restaurant operates during both lunch and dinner hours, and $x_{2i}=0$ indicates otherwise; (3) transaction amount per customer $X_3$, which is defined as the average amount of each transaction in the previous two months. All the continuous variables are standardized. For privacy protection, the platform has added Gaussian noise to the average transaction amount per customer for restaurants with $\lambda_x^2=0.2$, and added Gaussian noise to the response $\mathbb{Y}$ with $\lambda^2=0.25$. To analyze the effect of the network structure, we construct a network between restaurants based on the customers who have dined at them. Specifically, if two restaurants $i$ and $j$ have at least five shared customers, an edge is defined with $a_{ij}=1$; otherwise, $a_{ij}=0$. Furthermore, the network density (i.e., $\{N(N-1)\}^{-1}\sum_{i_1,i_2}a_{i_1i_2}$) is 0.31%, which implies that this is a relatively sparse network.

Then, we estimate the model using the CLE and CLS estimators. In addition to reporting the measurements from the simulation results, we also provide the root mean square error (RMSE) of each method for comparison. To be more precise, define $(\hat\rho,\hat\beta)$ as the estimators of $(\rho_0,\beta_0)$ (i.e., CLE or CLS). The RMSE is calculated by $ \mbox{RMSE} =\{N^{-1}\sum_{i=1}^N(\hat{Y_i}-Y_i)^2\}^{1/2}$, and $ (\hat{Y_1},\dots,\hat{Y}_N)^T=(I_N-\hat\rho W)^{-1}\mathbb{X}^{*}\hat\beta$. The estimation results are summarized in Table (ref).

From the results, we can draw the following conclusions. First, it can be observed that the estimation results for both the network effect and covariates using both methods are similar, where the estimated values differ by no larger than 0.04. Second, the network effects estimated by both methods are significantly negative. This is because the network construction is based on customers who have dined at the same restaurants. The more customers who have dined at both restaurants, the stronger the competition between those restaurants. Thus the negative network effect can be intuitively explained. Third, the CLS has slightly larger RMSE than that of the CLE. Lastly, the CLS only requires less than one tenth of the computational time of the CLE as it does not involve the calculation of $\Omega^{-1}$.

\csection{CONCLUDING REMARKS}

This study addressed the issue of privacy protection in the SAR model with noise-added response and covariates. We focused on achieving a consistent estimation of the model parameters and proposed two different estimators: the CLE and the CLS. Both of the estimators provide reliable solutions in privacy-protected scenarios. Theoretical properties have been carefully established and model generalizations are discussed.

For future research, we will discuss three potential directions. First, we focused on the classical form of the SAR model and some of its extensions. While these models are well-established, exploring the properties of estimators in dynamic and nonlinear network models with added noise remains an important area for future research. Second, Condition (C5) imposes relatively strict requirements. Investigating how to better protect the privacy of network edges by relaxing these conditions is essential and warrants further study. Finally, examining other methods of privacy protection for the spatial autoregressive model and its extensions presents an intriguing and valuable research direction.