Dynamic Network Quantile Regression Model

Dynamic Network Quantile Regression Model

{\em JEL classification}: C32, C51, G17 \\ {\em Keywords}: Dynamic Network Quantile Regression Model, Simultaneous Network Endogeneity, IVQR Estimator, Quantile Connectedness.

Introduction

The topology of financial networks is central to the study of financial contagion and systemic risk, see fafchamps2007, acemoglu2015, hautsch2015 among others. Given the relevance of tail dependence for financial supervision and risk management betz2016, a joint analysis of network effect and tail dependence becomes more important because the implications derived from network models evaluated by conventional conditional mean estimators cannot necessarily be generalized to the tails. ando2021connect show that major adverse events are associated with an increase in average connectedness but that their effects on the tails significantly differ.

Quantile regression (QR) has been a powerful tool for characterizing the heterogeneous policy impacts and measuring tail-event driven risk (e.g. hardle2016). Following a seminal work by bassett1978, QR can be used to evaluate the entire range of the conditional distribution. Recently, the literature on quantile time series regression has been rapidly growing. koenker2006quantile propose a quantile autoregressive model while galvao2013 study QR in an autoregressive dynamic framework with exogenous stationary covariates. Following the analysis of quantile cointegration in xiao2009quantile, cho2015quantile bring QR to the autoregressive distributed lag (ARDL) model literature. The quantile ARDL process captures both the long-run and short-run relationships at any desired location in the conditional distribution. engle2004 advance a conditional autoregressive value at risk model whilst white2015 propose a vector autoregressive (VAR) model for analyzing quantile dynamics.

However, in most financial systems, multiple entities are often intertwined and interacted with each other, which can be represented as networks hautsch2015, hardle2016, chen2019. In this context, zhu2017 develop a network autoregression (NAR) model, which has gained great popularity in social network analysis. A number of extensions have been developed. zhu2020multivariate consider a multivariate spatial autoregression model. zhu2019portal investigate a screening method to select influential nodes. zhu2020nonconcave studies nonconvex penalized estimation methods in high-dimensional VAR models while zhu2018grouped extend the network VAR model to allow for group-specific parameters.

In particular, zhu2019 extend the NAR model by zhu2017 and propose a network quantile autoregression (NQAR) model in order to analyze tail dependency in a dynamic network with a large number of nodes. The NQAR model consists of a system of equations, in which a continuous response is related to its lagged connected nodes, the response of the same node in the previous time period as well as node specific characteristics in a quantile autoregression process. However, main weakness of the NQAR model lies in that it does not accommodate the contemporaneous impact of connected nodes, even though the simultaneous network/peer effects are pervasive in empirical studies liu2014,forni2010. If they are statistically significant, the estimation of the NQAR model would become inconsistent and misleading.

In this paper, as a main contribution, we extend the NQAR model and propose a general dynamic network quantile regression model (DNQR) by explicitly incorporating contemporaneous and lagged network effects of connected nodes as well as the impacts of node-sepcific variables and observed common effects. Notice, however, that the simultaneous network effects are inherently endogenous to the system, which leads to inconsistent estimates at conditional mean regression as well as in QR, see wuthrich2019, wuthrich2020, chernozhukov2020network. To cope with the endogeneity issue in the different contexts, many studies have attempted to apply the instrumental variable quantile regression (IVQR) estimation advanced by Chernozhukov2006, e.g. frolich2013, su2016 and machado2019.

To deal with this challenging issue we adopt the IVQR approach. The social network data are similar to the spatial data, in the sense that observations from connected users are correlated. This makes the spatial autoregressive model good candidate for network data analysis. However, our work can be regarded as a nontrivial extension of Su2011, who apply the IVQR approach to analyzing the cross-section data using a linear spatial autoregressive model, to a dynamic network quantile model with nodal heterogeneity and common factors, which can shed further lights on uncovering the complex tail dependency in dynamic networks with a large number of nodes and time periods. Our study is also closely related to the growing literature on the tail comovements in a complex financial system, see diebold2014, hautsch2014, white2015 and ando2021connect.

More importantly, we follow Jenish2012 and xu2015maximum, and develop the general asymptotic theory for the IVQR estimator by applying the spatial near epoch dependence (NED) of the underlying network processes. We derive the detailed conditions on the network processes in order to establish the consistency and asymptotic normality of the IVQR estimator. Via Monte Carlo simulations, we confirm that the finite sample performance of the IVQR estimator is satisfactory across quantiles and the different error distributions under the different network structures.

Finally, we demonstrate the utility of our approach with an application to the dataset on the stocks traded in NYSE and NASDAQ in 2016, through the common shareholder network constructed using the information on the common mutual fund ownership, e.g. Anton2014. We find that the contemporary network effects (measured by returns of connected stocks in the same period) are positive and significant and dominate all the other effects across all quantiles. Furthermore, they are stronger at the lower tails than at the upper tails, suggesting that the contemporaneous network effects should be explicitly and carefully analysed in the dynamic network quantile model.

This paper proceeds as follows. In Section (ref) we outline the DNQR model and derive the stationarity condition of the underlying network process. Section (ref) introduces the IVQR estimation and develops its asymptotic properties using the spatial NED approach. In Section (ref) we provide simulation results, showing that the IVQR estimation performs satisfactory. In Section (ref) the DNQR model is applied to the US financial market data. Section (ref) provides concluding remarks. The mathematical proofs and the additional simulation and empirical results are presented in the Online Appendix. The replication code can be found \href{https://github.com/ChaowenZheng/Dynamic-Network-Quantile-Model}{here} on GitHub.

Notations: For a vector $v = (v_1, \ldots, v_m)^\top \in \mathbb{R}^m$, we denote $|v|_k = (\sum_{i=1}^m |v_i|^k)^{1/k}$, $\|v\|_k= (\sum_{i=1}^m\mathop{\mbox{\sf E}}|v_{i}|^k)^{1/k}$, and $|v|_\infty = \max_{i \le m} |v_i|$, where $k$ is a positive integer, and $\mathop{\mbox{\sf E}}$ is the expectation operator. For any $n \times m$ matrix $A = (a_{i j})_{1\le i\le n, 1\le j \le m}$, we define the two norm and the max norm by $|A|_2 = \sup_{\{v\in \mathbb{R}^m ,|v|_2 = 1\}} |A v|_2$ and $|A|_{\max} = \max_{i,j}|a_{ij}|$, respectively. Define the column-sum and the row-sum by $\|A\|_1 = \displaystyle \max_{1\leqslant j\leqslant m} \displaystyle \sum^{n}_{i=1} |a_{ij}|$ and $\|A\|_{\infty} = \displaystyle \max_{1\leqslant i\leqslant n} \displaystyle \sum^{m}_{j=1} |a_{ij}|$. We write $a_n=O(b_n)$ or $a_n\lesssim b_n$ if there exists a positive constant $C$ such that $a_n/b_n\leqslant C$ for all large $n$, and denote $a_n=o(b_n)$ (resp. $a_n\sim b_n$), if $a_n/b_n\rightarrow 0$ (resp. $a_n/b_n \rightarrow c$ for a positive constant $c$). For two sequences of random variables $(X_n)$ and $(Y_n)$, we write $X_n=o_p( Y_n)$ if $X_n/Y_n \rightarrow 0$ in probability. Let $I_{N}$ be an $N \times N$ identity matrix, $\mathbf{I}(\cdot)$ the indicator function, $\mathbf{1}_{N}$ a $N \times 1$ vector with each element as one, and $\mathbb{N}$ the integer set.

The Model

Consider the large scale network time series data with $N$ nodes for $1 \leqslant i \leqslant N$, and $T$ time periods for $1\leqslant t\leqslant T$, which is observationally equivalent to a regular panel data. To describe their relationship, we construct an adjacency matrix, $A = (a_{ij}) \in \mathbb{R}^{N \times N}$, where $a_{ij} = 1$ if the node $i$ follows the node $j$, and $a_{ij} = 0$ otherwise. We do not allow the self-following relation, $a_{ii} = 0$. Define the row-normalised network matrix as $W = (w_{ij}) \in \mathbb{R}^{N \times N}$, where $w_{ij} = n_i^{-1} a_{ij}$ and $n_i = \sum_{j=1}^N a_{ij}$. Let $\mathbb{Y}_{t}=(Y_{1t}, \cdots, Y_{Nt})^{\top} \in \mathbb{R}^{N}$ be the continuous response (e.g., tweet length) collected at time $t$, and $U_{it}$ be a sequence of $i.i.d.$ uniform random variables on the set $[0,1]$.

We then consider the following DNQR model:

for $i = 1, \cdots, N$ and $t = 1, \cdots, T$, where $\gamma_{j}^0(\cdot)$ for $j=0, 1, 2, 3$, $\alpha_{l}^0$ for $l=1, \cdots, q$, and each elements in $\beta_{k}^0 \in \mathbb{R}^{m} $ for $k = 0, 1, \cdots, p$ are unknown parameter functions from $[0,1]$ to $\mathbb{R}$, and the superscript $0$ is used to denote the true value of parameters. $Z_{i} = (Z_{i1}, \cdots, Z_{iq})^{\top} \in \mathbb{R}^{q}$ is a $q \times 1$ vector of time-invariant node-specific covariates, and $F_{t} = (F_{t1}, \cdots, F_{tm})^{\top} \in \mathbb{R}^{m}$ is an $m \times 1$ vector of time-varying common covariates that capture the systematic influences on response variable, $Y_{it}$.

If the right hand side of the DNQR model (ref) is monotonically increasing in $U_{it}$, then we can write the $\tau$-th conditional quantile function of $Y_{it}$ as

where $\mathcal{F}_{t} = \{Z_{1}, \cdots, Z_{N}, \mathbb{Y}_{t-1}, \mathbb{Y}_{t}, F_{t}, F_{t-1},\cdots, F_{t-p}\}$ is the information set. The first component, $\gamma_{0}^0(\tau) + \sum_{l=1}^{q} \alpha _{l}^0(\tau ) Z_{il}$ is the quantile-specific nodal impact of the node $i$, where $\gamma_{0}^0(\tau)$ is the baseline function and $Z_{il}$s are assumed to be independent from $U_{it}s$. Next, network interactions between nodes are captured via both contemporaneous and lagged network variables, $\sum_{j=1}^{N}w_{ij}Y_{jt}$ and $\sum_{j=1}^{N}w_{ij}Y_{j,t-1}$, with $\gamma_{1}^0(\tau)$ capturing the quantile-specific simultaneous network effects and $\gamma_{2}^0(\tau)$ measuring the lagged diffusion network effects. $\gamma_{3}^0(\tau)$ is the quantile-specific momentum function, capturing the temporal dynamics for the same node. Furthermore, we control for the dynamic impacts of the (observed) common macroeconomic and financial factors, $F_{t}$, which can mitigate any remaining {common shock effect} in the data.

Let $\mathbb{F}_{t} =(F_{t}^{\top},\cdots, F_{t-p}^{\top })^{\top} \in \mathbb{R}^{(p+1)m}$. Define $\mathbf{A}_{0t} = (\gamma_{0}^0(U_{it}) + \sum_{l=1}^{q} \alpha^0_{l}(U_{it}) Z_{il}, 1\leqslant i\leqslant N)^{\top }\in \mathbb{R}^{N}$, $\mathbf{A} _{1t}=diag\{\gamma _{1}^0(U_{it}),1\leqslant i\leqslant N\} \in \mathbb{R}^{N\times N}$, $\mathbf{A}_{2t}=diag\{\gamma^0 _{2}(U_{it}),1\leqslant i\leqslant N\} \in \mathbb{R}^{N\times N}$, $\mathbf{A}_{3t}=diag\{\gamma^0 _{3}(U_{it}),1\leqslant i\leqslant N\} \in \mathbb{R}^{N \times N}$, and $\mathbf{B}_{t} = ((\beta_{0}^{0\top}(U_{it}), \cdots ,\beta _{p}^{0\top }(U_{it}))^{\top}, 1\leqslant i\leqslant N )^{\top} \in \mathbb{R}^{N\times (p+1)m}$. The DNQR model, (ref) can be expressed compactly in a matrix form:

where $\mathbf{H}_{t}=\mathbf{A}_{2t}W+\mathbf{A} _{3t}\in \mathbb{R}^{N\times N}$, $\Gamma =\mathop{\mbox{\sf E}}\left( \mathbf{A}_{0t}\right)$, and $V_{t}=\mathbf{A}_{0t}-\Gamma \in \mathbb{R}^{N}$ is $i.i.d.$ over $t$ with mean $\mathbf{0}$ and variance-covariance matrix, $\Sigma_{V} = \sigma_{V}^{2} I_{N} \in \mathbb{R}^{N \times N}$.

Notice that the DNQR model can be regarded as a substantial extension of koenker2006quantile, who provide a classic framework for the analysis of the random-coefficient model in the quantile autoregression. Moreover, the DNQR model encompasses the NAR model by zhu2017 and the NQAR model by zhu2019, through jointly incorporating contemporaneous and lagged network effects of connected nodes as well as exogenous common effects.

We show that the DNQR model is subject to the endogeneity issue due to contemporaneous network spillovers across nodes. Consider a simple two-equation system:

Assuming that $1 - a_{21}a_{12} \gamma_1^0(U_{1t}) \gamma_1^0(U_{2t}) \neq 0$, we obtain the following solutions:

As $Y_{1t}$ is a function of $U_{1t}$ and $U_{2t}$, the monotone argument cannot be applied because

where $\overline{Y}_{1t} = a_{12} Y_{2t}$ and $\overline{Y}_{2t} = a_{21} Y_{1t}$. This shows that the endogeneity is caused by the contemporaneous network term, $\overline{Y}_{it}$.

The simultaneous network spillover would cause inconsistency. Consider the simple mean regression, $Y_{it} = \lambda \sum_{j \neq i} w_{ij} Y_{jt} + \varepsilon_{it}$. Let $w_{i} = (w_{i1}, ..., w_{ij}, ..., w_{iN})^{\top} \in \mathbb{R}^{N} $, and $\tilde{w}_{ij}$ as the $(i,j)$ element of the matrix $(I_{N} - \lambda W)^{-1}$. Assuming that $\mathop{\mbox{\sf E}} (\varepsilon_{it} \varepsilon_{jt})= 0$ if $i\neq j$ and $\mathop{\mbox{\sf E}}(\varepsilon_{it}^2) = \sigma_i^{2}$, then the bias term (the average correlation between the endogenous variable and the error term) will be of the order, $\mbox{lim}_{N,T \to \infty} (NT)^{-1}\sum_i\sum_t \mathop{\mbox{\sf E}} (w_i^{\top} \mathbb{Y}_{t}\varepsilon_{it} ) =\mbox{lim}_{N \to\infty} N^{-1} \sum_i\sum_{j\neq i} w_{ij} \tilde{w}_{ji} \sigma_i^2 \approx c$, where $c$ is a constant. This is not equal to zero unless $\mbox{lim}_{N \to \infty} N^{-1}\sum_i\sum_{j\neq i} w_{ij} \tilde{w}_{ji} \sigma_i^2 = o(1)$. Thus, the estimation is likely to be biased unless the link of the network is very weak. In the quantile case, the leading bias term will be of the order:

which does not tend to zero, where $\varepsilon_{it,\tau}=Y_{it}- Q_{Y_{it}}(\tau|\mathbb{Y}_{t})$ is the $\tau$-th QR error.

The nontrivial estimation issue for the DNQR model lies in that the endogeneity caused by contemporaneous network spillovers renders the ordinary QR estimator to be inconsistent. Chernozhukov2005,Chernozhukov2006, Chernozhukov2008 propose the IVQR approach to estimating quantile treatment effects and develop the robust inference. chernozhukov2020network develop a novel technique to constructing simultaneous confidence bands for quantile functions and quantile effects in nonlinear network panels.

We follow the IVQR approach to cope with the simultaneous network endogeneity. Notice that Su2011 apply the IVQR approach to an analysis of the cross-section data using a linear spatial autoregressive model. However, our work can be regarded as a nontrivial extension of Su2011 to a dynamic network quantile model with nodal heterogeneity and common factors, which can shed further lights on uncovering the complex tail dependency in dynamic networks with a large number of nodes and time periods. More importantly, we derive a general asymptotic theory by using the spatial NED property of the network process in Section (ref).

Stationarity

In this section we derive the stationarity conditions for $\mathbb{Y}_{t}$ in (ref), and its asymptotic distribution. Notice that the DNQR model can easily produce predictions of quantiles, $\widehat{Q}_{Y_{it}}(\tau|\mathcal{F}_{t})$, given the network structure and the data history, by plugging the estimated parameters into (ref). To this end, it is important to derive the conditions under which the network process is stationary. Further, stationarity may be required to identify some parameters. For example, if we wish to uncover the variance structure of series of interest, it would be crucial to check whether it changes over time or not.

Define $S_{t} = I_N -\mathbf{A}_{1t}W$. Then, we make the following assumptions:

Assumption (ref)(1) assures the invertibility of $S_{t}$. The model (ref) has a unique solution if and only if every principal minor of $I_N -\mathbf{A}_{1t}W$ is positive, which is met by Assumption (ref)(1), though it is only a sufficient condition. Assumption (ref)(2) is necessary to obtain the strict stationarity of $\{\mathbb{Y}_{t}\}_t$. Under Assumptions (ref)(2) and (3), the covariance stationarity can be achieved. Then, we have the following lemma.

We introduce the NED concept in Section (ref) to ensure that the dependency of the processes is decaying appropriately, which is the key in proving the consistency and asymptotic normality of the proposed IVQR estimator. In sum, stationarity is required for moment estimation and forecasting while the NED property is utilized to prove the parameter consistency and asymptotic normality.

Once $\mathbb{Y}_{t}$ is shown to be strictly stationary, $\mathbb{Y}_{t}$ is covariance stationary if $\mathop{\mbox{\sf Var}}(\mathbb{Y}_{t})$ and $\Gamma_{l} = \mathop{\hbox{Cov}}(\mathbb{Y}_{t},\mathbb{Y}_{t-l})$ exist. Rewrite the model (ref) as

where $S_{t}=I_{N}-\mathbf{A}_{1t}W$. Then, we have the following covariance stationary solution:

where $\mathbb{D} _{t}=S_{t}^{-1}(\mathbf{B}_{t}\mathbb{F}_{t}+\mathbf{A}_{0t})$, $M_{t}=S_{t}^{-1}\mathbf{H}_{t}$ and $\Pi_{l}=M_{t}\times \cdots \times M_{t-l+1}$ for $l>1$ with $\Pi_{0}=I_N$ and $\Pi_{1}=M_{t}$. In the Online Appendix (ref) we prove that the covariance of $\mathbb{Y}_{t}$ exist under Assumption (ref).

Asymptotic Stationary Distribution

Define any vector $a\in \mathbb{R}^{N}$ with $|a|_{2}=1$ and fixed $d$ number of non zero elements. Let $\tilde{\mathbb{Y}}_{t} =\mathbb{Y}_{t}-\mu _{\mathbb{Y}}$, $L_{T}=\sum_{t=1}^{T}a^{\top}\tilde{\mathbb{Y}}_{t}$, and $L_{t}=L_{\lfloor t \rfloor } + (t - \lfloor t \rfloor)a^{\top} \tilde{\mathbb{Y}}_{\lfloor t \rfloor + 1}, t \geqslant 1$, where $\mu _{\mathbb{Y}} = \sf E(\mathbb{Y}_{t})$ and $\lfloor t \rfloor = \max \{k\in \mathbb{Z}:k\geqslant t\}$ is the floor function. We then show that the average response is asymptotically normally distributed.

Remark For $u =1$, Theorem (ref) implies:

where $\overline{\mathbb{Y}}=T^{-1}\sum_{t=1}^{T}\mathbb{Y}_{t}$. Thus, the mean of $\mathbb{Y}_{t}$ converges in law to a normal distribution.

The IVQR Estimation

We first introduce the estimation algorithms of the IVQR approach. We then discuss the underlying assumptions and develop the asymptotic theory.

IVQR Estimator

Suppose that there exists an $N \times \ell$ matrix of instrumental variables (IV), denoted $\mathbf{R}_{t}=(R_{1t}, \cdots ,R_{Nt})^{\top} \in \mathbb{R}^{N\times \ell}$, which is assumed to be independent of $U_{it}$. Then, we have the following quantile conditions:

where $\overline{Y}_{it} = \sum_{j=1}^{N} w_{ij} Y_{jt}$ and $X_{it}=\left( 1,Z_{i}^{\top}, \overline{Y}_{i,t-1}, Y_{i,t-1}, F_{t}^{\top}, \cdots, F_{t-p}^{\top}\right)^{\top}$ with $\phi^{0}(\tau) = [\gamma^{0}_{0}(\tau), \alpha_{1}^{0}(\tau), ...,\alpha_{q}^{0}(\tau), \gamma^{0}_{2}(\tau), \gamma^{0}_{3}(\tau), \beta_{0}^{0 \top}(\tau),..., \beta_{p}^{0 \top}(\tau)]^{\top} \in \mathbb{R}^{3+q+(p+1)m}$. The above conditional probability is a measurable function of $(X_{it},R_{it})$.

In general, the valid IVs should satisfy the quantile conditions in (ref), and do not lead to collinearity among $R_{it}$ and $X_{it}$. See Theorem (ref) for the asymptotic formula of the variance matrix of the IVQR estimator. The estimation efficiency will be improved by choosing $R_{it}$ appropriately. Following the literature, we may choose $R_{it}$ to be the higher network orders of lagged dependent variables such as $e_i^{\top}W^2\mathbb{Y}_{t-1}$, $[e_i^{\top}W^2\mathbb{Y}_{t-1}, e_i^{\top}W\mathbb{Y}_{t-2}]$ and so on, where $e_i$ is a vector with unity on the $i$-th element and zeros otherwise. Based on the satisfactory simulation evidence reported in Section (ref), we suggest using $[e_i^{\top}W^2\mathbb{Y}_{t-1}, e_i^{\top}W^3\mathbb{Y}_{t-1}]$ as IVs.

To solve (ref) we need to find the unknown true parameters $(\gamma^{0}_{1}(\tau ),\phi^{0\top}(\tau ))^{\top}$ such that $\boldsymbol{0}$ is a solution to the quantile estimation of $Y_{it} - \gamma^{0}_{1}(\tau ) \overline{Y}_{it}-X_{it}^{\top} \phi^{0}(\tau)$ on $(X_{it},R_{it})$:

where $\mathscr{G}$ is the class of measurable functions of $(X_{it},R_{it})$ and $\rho_{\tau }(u)=u\{\tau - \mathbf{I}(u<0)\}$ is the check function with $\mathbf{I}(\cdot)$ the indicator function. We then restrict $\mathscr{G}$ to the class of linear-in-parameter functions:

where $\Lambda $ is a compact set in $\mathbb{R}^{\ell }$. Alternatively, we may construct the transformed IVs by the least squares projection of $\overline{Y}_{it}$ on $(X_{it},R_{it})$ as in Chernozhukov2005,Chernozhukov2006. Then, we obtain the sample analogue of the objective function:

Let $\theta(\tau ) =(\gamma_{1}(\tau ),\phi^{\top}(\tau ))^{\top}$ and $\eta(\tau)=(\phi^{\top} (\tau),\lambda^{\top} (\tau))^{\top}$. The IVQR estimator, $(\widehat{\gamma}_{1}(\tau),\widehat{\phi}^{\top}(\tau),\widehat{\lambda}^{\top}(\tau))^{\top}$, obtained by minimizing (ref), is expected to converge to the true parameters, $(\gamma^{0}_{1}(\tau), \phi^{0\top}(\tau), \mathbf{0}^{\top})^{\top}$. For a given value of endogenous parameter, $\tilde{\gamma}_{1}(\tau)$, over a grid set of the interval $(-1,1)$, we first run the ordinary QR of $Y_{it}-\tilde{\gamma} _{1}(\tau )\overline{Y}_{it}$ on $ (X_{it},R _{it})$ and obtain the corresponding estimator, denoted $\widehat{\eta}(\tilde{\gamma}_{1}(\tau),\tau) = \left[ \widehat{\phi}^{\top}(\tilde{\gamma}_{1}(\tau),\tau), \widehat{\lambda}^{\top} (\tilde{\gamma}_{1}(\tau),\tau)\right]^{\top}$. Next, we select $\tilde{\gamma}_{1}(\tau)$ which minimizes $|\widehat{\lambda} (\tilde{\gamma}_{1}(\tau),\tau) |_2^2$ over the interval $(-1,1)$, denoted as $\widehat{\gamma}_{1}(\tau)$. The IVQR estimator of $\theta (\tau)$ is then obtained by $\widehat{\theta}(\tau) = (\widehat{\gamma}_{1}(\tau), \widehat{\phi}^{\top}(\widehat{\gamma}_{1}(\tau),\tau))^{\top}$. For a given quantile index $\tau $, the IVQR estimation can proceed as follows:

Step 1. For a given value of $\tilde{\gamma}_{1}(\tau)$, run the QR of $Y_{it} - \tilde{\gamma}_{1}(\tau) \overline{Y}_{it}$ on $(X_{it},R_{it})$ and obtain:

Step 2. Minimize a weighted norm of $\widehat{\lambda}(\tilde{\gamma}_{1}(\tau),\tau)$ over $\tilde{\gamma}_{1}(\tau)$ to obtain the IVQR estimator of $\gamma_{1}(\tau)$:

where $\mathrm{A}$ is some positive definite matrix. Without loss of generality we set $\mathrm{A} = I$ throughout the paper.

Step 3. Run the QR of $Y_{it}-\widehat{\gamma}_{1}(\tau) \overline{Y}_{it}$ on $X_{it}$, and obtain the estimator of $\phi(\tau)$, denoted $\widehat{\phi}(\tau) = \widehat{\phi}(\widehat{\gamma}_{1}(\tau),\tau)$. Finally, we obtain the IVQR estimator by

Asymptotic Theory

To develop the asymptotic theory for the IVQR estimator, we need to deal with some topological properties of $\mathbb{Y}_{t}$ that are spatially and temporally dependent. We follow Jenish2009,Jenish2012 and utilize NED to address the spatial dependence of the statistics. The derivation of the asymptotic property follows from the standard M-estimation, including the quantile loss function as a special case. First, conditional on the common factors, we show in Section (ref) that the elements of $\{\mathbb{Y}_{t}\}_t$ is an NED process. Then, in Section (ref), we derive the asymptotic distribution of the IVQR estimator under certain regularity conditions. As we aim to apply the DNQR model to a network dataset with the large number of nodes and time periods, we mainly employ the large $N$ and large $T$ asymptotics, though the asymptotic theory can be equally developed for {large $T$ and fixed $N$ or fixed $N$ and large $T$ (as pointed out by an anonymous referee.)}

NED Properties of the Network Processes

Jenish2009,Jenish2012 extend the notion of NED processes used in the time series to random fields. This class of NED processes can accommodate a wide range of models with a spatial dependence. They derive the central limit theorem and the law of large numbers for NED random fields. We first review the definition and some properties of NED random fields in Jenish2012.

The observations for each node can be modeled as a realization of a dependent heterogeneous process indexed by a point in $\mathbb{R}^{d} $ with $d \geqslant 1$. We consider a random field $D \subseteq \mathbb{R}^{d}$. The space $\mathbb{R}^{d}$ is endowed with the metric $\rho(j,j') = \max_{1 \leqslant l \leqslant d} |j_{l}-j'_{l}|$ with the corresponding norm, $|j| = \max_{1 \leqslant l \leqslant d}|j_{l}|$, where $j_{l}$ is the $l$-th element of $j$. The distance between any subsets $U,V\subseteq D$ is defined as $\rho (U,V) = \inf \{\rho (j,j'): j\in U \text{ and } j' \in V\}$. Let $|U|$ denote the cardinality of a finite subset, $U$. In the two dimensional case with $d =2 $ and $j = j(i,t)$, we have: $\rho((i,t), (i^{\prime}, t^{\prime })) = \max(|i-i^{\prime}|, |t-t^{\prime}|)$.

The minimum distance assumption in Assumption (ref) is used for increasing domain asymptotics. It ensures the growth of the sample size as the sample regions $D_{NT}$ expands. The setting is introduced in Jenish2012 for spatial mixing and NED processes. Note that the space $D$ can be a space of socio-economomic characteristics or geographical space, and the metric is not restricted to physical distance.

Next, we present the $L_{2}$-NED properties of random field $\mathcal {Z}$ on some $\alpha$-mixing random field $\zeta$. The $\alpha$-mixing coefficient is defined below.

Unlike the standard mixing time-series processes, the mixing coefficients for random fields depend not only on the distance between two sets but also on their sizes. We further assume that $\overline{\alpha }(u,v,h)\leqslant \varphi (u,v)\widehat{\alpha }(h)$, where the function $\varphi (u,v)$ is nondecreasing with $u$ and $v$, and $\widehat{\alpha }(h)\rightarrow 0$ as $h\rightarrow \infty $. This implies the two different sources of dependence: (i) the decay of dependence with the distance, and (ii) the accumulation of dependence as the sample region expands. In the random field literature, $\varphi (u,v)=(u+v)^{a},a\geqslant 0$ or $\varphi (u,v)=\min (u,v)$, can be commonly selected, see Jenish2012.

Define $\mathcal{C}_f \stackrel{\mathrm{def}}{=} \sigma ( \mathbb{F}_{t}, \cdots, \mathbb{F}_{t-p}) $, $\mathcal{C}_z \stackrel{\mathrm{def}}{=} \sigma ( Z_{1}, \cdots, Z_{N} )$, where $\sigma(\cdot)$ denotes a sigma field, and $\mathcal{C} = \mathcal{C}_f \cup \mathcal{C}_z $. We now discuss the NED properties of $\{Y_{it}\}_{i,t}$ on the base $\{U_{it}\}_{i,t}$; $\mathscr{F}_{it}(s) = \sigma(U_{i't'}, \mathcal{C}: (i',t')\in D_{NT}, \rho((i',t'),(i,t) )\leqslant s)$ is the $\sigma$-field generated by random vectors, $ U_{i't'}$ located within distance $s$ from $(i,t)$.

Notice that the innovation $U_{it}$ is assumed to be $i.i.d.$ over $i$ and $t$ in Assumption (ref)(1), though it is well-known that $i.i.d.$ is a special case of $\alpha$-mixing. The above condition implies that $\{U_{it}\}_{i,t}$ is an $\alpha$-mixing random field with an $\alpha$-mixing coefficient, $\overline{\alpha}(u,v,h) \leqslant (u+v)^{\varsigma} \widehat{\alpha}(h)$, $\varsigma \geqslant 0$, where $\widehat{\alpha}(h)$ satisfies $\displaystyle \sum_{h=1}^{\infty} h^{d(\varsigma_{0}+1) - 1} \widehat{\alpha}^{\delta /(4+2\delta)}(h) <\infty$ with $\varsigma_{0}=\delta \varsigma /(2+\delta)$ and some constant $\delta >0$.

Following xu2015maximum, we outline some conditions on NED properties of $\mathbb{Y}_{t}$.

Assumption (ref) is mainly used to restrict the NED coefficients, $\varphi (s)\rightarrow 0$ as $s\rightarrow \infty $. Assumption (ref) (1) allows two individuals to have direct interaction even though their locations are far away from each other, with the requirement that the strength of interaction $w_{ij}$ declines with the distance $\rho (i,j)$ in the power of $c_{w}$. This assumption is in line with xu2015maximum. By excluding a limited number of nodes $K(\geqslant 1)$, the total effects on other units from each node should be bounded, i.e., we assume that ${\sup }|\gamma_{1}^0(u)|\sup_{j}\sum_{i=1}^{N}|w_{ij}|<\Upsilon $ or $\mbox{sup}|\gamma_{1}^0(u)| \sup_{j}\sum_{i=1}^{N}|w_{ij}|<1$ w.l.o.g. This corresponds to the existence of a narrow number of units with large aggregate effects on others even as the total number of nodes rises. Assumption (ref)(2) allows two individuals to have direct interaction only if they are located within a specific distance. Notice that this assumption does not allow star nodes, but one can see pesaran2020econometric,pesaran2021estimation and kapetanios2021detection for some extensions.

Let $u_{it} = u_{it}(\gamma_{1},\phi ,\lambda ,\tau ) = Y_{it} -\gamma_{1}(\tau) \overline{Y}_{it} - X_{it}^{\top}\phi(\tau) -{R}_{it}^{\top}\lambda (\tau)$ with the check function, $\rho_{\tau}(u)=(\tau -\mathbf{I}(u \leqslant 0))u$ and $\psi_{\tau}(u) =\tau -\mathbf{I}(u \leqslant 0)$ (the directional derivative of $\rho_{\tau}(u)$). Proposition (ref) provides the NED properties of $\{Y_{it}\}_{i,t}$, and its transformations $\{\rho_{\tau}(u_{it})\}_{i,t}$, $\{\psi_{\tau}(u_{it})\}_{i,t}$ on the base $\{U_{it}\}_{i,t}$.

Define $s_{it}(\gamma _{1}^{0},\eta ^{0}(\gamma _{1}^{0},\tau ),\tau )=\psi _{\tau }\left\{ Y_{it}-\gamma _{1}^{0}(\tau )\overline{Y}_{it}-\Psi _{it}^{\top }\eta ^{0}(\gamma _{1}^{0},\tau )\right\} \Psi _{it}$ where $\Psi_{it} = (R^{\top}_{it}, X^{\top}_{it})^{\top}$, $\check{s}_{it}= \check{s}_{it}(\gamma _{1}^{0},\eta ^{0}(\gamma _{1}^{0},\tau ),\tau )=s_{it}(\gamma_{1}^{0},\eta ^{0}(\gamma_{1}^{0},\tau),\tau ) - \mathop{\mbox{\sf E}}s_{it}(\gamma_{1}^{0},\eta ^{0}(\gamma _{1}^{0},\tau ),\tau )$ and $\mathop{\mbox{\sf E}}s_{it}(\gamma_{1}^{0},\eta ^{0}(\gamma _{1}^{0},\tau ),\tau ) = 0$ . Conditioning on $\mathcal{C}$, it is easily seen that the process $\{\check{s}_{it}\}_{i,t}$ is NED. Define an array of positive constants by $\{c_{it},(i,t) \in D_{NT},NT\geqslant 1\}$. To derive the central limit theorem for $G_{NT}^{0} = \frac{1}{\sqrt{NT}} \sum_{i=1}^{N}\sum_{t=1}^{T}\check{s}_{it} = \frac{1}{\sqrt{NT}}\sum_{i=1}^{N}\sum_{t=1}^{T}[s_{it}(\gamma _{1}^{0},\eta^{0}(\gamma _{1}^{0},\tau ),\tau )-\mathop{\mbox{\sf E}}s_{it}(\gamma_{1}^{0},\eta ^{0}(\gamma _{1}^{0},\tau ),\tau )]$, where the variance of $G_{NT}^{0}$ is given by $\displaystyle \Omega_{0} = \tau(1- \tau) \lim_{N, T\to\infty} (NT)^{-1}\sum_{i,t}\mathop{\mbox{\sf E}} (\Psi_{it} \Psi_{it}^{\top}|\mathcal{C}) $, we make the following assumptions:

Assumption (ref) imposes the moment conditions of $\rho_{\tau}(u_{it})$ while Assumption (ref) sets the mixing coefficients of the underlying mixing fields. Since $U_{it}$ is assumed to be $i.i.d.$, these conditions are automatically satisfied.

Asymptotic Distribution of the IVQR Estimator

The compactness of the parameter space in Assumption (ref)(1) is needed for $\gamma_{1}(\tau)$ due to the non-convex objective function. Assumption (ref)(2) implies the global identification of the parameters while the continuity condition is required for deriving the asymptotic normality. Assumption (ref)(3) requires that $\theta^{0}(\tau) = (\gamma_1^{0}(\tau), \phi^{0\top}(\tau))^{\top}$ to be the unique solution to ${S}_{\infty}(\theta, \tau) =0$, which is necessary for consistency of the estimator.

Let $\widehat{\theta}(\tau) = (\widehat{\gamma}_{1}(\tau),\widehat{\phi}^{\top}(\tau))^{\top}$ be the IVQR estimator of $\theta^0 (\tau) =(\gamma_{1}^0(\tau ),\phi^{0\top}(\tau ))^{\top}$, where $\widehat{\phi}(\tau) = \widehat{\phi}(\widehat{\gamma}_{1}(\tau),\tau)$. Define the $(4+q+(p+1)m)\times (4+q+(p+1)m)$ matrices:

Under Assumption (ref), the NED process, $\{\check{s}_{it}\}_{i,t}$ satisfies the central limit theorem: $G_{NT}^{0}(\theta^{0},\tau) = \frac{1}{\sqrt{NT}} \sum_{i=1}^{N}\sum_{t=1}^{T} \check{s}_{it}$ follows a zero mean Gaussian process with covariance function, $\displaystyle \Omega_{0} = \tau(1- \tau) \lim_{N, T\to\infty} (NT)^{-1}\sum_{i,t}\mathop{\mbox{\sf E}} (\Psi_{it} \Psi_{it}^{\top}) $.

We estimate $\Omega_{0}^*$ and $J^{*}(\tau)$ consistently by

where $\widehat{u}_{it} = u_{it}(\widehat{\theta}(\tau), \tau )$, and $h_{b}$ is the bandwidth (see (ref) below for details). Notice that the use of valid IVs do not result in $\widehat{J}^{*}(\tau)$ and $\widehat\Omega_{0}^*$ having or being closer to singularities, which causes the variance of the estimator $\widehat{\theta}(\tau)$ to be unreliably large.

Monte Carlo Simulations

In this section, we examine finite sample properties of the IVQR estimator via a Monte Carlo simulation study using three different network structures.

The Setup

We construct the data generating process based on the DNQR model as follows: First, we generate the five nodal covariates, $Z_{i}=(Z_{i1},\cdots ,Z_{i5})^{\top} \in \mathbb{R}^{5}, (q=5)$ from a multivariate normal distribution $\mathop{\mbox{\sf N}} (\mathbf{0},\Sigma _{z})$, where $\Sigma _{z}=(\sigma _{j_{1}j_{2}})$ and $\sigma _{j_{1}j_{2}}=0.5^{|j_{1}-j_{2}|}$. Then, we construct the two common covariates, $F_{t}=(F_{1t},F_{2t})^{\top }\in \mathbb{R}^{2} , (m=2)$ from the i.i.d standard normal distribution. Let the true parameters $\gamma ^{0} _{j,it}=\gamma ^{0}_{j}(U_{it})$ for $ j=0,1,2,3$, $\alpha^{0} _{j,it}=\alpha^{0} _{j}(U_{it})$ for $j=1,...,5$, and $\beta^{0} _{jk,it}=\beta^{0} _{jk}(U_{it})$ for $j=1,2$ and $k=0,1,$ where we set the lag of two common covariates to 1 ($p=1)$. We then generate the random coefficients by

where $\Phi (\cdot )$ is the standard normal distribution function, $\mathbf{G}(\cdot ,a,b)$ is the Gamma distribution function with shape parameter $a$ and scale parameter $b$, and $\mathfrak{u}_{it}$s are i.i.d random variables, generated either from (a) the standard normal distribution or from (b) the $t$-distribution with $5$ degrees of freedom. Notice that $U_{it}$ can be generated by $U_{it}=F(\mathfrak{u}_{it})$, where $F(\cdot)$ is cumulative distribution function of $\mathfrak{u}_{it}$. Finally, $\mathbb{Y}_{t}$s are generated by (ref).

To check the robustness of the finite sample performance of the IVQR estimator, we consider the following three different adjacency matrices, e.g. zhu2019.

Type 1. (Dyad Independence Model) holland1981exponential introduce this model with a dyad, $D_{ij}=(a_{ij},a_{ji})$ for $1\leqslant i<j\leqslant N$, where $D_{ij}$s are assumed to be independent. We set the probability of dyad being mutually connected to $\operatorname{P}(D_{ij}=(1,1))=2N^{-1}$ to ensure the network sparsity. Then, we set $\operatorname{P}(D_{ij}=(1,0))=\operatorname{P}(D_{ij}=(0,1))=0.5N^{-0.8}$, which implies that the expected degree for each node is $O(N^{0.2})$. Accordingly, we have $\operatorname{P}(D_{ij}=(0,0))=1-2N^{-1}-N^{-0.8}$, which tends to 1 as $N \rightarrow \infty $.

Type 2. (Stochastic Block Model) We first consider the Stochastic Block Model with an important application in community detection by zhao2012consistency. We follow nowicki2001estimation and randomly assign each node a block label index from 1 to $L$, where $L\in \{5,10,20\}$. We then set $ \operatorname{P}(a_{ij}=1)=0.3N^{-0.3}$ if $i$ and $j$ are in the same block, and $\operatorname{P}(a_{ij}=1)=0.3N^{-1}$ otherwise. Thus, the nodes within the same block have higher probability of connecting with each other than the nodes between blocks.

Type 3. (Power-law Distribution Network) In practice, the majority of nodes in the network have a small number of links while a small number of nodes have a large number of links, see barabasi1999emergence. In this case the degrees of nodes can be characterized by the power-law distribution. We simulate the adjacency matrix as follows: For each node, we generate the in-degree, $d_{i}=\sum_{j}a_{ji}$ according to the discrete power-law distribution such as $\operatorname{P}(d_{i}=\check{k})=c\check{k}^{-\beta }$, where $c$ is a normalizing constant and the exponent parameter $\beta $ is set at $2.5$ as in clauset2009power. Finally, for the $i$-th node, we randomly select $d_{i}$ nodes as its followers.

To estimate the variance of the IVQR estimator, we follow koenker2006quantile and select the bandwidth $h_{b}$ in (ref) as follows (hall1988distribution):

where $\bar{\varphi}(\cdot)$ and $\Phi(\cdot)$ are the probability density and distribution function of standard normal distribution and $\varrho_{\alpha}$ satisfies $\Phi(\varrho_{\alpha})=1-\alpha/2$ for the construction of $1-\alpha$ confidence intervals. We have also considered an alternative selection criterion by bofinger1975estimation, and obtained qualitatively similar results that are available upon request.

For the IVQR estimation, we suggest using $\mathbf{R}_{t}=[W^{2}\mathbb{Y}_{t-1},W^{3}\mathbb{Y}_{t-1}]$, where $\mathbf{R}_{t}=(R_{1t},\cdots ,R_{Nt})^{\top}$ and $W$ is the row-sum normalized network matrix. Although we may select the higher network orders such as $[W^{2} \mathbb{Y}_{t-1},W^{3}\mathbb{Y}_{t-1},W^{2}\mathbb{Y}_{t-2},W^{3}\mathbb{Y}_{t-2},...]$, we find that these two instruments are often the best choice.

Simulation Results

Using $1000$ replications, we evaluate the finite sample performance of the IVQR estimator at the different quantiles, $\tau = 0.1,0.5,0.9$ for the $(N,T)$ pairs with $N,T = 100,200,500$. Table (ref) presents the simulation results for Type 1 network in terms of RMSE. Overall, RMSEs of all the parameters decrease monotonically as $N$ or $T$ increases across the different quantiles and for the different distributions of $\mathfrak{u}_{it}$, which is in line with the asymptotic theory. But, RMSEs of $\gamma_{1}$ are larger than those of other parameters, especially in small samples ($N,T =100$), which mainly reflects uncertainty associated with the selection of the IVs. As the sample size grows ($N=500$ or $T=500$), all RMSEs decline sharply. Turning to the simulation results for Type 2 and Type 3 networks, reported in (ref) and (ref) in the Online Appendix, we observe qualitatively similar findings. In the Online Appendix, we also provide the simulation results for biases, see Tables (ref)--(ref). The biases of the IVQR estimators are mostly negligible across the different quantiles, for the different distributions of $\mathfrak{u}_{it}$ and for all the sample pairs of $(N, T)$.

{[Insert Table (ref) here]}

Next, in Table (ref), we report the coverage probability for the Type 1 network by evaluating if the estimates fall into the 95% confidence interval at each replication. Overall, we find that the coverage probabilities of all the parameters are close to the nominal 95% level across the different quantiles and the different distributions of $\mathfrak{u}_{it}$, and for all the sample pairs $(N,T)$. This implies the accuracy of the inference for the IVQR estimator. From the results for Type 2 and Type 3 networks, reported in Tables (ref) and (ref) in the Online Appendix, we also find that the coverage probabilities of all parameters are close to the nominal 95% level.

[Insert Table (ref) here]

For comparison, we report the simulation results by applying the ordinary QR estimator in Tables (ref)--(ref) in the Online Appendix (ref). We observe that RMSEs are larger. More importantly, RMSEs barely decrease with the sample size, especially for $\gamma_{1}$. Furthermore, coverage probabilities are well below the nominal $95\%$ level, and the biases are large especially for $\gamma_{1}$, which remain substantial as the sample size increases. This clearly demonstrates the importance of using the IVQR estimator for the DNQR model.

Application

We explore the financial network quantile connectedness among the stock returns. Anton2014 find that stock returns tend to display significant comovements due to common active mutual fund owners. In addition, Pirinsky2006 document strong comovements in the stock returns of firms headquartered in the same geographic area. Garcia2012 point out that the firms headquartered in the same geographic area have achieved uniformly excessive returns compared to geographically dispersed firms (the return local bias).

We consider the two different financial network structures: the common shareholder based network and the headquarter location based network. Notice that the order of network nodes can be randomly pre-determined and has no effect on the estimation results. For simplicity, we set up $\mathbb{Y}_{t} \in \mathbb{R}^{N}$ in which the row elements are arranged in alphabetical order by the unique trading code. Then the order of nodes in the pre-determined networks are the same as the individual order in $\mathbb{Y}_{t}$, that is alphabetically ordered by the unique stock trading code. The pre-determined networks are constructed by using information on the common mutual fund ownership and the uniform headquarter location. In particular, we let the stocks be connected if they are invested in by at least five common shareholders ($W_{CS5}$) while the companies with headquarters located in the same city are treated as connected ($W_{HQ}$).

We collect the data on all the stocks traded in NYSE and NASDAQ in 2016 from Datastream. The dataset on mutual fund holdings are downloaded from Thomson Reuters whilst the addresses of firms' headquarters are collected from COMPUSTAT. After merging these data from the databases according to the unique trading code and removing the stocks with missing values, we finally obtain 943 stocks ($N = 943$) over the whole time period $T = 252$. We then collect these stock return data from Datastream. We also obtain node-specific covariates such as market capitalization, book value per share, cash flow and price-earning ratio from Datastream, which are then standardized. Finally, we collect VIX from Datastream and the Fama-French three factors (excess market return, SMB, HML) from the website of French's homepage as the common covariates.

The network density is 3.24% for $W_{CS5}$ and 0.63 % for $W_{HQ}$, respectively. In Figure (ref) we display the topology of two networks for the top 100 market-value stocks only for visualization convenience. The larger nodes imply the higher market capitalization while the darker nodes present the higher connectedness especially for the network with $W_{CS5}$. Here we observe quite different network structures. There is a large connected group in Figure (ref), showing that the stocks are more centrally connected by common investors. On the contrary Figure (ref) displays more small groups, implying that the stocks are more locally connected when the network is measured by uniform headquarter locations.

[Insert Figure (ref) here]

Following the simulation results, we select IVs as $\mathbf{R}_{t}=[W^{2}\mathbb{Y}_{t-1},W^{3}\mathbb{Y}_{t-1}]$. We present the estimation results for the proposed DNQR model together with the two alternative models: (i) the original NQAR model without contemporaneous network effects and common covariates, and (ii) the factor-augmented NQAR model without contemporaneous network effects, denoted NQARF. To compare the relative performance of the alternative models, we follow Koenker1999 and evaluate the goodness of fits across the different quantiles. Consider a linear model for the conditional quantile function,

where $X_{it}=(X^{\top}_{1it},X^{\top}_{2it})^{\top}$. Let $\widehat{\theta}(\tau)=(\widehat{\theta}_{1}^{\top}(\tau),\widehat{\theta}_{2}^{\top}(\tau))^{\top}$ be the unrestricted estimator, which is the minimizer of $\widehat{V}(\tau) = \min \displaystyle \sum_{i=1}^{N} \displaystyle \sum_{t=1}^{T} \rho_{\tau} \left\{Y_{it} - X^{\top}_{it} \theta \right\}$ while $\widetilde{\theta}(\tau) = (\widetilde{\theta}^{\top}_{1}(\tau), \boldsymbol{0}^{\top})^{\top} $ denotes the minimizer for the constrained model, $\widetilde{V}(\tau) = \min \displaystyle \sum_{i=1}^{N} \displaystyle \sum_{t=1}^{T} \rho_{\tau} \left\{Y_{it} - X^{\top}_{1it} \theta_{1} \right\}$. Define the goodness-of-fit criterion as

which measures the overall decreased percentage of the quantile loss function of the unrestricted model with respect to the restricted model.

The estimation results for the network $W_{CS5}$ are presented in Table (ref). For convenience we present the coefficients and the standard errors multiplied by $10^{2}$. We find that the estimated contemporaneous network effects ($\gamma_1$) by the DNQR model, are significantly positive and dominate all other effects across all quantiles. The lagged network coefficient ($\gamma_2$) and the dynamic coefficient ($\gamma_3$) are also significant across quantiles with relatively smaller magnitudes. Furthermore, the goodness of fit, $R^2(\tau)$ reported in the last row, shows that the overall loss function of the DNQR model drops about 7%--9.5% relative to the NQAR and 6.9%--9.5% relative to the NQARF, respectively, suggesting that the contemporaneous network effects should be explicitly accommodated in the dynamic network quantile model. Moreover, the effects of node-specific covariates are significant across quantiles (except Cash) while those of common covariates are all significant across quantiles.

[Insert Table (ref) here]

Next, we display the QR coefficients across the different quantiles in Figure (ref). The dashed line is the QR coefficient while the grey area indicates a kernel density based confidence band advanced by powell1991. The contemporaneous network and dynamic coefficients, $\widehat{\gamma}_1(\cdot)$ and $\widehat{\gamma}_3(\cdot)$, are significant (as their bands exclude the null effect), while the lagged diffusion network coefficient, $\widehat{\gamma}_2(\cdot)$ tends to be insignificant only at the middle. $\widehat{\gamma}_1(\cdot)$ shows a downward trend with the quantile level, suggesting that the contemporaneous quantile connectedness is stronger at the lower tails (mainly characterised with the market turmoils). On the other hand, both $\widehat{\gamma}_2(\cdot)$ and $\widehat{\gamma}_3(\cdot)$ display the $U$-shaped pattern, implying that their effects are stronger at the tails than at the median. The QR effects of node-sepcific covariates mostly display a downward trend with the quantile level (except for insignificant Cash), suggesting their effects are stronger at the lower tails than at the upper tails (mainly chracterised with the bulls market). Turning to the QR effects of common factors, we observe a mixed finding: the impacts of VIX and the market factor increase with quantiles whilst those of SMB and HML factors decrease with quantiles.

[Insert Figure (ref) here]

Finally, as the robustness check, we provide the two additional estimation results in the Online Appendix. First, we reconstruct the network matrix by changing the number of common shareholders to $CS = 3$ $(W_{CS3})$ and $CS = 7$ $(W_{CS7})$, and the results are reported in Tables (ref)--(ref) and Figures (ref)--(ref). Notice that the network density of $W_{CS3}$ is dense at 25.25% and relatively sparse at 0.41% for $W_{CS7}$. Overall, we find qualitatively similar results to those reported for $W_{CS5}$. One notable observation is that the contemporaneous network effects measured by $\gamma_1$ tend to decrease monotonically as the network becomes more sparse. For example, at $\tau = 0.1$, $\widehat{\gamma}_1$ is estimated at 0.69 for $W_{CS3}$, 0.54 for $W_{CS5}$, and 0.35 $W_{CS7}$, respectively. Still, we find that the patterns of the quantile specific coefficients reported in Figures (ref)--(ref) are qualitatively similar to those displayed in Figure (ref).

Next, we estimate the models using the headquarter location network $W_{HQ}$, and present the estimation results in Table (ref) and Figure (ref) in the Online Appendix. Again, we find the qualitatively similar results, highlighting the importance of the contemporaneous network effect, which is also stronger at the lower tails than at the upper tails.

Conclusion

We develop a dynamic network quantile model that accommodates both temporal and cross-sectional dependence. Using the predetermined network information, we analyse the dynamic quantile connectedness within a network topology. The distinguishing feature of the DNQR model lies in that the behavior/response of a given node is not only influenced by its previous behavior/response, but also connected with a weighted average of contemporaneous and lagged behaviors/responses from peers.

The main challenge associated with the DNQR model is the presence of endogeneity stemming from the simultaneous network effect. In this regard, we develop the IVQR estimation, and derive the consistency and asymptotic normality of the IVQR estimator using the NED property of the network process. Monte Carlo exercises confirm the satisfactory performance of the IVQR estimator with the predetermined internal IVs across different quantiles under the different network structures.

Finally, we demonstrate the usefulness of our proposed approach with an application to the dataset on the stocks traded in NYSE and NASDAQ in 2016. In particular, we find that the contemporary network effects are significant and dominant across all quantiles. Furthermore, their effects display a downward trend with the quantile level, suggesting that the contemporaneous quantile connectedness is stronger at the lower tails.

Acknowledgements

We are mostly grateful for the insightful comments by the editor, the associated editor and three anonymous referees. Xu acknowledges partial financial support of the Natural Science Foundation of China (Grant No. 71803140). Shin and Wang gratefully acknowledge partial financial support from the ESRC (Grant Reference: ES/T01573X/1). Zheng gratefully acknowledges partial financial support from the Royal Economic Society. The usual disclaimer applies.

\numberwithin{equation}{section}

Section (ref) provides the proofs for Theorems (ref)--(ref). Section (ref) presents the additional simulation results on the performance of the IVQR and the ordinary QR estimators under the different network structures. Section (ref) reports the additional empirical results by employing the alternative common shareholder network structures constructed by imposing the different number of common shareholders, and the uniform headquarter location network.