*

#  Robust Estimation in Network Vector Autoregression with Nonstationary
Regressors††thanks: Article history: First draft: 21 April 2021. Second draft:
11 November 2023. The first draft of this paper was prepared after PhD
Candidature was confirmed at the Department of Economics, University of
Southampton. Relevant papers prepared during my PhD studies include Katsouris,
(2021); Katsouris, 2023a ; Katsouris, 2023b ; Katsouris, 2023d . The main
ideas presented in the original draft are extended to Chapter 4 of my PhD
thesis submitted at the University of Southampton titled: "Estimation and
Inference in Seemingly Unrelated Systems of Nonstationary Quantile Predictive
Regression Models". The author is currently working on a follow-up paper
titled: "Robust Identification and Estimation in Non-Gaussian Structural VARs
with Near Unit Roots". Financial support from the Research Council of Finland
(grant 347986) is gratefully acknowledged. Address correspondence to Christis
Katsouris, Faculty of Social Sciences, University of Southampton, United
Kingdom. E-mail<EMAIL_ADDRESS>
Christis Katsouris
University of Southampton
$\&$ University of Helsinki

Job Market Paper I Dr. Christis Katsouris is a Postdoctoral Researcher at the
Faculty of Social Sciences, University of Helsinki.

(August 28, 2024)

###### Abstract

This article studies identification and estimation for the network vector
autoregressive model with nonstationary regressors. In particular, network
dependence is characterized by a nonstochastic adjacency matrix. The
information set includes a stationary regressand and a node-specific vector of
nonstationary regressors, both observed at the same equally spaced time
frequencies. Our proposed econometric specification correponds to the NVAR
model under time series nonstationarity which relies on the local-to-unity
parametrization for capturing the unknown form of persistence of these node-
specific regressors. Robust econometric estimation is achieved using an IVX-
type estimator and the asymptotic theory analysis for the augmented vector of
regressors is studied based on a double asymptotic regime where both the
network size and the time dimension tend to infinity.

Keywords: Network dependence; Local-to-unity; persistence; IVX; Vector
Autoregression.

JEL Classification: C12, C22

##  Introduction

Network connectivity has important implications for the risk management of
economic, financial and societal events such as the diffusion of spillover
effects across large networks (e.g., He and Song,  (2018)), the financial
contagion in stock markets (e.g., Härdle et al.,  (2016), Chen et al.,
(2019), Mitchener and Richardson,  (2019)), the detection of market exuberance
(e.g., Magdalinos and Phillips, 2009a ) as well as the spread of epidemic
diseases (e.g., Keeling and Eames,  (2005)). Time series network-driven models
have seen growing attention in the literature as a statistical mechanism for
identifying underline network structures based on methods commonly used for
Vector Autoregression Models (e.g., Zhu et al.,  (2017, 2019) and Zhu and Pan,
(2020)). Our study is about the aspects of identification and estimation of
cointegration dynamics and the stability of autoregressive processes under
network dependence. We focus on robust econometric estimation of the Network
Vector Autoregression (NVAR), with nonstationary (near-unit-root) node-
specific regressors.

We consider a network with $N$ nodes (possibly high dimensional) such as a
social or a financial network which is indexed by
$i\in\left\\{1,...,N\right\\}$. The network structure is characterized by a
binary nonstochastic adjacency matrix
$\boldsymbol{\Omega}=\left(\omega_{ij}\right)\in\mathbb{R}^{N\times N}$ such
that $\omega_{ij}=1$ if a direct link between the pair $(i,j)$ exists and
$\omega_{ij}=0$ otherwise. We denote with $Y_{i(t)}\in\mathbb{R}$ the
continuous response variable obtained from node $i$ at time point $t$ and
respectively
$\mathbb{Y}_{t}=\left(Y_{1t},...,Y_{Nt}\right)^{\top}\in\mathbb{R}^{N}$ the
possibly high dimensional vector or responses for large $N$. Under the
assumption of identifiable network dependence we assume that the response
variable of the node $i$, $Y_{i(t)}$, is affected by its lag value
$Y_{i,(t-1)}$, as well as by its connected nodes which are collected via
$\left\\{j:\omega_{ij}=1\right\\}$. Moreover, the notion of network-driven
predictability is captured by a set of node-specific variables, denoted by
$X_{i(t)}\in\mathbb{R}^{p}$ which represent regressors of abstract degree of
persistence generated by local-to-unit root processes.

Our research objective is to study the relationship between the stationary
node-specific regressands and a set of node-specific nonstationary regressors
that capture cointegration dynamics under network dependence. Our proposed
framework aims to study economic phenomena under network dependence which
implies that time series observations from the infinite past evolve
conditioning on initial values. Thus we provide an analysis of discrete time,
time invariant, and causal dynamic systems where time, following a finite
sequence of initial values, is explicitly confined to positive integers. In
other words, the unobserved random disturbances determine how random
influences enter the multivariate system. We can think of the random
disturbances of our system to be determining the shocks that enter the system
under the simultaneous presence of cointegration dynamics and network
dependence111Regardless of the presence of network dependence we don’t study
temporal network dynamics under nonstationarity but the simultaneous presence
of nonstationarity and network dependence. Although, we don’t extend our
framework into a high dimensional setting, some relevant studies include
Adamek et al.,  (2022), Chen et al.,  (2023), Krampe et al.,  (2023) and
Zhang,  (2023) as well as Barigozzi et al.,  (2023), Cho et al.,  (2023) and
Fang et al.,  (2023). (see, also Basu and Subba Rao,  (2023)). Recently,
Bykhovskaya,  (2022) consider an estimation and prediction framework where the
model specifies the temporal evolution of a weighted network that combines
classical autoregression with non-negativity, a positive probability of
vanishing and peer effect interactions between weights assigned to edges in
the process.

Although nonstationarity in time series can be interpreted as time-varying
model parameters, we shall focus on cointegrating and unit root dynamics when
the underline data structure is defined across nodes of a single network. We
are interested for conditions of network stationarity against explosiveness of
the underline network evolution process. Regarding suitable mixing condition,
we shall conjecture whether the network effect dominate the nonstationary
property of regressors or whether the network effects appear in the limiting
distribution of the $\mathsf{NVAR}$ estimator.

Generally, we aim to investigate aspects related to the relation between the
network topology and the nonstationarity properties of regressors such as
under which conditions the degree of centrality as captured by the network-
covariates and the nonstationarity as captured by the LUR parametrization
asymptotically dominates when establishing the asymptotic theory of model
parameters. In particular, we shall aim to derive related stability conditions
of the system representation based on the network topology and the properties
of the nonstationary time series regressors. In other words, some of the
challenges that we need to address include the presence of network dependence
and the long memory properties (e.g., see Schennach,  (2018)) of the system
under time series nonstationarity. Some specific questions of interest include
the study of asymptotics and stochastic behaviour of processes in the case of
explosive regressors as well as the asymptotic theory analysis of the IVX
filtration under the presence of network dependence.

To examine both aspects in a unified framework, we begin by studying the
limiting behaviour of the classical least squares estimator within the
stationary framework under the assumption of network dependence. In
particular, the presence of network effects (e.g., the centrality measure as a
covariate), the inclusion of model intercepts (to capture mean nodal effects),
the lagged response variable as well as the nonstationary regressors (such as
persistent regressors), can induce different convergence rates which require
to develop suitable asymptotic theory that includes these features. Firstly,
we need to examine whether the network dependence affects the limiting
behaviour of the OLS type estimator. A starting point is to verify a Gaussian
random variant as a limiting distribution for the stationary framework, in a
similar fashion as in Zhu et al.,  (2017)). The stationary framework, allows
to simplify the econometric specification by assuming that the nodes’
regressors are stationary, so that the LUR process is not employed to model
the regressors.

Secondly, for the nonstationary framework, which is the case that the nodes’
regressors are generated by the LUR specification we need to examine how the
presence of network dependence affects the asymptotic theory of the OLS
estimator. A non-Gaussian limiting distribution can indicate that the presence
of both network effects and nonstationarity requires a different methodology
to robustify inference. To this direction, we can further examine the
implementation of the IVX instrumentation of Phillips and Magdalinos,  (2009),
which has been proven to be robust under the assumption of temporal
dependence, and provides a statistical methodology for filtering the unknown
degree of persistence. Thus, by demonstrating that the asymptotic distribution
of an IVX-type estimator under the assumption of both network dependence and
nonstationarity, converges to a mixed Gaussian distribution, then this
provides a robust methodology for estimation and inference for the NVAR model
within the aforementioned framework.

###  Contributions and Outline of the paper

To the best of our knowledge our study is the first to consider a network-type
of predictability in the same spirit as the conventional predictability
literature based on nonstationary time series regressions (see, Kostakis et
al.,  (2015, 2018). The LUR parametrization is employed to capture the unknown
form of persistence that these nonstationary regressors exhibit and along with
the network-dependent covariates, our functional form specification allows to
model cointegration dynamics in such settings. Our proposed system
representation is novel and thus the closer to our identification and
estimation strategy is the framework presented in the study of Magdalinos,
(2021). In particular, within the predictive regression literature the
regressand corresponds to a continuous stationary covariate such as stock
returns while regressors correspond to persistent data222 A different stream
of literature considers the case in which the response variable represents a
time series of counts, which is useful when modelling corporate defaults as in
Agosto et al.,  (2016)) and Armillotta and Fokianos,  (2022) who consider
modelling count data time series reponsens using nonlinear network vector
autoregressions.. Motivated by these observations we propose a network-
dependent autoregressive representation under the presence of time-series
nonstationarity while we consider suitable stability conditions for an
increasing network size. We believe that our proposed system representation as
well as identification and estimation strategy are of relevance both from the
theoretical and applied econometrics perspective.

Extending existing econometric specifications and estimation approaches to the
case of $\mathsf{NVAR}(1)$ with nonstationary regressors implies overcoming
several challenges to ensure robust implementation approaches. First, it is
not immediately clear which is the most suitable system representation that
incorporates both nonstationary regressors and network dynamics, in the form
of network covariates induced from an non-stochastic adjacency matrix. Second,
while asymptotic theory results for the case of nonstationary predictive
regression models are well developed (see Jansson and Moreira,  (2006),
Phillips and Lee,  (2013), Breitung and Demetrescu,  (2015), Lee,  (2016),
Kostakis et al.,  (2015), Kasparis et al.,  (2015), Andersen and Varneskov,
(2021) and Liu and Phillips,  (2023) among others) and systems of predictive
and cointegrated regressors (see Phillips and Magdalinos,  (2008), Magdalinos
and Phillips, 2009b and Magdalinos,  (2021)), as well as limit results for the
case of network vector autoregressive models under time series stationarity
(see Zhu et al.,  (2017)) and inference techniques under cluster dependence in
two or more dimensions (see Menzel,  (2021), Olmo and Sanso-Navarro,  (2023));
there is no related asymptotic theory that combines these features in a
unified framework (see, also discussion in Katsouris, 2023b ; Katsouris, 2023c
). Third, and perhaps the most profound challenge that would need to be
addressed is the fact that the presence of stationary and nonstationary
regressors, in the form of partially nonstationary and cointegration dynamics
(see, Ahn and Reinsel,  (1990), Toda and Yamamoto,  (1995), Cavanagh et al.,
(1995), Paruolo,  (1997) and Poskitt,  (2006)) along with near unit roots
(see, Phillips and Magdalinos,  (2007)), satisfy certain properties and
regularity conditions which would otherwise be violated under the presence of
network dependence, if these combined features are not being taken care
properly333Suitable markovian conditions are commonly used to preserve the
time series structure when network dependence is jointly modelled in the time
series dimension. Although the markovian property, as captured by the
transition matrix. can be expressed in terms of spatial dependence, in this
paper we consider a more general form of dependence without restricting the
possible dependence structures only to spatio-temporal processes. (see, also
White,  (2000)).

The aforementioned challenges have an impact on inferential procedures and
especially due to the well-known problems with obtaining uniform inference
results to the conventional literature of predictive and cointegrating
regression models (see, Mikusheva,  (2007) and Phillips,  (2014)). Thus, our
contributions in this paper are summarized as below:

  * •

Our econometric environment corresponds to multivariate time series processes
under network dependence. We develop asymptotic theory and estimation
techniques for inference under two large-sample regimes such that $(i)$ with
increasing time sample size, $T\to\infty$, and fixed network dimension,
denoted by $N$, and $(ii)$ with $N\to\infty$ and $T_{N}\to\infty$, where the
temporal size depends on $N$. In other words, the case in which both indices
tend to infinity has certain challenges when deriving the asymptotic
properties or such complex tail dependent processes.

  * •

Our asymptotic theory analysis is developed based on the classical joint weak
convergence arguments of functionals of Brownian motion equipped with $J_{1}$
topology of $\mathcal{D}_{\mathbb{R}^{p}}\left([0,1]\right)$, which holds
regardless of the presence of network dependence. Moreover, we carefully
explain the additional necessary conditions to ensure such that such
topological convergence results are still valid in a more general
autoregressive environment.

The organization of the paper is as follows. Section 2 introduces the Network
Vector Autoregressive Model with nonstationary Regressors. Section 3 develops
large-sample theory for the proposed modelling environment. Section 4 explains
the Monte Carlo Simulation study of the paper and Section 5 explains the
empirical application. Section 6 concludes. Technical proofs can be in the
Appendix of the paper. For any real arbitrary matrix $\boldsymbol{A}$, the
norm is denoted by $\left\lVert\boldsymbol{A}\right\rVert$ and corresponds to
the Frobenius norm defined by
$\left\lVert\boldsymbol{A}\right\rVert=\sqrt{\mathsf{trace}(\boldsymbol{A}^{\prime}\boldsymbol{A})}$.
Let $\lambda_{i}(\boldsymbol{M})$ denote the $i-$th largest eigenvalue of an
$(n\times n)$ symmetric matrix $\boldsymbol{M}$ with its eigenvalues such that
$\lambda_{1}(\boldsymbol{M})\geq...\geq\lambda_{n}(\boldsymbol{M})$. The
spectral norm of $\boldsymbol{A}$ is denoted by
$\left\lVert\boldsymbol{A}\right\rVert_{2}$, such that,
$\left\lVert\boldsymbol{A}\right\rVert_{2}=\sqrt{\lambda_{1}(\boldsymbol{A}^{\prime}\boldsymbol{A})}$,
is maximum column sum norm is denoted by
$\left\lVert\boldsymbol{A}\right\rVert_{1}$, such that,
$\left\lVert\boldsymbol{A}\right\rVert_{1}=\mathsf{max}_{1\leq j\leq
n}\sum_{i=1}^{m}|A_{ij}|$ and its maximum row sum norm is denoted by
$\left\lVert\boldsymbol{A}\right\rVert_{\infty}$, such that,
$\left\lVert\boldsymbol{A}\right\rVert_{\infty}=\mathsf{max}_{1\leq i\leq
n}\sum_{i=1}^{m}|A_{ij}|$. Moreover, the operator $\overset{P}{\to}$ denotes
convergence in probability, and $\overset{D}{\to}$ denotes convergence in
distribution.

###  Prior and related work

A large stream of literature methodologies for modelling cross sectional
dependence and heterogeneity via the use of dynamic panel models. In
particular, Kapetanios et al.,  (2014) develop asymptotic theory for nonlinear
panel models with cross-sectional dependence. Moreover, Huang et al.,  (2020)
propose a network autoregressive model for two-mode networks where the
econometric specification allows for different network autocorrelation
coefficients such that

(1.1) $\displaystyle Y_{1}$
$\displaystyle=\rho_{12}W_{12}Y_{2}+X_{1}\beta_{1}+\varepsilon_{1},$ (1.2)
$\displaystyle Y_{2}$
$\displaystyle=\rho_{21}W_{21}Y_{1}+X_{2}\beta_{2}+\varepsilon_{2},$

where $\varepsilon_{1}$ and $\varepsilon_{2}$ are assumed to be independent.
In other words, the particular structure measures the response collected for
the $k-$th group of nodes. Then, it can be shown that an equivalent model
representation is given by
$Y=\left(I_{N}-\mathbb{W}_{\rho}\right)^{-1}(\mathbb{X}\beta+\varepsilon)$,
where $N=n_{1}+n_{2}$ corresponds to the network size as the sum of the sample
size across the two groups. The estimation of the model parameters is achieved
using the quasi-maximum likelihood approach although robustness is achieved
when heteroscedasticity is properly modelled (see, also Anufriev and
Panchenko,  (2015)). On the other hand, Huang et al.,  (2020) propose an
approximation obtained via the least squares estimator to handle the
computational complexity due to the large-scale network while in our case the
IVX estimator is known to handle the aspects of endogeneity and unknown form
of persistence.

Further examples include the multivariate spatial autoregressive model for
large scale social networks proposed by Zhu et al.,  (2020). All
aforementioned studies correspond to modelling the relationship between the
node-specific responses and a set of exogenous covariates for each node.
Recently, Meitz and Saikkonen,  (2022) propose subgeometrically ergodic
autoregressions with autoregressive conditional heteroscedasticity. The
particular time series approach is applicable for modelling univariate
nonlinear autoregressions with autoregressive conditional heteroscedasticity
of the form:
$y_{t}=\alpha_{1}y_{t-1}+...+\alpha_{p}y_{t-p}+g(u_{t-1})+\sigma_{t}\varepsilon_{t}$.
The types of models we consider in this paper are of the form:
$y^{(i)}_{t}=\beta_{0}+\beta_{1}n_{i}^{-1}\sum_{j=1}^{N}\alpha_{ij}y^{(j)}_{t-1}+\beta_{2}y_{t-1}^{(i)}+\varepsilon_{t}^{(i)}$
as well as
$y^{(i)}_{t}=\beta_{0}+\beta_{1}n_{i}^{-1}\sum_{j=1}^{N}\alpha_{ij}y^{(j)}_{t-1}+\beta_{2}y_{t-1}^{(i)}+X^{\prime}_{i}\gamma+\varepsilon_{t}^{(i)}$,
with the only difference that any node-specific covariates, as specified by
Zhu et al.,  (2017) is replaced by a set of nonstationary regressors. Another
relevant framework is proposed by Nicholson et al.,  (2017) who consider a
high-dimensional VAR model with exogenous variables with a vector of
parameters estimated based on the objective function
$\underset{\boldsymbol{\mu},\boldsymbol{\Phi},\boldsymbol{\beta}}{\mathsf{arg\
min}}\sum_{t=1}^{T}\left\lVert\boldsymbol{y}_{t}-\boldsymbol{\mu}-\sum_{i=1}^{p}\boldsymbol{\Phi}_{(i)}\boldsymbol{y}_{t-j}-\sum_{j=1}^{q}\boldsymbol{\beta}_{(j)}\boldsymbol{x}_{t-j}\right\rVert$.
Moreover, the authors consider an extension to unit-root dynamics but they
tackle the presence of nonstationarity using stationarity transformations
which have the disadvantage of destroying information about the long-run
dynamics of regressors with possible cointegration and nonstationarity.

In our study we study the presence of nonstationarity via the local-to-unity
parametrization to model persistent data which we convert to mildly integrated
using the IVX filtration which controls the degree of endogeneity between the
innovation term of predictors and the innovation term of the predictive
regression. Our proposed approach aims to unify these two aspects; (i) the
network structure, characterized by the adjacency matrix, which allows network
dependence in Vector Autoregression models, and (ii) the persistence
properties of regressors pioneered with the seminal work of Phillips, 1987a ;
Phillips, 1987b and Phillips and Magdalinos,  (2007). To the best of our
knowledge, the inclusion of persistence properties of regressors for the
development of the asymptotic theory of the NVAR model its a novel
contribution to the literature regardless of employing limit theory and
stochastic calculus algebra from existing studies for the development of our
asymptotic theory. The studies of Katsouris,  (2021) and Katsouris, 2023a ;
Katsouris, 2023d , motivated from the aspects of financial connectedness (see,
Diebold and Yılmaz,  (2014) and Baruník and Křehlík,  (2018)) and systemic
risk were the first to bridge the gap in these two streams of literature by
modelling the persistence when estimating systemic risk measures as the tail
risk measures in Adrian and Brunnermeier,  (2016) and Härdle et al.,  (2016).

##  Econometric Identification and Estimation

###  Network Vector Autoregression

Consider a high dimensional network of size $N$ and $Y_{i(t)}$ to be the
response variable which corresponds to node $i$ obtained at time $t$. For each
node $i$, we assume the existence of a $p-$dimensional node specific random
vector of regressors
$X_{i(t)}=\left(X_{i1},...,X_{ip}\right)^{\prime}\in\mathbb{R}^{p}$.
Therefore, to model $Y_{i(t)}$, we propose the following $\mathsf{NVAR}(1)$
model

(2.1) $\displaystyle y_{i(t)}$
$\displaystyle=\beta_{0}+\underbrace{\beta_{1}n_{i}^{-1}\sum_{j=1}^{N}\omega_{ij}y_{j(t-1)}+\beta_{2}y_{i(t-1)}}_{(\text{network
effect})}+AX_{i(t-1)}+\epsilon_{i(t)}$ (2.2)
$\displaystyle\boldsymbol{X}_{i(t)}$
$\displaystyle=\boldsymbol{R}_{n}\boldsymbol{X}_{i(t-1)}+\boldsymbol{U}_{i(t)}$

with $1\leq t\leq T$, where the autoregressive coefficient matrix
$\boldsymbol{R}_{n}$ is expressed as below

(2.3) $\displaystyle\boldsymbol{R}_{n}$
$\displaystyle=\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}_{p}}{T}\right),$

with $\boldsymbol{C}_{p}=\text{diag}\left\\{c_{1},...,c_{p}\right\\}$ such
that $c_{i}$ denotes the unknown coefficient of persistence such that
$n_{i}=\sum_{j\neq i}\omega_{ij}$ is the total number of nodes that $i$
follows, which is the out-degree measure.

(2.4) $\displaystyle\begin{bmatrix}Y_{1}^{i}\\\ \vdots\\\ \vdots\\\
Y_{T}^{i}\end{bmatrix}=\begin{bmatrix}1&1&&...&0\\\ ...&...&&...&...\\\
...&...&&...&...\\\ 1&1&&...&0\end{bmatrix}\begin{bmatrix}X_{1(t-1)}^{i}\\\
\vdots\\\ \vdots\\\
X_{p(t-1)}^{i}\end{bmatrix}^{\prime}+\begin{bmatrix}U_{1(t)}^{i}\\\ \vdots\\\
\vdots\\\ U_{p(t)}^{i}\end{bmatrix}$

In matrix form we have that the regressors are generated via

(2.5) $\displaystyle\begin{bmatrix}X_{1(t)}^{i}\\\ \vdots\\\ \vdots\\\
X_{p(t)}^{i}\end{bmatrix}^{\prime}=\begin{bmatrix}\left(1-\frac{c_{1}}{T}\right)&0&...&...&0\\\
0&\left(1-\frac{c_{2}}{T}\right)&0&...&0\\\ 0&...&....&...&0\\\
0&...&....&...&\left(1-\frac{c_{p}}{T}\right)\end{bmatrix}\begin{bmatrix}X_{1(t-1)}^{i}\\\
\vdots\\\ \vdots\\\
X_{p(t-1)}^{i}\end{bmatrix}^{\prime}+\begin{bmatrix}U_{1(t)}^{i}\\\ \vdots\\\
\vdots\\\ U_{p(t)}^{i}\end{bmatrix}$

for $i\in\left\\{1,...,N\right\\}$. The NVAR(1) model given by (2.1)-(2.3)
implicitly assumes that a particular node $i$ can be affected by another node
$j$, if and only if the pair $(i,j)$ are interconnected as described by the
nonstochastic adjacency matrix $\left[\boldsymbol{\Omega}\right]_{ij}$.
Therefore, in this paper we propose a network vector autoregression (NVAR)
model with nonstationary regressors. The $\mathsf{NVAR}(1)$ model assumes that
each node’s response at a given time time point is a linear combination of (a)
its lag value, (b) the average of its connected neighbours, (c) a set of node-
specific regressors; and (d) an independent noise.

###  Statistical Framework

Recall that $N$ is the network size and $Y_{it}$ is the stationary regressand
that corresponds to the $i-$th subject at time point $t$. In particular, we
focus on the statistical estimation of $\mathsf{NVAR}$ models with
nonstationary regressors for the case that their asymptotic properties depend
on both $N$ (size) and $T$ (time). Therefore, when the dependence of
asymptotic approximations to sample moments relies on both indices, such that
$N$ and $T$, tend to infinity, it requires us to employ a double asymptotic
regime. At the same time an additional challenge we need to address is the
presence of nonstationary regressors when deriving asymptotic theory. However,
within our proposed econometric framework and in contrast to the study of Zhu
et al.,  (2017) and Zhu and Pan,  (2020), we replace the $p-$dimensional node-
specific random vector with a $p-$dimensional node-specific time indexed
vector of nonstationary regressors, which is assumed to be generated using a
local-to-unity parametrization. Specifically, the LUR parametrization for the
nonstationary regressors of the model along with the network-dependent
covariates aims to model cointegration dynamics.

To the best of our knowledge our study is the first to consider a network-type
of predictability in the same spirit as the conventional predictability
literature based on nonstationary time series regressions. In particular,
within the predictive regression literature the regressand corresponds to a
continuous stationary covariate such as stock returns while regressors
correspond to persistent data444 A different stream of literature considers
the case in which the response variable represents a time series of counts,
which is useful when modeling for example corporate defaults as in Agosto et
al.,  (2016)) as well as Armillotta and Fokianos,  (2022) who consider
modelling count data time series reponsens using nonlinear NVARs.. Motivated
from both of these observations we propose a network-dependent representation
of network interactions under the presence of time-series nonstationarity
while we consider suitable stability conditions for increasing network size.
We believe that our proposed system representation as well as identification
and estimation strategy are of relevance both from the theoretical and applied
econometrics perspective. In terms of asymptotic theory properties we consider
sufficient conditions for the information matrix to be stochastic
equicontinuous. Moreover, the main challenge is that the presence of the
nuisance parameter of persistence cannot be consistently estimated which is an
aspect discussed in several papers found in the predictive regression
literature. In particular, for a sequence of parameters
$\boldsymbol{R}(n)=\boldsymbol{I}_{p}-\boldsymbol{C}_{p}/n$, where the real
part of the eigenvalues of $\boldsymbol{C}_{p}\in\mathbb{R}^{p\times p}$ are
all strictly negative, then the statistical problem is equivalent to
estimating
$\boldsymbol{C}_{p}=n\big{(}\boldsymbol{I}_{p}-\boldsymbol{R}(n)\big{)}$ which
implies that in this setting the matrix $\boldsymbol{R}(n)$ can only be
estimated at a rate of $O(n^{-1})$.

Within our framework we consider the uniform convergence of random variables.
Then, for a $\mathsf{VAR}(1)$ process of the form
$\boldsymbol{X}_{t}=\boldsymbol{R}(n)\boldsymbol{X}_{t-1}+\boldsymbol{\varepsilon}_{t}$,
we denote with $\boldsymbol{\vartheta}:=(\boldsymbol{R},\boldsymbol{\Sigma})$
where
$\boldsymbol{\Sigma}:=\mathbb{E}\left(\boldsymbol{\varepsilon}_{t}\boldsymbol{\varepsilon}_{t}^{\prime}|\mathcal{F}_{t-1}\right)$.
Notice that for for every set $\boldsymbol{\vartheta}$, there are infinitely
many processes $\boldsymbol{X}_{t}$ and $\boldsymbol{\varepsilon}_{t}$
satisfying the assumptions above. Consequently, by imposing certain
requirements on the parameter $\boldsymbol{\vartheta}$ and specifically on the
real-valued matrix $\mathbb{R}$, then all these processes are ensured to
exhibit the same limiting behaviour.

We shall clarify that this is not a parameter restriction (as commonly done in
the SVAR literature), but simply an a priori assumption regarding the
asymptotic behaviour of certain stochastic processes (e.g., see ). Moreover,
for every $\boldsymbol{\vartheta}\in\Theta$, suppose that there exists
matrices $\boldsymbol{\Gamma}\in\mathbb{C}^{p\times p}$ and
$\boldsymbol{J}\in\mathbb{C}^{p\times p}$ such that $\boldsymbol{J}$ is a
Jacobian matrix and it holds that
$\boldsymbol{R}(n)\equiv\boldsymbol{\Gamma}\boldsymbol{J}\boldsymbol{\Gamma}^{-1}$.
Then, up to the reordering of the eigenvalues, the matrix $\boldsymbol{J}$ is
unique and satisfies $\boldsymbol{J}\in\mathcal{J}\big{(}|\lambda_{1}|\big{)}$
which is called the Jordan canonical form of the LUR matrix
$\boldsymbol{R}(n)$ (see, Dou and Müller,  (2021)).

Consider the LUR parametrization of a $\mathsf{VAR}(1)$ model expressed as
below

$\displaystyle\boldsymbol{X}_{t}=\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}}{n}\right)\boldsymbol{X}_{t-1}+\boldsymbol{u}_{t}$

Notice that we denote with
$\boldsymbol{\Xi}_{xx}:=\mathbb{E}\big{[}\boldsymbol{X}_{t}\boldsymbol{X}_{t}^{\prime}|\mathcal{F}_{t-1}\big{]}$,
which satisfies

$\displaystyle\boldsymbol{\Xi}_{xx}=\xi^{2}\left(\boldsymbol{I}_{p^{2}}-\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}}{n}\right)\otimes\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}}{n}\right)\right)^{-1}.$

Moreover, we can express the diagonal matrix that includes the nuisance
parameters of persistence

$\displaystyle\boldsymbol{C}$
$\displaystyle\equiv\boldsymbol{\Gamma}\boldsymbol{J}\boldsymbol{\Gamma}^{-1}$
$\displaystyle\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}}{n}\right)$
$\displaystyle\equiv\boldsymbol{\Gamma}\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{J}_{p}}{n}\right)\boldsymbol{\Gamma}^{-1}$

In a similar fashion, lets suppose that we have the local-to-unity model of a
univariate time series such that
$x_{t}=\left(1-\frac{c}{n}\right)x_{t-1}+u_{t}$, then if we assume that there
exists $p$ autoregressive roots such that

$\displaystyle\left(1-\rho_{1}L\right)\left(1-\rho_{2}L\right)...\left(1-\rho_{p}L\right)\big{(}x_{t}-\mu\big{)}\equiv
u_{t}.$

where there is a sequence of nuisance parameters of persistence
$\left\\{c_{j}\right\\}_{j=1}^{p}$. Notice that the above expression
corresponds to the case of a univariate time series which is expressed in the
form of $p$ autoregressive near unit roots in the similar spirit as an $AR(p)$
model has $p$ stationary roots. On the other hand, a VAR$(1)-$model
corresponds to the case where we have, say, a $d-$dimensional (multivariate)
series $\boldsymbol{x}_{t}$, which has near unit roots. Similarly, a
$\mathsf{VAR}(p)-$model implies the presence of $p-$LUR across the
$d-$dimensional vector, which is again different than a multivariate
$\mathsf{VAR}(p)-$model with stationary roots (i.e., meaning stable without
any LUR parametrization). There are of course the econometric specifications
of a $\mathsf{SVAR}(1)$, $\mathsf{SVAR}(p)$ under time series stationarity and
$\mathsf{SVAR}(p)$ under time series nonstationarity which are all
representations beyond the scope of this paper. Next, we shall focus on the
required assumptions that we need to impose on the eigenvalues of the sample
matrices. Notice that we assume that uniform convergence of random matrices
means uniform convergence of the vectorization of these matrices. Although it
is not trivial to conduct inference on the matrix $\boldsymbol{R}_{n}$ due to
the presence of the nuisance parameter $C_{n}$, which cannot be uniformly
consistently estimated, the IVX estimator of the matrix induces mixed Gaussian
asymptotics.

###  Main Assumptions

######  Assumption 1.

Within the stationary framework we impose the following assumptions:

  * (C1)

(Nodal Assumption) Assume that $Z_{i}$’s are i.i.d random variables, with mean
0 and covariance $\Sigma_{z}\in\mathbb{R}^{p\times p}$ and finite fourth order
moments. Moreover, for this case $\left\\{Z_{i}\right\\}$ and
$\left\\{\epsilon_{it}\right\\}$ are assumed to be mutually independent.

  * (C2)

(Network Structure) Assume that $\left\\{\boldsymbol{\Omega}_{i}\right\\}$ is
a sequence of nonstochastic matrices such that $i\in\left\\{1,...,N\right\\}$.

    * (C2.1)

(Connectivity) Treat $\boldsymbol{\Omega}$ as a transition probability matrix
of a Markov Chain, whose state space is defined as the set of all the nodes in
the network, for $i\in\left\\{1,...,N\right\\}$. We assume that the Markov
chain is irreducible and aperiodic. Furthermore, we define
$\pi=\left(\pi_{1},...,\pi_{N}\right)^{\prime}\in\mathbb{R}^{N}$ as the
stationary distribution of the Markov chain, such that (a) $\pi\geq 0$ and
$\sum_{i=1}^{N}\pi_{i}=1$, (b) $\pi=W^{\prime}\pi$. Moreover,
$\sum_{i=1}^{N}\pi_{i}^{2}$ is assumed to converge to 0 as $N\to\infty$.

    * (C2.2)

(Uniformity) Define $W^{*}=W+W^{*}$ as a symmetric matrix. Assume that
$\lambda_{\text{max}}\left(W^{*}\right)=\mathcal{O}\left(\text{log}N\right)$,
where $\lambda_{\text{max}}\left(A\right)$ denotes the largest absolute
eigenvalue of an arbitrary symmetric matrix $M$.

  * (C3)

(Law of Large Numbers) Define $Q:=(I-G)^{-1}(I-G^{\prime})^{-1}$, and recall
$G=\beta_{1}W+\beta_{2}I$. Then, assume that the following limit exists:
$\kappa_{1}=\text{lim}_{N\to\infty}N^{-1}\text{trace}\left\\{\Gamma(0)\right\\}$,
$\kappa_{2}=\text{lim}_{N\to\infty}N^{-1}\text{trace}\left\\{W\Gamma(0)\right\\}$,
$\kappa_{3}=\text{lim}_{N\to\infty}N^{-1}\text{trace}\left\\{(I-G)^{-1}\right\\}$
and $\kappa_{4}=\text{lim}_{N\to\infty}N^{-1}\text{trace}\left\\{Q\right\\}$,
where $\kappa_{1},\kappa_{2},\kappa_{3}$ and $\kappa_{4}$ are fixed constants.

######  Remark 1.

Some relevant remarks regarding our conditions are summarized below:

  * •

Condition (C1) provides a basic assumption regarding the nodal regressors and
the innovation sequence $\epsilon_{it}$. Moreover, we can check how to
correctly impose the covariance structure and the estimation of the
corresponding covariance matrices as in PM. The CLT above does not apply in
the case of stochastic integrals (only valid for the stationary framework).

  * •

Condition (C2) is related to the network structure, as characterized by the
adjacency matrix $W$. Firstly, condition (C2.1) assumes that all nodes are
reachable to each other, that is, irreducibility property. More specifically,
for two arbitrary nodes $i$ and $j$ a path of finite length connecting $i$ and
$j$ should exist. A simple sufficient condition for both irreducibility and
aperiodicity is that the network is always fully connected after a finite
number of steps. That is, there exists an $n^{*}$ such that , for any $n\geq
n^{*}$, each component in $W^{n}$ is always positive. Secondly, condition
(C2.2) requires that the network structure should admit certain uniformity
property so that the diverging rate $\lambda_{\text{max}}\left(W^{*}\right)$
should be sufficiently slow.

For stationary and ergodic time series models the common practice is to
establish the existence and uniqueness of solutions of related recurrent
equations of the form
$\boldsymbol{X}_{t}=\boldsymbol{A}_{t}\boldsymbol{X}_{t-1}+\boldsymbol{B}_{t},t\in\mathbb{Z}$,
as in Fort and Roberts,  (2005), Meitz and Saikkonen,  (2021, 2022) and Matsui
and Pedersen,  (2022) who are exploiting the property that a class of BEKK-
ARCH processes have multivariate stochastic recurrence equation
representations to show the existence of strictly stationary solutions under
mild conditions (see, also Meitz and Saikkonen, 2008a ; Meitz and Saikkonen,
2008b ). Moreover, Doukhan et al.,  (2023) study the stationarity and ergodic
properties for observation-driven models in random environments using strict
exogeneity assumptions. Specifically, the existence of stationary solutions in
the sequential exogeneity framework often relies on additional Lipschitz type
properties. However, in our study we consider a specific functional form which
includes both nonstationary regressors as well as network-dependent
covariates, thereby making it more challenging to either establish a direct
recurrent equation representation or employing other VAR representations as
done in the SVAR literature.

####  Stability Conditions

An alternative way for classifying the stability conditions of the
autoregressive processes which is commonly done by the assumption on the
persistence properties of regressors as in Kostakis et al.,  (2015) and
Magdalinos,  (2021), is to consider relevant eigenvalue conditions on the LUR
coefficient matrix which are based on whether the ordered eigenvalues are
below unity or not (e.g., see Holberg and Ditlevsen,  (2023)). Nevertheless,
such conditions can be also thought as a direct implication of the persistence
class given by Phillips and Magdalinos,  (2009) conditions and the spectral
radius conditions in M $\&$ PCB (2020).

Consider the case where we split the parameter space $\Theta$ into two
overlapping regions with respect to the eigenvalues such that

$\displaystyle R_{n}^{(1)}:=\left\\{\theta\in\Theta:|\lambda_{1}|\leq
1-\frac{\mathsf{log}}{n}\right\\},\ \ \
R_{n}^{(0)}:=\left\\{\theta\in\Theta:|\lambda_{1}|\leq
1-\frac{\mathsf{log}}{n}\right\\}$

Notice that for multivariate time series settings cointegrated systems or
almost cointegrated systems in the case where the roots are only close to
unity. Therefore, due to the presence of asymptotic distribution
discontinuities in order to ensure limit results of known form proper
normalizations of the sample covariances are required555In particular, in
Structural Vector Autoregressive models with nonstationary regressors in order
to ensure valid uniformly inference over a set of parameters including
processes that are cointegrated and with roots arbitrarily close to the unit
circle, adaptive estimation and inference procedures are needed, although that
case is beyond the scope of this paper which are currently investigated in a
different paper.. Usually in the literature the assumption regarding the LUR
matrix is that an assumption is imposed that $\boldsymbol{R}_{n}$ is either a
drifting sequence of diagonal matrices or a fixed diagonal matrix. We focus on
proving that the asymptotic distributions of the relevant sample covariances
can be approximated by stochastic integrals of Orneisten-Uhlenbeck processes.
Therefore, we aim to better understand what is the relation between network
dependence and the unknown degree of persistence.

In particular, it is the case that the OLS estimator is not standard normal
distribution asymptotically, then we can examine the implementation of the IVX
estimator as a filtration method for both the persistence and the network
dependence. In other words, we consider robust estimation and inference
methods for a Network Vector Autoregression Model with LUR regressors. The
model we consider is closely related to the NAR model proposed by Zhu et al.,
(2017) and Zhu et al.,  (2019) as well as the frameworks proposed by
Magdalinos and Phillips, 2009b .

###  Limit theory for Vector Autoregression

####  Weak Convergence Results

To derive the asymptotic theory we employ the local-unit-root specification
proposed by phillips1987time. The regressors are assumed to be generated via
the LUR process $X_{t}=\left(I_{p}-\frac{C_{p}}{T^{\lambda}}\right)$
$X_{t-1}+U_{t}$, with $\lambda<1,\lambda\in(0,1)\ \text{or}\ \lambda>1$. For
example, for the case $\lambda=1$, we consider the following $p-$dimensional
Gaussian process

(2.6) $\displaystyle\displaystyle J_{C}(r)=\int_{0}^{r}e^{C(r-s)}dB_{u}(s),\ \
\ \ r\in(0,1).$

which satisfies the Black-Scholes differential equation $dJ_{c}(r)\equiv
CJ_{C}(r)+dB_{u}(r)$, with $J_{C}(r)=0$, implying also that
$J_{C}(r)\equiv\sigma_{v}J_{C}(r)$, where $\displaystyle
J_{C}(r)=\int_{0}^{r}e^{C(r-s)}dW_{u}(s)$ and $J_{C}(r)$ the Ornstein-
Uhlenbeck, (OU) process,666The OU is a stationary Gaussian process with an
autocorrelation function that decays exponentially over time. Moreover, the
continuous time OU diffusion process has a unique solution. which encompasses
the unit root case such that $J_{C}(r)\equiv B_{u}(r)$, for $C=0$.

######  Theorem 1 (Uniform Convergence of Covariance Matrices (see, Holberg
and Ditlevsen,  (2023))).

Under assumptions above, and for the enlarged probability space
$(\Omega,\mathcal{F},\mathbb{P})$, there exists a standard $d-$dimensional
Brownian motion, denoted by $\left\\{\boldsymbol{W}(t)\right\\}_{t\in[0,1]}$,
and a family of stochastic processes
$\left\\{\boldsymbol{J}_{\boldsymbol{C}}(t)\right\\}_{t\in[0,1],n\in\mathbb{N}}$
such that

$\displaystyle\boldsymbol{J}_{\boldsymbol{C}}(t)=\int_{0}^{t}e^{(t-s)\boldsymbol{C}}\boldsymbol{\Gamma}\boldsymbol{\Sigma}^{1/2}d\boldsymbol{W}(s),\
\ \text{with}\ \ \boldsymbol{J}_{\boldsymbol{C}}(0)=\boldsymbol{0}.$

However, for the remaining of this paper, we consider the special case in
which $\boldsymbol{\Gamma}_{p}\equiv\boldsymbol{I}_{p}$, is the identify
matrix. Moreover, suppose that the LUR matrix $\boldsymbol{R}(n)$ can be
decomposed as $\boldsymbol{R}(n)=\begin{bmatrix}\lambda&1\\\
0&\lambda^{\prime}\end{bmatrix}$. As we have previously mentioned, for every
$\boldsymbol{X}_{t}$ generated by $\boldsymbol{R}(n)\in\mathbb{R}^{p\times
p}$, there exists $\boldsymbol{\Gamma}\in\mathbb{C}^{p\times p}$ and
$\boldsymbol{J}\in\mathcal{J}$ such that $\boldsymbol{\Gamma}$ is invertible
with
$\boldsymbol{\Gamma}^{-1}\boldsymbol{J}\boldsymbol{\Gamma}^{-1}=\boldsymbol{R}(n)$.
We begin our asymptotic theory analysis by considering sequence of parameters
in the stationary region of the LUR matrix $\boldsymbol{R}(n)$.

######  Theorem 2.

Suppose that
$\boldsymbol{V}\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{I}_{p^{2}}\right)$.
For all $\delta>0$ and $r\in[0,1]$, it holds that

$\displaystyle\underset{n\to\infty}{\mathsf{lim}}\
\underset{\theta\in\Theta}{\mathsf{sup}}\
\mathbb{P}\left(\left\lVert\frac{1}{n}\boldsymbol{M}^{-1/2}\left(\sum_{t=1}^{\lfloor
nr\rfloor}\boldsymbol{X}_{t-1}\boldsymbol{X}_{t-1}^{\prime}\right)\boldsymbol{M}^{-1/2}-r\boldsymbol{I}_{p}\right\rVert>\delta\right)=0.$

In other words, for the special case that $r=1$ (full sample sum), then the
above equation shows that the matrix $\boldsymbol{S}_{xx}$ converges in
probability to the identity matrix.

Some further useful results include the following

$\displaystyle\underset{t\geq 1}{\mathsf{sup}}\
\underset{\theta\in\Theta}{\mathsf{sup}}\
\left\lVert\mathbb{E}\big{[}(\boldsymbol{I}-\boldsymbol{R}_{n})^{1/2}\boldsymbol{X}_{t}\boldsymbol{X}_{t}^{\prime}(\boldsymbol{I}-\boldsymbol{R}_{n})^{1/2}\big{]}\right\rVert<\infty.$

Furthermore, it holds that

$\displaystyle\underset{t\geq 1}{\mathsf{sup}}\ \underset{\theta\in
R}{\mathsf{sup}}\
\left\lVert(\boldsymbol{I}-\boldsymbol{R}_{n})^{1/2}\boldsymbol{M}(\boldsymbol{I}-\boldsymbol{R}_{n})^{1/2}-\boldsymbol{\Sigma}_{x}\right\rVert=0.$

where
$\mathsf{vec}(\boldsymbol{\Sigma}_{x})=(\boldsymbol{I}-\boldsymbol{R}_{n})^{1/2}\otimes(\boldsymbol{I}-\boldsymbol{R}_{n})^{1/2\prime}\big{(}\boldsymbol{I}-\boldsymbol{\Gamma}\otimes\boldsymbol{\Gamma}^{\prime}\big{)}^{-1}\mathsf{vec}(\boldsymbol{\Sigma})$.

####  IVX instrumentation

######  Example 1.

Suppose that $\boldsymbol{X}_{t}$ is a $\mathsf{VAR}(1)-$process, then we can
split $X_{t}=\left(Y_{t},\tilde{X}_{t}\right)^{\prime}$ (see also Magdalinos,
(2021)) and the error term $u_{t}=(u_{yt},u_{xt})^{\prime}$ into the first
coordinate and their last $(p-1)$ coordinates

$\displaystyle\boldsymbol{Y}_{t}$
$\displaystyle=\boldsymbol{\beta}\tilde{\boldsymbol{X}}_{t-1}+u_{yt},\tilde{\boldsymbol{X}}_{t}:=\boldsymbol{\Gamma}\boldsymbol{X}_{t}$
$\displaystyle\tilde{\boldsymbol{X}}_{t}$
$\displaystyle=\boldsymbol{R}_{n}\tilde{\boldsymbol{X}}_{t-1}+\boldsymbol{u}_{xt}$

Furthermore, we consider the IVX instrumentation procedure which requires that

(2.7)
$\displaystyle\widetilde{\boldsymbol{Z}}_{i(t)}=\sum_{j=0}^{t-1}\left(\boldsymbol{I}_{p}+\frac{\boldsymbol{C}_{z}}{T^{\lambda}}\right)\big{(}\boldsymbol{X}_{i(t-j)}-\boldsymbol{X}_{i(t-j-1)}\big{)}$

for some $\boldsymbol{C}_{z}=\text{diag}\left\\{c_{z1},...,c_{zp}\right\\}$
with $c_{zi}<0$ and $0<\kappa<1$. Notice that the particular instrumentation
method is considered to be a linear filtering transformation of the regressor
$X_{i(t)}$ into mildly integrated process, as it is explained in Phillips and
Magdalinos,  (2007) and Phillips and Magdalinos,  (2009). However, the novelty
here is the fact that we apply the IVX instrument across the nodes of the
network, which ensures that the regressos of the NVAR(1) model have less
degree of persistence than the original ones. Therefore, the development of a
suitable FCLT requires to consider the weakly convergence of the matrix
moments which are based on the Network Vector Autoregression.

Furthermore, due to the presence of the nuisance parameter $c_{i}$ which
appear in the LUR process the OLS estimator $\widehat{\theta}$ will be biased
(i.e., second degree bias), especially under the assumption of nonzero
covariance terms in the covariance matrix of the innovations of the system. By
substituting the IVX instrument, estimated in the second estimation step
before obtaining the OLS counterpart in the first step, we allow for the
abstract degree of persistence to be filtered out. Moreover, for simplicity in
the derivations of the asymptotic theory we assume that the regressors for all
nodes in the network are identical and belong to the same persistence class as
defined by PM.In particular, the IVX instrumentation corresponds to
endogeneously generated instruments to slow down the rate of convergence of
the estimator enough to ensure mixed Gaussian limiting distributions, which is
based on the assumption that all roots converge to unity at the same rate.
This significant restriction since it implies that it excludes cases where
parts of the process are stationary and other exhibit unit root behaviour such
as mixed integration order, but nevertheless there are some studies which
consider regressors of mixed integration order (see, Phillips and Lee,  (2013)
and M $\&$ PCB (2022)).

###  Limit theory for Network Vector Autoregression

####  Time Series Stationarity

Below we follow the framework proposed by Zhu et al.,  (2017) to motivate
further our study. Therefore, we begin by assuming that the NVAR model has
exogenous regressors which are time-invariant; then the estimation problem
reduces to fitting a Vector Autoregression model with network dependence.
Within this setting, assume for simplicity that we have the following NVAR(1)
model:

(2.8) $\displaystyle Y_{i(t)}$
$\displaystyle=\beta_{0}+\beta_{1}n_{i}^{-1}\sum_{j=1}^{N}\omega_{ij}Y_{j(t-1)}+\beta_{2}Y_{i(t-1)}+\boldsymbol{Z}_{i}^{\prime}\boldsymbol{\xi}+\epsilon_{i(t)}$

Estimating the NVAR(1) model using the econometric specification (2.3),
implies that the pair $(Y_{i},X_{i})$ represents a sequence of stationary
random variables. In this case, we have node-specific covariates which are
time-invariant. This simplification allows to use the conventional central
limit theorem under network dependence for the development of the asymptotic
theory of the model estimates.

Let $\mathbb{Z}=\left(Z_{1},...,Z_{N}\right)^{\prime}\in\mathbb{R}^{N\times
p}$ and
$\boldsymbol{\mathcal{B}}_{0}=\left(\beta_{01},...,\beta_{0N}\right)^{\prime}=\beta_{0}\mathbf{1}+\mathbb{Z}\boldsymbol{\xi}\in\mathbb{R}^{N}$,
where $\textbf{1}=\left(1,...,1\right)^{\prime}$ is the unit vector of the
same dimension,
$\boldsymbol{\xi}=\left(\xi_{1},...,\xi_{p}\right)^{\prime}\in\mathbb{R}^{p}$
and the vector of response variables across the nodes is denoted with
$\mathbb{Y}_{t}=\left(Y_{1t},...,Y_{Nt}\right)^{\prime}\in\mathbb{R}^{N}$.
Then, we rewrite the model (2.3) in the matrix companion form as below

(2.9)
$\displaystyle\mathbb{Y}_{t}=\boldsymbol{\mathcal{B}}_{0}+\boldsymbol{G}\mathbb{Y}_{t-1}+\boldsymbol{\mathcal{E}}_{t}$

where
$\boldsymbol{G}:=\beta_{1}\widetilde{\boldsymbol{\Omega}}+\beta_{2}\boldsymbol{I}$,
with
$\widetilde{\boldsymbol{\Omega}}:=\mathsf{diag}\left\\{n_{1}^{-1},...,n_{N}^{-1}\right\\}\otimes\boldsymbol{\Omega}$,
is the row-normalized adjacency matrix, $\boldsymbol{I}$ is the identity
matrix with compatible dimension and
$\boldsymbol{\mathcal{E}}_{t}=\left(\epsilon_{1t},...,\epsilon_{Nt}\right)\in\mathbb{R}^{N}$
is the innovation vector. We assume the existence of a non-random adjacency
matrix, that captures the network dependence, therefore both $\boldsymbol{G}$
and $\widetilde{\boldsymbol{\Omega}}$ are nonstochastic. However, the matrix
of intercepts $\boldsymbol{\mathcal{B}}_{0}$ is a stochastic quantity which
has to be estimated.

####  Strict Stationarity

We consider the stationary properties of the time series $\mathbb{Y}_{t}$ for
both the cases where the network has a fixed structure, that is $N$ is assumed
to be fixed, as well as the case where the network has an unbounded size, such
that $N\to\infty$.

######  Theorem 3.

Suppose that $\mathbb{E}\left\lVert Z_{i}\right\rVert<\infty$ and $N$ is
fixed. If $|\beta_{1}|+|\beta_{2}|<1$, then there exists a unique strictly
stationary solution with a finite first-order moment to the NVAR(1) model
(2.3). The solution has the following form

(2.10)
$\displaystyle\mathbb{Y}_{t}=\left(\boldsymbol{I}-\boldsymbol{G}\right)^{-1}\boldsymbol{\mathcal{B}}_{0}+\sum_{j=0}^{\infty}\boldsymbol{G}^{j}\boldsymbol{\mathcal{E}}_{t-j}.$

Assuming the existence of the strictly stationary solution (2.23), we are
interested to obtain its conditional distribution given the nodal information
set that characterizes the nodes, denoted with $\mathbb{Z}$. Define with
$\mathbb{E}^{*}(.)=\mathbb{E}\big{(}.|\mathbb{Z}\big{)}$ and
$\mathsf{cov}^{*}=\mathsf{cov}(.|\mathbb{Z})$. For any integer $h$, we denote
the conditional auto-covariance function as
$\boldsymbol{\Gamma}(h)=\mathsf{cov}^{*}\left(\mathbb{Y}_{t},\mathbb{Y}_{t-h}\right)$.
Moreover, it holds that
$\boldsymbol{\Gamma}(h)=\boldsymbol{\Gamma}^{h}(0)\boldsymbol{\Gamma}(0)$ for
$h>0$, and
$\boldsymbol{\Gamma}(h)=\boldsymbol{\Gamma}(0)\left(\boldsymbol{\Gamma}^{\prime}\right)^{-h}$
for $h<0$. Then, the conditional mean and covariance of $\mathbb{Y}_{t}$ are
obtained by Proposition 1.

######  Proposition 1.

Assume the same conditions as in Theorem 3. Then, given $\mathbb{Z}$, the
strictly stationary solution of (2.10) converges to a normal distribution with
mean and covariance given as below

(2.11) $\displaystyle\boldsymbol{\mu}$
$\displaystyle=\left(\boldsymbol{I}-\boldsymbol{G}\right)^{-1}\mathcal{B}_{0}=\left(\boldsymbol{I}-\beta_{1}\boldsymbol{\Omega}-\beta_{2}I\right)^{-1}\mathcal{B}_{0}$
(2.12) $\displaystyle\mathsf{vec}\big{\\{}\boldsymbol{\Gamma}(0)\big{\\}}$
$\displaystyle=\sigma^{2}\big{(}\boldsymbol{I}-\boldsymbol{G}\otimes\boldsymbol{G}\big{)}^{-1}\mathsf{vec}(\boldsymbol{I}).$

The implications of Proposition 1 is that we can determine the factors that
affect the conditional mean of $\mathbb{Y}_{t}$, which are: (i) the nodal
impact $\boldsymbol{\mathcal{B}}_{0}$, (ii) the network effect $\beta_{1}$,
(iii) the momentum effect $\beta_{2}$, and (iv) the network structure777Note
that we assume that the adjacency matrix can take either binary values
indicating the existence of a link between the pair $(i,j)$ or represent a
weight of the strength of connection. given by the adjacency matrix
$\widetilde{\boldsymbol{\Omega}}$.

####  Parameter estimation

The parameters of the NVAR(1) model can be estimated by assuming a Gaussian
innovation process. Let
$\boldsymbol{\beta}=\left(\beta_{0},\beta_{1},\beta_{2}\right)^{\prime}\in\mathbb{R}^{3}$
and
$\boldsymbol{\theta}=(\boldsymbol{\theta}_{j})^{\prime}=\big{(}\boldsymbol{\beta}^{\prime},\boldsymbol{\xi}^{\prime}\big{)}^{\prime}\in\mathbb{R}^{p+3}$.
To estimate the unknown parameter $\boldsymbol{\theta}$, we rewrite the
$\mathsf{NVAR}(1)$ model as below

(2.13) $\displaystyle
Y_{i(t)}=\beta_{0}+\beta_{1}w_{i}^{\prime}\mathbb{Y}_{t-1}+\beta_{2}Y_{i(t-1)}+\boldsymbol{Z}_{i}^{\prime}\boldsymbol{\xi}+\epsilon_{i(t)}=\boldsymbol{\mathcal{X}}_{i(t-1)}^{\prime}\boldsymbol{\theta}+\boldsymbol{\epsilon}_{i(t)},$

where
$\boldsymbol{\mathcal{X}}_{i(t-1)}=\left(1,w_{i}^{\prime}\mathbb{Y}_{t-1},Y_{i(t-1)},Z_{i}^{\prime}\right)^{\prime}\in\mathbb{R}^{p+3}$,
and $w_{i}=\left(\omega_{ij}/n_{i}\right)^{\prime}\in\mathbb{R}^{N}$ for
$1\leq j\leq N$ is the $i-$ row vector of $W$.

Moreover, denote with
$\mathbb{X}_{t}=\left(\boldsymbol{\mathcal{X}}_{1t},\boldsymbol{\mathcal{X}}_{2t},...,\boldsymbol{\mathcal{X}}_{Nt}\right)^{\prime}\in\mathbb{R}^{N\times(p+3)}$.
Then, the NVAR(1) model (2.13) can be rewritten in vector form
$\mathbb{Y}_{t}=\mathbb{X}_{t-1}\boldsymbol{\theta}+\boldsymbol{\mathcal{E}}_{t}$.
Therefore, an ordinary least squares type estimator can be obtained by the
following expression

(2.14)
$\displaystyle\widehat{\boldsymbol{\theta}}_{ols}=\left(\sum_{t=1}^{T}\mathbb{X}_{t-1}^{\prime}\mathbb{X}_{t-1}\right)^{-1}\left(\sum_{t=1}^{T}\mathbb{X}_{t-1}^{\prime}\mathbb{Y}_{t}\right),$

whose asymptotic properties are to be investigated subsequently. The following
conditions are imposed to allow the development of the asymptotic theory.
Within the stationary framework, the asymptotic behaviour of the OLS estimator
for the $\mathsf{NVAR}(1)$ model is only affected by the given structure of
the model, which implies no presence of nuisance parameters such as the
unknown coefficient of persistence.

######  Theorem 4.

Assume that the stationary condition $|\beta_{1}|+|\beta_{2}|<1$ and technical
conditions (C1)-(C3) hold, we then have that

$\displaystyle\sqrt{NT}\left(\widehat{\theta}-\theta\right)\to\mathcal{N}\left(0,\sigma^{2}\Sigma^{-1}\right)$

as min$\left\\{N,T\right\\}\to\infty$, where $\Sigma$.

######  Proposition 2.

Assume that $T$ is fixed and conditions in Theorem 3 hold. Then, we have that

$\displaystyle\sqrt{N}\left(\widehat{\theta}-\theta\right)\to\mathcal{N}\left(0,\sigma^{2}T^{-1}\Sigma^{-1}\right)$

as $N\to\infty$.

######  Remark 2.

For example, recently Fan et al.,  (2023) consider conditions for covariance
stationarity of a quantile vector autoregressive model such that

(2.15)
$\displaystyle\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\big{(}Y_{t}-\mu_{Y}\big{)}\sim\mathcal{N}\left(0,\underset{T\to\infty}{\mathsf{lim}}\
\sum_{t=1}^{T}\mathbb{E}\big{(}Y_{t}-\mu_{Y}\big{)}\big{(}Y_{t}-\mu_{Y}\big{)}^{\prime}\right).$

Furthermore, in the absence of network dependence our econometric
specification reduces to the predictive regression model around the general
vicinity of unity as in Phillips and Magdalinos,  (2009). Usually, in those
frameworks the study of vector autoregressive processes of order 1 implies
that are integrated of order 1 and cointegrated, that is, processes for which
the first difference is stationary and there exists some linear combinations
of the coordinate processes that are stationary.

####  Time Series Nonstationarity

We are interested to develop the asymptotic theory for the case where we have
time varying regressors with certain nonstationary properties. Therefore, we
aim to consider a NVAR model with exogenous nonstationary regressors under a
known network structure. In our case, we are interested to examine the
conditional distribution given the nonstationary regressors of the nodes.

Consider again the model specification with the LUR process as below

(2.16) $\displaystyle Y_{i(t)}$
$\displaystyle=\alpha_{0}+\alpha_{1}n_{i}^{-1}\sum_{j=1}^{N}\omega_{ij}Y_{j(t-1)}+\alpha_{2}Y_{i(t-1)}+BX_{i(t-1)}+\epsilon_{i(t)}$
(2.17) $\displaystyle X_{i(t)}$
$\displaystyle=\boldsymbol{R}X_{i(t-1)}+U_{i(t)}$

where
$\boldsymbol{R}_{T}=\left(\boldsymbol{I}_{p}+\frac{\boldsymbol{C}_{p}}{T^{\lambda}}\right)$,
the autocorrelation coefficient matrix,
$\boldsymbol{C}_{p}=\mathsf{diag}\left\\{c_{1},...,c_{p}\right\\}$ with
$c_{i}$ the unknown persistence coefficient, $\lambda<1,\lambda\in(0,1)$ or
$\lambda>1$, the exponent rate and $\alpha_{0},\alpha_{1},\alpha_{2}$ and
$B=\text{diag}\left\\{\beta_{1},...,\beta_{p}\right\\}$ are the model
parameters to be estimated.

####  Parameter estimation

Denote with
$\widetilde{X}_{i(t-1)}=\left(1,w_{i}^{\prime}\mathbb{Y}_{t-1},Y_{i(t-1)},\widetilde{Z}_{i(t)}^{\prime}\right)^{\prime}\in\mathbb{R}^{p+3}$,
and $w_{i}=\left(\omega_{ij}/n_{i}\right)^{\prime}\in\mathbb{R}^{N}$ for
$1\leq j\leq N$ is the $i-$ row vector of $W$. Moreover, denote with
$\widetilde{\mathbb{X}}_{t}=\left(\widetilde{X}_{1t},\widetilde{X}_{2t},...,\widetilde{X}_{Nt}\right)^{\prime}\in\mathbb{R}^{N\times(p+3)}$
the vector composing the design matrix that aligns with $T$ time observations.
Notice also that in the case where the NVAR model is expressed in a similar
manner as the predictive regression system the design matrix includes the IVX
instrument which is mildly integrated version of the original time varying
regressor.

Then, the NVAR(1) model (2.16) can be rewritten in vector form
$\mathbb{Y}_{t}=\widetilde{\mathbb{X}}_{t-1}\widetilde{\theta}+\mathcal{E}_{t}$.
Therefore, the IVX type estimator of the NVAR(1) model can be obtained by

(2.18)
$\displaystyle\widetilde{\theta}_{\text{IVX}}=\left(\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}^{\prime}\widetilde{\mathbb{X}}_{t-1}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}^{\prime}\mathbb{Y}_{t},$

We aim to show that the asymptotic distribution of the
$\widetilde{\theta}_{\text{IVX}}$ is mixed Gaussian and to determine the
stochastic variance term of its limiting distribution.

######  Theorem 5.

Under Assumption 2 (innovation covariance structure), we then have that

$\displaystyle\sqrt{NT}\big{(}\widehat{\widetilde{\theta}}_{\text{IVX}}-\widetilde{\theta}_{\text{IVX}}\big{)}\to\mathcal{MN}\left(0,\sigma^{2}\mathbb{V}^{-1}\right)$

as min$\left\\{N,T\right\\}\to\infty$, where $\mathbb{V}$ a positive
covariance matrix.

###  General NVAR(p) Model

In the previous sections, we consider the NVAR(1) model, however we can easily
extend the model in the case where $p>1$ while concentrating on the
nonstationary case. Therefore, the NVAR$(p)$ model is expressed as below:

(2.19) $\displaystyle Y_{i(t)}$
$\displaystyle=\beta_{0}+\sum_{s=1}^{p}\alpha_{s}n_{i}^{-1}\sum_{j=1}^{N}\omega_{ij}Y_{j(t-s)}+\sum_{j=1}^{p}\beta_{s}Y_{i(t-s)}+\Xi
X_{i(t-1)}+\epsilon_{i(t)}$ (2.20) $\displaystyle X_{i(t)}$
$\displaystyle=\mathcal{R}_{T}X_{i(t-1)}+U_{i(t)}$

where $\mathcal{R}_{T}=\left(I_{p}+\frac{C_{p}}{T^{\lambda}}\right)$, with
$\lambda<1,\lambda\in(0,1)$ or $\lambda>1$ and
$C_{p}=\text{diag}\left\\{c_{1},...,c_{p}\right\\}$ and
$\Xi=\text{diag}\left\\{\xi_{1},...,\xi_{p}\right\\}$ is the coefficient
matrix for the vector of regressors. Moreover, we define with
$\mathbb{Y}_{t}=\left(\mathbb{Y}_{t}^{\prime},\mathbb{Y}_{t-1}^{\prime},....,\mathbb{Y}_{t-p+1}^{\prime}\right)^{\prime}\in\mathbb{R}^{Np}$.
Then, we express the NVAR$(p)$ model in the the matrix companion form as below

(2.21)
$\displaystyle\mathbb{Y}_{t}^{*}=\mathcal{B}_{0}^{*}+G^{*}\mathbb{Y}_{t-1}^{*}+\mathcal{E}^{*}_{t},$

with
$\mathcal{B}_{0}^{*}=\left(\mathcal{B}_{0}^{\prime},\mathbf{0}^{\prime}_{N(p-1)}\right)\in\mathbb{R}^{Np}$,
$\mathcal{E}^{*}_{t}=\left(\mathcal{E}_{t}^{\prime},\mathbf{0}^{\prime}_{N(p-1)}\right)^{\prime}\in\mathbb{R}^{Np}$
where $G^{*}$ is defined below

(2.22) $\displaystyle G^{*}=\begin{pmatrix}\mathcal{K}\ \ &\ \
\alpha_{p}W+\beta_{p}I_{N}\\\ I_{N(p-1)}\ \ &\ \ \mathbf{0}_{N(p-1),N}\\\
\end{pmatrix},$

such that
$\mathcal{K}=\left(\alpha_{1}W+\beta_{1}I_{N},....,\alpha_{p-1}W+\beta_{p-1}I_{N}\right)\in\mathbb{R}^{N\times
N(p-1)}$, $\mathbf{0}_{n}$ the n-dimensional zero vector and $I_{n}$ the
identity matrix. Furthermore, in the case that we replace the vector of time-
varying regressors $X_{i(t-1)}$ with $Z_{i}$ some time-invariant node-specific
characteristic, then we can consider the strictly stationary solution of the
NVAR$(p)$ model as we presented with Theorem 3 above in the case of the
NVAR(1) model.

######  Theorem 6.

Consider the case where $X_{i(t-1)}\equiv Z_{i}$. Assume that
$\textbf{E}\left\lVert Z_{i}\right\rVert<\infty$ and $N$ is fixed. If
$\sum_{s=1}^{p}\left(|\alpha_{s}|+|\beta_{s}|\right)<1$, then there exists a
unique strictly stationary solution with a finite first-order moment to the
NVAR(p) model (2.19) of the form:

(2.23) $\displaystyle\mathbb{Y}_{t}=\mathcal{J}\mathbb{Y}_{t}^{*},\ \
\text{where}\ \
\mathbb{Y}_{t}^{*}=\left(I_{p}-G^{*}\right)^{-1}\mathcal{B}_{0}^{*}+\sum_{j=0}^{\infty}G^{*j}\mathcal{E}_{t-j}^{*}$

such that the following hold

$\displaystyle\mathcal{J}\left(I-G^{*}\right)^{-1}\mathcal{B}_{0}^{*}=\left(I-\widetilde{G}\right)^{-1}\mathcal{B}_{0},\
\widetilde{G}=\sum_{s=1}^{p}\left(\alpha_{s}W+\beta_{s}I_{N}\right)\
\text{and}\ \mathcal{J}=\left[I_{N}\ \mathbf{0}_{N(p-1),N}\right].$

####  Parameter Estimation for Stationary framework

Assume that the node-specific covariates $Z_{i}$ is a $d-$dimensional vector.
Then, we write
$\mathcal{X}^{*}_{i(t-1)}=\left(1,w_{i}^{\prime}\mathbb{Y}_{t-1},...,\mathbb{Y}_{t-p},Y_{i(t-p)},Z_{i}^{\prime}\right)^{\prime}\in\mathbb{R}^{2p+d+1}$
and
$\mathbb{X}^{*}_{t-1}=\left(\mathcal{X}^{*}_{1(t-1)},...,\mathcal{X}^{*}_{N(t-1)}\right)\in\mathbb{R}^{N\times(2p+d+1)}$.
Moreover, denote the parameter vector with
$\theta^{*}=\left(\beta_{0},\alpha^{\prime},\beta^{\prime},\xi^{\prime}\right)^{\prime}\in\mathbb{R}^{2p+d+1}$,
where $\alpha=\left(\alpha_{1},...,\alpha_{p}\right)^{\prime}$ and
$\beta=\left(\beta_{1},...,\beta_{p}\right)^{\prime}$. Then, model (2.21) can
be written as $\mathbb{Y}_{t}=\mathbb{X}^{*}_{t-1}\theta^{*}+\mathcal{E}_{t}$.
Then, the ordinary least squares type estimator can be obtained by

(2.24)
$\displaystyle\widehat{\theta}_{\text{OLS}}^{*}=\left(\sum_{t=p+1}^{T}\mathbb{X}_{t-1}^{*\prime}\mathbb{X}^{*}_{t-1}\right)^{-1}\sum_{t=p+1}^{T}\mathbb{X}_{t-1}^{*\prime}\mathbb{Y}_{t},$

In order to investigate the asymptotic properties of
$\widehat{\theta}_{\text{OLS}}^{*}$, we define
$\Gamma^{*}(h)=\text{cov}^{*}\left(\mathbb{Y}_{t},\mathbb{Y}_{t-h}\right)$ to
be the conditional auto-covariance function for the NVAR$(p)$ model under the
assumption of strict stationarity. Theorem 7 gives the limiting distribution
of the estimator.

######  Theorem 7.

Assume that $\sum_{s=1}^{p}\left(|\alpha_{s}|+|\beta_{s}|\right)<1$ and
technical conditions (C1),(C2), (C4) hold. We then have that

$\displaystyle\sqrt{NT}\left(\widehat{\theta}^{*}_{\text{OLS}}-\theta^{*}_{\text{OLS}}\right)\to\mathcal{N}\left(0,\sigma^{2}\Sigma^{*-1}\right)$

as min$\left\\{N,T\right\\}\to\infty$, where $\Sigma^{*}$.

####  Parameter Estimation for Nonstationary framework

Under the assumptions we impose for the nonstationary framework, we utilize
the models (2.19)-(2.21) and in particularly we assume the existence of a
vector of nonstationary regressors generated by the LUR specification.
Similarly, as in the previous section we aim to investigate the asymptotic
distribution of the IVX estimator for the NVAR$(p)$.

Denote with
$\widetilde{X}_{i(t-1)}^{*}=\left(1,w_{i}^{\prime}\mathbb{Y}_{t-1},...,\mathbb{Y}_{t-p},Y_{i(t-p)},\widetilde{Z}_{i}^{\prime}\right)^{\prime}\in\mathbb{R}^{2p+d+1}$
where the elements
$w_{i}=\left(\omega_{ij}/n_{i}\right)^{\prime}\in\mathbb{R}^{N}$ for $1\leq
j\leq N$ represent the $i-$ row vector of $W$. Moreover, denote with

$\displaystyle\widetilde{\mathbb{X}}^{*}_{t-1}=\left(\widetilde{X}^{*}_{1(t-1)},\widetilde{X}^{*}_{2(t-1)},...,\widetilde{X}^{*}_{N(t-1)}\right)^{\prime}\in\mathbb{R}^{N\times(2p+d+1)},$

the vector that corresponds to the design matrix of the model that aligns with
$T$ time observations. We also, denote the parameter vector with
$\theta^{*}=\left(\beta_{0},\alpha^{\prime},\beta^{\prime},\Xi\right)^{\prime}\in\mathbb{R}^{2p+d+1}$,
where $\alpha=\left(\alpha_{1},...,\alpha_{p}\right)^{\prime}$ and
$\beta=\left(\beta_{1},...,\beta_{p}\right)^{\prime}$. Therefore, the
NVAR$(p)$ model (2.8) within the nostationary framework can be rewritten in
vector form
$\mathbb{Y}_{t}=\widetilde{\mathbb{X}}^{*}_{t-1}\widetilde{\theta}^{*}+\mathcal{E}_{t}$.
Therefore, the IVX type estimator of the NVAR(1) model can be obtained by

(2.25)
$\displaystyle\widehat{\widetilde{\theta}}^{*}_{\text{IVX}}=\left(\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}^{*\prime}\widetilde{\mathbb{X}}^{*}_{t-1}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}^{*\prime}\mathbb{Y}_{t},$

######  Theorem 8.

Under Assumption 2 (innovation covariance structure), we then have that

$\displaystyle\sqrt{NT}\big{(}\widehat{\widetilde{\theta}}^{*}_{\text{IVX}}-\widetilde{\theta}^{*}_{\text{IVX}}\big{)}\to\mathcal{MN}\left(0,\sigma^{2}\mathbb{V}^{*-1}\right)$

as min$\left\\{N,T\right\\}\to\infty$, where $\mathbb{V}^{*}$ is defined to
be…

The proof of this theorem should be the second theoretical contribution of the
paper to the literature. The proof should be similar to the proof of Theorem
5.

##  Asymptotic Theory

###  Assumptions on Network Dependence

To develop the asymptotic theory of our framework we shall discuss some
relevant assumptions and limit theorems from the network analysis perspective.

####  Asymptotic Uncorrelation

Let $G_{n}=\left(V_{n},E_{n}\right)$ be an undirected network and let
$\left(Z_{i}\right)_{i\in G_{n}}$ be a set of random variables indexed by the
edges. We assume that $\left(Z_{i}\right)_{i\in G_{n}}$ is jointly
exchangeable in the vertices, that is, $\sigma:V_{n}\to V_{n}$ be a
permutation of the vertex set. For an edge $i=(v_{1},v_{2})\in G_{n}$, we
denote by $\sigma_{(i)}:=\left(\sigma_{(v_{1})},\sigma_{(v_{2})}\right)$ the
permuted edge. The network is called exchangeable if $\left(Z_{i}\right)_{i\in
G_{n}}$ and $\left(Z_{\sigma(i)}\right)_{i\in G_{n}}$ have the same joint
distribution for all permutations of the vertex set $\sigma$. In particular,
this means that $Z_{i}$ and $Z_{j}$ are identically distributed if there is a
permutation $\sigma$ such that $i=\sigma(j)$. Notice that this is always the
case because we have defined networks as having no loops, that is, for
$i=(v,v^{\prime})$ we always assume that $v\neq v^{\prime}$. Next, we want to
describe under which circumstances two pairs
$\left(Z_{i_{1}},Z_{j_{1}}\right)$ and $\left(Z_{i_{2}},Z_{j_{2}}\right)$ for
$i_{1},j_{1},i_{2},j_{2}\in G_{n}$ have the same distribution. Denote with
$\kappa(i,j):=|e_{i}\cap e_{j}|\in\left\\{0,1,2\right\\}$ be the number of
common vertices of $i$ and $j$. The following lemma is easy to prove.

######  Lemma 1.

Let $i_{1},j_{1},i_{2},j_{2}\in G_{n}$. There is a permutation of the vertices
$\sigma$ such that $(i_{1},j_{1})=\left(\sigma(i_{2}),\sigma(j_{2})\right)$ if
and only if $\kappa(i_{1},j_{1})=\kappa(i_{2},j_{2})$.

###### Proof.

Let $i=(v_{i_{1}},v^{\prime}_{i_{1}})$ and analogously for $i_{2},j_{1}$ and
$j_{2}$. Let $\sigma$ be a permutation such that
$(i_{1},j_{1})=\left(\sigma(i_{2}),\sigma(j_{2})\right)$. Then, we have that

$\displaystyle\kappa(i_{1},j_{1})=\kappa\left(\sigma(i_{2}),\sigma(j_{2})\right)=\big{|}\left\\{\sigma(v_{i_{2}}),\sigma(v^{\prime}_{i_{2}})\right\\}\cap\left\\{\sigma(v_{j_{2}}),\sigma(v^{\prime}_{j_{2}})\right\\}\big{|}=|e_{i_{2}}\cap
e_{j_{2}}|=\kappa(i_{2},j_{2}).$

Thus, if $\kappa(i_{1},j_{1})=\kappa(i_{2},j_{2})$ we can easily construct
$\sigma$ with $(i_{1},j_{1})=\left(\sigma(i_{2}),\sigma(j_{2})\right)$ just by
mapping the corresponding vertices onto each other and letting the other
vertices unchanged. ∎

Next, we obtain the following corollary which characterizes when two pairs of
random variables are identically distributed.

######  Corollary 1.

Let $\left(Z_{\sigma(i)}\right)_{i\in G_{n}}$ be a sequence of exchangeable
random variables indexed by the edges of a network $G_{n}$. For any vertices
$i_{1},j_{1},i_{2},j_{2}\in G_{n}$, we have that
$\left(Z_{i_{1}},Z_{j_{1}}\right)$ and $\left(Z_{i_{2}},Z_{j_{2}}\right)$ are
identically distributed if $\kappa(i_{1},j_{1})=\kappa(i_{2},j_{2})$.

######  Corollary 2.

For all $n\in\mathbb{N}$, let $G_{n}=(V_{n},E_{n})$ be undirected and complete
networks and assume that $\left(Z_{\sigma(i)}\right)_{i\in G_{n}}$ are
interchangeable and square integrable. Recall that
$r_{n}=|G_{n}|=\frac{n(n-1)}{2}$ is the number of edges. Then, for pairwise
different vertices $v_{1},v_{2},v_{3},v_{4}\in V_{n}$,

$\displaystyle Var\left(\frac{1}{v_{n}}\sum_{i\in G_{n}}Z_{n,i}\right)$
$\displaystyle=\frac{1}{v_{n}^{2}}\sum_{i\in
G_{n}}\text{Var}\left(Z_{n,i}\right)+\frac{1}{v_{n}^{2}}\sum_{\begin{subarray}{c}i,j\in
G_{n}\\\
\kappa(i,j)=1\end{subarray}}Cov\left(Z_{n,i},Z_{n,j}\right)+\frac{1}{v_{n}^{2}}\sum_{\begin{subarray}{c}i,j\in
G_{n}\\\ \kappa(i,j)=0\end{subarray}}Cov\left(Z_{n,i},Z_{n,j}\right)$
$\displaystyle=r_{n}^{-1}Var\left(Z_{n,v_{1}v_{2}}\right)+\mathcal{O}\left(r_{n}^{-1/2}\right)Cov\left(Z_{n,v_{1}v_{2}},Z_{n,v_{2}v_{3}}\right)+\mathcal{O}(1)Cov\left(Z_{n,v_{1}v_{2}},Z_{n,v_{3}v_{4}}\right)$

######  Remark 3.

For instance, if we assume that $Z_{n,i}$ and $Z_{n,j}$ are uncorrelated when
$i\neq j$, the covariances vanish and we see that
$Var\left(\frac{1}{v_{n}}\sum_{i\in G_{n}}Z_{n,i}\right)\to 0$ as $n\to\infty$
if $r_{n}^{-1}Var\left(Z_{n,v_{1}v_{2}}\right)\to 0$ as $n\to\infty$. Notice
that as we have motivated in the previous corollary we do not need asymptotic
uncorrelation for all edges, we need it only for edges $i$ and $j$ with
$\kappa(i,j)=0$. Furthermore, for edges $i$ and $j$ with $\kappa(i,j)\leq 1$
we merely need that the covariances do not grow too fast. Thus, we make an
assumption on the average behaviour of disjoint edges, For $\kappa(i,j)=0$, we
argue that $Cov\left(Z_{n,v_{1}v_{2}},Z_{n,v_{3}v_{4}}\right)\to 0$ is a
reasonable assumption, because we believ that most edges are separated and do
not strongly influence each other.

Consequently, to correctly establish the asymptotic theory within our proposed
econometric environment we shall consider the implementation of suitable
stochastic approximations for functionals of Brownian motion under the
presence of network (graph) dependence. We discuss two relevant studies,
namely the framework proposed by Laurent and Shi,  (2022) as well as the
framework of Gobet and Matulewicz,  (2017). In particular, the latter
considers parameter estimation of Orneisten-Uhlenbeck processes generating a
stochastic graph while the former considers the development of unit root
testing for high-dimensional data. Notice that there are several issues we
need to make sure are clarified. In particular, a more complex scenario would
be the case we allow for an increasing sequence of graphs which have some form
of temporal dependence. On the other hand, to simplify our setting we assume
that we have a fix time length $T$ and consider a random graph $A_{T}$. For
instance, if we would allow for the evolution of the network through time then
we will need to impose related conditions on the correlation structure between
graphs at fixed time length increments.

Suppose that we have a set of stochastic graphs such as the adjacency value
between vertices $i$ and $j$ is given by
$\boldsymbol{\Omega}_{ij}:=\boldsymbol{1}\left\\{\boldsymbol{J}_{C}(t)\in
S_{ij}\right\\}$. In other words, we can consider those network-dependent
covariates as the jump-components in a time series regression for high
dimensional data formulated to represent a stochastic process with unit roots.
There is also a growing literature which considers unstable networks when
modeling interactions of economic agents as in the study of Badev,  (2021).

###  Nonstationary Framework

In order to establish the asymptotic theory in our framework we shall
establish the weak convergence of functionls of Brownian motion in graphs into
OU processes under graph dependence. In this Section we focus on the proofs of
the related theorems that give the asymptotic distribution of the estimators
based on the IVX instrumentation. We follow similar derivations as in Section
LABEL:Section3.1, but for the nonstationary framework we assume that the
vector of regressors $X_{i(t)}$ is generated by the LUR process.

Let
$\widetilde{X}_{i(t-1)}=\left(1,w_{i}^{\prime}\mathbb{Y}_{t-1},Y_{i(t-1)},\widetilde{Z}_{i(t)}^{\prime}\right)^{\prime}\in\mathbb{R}^{p+3}$,
and $w_{i}=\left(\omega_{ij}/n_{i}\right)^{\prime}\in\mathbb{R}^{N}$ for
$1\leq j\leq N$ is the $i-$ row vector of $W$. Moreover, denote with
$\widetilde{\mathbb{X}}_{t}=\left(\widetilde{X}_{1t},\widetilde{X}_{2t},...,\widetilde{X}_{Nt}\right)^{\prime}\in\mathbb{R}^{N\times(p+3)}$.
Then, the NVAR(1) model (2.8) can be rewritten in vector form
$\mathbb{Y}_{t}=\widetilde{\mathbb{X}}_{t-1}\widetilde{\theta}+\mathcal{E}_{t}$.
Therefore, the IVX type estimator of the NVAR(1) model can be obtained by

$\displaystyle\widetilde{\theta}_{\text{IVX}}=\left(\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}^{\prime}\widetilde{\mathbb{X}}_{t-1}\right)^{-1}\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}^{\prime}\mathbb{Y}_{t},$

##### Sketch Proof of Theorem 3

We define with

(3.1)
$\displaystyle\widehat{\Sigma}:=\frac{1}{NT}\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}\widetilde{\mathbb{X}}_{t-1}^{\prime}\
\ \text{and}\ \
\widehat{\Sigma}_{xe}=\frac{1}{NT}\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}\mathcal{E}_{t}$

We have that

(3.2)
$\displaystyle\widehat{\Sigma}=\frac{1}{NT}\sum_{t=1}^{T}\widetilde{\mathbb{X}}_{t-1}\widetilde{\mathbb{X}}_{t-1}^{\prime}=\begin{bmatrix}1&\mathcal{S}_{12}&\mathcal{S}_{13}&\mathcal{S}_{14}\\\
&\mathcal{S}_{22}&\mathcal{S}_{23}&\mathcal{S}_{24}\\\
&&\mathcal{S}_{33}&\mathcal{S}_{34}\\\ &&&\mathcal{S}_{44}\\\ \end{bmatrix}$

where

$\displaystyle\mathcal{S}_{12}$
$\displaystyle=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}w_{i}^{\prime}\mathbb{Y}_{t-1},\
\ \mathcal{S}_{13}=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}Y_{i(t-1)},\ \
\mathcal{S}_{14}=\frac{1}{N}\sum_{t=1}^{T}\widetilde{Z}_{t}^{\prime}$
$\displaystyle\mathcal{S}_{22}$
$\displaystyle=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}\left(w_{i}^{\prime}\mathbb{Y}_{t-1}\right)^{2},\
\
\mathcal{S}_{23}=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}w_{i}^{\prime}\mathbb{Y}_{t-1}Y_{i(t-1)}$
$\displaystyle\mathcal{S}_{24}$
$\displaystyle=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}w_{i}^{\prime}\mathbb{Y}_{t-1}\widetilde{Z}_{t}^{\prime},\
\ \mathcal{S}_{33}=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}Y^{2}_{i(t-1)}$
(3.3) $\displaystyle\mathcal{S}_{34}$
$\displaystyle=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}Y_{i(t-1)}\widetilde{Z}_{t}^{\prime},\
\
\mathcal{S}_{44}=\frac{1}{T}\sum_{t=1}^{T}\widetilde{Z}_{t}\widetilde{Z}_{t}^{\prime}$

where $\widetilde{Z}_{t}$ denotes the IVX instrument.

Therefore, we need to examine the weakly convergence for each of the above
terms in order to determine the covariance of the limiting distribution of the
estimator $\widetilde{\theta}_{\text{IVX}}$.  Convergence of
$\mathcal{S}_{12}$:

(3.4)
$\displaystyle\mathcal{S}_{12}=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}w_{i}^{\prime}\mathbb{Y}_{t-1}=\frac{1}{NT}\sum_{t=1}^{T}\mathbf{1}^{\prime}W\mathbb{Y}_{t-1}=\frac{\beta_{0}}{1-\beta_{1}-\beta_{2}}+\mathcal{S}_{12}^{A}+\mathcal{S}_{12}^{B}$

where

(3.5)
$\displaystyle\mathcal{S}_{12}^{A}=\frac{1}{N}\mathbf{1}^{\prime}W\left(I-G\right)^{-1}\mathbb{Z}\Xi\
\ \text{and}\ \
\mathcal{S}_{12}^{B}=\frac{1}{NT}\sum_{t=1}^{T}\mathbf{1}^{\prime}W\widetilde{\mathbb{Y}}_{t-1}$

To check whether the above terms converge to zero, that is,
$\mathcal{S}_{12}^{A}\to 0$ and $\mathcal{S}_{12}^{B}\to 0$.

Convergence of $\mathcal{S}_{13}$:

(3.6)
$\displaystyle\mathcal{S}_{13}=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}Y_{i(t-1)}=\frac{1}{NT}\sum_{t=1}^{T}\mathbf{1}^{\prime}\mathbb{Y}_{t-1}=\frac{\beta_{0}}{1-\beta_{1}-\beta_{2}}+\mathcal{S}_{13}^{A}+\mathcal{S}_{13}^{B}$

where

(3.7)
$\displaystyle\mathcal{S}_{13}^{A}=\frac{1}{N}\mathbf{1}^{\prime}\left(I-G\right)^{-1}\mathbb{Z}\Xi\
\ \text{and}\ \
\mathcal{S}_{13}^{B}=\frac{1}{NT}\sum_{t=1}^{T}\mathbf{1}^{\prime}\widetilde{\mathbb{Y}}_{t-1}$

We have that $N^{-2}\mathbf{1}^{\prime}Q\mathbf{1}\to 0$ and
$N^{-1}\sum_{j=0}^{\infty}\left\\{\mathbf{1}^{\prime}G^{j}\mathbf{1}\right\\}^{1/2}\to
0$, as $N\to\infty$, which implies that $\mathcal{S}_{13}^{A}\to 0$ and
$\mathcal{S}_{14}^{B}\to 0$. Therefore,
$\mathcal{S}_{13}\to\frac{\beta_{0}}{1-\beta_{1}-\beta_{2}}$.

###  Statistical Hypothesis Testing

In this Section, we test hypotheses of interest based on the IVX-Wald test for
the NVAR(1) model. For instance, we can examine whether a subset of parameters
of the NVAR(1) model is zero. We begin, with a more simplified example since
we are particularly interested to check the presence of both predictability
and causality of certain regressors. The first example, we consider is to test
whether a subset of parameters of the NVAR(1) model are simultaneously zero
across the nodes.

To construct such a test we define
$\mathcal{I}\subset\left\\{(j,r,s):j,r\in\left\\{1,...,p\right\\}\ \text{and}\
s\in\left\\{1,...,d\right\\}\right\\}$ be a subset of indices due to the
presence of both nonstationary regressors and lagged $Y_{i,(t)}$ variables as
predictors. We consider the following testing hypothesis:

$\displaystyle H_{0}$ $\displaystyle:\theta^{(s)}_{jr}=0,\text{for all}\
(j,r,s)\in\mathcal{I}$ $\displaystyle H_{1}$ $\displaystyle:\text{There exists
at least one pair}\ (j,r,s)\in\mathcal{I}\ \text{such that}\
\theta^{(s)}_{jr}\neq 0.$

Here we need to check whether the linear restrictions imposed on the parameter
space of the NVAR(1) model under both the null and the alternative hypothesis
do not violate the stability of the model. For example, whether we can
establish that det$\left(\theta^{*}(z)\right)\neq 0$, for all $z\leq 1$.
Furthermore, although the presence effects and nonlinearities in network
dependence is an interesting research avenue (e.g., see Armillotta and
Fokianos,  (2022), Teräsvirta et al.,  (1994)), we focus on testing the
simultaneous presence of nonstationarity and network dependence, which is a
novel test on its own right in the existing literature, using our proposed
estimation and inference methodology for NVAR models with nonstationary
regressors.

##  Empirical Application

We present an empirical implementation based on the dataset of Härdle et al.,
(2016). The particular dataset contains the stock returns of the top 100 US
financial institutions by market capitalization with weekly time observations
spanning the trading period from 2007 to 2013 as well as a set of
macroeconomic and financial variables with unknown persistence properties. We
begin by explaining the procedure to estimate the nonstochastic adjacency
matrix for the full sampling period.

In particular, we construct a tail network which is a suitable to capture tail
interconnectedness (see, Chen et al.,  (2019)), that is, financial
conntectedness estimated with tail risk measures (see, also Daouia et al.,
(2022)). Suppose that $t\in\left\\{1,...,T\right\\}$, then by employing a
quantile optimization function888Consider the conditional quantile function
estimated via
$Q_{y_{i}}\left(\tau|x_{i}\right)=F_{y_{i}}^{-1}\left(\tau|x_{i}\right)$. The
optimization function to obtain the model estimates is expressed as
$Q_{\tau}\left(y_{i}|x_{i}\right)=argmin_{q(x)}\
\mathbb{E}\big{[}g_{\tau}\left(y_{i}-q(x_{i})\right)\big{]}$, where
$\tau\in(0,1)$ is a specific quantile level, and
$g_{\tau}(u)=u\left(\tau-\mathbf{1}{\left\\{u<0\right\\}}\right)$ is the check
function. we can estimate the risk measures of the VaR and CoVaR via the
following

(4.1) $\displaystyle R_{i,t}$
$\displaystyle=a_{i}+b_{i}^{\prime}X_{t-1}+u_{i,t}$ (4.2) $\displaystyle
R_{j,t}$
$\displaystyle=a_{j|i}+b_{j|i}^{\prime}X_{t-1}+\gamma_{j|i}R_{i,t}+u_{j|i,t}$

where $\\{R_{i,t}\\}_{i=1,...,N}$ is the vector of portfolio returns at time
$t$, and $X_{t-1}$ is a vector of exogenous regressors containing
macroeconomic characteristics common across assets. Let $Q_{\tau}\left(\cdot\
|\ \mathcal{F}_{t-1}\right)$ denote the quantile operator for $\tau\in(0,1)$
conditional on an information set $\mathcal{F}_{t-1}$. Then, the error terms
satisfy $Q_{\tau}\left(u_{i,t}\ |\ X_{t-1}\right)=0$ for the quantile
regression (4.1), and $Q_{\tau}\left(u_{j|i,t}\ |\ R_{it},X_{t-1}\right)=0$
for the quantile regression (4.2).

Denote $\theta_{(1)}:=\left(a_{i},b_{i}^{\prime}\right)$ and
$\check{X}_{t-1}^{(1)}:=\left(1,X_{t-1}^{\prime}\right)^{\prime}$ for (4.1)
and $\theta_{(2)}:=\left(a_{j|i},b_{j|i}^{\prime},\gamma_{j|i}\right)$ and
$\check{X}_{t-1}^{(2)}:=\left(1,X_{t-1}^{\prime},R_{i,t}\right)^{\prime}$ for
(4.2). Then, the QR estimator for model (4.1) is given by

(4.3)
$\displaystyle\widehat{\theta}_{(1)}\left(\tau\right):=\underset{\theta_{(1)}\in\mathbb{R}^{p+2}}{\text{argmin}}\sum_{t=1}^{T}g_{\tau}\left(R_{i,t}-\theta_{(1)}^{\prime}\check{X}_{t-1}^{(1)}\right)$

Similarly, the quantile estimator of (4.2) is obtained by

(4.4)
$\displaystyle\widehat{\theta}_{(2)}\left(\tau\right):=\underset{\theta_{(2)}\in\mathbb{R}^{p+1}}{\text{argmin}}\sum_{t=1}^{T}g_{\tau}\left(R_{j,t}-\theta_{(2)}^{\prime}\check{X}_{t-1}^{(2)}\right)$

The QR estimator from the optimization function (4.3) allows to construct the
one-period ahead forecast for the Value-at-Risk of firm $i$ such that
$\hat{\text{VaR}}_{i,t+1}(\tau)=\hat{a}_{i}+\hat{b}_{i}X_{t}$. Similarly, we
can obtain the one-period ahead forecast for the Conditional-Value-at-Risk of
firm $i$ such that
$\hat{\text{CoVaR}}_{j|i,t+1}(\tau)=\hat{a}_{j|i}+\hat{b}_{j|i}X_{t}+\hat{\gamma}_{j|i}\hat{VaR}_{i,t+1}(\tau)$.
This procedure allows to construct the VaR-$\Delta$CoVaR risk matrix proposed
by katsourisOlmo20. In this paper, we consider the construction of a binary
adjacency matrix $\Omega_{ij}$, with elements such that $\omega_{ij}=1$, if
$\gamma_{i|j}$ is individually statistical significant and $\omega_{ij}=0$
otherwise.

The above approach allows to obtain a non-random adjacency matrix which
captures the tail connectivity of the nodes for the whole sampling period. The
second step, is to estimate the NVAR(1) model assuming the regressors are
generated by the local-to-unit root process as described in Section 2.1. The
testing hypothesis we are interested to examine is whether there is joint
predictability and causality, that is, there exists some $\xi_{i}\neq 0$, for
$i\in\left\\{1,...,p\right\\}$ which is a subset of the parameter space of
$\xi$ and $\beta_{2}\neq 0$. Under the null hypothesis there is no joint
predictability and no causality while under the alternative hypothesis there
is simultaneously joint predictability and causality.

##  Conclusion

In this paper, we consider a high-dimensional network and propose the Network
Vector Autoregression model as a suitable statistical methodology to capture
the network dynamics. The particular NVAR model we propose, allows to
incorporate a nonstochastic adjacency matrix as well as it captures the time
series properties of regressors via the local-to-unit root specification.
Moreover, we examine the asymptotic theory of the model within the framework
of network dependence and conditional neighborhood dependence presented by
Kojevnikov et al.,  (2021) and Lee and Song,  (2019) respectively. More
specifically, we develop the asymptotic theory for the IVX estimator of the
NVAR model and show that the estimator is robust to the degree of persistence
of regressors, and specify the matrix moments of its limiting distribution. To
summarize, in this paper we consider the persistence properties of the
regressors which correspond to each node in the network and develop an
asymptotic theory based on the IVX instrumentation. An empirical application
and two additional applications demonstrate the usefulness of our framework.

##### Conflicts of interest

The author declares that there are no known conflicts of interest.

##### Data availability

No data was used for the research described in this article.

##### Acknowledgements

I wish to thank Professor Jose Olmo and Professor Tassos Magdalinos from the
Department of Economics, University of Southampton for helpful discussions as
well as Dr. Julius Vainora from the Faculty of Economics, University of
Cambridge. Moreover, I am grateful to Professor Markku Lanne and Professor
Mika Meitz from the Faculty of Social Sciences, University of Helsinki for
helpful conversations. Financial support from the Research Council of Finland
(grant 347986) is gratefully acknowledged. All remaining errors are my own
responsibility.

##  Appendix

###  Illustrative Example of VAR Representations

######  Example 2.

We consider the multivariate autoregressive index model representation
proposed by Reinsel,  (1983), which provide a suitable parametrization for
demensionality reduction. Specifically, consider an $n-$dimensional stationary
vector autoregressive time series denoted by
$\boldsymbol{Y}_{t}=(y_{1t},...,y_{nt})^{\prime}$

(6.1)
$\displaystyle\boldsymbol{Y}_{t}-\sum_{j=1}^{p}\boldsymbol{\Phi}_{j}\boldsymbol{Y}_{t-j}=\boldsymbol{\varepsilon}_{t},$

where $\boldsymbol{\varepsilon}_{t}$ are i.i.d
$\mathcal{N}(0,\boldsymbol{\Omega})$. Define with
$\boldsymbol{\Phi}(L)=\boldsymbol{I}-\boldsymbol{\Phi}_{1}L-...-\boldsymbol{\Phi}_{p}L^{p}$,
where $L$ denotes the lag operator, and assume that
$\mathsf{det}[\boldsymbol{\Phi}(z)]\neq 0$ for all complex numbers $|z|\leq
1$. In addition, we consider that
$\boldsymbol{I}-\boldsymbol{\Phi}(L)=\sum_{j=1}^{p}\boldsymbol{\Phi}_{j}L^{j}$
can be factorized as
$\boldsymbol{I}-\boldsymbol{\Phi}(L)\equiv\boldsymbol{A}(L)\boldsymbol{B}(L)$
such that

(6.2) $\displaystyle\boldsymbol{A}(L)_{(m\times r)}$
$\displaystyle=\boldsymbol{A}_{1}L+...+\boldsymbol{A}_{p_{1}}L^{p_{1}}$ (6.3)
$\displaystyle\boldsymbol{B}(L)_{(r\times m)}$
$\displaystyle=\boldsymbol{B}_{0}+\boldsymbol{B}_{1}L+...+\boldsymbol{B}_{p_{2}}L^{p_{2}}$

Therefore, we have that

(6.4)
$\displaystyle\boldsymbol{Y}_{t}=\boldsymbol{A}(L)\boldsymbol{B}(L)+\boldsymbol{\varepsilon}_{t}\equiv\sum_{\mathsf{u}=1}^{p_{1}}\sum_{\mathsf{v}=1}^{p_{2}}\boldsymbol{A}_{\mathsf{u}}\boldsymbol{B}_{\mathsf{v}}\boldsymbol{Y}_{t-\mathsf{u}-\mathsf{v}}+\boldsymbol{\varepsilon}_{t}.$

where $p_{1}+p_{2}=p$. Suppose that $\boldsymbol{\eta}_{t}$ is an
$r-$dimensional series such that
$\boldsymbol{\eta}_{t}=\boldsymbol{B}(L)\boldsymbol{Y}_{t}\equiv\boldsymbol{B}_{0}\boldsymbol{Y}_{t}+\boldsymbol{B}_{1}\boldsymbol{Y}_{t-1}+...+\boldsymbol{B}_{p_{2}}\boldsymbol{Y}_{t-p_{2}}$
then the variables
$\left\\{\boldsymbol{\eta}_{t-1},...,\boldsymbol{\eta}_{t-p_{1}}\right\\}$,
provide an adaptive filtration of all past information required for
prediction. In particular, we have that

(6.5)
$\displaystyle\boldsymbol{Y}_{t}=\boldsymbol{A}(L)\boldsymbol{\eta}_{t}+\boldsymbol{\varepsilon}_{t}\equiv\sum_{\mathsf{u}=1}^{p_{1}}\boldsymbol{A}_{\mathsf{u}}\boldsymbol{\eta}_{t-\mathsf{u}}+\boldsymbol{\varepsilon}_{t},$

Assume that $p_{2}=0$ and $p_{1}=p$, then the model can be formulated as below

(6.6)
$\displaystyle\boldsymbol{Y}_{t}=\sum_{j=1}^{p}\boldsymbol{A}_{j}\boldsymbol{B}_{0}\boldsymbol{Y}_{t-j}+\boldsymbol{\varepsilon}_{t}\equiv\boldsymbol{\mathcal{A}}\boldsymbol{X}_{t-1}+\boldsymbol{\varepsilon}_{t}$

where
$\boldsymbol{X}_{t-1}^{\prime}=\big{(}\boldsymbol{Y}_{t-1}^{\prime}\boldsymbol{B}_{0}^{\prime},...,\boldsymbol{Y}_{t-p}^{\prime}\boldsymbol{B}_{p}^{\prime}\big{)}=\big{(}\boldsymbol{Y}_{t-1}^{\prime},...,\boldsymbol{Y}_{t-p}^{\prime}\big{)}\big{(}\boldsymbol{I}_{p}\otimes\boldsymbol{B}_{0}^{\prime}\big{)}$
and $\boldsymbol{\mathcal{A}}=(\boldsymbol{A}_{1},....,\boldsymbol{A}_{p})$.

In addition, full rank $r$ conditions for these matrices are assumed to hold.
Since the elements of $\boldsymbol{A}_{j}^{\prime}$ and $\boldsymbol{B}_{0}$
are determined only up to nonsingular linear transformations, that is,
$\boldsymbol{A}_{j}\boldsymbol{B}_{0}\equiv\boldsymbol{A}_{j}\boldsymbol{P}^{-1}\boldsymbol{P}\boldsymbol{B}_{0}$,
for $j=1,...,p$ and any $(r\times r)$ nonsingular matrix $\boldsymbol{P}$, we
must impose some normalization conditions to ensure uniqueness of the
parameters. This example can be extended to more complex functional forms
provided that we can obtain observationally equivalent processes.

######  Example 3 (see, Cubadda and Hecq,  (2022)).

Consider the Dimension-Reducible VAR model as

(6.7)
$\displaystyle\boldsymbol{Y}_{t}=\sum_{j=1}^{p}\boldsymbol{A}\alpha_{j}\boldsymbol{A}^{\prime}\boldsymbol{Y}_{t-j}+\boldsymbol{u}_{t},$

Under the assumption that joint DGP of the observed variables
$\boldsymbol{Y}_{t}$ follows that of a FAVAR then the coefficient matrix has
the following structure

(6.8)
$\displaystyle\boldsymbol{A}=\begin{bmatrix}\boldsymbol{I}_{m}&\boldsymbol{0}_{(n-m)\times
m}\\\
\boldsymbol{0}_{m\times(n-m)}&\boldsymbol{B}_{(n-m)\times(r-m)}\end{bmatrix}$

Therefore, to perform structural analysis through the DRVAR, then the
invertability condition for the polynomial VAR coefficient matrix should hold
to obtain the Wold representation of the time series $\boldsymbol{Y}_{t}$.
However, we can invert the polynomial coefficient matrix of
$x_{t}=\sum_{j=1}^{p}\alpha_{j}x_{t-j}+\xi_{t}$ and insert the Wold
representation of the dynamic components $x_{t}$ in expression
$Y_{t}=Ax_{t}+\varepsilon_{t}$ to obtain

(6.9)
$\displaystyle\boldsymbol{Y}_{t}=\boldsymbol{A}\boldsymbol{\gamma}(L)\boldsymbol{\xi}_{t}+\boldsymbol{\varepsilon}_{t},$

where
$\boldsymbol{\gamma}(L)^{-1}=\boldsymbol{I}_{n}-\sum_{j=1}^{p}\alpha_{j}L^{j}$.
Lastly, by linearly projecting $\varepsilon_{t}$ on $\xi_{t}$, we can
decompose the static component as $\varepsilon_{t}=\rho\xi_{t}+v_{t}$, where
$\rho=A_{\perp}A_{\perp}^{\prime}\Sigma_{u}\boldsymbol{A}(\boldsymbol{A}^{\prime}\Sigma_{u}\boldsymbol{A})^{-1}$

(6.10)
$\displaystyle\boldsymbol{Y}_{t}=\underbrace{\boldsymbol{C}(L)\boldsymbol{\xi}_{t}}_{{\color[rgb]{0,0,1}\chi_{t}}}+\boldsymbol{v}_{t},$

where $C_{0}=(A+\rho)$ and $C_{j}=A\gamma_{j}$ for some $j>0$.

The above derivations show how to decompose the underline dynamics of the
observable series $\boldsymbol{Y}_{t}$ into the common component
$\boldsymbol{\chi}_{t}$, and the ignorable errors, $\boldsymbol{v}_{t}$.
Moreover, since the errors $\boldsymbol{\xi}_{t}$ and $\boldsymbol{v}_{t}$ are
uncorrelated at any lead and lags, we can recover the structural shocks solely
by the reduced form errors $\boldsymbol{\xi}_{t}$ of the common component
$\boldsymbol{\chi}_{t}$ using procedures commonly used in structural VAR
analysis. In particular, we can obtain the structural shocks as
$\boldsymbol{u}_{t}=\boldsymbol{C}^{-1}\boldsymbol{D}\boldsymbol{\xi}_{t}$ and
the impulse response functions from
$\boldsymbol{\Psi}(L)=\boldsymbol{C}(L)\boldsymbol{D}^{-1}\boldsymbol{C}$,
where $\boldsymbol{D}$ is the matrix formed by the first $r$ rows of
$\boldsymbol{C}_{0}$ and $\boldsymbol{C}$ is the lower triangular matrix such
that
$\boldsymbol{C}\boldsymbol{C}^{\prime}=\boldsymbol{D}\boldsymbol{A}^{\prime}\boldsymbol{\Sigma}_{u}\boldsymbol{A}\boldsymbol{D}^{\prime}$.
Notice that such identification strategy is based on a unique rotation of the
reduced form common shocks $\boldsymbol{\xi}_{t}$, and hence it does not
require to endow the dynamic component $\boldsymbol{\chi}_{t}$ with an
economic interpretation.

###  Illustrative Examples

######  Example 4 (NVAR$(p,q)$).

The main idea of a NVAR-LUR$(p,q)$ model is as follows. To begin with,
innovation transmission across bilateral links takes place with a lag.
Moreover, uncertainty occurs due to unequal variation in the frequency of
network interactions in comparison to the time series observations. Thus, the
proposed econometric functional form specification accommodates
characteristics on how innovations transmit through the graph (network) over
time.

(6.11) $\displaystyle
x_{t}=\alpha_{1}Ax_{t-1}+...+\alpha_{p}Ax_{t-p}+v_{t}\equiv
A\sum_{j=1}^{p}\theta_{j}x_{t-j}+v_{t},$

######  Proposition 3 (Granger-Causality in NVAR$(p,1)$).

Suppose that $x_{t}$ is generated by (6.11) and assume that $\alpha_{\ell}\neq
0\ \forall\ \ell\in\left\\{1,...,p\right\\}$, $x_{j}$ Granger-causes $x_{i}$
at horizon $h$ iff there exists a connection from $i$ to $j$ of at least one
order $k\in\left\\{k^{\star},k^{\star}+1,...,h\right\\}$, where
$k^{\star}=\mathsf{ceil}(h/p)$.

Consequently, the GIRF has the following form:

(6.12) $\displaystyle\frac{\partial x_{i,t+h}}{\partial
v_{j,t}}\bigg{|}\mathcal{F}_{t}:=\vartheta_{t=k^{\star}}(\alpha)[A_{ij}]|_{t=k^{\star}}+...+\vartheta_{t=k^{\star}}(\alpha)[A_{ij}]|_{t=h}$

The coefficients $\left\\{\vartheta_{t}(\alpha)\right\\}_{t=k^{\star}}^{h}$
are polynomials of $\left\\{\alpha_{j}\right\\}_{j=1}^{p}$. Our aim is to
bridge the interaction between network connectedness and transmission of shock
dynamics under the presence of possibly nonstationarity.

###### Proof.

Consider again the functional form of the NVAR$(p,1)$ model as below:

(6.13) $\displaystyle y_{t}=\alpha_{1}Ay_{t-1}+...+\alpha_{p}Ay_{t-p}+u_{t},$

Let
$\boldsymbol{\alpha}=(\alpha_{1},...,\alpha_{p})^{\prime}\in\mathbb{R}^{p}$
and denote with $\boldsymbol{X}_{t}$ the covariates matrix which includes
information in the lags 1 to $p$ of $y_{t}$ using first-order network
connections,

(6.14) $\displaystyle\boldsymbol{X}_{t}=\big{[}Ay_{t-1},...,Ay_{t-p}\big{]}$

Since the network adjacency matrix is assumed to be fixed (time invariance in
social interactions), then this implies that the dependence of the regressors
on the network dependence can suppressed. This implies that given a known
adjacency matrix $A$, then the unknown parameter vector $\boldsymbol{\alpha}$
can be estimated by OLS, especially due to the absence of exogenous regressors
or nonstationary regressors in the system.

This yields the following optimization problem:

(6.15)
$\displaystyle\underset{\boldsymbol{\alpha}\in\mathbb{R}^{p+1}}{\mathsf{min}}\
\frac{1}{NT}\sum_{t=1}^{T}\big{(}y_{t}-\boldsymbol{X}_{t}\boldsymbol{\alpha}\big{)}^{\prime}\boldsymbol{\Sigma}\big{(}y_{t}-\boldsymbol{X}_{t}\boldsymbol{\alpha}\big{)},$

where $\boldsymbol{\Sigma}=\mathsf{Var}(u_{t})$. Therefore, this implies that

(6.16)
$\displaystyle\hat{\boldsymbol{\alpha}}_{ols}=\left(\sum_{t=1}^{T}\boldsymbol{X}_{t}^{\prime}\boldsymbol{\Sigma}^{-1}\boldsymbol{X}_{t}\right)^{-1}\left(\sum_{t=1}^{T}\boldsymbol{X}_{t}^{\prime}\boldsymbol{\Sigma}^{-1}y_{t}\right)$

Specifically, when we assume that the covariance matrix of the disturbances is
the identity matrix, that is, $\boldsymbol{\Sigma}=\boldsymbol{I}$, then the
OLS estimator takes the form of a pooled OLS estimator as below:

(6.17)
$\displaystyle\hat{\boldsymbol{\alpha}}_{ols}=\left(\sum_{t=1}^{N}\sum_{t=1}^{T}x_{it}x_{it}^{\prime}\right)^{-1}\left(\sum_{t=1}^{N}\sum_{t=1}^{T}x_{it}y_{t}\right)$

As $n\to\infty$, then it holds that

(6.18)
$\displaystyle\sqrt{N}\left(\hat{\boldsymbol{\alpha}}_{ols}-\boldsymbol{\alpha}\right)\Rightarrow\mathcal{N}\left(0,\frac{\sigma^{2}}{T}\mathbb{E}\left[x_{it}x_{it}^{\prime}\right]\right)$

On the other hand, as $T\to\infty$ then the $\hat{\boldsymbol{\alpha}}_{ols}$
estimator is consistent if the model is specified correctly and $y_{t}$ is
ergodic and strictly stationary which implies that

(6.19)
$\displaystyle\sqrt{T}\left(\hat{\boldsymbol{\alpha}}_{ols}-\boldsymbol{\alpha}\right)\Rightarrow\mathcal{N}\big{(}0,\mathbb{E}\left[\boldsymbol{X}_{t}\boldsymbol{X}_{t}^{\prime}\right]^{-1}\mathbb{E}\left[\boldsymbol{X}_{t}^{\prime}\boldsymbol{\Sigma}\boldsymbol{X}_{t}^{\prime}\right]\mathbb{E}\left[\boldsymbol{X}_{t}\boldsymbol{X}_{t}^{\prime}\right]^{-1\prime}\big{)}$

######  Remark 4.

Now, there are certain cases in which although the estimator
$\boldsymbol{\alpha}$ that corresponds to the NVAR$(p,1)$ model can be
obtained via the Expectation-Maximization (EM) algorithm, point identification
of the estimator is not guaranteed. This might suggest that the mapping
between parameters in the process for $\left\\{y_{t}\right\\}_{t=1}^{T}$ and
$\alpha$ is not bijective, similar to the statistical problem of estimating
continuous time models using discrete time data. Furthermore, Bayesian methods
can be employed especially in cases of lack of point identification by
imposing specific distributional assumptions on the Prior and Posterior
distributions corresponding to the solution of the EM algorithm for the model
estimator.

Notice that the process $\left\\{y_{t}\right\\}$ is weakly stationary iff for
all eigenvalues $\lambda_{i}$ of $A$ it holds that $|\lambda_{i}|<1/|a|$.

∎

###  Impulse Response Functions

Assume that $y_{t}$ is stationary, then the long-term response of $y_{t}$ to a
permanent increase in $u_{t}$ is equivalent to the (contemporaneous) response
of $y$ to a disturbance in $\epsilon$, $\partial y/\partial\epsilon$, such
that

(6.20) $\displaystyle R=\underset{h\to\infty}{\mathsf{lim}}\
\left(\frac{\partial y_{t+h}}{\partial u_{t}}+\frac{\partial y_{t+h}}{\partial
u_{t+1}}+....+\frac{\partial y_{t+h}}{\partial
u_{t+h}}\right)\equiv\frac{\partial}{\partial\varepsilon}$

Notice that it holds that $y=\big{(}I-\alpha A\big{)}^{-1}\varepsilon$.
Moreover, since $y_{t}$ is assumed to be a stationary process then it holds
that

(6.21) $\displaystyle R=\underset{h\to\infty}{\mathsf{lim}}\
\sum_{j=0}^{h+1}\frac{\partial y_{t+h}}{\partial
u_{t+h-j}}=\underset{h\to\infty}{\mathsf{lim}}\ \sum_{j=0}^{h+1}\frac{\partial
y_{t+h}}{\partial u_{t}}=\sum_{j=0}^{\infty}\frac{\partial y_{t+h}}{\partial
u_{t}}.$

Therefore, to obtain the impulse response for $x_{t}$, one can write the
econometric specification in the following companion form:

(6.22)
$\displaystyle\boldsymbol{z}_{t}=\boldsymbol{F}\boldsymbol{z}_{t-1}+\boldsymbol{e}_{t},$

where the $(N\times N)$ $\boldsymbol{F}$ matrix is expressed as below:

(6.23) $\displaystyle\boldsymbol{F}=\begin{bmatrix}\alpha_{1}A\ &\
\alpha_{2}A\ &\ldots\ &\ \alpha_{p-1}A\ &\ \alpha_{p}A\\\ I_{N}\ &\ 0_{N}\ &\
\ldots\ &\ 0_{N}\ &\ 0_{N}\ \\\ 0_{N}&I_{N}&\ldots&0_{N}&0_{N}\\\ \vdots\ &\
\vdots\ &\ \ddots\ &\ \vdots\ &\ \vdots\\\ 0_{N}\ &\ 0_{N}\ &\ \ldots\ &\
I_{N}\ &0_{N}\end{bmatrix}$

Therefore, the impulse response of $x_{t}$ to a disturbance in $v_{t}$ is then
given by $(n\times n)$ upper left block in $F^{h}$ such that:

(6.24) $\displaystyle\frac{\partial y_{t+h}}{\partial u_{t}}=\frac{\partial
y_{t+h}}{\partial z_{t+h}}\cdot\frac{\partial z_{t+h}}{\partial
e_{t}}\cdot\frac{\partial e_{t}}{\partial u_{t}}$

SUPPLEMENTARY APPENDIX

"Robust Estimation in Network Vector Autoregression with Nonstationary
Regressors"

Christis Katsouris

University of Southampton $\&$ University of Helsinki

##  Appendix A Supplementary Results

Let $\mathcal{D}$ be the space of càdlàg functions $f:[0,1]\to\mathbb{R}$
equipped with the Skorokhod topology. As a measurable structure on
$\mathcal{D}$ we consider the corresponding Borel $\sigma-$algebra
$\mathcal{B}(D)$. We write $X_{n}\overset{\mathcal{D}}{\to}X$ whenever
$X_{n}$, $X$ are random variables taking values in
$\big{(}\mathcal{D},\mathcal{B}(\mathcal{D})\big{)}$ such that $X_{n}$
converges weakly to $X$ in $\mathcal{D}$ as $n\to\infty$, that is,
$\underset{n\to\infty}{\mathsf{lim}}\mathbb{E}f(X_{n})\to\mathbb{E}f(X)$ for
all bounded continuous functions $f:\mathcal{D}\to\mathbb{R}$.

###  Matrix Representation of multivariate LUR process

Consider that the $p-$dimensional vector of the system’s predictors
$\boldsymbol{X}_{t}$ is generated via the following multivariate LUR process
(see, also Katsouris, 2023d )

(A.1)
$\displaystyle\boldsymbol{X}_{t}=\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}_{p}}{T^{\upgamma_{x}}}\right)\boldsymbol{X}_{t-1}+\boldsymbol{u}_{t},\
\ \ \ t\in\left\\{1,...,T\right\\}$

Let
$\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x}):=\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}_{p}}{T^{\upgamma_{x}}}\right)$,
then under the assumption of homogeneous coefficients of persistence, that is,
$c_{i}\equiv c$ for all $i\in\left\\{1,...,p\right\\}$, then (A.1) can be
written as below

(A.2)
$\displaystyle\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})\boldsymbol{X}=\boldsymbol{U}$

where the non-stochastic $T\times T$ matrix
$\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})$ takes the bi-diagonal form as
below

(A.3)
$\displaystyle\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})=\begin{bmatrix}{\color[rgb]{1,0,0}1}&0&0&\ldots&0&0\\\
{\color[rgb]{0,0,1}-\varphi_{T}(c,\upgamma_{x})}&{\color[rgb]{1,0,0}1}&0&\ldots&0&0\\\
0&{\color[rgb]{0,0,1}-\varphi_{T}(c,\upgamma_{x})}&{\color[rgb]{1,0,0}1}&\vdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&0&{\color[rgb]{0,0,1}-\varphi_{T}(c,\upgamma_{x})}&{\color[rgb]{1,0,0}1}\end{bmatrix}$

with
${\color[rgb]{0,0,1}-\varphi_{T}(c_{i},\upgamma_{x})=(1-\frac{c_{i}}{T^{\upgamma_{x}}})}$.

Denote with
$\boldsymbol{K}_{T}=\boldsymbol{I}_{T}-\frac{1}{T}\mathbf{1}\mathbf{1}^{\top}$
denote the $T\times T$ centering matrix and
$\boldsymbol{X}_{K}=\boldsymbol{K}_{T}\boldsymbol{X}$ the corresponding
centred big data matrix. Therefore, the corresponding centered covariance
matrix can be written in the following form

(A.4)
$\displaystyle\boldsymbol{\mathcal{S}}_{T}=\frac{1}{p}\boldsymbol{X}_{K}\boldsymbol{X}_{K}^{\top}.$

######  Remark 5.

  * •

Firstly, notice that the centering matrix $\boldsymbol{K}_{T}$ is employed for
demeaning the original data in high dimensional settings. One application can
be found in the representation of the dynamic time series model (see, Handbook
of Econometrics). A different implementation of this property is employed by
Maillet et al.,  (2015) as well as by Davis et al.,  (2016) who study the
asymptotic theory for the sample covariance matrix of a heavy-tailed
multivariate matrix. Another application is presented by Phillips and
Magdalinos,  (2013) who derive an expression of the condition number in
explosive VAR models.

  * •

Secondly, notice that the covariance matrix representation given by expression
(A.4) becomes a functional of the nuisance coefficient of persistence, which
is a data specific feature. However considering the structure of our proposed
novel VaR-CoVaR risk matrix, the asymptotic theory of its largest eigenvalues
will be more challenging than conventional approaches. In this case, we will
need to consider both the framework proposed by Davis et al.,  (2016) as well
as the proposed regression-based estimation methodology for the risk matrix,
the quantile predictive regressions which capture the possibly nonstationary
nature of predictors; in order to derive an analytical tractable limiting
distribution with desirable properties for statistical inference.

Moreover, we have that
$\boldsymbol{X}=\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})^{-1}\boldsymbol{U}$
which implies that the covariance type matrix $\boldsymbol{\mathcal{S}}_{T}$
can be written in the following form

(A.5)
$\displaystyle\boldsymbol{\mathcal{S}}_{T}=\frac{1}{p}\left(\boldsymbol{K}_{T}\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})^{-1}\boldsymbol{U}\right)\left(\boldsymbol{K}_{T}\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})^{-1}\boldsymbol{U}\right)^{\top}$

Then, the uncentered matrix can be written as below

(A.6)
$\displaystyle\boldsymbol{\mathcal{S}}^{o}_{T}=\frac{1}{p}\boldsymbol{X}\boldsymbol{X}^{\top}.$

and using the multivariate LUR process representation (A.2) we obtain

(A.7)
$\displaystyle\boldsymbol{\mathcal{S}}^{o}_{T}=\left(\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})^{-1}\boldsymbol{U}\right)\left(\boldsymbol{\Phi}_{T}(c_{i},\upgamma_{x})^{-1}\boldsymbol{U}\right)^{\top}$

Denote with $\hat{\lambda}_{k}$ and $\hat{\lambda}^{o}_{k}$ the $k-$th largest
eigenvalues of $\boldsymbol{\mathcal{S}}_{T}$ and
$\boldsymbol{\mathcal{S}}^{o}_{T}$ respectively. Then, one would be interested
to examine the asymptotic behaviour of the largest eigenvalues of the matrices
$\boldsymbol{\mathcal{S}}_{T}$ and $\boldsymbol{\mathcal{S}}^{o}_{T}$ in a
high-dimensional setting.

###  Predictive Regression with variable addition

Consider the model

(A.8) $\displaystyle y_{t}$ $\displaystyle=\beta x_{t-1}+u_{t}$ (A.9)
$\displaystyle x_{t}$ $\displaystyle=\rho x_{t-1}+v_{t}$

with the IVX instrument defined as

(A.10) $\displaystyle z_{t-1}=(1-\rho_{z}L)^{-1}\Delta x_{t-1},\ \
\rho_{z}=\left(1-\frac{c_{z}}{T^{\lambda}}\right)$

For the IVX estimator it holds that
$\hat{\beta}_{ivx}-\beta=\mathcal{O}_{p}\left(T^{-\frac{\lambda+1}{2}}\right)$.
Moreover, the limiting distribution of the IVX-based t type statistic
$t_{ivx}$ for the null $\beta=0$, under the sequences of local alternatives of
the type $\beta=\beta_{T}=bT^{-\frac{\lambda+1}{2}}$, is shown to be

(A.11) $\displaystyle
t_{ivx}\Rightarrow\mathcal{Z}+b\frac{\sigma_{v}\sqrt{2}}{\sigma_{u}\sqrt{c_{z}}}\left(B_{c}^{2}(1)-\int_{0}^{1}B_{c}(s)dB_{c}(s)\right)$

where $dB_{c}(s):=-cB_{c}(s)ds+dV(s)$ and $\mathcal{Z}$ is a standard normal
variate independent of $V(t)$ and $B_{c}(s)$.

######  Proposition 4.

Denote by $\hat{\beta}_{ivx}^{*}$ the IVX estimator of $\beta$ with $y_{t-1}$
instrumented by itself and by $t_{ivx}^{*}$ the corresponding t statistic. Let
$\sigma_{uv}\neq 0$, then we have that

Under the assumptions of the proposition the following hold

(A.12)
$\displaystyle\frac{1}{T^{\frac{1+\lambda}{2}}}\sum_{t=2}^{T}z_{t-1}u_{t}$
$\displaystyle\Rightarrow\mathcal{N}\left(0,\frac{\sigma_{u}^{2}\sigma_{v}^{2}}{2c_{z}}\right)$
(A.13) $\displaystyle\frac{1}{T^{1+\lambda}}\sum_{t=2}^{T}z^{2}_{t-1}$
$\displaystyle\to\frac{\sigma_{v}^{2}}{2c_{z}}$ (A.14)
$\displaystyle\frac{1}{T^{\frac{1+\lambda}{2}}}\sum_{t=2}^{T}z_{t-1}x_{t-1}$
$\displaystyle\Rightarrow\frac{\sigma_{v}^{2}}{c_{z}}\left(B_{c}^{2}(1)-\int_{0}^{1}B_{c}(s)dB_{c}(s)\right),$

###  Limit Theory for Moderately Explosive Systems

##### Case II: C with not distinct diagonal elements, $c_{i}=c_{j}$ for some
$i\neq j$

We now consider case (II), where the matrix $C$ does not have distinct
diagonal elements. In this case two or more elements of $x_{t}$ have
comparable moderately explosive behaviour governed by a common autoregressive
root of the form $\rho_{nj}=1+c_{j}/n^{\alpha}$. We know that under such
conditions the second moment matrix $\sum_{t=1}^{n}x_{t}x_{t}^{\prime}$ in the
regression model is asymptotically singular, which is explained by the fact
that some elements of $x_{t}$ have common moderately explosive behaviour.
These elements of $x_{t}$ are then asymptotically multicollinear in much the
same way as regressors that are cointegrated or have common deterministic
trends. To deal with this singularity in the regression model, we can develop
an asymptotic theory for the regression in a similar way by rotating the
regression coordinates in the direction of the common explosive behaviour and
in an orthogonal direction. Here, however, the rotation is a random process
determined by the regressor vector $x_{n}$. The randomness and the sample size
dependence in the rotation present further complications in the development of
the asymptotics because the limit theory segementation then depends on weak
convergence of the rotation matrix. In what follows, we develop the related
asymptotic theory for the special case where $C$ is the scalar matrix
$C=cI_{K}$ and $\rho_{N}=1+c/n^{\alpha}$.

For $\left(k_{n}\right)_{n\in\mathbb{N}}$ any sequence increasing to infinity,
we define

(A.15) $\displaystyle
Y_{C_{n}}:=\frac{1}{n^{\alpha/2}}\sum_{j=1}^{k_{n}}R_{n}^{-j}F_{x}(1)\epsilon_{j}.$

The stochastic sequence of $Y_{C_{n}}$ plays an important role in determining
the asymptotic behaviour of the explosive systems.

##### Proof of Lemma 4.2 MP (2009)

######  Lemma 2.

For each $\alpha\in(0,1)$, $C>0$ and a sequence $k_{n}$ satisfying
$\left\lVert R_{n}\right\rVert^{-k_{n}}\to 0$ and $n^{\alpha}\left\lVert
R_{n}\right\rVert^{-(n-k_{n})}\to 0$, we have that, as $n\to\infty$,

  * (i)

(A.16) $\displaystyle
n^{-\alpha}\mathbb{E}\left\lVert\text{vec}\sum_{t=1}^{n}u_{0t}\left(\sum_{j=t+1}^{n}R_{n}^{t-j}u_{xj}\right)^{\prime}R_{n}^{-n}\right\rVert\to
0,$

  * (ii)

(A.17) $\displaystyle
n^{-\alpha}\mathbb{E}\left\lVert\text{vec}\sum_{t=1}^{k_{n}}u_{0t}\left(\sum_{j=t+1}^{n}R_{n}^{t-j}u_{xj}\right)^{\prime}R_{n}^{-n}\right\rVert\to
0.$

###### Proof.

Notice that the sample covariance can be written as below

$\displaystyle\text{vec}\frac{1}{n^{\alpha}}\sum_{t=1}^{n}u_{0t}x_{t}^{\prime}R_{n}^{-n}$
$\displaystyle=\frac{1}{n^{\alpha}}\text{vec}\left\\{\sum_{t=1}^{n}u_{0t}\left(\sum_{j=1}^{t}R_{n}^{t-j}u_{xj}\right)^{\prime}R_{n}^{-n}\right\\}+o_{p}(1)$
$\displaystyle=\frac{1}{n^{\alpha}}\text{vec}\left\\{\sum_{t=k_{n}+1}^{n}u_{0t}\left(\sum_{j=1}^{n}R_{n}^{t-j}u_{xj}\right)^{\prime}R_{n}^{-n}\right\\}+o_{p}(1)$
$\displaystyle=\frac{1}{n^{\alpha/2}}\sum_{t=k_{n}+1}^{n}\left(R_{n}^{-(n-t)}\otimes
u_{0t}\right)\left(\frac{1}{n^{\alpha/2}}\sum_{j=1}^{n}R_{n}^{t-j}u_{xj}\right)^{\prime}R_{n}^{-n}+o_{p}(1)$
$\displaystyle=\big{[}I_{mK}+o_{p}(1)\big{]}\frac{1}{n^{\alpha/2}}\sum_{t=k_{n}+1}^{n}\left(R_{n}^{-(n-t)}\otimes
u_{0t}\right)Y_{C_{n}}+o_{p}(1).$

Applying the BN decomposition on the proceeding expression and making use of

(A.18)
$\displaystyle\frac{1}{n^{\alpha/2}}\sum_{t=k_{n}+1}^{n}\left(R_{n}^{-(n-t)}\otimes\Delta\tilde{\epsilon}_{0t}\right)=o_{p}(1),\
\ \text{as}\ \ n\to\infty,$

we obtain the following expression for the sample covariance as $n\to\infty$

vec
$\displaystyle\left\\{\frac{1}{n^{\alpha}}\sum_{t=1}^{n}u_{0t}x_{t}^{\prime}R_{n}^{-n}\right\\}$
$\displaystyle=\big{[}I_{mK}+o_{p}(1)\big{]}\frac{1}{n^{\alpha/2}}\sum_{t=1}^{n-k_{n}}\left[R_{n}^{-(n-k_{n}-t)}Y_{C_{n}}\otimes
F_{0}(1)\epsilon_{t+\kappa_{n}}\right].$

Furthermore, letting
$\mathcal{F}_{n,i}:=\sigma\left(x_{0},\epsilon_{i},\epsilon_{i-1},...\right)$
since $Y_{C_{n}}$ is $\mathcal{F}_{n,k_{n}}$ measurable, then

(A.19)
$\displaystyle\xi_{n,t+k+n}:=\frac{1}{n^{\alpha/2}}\left[R_{n}^{-(n-k_{n}-t)}Y_{C_{n}}\otimes
F_{0}(1)\epsilon_{t+\kappa_{n}}\right].$

is an $\mathbb{R}^{mK}-$valued martingale difference array with respect to
$\mathcal{F}_{n,t+k_{n}}$. Moreover, we denote by
$M_{n,k}:=\sum_{t=1}^{k}\xi_{n,t+k_{n}}$ the martingale array corresponding to
$\xi_{n,t+k_{n}}$ and by $\langle
M_{n,k}\rangle:=\sum_{t=1}^{k}\mathbb{E}_{\mathcal{F}_{n,t+k_{n}}}\left[\xi_{n,t+k_{n}}\xi_{n,t+k_{n}}^{\prime}\right]$
the predictable quadratic variation of $M_{n,k}$. Since we have that

$\displaystyle\langle M_{n,k}\rangle_{n-k_{n}}$
$\displaystyle=\left[\frac{1}{n^{\alpha}}\sum_{t=1}^{n-k_{n}}R_{n}^{-(n-k_{n}-t)}Y_{C_{n}}Y_{C_{n}}^{\prime}R_{n}^{-(n-k_{n}-t)}\right]\otimes\Omega_{00}$
$\displaystyle=\big{[}I_{mK}+o_{p}(1)\big{]}\left(\int_{0}^{\infty}e^{-pC}Y_{C_{n}}Y_{C_{n}}^{\prime}e^{-pC}\otimes\Omega_{00}\right)$

and $\langle M_{n,k}\rangle_{n-k_{n}}$ is $\mathcal{F}_{n,k_{n}}-$measurable
with $\mathcal{F}_{n,k_{n}}\subset\mathcal{F}_{n,k_{n}+1}$, that is, the
$\sigma-$algebra supporting the predictable quadratic variation is smaller
than each of the elements of the filtration supporting the array
$\xi_{n,t+k_{n}}$.

Therefore, we obtain that

$\displaystyle\text{ vec
}\left\\{\frac{1}{n^{\alpha}}\sum_{t=1}^{n}u_{0t}x_{t}^{\prime}R_{n}^{-n}\right\\}$
$\displaystyle=M_{n,n-k_{n}}+o_{p}(1)$
$\displaystyle\Rightarrow\mathcal{MN}\left(0,\int_{0}^{\infty}e^{pC}Y_{C}Y_{C}^{\prime}e^{pC}\otimes\Omega_{00}\right)$

as $n\to\infty$. The verification of the Lindeberg and tightness conditions
needed for this result is provided by Proposition A1 of MP.

Therefore, we can now deduce the limit behavior of the regression coefficient
for the case (I), where the localizing coefficients are distinct. This is
given by

$\displaystyle\text{vec}\left\\{n^{\alpha}\left(\hat{A}_{n}-A\right)R_{n}^{n}\right\\}$
$\displaystyle=\left[\left(\frac{1}{n^{2\alpha}}\sum_{t=1}^{n}R_{n}^{-n}x_{t}x_{t}^{\prime}\right)^{-1}\otimes
I_{m}\right]\text{vec}\left(\frac{1}{n^{\alpha}}\sum_{t=1}^{n}u_{0t}x_{t}^{\prime}R_{n}^{-n}\right)$
$\displaystyle=\big{[}I_{mK}+o_{p}(1)\big{]}\left[\left(\int_{0}^{\infty}e^{-pC}Y_{C_{n}}Y_{C_{n}}^{\prime}e^{-pC}dp\right)^{-1}\otimes
I_{m}\right]\sum_{t=1}^{n-k_{n}}\xi_{n,t}$
$\displaystyle=\big{[}I_{mK}+o_{p}(1)\big{]}\left[\left(\int_{0}^{\infty}e^{-pC}Y_{C_{n}}Y_{C_{n}}^{\prime}e^{-pC}dp\right)^{-1}\otimes\Omega_{00}^{-1}\right]\left(I_{mK}\otimes\Omega_{00}\right)\sum_{t=1}^{n-k_{n}}\xi_{n,t}$
$\displaystyle=\big{[}I_{mK}+o_{p}(1)\big{]}\langle
M_{n}\rangle_{n-k_{n}}^{-1}\left(I_{mK}\otimes\Omega_{00}\right)M_{n,n-k_{n}}.$

The limiting distribution of $M_{n,n-k_{n}}$ is established in (26). We also
show in the Appendix that $M_{n,n-k_{n}}$ satisfies the requirements of
Proposition A1 (iii), so that joint convergence of $M_{n,n-k_{n}}$ and
$\langle M\rangle_{n-k_{n}}$ applies. Therefore, with the next theorem we
prove that the regression coefficient has a mixed normal limiting
distribution. ∎

######  Theorem 9.

For the model (1)-(2) with $R_{n}=I_{K}+C/n^{\alpha}$, $c_{i}>0$ for all
$i,c_{i}\neq c_{j}$ for all $i\neq j,\alpha\in(0,1)$, and weakly dependent
errors satisfying Assumption LP, we have that

(A.20) $\displaystyle
n^{\alpha}\left(\hat{A}_{n}-A\right)R_{n}^{n}\Rightarrow\mathcal{MN}\left(0,\left(\int_{0}^{\infty}e^{-pC}Y_{C_{n}}Y_{C_{n}}^{\prime}e^{-pC}dp\right)^{-1}\otimes\Omega_{00}\right)$

###### Proof.

Start by defining the orthogonal random matrix $H_{n}=[H_{cn},H_{\perp n}]$,
where

(A.21) $\displaystyle
H_{cn}=\frac{x_{n}}{\left(x_{n}^{\prime}x_{n}\right)^{1/2}},\ \ H_{\perp
n}^{\prime}H_{cn}=0\ \ \text{almost surely},$

and the $K\times(K-1)$ random matrix $H_{\perp n}$ is an orthogonal complement
to $H_{cn}$ satisfying $H_{\perp n}^{\prime}H_{\perp n}=I_{K-1}$ and $H_{\perp
n}^{\prime}H_{\perp n}=I_{K}-H_{cn}^{\prime}H_{cn}$ almost surely. Therefore,
the asymptotic behaviour of $H_{\perp n}H_{\perp n}^{\prime}$ is given by

(A.22) $\displaystyle H_{\perp n}H_{\perp
n}^{\prime}=I_{K}-\frac{Y_{cn}Y_{cn}^{\prime}}{Y_{cn}^{\prime}Y_{cn}}+o_{p}(1)\Rightarrow
I_{K}-\frac{Y_{c}Y_{c}^{\prime}}{Y_{c}^{\prime}Y_{c}}\ \ \text{as}\
n\to\infty,$

where $Y_{cn}$ and $Y_{c}$ are random vectors $Y_{C_{n}}$ and $Y_{C}$ of Lemma
4.1 with $C=cI_{K}$. Then, applying this rotation to the moderately explosive
regressor vector yields

(A.23) $\displaystyle
z_{t}=H_{n}^{\prime}x_{t}=\begin{bmatrix}H_{cn}^{\prime}x_{t}\\\ H_{\perp
n}^{\prime}x_{t}\end{bmatrix}:=\begin{bmatrix}z_{1t}\\\ z_{2t}\end{bmatrix}$

with $z_{2t}$ satisfying the reverse autoregression
$z_{2t}=\rho_{n}^{-1}z_{2t+1}-\rho_{n}^{-1}H_{\perp n}^{\prime}u_{xt}$, which
gives rise to the following

(A.24) $\displaystyle z_{2t}=-H_{\perp
n}^{\prime}\sum_{j=1}^{n-t}\rho_{n}^{-j}u_{xt+j},$

because $z_{2n}=H_{\perp n}^{\prime}x_{n}=0$. Using the orthogonality
condition of $H_{n}$, we obtain the following expression for the least squares
estimator after rotation of the regression space

$\displaystyle n^{\frac{1+\alpha}{2}}\left(\hat{A}_{n}-A\right)$
$\displaystyle=\left(\frac{1}{n^{(1+\alpha)/2}}\sum_{t=1}^{n}u_{0t}z_{t}^{\prime}\right)\left(\frac{1}{n^{1+\alpha}}\sum_{t=1}^{n}z_{t}z_{t}^{\prime}\right)^{-1}H_{n}^{\prime}$
$\displaystyle=\left(\frac{1}{n^{(1+\alpha)/2}}\sum_{t=1}^{n}u_{0t}z_{t}^{\prime}\right)$

∎

## References

  * Adamek et al., (2022) Adamek, R., Smeekes, S., and Wilms, I. (2022).  Local projection inference in high dimensions.  arXiv preprint arXiv:2209.03218.
  * Adrian and Brunnermeier, (2016) Adrian, T. and Brunnermeier, M. K. (2016).  Covar.  The American Economic Review, 106(7):1705.
  * Agosto et al., (2016) Agosto, A., Cavaliere, G., Kristensen, D., and Rahbek, A. (2016).  Modeling corporate defaults: Poisson autoregressions with exogenous covariates (parx).  Journal of Empirical Finance, 38:640–663.
  * Ahn and Reinsel, (1990) Ahn, S. K. and Reinsel, G. C. (1990).  Estimation for partially nonstationary multivariate autoregressive models.  Journal of the American statistical association, 85(411):813–823.
  * Andersen and Varneskov, (2021) Andersen, T. G. and Varneskov, R. T. (2021).  Consistent inference for predictive regressions in persistent economic systems.  Journal of Econometrics, 224(1):215–244.
  * Anufriev and Panchenko, (2015) Anufriev, M. and Panchenko, V. (2015).  Connecting the dots: Econometric methods for uncovering networks with an application to the australian financial institutions.  Journal of Banking & Finance, 61:S241–S255.
  * Armillotta and Fokianos, (2022) Armillotta, M. and Fokianos, K. (2022).  Nonlinear network autoregression.  arXiv preprint arXiv:2202.03852.
  * Badev, (2021) Badev, A. (2021).  Nash equilibria on (un) stable networks.  Econometrica, 89(3):1179–1206.
  * Barigozzi et al., (2023) Barigozzi, M., Cho, H., and Owens, D. (2023).  Fnets: Factor-adjusted network estimation and forecasting for high-dimensional time series.  Journal of Business & Economic Statistics, pages 1–13.
  * Baruník and Křehlík, (2018) Baruník, J. and Křehlík, T. (2018).  Measuring the frequency dynamics of financial connectedness and systemic risk.  Journal of Financial Econometrics, 16(2):271–296.
  * Basu and Subba Rao, (2023) Basu, S. and Subba Rao, S. (2023).  Graphical models for nonstationary time series.  The Annals of Statistics, 51(4):1453–1483.
  * Breitung and Demetrescu, (2015) Breitung, J. and Demetrescu, M. (2015).  Instrumental variable and variable addition based inference in predictive regressions.  Journal of Econometrics, 187(1):358–375.
  * Bykhovskaya, (2022) Bykhovskaya, A. (2022).  Time series approach to the evolution of networks: Prediction and estimation.  Journal of Business & Economic Statistics, 41(1):170–183.
  * Cavanagh et al., (1995) Cavanagh, C. L., Elliott, G., and Stock, J. H. (1995).  Inference in models with nearly integrated regressors.  Econometric theory, 11(5):1131–1147.
  * Chen et al., (2019) Chen, C. Y.-H., Härdle, W. K., and Okhrin, Y. (2019).  Tail event driven networks of sifis.  Journal of Econometrics, 208(1):282–298.
  * Chen et al., (2023) Chen, E. Y., Fan, J., and Zhu, X. (2023).  Community network auto-regression for high-dimensional time series.  Journal of Econometrics, 235(2):1239–1256.
  * Cho et al., (2023) Cho, H., Maeng, H., Eckley, I. A., and Fearnhead, P. (2023).  High-dimensional time series segmentation via factor-adjusted vector autoregressive modeling.  Journal of the American Statistical Association, pages 1–13.
  * Cubadda and Hecq, (2022) Cubadda, G. and Hecq, A. (2022).  Dimension reduction for high-dimensional vector autoregressive models.  Oxford Bulletin of Economics and Statistics, 84(5):1123–1152.
  * Daouia et al., (2022) Daouia, A., Gijbels, I., and Stupfler, G. (2022).  Extremile regression.  Journal of the American Statistical Association, 117(539):1579–1586.
  * Davis et al., (2016) Davis, R. A., Mikosch, T., and Pfaffel, O. (2016).  Asymptotic theory for the sample covariance matrix of a heavy-tailed multivariate time series.  Stochastic Processes and their Applications, 126(3):767–799.
  * Diebold and Yılmaz, (2014) Diebold, F. X. and Yılmaz, K. (2014).  On the network topology of variance decompositions: Measuring the connectedness of financial firms.  Journal of Econometrics, 182(1):119–134.
  * Dou and Müller, (2021) Dou, L. and Müller, U. K. (2021).  Generalized local-to-unity models.  Econometrica, 89(4):1825–1854.
  * Doukhan et al., (2023) Doukhan, P., Neumann, M. H., and Truquet, L. (2023).  Stationarity and ergodic properties for some observation-driven models in random environments.  The Annals of Applied Probability, 33(6B):5145–5170.
  * Fan et al., (2023) Fan, Y., Han, F., and Park, H. (2023).  Estimation and inference in a high-dimensional semiparametric gaussian copula vector autoregressive model.  Journal of Econometrics, 237(1):105513.
  * Fang et al., (2023) Fang, P., Gao, Z., and Tsay, R. S. (2023).  Determination of the effective cointegration rank in high-dimensional time-series predictive regressions.  arXiv preprint arXiv:2304.12134.
  * Fort and Roberts, (2005) Fort, G. and Roberts, G. O. (2005).  Subgeometric ergodicity of strong markov processes.  Annals of Applied Probability, 15(2):1565–1589.
  * Gobet and Matulewicz, (2017) Gobet, E. and Matulewicz, G. (2017).  Parameter estimation of ornstein–uhlenbeck process generating a stochastic graph.  Statistical Inference for Stochastic Processes, 20:211–235.
  * Härdle et al., (2016) Härdle, W. K., Wang, W., and Yu, L. (2016).  Tenet: Tail-event driven network risk.  Journal of Econometrics, 192(2):499–513.
  * He and Song, (2018) He, X. and Song, K. (2018).  Measuring diffusion over a large network.  arXiv preprint arXiv:1812.04195.
  * Holberg and Ditlevsen, (2023) Holberg, C. and Ditlevsen, S. (2023).  Uniform inference for cointegrated vector autoregressive processes.  arXiv preprint arXiv:2306.03632.
  * Huang et al., (2020) Huang, D., Wang, F., Zhu, X., and Wang, H. (2020).  Two-mode network autoregressive model for large-scale networks.  Journal of Econometrics, 216(1):203–219.
  * Jansson and Moreira, (2006) Jansson, M. and Moreira, M. J. (2006).  Optimal inference in regression models with nearly integrated regressors.  Econometrica, 74(3):681–714.
  * Kapetanios et al., (2014) Kapetanios, G., Mitchell, J., and Shin, Y. (2014).  A nonlinear panel data model of cross-sectional dependence.  Journal of Econometrics, 179(2):134–157.
  * Kasparis et al., (2015) Kasparis, I., Andreou, E., and Phillips, P. C. B. (2015).  Nonparametric predictive regression.  Journal of Econometrics, 185(2):468–494.
  * Katsouris, (2021) Katsouris, C. (2021).  Optimal portfolio choice and stock centrality for tail risk events.  arXiv preprint arXiv:2112.12031.
  * (36) Katsouris, C. (2023a).  Estimating conditional value-at-risk with nonstationary quantile predictive regression models.  arXiv preprint arXiv:2311.08218.
  * (37) Katsouris, C. (2023b).  Limit theory under network dependence and nonstationarity.  arXiv preprint arXiv:2308.01418.
  * (38) Katsouris, C. (2023c).  Optimal estimation methodologies for panel data regression models.  arXiv preprint arXiv:2311.03471.
  * (39) Katsouris, C. (2023d).  Statistical estimation for covariance structures with tail estimates using nodewise quantile predictive regression models.  arXiv preprint arXiv:2305.11282.
  * Keeling and Eames, (2005) Keeling, M. J. and Eames, K. T. (2005).  Networks and epidemic models.  Journal of the royal society interface, 2(4):295–307.
  * Kojevnikov et al., (2021) Kojevnikov, D., Marmer, V., and Song, K. (2021).  Limit theorems for network dependent random variables.  Journal of Econometrics, 222(2):882–908.
  * Kostakis et al., (2015) Kostakis, A., Magdalinos, T., and Stamatogiannis, M. P. (2015).  Robust econometric inference for stock return predictability.  The Review of Financial Studies, 28(5):1506–1553.
  * Kostakis et al., (2018) Kostakis, A., Magdalinos, T., and Stamatogiannis, M. P. (2018).  Taking stock of long-horizon predictability tests: Are factor returns predictable?  Available at SSRN 3284149.
  * Krampe et al., (2023) Krampe, J., Paparoditis, E., and Trenkler, C. (2023).  Structural inference in sparse high-dimensional vector autoregressions.  Journal of Econometrics, 234(1):276–300.
  * Laurent and Shi, (2022) Laurent, S. and Shi, S. (2022).  Unit root test with high-frequency data.  Econometric Theory, 38(1):113–171.
  * Lee, (2016) Lee, J. H. (2016).  Predictive quantile regression with persistent covariates: Ivx-qr approach.  Journal of Econometrics, 192(1):105–118.
  * Lee and Song, (2019) Lee, J. H. and Song, K. (2019).  Stable limit theorems for empirical processes under conditional neighborhood dependence.  Bernoulli, 25(2):1189–1224.
  * Liu and Phillips, (2023) Liu, Y. and Phillips, P. C. (2023).  Robust inference with stochastic local unit root regressors in predictive regressions.  Journal of Econometrics, 235(2):563–591.
  * Magdalinos, (2021) Magdalinos, T. (2021).  Least squares and ivx limit theory in systems of predictive regressions with garch innovations.  Econometric Theory, pages 1–38.
  * (50) Magdalinos, T. and Phillips, P. C. B. (2009a).  Limit theory for cointegrated systems with moderately integrated and moderately explosive regressors.  Econometric Theory, 25(2):482–526.
  * (51) Magdalinos, T. and Phillips, P. C. B. (2009b).  Limit theory for cointegrated systems with moderately integrated and moderately explosive regressors.  Econometric Theory, 25(2):482–526.
  * Maillet et al., (2015) Maillet, B., Tokpavi, S., and Vaucher, B. (2015).  Global minimum variance portfolio optimisation under some model risk: A robust regression-based approach.  European Journal of Operational Research, 244(1):289–299.
  * Matsui and Pedersen, (2022) Matsui, M. and Pedersen, R. S. (2022).  Characterization of the tail behavior of a class of bekk processes: A stochastic recurrence equation approach.  Econometric Theory, 38(1):1–34.
  * (54) Meitz, M. and Saikkonen, P. (2008a).  Ergodicity, mixing, and existence of moments of a class of markov models with applications to garch and acd models.  Econometric Theory, 24(5):1291–1320.
  * (55) Meitz, M. and Saikkonen, P. (2008b).  Stability of nonlinear ar-garch models.  Journal of Time Series Analysis, 29(3):453–475.
  * Meitz and Saikkonen, (2021) Meitz, M. and Saikkonen, P. (2021).  Subgeometric ergodicity and $\beta$-mixing.  Journal of Applied Probability, 58(3):594–608.
  * Meitz and Saikkonen, (2022) Meitz, M. and Saikkonen, P. (2022).  Subgeometrically ergodic autoregressions with autoregressive conditional heteroskedasticity.  arXiv preprint arXiv:2205.11953.
  * Menzel, (2021) Menzel, K. (2021).  Bootstrap with cluster-dependence in two or more dimensions.  Econometrica, 89(5):2143–2188.
  * Mikusheva, (2007) Mikusheva, A. (2007).  Uniform inference in autoregressive models.  Econometrica, 75(5):1411–1452.
  * Mitchener and Richardson, (2019) Mitchener, K. J. and Richardson, G. (2019).  Network contagion and interbank amplification during the great depression.  Journal of Political Economy, 127(2):000–000.
  * Nicholson et al., (2017) Nicholson, W. B., Matteson, D. S., and Bien, J. (2017).  Varx-l: Structured regularization for large vector autoregressions with exogenous variables.  International Journal of Forecasting, 33(3):627–651.
  * Olmo and Sanso-Navarro, (2023) Olmo, J. and Sanso-Navarro, M. (2023).  A nonparametric spatial regression model using partitioning estimators.  Econometrics and Statistics.
  * Paruolo, (1997) Paruolo, P. (1997).  Asymptotic inference on the moving average impact matrix in cointegrated 1 (1) var systems.  Econometric Theory, 13(1):79–118.
  * (64) Phillips, P. C. B. (1987a).  Time series regression with a unit root.  Econometrica: Journal of the Econometric Society, pages 277–301.
  * (65) Phillips, P. C. B. (1987b).  Towards a unified asymptotic theory for autoregression.  Biometrika, 74(3):535–547.
  * Phillips, (2014) Phillips, P. C. B. (2014).  On confidence intervals for autoregressive roots and predictive regression.  Econometrica, 82(3):1177–1195.