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.
89,354 characters · 15 sections · 76 citation commands
Estimation of Grouped Time-Varying Network Vector Autoregression Models
\newtheorem{corollary}{Corollary} \newtheorem{definition}{Definition} \newtheorem{lemma}{Lemma} \newtheorem{proposition}{Proposition} \newtheorem{remark}{Remark} \newtheorem{theorem}{Theorem} \newtheorem{example}{Example} \newtheorem{assumption}{Assumption} \newtheorem{prop}{Proposition}
\numberwithin{corollary}{section} \numberwithin{definition}{section} \numberwithin{equation}{section} \numberwithin{lemma}{section} \numberwithin{proposition}{section} \numberwithin{remark}{section} \numberwithin{theorem}{section}
\allowdisplaybreaks[4]
\setcounter{equation}{0}
Modeling large-scale time series has been the main frontier of recent advances of time series analysis and is of fundamental importance in various fields of applications such as climatology, economics, finance and social networks. Since Si80's seminal work, the vector autoregressive (VAR) model has become a commonly-used statistical tool to tackle multivariate time series, see Lu06 and KL17 for a comprehensive review of classic estimation and forecasting techniques. However, the VAR-based estimation and forecasting are challenging when the number of time series sequences $N$ diverges to infinity. In this case, the number of unknown parameters in VAR transition matrices is of order $O(N^2)$, which may be substantially larger than the time series length $T$. In order to construct sensible estimate and forecast, two dimension-reduction approaches are often employed: VAR with sparse transition matrices and regularised estimation SB11, BM15, KC15, DZZ16, MPS23 and factor-augmented VAR BBE05, BN06, BLL16. Although some sound asymptotic properties have been developed for the sparse or factor-augmented VAR estimates, they often neglect possible network structures in large-scale time series and cannot directly capture dynamic network effects.
Consider time series observation vectors $X_{t}=\left(x_{1,t},\cdots,x_{N,t}\right)^{^\intercal}$ with $N$ being the number of nodes in the large-scale network, and denote an adjacency matrix by ${\mathbf W}=(w_{ij})_{N\times N}$, where $w_{ii}=0$, $w_{ij}=1$ for $i\neq j$ if there exists a directed edge from $i$ to $j$ and $w_{ij}=0$ otherwise. The matrix ${\mathbf W}$ is assumed to be observable and can be either directed (${\mathbf W}^{^\intercal}\neq {\mathbf W}$) or undirected (${\mathbf W}^{^\intercal}={\mathbf W}$). The classic network VAR model is defined by
where $\widetilde{w}_{ij}=w_{ij}/n_i$ with $n_i=\sum_{j\neq i} w_{ij}$, $\beta_1$ and $\beta_2$ are unknown parameters, and $\varepsilon_t=(\varepsilon_{1,t},\cdots,\varepsilon_{N,t})^{^\intercal}$ is a sequence of independent and identically distributed (i.i.d.) random vectors. The above network VAR model formulation contains two regression components: $\beta_1 \sum_{j\neq i}\widetilde{w}_{ij}x_{j,t-1}$ and $\beta_2x_{i,t-1}$, corresponding to the network (cross-lag) and momentum (own-lag) effects, respectively. ZPLLW17 discuss stationarity conditions for an extended network VAR model with extra nodal effects, propose the least squares estimation method and derive the relevant asymptotic theory. Although the classic linear network VAR model is easy to interpret and implement, it may be invalid in empirical applications. In particular, there exist two practical issues: (i) the stable network VAR model cannot capture smooth structural changes in the underlying data generating process over a long time span; and (ii) it is often too restrictive to impose the homogeneity assumption on the autoregressive coefficients over $N$ nodes. Consequently, the homogenous linear network VAR ((ref)) may suffer from the model misspecification problem, resulting in biased estimates, inaccurate forecast and misleading inference. There have been some attempts in recent years to address one of the aforementioned two issues. To incorporate structural changes in the autoregressive structures, Su16, SM18, Wu19, CLLL23 and YSM24 extend the linear network VAR models, allowing the coefficients to vary smoothly over time or with a stationary index variable. On the other hand, to relax the homogeneity restriction in linear VAR models, YSM23, YSM24 introduce fully heterogenous network VAR models with the momentum and network spillover effects varying over nodes. However, the number of unknown coefficients in the latter models grow with the number of nodes, resulting in slow estimation convergence when the time span is not sufficiently long.
This paper aims to jointly tackle the aforementioned two issues by introducing a general time-varying network VAR model framework satisfying a latent group structure, i.e., the time-varying network autoregressive relationships are invariant within a group of nodes, but change over different groups. The grouped time-varying network VAR model achieves a good balance between model flexibility and parsimony. It not only covers the homogenous network VAR models ZPLLW17, Wu19 as a special case, but also provides a more parsimonious model formulation than the fully heterogenous network VAR models YSM23, YSM24, achieving dimension reduction in estimation and improving the subsequent out-of-sample forecasting performance. The main methodological and theoretical contributions of our paper with connection to the existing literature are summarized as follows.
The finite-sample Monte-Carlo simulation study shows that the proposed clustering algorithm and ratio criterion can consistently estimate the latent group membership and number as long as the time series length $T$ is moderately large; the post-grouping local linear estimates perform significantly better than the naive heterogeneous local linear estimates; and the developed two-stage estimation procedure can precisely locate the break point and estimate the group structure before and after the break. The developed model and methodology are further applied to analyze the monthly temperature data collected in $37$ UK weather stations over the time period between January 1950 and February 2023. The empirical result reveals that there exist two groups over the node-specific time-varying coefficients: weather stations in Northern Ireland and Wales tend to form one group and those in England and Scotland form the other. In addition, our proposed method outperforms other competing methods in terms of out-of-sample prediction.
The rest of the paper is organized as follows. Section (ref) introduces the time-varying network VAR model together with the latent group structure and provides some fundamental model assumptions. Section (ref) describes the estimation methodology for the group structure and establishes the consistency property for the resulting estimation. Section (ref) proposes the post-grouping nonparametric estimation and derives the relevant asymptotic theory. Section (ref) extends the model, method and theory to the case with breaks in either the group structure or group-specific coefficient functions. Section (ref) reports both the simulation and empirical studies. Section (ref) concludes the paper. Appendix A introduces the clustering algorithm to identify the homogeneity structure on the network effects. The online supplement LPTW24 contains proofs of the main asymptotic theorems, some technical lemmas with their proofs, discussion on the refined break point estimation, and extra numerical results. Throughout the paper, we let $\lfloor\cdot\rfloor$ and $\lceil\cdot\rceil$ be the floor and ceiling functions, respectively. For a vector $x=(x_1,\cdots,x_p)^{^\intercal}\in\mathscr{R}^p$, we write $\vert x\vert_q=(\sum_{i=1}^p|x_i|^q)^{1/q}$ for $q\geq1$ and $\vert x\vert_\infty=\max_{1\leq i\leq p}|x_i|$; for a matrix $\boldsymbol{\Sigma}=(\sigma_{ij})_{p\times p}$, write $\vert\boldsymbol{\Sigma}\vert_{\sf F}=(\sum_{i=1}^{p}\sum_{j=1}^{p}\sigma_{ij}^2)^{1/2}$, $\vert\boldsymbol{\Sigma}\vert_\infty=\max_{1\leq i,j\leq p}\vert\sigma_{ij}\vert$ and $\vert\boldsymbol{\Sigma}\vert_{\sf O}=\max_{x\in\mathscr{R}^p: \vert x\vert_2=1}\vert\boldsymbol{\Sigma} x\vert$; and for a $p$-dimensional random vector $Z$, write $Z\in\mathscr{L}^\kappa$ if $\Vert Z\Vert_\kappa:=[{\sf E}(\vert Z\vert_2^\kappa)]^{1/\kappa}<\infty, \kappa\ge1$. Let $\mathbf{I}_k$ be a $k\times k$ identity matrix and $\mathbf{O}_{k\times l}$ a $k\times l$ zero matrix. For a square matrix, $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ denote the minimum and maximum eigenvalues and ${\sf det}(\cdot)$ denotes its determinant. Let $a_n=o(b_n)$, $a_n=o_P(b_n)$ and $a_n\propto b_n$ denote that $a_n/b_n\to 0$ as $n\to\infty$, $a_n/b_n\to 0$ with probability approaching one ({\em w.p.a.1}), and $0<\underline{c}\leq a_n/b_n\leq \overline{c}<\infty$, respectively.
\setcounter{equation}{0}
In this section, we introduce the main model framework, i.e., the time-varying network VAR model with a latent group structure, and impose some fundamental assumptions, ensuring the network time series are locally stable.
Suppose that there exists a partition of the node index set $\{1,2,\cdots,N\}$, denoted by ${\mathscr{G}}=\{{\mathscr G}_1,{\mathscr G}_2,\cdots,{\mathscr G}_{K_0}\}$, such that ${\mathscr G}_i\cap {\mathscr G}_j=\emptyset$ for $1\leq i\neq j\leq K_0$. Let $g_i\in\{1,\cdots,K_0\}$ be the group membership of the $i$-th node, i.e., $g_i=k$ is equivalent to $i\in{\mathscr G}_k$. Neither the group membership nor the group number is known a priori. Consider the following grouped time-varying network VAR model:
where $\tau_t=t/T$ denotes the scaled time point, $\alpha_{g_ig_j}(\cdot)$ and $\alpha_{g_i}(\cdot)$ are smooth coefficient functions, and the remaining elements are defined as those in ((ref)). In contrast to the linear network VAR model ((ref)), the grouped time-varying network VAR model ((ref)) provides a much more flexible framework, allowing the network and momentum effects to change over time and nodes. In particular, the interaction between nodes from two groups share the same time-varying network effects which is appealing for modeling social networks with smooth structural changes. Although we only consider one lag in model ((ref)) for simplicity of exposition, the method and theory developed in Sections (ref) and (ref) below can be easily extended to the model setting with finite lags. It is worth pointing out that our network structure is deterministic ZPLLW17 and the group membership is determined by node-specific time-varying momentum and network spillover effects. This is substantially different from the community network VAR model studied by CFZ23, where the network structure is random and the community structure is used for the network generating mechanism.
Let $\widetilde{\mathbf W}$ be the row-normalized adjacency matrix with the $(i,j)$-entry being $\widetilde{w}_{ij}$, ${\mathbf B}_1(\cdot)$ be an $N\times N$ matrix being the diagonal entries being zeros and the off-diagonal $(i,j)$-entry being $\alpha_{g_ig_j}(\cdot)$ and ${\mathbf B}_2(\cdot)={\sf diag}\{\alpha_{g_1}(\cdot),\cdots,\alpha_{g_N}(\cdot)\}$. Then, we may rewrite model ((ref)) as
where $\odot$ denotes the Hadamard product between matrices. Model ((ref)) thus falls within the high-dimensional time-varying VAR model framework which has received increasing attention in recent years. For instance, DQC17 propose a kernel-weighted $\ell_1$-regularised estimation for a time-varying VAR model; XCW20 study a high-dimensional VAR model with multiple breaks and estimate smooth time-varying covariance and precision matrices between the break points; CLLL23 estimate dual network structures via directed Granger causality and undirected partial correlation linkages within the high-dimensional time-varying VAR framework. However, the aforementioned literature often assumes a sparsity condition on the time-varying VAR transition matrices to facilitate the use of the regularised estimation techniques and cannot directly capture possible time-varying network effects. In this paper, we decompose the time-varying transition matrix into two components: ${\mathbf B}_2(\tau_t)$ capturing the momentum effects and ${\mathbf B}_1(\tau_t)\odot\widetilde{\mathbf W}$ capturing the dynamic network spillover effects.
The proposed model ((ref)) contains the homogenous time-varying network VAR model Wu19 as a special case, i.e., \[ x_{i,t}=\alpha_{1}^\dagger(\tau_t)\sum_{j\neq i}\widetilde{w}_{ij}x_{j,t-1}+\alpha_{2}^\dagger(\tau_t)x_{i,t-1}+\varepsilon_{i,t}. \] In contrast to the fully heterogenous network VAR model YSM24, our model achieves substantial dimension reduction, reducing the number of unknown coefficient functions to $K_0^2+K_0$ which is finite (as $K_0$ is assumed to be fixed). In the homogenous or grouped linear network VAR ZPLLW17, CFZ23, ZXF23, the so-called nodal effect is often incorporated in the model formulation. Hence, we may further extend model ((ref)) to
where $Z_i$ is a $p$-dimensional vector of node-specific exogenous covariates and $\gamma_{g_i}(\cdot)$ is a vector of smooth coefficient functions. Letting $\{Z_i\}$ be independent of $\{\varepsilon_{i,t}\}$, the model framework and methodology developed in the sequel can be extended to tackle ((ref)) with slight modification. However, for notational simplicity, we mainly focus on model ((ref)) without the nodal effect throughout the paper.
Let $f^{\prime}(\cdot)$ and $f^{\prime\prime}(\cdot)$ be the first- and second-order derivatives of $f(\cdot)$. We impose the following assumption on model ((ref)).
\setcounter{equation}{0} \setcounter{prop}{0}
In this section, we introduce the methodology for estimating the latent group membership and number, and present the relevant asymptotic properties.
We next introduce the nonparametric estimation and clustering methods and obtain the consistent estimation of the group membership ${\mathscr G}$ when $K_0$ is known a priori. The methodology can be split into the following three stages.
{\bf Stage 1}:\ \ As there is no prior information on the latent group structure, we start with the fully heterogenous time-varying network VAR model:
where ${\mathscr N}_i=\{j\neq i:\ w_{ij}=1\}$ is the index set of nodes which the $i$-th node follows, $\beta_{ij}(\cdot)=\alpha_{g_ig_j}(\cdot)$ and $\beta_i(\cdot)=\alpha_{g_i}(\cdot)$. Model ((ref)) is similar to that in YSM24. Note that $\beta_{ij}(\cdot)$, $j\notin{\mathscr N}_i$, are unidentifiable and thus not estimable. With the smoothness condition in Assumption (ref)(i), we adopt the local linear smoothing FG96 to estimate the heterogenous time-varying coefficient functions, only using the sample information from the $i$-th node and its direct neighbors, i.e., $\widetilde{w}_{ij}\neq0$.
Define \[ \widetilde{X}_{i,t-1}=\left[\left(\widetilde{w}_{ij}x_{j,t-1}:\ j\in{\mathscr N}_i\right)^{^\intercal}, x_{i,t-1}\right]^{^\intercal}, \] which is a random vector with dimension $n_i+1$, where $n_i={\rm card}({\mathscr N}_i)$ is allowed to diverge slowly to infinity and ${\rm card}(\cdot)$ denotes the cardinality of a set. Letting \[ \beta_{i\bullet}(\tau)=\left[\left(\beta_{ij}(\tau):\ j\in{\mathscr N}_i\right)^{^\intercal}, \beta_i(\tau)\right]^{^\intercal}, \] with Assumption (ref)(i), we have the following Taylor expansion: \[ \beta_{i\bullet}(\tau_t)\approx \beta_{i\bullet}(\tau)+\beta_{i\bullet}^{\prime}(\tau)(\tau_t-\tau) \] when $\tau_t$ falls in a small neighborhood of $\tau$. Define the node-specific local linear weighted objective function:
where $a$ and $b$ are $(n_i+1)$-dimensional vectors, $K_h(\cdot)=\frac{1}{h}K(\cdot/h)$, $K(\cdot)$ is a kernel function and $h$ is a bandwidth. Minimizing ${\cal L}_i(a,b)$ with respect to the vectors $a$ and $b$, we obtain the solution denoted as $\widehat a$ and $\widehat b$, and then the local linear estimate of $\beta_{i\bullet}(\tau)$ as
In practice, we obtain the local linear estimates at $\tau_l^\ast$, $l=1,\cdots,L$, a sequence of user-specified equidistant grid points between $0$ and $1$ satisfying $L\rightarrow\infty$ and $L=O(T)$.
{\bf Stage 2}:\ \ Let $\bar{\mathscr N}_N=\{(i,j): 1\leq i\leq N, j\in{\mathscr N}_i\}$. It follows from ((ref)) that there exists a latent homogeneity structure for $\beta_{ij}(\cdot)$, $(i,j)\in\bar{\mathscr N}_N$, and the number of distinct time-varying coefficient functions is at most $K_0^2$. Let $\beta_m^\circ(\cdot)$, $m=1,\cdots,M_0$, denote the true distinct time-varying coefficient functions for network spillover effects, $M_0\leq K_0^2$, and $g_{ij}\in\{1,\cdots,M_0\}$ be the group membership for the index pair $(i,j)\in\bar{\mathscr N}_N$. With $\{\widehat{\beta}_{ij}(\tau_l^\ast):\ 1\leq l\leq L, (i,j)\in\bar{\mathscr N}_N\}$ obtained in Stage 1, combining the clustering algorithm and the ratio criterion with details provided in Appendix A (see also Stage 3 and Section (ref)), we may obtain a consistent estimate of $M_0$ denoted by $\widehat{M}$, and the estimated membership $\widehat{g}_{ij}$. For the $i$-th node, we construct
where $\tau$ is chosen as the grid points $\tau_{l}^\ast$ defined in Stage 1, \[ \widehat{\omega}_{ij,m}=\left\{
\right. \] and $I(\cdot)$ denotes the indicator function.
{\bf Stage 3}:\ \ With the estimates $\widehat{\beta}_i(\cdot)$ and $\widehat{\beta}_{i\bullet}^\circ(\cdot)$ defined in Stages 1 and 2, respectively, we may compute the point-wise distance between nodes $i$ and $j$:
and subsequently define the distance matrix: \[ \widehat{\mathbf D}=\left\{\widehat{D}_{ij}\right\}_{N\times N},\quad \widehat{D}_{ij}=\frac{1}{L}\sum_{l=1}^L\widehat{d}_{ij}(\tau_l^\ast). \] It is clear that the diagonal elements of $\widehat{\mathbf D}$ are zeros. With the distance matrix $\widehat{\mathbf D}$, we may adopt the agglomerative hierarchical clustering algorithm which is commonly used in unsupervised cluster analysis HTF09, ELLS11. This clustering algorithm has been recently combined with the kernel-based estimation technique to identify the homogeneity/group structure in nonparametric panel regression models. For instance, Ch19 constructs a similar distance matrix and further estimates the latent group structure in time-varying coefficient panel data models; and VL20 introduce a bandwidth-free normalized distance measure in the clustering algorithm but assume the panel observations are independent over subjects. The latter assumption may be too restrictive for large-scale network time series data and is thus removed in this paper. Another relevant paper is Z13 which clusters nonlinear trend functions based on parallelism and allows the number of time series to grow at a slow polynomial rate of $T$. In contrast, the number of nodes $N$ can be much larger than $T$ in this paper, see ((ref)) in Assumption (ref)(iii).
Assuming the true group number $K_0$ is known a priori, we start with $N$ clusters each of which corresponds to one node, search for the smallest off-diagonal entry in $\widehat{\mathbf D}$ (which is the smallest estimated distance between nodes), and merge the two corresponding nodes. Consequently the cluster number reduces from $N$ to $N-1$. Use a linkage technique (such as the single or complete linkage) to calculate the distance between the merged cluster and the remaining ones and update the estimated distance matrix with size $(N-1)\times (N-1)$. Repeat the previous steps with the updated distance matrix, and stop the algorithm when the number of clusters reaches $K_0$. We denote the estimated clusters by $\widehat{\mathscr G}_k$, $k=1, \cdots,K_0$.
Let ${\boldsymbol\Delta}_{i,t}={\sf E}[\widetilde{X}_{i,t}\widetilde{X}_{i,t}^{^\intercal}]$ with $\widetilde{X}_{i,t}$ defined in Stage 1, and write
and \[ \omega_{ij,m}=\left\{
\right. \] The latter is estimated by $\widehat{\beta}_{i\bullet}^\circ(\tau)$ defined in ((ref)) (up to permutation). The following conditions are required to derive the consistency property of $\widehat{\mathscr G}_k$, $k=1, \cdots,K_0$.
Theorem (ref) below establishes the consistency property of the group membership estimation when $K_0$ is pre-specified.
\setcounter{theorem}{0}
In practice, the true number of groups is unknown and a data-driven criterion is thus required to obtain its consistent estimation. We next introduce an easy-to-implement ratio criterion to estimate $K_0$ consistently. Assuming the group number to be $K$, we may terminate the clustering algorithm in Stage 3 when the cluster number reaches $K$ and obtain the estimated clusters denoted by $\widehat{\mathscr G}_{k|K}$, $k=1,\cdots,K$. With these estimated clusters, we pool the heterogenous time-varying coefficient estimates $\widehat{\beta}_i(\cdot)$ and $\widehat{\beta}_{i\bullet}^\circ(\cdot)$ over $i\in\widehat{\mathscr G}_{k|K}$, and obtain
Then, we define the average deviation:
The grouped time-varying network VAR model is either correctly- or over-fitted when $K\geq K_0$, indicating that $\widehat{R}(K)$ converges to zero, and is under-fitted when $K<K_0$. For the latter scenario, at least two groups would be falsely merged, leading to biased estimation of some group-specific time-varying coefficient functions and a relatively large value of $\widehat{R}(K)$. Hence, it is sensible to estimate $K_0$ by
where $\overline{K}$ is a pre-specified positive integer larger than $K_0$. In practical implementation, we set $\widehat{R}(1)/\widehat{R}(0)\equiv1$, $\widehat{R}(K)=0$ if it is smaller than $\rho_{NT}$, a user-specified tuning parameter, and define $0/0\equiv1$. A similar ratio criterion is adopted by YCLL23 to consistently estimate the group number in nonparametric grouped panel quantile regression models. Other applications of the ratio criterion can be found in LY12 and LRS20.
We require some further conditions to derive the consistency property of $\widehat K$.
We establish the consistency property of the ratio criterion in the following theorem.
Finally, we use the estimated group number $\widehat K$ and terminate the clustering algorithm in Section (ref) when the cluster number reaches $\widehat K$. In order to avoid unnecessary notational burden, we still denote the estimated groups as $\widehat{\mathscr G}_k$, $k=1, \cdots,\widehat{K}$. Combining Theorems (ref) and (ref), we readily have the following result.
\setcounter{corollary}{2}
\setcounter{equation}{0}
The heterogenous local linear estimation defined in ((ref)) only makes use of the sample information from node $i$ and its direct neighbors, resulting in rather slow convergence rates (see Lemma D.3 in the supplement) and unstable numerical performance if $T$ is not sufficiently large in finite samples. We next aim to address this issue by pooling the sample information over nodes in the same cluster and proposing a post-grouping local linear estimation.
It follows from Corollary (ref) that, for any $k=1,\cdots,K_0$, there exists $1\leq k^\dag\leq \widehat{K}$ such that ${\mathscr G}_k=\widehat{\mathscr G}_{k^\dag}$ {\em w.p.a.1}. Without loss of generality, we may consider ${\mathscr G}_k=\widehat{\mathscr G}_{k}$ (conditioning on $\widehat{K}=K_0$) throughout this section. For $i\in\widehat{\mathscr G}_k$, define \[ \check{X}_{i,t-1}=\left[\sum_{j\in\widehat{\mathscr G}_1}\widetilde w_{ij} x_{j,t-1},\cdots,\sum_{j\in\widehat{\mathscr G}_{\widehat{K}}}\widetilde w_{ij} x_{j,t-1}, x_{i,t-1}\right]^{^\intercal}, \] which is a random vector with dimension $\widehat{K}+1$. For the $k$-th group, let \[ \alpha_{k\bullet}(\tau)=\left[\alpha_{k1}(\tau),\cdots,\alpha_{kK_0}(\tau), \alpha_k(\tau)\right]^{^\intercal} \] be a vector of true group-specific time-varying coefficient functions to be estimated. For each $k$ and given $\tau\in(0,1)$, we define the following post-grouping local linear objective function:
where $h_\dag$ is a bandwidth which may be different from $h$ used in the heterogenous local linear estimation ((ref)) and ((ref)). Minimizing the post-grouping objective function with respect to the vectors $a$ and $b$, we obtain the solutions denoted by $\check a$ and $\check b$, and construct the post-grouping local linear estimation as
Let $\sigma_{ij}={\sf E}\left(\varepsilon_{i,t}\varepsilon_{j,t}\right)$ and \[ {\boldsymbol\Delta}_{ij}^\diamond(\tau)={\sf E}\left[X_{i,t}^\diamond(\tau) X_{j,t}^{\diamond^\intercal}(\tau)\right],\ \ 1\leq i,j\leq N, \] where \[ X_{i,t}^\diamond=\left[\sum_{j\in {\mathscr G}_1}\widetilde w_{ij} x_{j,t}^\circ(\tau),\cdots,\sum_{j\in {\mathscr G}_{K_0}}\widetilde w_{ij} x_{j,t}^\circ(\tau), x_{i,t}^\circ(\tau)\right]^{^\intercal} \] with $x_{i,t}^\circ(\tau)$ being the $i$-th element of $X_{t}^\circ(\tau)$ defined in ((ref)).
\setcounter{theorem}{0}
\setcounter{equation}{0}
The model, methodology and theory developed in Sections (ref)--(ref) rely on the assumption that the latent group structure is time invariant and the group-specific coefficient functions are smooth over the entire time span. As discussed in the introductory section, this assumption may be restrictive for some empirical applications. Hence, we next make a further extension of the model, methodology and theory, allowing structural breaks in either the group membership, group number or group-specific coefficient functions. Our main interest lies in locating the break point and estimating the group structure before and after the break. We mainly consider the case of a single break for notational brevity and will briefly discuss its extension to the case of multiple breaks later in Remark (ref)(ii).
Assume that the break occurs at an unknown time point $t_0$. Let ${\mathscr{G}}^1=\{{\mathscr G}_1^1,\cdots,{\mathscr G}_{K_1}^1\}$ and $g_i^1\in\{1,\cdots,K_1\}$ be the group structure and membership label (for node $i$) before the break, whereas let ${\mathscr{G}}^2=\{{\mathscr G}_1^2,\cdots,{\mathscr G}_{K_2}^2\}$ and $g_i^2\in\{1,\cdots,K_2\}$ be defined similarly for those after the break. Consider the time-varying network VAR model with break in the group structure:
where $\alpha_{g_i^1g_j^1}^1(\cdot)$ (or $\alpha_{g_i^2g_j^2}^2(\cdot)$) and $\alpha_{g_i^1}^1(\cdot)$ (or $\alpha_{g_i^2}^2(\cdot)$) are the smooth time-varying network spillover and momentum effects before (or after) the break. Model ((ref)) can be seen as an extension of the linear panel model framework (with break in the group structure) in LOW23 and WPS23, taking into account the smooth time-varying feature and network structure.
Throughout this section, assume that $\bar n$ is bounded. We next introduce a two-stage estimation procedure with break location estimation in Stage 1 and then group estimation in Stage 2.
{\bf Stage 1}:\ \ As the time-varying group structure is latent, similar to ((ref)), we first consider the fully heterogenous time-varying network VAR model with break at $t_0$:
where
Write \[ \beta_{i\bullet}^\ddag(\tau_t)=\left[\left(\beta_{ij}^\ddag(\tau_t):\ j\in{\mathscr N}_i\right)^{^\intercal}, \beta_i^\ddag(\tau_t)\right]^{^\intercal} \] as in Section (ref), and let \[ \beta_{i\bullet}^{\ddag,{\sf l}}(\tau)=\lim_{x\uparrow\tau} \beta_{i\bullet}^\ddag(x)\quad {\rm and}\quad \beta_{i\bullet}^{\ddag,{\sf r}}(\tau)=\lim_{x\downarrow\tau} \beta_{i\bullet}^\ddag(x) \] denote the left and right limits of $\beta_{i\bullet}^\ddag(\tau)$, respectively. Define
It follows from ((ref))--((ref)) that $\delta_\beta(t)$, $t=1,\cdots,T$, achieve the maximum at $t=t_0$, which motivates the subsequent estimation procedure. Specifically, we estimate $\delta_\beta(t)$ by a one-sided kernel smoothing method\footnote{The extension to the one-sided local polynomial smoothing is straightforward, see CWW22.} and then locate the break point by maximizing the estimate of $\delta_\beta(t)$ over $t$.
Let $K^\ddag(\cdot)$ be a one-sided kernel function with a compact support $[0,1]$, say, the one-sided version of the Epanechnikov kernel. Define
where $h_\ddag$ is the bandwidth and $\widetilde{X}_{i,s-1}$ is defined as in Section (ref). We estimate $\beta_{i\bullet}^{\ddag,{\sf l}}(\tau_t)$ and $\beta_{i\bullet}^{\ddag,{\sf r}}(\tau_t)$ by
respectively, and then construct
The estimation of $t_0$ is defined as
{\bf Stage 2}:\ \ Let \[ {\mathscr T}_1=\{1,2,\cdots,\widehat{t}-\lfloor\epsilon_0T\rfloor\}\quad {\rm and}\quad {\mathscr T}_1=\{\widehat{t}+\lfloor\epsilon_0T\rfloor,\cdots,T-1,T\}, \] where $\epsilon_0$ is an arbitrarily small positive number (say $0.01$). From Theorem (ref)(i), the group membership and number are time invariant {\em w.p.a.1} over the two time periods ${\mathscr T}_1$ and ${\mathscr T}_2$. Hence, we may adopt the clustering algorithm and ratio criterion developed in Section (ref) to estimate ${\mathscr G}_k^1$ and $K_1$ (using the network time series sample over ${\mathscr T}_1$) as well as ${\mathscr G}_k^2$ and $K_2$ (using the sample over ${\mathscr T}_2$). We denote the resulting estimates as $\widehat{\mathscr G}_k^1$, $\widehat K_1$, $\widehat{\mathscr G}_k^2$ and $\widehat K_2$, whose consistency property is derived in Theorem (ref)(ii) below.
The following conditions are required to derive the asymptotic property of the above two-stage estimation method.
Let $\beta_{i\bullet}^{\circ1}(\cdot)$ be defined similarly to $\beta_{i\bullet}^\circ(\cdot)$ in ((ref)) but with $\beta_{ij}(\cdot)$ replaced by $\beta_{ij}^1(\cdot)$. As $\beta_m^\circ(\cdot)$ defined in Stage 2 of Section (ref), we let $\beta_m^{\circ1}(\cdot)$, $m=1,\cdots,M_1$, denote the distinct coefficient functions for time-varying network effects before the break, where $M_1\leq K_1^2$. The definitions of $\beta_{i\bullet}^{\circ2}(\cdot)$ and $\beta_m^{\circ2}(\cdot)$, $m=1,\cdots,M_2$, are analogous. Define
Let ${\mathscr N}_i^1(j)=\{k\in{\mathscr N}_i: g_k^1=j\}$ and ${\mathscr N}_i^2(j)=\{k\in{\mathscr N}_i: g_k^2=j\}$.
\setcounter{theorem}{0}
\setcounter{equation}{0}
In this section, we conduct both the simulation and empirical studies. Sections (ref) and (ref) assess the finite-sample performance of the developed methodology and verify the main convergence properties via simulation, whereas Section (ref) reports the empirical application to a network time series data set for UK temperature.
We use the grouped network time-varying VAR model (ref) for data generation. The entries of the adjacency matrix ${\mathbf W}$ are defined by $w_{ij} = I(u_{ij}\le \bar{w})$, where $u_{ij}\sim {\sf U}(0,1)$ and $0<\bar{w}<1$, controlling the sparsity level of the network structure. The innovation vectors are generated by $\varepsilon_t\sim {\sf N}({\bf 0}, {\boldsymbol\Sigma}_\varepsilon)$ independently over $t$, where ${\boldsymbol\Sigma}_\varepsilon=\{\sigma_{ij}\}_{N\times N}$ with $\sigma_{ij}=0.1^{|i-j|}$, allowing cross-sectional dependence over components. Consider $K_0=2$ and define the group-specific coefficients as
The group membership is generated as follows: assign each node $i\in \{1,\ldots, N\}$ to $\mathscr{G}_1$ and $\mathscr{G}_2$ with respective probabilities $0.65$ and $0.35$. Note that there exists a further group structure on the time-varying spillover effects $\alpha_{g_ig_j}(\cdot)$ with $M_0=2$, as described in Stage 2 of Section (ref). In the simulation, we consider two scenarios for generating the group membership: (i) fixed group, i.e., the group membership is only generated once and remains the same over replications; and (ii) random group, i.e., the group membership is randomly generated for each replication. We conduct the simulation over $R=1000$ replications and set $N=100,200$, $T=300,600$, and $\bar{w} = 0.025,0.075$.
We use the Epanechnikov kernel in the local linear smoothing, where the bandwidth is determined by the rule of thumb (SW17): $h=(2.35/\sqrt{12})T^{-1/5}$ for the fully heterogenous local linear estimation ((ref)) and $h_\dag=(2.35/\sqrt{12})[\text{card}(\widehat{\mathscr{G}}_j)T]^{-1/5}$ for the post-grouping estimation ((ref)), where $\widehat{\mathscr{G}}_j$ denotes the estimated group. We notice that the clustering results are insensitive to the bandwidth choice in our simulation. For each simulated data set, we first estimate the group membership and number as in Section (ref), and then conduct the post-grouping estimation as in Section (ref). To evaluate the group structure estimation accuracy, we adopt the following two measurements:
where $\widehat{K}_r$ and $\widehat{\mathscr{G}}_{k, r}$ denote the estimates of the group number and membership in the $r$-th replication. The purity quantity ${\sf Purity}({\mathscr G})$ is a simple and transparent evaluation measure with value close to one when the clustering method is precise. To compare the estimation performance between the pre-grouping local linear estimation and the post-grouping one, we compute the root mean squared errors for the estimated momentum and spillover effects:
where $ \beta_{i\bullet}^\ast(\cdot)=\left(\beta_{ij}(\cdot):\ j\in{\mathscr N}_i\right)^{^\intercal}$, $\widetilde{\beta}_{r,i}(\cdot)$ and $\widetilde{\beta}_{r,i\bullet}^{*}(\cdot)$ stand for the estimates in the $r$-th replication, $\tau_s = 0.05, 0.1,\cdots, 0.95$ and $S=19$.
The simulation results are summarized in Tables (ref) and (ref). The numbers in parentheses of Table (ref) are standard deviations of $\sf{RMSE}_{M,r}$ and $\sf{RMSE}_{S,r}$ over 1000 replications. It follows from Table (ref) that both ${\sf AC}(K_0)$ and ${\sf Purity}({\mathscr G})$ converge to one as the time series length $T$ increases from $300$ to $600$, and the results remain stable when the sparsity level $\bar w$ changes from $0.025$ to $0.075$. Table (ref) shows that the post-grouping time-varying coefficient estimation substantially outperforms the pre-grouping one, confirming that the estimation accuracy is significantly improved by making use of the estimated group structure. The standard deviations are generally small, indicating that the nonparametric estimation performance is stable over replications. In addition, both the pre-grouping and post-grouping local linear estimates deteriorate when $\bar w$ increases from $0.025$ to $0.075$.
We next examine the numerical performance of the estimation method with a break in the group structure introduced in Section (ref). The break point is set at $t_0=\lfloor T/2\rfloor +1$ when we re-assign each node $i$ to $\mathscr{G}_1$ and $\mathscr{G}_2$ with respective probabilities $0.65$ and $0.35$. This results in a break in the group membership. Before the break time, the data generating process is the same as that in Section (ref), whereas, after the break, the group-specific time-varying coefficients are defined as \[ \alpha_{g_i} (\tau) = \left\{
\right. \quad\quad \alpha_{g_ig_j} (\tau) = \left\{
\right. \] As in Section (ref), we consider both the fixed and random groups when generating the group membership (with a break) over $R=1000$ replications. In order to obtain stable finite-sample performance, we slightly increase $T$ from $(300,600)$ to $(400,800)$. The number of nodes remains as $N=100$ and $200$.
The one-sided version of the Epanechnikov kernel function is adopted in our nonparametric method. We first estimate the break point via ((ref)), compute the (scaled) measurement $\Delta_{t_0} = (\widehat{t}-t_0)/T$, and then report the means and standard deviations (in parentheses) of $\Delta_{t_0}$ in Table (ref). It is clear that the scaled break point $t_0/T$ can be accurately detected, and its estimation accuracy is not sensitive to the sparsity level of the adjacency matrix ${\mathbf{W}}$. With the estimated break point, we may split the entire time period into the "pre-break" and "post-break" periods and compute their respective ${\sf AC}(K_0)$ and ${\sf Purity}({\mathscr G})$. We further take a simple average of those values over the two periods and report them in Table (ref). We note that both the group number and membership are estimated very accurately even when there exists a break in the group structure. We finally compare the estimation performance between the pre-grouping and post-grouping local linear estimation after the break point is detected. When computing ${\sf RMSE}_{\sf M}$ and ${\sf RMSE}_{\sf S}$, we choose $\tau_s = 0.05, 0.1,\ldots, 0.4$ for the “pre-break" period and $\tau_s = 0.6, 0.65,\ldots, 0.95$ for the “post-break" period, avoiding possible boundary effect in the estimation. The general pattern in Table (ref) is very similar to that in Table (ref), again confirming the significant advantage of the post-grouping estimation.
Appendix F in the online supplement LPTW24 contains extra simulation results: the finite-sample performance of the clustering algorithm introduced in Appendix (ref) on estimating the homogeneity structure for the network spillover effects, and the clustering result by using ZXF23's grouped network VAR model with constant coefficients.
There has been increasing interest in investigating the spatial pattern of climate data, see, for example, PSH2009, Kendon20, Hanlon21 and the references therein. We next apply the proposed model and methodology to analyze a set of UK climate data, exploring the network and latent group structures and allowing for smooth structural changes to account for possible climate changes in the past few decades. The data that we use are collected from the UK Meteorological Office,\footnote{https://www.metoffice.gov.uk/research/climate/maps-and-data/historic-station-data} containing temperature recordings (in Celsius) of 37 weather stations with their geographical locations presented in Figure (ref). The original dataset collects minimum and maximum temperatures per month over the period from January 1950 to February 2023. Hence, the time series length is $T=878$. We consider the following two scenarios when building the network VAR model: (i) both the minimum and maximum temperatures are used as elements of $X_t$; (ii) the averaged temperature per month is used as elements of $X_t$. In model (ii), each weather station is treated as a node and $N=37$; whereas in model (i), the minimum and maximum recordings in each weather station are treated as two nodes and $N=74$. The time series observations are standardized to have zero mean and unit standard deviation. The adjacency matrix ${\mathbf W}$ is constructed by following the UK climate region map\footnote{https://www.metoffice.gov.uk/research/climate/maps-and-data/about/regions-map}. Specifically, it is classified into the following five regions: Southern England, Northern England, Wales, Scotland, and Northern Ireland. When stations $i$ and $j$ are in the same region, we set $w_{ij}=1$, otherwise, $w_{ij}=0$. Consequently, the percentage of non-zero elements of the adjacency matrix is 0.0559 and 0.2267 for the two models.
The primary interest of this empirical study lies in identifying potential group structure over the weather stations to achieve dimension reduction in the subsequent network VAR model estimation. With the clustering algorithm in Section (ref), we obtain two estimated groups, i.e., $\widehat{K}=2$, for both models. Table (ref) reports the estimated group membership which is very similar between the two models. For model (i), all the stations of Group 1 are in Northern Ireland and Wales, while those of Group 2 are in England and Scotland. This indicates that Northern Ireland and Wales often have common patterns in terms of temperature change whereas the weather stations in the island of Great Britain usually have similar temperature recordings. It is noteworthy that, although each weather station in model (i) contains the minimum and maximum temperatures (as two nodes), both nodes are classified into the same group (thus we only report the station name in Table (ref)). For model (ii), three weather stations in the coastal towns (Eastbourne, Tiree and Whitby) move from Group 2 to Group 1. The estimated time-varying coefficient functions are plotted in the online supplement LPTW24.
To further show the necessity of accounting for the group structure and smooth structural changes, we compare the out-of-sample prediction performance among the following three methods: fully heterogeneous time-varying network VAR model YSM24; the proposed grouped time-varying network VAR model; and the grouped network VAR model with constant coefficients ZXF23, denoted by "fully heterogeneous", "grouped+TV" and "grouped+linear", respectively, in Table (ref). We leave the last $T_{\sf pre}$ observations out for prediction, where $T_{\sf pre}=12, 24$ and $36$, corresponding to one, two and three years, respectively. For any time point $t_\bullet$ in the prediction (or test) period, we use the observations over $t=1,\cdots, t_\bullet-1$ to get the estimates of time-varying or constant coefficients, which are subsequently used to forecast the value of $X_{t_\bullet}$. We conduct this expanding-window one-step ahead forecasting exercise for all the three models and report their out-of-sample RMSE in Table (ref). It follows from Table (ref) that our proposed model produces the most accurate out-of-sample forecasting results with much smaller RMSE than the other two competing methods.
In this paper we have introduced a general nonlinear network VAR model for high-dimensional time series, where the momentum and network spillover effects are allowed to change over time and nodes. To achieve dimension reduction and obtain satisfactory estimation convergence rates, we impose a latent group structure on time-varying coefficients in the heterogenous network VAR model. The unknown group number is determined by an easy-to-implement criterion whereas the group membership is estimated by the agglomerative clustering algorithm with the nonparametrically estimated distance matrix. Theorems (ref) and (ref) show that the developed methodology consistently estimates the latent group structure. To further improve the convergence rates of the time-varying coefficient estimation, we have proposed a post-grouping local linear smoothing to estimate the group-specific time-varying momentum and network effects. In addition, we further extend the model, methodology and theory to allow for structural breaks in either the group structure or group-specific coefficient functions. The simulation study demonstrates that (i) the developed method can accurately estimate the latent group structure in finite samples; (ii) the post-grouping local linear estimation significantly outperforms the naive heterogenous estimation which ignores the latent structure; and (iii) the developed two-stage method can accurately locate the break point and estimate the group structure before and after the break. The empirical study of the UK temperature time series shows that there exist two groups over the 37 UK weather stations and our proposed method has better out-of-sample prediction performance than the other two competing methods which ignore either the grouped or time-varying feature in network VAR model building.
The authors thank the Editor, an Associate Editor and two anonymous referees for their constructive comments, which helped to substantially improve an earlier version of the article. The first author is supported by the Leverhulme Research Fellowship (RF-2023-396). The second author is supported by the Australian Research Council Discovery Grants (DP210100476).