EconBase
← Back to paper

Robust Estimation in Network Vector Autoregression with Nonstationary Regressors

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

111,227 characters · 31 sections · 105 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Robust Estimation in Network Vector Autoregression with Nonstationary Regressors

\pagenumbering{roman}

abstractThis 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

\setcounter{page}{1} \pagenumbering{arabic}

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., he2018measuring), the financial contagion in stock markets (e.g., hardle2016tenet, chen2019tail, mitchener2019network), the detection of market exuberance (e.g., Magdal2009limit) as well as the spread of epidemic diseases (e.g., keeling2005networks). 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., zhu2017network, zhu2019network and zhu2020grouped). 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 dependence\footnote{Regardless 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 adamek2022local, chen2023community, krampe2023structural and zhang2023statistical as well as barigozzi2023fnets, cho2023high and fang2023determination.} (see, also basu2023graphical). Recently, bykhovskaya2022time 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 schennach2018long) 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 zhu2017network). 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 PM2009econometric, 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, kostakis2015Robust, kostakis2018taking. 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 magdalinos2021least. In particular, within the predictive regression literature the regressand corresponds to a continuous stationary covariate such as stock returns while regressors correspond to persistent data\footnote{ 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 agosto2016modeling) and armillotta2022nonlinear 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 jansson2006optimal, phillips2013predictive, breitung2015instrumental, lee2016predictive, kostakis2015Robust, kasparis2015nonparametric, andersen2021consistent and liu2023robust among others) and systems of predictive and cointegrated regressors (see Phillips2008limit, magdalinos2009limit and magdalinos2021least), as well as limit results for the case of network vector autoregressive models under time series stationarity (see zhu2017network) and inference techniques under cluster dependence in two or more dimensions (see menzel2021bootstrap, olmo2023nonparametric); there is no related asymptotic theory that combines these features in a unified framework (see, also discussion in katsouris2023limit, katsouris2023optimal). 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, ahn1990estimation, toda1995statistical, cavanagh1995inference, paruolo1997asymptotic and poskitt2006identification) along with near unit roots (see, phillips2007limit), 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 properly\footnote{Suitable 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 white2000asymptotic).

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, mikusheva2007uniform and phillips2014confidence). Thus, our contributions in this paper are summarized as below:

itemize• 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 \hyperref[appA]{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, kapetanios2014nonlinear develop asymptotic theory for nonlinear panel models with cross-sectional dependence. Moreover, huang2020two propose a network autoregressive model for two-mode networks where the econometric specification allows for different network autocorrelation coefficients such that

align[align omitted — 137 chars of source]

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 anufriev2015connecting). On the other hand, huang2020two 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 zhu2020multivariate. All aforementioned studies correspond to modelling the relationship between the node-specific responses and a set of exogenous covariates for each node. Recently, meitz2022subgeometrically 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 zhu2017network is replaced by a set of nonstationary regressors. Another relevant framework is proposed by nicholson2017varx 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 Phillips1987time, Phillips1987towards and phillips2007limit. 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 katsouris2021optimal and katsouris2023estimating, katsouris2023statistical, motivated from the aspects of financial connectedness (see, diebold2014network and barunik2018measuring) 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 Adrian2016covar and hardle2016tenet.

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

align[align omitted — 307 chars of source]

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

align[align omitted — 111 chars of source]

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.

align[align omitted — 429 chars of source]

In matrix form we have that the regressors are generated via

align[align omitted — 529 chars of source]

for $i \in \left\{ 1,..., N \right\}$. The NVAR(1) model given by (ref)-(ref) 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 zhu2017network and zhu2020grouped, 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 data\footnote{ 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 agosto2016modeling) as well as armillotta2022nonlinear 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, dou2021generalized).

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

align*[align* omitted — 137 chars of source]

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

align*[align* omitted — 227 chars of source]

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

align*[align* omitted — 290 chars of source]

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

align*[align* omitted — 139 chars of source]

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

assumptionWithin the stationary framework we impose the following assumptions: \begin{itemize} • (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. • (Network Structure) Assume that $\left\{ \boldsymbol{\Omega}_i \right\}$ is a sequence of nonstochastic matrices such that $i \in \left\{ 1,...,N \right\}$. \begin{itemize} • \textbf{(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$. • \textbf{(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$. \end{itemize} • \textbf{(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. \end{itemize}
remarkSome relevant remarks regarding our conditions are summarized below: \begin{itemize} • 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. \end{itemize}

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 fort2005subgeometric, meitz2021subgeometric, meitz2022subgeometrically and matsui2022characterization 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 meitz2008ergodicity, meitz2008stability). Moreover, doukhan2023stationarity 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 kostakis2015Robust and magdalinos2021least, 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 holberg2023uniform). Nevertheless, such conditions can be also thought as a direct implication of the persistence class given by PM2009econometric 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

align*[align* omitted — 229 chars of source]

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 required\footnote{In 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 zhu2017network and zhu2019network as well as the frameworks proposed by magdalinos2009limit.

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

align[align omitted — 86 chars of source]

which satisfies the Black-Scholes differential equation $d J_c(r) \equiv C J_C(r) + d B_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)} d W_u(s)$ and $J_C(r)$ the Ornstein-Uhlenbeck, (OU) process,\footnote{The 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[Uniform Convergence of Covariance Matrices (see, holberg2023uniform)] \ 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 \begin{align*} \boldsymbol{J}_{ \boldsymbol{C} } (t) = \int_0^t e^{ (t-s) \boldsymbol{C} } \boldsymbol{\Gamma} \boldsymbol{\Sigma}^{1/2} d \boldsymbol{W}(s), \ \ with \ \ \boldsymbol{J}_{ \boldsymbol{C} }(0) = \boldsymbol{0}. \end{align*}

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) =

bmatrix[bmatrix omitted — 51 chars of source]

$. 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)$.

theoremSuppose 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 \begin{align*} \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}^{ \floor{nr} } \boldsymbol{X}_{t-1} \boldsymbol{X}_{t-1}^{\prime} \right) \boldsymbol{M}^{-1/2} - r \boldsymbol{I}_p \right\rVert > \delta \right) = 0. \end{align*} 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

align*[align* omitted — 295 chars of source]

Furthermore, it holds that

align*[align* omitted — 261 chars of source]

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

exampleSuppose 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 magdalinos2021least) and the error term $u_t = ( u_{yt}, u_{xt} )^{\prime}$ into the first coordinate and their last $(p-1)$ coordinates \begin{align*} \boldsymbol{Y}_t &= \boldsymbol{\beta} \tilde{\boldsymbol{X}}_{t-1} + u_{yt}, \tilde{\boldsymbol{X}}_t := \boldsymbol{\Gamma} \boldsymbol{X}_t \\ \tilde{\boldsymbol{X}}_t &= \boldsymbol{R}_n \tilde{\boldsymbol{X}}_{t-1} + \boldsymbol{u}_{xt} \end{align*}

Furthermore, we consider the IVX instrumentation procedure which requires that

align[align omitted — 206 chars of source]

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 phillips2007limit and PM2009econometric. 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, phillips2013predictive and M $\&$ PCB (2022)).

Limit theory for Network Vector Autoregression

Time Series Stationarity

Below we follow the framework proposed by zhu2017network 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:

align[align omitted — 188 chars of source]

Estimating the NVAR(1) model using the econometric specification (ref), 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 (ref) in the matrix companion form as below

align[align omitted — 137 chars of source]

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$.

theoremSuppose 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 (ref). The solution has the following form \begin{align} \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}. \end{align}

Assuming the existence of the strictly stationary solution (ref), 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 (ref).

propositionAssume the same conditions as in Theorem (ref). Then, given $\mathbb{Z}$, the strictly stationary solution of (ref) converges to a normal distribution with mean and covariance given as below \begin{align} \boldsymbol{\mu} &= \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 \\ \mathsf{vec} \big\{ \boldsymbol{\Gamma}(0) \big\} &= \sigma^2 \big( \boldsymbol{I} - \boldsymbol{G} \otimes \boldsymbol{G} \big)^{-1} \mathsf{vec} ( \boldsymbol{I} ). \end{align}

The implications of Proposition (ref) 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 structure\footnote{Note 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

align[align omitted — 271 chars of source]

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 (ref) 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

align[align omitted — 198 chars of source]

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.

theoremAssume that the stationary condition $| \beta_1 | + | \beta_2 | < 1$ and technical conditions (C1)-(C3) hold, we then have that \begin{align*} \sqrt{NT} \left( \widehat{\theta} - \theta \right) \to \mathcal{N} \left( 0, \sigma^2 \Sigma^{-1} \right) \end{align*} as min$\left\{ N, T \right\} \to \infty$, where $\Sigma$.
propositionAssume that $T$ is fixed and conditions in Theorem 3 hold. Then, we have that \begin{align*} \sqrt{N} \left( \widehat{\theta} - \theta \right) \to \mathcal{N} \left( 0, \sigma^2 T^{-1} \Sigma^{-1} \right) \end{align*} as $N \to \infty$.
remarkFor example, recently fan2023estimation consider conditions for covariance stationarity of a quantile vector autoregressive model such that \begin{align} \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). \end{align}

Furthermore, in the absence of network dependence our econometric specification reduces to the predictive regression model around the general vicinity of unity as in PM2009econometric. 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

align[align omitted — 230 chars of source]

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 (ref) 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

align[align omitted — 213 chars of source]

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.

theoremUnder Assumption 2 (innovation covariance structure), we then have that \begin{align*} \sqrt{NT} \big( \widehat{\widetilde{\theta}}_{IVX} - \widetilde{\theta}_{IVX} \big) \to \mathcal{MN} \left( 0, \sigma^2 \mathbb{V}^{-1} \right) \end{align*} 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:

align[align omitted — 243 chars of source]

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

align[align omitted — 116 chars of source]

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

align[align omitted — 151 chars of source]

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 (ref) above in the case of the NVAR(1) model.

theoremConsider 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 (ref) of the form: \begin{align} \mathbb{Y}_t = \mathcal{J} \mathbb{Y}_t^{*}, \ \ where \ \ \mathbb{Y}_t^{*} = \left( I_p - G^{*} \right)^{-1} \mathcal{B}_0^{*} + \sum_{j=0}^{\infty} G^{*j} \mathcal{E}_{t-j}^{*} \end{align}

such that the following hold

align*[align* omitted — 273 chars of source]

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 (ref) 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

align[align omitted — 193 chars of source]

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 (ref) gives the limiting distribution of the estimator.

theoremAssume that $\sum_{s=1}^p \left( | \alpha_s | + | \beta_s | \right) < 1$ and technical conditions (C1),(C2), (C4) hold. We then have that \begin{align*} \sqrt{NT} \left( \widehat{\theta}^{*}_{OLS} - \theta^{*}_{OLS} \right) \to \mathcal{N} \left( 0, \sigma^2 \Sigma^{*-1} \right) \end{align*} 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 (ref)-(ref) 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

align*[align* omitted — 199 chars of source]

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 (ref) 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

align[align omitted — 235 chars of source]
theoremUnder Assumption 2 (innovation covariance structure), we then have that \begin{align*} \sqrt{NT} \big( \widehat{ \widetilde{\theta}}^{*}_{IVX} - \widetilde{\theta}^{*}_{IVX} \big) \to \mathcal{MN} \left( 0, \sigma^2 \mathbb{V}^{*-1} \right) \end{align*} 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 (ref).

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.

lemmaLet $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)$.
proofLet $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 \begin{align*} \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 ). \end{align*} 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.

corollaryLet $\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)$.
corollaryFor 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$, \begin{align*} Var \left( \frac{1}{v_n} \sum_{i \in G_n } Z_{n,i} \right) &= \frac{1}{v_n^2} \sum_{i \in G_n } Var \left( Z_{n,i} \right) + \frac{1}{v_n^2} \sum_{ \substack{ i, j \in G_n \\ \kappa(i,j) = 1} } Cov \left( Z_{n,i}, Z_{n,j} \right) + \frac{1}{v_n^2} \sum_{ \substack{ i, j \in G_n \\ \kappa(i,j) = 0} } Cov \left( Z_{n,i}, Z_{n,j} \right) \\ &= r_n^{-1} Var \left( Z_{n, v_1v_2 } \right) + \mathcal{O} \left( r_n^{ -1/2} \right) Cov \left( Z_{n,v_1v_2}, Z_{n,v_2v_3} \right) + \mathcal{O}(1) Cov \left( Z_{n,v_1v_2}, Z_{n,v_3v_4} \right) \end{align*}
remarkFor 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_1v_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_1v_2}, Z_{n,v_3v_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 laurent2022unit as well as the framework of gobet2017parameter. 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 badev2021nash.

wrap\subsection{Quasi-Maximum Likelihood Estimation Approach} Estimation for the unknown parameter vector $\theta$ is developed by means of QMLE. We define the quasi-log-likelihood function for $\boldsymbol{\theta}$ as below \begin{align} \ell_{NT} (\boldsymbol{\theta}) = \sum_{t=1}^T \sum_{i=1}^N \ell_{i,t} (\boldsymbol{\theta}), \end{align} where $\ell_{i,t} (\boldsymbol{\theta})$ is the log-likelihood contribution of a single network node which depends on the data, that is on the persistence properties of regressors in our case. We consider the following QMLE estimator \begin{align} \ell_{NT} (\boldsymbol{\theta}) = \sum_{t=1}^T \sum_{i=1}^N \big( Y_{i,t} \mathsf{log} \lambda_{i,t}(\boldsymbol{\theta}) - \lambda_{i,i}(\boldsymbol{\theta}) \big), \end{align} which is the log-likelihood obtained if all time series were contemporaneously independent. The particular property simplifies computations allowing to establish consistency and asymptotic normality of the resulting estimator. Some of the challenges include the fact that non standard proofs are required for establishing stationarity or infinite-dimensional processes in order to establish robust inference procedures within the double asymptotic regime. Suppose that the functional form of the model is described by $\lambda_t ( \theta)$, then $\hat{\theta}$ corresponds to the maximizer of the least squares criterion such that \begin{align} \ell_{NT} (\theta) = - \sum_{t=1}^{T} \big( Y_t - \lambda_t ( \boldsymbol{\theta} ) \big)^{\prime} \big( Y_t - \lambda_t ( \boldsymbol{\theta} ) \big). \end{align} Then, it follows that \begin{align} \mathcal{S}_{NT}(\theta) = \sum_{t=1}^T \frac{ \partial \lambda_t (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta} } \big( Y_t - \lambda_t ( \boldsymbol{\theta} ) \big) = \sum_{t=1}^T s_{NT} (\boldsymbol{\theta}). \end{align} Furthermore, the empirical Hessian and information matrices are respectively as below \begin{align} H_{NT}(\boldsymbol{\theta}) &= \sum_{t=1}^T \sum_{i=1}^N \frac{ \partial \lambda_{i,t} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta} } \cdot \frac{ \partial \lambda_{i,t} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta}^{\prime} } - \sum_{t=1}^T \sum_{i=1}^N \big( Y_{i,t} - \lambda_{i,t} ( \boldsymbol{\theta} ) \big) \frac{ \partial^2 \lambda_{i,t} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{\prime} } \\ \nonumber \\ B_{NT}(\boldsymbol{\theta}) &= \sum_{t=1}^T \frac{ \partial \lambda_{t}^{\prime} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta} } \boldsymbol{\Sigma}_t (\boldsymbol{\theta}) \frac{ \partial \lambda_{t}^{\prime} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta}^{\prime} }. \end{align} Moreover, we have that \begin{align} H_{N}(\theta) &= \mathbb{E} \left( \frac{ \partial \lambda_{t}^{\prime} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta} } \cdot \frac{ \partial \lambda_{t} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta}^{\prime} } \right) \\ B_N &= \mathbb{E} \left( \frac{ \partial \lambda_{t}^{\prime} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta} } \cdot \xi_t(\theta) \cdot \xi^{\prime}_t(\theta) \frac{ \partial \lambda_{t} (\boldsymbol{\theta}) }{ \partial \boldsymbol{\theta}^{\prime} } \right) \end{align} There, there exists a fixed open neighborhood $\mathcal{D} (\theta_0) = \big\{ \theta : \left| \theta - \theta \right|_2 \big\}$ of $\theta_0$ such that with probability tending to 1 as $\left\{ N, T_N \right\} \to \infty$, then the score function equation $\mathcal{S}_{NT} (\theta) = 0$, has a unique solution which is denoted by $\hat{\theta} \overset{p}{\to} \theta_0$ and $\sqrt{ N T_N } \left( \hat{\theta} - \theta_0 \right) \overset{d}{\to} \mathcal{N} \left( 0, B^{-1} \right)$, where $B = \sigma^2 H$ and $H$ is defined, with $\Sigma_t = D_t = I$.
wrap\subsection{Theoretical Proofs} \paragraph{Proof of Theorem 2 of zhu2017network} \ \color{blue} Write the estimator of $\widehat{\theta}$ such that $\widehat{\theta} = \theta + \widehat{\Sigma}^{-1} \widehat{\Sigma}_{xe}$, where \begin{align} \widehat{\Sigma} = (NT)^{-1} \sum_{t=1}^T \mathbb{X}_{t-1} \mathbb{X}_{t-1} ^{\prime} \ \ and \ \ \widehat{\Sigma}_{xe} = (NT)^{-1} \sum_{t=1}^T \mathbb{X}_{t-1} \mathcal{E}_{t} \end{align} As a result the conclusion of Theorem 3 holds if, \begin{align} \widehat{\Sigma} &\to \Sigma \\ \sqrt{NT} \widehat{\Sigma}_{xe} &\to \mathcal{N} \left( 0 , \sigma^2 \Sigma \right) \end{align} as min$\left\{ N, T \right\} \to \infty$. Below, we prove (ref) in Step 1 and (ref) in Step 2. Step 1. In this step, we prove the following result \begin{align} \widehat{\Sigma} = \frac{1}{NT} \sum_{t=1}^T \mathbb{X}_{t-1}^{\prime} \mathbb{X}_{t-1} = \begin{pmatrix} 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{pmatrix} \to \begin{pmatrix} 1 & c_{\beta} & c_{\beta} & 0 \\ & \Sigma_1 & \Sigma_2 & \kappa_3 \gamma^{\prime} \Sigma_z \\ & & \Sigma_3 & \kappa_8 \gamma^{\prime} \Sigma_z \\ & & & \Sigma_z \\ \end{pmatrix} \end{align} where \begin{align} \mathcal{S}_{12} &= \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_{i=1}^N Z_{i}^{\prime} \nonumber \\ \nonumber \\ \mathcal{S}_{22} &= \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)} \nonumber \\ \nonumber \\ \mathcal{S}_{24} &= \frac{1}{NT} \sum_{t=1}^T \sum_{i=1}^N w_i^{\prime} \mathbb{Y}_{t-1} Z_i^{\prime} , \ \ \mathcal{S}_{33} = \frac{1}{NT} \sum_{t=1}^T \sum_{i=1}^N Y^2_{i(t-1)} \nonumber \\ \nonumber \\ \mathcal{S}_{34} &= \frac{1}{NT} \sum_{t=1}^T \sum_{i=1}^N Y_{i(t-1)} Z_i^{\prime}, \ \ \mathcal{S}_{44} = \frac{1}{N} \sum_{i=1}^N Z_i Z_i^{\prime} \end{align} Therefore, we have that \begin{align} \mathbb{Y}_t = c_{\beta} \mathbf{1} + (I-G)^{-1} \mathbb{Z} \gamma + \widetilde{\mathbb{Y}}_t, \end{align} almost surely. where $\widetilde{\mathbb{Y}}_t = \sum_{j=0}^{\infty} G^j \mathcal{E}_{t-j}$, $\mathbb{Z} \gamma = \left( Z_1^{\prime}, ... , Z_N^{\prime} \right)^{\prime}$ For the case where no information regarding the persistence properties of regressors (that is, we assume that $X_{i(t)}$ are not generated via the LUR specification), then we see that the limiting distribution of the model estimate for the NVAR(1) model weakly converges to a Normal random variable.

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 (ref), 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 (ref) 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

align*[align* omitted — 214 chars of source]

\paragraph{Sketch Proof of Theorem 3}

We define with

align[align omitted — 240 chars of source]

We have that

align[align omitted — 354 chars of source]

where

align[align omitted — 916 chars of source]

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}}$. \color{red} Convergence of $\mathcal{S}_{12}$: \color{black}

align[align omitted — 253 chars of source]

where

align[align omitted — 230 chars of source]

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

\color{red} Convergence of $\mathcal{S}_{13}$: \color{black}

align[align omitted — 232 chars of source]

where

align[align omitted — 226 chars of source]

\color{blue} 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}$. \color{black}

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:

align*[align* omitted — 204 chars of source]

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 armillotta2022nonlinear, terasvirta1994aspects), 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 hardle2016tenet. 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, chen2019tail), that is, financial conntectedness estimated with tail risk measures (see, also daouia2022extremile). Suppose that $t \in \left\{ 1,...,T \right\}$, then by employing a quantile optimization function\footnote{Consider 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

align[align omitted — 173 chars of source]

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 (ref), and $Q_{\tau} \left( u_{j|i,t} \ | \ R_{it}, X_{t-1} \right) = 0$ for the quantile regression (ref).

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 (ref) 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 (ref). Then, the QR estimator for model (ref) is given by

align[align omitted — 232 chars of source]

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

align[align omitted — 234 chars of source]

The QR estimator from the optimization function (ref) 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 (ref). 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 kojevnikov2021limit and lee2019stable 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.

\paragraph{Conflicts of interest}

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

\paragraph{Data availability}

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

\paragraph{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

exampleWe consider the multivariate autoregressive index model representation proposed by reinsel1983some, 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}$ \begin{align} \boldsymbol{Y}_t - \sum_{j=1}^p \boldsymbol{\Phi}_j \boldsymbol{Y}_{t-j} = \boldsymbol{\varepsilon}_t, \end{align} 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 \begin{align} \boldsymbol{A}( L )_{( m \times r)} &= \boldsymbol{A}_1 L + ... + \boldsymbol{A}_{p_1} L^{p_1} \\ \boldsymbol{B}( L )_{( r \times m)} &= \boldsymbol{B}_0 + \boldsymbol{B}_1 L + ... + \boldsymbol{B}_{p_2} L^{p_2} \end{align} Therefore, we have that \begin{align} \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. \end{align} 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 \begin{align} \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, \end{align} Assume that $p_2 = 0$ and $p_1 = p$, then the model can be formulated as below \begin{align} \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 \end{align} 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[see, cubadda2022dimension] Consider the Dimension-Reducible VAR model as \begin{align} \boldsymbol{Y}_t = \sum_{j=1}^p \boldsymbol{A} \alpha_j \boldsymbol{A}^{\prime} \boldsymbol{Y}_{t-j} + \boldsymbol{u}_t, \end{align} 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 \begin{align} \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} \end{align} 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 = A x_t + \varepsilon_t$ to obtain \begin{align} \boldsymbol{Y}_t = \boldsymbol{A} \boldsymbol{\gamma} (L) \boldsymbol{\xi}_t + \boldsymbol{\varepsilon}_t, \end{align} 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}$ \begin{align} \boldsymbol{Y}_t = \underbrace{ \boldsymbol{C} (L) \boldsymbol{\xi}_t }_{ \textcolor{blue}{ \chi_t } } + \boldsymbol{v}_t, \end{align} 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[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. \begin{align} x_t = \alpha_1 A x_{t-1} + ... + \alpha_p A x_{t-p} + v_t \equiv A \sum_{j=1}^p \theta_j x_{t-j} + v_t, \end{align}
proposition[Granger-Causality in NVAR$(p,1)$] \ Suppose that $x_t$ is generated by (ref) 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:

align[align omitted — 225 chars of source]

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.

proofConsider again the functional form of the NVAR$( p, 1 )$ model as below: \begin{align} y_t = \alpha_1 A y_{t-1} + ... + \alpha_p A y_{t-p} + u_t, \end{align} 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, \begin{align} \boldsymbol{X}_t = \big[ A y_{t-1},..., A y_{t-p} \big] \end{align} 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: \begin{align} \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), \end{align} where $\boldsymbol{\Sigma} = \mathsf{Var} ( u_t )$. Therefore, this implies that \begin{align} \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) \end{align} 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: \begin{align} \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) \end{align} As $n \to \infty$, then it holds that \begin{align} \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) \end{align} 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 \begin{align} \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) \end{align} \begin{remark} 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. \end{remark} 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

align[align omitted — 279 chars of source]

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

align[align omitted — 339 chars of source]

Therefore, to obtain the impulse response for $x_t$, one can write the econometric specification in the following companion form:

align[align omitted — 87 chars of source]

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

align[align omitted — 341 chars of source]

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:

align[align omitted — 195 chars of source]
center[center omitted — 265 chars of source]