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.
136,348 characters · 24 sections · 39 citation commands
Estimation and Inference for Latent Dual Networks Using High-Dimensional IV Screening
Networks and strategic interactions play an important role in shaping economic outcomes across a wide range of settings, including innovation, finance, trade, and regional development. A fundamental feature of many economic environments is that interactions may reflect multiple and potentially opposing forces. The same network may simultaneously generate reinforcing effects through imitation, learning, and information transmission, as well as displacement effects arising from competition, market-share reallocation, and strategic differentiation. Consequently, spillovers may reflect qualitatively different forms of dependence.
Motivated by this observation, we consider an econometric framework in which interactions are generated by two latent networks, one capturing the reinforcing effects and the other capturing the displacement effects
Here, $y_{i,t}$ is the outcome of unit $i$, $\boldsymbol{x}_{i,t}$ is a vector of observed covariates, $\eta_i$ is a unit-specific effect, and $\varepsilon_{i,t}$ is an idiosyncratic disturbance. A key departure from conventional network and spatial models is that Eq. (ref) allows interactions to arise through two latent channels. These channels are represented by non-negative interaction/weighting matrices $\boldsymbol{W}_1=(w_{1,i,j})$ and $\boldsymbol{W}_2=(w_{2,i,j})$, with associated intensities $\rho_1$ and $\rho_2$. Together, they give rise to a dual-network framework.
Given that the two networks are latent, for the sake of identification and interpretation we assume $\rho_1>0$ and $\rho_2<0$,\footnote{Or, the observationally equivalent setting $\rho_1<0$ and $\rho_2>0$.} so that the two channels operate in opposite directions. The first channel captures positive spillovers, strategic complementarities, or clustering effects, whereas the second captures strategic substitution, competition, or spatial dispersion effects. Every nonzero bilateral interaction is assumed to belong to at most one of the reinforcing and displacement networks; a dyad may belong to neither. This restriction has a natural economic interpretation and is important for identification.\footnote{Without this restriction, the data would generally identify only the combined effect $\rho_1 w_{1,i,j}+\rho_2 w_{2,i,j}$ rather than the separate contributions of the two interaction channels.}
Models based on a single interaction network arise as special cases of Eq. (ref). In particular, when one interaction component is absent, or when $\rho_1$ and $\rho_2$ have the same sign, the model reduces to a conventional interaction structure in which all links generate effects of the same qualitative type. Such models are related to the literature on social interactions and the reflection problem, beginning with Manski1993, to identification through network structure in BramoulleEtAl2009, and to micro-founded linear social interaction models such as BlumeEtAl2015 and LewbelQuTang2023. They are also closely connected to spatial econometric models, see e.g. Anselin1988, Elhorst2003, ArbiaBaltagi2009, Baltagi2021Spatial, and Elhorst2014Book.
The majority of the existing literature assumes that the interaction matrix is known a priori and focuses on estimation and inference for the structural parameters conditional on the network structure, see e.g. KelejianPrucha2010,BaltagiEtAl2013,ShiLee2017,Yang2018,CuiEtAl2023,Elhorst2024TwoLaws.\footnote{Related contributions include HuangEtAl2020, who develop a two-mode network autoregressive model that allows interactions across different classes of nodes through observed network matrices; and ChenEtAl2025, who allow for heterogeneous interaction effects while maintaining a known interaction structure.} More recently, a growing literature has considered estimating sparse network links directly from the data; see, among others, AhrensBhattacharjee2015, Manresa2016, Rose2017, LamSouza2020JBES, dePaula2020, dePaulaRasulSouza2024, and KrisztinPiribauer2023. This literature builds upon available results on high-dimensional sparse estimation, where techniques such as the Lasso, the adaptive Lasso, and the elastic net regularisation are used to recover sparse dependence structures in large systems. These methods typically impose a common sign structure on network effects, effectively treating all links as manifestations of a single interaction mechanism.\footnote{Although penalised network-recovery procedures may permit negative first-stage interaction estimates, these are typically treated within a single-network framework rather than as a distinct interaction mechanism. For example, in the adaptive elastic-net procedure of dePaulaRasulSouza2024, the first-stage step imposes only $|w_{ij}|\le 1$, while the adaptive step sets $\widetilde w_{i,j}=0.05$ whenever $\widetilde w_{i,j}<0.05$ when constructing the adaptive penalty weights. Thus negative first-stage estimates are treated as small positive values for penalisation purposes, rather than as evidence of a separate displacement network.} Identification is then achieved through the reduced-form representation of the model.
The present paper is related to this literature, but differs in two important respects. First, we allow for two latent interaction channels with opposite signs, while treating both \(\boldsymbol W_1\) and \(\boldsymbol W_2\) as unknown. Therefore, the first object of estimation is the composite structural interaction matrix \[ \boldsymbol A = \rho_1\boldsymbol W_1+\rho_2\boldsymbol W_2 = (\alpha_{i,j}), \qquad \alpha_{i,j} = \rho_1w_{1,i,j}+\rho_2w_{2,i,j}. \]
Second, the proposed approach differs from existing data-driven network-recovery methods in the way $\boldsymbol A$ is recovered. Reduced-form approaches identify the network by first recovering the equilibrium mapping from covariates or shocks to outcomes and then using this mapping to infer the underlying interaction structure. Such approaches are necessarily global: the reduced form is an all-or-nothing object that depends on the joint behaviour of the entire system. As a result, conditions ensuring the existence, uniqueness, and stability of the global reduced-form representation play a central role. By contrast, the estimation approach developed here is constructive and equation-specific. For each equation $i$, it asks which candidate outcomes $y_{j,t}$ enter the structural equation for $y_{i,t}$; that is, which coefficients $\alpha_{i,j}$ are nonzero. Estimation is therefore based on local incremental contributions within each equation, using instrumental-variable screening statistics and multiple-testing control.
A central feature of this row-wise recovery strategy is that the nonzero entries of $\boldsymbol A$ may be positive or negative. The procedure does not impose a common sign restriction on the interaction matrix. This allows the recovered structural interaction matrix to contain both reinforcing and displacement forces. To this end, we develop a step-wise instrumental-variables (IV) screening procedure with multiple-testing control for recovering sparse latent interaction structures in panel network models. The procedure builds on the least-squares boosting framework of Kapetanios2026BMT, but adapts it to structural network recovery with endogenous contemporaneous outcomes. For each equation, every other unit is initially treated as a candidate link. The candidates are evaluated using IV screening statistics. At each step, the procedure selects at most one link: the candidate with the largest statistic, provided that it exceeds the multiple-testing threshold. The screening statistics are then recomputed conditional on the links selected in previous steps. The algorithm stops when no remaining candidate exceeds the threshold. We refer to the resulting procedure as Boosting One-Link-at-a-Time with Multiple Testing (BOLMT).
After $\boldsymbol A$ has been estimated, positive and negative entries are assigned to distinct latent channels. Together with a disjoint-support restriction and row normalisations, this yields a decomposition of the recovered composite matrix into reinforcing and displacement networks. Hence, the dual-network structure is not imposed through a pre-specified interaction matrix; it is obtained as a structured decomposition of the estimated structural interaction matrix.
The structural network recovery problem we consider raises several challenges that do not arise in standard model selection problems, e.g. as in Kapetanios2026BMT. First, the candidate regressors are endogenous, such that the algorithm must rely on fitted values obtained from the corresponding first-stage projections given some instruments. Second, network recovery requires solving \(N\) high-dimensional selection problems, one for each equation of the system. Finally, while in standard high-dimensional regression, the inclusion of a small number of irrelevant variables need not be particularly harmful (provided that the selected model delivers a good approximation to the underlying signal), in the present context the support itself has a structural interpretation. Consequently, distinguishing direct links from indirect or spurious links that arise through common sources of dependence becomes important for recovery of the underlying interaction structure and its economic interpretation. In this respect, step-wise conditioning is particularly valuable in the present framework because the explanatory power of proxy links is expected to diminish once the relevant direct links have been selected.
To the best of our knowledge, this is the first framework to recover latent reinforcing and displacement interaction networks and to establish exact support recovery together with post-selection IV inference in a high-dimensional panel setting.
The remainder of the paper is organised as follows. Section (ref) introduces the dual-network framework, presents the BOLMT procedure, and discusses post-selection estimation and network decomposition. Section (ref) establishes exact network recovery, recovery of the dual-network structure, and oracle-equivalent post-selection inference. Section (ref) discusses extensions, including block screening statistics and latent factor structures, as well as implementation issues. Section (ref) reports Monte Carlo evidence and Section (ref) provides an empirical illustration based on corporate financial policies. A final section concludes.
We rewrite the dual-network model in Eq. (ref) in the equivalent form
where we define \[ \alpha_{i,j} := \rho_1 w_{1,i,j} + \rho_2 w_{2,i,j}. \] The matrix of the corresponding interaction coefficients, $\boldsymbol{A}=(\alpha_{i,j})$, may be interpreted as a directed weighted network. There is an edge $j\rightarrow i$ whenever $\alpha_{i,j}\neq0$, since the outcome of unit $j$ enters the equation for unit $i$. The network need not be symmetric: the presence of $j\rightarrow i$ does not imply the presence of $i\rightarrow j$.
The representation in Eq. (ref) highlights two important features of the recovery problem with latent $\alpha_{i,j}$. First, the interaction intensities and the corresponding network weights are not separately identified from the composite coefficients $\alpha_{i,j}$. The decomposition into reinforcing $(\rho_1,\boldsymbol W_1)$ and displacement $(\rho_2,\boldsymbol W_2)$ networks is therefore a second-stage operation performed only after recovery and estimation of the structural interaction matrix $\boldsymbol A=(\alpha_{i,j})$. Second, the interaction structure is inherently heterogeneous across equations. The supports need not be symmetric, may differ in size, and may involve coefficients of different magnitudes. As a result, both the support and the coefficients of each row of $\boldsymbol A$ are individual-specific. The support-recovery problem therefore consists of $N$ separate high-dimensional selection problems, one for each row of $\boldsymbol A$.
For each unit $i$, denote by $\mathcal S_i^n$ the set of units whose outcomes enter the structural equation of unit $i$ through nonzero interaction coefficients
This support set is unobserved. The researcher observes only the set of potential candidate links, $\{1,\ldots,N\}\setminus\{i\}$, and seeks to determine which of these candidates belong to $\mathcal S_i^n$. The framework developed in this paper builds upon the decomposition
where the three (non-overlapping) sets are defined as:
The distinction between the three types of units is central for network recovery. In standard high-dimensional regression, the inclusion of a small number of proxy variables may still yield a useful approximating model. In the present context, however, the support itself has a structural interpretation. Selecting proxy links or irrelevant units may distort the estimated interaction structure and, consequently, the implied network propagation mechanism. Moreover, when the reduced form is well defined, reduced-form dependence may still be misleading about the sign and strength of the underlying structural interaction coefficient.\footnote{For example, suppose \(N=4\), and the reinforcing network contains links \(2\rightarrow 1\) and \(3\rightarrow 2\), while the displacement network contains links \(3\rightarrow 1\) and \(4\rightarrow 2\). Let \(\rho_1=0.5\) and \(\rho_2=-0.2\). Then unit 3 has a negative immediate structural interaction coefficient in the equation for unit 1 through the displacement link \(3\rightarrow 1\). At the same time, unit 3 also affects unit 1 indirectly through the path \(3\rightarrow 2\rightarrow 1\), generating a positive second-order effect. The latter dominates the former, so that reduced-form dependence between units 1 and 3 has the opposite sign from the immediate structural interaction coefficient. }
Figure (ref) provides a simple illustration of these concepts for equation $i$, where unit $j_1 \in \mathcal S_i^n$, unit $j_2 \in \mathcal S_i^p$, and unit $j_3 \in \mathcal S_i^d$.
In what follows, we formally introduce the BOLMT algorithm. The contemporaneous nature of the model implies that the network-related regressors are endogenous. In particular, if unit $j$ responds contemporaneously to shocks elsewhere in the system, then $y_{j,t}$ can be correlated with $\varepsilon_{i,t}$. Consequently, the screening statistics are constructed using instrumental variables and fitted regressors obtained from the corresponding first-stage projections. The instruments may consist of exogenous covariates, predetermined variables, or other application-specific instruments satisfying the usual relevance and exogeneity conditions.
At any stage of the algorithm, the objective is to test whether a remaining candidate unit \(j\) contributes additional information about the outcome of unit \(i\), conditional on the controls and the links selected in previous stages. Since the previously selected contemporaneous neighbour outcomes are endogenous in the network system, the stagewise statistic is computed as a partial IV statistic. Thus the previously selected link outcomes and the current candidate are instrumented jointly, and the candidate is tested after partialling out the fitted selected controls.
Let \(H_{i,m-1}\subseteq \{1,\ldots,N\}\setminus\{i\}\) denote the set of links already selected before a generic stage of the algorithm, $m$. Let \(\mathbf C_{i,H_{i,m-1}}\) denote the matrix containing all always-included regressors together with the outcomes of the units in \(H_{i,m-1}\).\footnote{Depending on the application, the always-included regressors may contain observed covariates, time lags of the dependent variable, deterministic trends, seasonal effects, fixed effects, factor proxies, or other controls that are not subject to selection.}
For a remaining candidate \(j\notin H_{i,m-1}\), let \(\mathbf v_j:= \mathbf y_j\) denote the candidate outcome vector. Let \(\boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}}\) denote the full stagewise instrument matrix used for both \(\mathbf C_{i,H_{i,m-1}}\) and \(\mathbf v_j\). Thus \(\boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}}\) contains instruments for the previously selected neighbour outcomes, instruments for the current candidate, and any exogenous controls used as their own instruments.
Define the fitted selected controls and fitted candidate by \[ \widehat{\mathbf C}_{i,j,H_{i,m-1}} = \boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}} (\boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}}'\boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}})^{-1} \boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}}'\mathbf C_{i,H_{i,m-1}}, \] and \[ \overline{\mathbf v}_{i,j,H_{i,m-1}} = \boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}} (\boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}}'\boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}})^{-1} \boldsymbol{\mathcal Z}_{i,j,H_{i,m-1}}'\mathbf v_j. \] Let
The partial IV/FWL residuals entering the scalar screening statistic are
Thus \(\widehat{\mathbf v}_{i,j,H_{i,m-1}}\) is the fitted candidate residual after projecting the fitted candidate on the fitted selected controls.
The corresponding IV estimator is given by define
with the corresponding screening statistic
Here \(\widehat{\sigma}_{i,j,H_{i,m-1}}^{2}\) is the residual variance estimator from the same full stagewise IV regression. Only the values of $j$ such that:
are included in the next stage of the algorithm. The full algorithm is summarized in Algorithm (ref), where procedure is repeated for every unit \(i=1,\ldots,N\).
In practice, we implement the threshold using
where \(p\in(0,1)\) is a tuning parameter and \(c,\delta>0\) are user-specified constants, and \(\Phi\) is the standard normal CDF. Under standard tail approximations, \(c_N\asymp\sqrt{\log N}\), so that \(c_p(N,\delta)\) provides a feasible tuning rule of the same order as the theoretical threshold sequence \(c_N\). The recovery theory covers this rule only when its leading constant is sufficiently large. Interpreting Eq. (ref) as a nominal Gaussian multiple-testing critical value would additionally require asymptotic standard normality, or appropriate long-run-variance studentisation, neither of which is imposed in the selection theory below.
Since the selected supports need not be symmetric, the resulting network is directed.
If all stagewise regressors are exogenous and the instrument matrix coincides with the regressor matrix, this construction reduces to ordinary least-squares partialling out, as in standard BMT of Kapetanios2026BMT.
Once BOLMT has been applied to all cross-sectional units, the estimated support matrix is constructed from the selected supports $\{\widehat{\mathcal{S}}_i\}_{i=1}^{N}$, with elements
Conditional on the selected support, we estimate model
for every $i$ using the corresponding IV regression, where $u_{i,t}(\widehat{\mathcal S}_i)$ includes any omitted interaction terms and equals $\varepsilon_{i,t}$ on the exact-recovery event. The post-selection instrument matrix is taken to be a fixed function of the selected support; the oracle estimator uses the same rule evaluated at $\mathcal S_i^n$.
Let $\widehat{\boldsymbol{\beta}}_i$ and $\widehat{\alpha}_{i,j}$ denote the post-selection IV estimates obtained from Eq. (ref), with $\widehat\alpha_{i,j}=0$ for $j\notin\widehat{\mathcal S}_i$. Define
Using this definition, the estimated interaction coefficients $\widehat{\alpha}_{i,j}$ are decomposed according to their sign. Define
and define the row-normalised weights, for $\ell=1,2$, by \[ \widehat w_{\ell,i,j} =
\] Thus an empty signed row of $\widehat{\mathbf W}_{\ell}$ is set equal to zero and no division by zero is made.
As a final step, we aggregate all unit $i$ specific estimates $\widehat{\boldsymbol{\theta}}_i = \left( \widehat{\boldsymbol{\beta}}_i', \widehat\rho_{1,i}, \widehat\rho_{2,i} \right)'$ using the Mean Group (MG) aggregation, as Pesaran199579:
The implementation is fundamentally local and proceeds equation by equation. It does not require global graph restrictions such as symmetry, connectedness, or particular topological features of the interaction network. Such restrictions may be useful in specific applications, but they are not needed for the row-wise selection arguments developed here.
The primary object recovered by BOLMT is the structural interaction matrix $\boldsymbol A=(\alpha_{i,j})$. The coefficient $\alpha_{i,j}$ indicates whether the outcome of unit $j$ enters the structural equation for unit $i$, and with what sign and magnitude. It is therefore an immediate interaction coefficient, not a direct effect in the usual spatial-econometric sense. Direct, indirect, and total effects are defined from the relevant multiplier matrix. For example, when the full system is stable, the effect matrix for covariate $\ell$ is based on \[ (\boldsymbol I_N-\boldsymbol A)^{-1} \operatorname{diag}(\beta_{1\ell},\ldots,\beta_{N\ell}), \] with direct effects obtained from the diagonal elements and indirect effects from the off-diagonal elements. If only a retained subsystem is stable, the same logic applies to the corresponding subsystem multiplier, conditional on the excluded outcomes being treated as observed external inputs.
The dual-network representation is therefore a structured decomposition of the recovered matrix $\boldsymbol A$. Without additional restrictions, the data identify the composite coefficients $\alpha_{i,j}$, not the separate objects $\rho_1w_{1,i,j}$ and $\rho_2w_{2,i,j}$. The sign decomposition, disjoint-support restriction, and row normalisations deliver the reinforcing and displacement weights and the corresponding row-specific interaction intensities.
This section collects the assumptions used for recovery of the structural interaction matrix $\boldsymbol{A}=(\alpha_{i,j})$ and for post-selection estimation. Afterwards, based on these assumptions consistency of the proposed screening procedure, exact recovery of the latent interaction network, and oracle properties of the resulting post-selection estimators is established. Throughout, this section all asymptotic statements are understood as jointly $N,T\to\infty$ unless otherwise stated.
The assumptions in this section are stated as population restrictions on a deterministic class of low-dimensional conditioning sets. They are not restrictions on a realised random path of the algorithm. The connection with the realised BOLMT path is made in the appendix by an induction argument: with probability tending to one, the selected set remains in the class of conditioning sets defined below.
Recall that \(\mathcal S_i^n\) denotes the set of true links. For the theoretical analysis, the remaining candidates are partitioned into two groups. The set \(\mathcal S_i^p\) contains proxy links, namely candidates for which \(\alpha_{i,j}=0\) but whose outcomes remain correlated with the structural component of \(y_{i,t}\) through indirect network propagation or common exposure. The remaining candidates are termed irrelevant links. Accordingly, define \[ \mathcal K_i := \mathcal S_i^n \cup \mathcal S_i^p \qquad k_i^K:=|\mathcal K_i|, \qquad \bar K:=\sup_i k_i^K, \] and \[ \mathcal S_i^d = \{1,\ldots,N\}\setminus(\{i\}\cup\mathcal K_i). \]
Let \[ \mathfrak U_i=\{U:U\subseteq\mathcal K_i\} \] be the deterministic class of low-dimensional conditioning sets. In the theorem statements and appendix we write \(H\) for a generic element of this class and set \[ \mathcal H_i:=\mathfrak U_i. \] Also define \[ \mathcal H_i^0=\{H\in\mathcal H_i:\mathcal S_i^n\not\subseteq H\}, \qquad \mathcal H_i^1=\{H\in\mathcal H_i:\mathcal S_i^n\subseteq H\}. \] We call \(H\in\mathcal H_i^0\) a pre-completion conditioning set and \(H\in\mathcal H_i^1\) a post-completion conditioning set. Completion means that all true links have been included in \(H\); it does not necessarily mean that the screening algorithm has stopped. Thus \(\mathcal H_i^0\) contains low-dimensional conditioning sets that do not yet contain all true links, while \(\mathcal H_i^1\) contains conditioning sets that contain the true support. Since \(\bar K<\infty\), \[ |\mathcal H_i|\le 2^{\bar K} \] uniformly in \(i\). Hence the high-dimensional probability bounds below are required to be uniform over equations and candidates, but not over all subsets of the \(N-1\) candidate set.
For any \(H\in\mathcal H_i\), let \(\mathbf C_{i,H,t}\) collect the always-included controls and the outcomes of units in \(H\). For a candidate \(j\notin H\), the scalar candidate outcome is \(y_{j,t}\). Let $\boldsymbol{\mathcal Z}_{i,j,H,t}$ denote the full stagewise instrument vector used jointly for \(\mathbf C_{i,H,t}\) and \(y_{j,t}\). It contains instruments for the conditioning variables, instruments for the current candidate, and any exogenous regressors used as their own instruments.
Define \[ \mathbf Q^{\boldsymbol{\mathcal Z}}_{i,j,H} := \operatorname{E}(\boldsymbol{\mathcal Z}_{i,j,H,t}\boldsymbol{\mathcal Z}_{i,j,H,t}'). \] The population fitted controls and fitted candidate are \[ \mathbf C^*_{i,j,H,t} := \boldsymbol\Pi_{C,i,j,H}'\boldsymbol{\mathcal Z}_{i,j,H,t}, \qquad \bar v^*_{i,j,H,t} := \boldsymbol\pi_{v,i,j,H}'\boldsymbol{\mathcal Z}_{i,j,H,t}, \] where \[ \boldsymbol\Pi_{C,i,j,H} := (\mathbf Q^{\boldsymbol{\mathcal Z}}_{i,j,H})^{-1} \operatorname{E}(\boldsymbol{\mathcal Z}_{i,j,H,t}\mathbf C_{i,H,t}'), \qquad \boldsymbol\pi_{v,i,j,H} := (\mathbf Q^{\boldsymbol{\mathcal Z}}_{i,j,H})^{-1} \operatorname{E}(\boldsymbol{\mathcal Z}_{i,j,H,t}y_{j,t}). \] Let \[ \boldsymbol\Gamma_{vC,i,j,H} := \{\operatorname{E}(\mathbf C^*_{i,j,H,t}\mathbf C_{i,j,H,t}^{*\prime})\}^{-1} \operatorname{E}(\mathbf C^*_{i,j,H,t}\bar v^*_{i,j,H,t}), \] and define the population partial fitted candidate
Here, \(\bar v^*_{i,j,H,t}\) is the population first-stage fitted value of \(y_{j,t}\), while \(v^*_{i,j,H,t}\) is the residual obtained after partialling this fitted value with respect to the population fitted controls.
For $H\in\mathcal H_i$, define the omitted structural component \[ \mu_{i,H,t} = \sum_{k\in\mathcal S_i^n\setminus H}\alpha_{i,k}y_{k,t}, \qquad \mu_{i,j,H,t}:=\mu_{i,H,t}, \] where the second notation is retained for compatibility with candidate-specific score arrays. Also let \[ \sigma_{i,j,H}^2 := \operatorname*{plim}_{T\to\infty}\widehat\sigma_{i,j,H}^2 \] denote the population limit of the residual-variance estimator from the corresponding stagewise IV regression.
The tail and dependence conditions below are imposed uniformly on the centred versions of the following scalar processes, over \(i\), \(j\), and \(H\in\mathcal H_i\): \[ v^*_{i,j,H,t}\varepsilon_{i,t}, \qquad v^*_{i,j,H,t}\mu_{i,j,H,t}, \qquad (v^*_{i,j,H,t})^2, \] \[ \boldsymbol{\mathcal Z}_{i,j,H,t,a} \boldsymbol{\mathcal Z}_{i,j,H,t,b}, \qquad \boldsymbol{\mathcal Z}_{i,j,H,t,a} y_{j,t}, \qquad \boldsymbol{\mathcal Z}_{i,j,H,t,a} \mathbf C_{i,H,t,b}, \qquad \varepsilon_{i,t}^{2}. \] These are the scalar score, signal, Gram-matrix, and first-stage arrays controlled by the probability bounds.
For a scalar array \(\xi_t\), write \(\xi_t\in E(s)\), \(s>0\), if \(\sup_t\operatorname{E} \exp(a|\xi_t|^s)<\infty\) for some \(a>0\). Write \(\xi_t\in H(\vartheta)\), \(\vartheta>2\), if \(\sup_t\operatorname{E}|\xi_t|^\vartheta<\infty\). For \(\zeta>0\), define
The functions \(f_T(\cdot)\) and \(g_T(\cdot)\) are the thin-tail and heavy-tail concentration envelopes of DendramisEtal2022. All population moments used below are assumed to be time invariant, unless stated otherwise.
For each equation \(i\), let \[ U_i(H)=\mathcal S_i^n\setminus H, \qquad P_i(H)=\mathcal S_i^p\setminus H, \] denote the sets of true and proxy links not yet selected under conditioning set \(H\). We use the convention that a maximum over an empty set is zero. For any remaining candidate \(j\notin H\), set \(\widetilde y_{j,H,t}:=v^*_{i,j,H,t}\), with the equation index $i$ suppressed on the left-hand side. Define \[ d_{i,j,H} = \sigma_{i,j,H} \left\{ \operatorname{E}[(\widetilde y_{j,H,t})^2] \right\}^{1/2}. \] Define the normalised IV alignment \[ \bar\gamma_{i,j,k}(H) = \frac{ \operatorname{E}(\widetilde y_{j,H,t}y_{k,t}) }{ d_{i,j,H} }, \qquad k\in U_i(H), \] and the normalised population IV screening signal \[ m_{i,j}(H) = \frac{ \operatorname{E}(\widetilde y_{j,H,t}\mu_{i,j,H,t}) }{ d_{i,j,H} }. \] Because $v^*_{i,j,H,t}$ is a linear combination of $\boldsymbol{\mathcal Z}_{i,j,H,t}$, \[ \operatorname{E}(v^*_{i,j,H,t}\mathbf C_{i,H,t}') =\operatorname{E}(v^*_{i,j,H,t}\mathbf C_{i,j,H,t}^{*\prime})=\mathbf0, \qquad \operatorname{E}(v^*_{i,j,H,t}\varepsilon_{i,t})=0. \] Consequently, $\operatorname{E}(\widetilde y_{j,H,t}y_{i,t})=\operatorname{E}(\widetilde y_{j,H,t}\mu_{i,j,H,t})$. The normalised population screening signal admits the decomposition
which follows directly from the structural equation. This decomposition links the primitive alignment quantities $\bar\gamma_{i,j,k}(H)$ to the population screening signal $m_{i,j}(H)$, and underlies the subsequent detectability and dominance assumptions.
Assumption (ref) formalises the population conditions under which true links remain distinguishable from irrelevant candidates throughout the screening procedure. Together with Assumption (ref), it implies the screening properties required for the subsequent selection proofs, as established in Lemma (ref).
The first part of Assumption (ref) places a primitive restriction on the population IV alignment of proxy links relative to true links. A true link contributes directly to the structural equation for \(y_i\), whereas a proxy link can generate a screening signal only through its association with omitted true links. The assumption therefore requires proxy links to be uniformly less aligned with the remaining structural component than the strongest remaining true link. The appendix shows that, together with Assumption (ref), this primitive restriction implies that, before completion, the strongest remaining true-link screening signal dominates the strongest remaining proxy-link screening signal. The second part imposes the corresponding post-completion requirement: once all true links have been selected, proxy links should have no remaining structural source of population screening signal and therefore should not cross the stopping threshold.
To see the intuition, suppose \(N=5\), \(i=1\), \(\mathcal S_1^n=\{2,3\}\), \(\mathcal S_1^p=\{4\}\), and \(\mathcal S_1^d=\{5\}\). Then \[ y_{1t} = \boldsymbol{\beta}_1'\mathbf x_{1t} + \alpha_{12}y_{2t} + \alpha_{13}y_{3t} + \varepsilon_{1t}. \] For convenience, define the normalised population IV screening signal
Also define the corresponding normalised IV alignment \[ \bar\gamma_{i,j,k}(H) = \frac{ E(\widetilde y_{j,H,t}y_{k,t}) }{ \sigma_{i,j,H} \left\{ E\!\left[ (\widetilde y_{j,H,t})^2 \right] \right\}^{1/2} }. \]
At the first stage, the normalised population IV screening signals are \[ m_{12} = \alpha_{12}\bar\gamma_{1,2,2} + \alpha_{13}\bar\gamma_{1,2,3}, \qquad m_{13} = \alpha_{12}\bar\gamma_{1,3,2} + \alpha_{13}\bar\gamma_{1,3,3}, \] whereas the proxy signal is \[ m_{14} = \alpha_{12}\bar\gamma_{1,4,2} + \alpha_{13}\bar\gamma_{1,4,3}. \] Thus the required dominance condition is \[ \max\{|m_{12}|,|m_{13}|\}>|m_{14}|. \] A proxy link may be correlated with several omitted true links simultaneously. When this occurs, its population screening signal combines information from multiple true neighbours and can, in principle, exceed the marginal signal of an individual omitted true link. If there is only one omitted true link, dominance reduces to the explicit normalised-alignment comparison below; it is not automatic with candidate-specific instruments and normalisations.\footnote{ To see this, suppose that the omitted true-link set is \(\mathcal S_i^n\setminus H=\{s\}\) and let \(p\in\mathcal S_i^p\setminus H\) be a proxy link. The structural equation is \[ y_{i,t} = \boldsymbol{\beta}_i'\mathbf x_{i,t} + \alpha_{i,s}y_{s,t} + \varepsilon_{i,t}. \] Hence the normalised population screening signals are \[ m_{i,s} = \frac{ \alpha_{i,s}E\!\left(\widetilde y_{s,H,t}y_{s,t}\right) }{ \sigma_{i,s,H} \{E[(\widetilde y_{s,H,t})^2]\}^{1/2} }, \qquad m_{i,p} = \frac{ \alpha_{i,s}E\!\left(\widetilde y_{p,H,t}y_{s,t}\right) }{ \sigma_{i,p,H} \{E[(\widetilde y_{p,H,t})^2]\}^{1/2} }. \] Therefore, \[ \frac{|m_{i,p}|}{|m_{i,s}|} = \frac{ \left| E\!\left(\widetilde y_{p,H,t}y_{s,t}\right) \right| }{ \left| E\!\left(\widetilde y_{s,H,t}y_{s,t}\right) \right| } \cdot \frac{ \sigma_{i,s,H} \left\{ E[(\widetilde y_{s,H,t})^2] \right\}^{1/2} }{ \sigma_{i,p,H} \left\{ E[(\widetilde y_{p,H,t})^2] \right\}^{1/2} }. \] Thus, in the single-link case, signal-to-proxy dominance reduces to comparing the normalised population IV alignment of the proxy with that of the true link. Dominance holds precisely when the displayed ratio is uniformly below one. }
Let $\widehat{\mathcal S}_i^t$ denote the set selected by the BOLMT procedure based on the scalar screening statistic $t_{i,j,H}$ and threshold \(c_N^t=C_t\sqrt{\log N}\), where $C_t$ is a sufficiently large fixed constant. Define the exact-recovery event
Theorem (ref) establishes exact recovery of the structural support of every unit. Assumption (ref) prevents premature stopping by requiring that, before the true support has been fully recovered, at least one omitted true link is sufficiently detectable to cross the high-dimensional screening threshold. Assumption (ref) ensures correct selection and stopping: true links dominate proxy candidates before completion, while remaining proxy links do not cross the threshold after the true support has been selected.
For Proposition (ref), assume that every augmented selected model satisfying \[ \mathcal S_i^n\subseteq\widehat{\mathcal S}_i\subseteq \mathcal S_i^n\cup\mathcal S_i^p \] satisfies, uniformly over this finite class, the same IV exogeneity, rank, sample-moment, and long-run variance conditions as the oracle model, with dimensions uniformly bounded.
The preceding results establish recovery of the structural interaction matrix \(\boldsymbol{A}\). As discussed previously, recovery of $\boldsymbol{A}$ is the primary objective of BOLMT and is sufficient in conventional single-network settings. To recover the reinforcing and displacement networks separately, additional structure is required.
Assumption (ref) ensures that the higher-order feedback effects generated by $(\mathbf I_N-\boldsymbol A)^{-1} = \mathbf I_N+\boldsymbol A+\boldsymbol A^2+\cdots$ are represented by a convergent Neumann series for each $N$. This assumption is not used for support recovery or for the sign decomposition below.
For $\ell=1,2$, define the population signed support and row intensity by \[ \mathcal S_{1,i}=\{j:\alpha_{i,j}>0\}, \qquad \mathcal S_{2,i}=\{j:\alpha_{i,j}<0\}, \qquad \rho_{\ell,i}=\sum_{j\in\mathcal S_{\ell,i}}\alpha_{i,j}, \] and set \[ w_{\ell,i,j} =
\] Under Assumption (ref), $\rho_{\ell,i}=\rho_\ell$ on every active row of $\mathbf W_\ell$, while inactive rows are zero.
Theorem (ref) shows that exact support recovery, together with uniform coefficient accuracy, implies exact recovery of the signed supports and consistent estimation of the reinforcing and displacement weights. The sign restrictions and disjoint-support conditions ensure that each recovered link can be assigned uniquely to one of the two interaction channels. The decomposition theorem is stated under common interaction intensities $(\rho_1,\rho_2)$ for expositional simplicity. In empirical applications, however, the interaction intensities and slope coefficients may vary across equations. The post-selection estimation step therefore permits equation-specific parameters and aggregates them using Mean Group estimation. Since support recovery depends only on the interaction matrix $\boldsymbol{A}$, this extension does not affect the preceding recovery theory.
The results above are asymptotic in their nature and are stated under exact sparsity of the relevant local interaction structure. Approximate sparsity can be accommodated by allowing sufficiently small interaction effects to enter the proxy set or approximation error. In that case, however, the target becomes an approximating interaction structure rather than exact support recovery. The present formulation keeps the structural target explicit.
Let $\widehat{\boldsymbol\theta}_i$ be the equation-level post-selection IV estimator and let $\widehat{\boldsymbol\theta}_{MG}$ denote the corresponding mean-group estimator, both defined in Eq. (ref).
Theorems (ref) and (ref) show that network selection does not affect first-order inference for the fixed finite-population target $\boldsymbol\theta_N$. Since exact recovery holds with probability approaching one, the feasible post-selection IV estimator coincides asymptotically with the oracle estimator, and the mean-group estimator has the same limiting distribution as if the true support sets were known. Feasible inference additionally requires a consistent estimator of $\mathbf V$ under the maintained cross-sectional dependence conditions.
The framework developed above extends naturally to settings in which each candidate link contributes multiple regressors. Such specifications arise in a variety of network and spatial models. For example, a single interaction link \((i,j)\) may generate both contemporaneous and lagged regressors, \[ \alpha_{i,j}y_{j,t} + \gamma_{i,j}y_{j,t-1}, \] while spatial Durbin specifications may augment interaction effects with neighbours' covariates, \[ \alpha_{i,j}y_{j,t} + \boldsymbol{\delta}_{i,j}'\mathbf x_{j,t}. \] More generally, a candidate link may contribute a fixed-dimensional block of regressors rather than a single interaction term. In such settings, the object of interest is not an individual coefficient but the presence or absence of a link. Consequently, screening is performed at the block level by jointly testing all regressors associated with a candidate link.
For candidate link \(j\), let \[ \mathbf v_{i,j,H,t} = (v_{i,j,H,1,t},\ldots,v_{i,j,H,d,t})', \] denote the corresponding block of regressors, where the block dimension \(d\) is fixed. Examples include \[ \mathbf v_{i,j,H,t} = (y_{j,t},y_{j,t-1})' \] in dynamic network models, and \[ \mathbf v_{i,j,H,t} = (y_{j,t},\mathbf x_{j,t}')' \] in spatial Durbin specifications.
The theoretical analysis closely parallels the scalar case. The key difference is that signal strength is now measured by a population Wald criterion rather than a scalar IV screening signal.
Define the population fitted block \[ \bar{\mathbf V}^{*}_{i,j,H,t} = \boldsymbol{\Pi}_{V,i,j,H}' \boldsymbol{\mathcal Z}_{i,j,H,t}, \qquad \boldsymbol{\Pi}_{V,i,j,H} = (\mathbf Q^{\boldsymbol{\mathcal Z}}_{i,j,H})^{-1} \operatorname{E}\!\left( \boldsymbol{\mathcal Z}_{i,j,H,t} \mathbf v_{i,j,H,t}' \right), \] and let \[ \boldsymbol{\Gamma}_{VC,i,j,H} = \left\{ \operatorname{E}\!\left( \mathbf C^*_{i,j,H,t} \mathbf C_{i,j,H,t}^{*\prime} \right) \right\}^{-1} \operatorname{E}\!\left( \mathbf C^*_{i,j,H,t} \bar{\mathbf V}_{i,j,H,t}^{*\prime} \right). \] The corresponding population partial fitted block is \[ \mathbf V^*_{i,j,H,t} = \bar{\mathbf V}^*_{i,j,H,t} - \boldsymbol{\Gamma}_{VC,i,j,H}' \mathbf C^*_{i,j,H,t}, \] namely the fitted block after projecting out the fitted selected controls.
In this subsection, let $\boldsymbol\vartheta_{i,k}$ denote the coefficient vector on the regressor block contributed by link $k$, and let $\mathcal S_i^n=\{k:\|\boldsymbol\vartheta_{i,k}\|>0\}$ denote the corresponding group support; the candidate partition and conditioning-set notation are understood groupwise. Reinterpret the omitted component as \[ \mu_{i,j,H,t}:=\mu^F_{i,H,t} = \sum_{k\in\mathcal S_i^n\setminus H} \boldsymbol\vartheta_{i,k}'\mathbf v_{i,k,H,t}. \] Likewise, $\widehat\sigma_{i,j,H}^2$ and $\sigma_{i,j,H}^2$ denote, respectively, the residual-variance estimator from the corresponding block-stage IV regression and its population limit.
Define \[ \mathbf G^*_{i,j,H} = \operatorname{E}\!\left( \mathbf V^*_{i,j,H,t} \mathbf V_{i,j,H,t}^{*\prime} \right), \qquad \mathbf m^F_{i,j,H} = \operatorname{E}\!\left( \mathbf V^*_{i,j,H,t} \mu_{i,j,H,t} \right), \] where \(\mathbf m^F_{i,j,H}\) is the population block screening score vector. Define the normalised population block screening signal
Let \[ \mathbf S_{i,j,H}=T^{-1/2}\widehat{\mathbf V}_{i,j,H}'\widetilde{\mathbf y}_{i,j,H}, \qquad \mathbf G_{i,j,H}=T^{-1}\widehat{\mathbf V}_{i,j,H}'\widehat{\mathbf V}_{i,j,H}. \] The sample block screening quadratic is
which corresponds to the homoskedastic fitted-regressor covariance normalisation. Under serial dependence it is used as a screening quadratic, not as a finite-sample Wald statistic.
For the block results, Assumptions (ref) and (ref) are imposed also on the centred scalar components of \[ V^*_{i,j,H,t,a}\varepsilon_{i,t}, \qquad V^*_{i,j,H,t,a}\mu_{i,j,H,t}, \qquad V^*_{i,j,H,t,a}V^*_{i,j,H,t,b}, \] where \(\mathbf V^*_{i,j,H,t}\) denotes the population partial fitted block.
For the block-screening results, Assumption (ref) is understood to apply also to the centred block score and Gram arrays entering the Wald statistic.
Let $\widehat{\mathcal S}_i^F$ denote the set selected by the block BOLMT procedure based on the statistic $F_{i,j,H}$ and threshold $c_N^F=C_F\log N$, where $C_F$ is a sufficiently large fixed constant. Define the exact-recovery event
Theorem (ref) establishes exact recovery of the structural support under block screening. The result applies to dynamic network models, spatial Durbin specifications, and more general settings in which each link contributes a fixed-dimensional block of regressors.
The block results can be viewed as tests of a candidate group. The population noncentrality \(Q^F_{i,j,H}\) aggregates information across the components of the block. A block may be powerful when the signal is spread across several variables or transformations, even if no single component has a dominant scalar signal. Because the statistic is quadratic, the threshold is of order \(\log N\), and the proxy-dominance condition must dominate both the \(\log N\) stochastic component and the quadratic cross term in the block expansion.
Many applications involve pervasive common shocks that generate cross-sectional dependence beyond the network structure itself. Suppose the disturbance admits the factor representation \[ \varepsilon_{i,t} = \boldsymbol{\lambda}_i' \mathbf f_t + u_{i,t}, \] where the number of factors is fixed.
A factor-augmented implementation includes observed factors or estimated factor proxies among the always-included regressors and in the first-stage projections. Let $\widehat{\mathbf F}$ denote the matrix of observed or estimated factors and let $\mathbf M_{\widehat F}$ denote the corresponding residual-maker matrix. The residual factor component must be negligible on the screening-score scale. For scalar screening this requires \[ \sup_{i,j,H} \left| T^{-1/2} \sum_{t=1}^T v^*_{i,j,H,t} (\mathbf M_{\widehat F}\mathbf F\boldsymbol{\lambda}_i)_t \right| = o_p(\sqrt{\log N}), \] and for block screening it requires \[ \sup_{i,j,H}\max_{1\le a\le d} \left| T^{-1/2} \sum_{t=1}^T V^*_{i,j,H,t,a} (\mathbf M_{\widehat F}\mathbf F\boldsymbol{\lambda}_i)_t \right| = o_p(\sqrt{\log N}). \] If, in addition, the factor-purged idiosyncratic errors satisfy the same mixing, tail, first-stage, and signal conditions as the baseline disturbances, and the feasible factor-replacement bounds in Lemma (ref) hold, then the scalar and block selection results apply to the purged variables.
The role of factor adjustment is particularly important in the present setting because latent common factors provide a natural source of proxy links. Even when a candidate unit does not enter the structural equation directly, it may exhibit a nonzero screening signal if its outcome co-moves with the true links through common macroeconomic, sectoral, or market-wide shocks. In this sense, factor-induced dependence can create spurious network connections that resemble genuine interaction effects.
The preceding argument establishes that the screening and recovery results continue to hold after factor adjustment, provided that the residual factor component is asymptotically negligible and the purged idiosyncratic errors satisfy the maintained regularity conditions. The extension is particularly relevant in economic and financial applications where latent network interactions coexist with pervasive common shocks.
Scalar screening is natural when each candidate link contributes a single interaction regressor. This includes the dual-network model developed in Sections 2--4, where each candidate contributes a contemporaneous interaction term. Block screening is appropriate when each candidate link contributes multiple regressors simultaneously. Examples include dynamic network models involving contemporaneous and lagged interaction effects, spatial Durbin specifications involving links' covariates, and more general contextual-effect models.
The two procedures need not select identical interaction structures in finite samples. A scalar statistic may select a link because one component generates a strong signal. By contrast, a block statistic may select a link whose joint contribution is substantial even when each individual component is only moderately informative. In empirical applications it may therefore be useful to examine the robustness of the selected interaction structure across scalar and block implementations.
The results are asymptotic and are stated under exact sparsity of the relevant local interaction structure. Approximate sparsity can be accommodated by allowing sufficiently small interaction effects to enter the proxy set or approximation error. In that case, however, the target becomes an approximating interaction structure rather than exact support recovery. The present formulation keeps the structural target explicit.
The theory is fundamentally local and proceeds equation by equation. It does not require global graph restrictions such as symmetry, bounded degree, connectedness, or particular topological features of the interaction network. Such restrictions may be useful in specific applications, but they are not needed for the row-wise selection arguments developed here.
Finally, the framework is agnostic regarding the economic mechanism generating the interaction structure. The same methodology applies to spatial, social, financial, production, trade, or other network environments, provided suitable instruments are available and the required signal conditions are satisfied.
In empirical applications, it is useful to report the estimated degree distribution, the sequence of selected links for representative equations, and diagnostics for the associated first-stage regressions. Weak instruments affect both screening and post-selection estimation and may compromise network recovery. In particular, if the fitted-regressor variance is close to zero, the relevance conditions underlying the screening statistics may not be credible for that candidate.
The scalar critical value can be implemented using the threshold $c_p(N,\delta)$ defined in Eq. (ref). For block screening, the analogous critical value is based on Wald-type quantiles with tail probability proportional to $N^{-\delta}$. These choices correspond to thresholds of order $\sqrt{\log N}$ for scalar statistics and $\log N$ for block statistics. The theory does not rely on exact finite-sample normal or Wald distributions. The thresholds should instead be interpreted as high-dimensional screening thresholds; uniform false-selection control requires the sufficiently large leading constants specified in the recovery theorems.
We examine the finite-sample performance of the proposed BOLMT procedure through Monte Carlo experiments based on the following dual-network interaction model:
for $i=1,\ldots,N$ and $t=1,\ldots,T$. The disturbance term satisfies $\varepsilon_{i,t} \sim i.i.d. N(0,\sigma_{\varepsilon}^{2})$, while the individual effects $\eta_i$ are generated independently from a standard normal distribution.
The covariates follow the processes
where the innovation vector $\boldsymbol v_{\ell,t} = (v_{\ell,1,t},\ldots,v_{\ell,N,t})'$ is cross-sectionally correlated according to \[ \boldsymbol v_{\ell,t} \sim N(\boldsymbol 0,\boldsymbol\Sigma_v), \qquad (\Sigma_v)_{ij} = \sigma_v^2\rho_x^{|i-j|}. \] Thus, $\rho_x$ controls both the persistence of the covariates through the autoregressive component and the cross-sectional correlation of the innovations through $\boldsymbol\Sigma_v$. The individual-specific covariate effects are correlated with the fixed effects according to
where $\xi_{\ell,i}\sim i.i.d.N(0,1)$. Accordingly, $\rho_{\eta}$ determines the degree of correlation between the covariates and the individual effects.
We consider three alternative dual-network designs.
Design 1: Circular network. Each unit is connected to four linking units arranged on a circular lattice (see e.g., KapoorKelejianPrucha2007). Two links are assigned to the displacement network $\boldsymbol W_1$ and two links to the reinforcing network $\boldsymbol W_2$. The displacement network is asymmetric, with row weights $(0.75,0.25)$, whereas the reinforcing network is symmetric with weights $(0.5,0.5)$.
Design 2: Block-diagonal network. Units are partitioned into independent groups of size five. Within each block, each unit interacts with four other units. Two links are assigned to $\boldsymbol W_1$ and two links to $\boldsymbol W_2$. As in Design 1, the displacement network is asymmetric and the reinforcing network is symmetric.
Design 3: Random network with a dominant unit. For each unit, links are generated randomly subject to sparsity. Each unit is connected to two displacement links and two reinforcing links. One unit acts as a dominant hub in the reinforcing network, thereby generating heterogeneous network centrality and stronger dependence propagation.
In all designs, the supports of $\boldsymbol W_1$ and $\boldsymbol W_2$ are disjoint and each row is normalised to sum to unity within each network separately.
Let \[ \boldsymbol A = \rho_1\boldsymbol W_1 + \rho_2\boldsymbol W_2. \] The outcome process is generated according to \[ \boldsymbol y_t = (\boldsymbol I-\boldsymbol A)^{-1} ( \boldsymbol\eta + \beta_1\boldsymbol x_{1,t} + \beta_2\boldsymbol x_{2,t} + \boldsymbol\varepsilon_t ). \] The signal-to-noise ratio is calibrated using the spectral radius of $\boldsymbol A$:
Throughout the experiments we fix $SNR=5$ and set $\sigma_{\varepsilon}^{2}=\psi_A/SNR$. The parameter values are fixed at $\beta_1=2, \beta_2=-1, \rho_1=0.5, \rho_2=-0.4$. We consider all combinations of $N\in\{25,50,100,200\}, T\in\{25,50,100,200\}$. Finally, for the critical value function in Eq. ((ref)), we set $\delta=1, c=1, p=0.05$. All experiments are based on $5,000$ Monte Carlo replications.
All variables are first transformed using the within transformation in order to eliminate the individual-specific fixed effects. For notational simplicity, the transformed variables are denoted using the original notation throughout.
In addition to the proposed BOLMT algorithm developed in this paper, we consider two alternative link-selection methods, based on the OCMT and MOCMT procedures of chu2018. These procedures have recently been adapted to network estimation by LiBhattacharjee2024. OCMT applies a single-stage multiple-testing selection rule, whereby all candidate links whose associated $t$-statistics exceed the critical value threshold are selected simultaneously. MOCMT extends this procedure to multiple iterations by updating the conditioning set after each selection round and then re-testing the remaining candidate links. Thus, both OCMT and MOCMT admit all statistically significant candidate links at a given iteration. By contrast, BOLMT proceeds in a more sequential and less greedy manner. At each step, if at least one remaining candidate exceeds the threshold, BOLMT selects only the single link associated with the largest absolute $t$-statistic before updating the conditioning set and continuing to the next step. If no remaining candidate exceeds the threshold, the procedure stops and the current selected set is retained. We index the three procedures by $h=1,2,3$, where $h=1$ corresponds to BOLMT, $h=2$ to OCMT, and $h=3$ to MOCMT.
For each unit $i$, the selected support set obtained from procedure $h$ is denoted by $\widehat{\mathcal S}^{(h)}_i$, where $h=1$ corresponds to BOLMT, $h=2$ to OCMT, and $h=3$ to MOCMT, with cardinality $\widehat{k}^{(h)}_i = |\widehat{\mathcal S}^{(h)}_i|$. Conditional on the selected support, we estimate the restricted structural equation
by instrumental variables.
Define the regressor matrix
and the instrument matrix
The corresponding individual-specific IV estimator is given by
where
The vector $\widehat{\boldsymbol\theta}^{(h)}_{IV,i}$ contains the IV estimates of the slope coefficients together with the estimated interaction coefficients associated with the selected links. Following the decomposition introduced in Section (ref), the estimated interaction coefficients are assigned to the reinforcing and displacement networks according to the sign restrictions imposed in the dual-network representation, yielding the corresponding estimates of $\widehat\rho^{(h)}_{1,i}$ and $\widehat\rho^{(h)}_{2,i}$, together with the associated row-normalised interaction weights.
Finally, the mean-group estimator is constructed by averaging the individual-specific estimates across cross-sectional units:
Table (ref) and Figure (ref) report the finite-sample performance of BLOMT and the aforementioned alternative link-selection procedures across the various network designs. Figures (ref)--(ref) in Appendix B provide additional visual evidence on the decomposition of network recovery performance into true and false link selection rates.
The recovery measures are defined in terms of the true and estimated support matrices. Let \[ s_{i,j} = \mathbf 1(\alpha_{i,j}\neq 0), \qquad \widehat s_{i,j} = \mathbf 1(\widehat\alpha_{i,j}\neq 0), \qquad i\neq j, \] and let $\boldsymbol S=(s_{i,j})$ and $\widehat{\boldsymbol S}=(\widehat s_{i,j})$ denote the corresponding true and estimated support matrices. The mean absolute deviation is \[ \mathrm{MAD} = \frac{1}{N(N-1)} \sum_{i=1}^N \sum_{j\neq i} \left| \widehat s_{i,j}-s_{i,j} \right|. \] Equivalently, MAD is the fraction of pairwise link decisions that are incorrectly classified.
Let \[ \mathrm{TP} = \sum_{i=1}^N\sum_{j\neq i} \mathbf 1(\widehat s_{i,j}=1,s_{i,j}=1), \qquad \mathrm{FP} = \sum_{i=1}^N\sum_{j\neq i} \mathbf 1(\widehat s_{i,j}=1,s_{i,j}=0), \] \[ \mathrm{FN} = \sum_{i=1}^N\sum_{j\neq i} \mathbf 1(\widehat s_{i,j}=0,s_{i,j}=1), \qquad \mathrm{TN} = \sum_{i=1}^N\sum_{j\neq i} \mathbf 1(\widehat s_{i,j}=0,s_{i,j}=0). \] The true positive rate, false positive rate, true discovery rate, and false discovery rate are then defined as \[ \mathrm{TPR} = \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}, \qquad \mathrm{FPR} = \frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}, \] and \[ \mathrm{TDR} = \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}, \qquad \mathrm{FDR} = \frac{\mathrm{FP}}{\mathrm{TP}+\mathrm{FP}}. \] Thus, TPR measures the fraction of true links that are recovered, FPR measures the fraction of non-links that are incorrectly selected, TDR measures the fraction of selected links that are true links, and FDR measures the fraction of selected links that are false links.
The proposed BOLMT procedure performs remarkably well across nearly all Monte Carlo designs. In particular, BOLMT achieves small misclassification rates, with median MAD values close to zero in both the full sample and the restricted sample with $T>25$. The procedure also delivers near-perfect true positive rates and true discovery rates, especially for moderate and large values of $T$. At the same time, both the false positive rate and the false discovery rate remain very small across almost all designs. These findings indicate that the proposed sequential selection mechanism is highly effective in recovering the underlying sparse interaction structure while avoiding spurious link selection.
By contrast, both OCMT and MOCMT exhibit substantially weaker recovery performance. Although the two procedures often attain relatively high true positive rates, this comes at the cost of substantially larger false positive and false discovery rates. In particular, the FDR results indicate that a considerable proportion of the selected links are spurious across most designs. This problem persists even in larger samples and is especially pronounced in the more challenging network structures.
Figure (ref) reports MAD values across all 48 Monte Carlo designs, grouped by network structure and different $(N,T)$ combinations. As shown, BOLMT performs consistently well across all three network structures. Recovery accuracy improves as $T$ increases, and by $T=50$ the procedure delivers very low mis-classification rates in most cases. This pattern remains stable even under the random network with a dominant unit, where higher-order interaction effects are expected to propagate more strongly through the network.
On the other hand, OCMT and MOCMT display noticeably higher MAD values, indicating less accurate support recovery. OCMT performs worst in most designs, especially when $T=25$, while MOCMT usually improves on OCMT but remains well above BOLMT. The deterioration is particularly pronounced under the random network with a dominant unit, where the denser selection rules appear more vulnerable to spurious links generated by indirect dependence propagation.
The differences between the procedures appear to reflect their distinct selection mechanisms. OCMT and MOCMT admit all candidates whose $t$-statistics exceed the multiple-testing threshold at a given iteration. This aggressive rule allows pseudo links generated by indirect dependence propagation to enter alongside true structural links and to retain spurious explanatory power in later iterations. BOLMT, by contrast, selects at most one link at each iteration and updates the conditioning set before reconsidering the remaining candidates. This less greedy strategy reduces premature admission of pseudo links and improves discrimination between direct and indirect dependence.
\begingroup \pagestyle{empty} \thispagestyle{empty}
\endgroup \pagestyle{plain}
Furthermore, for BOLMT, we also examined the accuracy with which the weights associated with the asymmetric displacement network are recovered. In the data generating process, the asymmetric component is constructed using row weights $(0.75,0.25)$, which determine the relative contribution of the two directional components within the displacement network. The BOLMT estimates are tightly concentrated around the true values across all network designs. The largest deviations arise when $T=25$, particularly under the circular network structure, but these discrepancies decline rapidly as the time dimension increases. Overall, the results indicate that the proposed procedure not only recovers the sparse network structure accurately, but also estimates the relative importance of the asymmetric displacement components with high precision.
Finally, Table (ref) reports median computation times per Monte Carlo replication for BOLMT and MOCMT across different network structures and cross-sectional dimensions. As expected, computation times increase substantially with $N$ for both procedures, reflecting the higher-dimensional nature of the network recovery problem. Nevertheless, BOLMT is consistently faster than MOCMT across all designs considered. This finding is noteworthy, since MOCMT admits multiple links simultaneously at each stage, whereas BOLMT proceeds sequentially by selecting only the single most significant link at a time. The computational advantage of BOLMT is particularly pronounced under the random network structure, where the BOLMT-to-MOCMT time ratio falls below 0.5 for all values of $N$.
The results reported in Table (ref) provide a concise overview of the finite sample performance of the competing estimators across the 48 Monte Carlo designs considered in the study. The table reports median absolute bias, median RMSE, and RMSE ratios relative to the Oracle IV estimator across all designs.
For each parameter, absolute bias is computed as the absolute difference between the Monte Carlo mean of the estimator and the true parameter value, multiplied by 100. In particular, for a vector of mean-group estimates $\widehat{\boldsymbol\beta}_{MG}$ and true values $\boldsymbol\beta$, this is computed as \[ \left| \operatorname{mean}\left(\widehat{\boldsymbol\beta}_{MG}\right) - \boldsymbol\beta' \right| \times 100 . \]
The RMSE ratio measures the efficiency loss relative to the infeasible benchmark that treats the true network structure as known. Detailed simulation results for all parameters and estimators are reported in Tables (ref)--(ref) in Appendix B.
First, the proposed BOLMT estimator performs well across all parameters. Median biases remain small for both slope coefficients and spatial parameters, while the corresponding RMSE values are consistently close to those of the Oracle IV benchmark. In particular, the RMSE ratios are close to unity for all parameters, ranging from 1.074 to 1.784. This indicates that the loss in overall estimation accuracy resulting from estimating the underlying network structure is relatively modest. The results therefore suggest that BOLMT is able to recover the relevant network information sufficiently accurately to deliver estimation performance close to the infeasible Oracle estimator that treats the true network structure as known.
Second, the results reveal substantial differences between BOLMT and the alternative model selection procedures. OCMT exhibits considerably larger median biases and RMSE values across all parameters. The deterioration is particularly pronounced for the spatial autoregressive coefficients $\rho_1$ and $\rho_2$, where RMSE ratios exceed 80 and 100, respectively. These findings indicate that OCMT frequently fails to recover the relevant network interactions in finite samples, leading to substantial estimation inaccuracies relative to the Oracle benchmark. MOCMT performs considerably better than OCMT and remains competitive for the slope coefficients $\beta_1$ and $\beta_2$. In particular, the RMSE ratios for these parameters remain close to one, indicating relatively small losses in overall estimation accuracy. This improvement is not surprising, since MOCMT applies the OCMT procedure iteratively, thereby allowing links omitted in earlier stages to be subsequently admitted. However, this additional flexibility may also increase the likelihood of selecting spurious links. While the iterative procedure improves performance relative to OCMT, the results indicate that MOCMT still encounters difficulties in accurately recovering the spatial dependence structure. In particular, the RMSE ratios for $\rho_1$ and $\rho_2$ remain substantially larger than those obtained by BOLMT, which is also reflected in the larger median biases for the spatial coefficients.
Figure (ref) reports the ratio of the RMSE of BOLMT relative to MOCMT for the estimation of the spatial parameter $\rho_2$ across all combinations of $(N,T)$ and the three network structures. Values below unity indicate that BOLMT outperforms MOCMT in terms of estimation accuracy. The solid, dotted, and dashed lines correspond to $\mathbf W_1$, $\mathbf W_2$, and $\mathbf W_3$, respectively.
Several patterns emerge. First, BOLMT consistently dominates MOCMT across all designs, with RMSE ratios well below one throughout. This finding complements the results reported in Table (ref) and suggests that the step-wise, one-link-at-a-time selection rule associated with BOLMT, which limits the inclusion of spurious interaction terms, contributes to the lower estimation error observed across the Monte Carlo designs.
Second, the gains from BOLMT are particularly pronounced in designs with relatively large cross-sectional dimensions. For a given value of $T$, the RMSE ratios generally decline as $N$ increases, indicating that BOLMT becomes especially advantageous when the candidate network is more high-dimensional. In several large-$N$ designs, the ratios fall below 0.1, implying that the RMSE of BOLMT is less than one tenth of that obtained by MOCMT.
Third, the size of the gains depends on network topology. The largest improvements are typically observed under the random network with a dominant unit, $\mathbf W_3$, where MOCMT appears most affected by highly asymmetric and heterogeneous interaction patterns. BOLMT also delivers sizeable gains under the block-diagonal network, $\mathbf W_2$, while the circular network, $\mathbf W_1$, produces smoother patterns and somewhat smaller relative improvements. This is consistent with regular network structures being easier to recover for both procedures.
This section illustrates the practical usefulness of BOLMT by examining network effects in corporate financial policies. The analysis is motivated by the corporate finance literature on peer effects, which shows that firms do not make financial decisions in isolation. Capital structure, liquidity management, investment, and other financial policies may be influenced by economically related firms operating in similar product markets, giving rise to strategic interactions and cross-firm dependence.
A substantial body of empirical work has documented such interactions. Early contributions include LearyRoberts2014, who provide evidence of peer effects in capital-structure decisions using instrumental-variable methods. Subsequent studies, including KaustiaRantala2015, AdhikariAgrawal2018, Grennan2019, and GrieserEtAl2022, employ richer peer-group or network structures to model strategic interactions among firms. In particular, GrieserEtAl2022 develop a spatial econometric framework based on product-market networks constructed from textual similarity measures and find economically important interactions in leverage decisions.
Although these studies differ in how peer groups are defined, they share the common feature that the interaction network is externally specified. In practice, networks are typically constructed using industry classifications, geographic proximity, social connections, or product-market similarity measures, and are then treated as fixed throughout the analysis. The empirical focus is therefore on estimating the strength and implications of interactions within a predetermined network, rather than recovering the network itself. Moreover, interactions are usually modelled through a single channel, so that all links are governed by the same qualitative mechanism. This restriction may be limiting: the relevant network is rarely observed directly, and different links may reflect distinct economic forces. Some interactions may reinforce firms' decisions through imitation, information transmission, learning, or common opportunities, whereas others may generate displacement effects through competition, market-share reallocation, or strategic differentiation.
The methodology proposed in this paper addresses both limitations. BOLMT estimates the network directly from the data by identifying the economically relevant links associated with each firm, rather than imposing a pre-specified interaction structure. The recovered interaction matrix is then decomposed into reinforcing and displacement components, allowing the structure of the network and the nature of the interactions to be determined empirically. In contrast to approaches where inference is conditional on a predetermined network, the present framework treats network recovery as part of the econometric problem and allows positive and negative interaction effects to coexist across different links.
The dataset consists of a panel of 55 publicly listed U.S. firms observed between 1990 and 2023. Firm-level accounting and market variables are obtained from the CRSP/Compustat Merged database, following AsimakopoulosEtAl2026. We consider the following dynamic network specification:
Equivalently, the procedure first recovers the composite interaction coefficients $\alpha_{i,j}=\rho_1w_{1,i,j}+\rho_2w_{2,i,j}$ and then assigns the estimated links to the reinforcing or displacement component according to the sign of $\widehat{\alpha}_{i,j}$, as described in Section (ref). The dependent variable $y_{i,t}$ denotes the firm's leverage ratio, and the vector \[ \mathbf x_{i,t} = \left( \text{Cash}_{i,t}, \text{MB}_{i,t}, \text{Tang}_{i,t}, \text{Prof}_{i,t}, \text{Size}_{i,t} \right)' \] contains the firm-specific characteristics commonly used in the capital-structure literature. Following GrieserEtAl2022, these variables measure cash holdings relative to assets, the market-to-book ratio, asset tangibility, profitability, and firm size, respectively. The parameters $\eta_i$ capture firm-specific fixed effects, while the slope coefficients $\boldsymbol{\beta}_i$ are allowed to differ across firms. The matrices $\mathbf W_1$ and $\mathbf W_2$ represent reinforcing and displacement interaction networks, respectively.
Table (ref) reports Mean Group IV estimates obtained using four alternative network-selection procedures. Our discussion focuses primarily on BOLMT, while the remaining methods are considered for comparison purposes.
Under BOLMT, the estimated autoregressive coefficient is $0.257$ and highly significant, indicating persistence in firms' leverage decisions. The magnitude suggests gradual adjustment over time, although the implied speed of adjustment remains relatively high.
The firm-specific controls generally have the expected signs. Cash holdings enter negatively and are statistically significant at the 10% level, suggesting that firms with greater internal liquidity rely less on leverage. Market-to-book ratios enter positively and are also significant at the 10% level, consistent with firms with stronger growth opportunities adopting more intensive financing policies. Tangibility has a positive and statistically significant coefficient, in line with the collateral role of tangible assets. Profitability enters negatively and is highly significant, consistent with the view that more profitable firms rely less on external financing. Firm size has a positive and significant effect, suggesting that larger firms may benefit from better access to capital markets and greater financial flexibility. These patterns are broadly stable across the alternative network-selection procedures.
The estimated network effects reveal the presence of two distinct interaction mechanisms. The reinforcing effect is positive and highly significant ($\widehat{\rho}_1=0.715$), while the displacement effect is negative and highly significant ($\widehat{\rho}_2=-0.606$). The comparable magnitudes of the two coefficients suggest that both mechanisms are quantitatively important in explaining cross-sectional dependence in firms' financial decisions.
The positive estimate of $\rho_1$ implies that firms connected through the reinforcing network exhibit positively associated leverage adjustments, after controlling for firm-specific characteristics and fixed effects. This is consistent with strategic complementarities arising from imitation, information transmission, learning, or exposure to common opportunities. By contrast, the negative estimate of $\rho_2$ points to a displacement mechanism, under which connected firms move in opposite directions. Such effects may reflect competitive pressures, market-share reallocation, or strategic differentiation. A conventional single-network specification would compress these opposing mechanisms into a single aggregate interaction effect, potentially masking important heterogeneity in the nature of cross-firm dependence. The dual-network framework therefore provides a richer description of strategic interactions in corporate financial policies.
The lower panel of Table (ref) reports network diagnostics. Under BOLMT, the estimated network contains 137 directed links, corresponding to a density of 4.6%. The average degree is 2.5 and the median degree equals 2, indicating a sparse network in which most firms interact with only a small number of economically relevant counterparts. The relatively low density suggests that only a small fraction of all potential interactions are supported by the data. The estimated spectral radius is 0.861, which is well below unity and therefore consistent with the stability condition imposed by the model. The estimated network is also balanced in terms of positive and negative links. Approximately 55% of connections belong to the reinforcing network and 45% belong to the displacement network. This balance provides additional evidence that both interaction mechanisms are empirically relevant.
The remaining columns of Table (ref) provide useful comparisons. BOLMTP is a pruning refinement of BOLMT in which the selected support is re-assessed during the boosting procedure. At each stage, previously selected links are re-tested conditional on the current selected set, and links that no longer remain significant are removed. BOLMTP produces results that are very similar to BOLMT, implying that the stage-wise pruning step has little effect on the recovered support. OCMT identifies a considerably denser network, with 228 links and an average degree exceeding four. MOCMT produces the densest network by a substantial margin, selecting 451 links and yielding an average degree above eight.
The greater density of the OCMT and MOCMT networks is accompanied by larger estimated interaction effects. In particular, MOCMT yields estimates of $\rho_1$ and $\rho_2$ whose absolute magnitudes exceed unity. Moreover, the estimated interaction matrix associated with MOCMT has a spectral radius of 1.334, violating the stability condition imposed by the model. By contrast, the BOLMT estimates imply a substantially sparser interaction structure and a spectral radius below unity, consistent with the theoretical framework developed in this paper.
Table (ref) reports the decomposition of covariate effects into direct, indirect, and total components. Because the model contains contemporaneous interactions, a change in a firm's characteristics affects not only its own outcome but also the outcomes of other firms connected through the estimated networks. These feedback effects propagate through the interaction structure and therefore modify the overall impact of firm-level characteristics. Since the reduced-form representation of the model can be written as $\mathbf y_t = (\mathbf I_N-\boldsymbol{A})^{-1} \left( \boldsymbol{\eta} + \mathbf X_t\boldsymbol{\beta} + \mathbf u_t \right)$, with $\boldsymbol A = \rho_1\mathbf W_1+\rho_2\mathbf W_2$, the matrix of marginal effects associated with covariate $k$ is \[ \mathbf S_k = (\mathbf I_N-\boldsymbol{A})^{-1}\beta_k. \] Following the spatial econometrics literature DebarsyLeSagePace2012,LeSagePace2009, the average direct effect (DE) is computed as the average of the diagonal elements of $\mathbf S_k$, measuring the impact of a change in a firm's own characteristic on its own outcome. The average indirect effect (IE) is computed as the average row sum of the off-diagonal elements of $\mathbf S_k$, capturing the impact transmitted through the network to other firms. The total effect (TE) is the sum of the direct and indirect effects. Because the interaction matrix is decomposed into reinforcing and displacement components, the indirect effect can be further separated into positive and negative contributions associated with $\mathbf W_1$ and $\mathbf W_2$, respectively.
Several findings emerge from Table (ref). First, the direct effects closely resemble the estimated slope coefficients reported in Table (ref), indicating that network feedback does not substantially alter the immediate impact of firm characteristics on financial policies. The indirect effects are nevertheless economically meaningful, increasing the magnitude of the direct effects by approximately 14% on average.
A notable feature of the results is that the indirect effects preserve the sign of the corresponding direct effects. Variables with positive direct effects, namely market-to-book ratios, tangibility, and firm size, also generate positive indirect effects. Similarly, cash holdings and profitability have negative direct and indirect effects. Thus, network propagation reinforces rather than offsets the underlying influence of firm characteristics. Profitability exhibits the largest effect in absolute value. Its direct effect equals $-0.348$, while the total effect reaches $-0.406$, indicating that network interactions strengthen the role of profitability in leverage decisions. More generally, the gap between direct and total effects is non-negligible across all variables, suggesting that part of the overall effect operates indirectly through the recovered network structure.
Overall, the results indicate that direct effects remain the dominant channel, but network propagation is quantitatively relevant. Approximately 14% of the total effect of firm characteristics operates indirectly through the estimated interaction networks.
To investigate whether the estimated network merely reflects similarities in firm size, we compare the distribution of pairwise asset distances for linked and unlinked firm pairs. Let \[ d_{i,j} = \left| \log(\text{Assets}_i) - \log(\text{Assets}_j) \right| \] denote the absolute difference in log assets between firms $i$ and $j$. If the estimated network simply connects firms of similar size, linked pairs should exhibit systematically smaller values of $d_{i,j}$ than unlinked pairs.
Following the nonparametric approach in AsimakopoulosEtAl2026 and KripfganzSarafidis2026, we employ a Wilcoxon rank-sum test to compare the two distributions. We also report the Relative Homophily Index (RHI), \[ \text{RHI} = \frac{ \operatorname{Median}(d_{i,j}\mid \widehat a_{i,j}\neq 0) }{ \operatorname{Median}(d_{i,j}\mid \widehat a_{i,j}=0) }, \] which compares the median asset distance among linked firms with the corresponding median distance among unlinked firms. Values below unity indicate that linked firms tend to be more similar in size than unlinked firms, whereas values close to unity suggest little evidence of size-based homophily.
The results provide little evidence that the recovered network is driven by firm size. The median asset distance among linked pairs equals 9.055, compared with 9.061 among unlinked pairs, yielding an RHI of 1.001. Moreover, the Wilcoxon rank-sum test also fails to reject equality of the two distributions ($p=0.680$). These findings suggest that the links identified by BOLMT do not simply connect firms of similar size and that the recovered network captures information beyond firm-size similarity.
This paper proposes a step-wise IV screening procedure for panel network models with unknown interaction matrices. The procedure estimates local supports one equation at a time, using scalar $t$-statistics or block $F$-statistics constructed from fitted regressors. This allows network links to be recovered directly from the structural equations, rather than being imposed a priori or inferred from reduced-form dependence.
A central feature of the framework is that it accommodates dual interaction structures. The recovered structural interaction matrix can be decomposed into reinforcing and displacement components, allowing positive and negative interaction mechanisms to coexist within the same system. This distinction is important in economic applications where some links may reflect imitation, learning, or common opportunities, while others may arise from competition, market-share reallocation, or strategic differentiation.
The theory is centred on exact structural recovery, but it also clarifies the weaker support-containment property that obtains when the post-completion proxy-control condition is relaxed. This distinction is useful because non-link proxies may generate nonzero population screening criteria through indirect network propagation or common sources of dependence. Exact recovery requires both detectability of omitted true links and sufficient separation between true links and proxy candidates along the selection path. Once exact recovery holds, post-selection IV estimators are asymptotically equivalent to oracle IV estimators that know the true supports.
The probability theory is developed under weak temporal dependence and allows either thin-tailed or heavy-tailed score arrays. The distinction affects the admissible growth of $N$ relative to $T$: thin tails support substantially higher-dimensional candidate sets, whereas heavy tails require polynomial growth restrictions. The framework also extends naturally to block screening, covering dynamic network specifications and spatial Durbin-type models in which each candidate link contributes multiple regressors.
The Monte Carlo experiments show that the proposed BOLMT procedure recovers sparse network structures accurately and delivers post-selection estimates close to the oracle benchmark. Compared with more aggressive selection rules such as OCMT and MOCMT, the step-wise, one-link-at-a-time strategy substantially reduces the inclusion of spurious links and leads to lower estimation error in the designs considered.
The empirical application to U.S. corporate leverage decisions illustrates the practical value of the approach. The results reveal the coexistence of reinforcing and displacement effects in corporate financial policies, suggesting that firms' leverage decisions are shaped by both complementarities and competitive forces. These findings highlight the usefulness of estimating the interaction structure directly and allowing different types of network effects to operate simultaneously.
Overall, the proposed methodology provides a transparent and flexible approach for recovering latent network interactions in high-dimensional panel settings. It offers a way to move beyond pre-specified interaction matrices and single-channel network models, while retaining support-recovery guarantees and oracle-equivalent post-selection inference.