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.
74,786 characters · 12 sections · 31 citation commands
Uniform Inference on High-dimensional Spatial Panel Networks
Network analysis has gained significant interest in recent years. In particular, measuring connectedness within a complex system has become a central task in learning networks. Various forms of regression, where the dependent variables are affected by the outcomes and characteristics of network members, have been formulated for that purpose. The established literature on social network analysis favors using a predetermined network structure, which is fully characterized by a specified adjacency matrix, to study peer effects in social networks; see, for example, lee2007identification,bramoulle2009identification,lee2010specification,yang2017identification,zhu2020multivariate. As for spatial panel networks, kuersteiner2020dynamic consider a class of GMM estimators for general dynamic panel models that allow for potential endogeneity and cross-sectional dependence. An alternative to imposing a known network structure is to estimate the adjacency matrix, provided that the structural parameters are already identified. Examples of related studies include blume2015linear,de2018recovering,lewbel2019social.
With the rise of big data availability, many applications are concerned with large-scale networks consisting of a large number of individuals. In particular, spatial panel data involving high-dimensional time series are observed in many financial and economic network analyses. This poses the challenge of estimating too many unknown parameters. To reduce the dimensionality, various machine learning methods based on sparsity and penalization are employed to shrink the parameters. manresa2013estimating uses LASSO (Least Absolute Shrinkage and Selection Operator) to quantify the spillover effects in social networks, where the endogenous interactions are not taken into consideration. de2018recovering apply Adaptive Elastic Net GMM to estimate the interaction model with important contributions to the identification of the structural parameters. ata2018latent consider a reduced-form estimation with the innovative discovery of the algebraic results on how the sparsity of the structural parameters relates to that of the parameters in the reduced form. lam2014regularization study the penalized estimation of the spatial weight matrix in a spatial lag model through adaptive LASSO and show the oracle properties of the sparse estimator. wang2024panel develop a high-dimensional interactive fixed effects estimator that allows for a growing number of latent factors and apply it to peer-effects analysis in networks with sparse links. They demonstrate the consistency of the new estimator and the asymptotic normality of the post-selection estimator of the slope coefficients. In this paper, we also aim to conduct inference on the network structure.
Machine learning methods are notably effective in improving prediction performance. However, statistical inference may suffer from substantial bias due to omitted variables. Debiasing is necessary to construct high-quality point and interval estimates. Taking LASSO-type methodologies as example, lam2014regularization establish the asymptotic normality of non-zero elements in the network structure. However, in practice, we often lack prior information about whether parameters are truly non-zero, necessitating a uniform inference theory that allows testing any parameters of interest. For independent and identically distributed (i.i.d.) data, extensive research explores uniform inference in high-dimensional regression settings under exogeneity conditions (e.g., BCH2014,zhang2014debiased,BCK15Bio,DML) and, more generally, considers GMM frameworks that allow for endogeneity (e.g., belloni2018high,belloni2017simultaneous,caner2018high), through various de-biasing and orthogonalization techniques. Building on the idea of orthogonality, ata2018latent present an algorithm incorporating bias-corrected Dantzig selector estimator to investigate large networks with latent agents, though without accounting for temporal dependence. Addressing data-generating processes exhibiting dependency, lasso2018 study LASSO-based inference for exogenous regression under general temporal and cross-sectional dependence.
In this paper, we are motivated by the need to understand the connectedness within a complex spatial panel network. Our focus is on exploring network structures, which need not be sparse, while allowing for flexible sparse deviations. These deviations can be viewed as either latent or misspecified relative to a predetermined adjacency matrix (e.g., credit chains or common ownership information in a financial system). Specifically, we examine network formation by framing the problem as a general system of dynamic regression equations, considering both temporal and spatial dependencies inherent in the data-generating processes. Methodologically, we extend the model setting in lasso2018 by allowing for endogeneity in the covariates, which is a natural concern when the regression system is featured with simultaneity by incorporating contemporaneous lags. As a result, sufficiently many moment conditions involving instrumental variables (IV) are needed and we build a debiased-regularized, high-dimensional GMM estimator to facilitate valid inference. Notably, the double LASSO estimation steps used in lasso2018 for debiasing are unsuitable in our case due to the endogeneity issue. This necessitates the identification of an appropriate moment selection matrix to achieve the desired orthogonality for valid inference. Given the high-dimensional nature of the covariance matrix and its inverse, a unified regularized estimation framework is required to ensure the consistency of both the preliminary estimator and the matrices involved in the debiasing step.
For implementation, we propose employing a Generalized Dantzig Selector (GDS) as an initial step, followed by a debiasing step. Theoretically, we establish the consistency of the GDS estimator and derive the linearization of the debiased estimator to enable the application of the central limit theorem for uniform inference on the parameters of interest (whether of fixed or growing dimension). In particular, we show the asymptotic properties of the debiased-regularized GMM (DRGMM) estimator for both linear and nonlinear moments cases. Moreover, we discuss the connection to the semiparametric efficiency literature, particularly in relation to the construction of our estimator when the dimension of the parameters of interest is fixed.
We contribute to the literature in four respects. First, we develop a method for estimating parameters in a high-dimensional endogenous equation system that incorporates both spatial and temporal dynamics. Our theoretical framework accords with general dynamic panel models, capturing heterogeneity through individual-specific parameters. Second, we propose a latent model that shrinks toward a pre-specified network structure. In particular, we provide theoretical insights into how the restricted eigenvalue conditions on the design matrix adapt to the transformation of the covariates. Third, we employ a debiased machine learning approach to conduct simultaneous hypothesis testing on high-dimensional parameters. Finally, we demonstrate the practical utility of our method through an empirical application in a financial network context.
Compared to the high-dimensional GMM estimator developed in belloni2018high, this study involves a spatial panel model setup, rather than i.i.d. data, introducing several technical challenges. First, to prove consistency, the verification of certain high-level assumptions requires significantly different steps. We demonstrate the validity of concentration under spatial-temporal dependent processes, ensuring that panel data with a network structure can be properly handled. Furthermore, to extend the framework to nonlinear and even non-smooth moments, we employ different techniques for proving tail probabilities and concentration inequalities, as detailed in Appendix (ref).
The following notations are adopted throughout the paper. For a vector $v = (v_1, \ldots, v_p)^\top$, let $|v |_k = (\sum_{i=1}^p |v_i|^k)^{1/k}$ with $k\geq1$, $|v|_\infty = \max\limits_{1\leq i\leq p} |v_i|$, and $|v|_0$ denote the number of nonzero components of the vector. For a random variable $X$, let $\|X\|_r\stackrel{\mathrm{def}}{=}(\mathop{\mbox{\sf E}}|X|^r)^{1/r}$, with $r>0$. For a matrix $A = (a_{ij})\in\mathbb{R}^{p\times q}$, we define $|A|_1 = \max\limits_{1\leq j\leq q} \sum_{i=1}^p |a_{ij}|$, $|A|_{\infty} = \max\limits_{1\leq i\leq p} \sum_{j=1}^q|a_{ij}|$, $|A|_{\max} = \max\limits_{1\leq i \leq p,1\leq j\leq q}|a_{ij}|$, and the spectral norm $|A|_2 = \sup_{|v|_2\leq1} |Av|_2$. Moreover, let $\lambda_i(A)$ denote the $i$-th largest eigenvalue of a square matrix $A$, and let $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the minimal and maximal eigenvalues of $A$, respectively. Similarly, let $\sigma_i(A)$ denote the $i$-th largest singular value of $A$, with $\sigma_{\min}(A)$ and $\sigma_{\max}(A)$ representing the minimal and maximal singular values of $A$, respectively. Let $\mathbf I_{p}$ denote the identity matrix of size $p\times p$. For any measurable function on a measurable space $g:\mathcal{W}\rightarrow\mathbb{R}$, define the sample average over the indices $t=1,\ldots,n$ as $\mathop{\mbox{\sf E}}_n(g(\omega_t))\stackrel{\mathrm{def}}{=} n^{-1}\sum_{t=1}^n g(\omega_t)$. Given two sequences of positive numbers $a_n$ and $b_n$, write $a_n\lesssim b_n$ (resp. $a_n\asymp b_n$) if there exists constant $C>0$ (independent of $n$) such that $a_n/b_n\leq C$ (resp. $1/C \le a_n / b_n\leq C$) for all large $n$. For a sequence of random variables $x_n$, we use the notation $x_n\lesssim_{\mathrm{P}} b_n$ to denote $x_n=\mathcal{O}_{\mathrm{P}}(b_n)$.
The rest of the article is organized as follows: Section (ref) outlines the model specification and estimation steps. Section (ref) presents the main theoretical results for the case of linear moments. Sections (ref) and (ref) provide simulation studies and an empirical application on financial network analysis with potential misspecification. The technical proofs and additional details--including extension to nonlinear moments, connection to semiparametric efficiency, and supplementary discussions--are provided in the Online Appendix. The codes to implement the algorithms are publicly accessible via the GitHub repository: \href{https://github.com/huangche/Uniform-Inference-on-High-dimensional-Spatial-Panel-Networks}{Uniform-Inference-on-High-dimensional-Spatial-Panel-Networks}.
For time points $t=1,\ldots,n$ and individual entities $j=1,\ldots,p$ (both $n,p$ tend to infinity), we consider a spatial panel network model for the nodal response $y_{j,t}$:
where we have an observed network structure $w_j=(w_{j1},\ldots,w_{jp})^\top$ for all $j=1,\ldots,p$, and $\rho^0$ is the spatial autoregressive parameter. In particular, $w_j^\top y_t$ is an observed weighted variable, and vectors $\delta^0_j=(\delta^0_{j1},\ldots,\delta^0_{jp})$, $j=1,\ldots,p$, denote approximately sparse misspecification errors of the network structure. Estimation and inference of $\delta_j^0$ and $\rho^0$ are of interest in analyzing both the actual connectedness among individuals and the joint network effect.
We let $w_{jj}=0$ and assume $\delta^0_{jj}=0$ for all $j$. It is worth noting that endogeneity is a concern, since the inclusion of $y_{k,t}$ ($k\neq j$) induces simultaneity in the structural equation system. To handle the simultaneity bias, instrumental variables (denoted by $z_{j,t}$) are needed. For example, lags $y_{j,t-1},y_{j,t-2}\ldots$ are commonly used in practice. We shall further assume that ${\varepsilon}_{j,t}$ are martingale difference sequences with respect to a suitable filtration, as defined below, and allow for temporal and spatial dependencies in the observed data sample (see \hyperref[A_dgp1]{(A1)(i)}, \hyperref[A_dan]{(A2)} and \hyperref[A_error]{(A3)}).
As a practical example, in de2018recovering, $y_{j,t}$ refers to the state tax liabilities for state $j$ in year $t$, $w_{jk}$ is observed as some known geographic measurement of neighborhood, and $\delta^0_{jk}$ contributes to the measurement deviations. In this case, the overall network effect, i.e., $\rho^0 w_j + \delta^0_j$ is interpreted as an overall economic measurement of the connections. On this basis, the social network effect of tax competition is analyzed.
In addition, we can expand the model by including equation-specific covariates $u_{j,t}\in\mathbb{R}^{d_j}$ whose dimension may grow with the sample size:
The compact form of the model is given by:
where $y_t=(y_{1,t}, \ldots, y_{p,t})^\top$, $u_t=(u_{1,t}^\top,\ldots,u_{p,t}^\top)^\top$, and ${\varepsilon}_t=({\varepsilon}_{1,t}, \ldots, {\varepsilon}_{p,t})^\top$. In this expression, $W$ and $\Delta^0$ are $p\times p$ matrices, with the $j$-th row of $W$ being $w_j^\top$ and the $j$-th row of $\Delta^0$ being $\delta_j^{0\top}$. If the covariates $u_t$ are exogenous, the transformed covariates $Wu_t$ and $W^2u_t$ are also commonly used as instrumental variables.
Following the spatial econometrics literature, we assume that $|(\rho^0W+\Delta^0)^t|_\infty\leq|c|^t$ with some $|c|<1$ to ensure the stationarity of the model. Without loss of generality, for identification purposes, suppose it is known that there exist $j^*,k^*$ ($k^*\neq j^*$) such that $W_{j^*k^*}\neq0$ and $\Delta^0_{j^*k^*}=0$. This implies that at least one of the non-zero actual links can be correctly specified by the observed linkage. This assumption ensures that the regression does not suffer from multicollinearity.
In reality, $W$ might be either sparse or dense. On the other hand, it is noted in the literature that the classical spatial estimator for $\rho^0$, such as the IV estimator, would not be consistent if the misspecification error $\Delta^0$ is too dense; see recent works by lewbel2023estimating,lewbel2023ignoring. We therefore posit that $\Delta^0$ is approximately sparse, though the observed or actual network structure may not necessarily be sparse.
When multiple options for the pre-specified matrix $W$ are available, a linear combination of the potential matrices $W_i$, $i=1,\ldots,M$, can be incorporated into the model. Such a generalization has been considered in articles such as lam2014regularization,higgins2023shrinkage, with $M$ increasing as $n$ grows. In this case, a regularized estimation can be performed on the weights associated with $W_i$'s, and the sparse weights would be included as part of the unknown parameters in our framework.
In our empirical section (ref), we attempt to quantify the spillover effect among individual stocks, where $y_t$ denotes a vector of stock returns, $W$ is a network matrix corresponding to the common shareholder information, and $\rho^0$ measures the joint network effect. The purpose of this application is to understand the overall network effect among firms and to uncover the latent links.
It is worth noting that the spatial panel network model we have discussed fits within the framework of high-dimensional regression equations, potentially involving endogeneity. In Appendix (ref), we present a general model framework that encompasses many examples in panel or longitudinal data analysis. For instance, the general model can be dynamic, allowing for the inclusion of lagged values of $y_{j,t}$ in the covariates. The primary theorems presented in Section (ref) and Appendix (ref) are applicable to the estimator of the general model when using linear or nonlinear moments.
In this subsection, we outline the estimation steps for the DRGMM estimator, which include obtaining a preliminary estimator using the Dantzig selector and the subsequent debiasing procedure, allowing us to perform inference on the parameters of interest.
For each equation $j=1,\ldots,p$, let $x_{j,t}$ and $\vartheta_j^0$ collect the regressors and the corresponding coefficients respectively. Recall the existence of indices $(j^*,k^*)$, where $j^*\neq k^*$. Specifically, for $j\neq j^*$, we have: $$x_{j,t}=(w_j^\top y_t, y_{t}^\top, u_{j,t}^\top)^\top, \quad \vartheta_j^0=(\rho^0,\delta_j^{0\top},\beta^{0\top})^\top;$$ for $j=j^*$, we have: $$x_{j,t}=(w_j^\top y_t, y_{t,-k^*}^\top, u_{j,t}^\top)^\top,\quad \vartheta^0_j=(\rho^0,\delta_{j,-k^*}^{0\top},\beta^{0\top})^\top,$$ where $y_{t,-k^*}$ denotes the subvector of $y_t$ obtained by excluding the $k^*$th element $y_{t,k^*}$, and similarly for $\delta^0_{j,-k^*}$. With these notations, we can rewrite the model in (ref) in the form of $y_{j,t}=x_{j,t}^\top\vartheta_j^{0}+{\varepsilon}_{j,t}$. Let $K_j$ denote the dimension of $x_{j,t}$. Define $\theta^0 = (\rho^0, \delta_1^{0\top},\ldots,\delta_p^{0\top},\beta^{0\top})^\top\in\mathbb{R}^K$ to collect all the parameters in the model. Note that, given $\delta_{j^*k^*}^0=0$, the parameter sets $(\vartheta_1^0,\ldots,\vartheta_p^0)$ and $\theta^0$ contain the same unknown parameters. We shall estimate $\theta^0$ under the assumption that it is sparse.
Due to the endogeneity in the structural model, we introduce the instrumental variables $z_t=[z_{j,t}]_{j=1}^p\in\mathbb{R}^q$, where $q=\sum_{j=1}^p q_j\geq K$, to construct the moments. Specifically, $z_{j,t}\in\mathbb{R}^{q_j}$ contains the instrumental variables for the $j$-th equation, ensuring that $\mathop{\mbox{\sf E}}({\varepsilon}_{j,t}|z_{j,t})=0$. Here, the notation $[A_j]_{j=1}^p$ indicates that we stack $A_j$ by rows over $j=1,\ldots,p$.
For each $j=1,\ldots,p$, we define a vector-valued score function $g_j(D_{j,t},\theta)$ that maps $\mathbb{R}^{K_j+q_j}\times\mathbb{R}^K$ to $\mathbb{R}^{q_j}$, where $D_{j,t}\stackrel{\mathrm{def}}{=}(x_{j,t}^\top, z_{j,t}^\top)^\top$. For the case with linear moments, the score function is given by $g_j(D_{j,t},\theta)=z_{j,t}{\varepsilon}_j(D_{j,t},\theta)$, where $ {\varepsilon}_j(D_{j,t},\theta)=y_{j,t}-x_{j,t}^\top\vartheta_j$. Thus, the moment functions mapping $\Theta\subseteq\mathbb{R}^K$ to $\mathbb{R}^{q_j}$ are: $$g_j(\theta)=\mathop{\mbox{\sf E}} g_j(D_{j,t},\theta)= \mathop{\mbox{\sf E}} [z_{j,t}(y_{j,t}-x_{j,t}^\top\vartheta_j)],$$ and we have $g_j(\theta^0)=0$. By stacking the moment functions across equations, we get $g(\theta)=[g_j(\theta)]_{j=1}^p$. The empirical counterpart is computed as: $$\hat g(\theta) =[\mathop{\mbox{\sf E}}{_n}g_j(D_{j,t},\theta)]_{j=1}^p=[\mathop{\mbox{\sf E}}{_n}\{z_{j,t}(y_{j,t}-x_{j,t}^\top\vartheta_j)\}]_{j=1}^p.$$ Additionally, the covariance matrix of the score functions is defined as $$\Omega_{q\times q} \stackrel{\mathrm{def}}{=} \frac{1}{n}\mathop{\mbox{\sf E}}\Big[\Big\{\sum_{t=1}^n g(D_t,\theta^0)\Big\}\Big\{\sum_{t=1}^n g(D_t,\theta^0)\Big\}^\top\Big],$$ where $D_t=[D_{j,t}]_{j=1}^p$ and $g(D_t,\theta)=[g_j(D_{j,t},\theta)]_{j=1}^p\in\mathbb{R}^{q}$. In our case, this simplifies to $\Omega=\mathop{\mbox{\sf E}}\big[[z_{j,t}{\varepsilon}_{j,t}]_{j=1}^p([z_{j,t}{\varepsilon}_{j,t}]_{j=1}^p)^\top\big]$.
Suppose the parameter vector $\theta^0\in\mathbb{R}^K$ is partitioned into two parts: the parameters of interest $\theta_1^0\in\mathbb{R}^{K^{(1)}}$ and the nuisance parameters $\theta_2^0\in\mathbb{R}^{K^{(2)}}$, where $K^{(1)}+K^{(2)} =K$. In this context, we are primarily interested in $\theta_1^0=(\rho^0, \delta_1^{0\top},\ldots,\delta_p^{0\top})^\top$, which includes the spatial autoregressive parameter and the misspecification errors of the network structure. Meanwhile, the coefficients on the control variables, denoted by $\theta_2^0=\beta^0$, are treated as nuisance parameters. Let $G_1$ and $G_2$ denote the Jacobian matrices of the moment function $g(\theta)$ with respect to $\theta_1$ and $\theta_2$, respectively. Specifically, since $\theta_1$ contains both common and equation-specific parameters, $G_1\stackrel{\mathrm{def}}{=}\partial_{\theta_1^\top}g(\theta_1,\theta_2^0)|_{\theta_1=\theta_1^0}$ can be decomposed as $G_1=(G_{11},G_{12})$, where $G_{11}$ is a $q\times1$ vector given by $-[\mathop{\mbox{\sf E}}(z_{j,t}w_j^\top y_t)]_{j=1}^p$, and $G_{12}$ is a $q\times (K^{(1)}-1)$ block diagonal matrix whose $j$th block is $-\mathop{\mbox{\sf E}}(z_{j,t}y_t^\top)$.{\linespread{1}\footnote{Under the assumption that $\delta^0_{j^*k^*}$ is known to be zero, we can simplify the parameter vector $\theta^0_1$ by excluding $\delta^0_{j^*k^*}$. Consequently, the corresponding column in the Jacobian matrix $G_{12}$ associated with $\delta^0_{j^*k^*}$ (i.e., the ($pj^*+k^*$)-th colum) should also be removed.}} Additionally, $G_2\stackrel{\mathrm{def}}{=}\partial_{\theta_2^\top}g(\theta_1^0,\theta_2)|_{\theta_2=\theta_2^0}$ is a $q\times K^{(2)}$ matrix given by $-[\mathop{\mbox{\sf E}}(z_{j,t}u_{j,t}^\top)]_{j=1}^p$. Other ways to partition the parameters are also possible, with the expressions for $G_1$ and $G_2$ adjusted accordingly.
The DRGMM estimator procedure will be carried out in two steps:
Our estimation procedure is designed for settings where $n,p,K_j$, and $q_j$ (and thus $K$ and $q$) can all diverge. It includes a special case of many IV problems with $n,q\to\infty$ while the number of unknown parameters is fixed. In this scenario, regularization on the parameters is not required in the first estimation step. For instance, in our supplementary simulation study in Appendix (ref), we consider the Arellano-Bond (AB) estimator for dynamic panel models, where an excessive number of instruments is used to estimate two parameters. We use the conventional AB estimator as the preliminary estimator, which is then refined in a subsequent debiasing step to address overidentification with optimal moment selection.
It is worth noting that in the high-dimensional setting ($q>n$), $\hat\Omega$ is singular due to the rank deficiency, necessitating the use of a regularized estimator for the precision matrix. Specifically, we consider the constrained $\ell_1$-minimization for inverse matrix estimation (CLIME; see cai2011constrained). In Appendix (ref), we will present a feasible debiased estimator $\check\theta_1$ that incorporates approximate inverse matrices. The convergence rates of the estimators involved in addressing the rank deficiency are analyzed in several auxiliary lemmas in the same appendix.
To provide additional clarity on the debiasing step, in Appendix (ref), we establish a link between our debiased estimator and the Two-Stage Least Squares (2SLS) estimator in a low-dimensional framework, where the number of unknown parameters and moment conditions remain fixed.
In Section (ref), we will demonstrate that the debiased estimator $\check\theta_1$ is asymptotically unbiased and Gaussian. This allows us to perform simultaneous inference on the parameters of interest.
In this section, we present the theoretical foundation of the proposed estimator for the case of linear moments. Specifically, Section (ref) focuses on the consistency of the preliminary GDS estimator, $\hat\theta$, in step 1, while Section (ref) examines the inference procedure for the final DRGMM estimator, $\check\theta_1$, for the parameters of interest. Extensions of the main theory to the case of nonlinear moments are discussed in Section (ref) of the Appendix.
Throughout this section, we impose the following assumptions and definitions:
In \hyperref[A_dgp1]{(A1)(i)}, we allow for overlap in the innovations $\zeta_{jk,t},\xi_{jm,t},\eta_{j,t}$ as long as the exogeneity condition $\mathop{\mbox{\sf E}}(z_{j,t}{\varepsilon}_{j,t})=0$ is satisfied. \hyperref[A_dan]{(A2)} assumes a sufficient decay rate of dependency. In the main text of this paper, we focus on the weak dependence case with $\varsigma>1/2-1/r$. In the detailed proofs in the appendix, we will discuss how the rates adapt to the case of strong dependence. It is worth noting that \hyperref[A_dan]{(A2)}, together with the stationary condition \hyperref[A_dgp2]{(A1)(ii)}, implies that the dependence adjusted norm for the transformed covariates $\|x_{jk,\cdot}\|_{r,\varsigma}$ is also finite.
Assumption \hyperref[A_error]{(A3)} restricts the dependence structure of the error term by assuming it follows a martingale difference sequence (m.d.s.) with respect to the filtration $\mathcal F_{t-1}$. While this rules out serial correlation, it remains reasonable as our general modeling framework accommodates dynamics through the inclusion of sufficiently many lags. Due to the m.d.s. nature of the error term, the long-run variance of the score functions need not be considered in forming $\Omega$ for debiasing. Additionally, we impose some structure on the conditional variance-covariance matrix to simplify the derivation. However, this setting could be extended to handle more complex structures, such as serial correlations, unobserved heterogeneity, and factor structures. See Appendix (ref) for further discussion.
\hyperref[ES]{(A4)} focuses on the sparsity of the true parameter $\theta^0$, which is the assumption we primarily rely on in demonstrating the main theorems. This condition can be extended to the case of approximate sparsity, a more general assumption in the literature on high-dimensional data analysis. In Appendix (ref), we will derive the estimation error bounds under the approximate sparsity assumption, taking into account the approximation error.
\hyperref[A_tune]{(A5)} ensures that $\theta^0$ is feasible for the problem in (ref) with probability at least $1-\alpha$. \hyperref[A_id]{(A6)} is an identification assumption that is crucial for ensuring consistency. In Appendix (ref), we discuss the conditions required to validate this assumption.
In order to establish the consistency of the GDS estimator $\hat\theta$, we need to derive the error bound for $|\hat{\theta}- \theta^0|_a$ for $a=1$ or $2$, and analyze the convergence rate. Under the identification condition \hyperref[A_id]{(A6)}, the error bound for $|\hat{\theta}- \theta^0|_a$ follows from the error bound for $| g(\hat\theta) - g(\theta^0)|_{\infty}$ (we will elaborate on this argument in the proof of Theorem (ref)). Using the identity $g(\theta^0) = 0$, we can bound $|g(\hat{\theta})- g(\theta^0)|_{\infty}$ as follows:
Recalling the definition of the GDS estimator, we have $|\hat g(\hat{\theta})|_{\infty}\leq\lambda_n$. Let $\mathcal{R}(\theta^0) \stackrel{\mathrm{def}}{=} \{\theta \in \Theta: |\theta|_1 \leq |\theta^0|_1\}$ denote the restricted set. As a consequence of \hyperref[A_tune]{(A5)}, we could have $\hat\theta\in\mathcal R(\theta^0)$ with probability at least $1-\alpha$. The remaining task is to demonstrate the concentration result, i.e., to show that: $$\sup\limits_{\theta \in \mathcal{R}(\theta^0)}|\hat g (\theta)- g(\theta)|_{\infty}\leq\epsilon_n$$ holds with probability approaching 1, for a sequence of positive constants $\epsilon_n\downarrow0$ as $n\to\infty$.
We focus on cases with linear moments, where $g(\theta)=G\theta+g(0)$ and $\hat g(\theta)=\hat G\theta +\hat g(0)$, with $G=\partial_{\theta^\top}g(\theta)$ and $\hat G=\partial_{\theta^\top}\hat g(\theta)$ being independent of $\theta$. It follows that
where $|\cdot|_{\max}$ denotes the element-wise max norm of a matrix.
To derive the convergence rate, we need to analyze the rates of $|\hat G-G|_{\max}$ and $|\hat g(0)- g(0)|_{\infty}$ by applying the concentration inequality in Lemma (ref). For this purpose, we define the following quantities:
which are all bounded by constants for some $r\geq4$ and $\varsigma>0$ according to \hyperref[A_dan]{(A2)}. Additionally, we define $\Phi_{r,\varsigma}^{yz} = \max\limits_{1\leq j\leq p,1\leq m\leq q_j} \|y_{j,\cdot} z_{jm,\cdot}\|_{r,\varsigma}$.
For each equation $j$, we aggregate the dependence adjusted norm of the vector of processes $x_{j,t}$ as follows: $$\||x_{j,\cdot}|_{\infty}\|_{r,\varsigma} = \sup_{d\geq 0}(d+1)^{\varsigma}\Delta_{d,r,j},\quad\Delta_{d,r,j}=\sum_{t=d}^\infty\||x_{j,t}- x_{j,t}^{\ast}|_\infty\|_r.$$ This is in comparison to the dependence adjusted norm for a univariate process as in Definition (ref). Similarly, we define $\||x_{j,\cdot}z_{jm,\cdot}|_{\infty}\|_{r,\varsigma}$. Additionally, we aggregate over $j=1,\ldots,p$ by: $$\Big\|\max_{1\leq j\leq p}|x_{j,\cdot}|_{\infty}\Big\|_{r,\varsigma} = \sup_{d\geq 0}(d+1)^{\varsigma}\Delta_{d,r},\quad \Delta_{d,r}=\sum_{t=d}^\infty\Big\|\max_{1\leq j\leq p}|x_{j,t}- x_{j,t}^{\ast}|_\infty\Big\|_r.$$ The definition for $\Big\|\max\limits_{1\leq j\leq p,1\leq m\leq q_j}|x_{j,\cdot}z_{jm,\cdot}|_{\infty}\Big\|_{r,\varsigma}$ follows similarly.
In the case where the dependence adjusted norms involved in $b_n$ and $b'_n$ are bounded by constants, and assuming that $n^{-1/2+1/r}(\log P_n) =\mathcal{O}(1)$ for sufficiently large $r$, we have the concentration rate $$\epsilon_n \lesssim (s+1)n^{-1/2}(\log P_n)^{1/2},$$ which matches the rate shown in Lemma 3.3 of belloni2018high for i.i.d. data.
Combining the results from Lemma (ref) with the identification condition \hyperref[A_id]{(A6)}, we obtain the bound on the estimation error of the GDS estimator. The rate of consistency is stated in the following theorem.
According to Corollary 5.1 of lasso2018, the order of $\lambda_n$ is given by $$n^{-1}\max_{1\leq j\leq p,1\leq m\leq q_j}\Big(\|z_{jm,\cdot}{\varepsilon}_{j,\cdot}\|_{2,\varsigma}(n\log q)^{1/2}\vee\|z_{jm,\cdot}{\varepsilon}_{j,\cdot}\|_{r,\varsigma}(nq)^{1/r}\Big).$$ In the case where$(nq)^{1/r} \lesssim (n\log q)^{1/2}$, we have $\lambda_n \lesssim n^{-1/2}(\log q)^{1/2}$. This implies that if $r$ is sufficiently large, $q$ can diverge as a polynomial rate of $n$ (a better dimension allowance for $q$ is possible under stronger exponential moment conditions; see Comment 5.5 in lasso2018). Consequently, assuming $\max\limits_{1\leq j\leq p,1\leq m\leq q_j}\|z_{jm,\cdot}{\varepsilon}_{j,\cdot}\|_{r,\varsigma}$ is bounded by a constant, we have: $$d_{n,a} \lesssim (s+2)s^{1/a}n^{-1/2} (\log P_n)^{1/2},$$ which is of the same order as the rate for the i.i.d. case studied in Theorem 3.1 of belloni2018high.
In this subsection we show the asymptotic properties of the debiased estimator $\check\theta_1$ obtained in the second step. In particular, we provide a key representation that linearizes the estimator, facilitating the application of a high-dimensional Gaussian approximation theorem for inference.
Define ${A}\stackrel{\mathrm{def}}{=}{G}_1^{\top}{\Omega}^{-1}(\mathbf I_q - {G}_2P({\Omega},{G}_2))$ and $B\stackrel{\mathrm{def}}{=}(AG_1)^{-1}$, where $P(\Omega, G_2) \stackrel{\mathrm{def}}{=} (G_2^{\top}\Omega^{-1} G_2)^{-1}G_2^{\top}\Omega^{-1}$. Consider estimators of $A$ and $B$, denoted by $\hat{A}$ and $\hat{B}$. More details about the construction of these estimators are discussed in Section (ref) and Appendix (ref).
We shall analyze the accuracy of estimator $\check\theta_1$. Observe that
where $r_n=r_{n,1}+r_{n,2}$, and $$ r_{n,1}=(\mathbf I - \hat B\hat A\hat G_1)(\hat\theta_1-\theta_1^0),\, r_{n,2}=(BA - \hat B\hat A)\hat g(\theta^0). $$ Note that, due to the Neyman orthogonality property, the term $r_{n,1}$ is expected to be small. Under mild conditions, the last term $r_{n,2}$ is also expected to vanish. By applying the triangle inequality and H\"{o}lder's inequality, we have the following bounds for the two terms $r_{n,1}$ and $r_{n,2}$, respectively:
The linearized representation in (ref) shows that the debiased estimator $\check\theta_1$ can be expressed by the true parameter $\theta_0$ plus a weighted empirical moment function evaluated at $\theta_0$, along with an approximation error. Consequently, relying on a high-dimensional Gaussian approximation of the leading term $BA\hat{g}({{\theta}^0})=(A{G}_1)^{-1}A \hat{g}({{\theta}^0})$, as will be discussed in Section (ref), valid inference can be conducted, provided that the linearization errors are asymptotically negligible in the sense that $|r_n|_{\infty}$ is of small order. We now present a theorem for the linearization of the debiased estimator.
The proof of this theorem and the detailed rate of $|r_n|_\infty$ are deferred to Appendix (ref). In particular, we will discuss the rate specifically under the special case where all the dependence adjusted norms involved are bounded by constants in Remark (ref). To enable valid inference through the Gaussian approximation on the leading term, we require that the linearization errors be sufficiently small, ensuring that $\sqrt{n}|r_n|_\infty=\mbox{\tiny $\mathcal{O}$}_\mathrm{P}(1)$. This condition imposes restrictions on the allowed dimensionality and sparsity relative to $n$, under mild assumptions.
In this subsection, we cite a high-dimensional Gaussian approximation theorem to facilitate the simultaneous inference of the parameters. The theorem is adapted from ZW15gaussian. Specifically, we focus on testing the hypothesis $H_0:\theta_{1,k}^0=0,\forall k\in \mathcal S$, where $\mathcal S\subseteq \{1,\ldots,K^{(1)}\}$, and $\theta_{1,k}^0$ denotes the $k$-th element of the vector $\theta_1^0$. To proceed with this inference, we first revisit some key definitions from Section (ref).
For the case of linear moments, the score functions evaluated at the true parameters are given by $g_j(D_{j,t},\theta^0)= z_{j,t}{\varepsilon}_{j,t}$, where $D_{j,t}\stackrel{\mathrm{def}}{=}(x_{j,t}^\top, z_{j,t}^\top)^\top$. Let $D_t=[D_{j,t}]_{j=1}^p$ and $g(D_t,\theta)=[g_j(D_{j,t},\theta)]_{j=1}^p$. Define the vector $\mathcal{G}_t=(\mathcal{G}_{k,t})_{k\in\mathcal S}$, where $\mathcal{G}_{k,t}=-\zeta_kg(D_t,\theta^0)$, and $\zeta_k$ is the $k$-th row of the matrix $(AG_1)^{-1}A$. Assuming $|(AG_1)^{-1}A|_\infty$ is bounded by a constant, for any $k\in\mathcal S$, the dependence adjusted norm of $\mathcal{G}_{k,t}$ is bounded by $$\|\mathcal{G}_{k,\cdot}\|_{r,\varsigma}\lesssim\max_{1\leq j\leq p,1\leq m\leq q_j}\|z_{jm,\cdot}\|_{2r,\varsigma}\|{\varepsilon}_{j,\cdot}\|_{2r,\varsigma}.$$
For simultaneous inference, we allow the number of parameters being tested, i.e., the cardinality $|\mathcal S|$, to increase as $n\to\infty$. Specifically, we consider a polynomial growth rate, $|\mathcal S|=n^c$ for some $c>0$. The admissible growth rate is specified in the following assumption:
We now state the Gaussian approximation results as follows. Denote by $c_\alpha$ the $(1-\alpha)$ quantile of the $\max_{k\in S}|\mathcal Z_k|$, where $\mathcal Z_k$ are the standard normal random variables. Let $\sigma_k^2$ be the $k$-th diagonal element of the covariance matrix $(AG_1)^{-1}A\Omega A^\top\{(AG_1)^{-1}\}^\top=(AG_1)^{-1}$. Under \hyperref[A_Srate]{(A7)} and the same conditions as in Theorem (ref), assume that there exists a constant $C>0$ such that {$\min\limits_{k\in \mathcal S}\operatorname{Var}\big(n^{-1/2}\sum_{t=1}^n\mathcal{G}_{k,t}\big)\geq C$}. Then, we have
The conclusion also holds when $\sigma_k$ is replaced by a consistent estimator $\hat\sigma_k$. Consequently, for each $k\in\mathcal S$, we can construct the two-sided $(1-\alpha)$ confidence interval using asymptotic normality as:
Based on the Gaussian approximation results, the multiplier bootstrap can be employed to determine the critical value. To account for temporal dependence, we adopt a block multiplier bootstrap procedure using non-overlapping blocks. Define the vector $\widehat{\mathcal T}=(\widehat{\mathcal T})_{k\in\mathcal S}$, where
$\hat\zeta_k$ is the $k$-th row of the matrix $(\hat A\hat G_1)^{-1}\hat A$ and $e_i$ are independently drawn from $\mathop{\mbox{\sf N}}(0, 1)$. Here, $l_n$ and $b_n$ denote the numbers of blocks and block size, respectively, with $b_n=\lfloor n/l_n\rfloor$. To ensure the validity of the multiplier bootstrap, as shown in the following theorem, we assume that the block size grows at a polynomial rate such that $b_n=\mathcal{O}(n^\eta)$ for some $0<\eta<1$. Intuitively, a larger block size is needed to effectively capture the dependency structure, while sufficiently many blocks are required for robust approximation of the bootstrapped statistics. To address this trade-off, a set of accompanying conditions is imposed, further narrowing the admissible range of $\eta$ to determine the optimal $b_n$ in (ref), as detailed in the proof. This range is influenced by its interplay with $r$ and $\varsigma$ (satisfying \hyperref[A_dan]{(A2)}) and the size of $|\mathcal S|$.
As a result of Theorem (ref), for each $k\in\mathcal S$, we can construct the two-sided $(1-\alpha)$ bootstrap confidence interval as:
In this section, we illustrate the finite sample properties of our proposed methodology across different simulation scenarios. We first present results for the primary example of spatial panel networks discussed in Section (ref), while Appendix (ref) focuses on dynamic linear panel models.
Consider the spatial panel network model with covariates:
where $h^0_j=(h^0_{j1},\ldots,h^0_{jp})^\top$ and $h^0_{jk}$ ($k\neq j$) represents the actual, unobserved peer effect of unit $k$ on unit $j$. Our objective is to estimate the joint network effect $\rho^0$, recognizing that $h_j^0$ may be misspecified as an observed network structure $w_j=(w_{j1},\ldots,w_{jp})^\top$ for all $j=1,\ldots,p$. The model can then be rewritten as: $$y_{j,t}=\rho^0 w_j^\top y_{t} + \rho^0 \delta_j^{0\top}y_t + \beta^{0\top} u_{j,t}+{\varepsilon}_{j,t},$$ where the vectors $\delta^0_j=(h_j^{0\top} - w_j^\top)$, $j=1,\ldots,p$, capture the misspecification errors of the network structure.
We randomly generate the actual links by using independent Bernoulli random variables, each with a probability of 0.5 of equaling one. Additionally, we set $h^0_{jj}=0$ and apply normalization to $h^0_j$ for each $j=1,\ldots,p$. We assume that misspecification occurs randomly with a probability of 0.2 when an actual link is non-zero; that is, $h^0_{jk}\neq0$ but $w_{jk}=0$.
To incorporate the dependency, we let the instrumental variables $Z_{j,t}\in\mathbb{R}^{q_j}$ for $j=1,\ldots,p$, follow a linear process such that $Z_{j,t}=\sum_{\ell=0}^{\infty}A^j_\ell\xi_{j,t-\ell}$, where $A^j_\ell=(\ell+1)^{-\tau-1}M^j_\ell$, and $M^j_\ell$ are independently drawn from Ginibre matrices (i.e., all entries of $M^j_\ell$ are i.i.d. $\operatorname{N}(0,1)$). In practice, the sum is truncated to $\sum_{\ell=0}^{500}A^j_\ell\xi_{j,t-\ell}$. We set $\tau=1.0$ for weaker dependence and $\tau=0.1$ for stronger dependence. For the $q_j$-dimensional vector $\xi_{j,t}$, we define each element as $\xi_{jk,t}=e_{jk,t}(0.8e_{jk,t-1}^2+0.2)^{1/2}$ for $k=1,\ldots,q_j$, where $e_{jk,t}$ are i.i.d. and follow a scaled $t(8)$-distribution: $t(8)/\sqrt{8/(8-2)}$, with $t(8)$ being the Student's $t$-distribution with 8 degrees of freedom.
Then, for each $j=1,\ldots,p$, we generate the $d$-dimensional covariates $u_{j,t}$ as follows: $$u_{j,t}=\pi^\top Z_{j,t} + v_{j,t},$$ where the $q_j\times d$ matrix $\pi$ is defined as $\pi=(3+3\kappa^{q_j/3})^{-1/3}(\iota_3\otimes \mathbf{I}_{q_j/3})$, with $d=q_j/3$, $\iota_3$ being a $3\times1$ vector of ones, and $\kappa=0.5$. We let $\beta^0=(10,10,10,10,10,5,5,5,1,1,0_{d-10}^\top)^\top$. The errors ${\varepsilon}_{j,t}$ and $v_{j,t}$ are generated independently from standard normal distribution.
We consider two cases: $p=30,d=30,n=100$ and $p=50,d=50,n=200$, where the total number of parameters, $p^2+d+1$, amounts to 931 and 2,551, respectively. The total number of moment conditions, $q=\sum_{j=1}^pq_j$, is 2,700 for the first case and 7,500 for the second. Specifically, we focus on $\rho^0$, $\beta^0$, and $\tilde\delta^0$, which includes the first 50 elements of the stacked vector $[\tilde\delta^0_j]_{j=1}^p$ as parameters of interest. Here, $\tilde\delta^0_j$ is defined as a subvector of $\delta^0_j$ with elements known to be zero removed. These removed zero elements correspond to non-zero $w_{jk}$ values, which are assumed to be correctly specified in this setting.
To assess the estimation accuracy of our proposed two-step method, we compute the absolute deviation for estimating $\rho^0$ and the estimation error for the vectors $\beta^0$ and $\tilde\delta^0$, measured by the Euclidean norm. These calculations are performed on estimators both with and without applying the debiasing step, namely, the DRGMM and GDS estimators. In the first step, we use penalty that is independent of the design matrix. Specifically, we set $\lambda_n=\Phi^{-1}(1-0.1/(2q))\max\limits_{1\leq j\leq q}\hat\sigma_j^2/\sqrt{n}$, given that $\hat g_j(\theta)$ asymptotically follows $\operatorname{N}(0,\sigma^2_j/n)$ for $j=1,\ldots,q$. This choice is intentionally conservative to mitigate the risk of overfitting, following belloni2018high. When the debiasing step is applied, we treat $\rho^0$, $\beta^0$, and $\tilde\delta^0$ as the parameters of interest respectively. For the convenience of comparison, we present the estimation errors as ratios, which measure the relative difference between the results obtained using the DRGMM and GDS estimator. In particular, a ratio smaller than 1 indicates better performance when the debiasing step is applied. The results, summarized in Tables (ref) and (ref), are aggregated over 500 replications using both the mean and the median.
Additionally, we evaluate the inference performance by examining the empirical power and size of the confidence intervals (with a nominal confidence level of 95%) constructed using the limiting distribution theory outlined in Section (ref). Specifically, the average rejection rate of the null hypotheses for the truly zero components reflects size performance, while the testing power is evaluated for the truly non-zero components. Inference results are reported separately for the structural parameters $(\rho^0,\beta^0)$ and for the network structure $\tilde\delta^0$. For comparison, the average false positive rate for truly zero parameters and the average true positive rate for truly non-zero parameters under the GDS estimator are also reported to assess the necessity of uniform inference via debiasing. The results, based on 500 simulations, are presented in Tables (ref) and (ref).
From Tables (ref) and (ref), it is evident that debiased regularization significantly outperforms the one-step GDS estimator in estimating the structural parameters, particularly when a stronger network effect is observed in $\rho^0$. Our proposed DRGMM estimator performs well in recovering the true network structure, with its superiority becoming more pronounced for larger-scale networks and under stronger dependency. Overall, we observe that the estimation errors are robust across simulations.
From Tables (ref) and (ref), we find that inference after applying the debiasing step provides size control close to the nominal level and high empirical power in most cases. While the Dantizig selection successfully detects the truly non-zero structural parameters, it is not reliable for recovering the latent network structure. Notably, our proposed method effectively avoids an excess of false positives, which can occur with the one-step regularized selection. These results confirm the necessity of uniform inference on parameters, including both truly zero and non-zero elements in the network structure.
In this section, our proposed methodology is employed to study the spatial network effect of stock returns. We use the public cross-ownership information as the pre-specified social network structure zhu2019network; however, there might be misspecification in the network since some of the cross-shareholder information is not published. Our purpose is to analyze the network effect and simultaneously recover the unobserved linkages.
Our empirical illustration is carried out on a dataset consisting of 100 individual stocks traded on the Chinese A-share market (Shanghai Stock Exchange and Shenzhen Stock Exchange), spanning 14 sectors as defined by the Industry Classification Guidelines of the China Securities Regulatory Commission. The data covers the period from January 2, 2019 to December 31, 2019 (i.e., $244$ trading days). Daily stock returns and annual cross-ownership data were sourced from the \href{https://www.wind.com.cn}{Wind Data Service}.
The spatial network model is specified as follows:
where $j=1,\ldots,p$ indexes the individual stocks, $y_t=(y_{1,t},\ldots,y_{p,t})^\top$ represents the daily log returns, and $u_{j,t}$ denotes the daily turnover ratio (trading volume divided by shares outstanding), which is used as a firm-specific control variable. We assume that ${\varepsilon}_{j,t}$ satisfies assumption \hyperref[A_error]{(A3)}, with $\mathop{\mbox{\sf E}}({\varepsilon}_{j,t}{\varepsilon}_{j',t}|\mathcal F_{t-1})=0$ for $j\neq j'$. An unobserved individual effect, $\alpha_j$, is included to account for potential serial correlation in the error term. Within our estimation framework, these fixed effects can be treated as equation-specific intercepts during estimation.
The term $w_{jk}$ represents the public cross-ownership between stock $k$ and $j$, defined as $w_{jk}=1$ if company $j$ holds shares in company $k$ based on available information, and $w_{jk}=0$ otherwise. The resulting network structure for $w_{jk}$ ($j,k=1,\ldots,p$) is illustrated in Figure (ref), where the stocks are grouped by sector. Notably, cross-ownership relationships are observed across sectors.
It is possible for $w_{jk}=0$ while $h^0_{jk}\neq0$ if some shareholders of company $j$ are not publicly disclosed. We set $w_{jj}=h^0_{jj}=0$. Without loss of generality, we assume that misspecification errors occur only when the actual link is non-zero; specifically, cases where $w_{jk}\neq0$ while $h^0_{jk}=0$ are excluded. Our goal is to estimate the network effect $\rho^0$ and the misspecification errors $\delta^0_{jk}=h^0_{jk}-w_{jk}$ using our proposed approach. Ultimately, we aim to recover the latent linkages $h^0_{jk}$ based on inference results for the deviation $\delta^0_{jk}$, particularly in cases where $w_{jk}$ is observed to be zero.
In particular, the two-step DRGMM estimation procedure described in Section (ref) is applied, with $y_{t-1},y_{t-2}$ chosen as the instrumental variables. The resulting debiased estimators are $\check\rho=0.2214$ and $\check\beta =0.0012$, with standard errors of 0.0061 and 0.0001. Both $\rho^0$ and $\beta^0$ are found to be statistically significant. For comparison, we also fit the spatial autoregressive (SAR) model based solely on the observed network structure, using the same moment conditions (i.e., same instruments) for estimation. The estimated structural parameters $(\rho^0,\beta^0)$ are 0.3188 and 0.0014, with standard errors of 0.2181 and 0.0001, respectively. Notably, our proposed approach identifies a significantly stronger network effect.
Furthermore, it is of interest to test the latent network structure, and the inference theory based on DRGMM provides a formal framework for doing so. Following the discussion in Section (ref), we perform individual hypothesis tests on $H_0^{(j,k)}:\delta_{jk}^0=0$ if the preliminary estimator in step 1, $\hat\delta_{jk}$, is found to be non-zero and $w_{jk}$ is observed to be zero. A total of 60 edges are considered, and debiasing is applied to the entire vector. The recovered network structure after hypothesis testings is shown in Figure (ref), where a directed link from $k$ to $j$ indicates that $\delta^0_{jk}$ is significantly non-zero, implying that $h^0_{jk}$ should also be non-zero.
We find that the recovered network, which accounts for latent link structures, differs substantially from the pre-specified network. Notably, the finance and insurance sector emerges as the one with the highest outbound degree centrality, while the most intensive connections are directed towards the manufacturing sector. At the individual stock level, Ping An Bank Co., Ltd. (000001.SZ) from the finance and insurance sector has the highest outbound degree centrality, with a value of 23, while ZTE Corporation (000063.SZ) from the manufacturing sector exhibits the highest inbound degree centrality, with a value of 5. These results highlight the importance of addressing misspecified network links when analyzing risk channels and financial stability within a financial system.
We thank Lung-Fei Lee for prompting the impetus to explore this topic. We are also grateful to Tim Christensen, Aureo de Paula, Wolfgang H\"ardle, Elena Manresa, Gerard van den Berg, and Jeffrey Wooldridge for helpful discussions. In addition, we thank the Editor, one Associate Editor, and two referees for their valuable comments, which have significantly improved the paper. We remain responsible for any errors or omissions. Chen Huang acknowledges financial support from the Independent Research Fund Denmark through the Inge Lehmann Grant (1132-00019B). Weining Wang is supported through the project “IDA Institute of Digital Assets”, CF166/15.11.2022, financed under the Romania's National Recovery and Resilience Plan; and the Marie Sk\lodowska-Curie Actions under the European Union's Horizon Europe research and innovation program for the Industrial Doctoral Network on Digital Finance, Project No. 101119635. Lastly, we thank GPT-4 for proofreading assistance; all content was reviewed and edited by the authors, who take full responsibility for the final version of the manuscript.
\vskip 2em \centerline{ \bf Appendix} \vskip -1em \setcounter{subsection}{0} \vskip 2em