Estimating Time-Varying Networks for High-Dimensional Time Series

Estimating Time-Varying Networks for High-Dimensional Time Series

\centerline{\bf Abstract}

We explore time-varying networks for high-dimensional locally stationary time series, using the large VAR model framework with both the transition and (error) precision matrices evolving smoothly over time. Two types of time-varying graphs are investigated: one containing directed edges of Granger causality linkages, and the other containing undirected edges of partial correlation linkages. Under the sparse structural assumption, we propose a penalised local linear method with time-varying weighted group LASSO to jointly estimate the transition matrices and identify their significant entries, and a time-varying CLIME method to estimate the precision matrices. The estimated transition and precision matrices are then used to determine the time-varying network structures. Under some mild conditions, we derive the theoretical properties of the proposed estimates including the consistency and oracle properties. In addition, we extend the methodology and theory to cover highly-correlated large-scale time series, for which the sparsity assumption becomes invalid and we allow for common factors before estimating the factor-adjusted time-varying networks. We provide extensive simulation studies and an empirical application to a large U.S. macroeconomic dataset to illustrate the finite-sample performance of our methods.

\noindentKeywords: CLIME, factor model, Granger causality, LASSO, local linear smoothing, partial correlation, time-varying network, VAR.

Introduction

\setcounter{equation}{0}

In recent years, the network analysis has become an effective tool to explore inter-connections among a large number of variables, with applications to various disciplines such as: epidemiology, economics, finance, and social networks N02, BKT13, DY14, DY15, HSS14, S17, BB19, ZCLW19. The so-called graphical model is commonly used in the network analysis to visualise the connectedness of a large panel with vertices representing variables in the panel and the presence of an edge indicating appropriate (conditional) dependence between the variables. In the past decades, most of the existing literature on statistical estimation and inference of network data limits attention to the static network, which is assumed to be invariant over time YL07, FFW09, LW13, BSM15, ZLWL22. However, such an assumption may be too restrictive and often fails in practical applications where the underlying data generating mechanism is dynamic. There have been some attempts in the recent literature to relax the static network assumption, allowing the connectivity structure to exhibit time-varying features. For example, KSAX10 and ZLW10 study dynamic network models with smooth time-varying structural changes; whereas WYR21 consider change-point detection and estimation in dynamic networks. However, most of the aforementioned literature typically assumes that the network data are independent, which often becomes invalid in practice. We aim to relax this restrictive assumption and model large-scale network data under a general temporal dependence structure.

Vector autoregression (VAR) is a fundamental modelling tool for multivariate time series data Lu06. In recent years, there has been increasing interest in extending the finite-dimensional VAR to the high-dimensional setting. Under appropriate sparsity restrictions on the transition (or autoregressive coefficient) matrices, various regularised methods have been proposed to estimate high-dimensional VAR models and identify non-zero entries in the transition matrices BM15, HLL15, KC15, DZZ16. ZPLLW17 introduce a network VAR model by incorporating the adjacency matrix to capture the network effect and estimate the model via ordinary least squares. More recently, CFZ20 and MPS22 further study high-dimensional VAR and network VAR with latent common factors, allowing strong cross-sectional dependence in large panel time series. The methodology and theory developed in these papers heavily rely on the stationarity assumption with both transition and volatility matrices being time-invariant.

The stable VAR model cannot capture smooth structural changes and breaks in the underlying data generating process, two typical dynamic features in time series data collected over a long time span. To address this problem, DQC17 consider a time-varying VAR model for high-dimensional time series (allowing the number of variables to diverge at a sub-exponential rate of the sample size), and estimate the time-varying transition matrices by combining the kernel smoothing with $\ell_{1}$-regularisation, whereas SS22 simultaneously detect breaks and estimate transition matrices in high-dimensional VAR via a three-stage procedure using the total variation penalty. XCW20 detect structural breaks and estimate smooth changes (between breaks) in the covariance and precision matrices of high-dimensional time series (covering VAR as a special case). In the present paper, we aim to jointly estimate the time-varying transition and precision matrices in the high-dimensional sparse VAR under the local stationarity framework. Motivated by the stable network time series analysis in BB19, we use the estimated transition and precision matrices to further construct two time-varying networks: one containing directed edges of Granger causality linkages, and the other containing undirected edges of partial correlation linkages.

The proposed time-varying network via VAR is naturally connected to the locally stationary models, which have been systematically studied in the literature for low-dimensional time series. D97 is among the first to introduce a locally stationary time series model via a time-varying spectral representation. DS06 study a time-varying ARCH model and propose a kernel-weighted quasi-maximum likelihood estimation method. HL10 further consider a time-varying version of GARCH model and introduce a semiparametric method to estimate both the parametric and nonparametric components involved. Vo12 and ZW12 study nonparametric kernel-based estimation and inference in a general class of locally stationary time series. KL12 extend the locally stationary model framework to the diffusion process. YGP20 develop a kernel estimation method and theory for time-varying vector moving average models. The present paper complements the locally stationary time series literature by further exploring the high-dimensional dynamic network structure.

We study the time-varying VAR and network models for large-scale time series, allowing the number of variables to be much larger than the time series length. Under the sparsity assumption on the transition and precision matrices with smooth structural changes, we introduce a three-stage estimation procedure: (i) preliminary local linear estimation of the transition matrices and their derivatives with time-varying LASSO; (ii) joint local linear estimation and feature selection of the time-varying transition matrices with weighted group LASSO; (iii) estimation of the precision matrix via time-varying CLIME. To guarantee the oracle property, the weights of LASSO in the second estimation stage are constructed via a local linear approximation to the SCAD penalty ZL08 using the consistent preliminary estimates obtained in the first stage. Our penalised estimation methodology for the time-varying transition matrices is connected to various nonparametric screening and shrinkage methods developed for high-dimensional functional-coefficient models WX09, L12, FMD14, LLW14, LKZ15, whereas the time-varying CLIME is a natural extension of the conventional CLIME for static precision matrix estimation CLL11. The theoretical properties of the techniques developed in the aforementioned literature (such as the oracle property and minimax optimal convergence rates) rely on the independent data assumption. Extension of the methodology and theory to the high-dimensional locally stationary time series is non-trivial, requiring new technical tools such as the concentration inequality for time-varying VAR. Under some regularity conditions, we show that the proposed local linear estimates with weighted group LASSO equal to the infeasible oracle estimates with prior information on the significant entries of time-varying transition matrices, and the precision matrix estimate with time-varying CLIME is uniformly consistent with sensible convergence rates under various matrix norms. The estimated transition matrices are used to consistently estimate the uniform network structure with directed Granger causality linkages, whereas the estimated precision matrix is used to construct the network structure with undirected partial correlation linkages.

We further consider highly-correlated large-scale time series, for which the sparsity model assumption is no longer valid in which case the methodology and theory need to be substantially modified. The approximate factor model CR83 or its time-varying version SW17 is employed to accommodate the strong cross-sectional dependence among a large number of time series. In particular, we assume that the high-dimensional idiosyncratic error process in the approximate factor model satisfies the time-varying VAR structure with the sparsity restriction imposed on its transition and precision matrices. The latent common and idiosyncratic components need to be estimated consistently. With the approximated idiosyncratic error vectors, the penalised local linear estimation method with weighted group LASSO and time-varying CLIME are applied to estimate the time-varying transition and precision matrices. Subsequently, the factor-adjusted time-varying network estimates with directed Granger causality and undirected partial correlation linkages are obtained. Our paper thus substantially extends the recent work on the factor-adjusted stable VAR model estimation FMM21, BCO22, KM22.

Our simulation studies demonstrate that the proposed methodology can accurately estimate the time-varying Granger and partial correlation networks when the number of time series variables is comparable to the sample size. In particular, for the time-varying transition matrix estimation, the penalised local linear method with weighted group LASSO outperforms the conventional local linear method (which often fails in the high-dimensional time series setting) and produces numerical results similar to those of the oracle estimation. For the time-varying error precision matrix estimation, the numerical performance of the proposed time-varying CLIME is comparable to that of the time-varying graphical LASSO. We further apply the developed methodology to the FRED-MD macroeconomic dataset and estimate both the Granger causality and partial correlation networks via the proposed time-varying VAR model.

The rest of the paper is organised as follows. Section (ref) introduces the time-varying VAR and network model structures. Section (ref) presents the estimation procedures for the time-varying transition and precision matrices and Section (ref) gives the asymptotic properties of the developed estimates. Section (ref) considers the factor-adjusted time-varying VAR model and network estimation. Sections (ref) and (ref) report simulation studies and an empirical application, respectively. Section (ref) concludes the paper. A supplemental document contains proofs of the main theorems, some technical lemmas with proofs, verification of a key assumption and discussions on tuning parameter selection. Throughout the paper, we let $\vert\cdot\vert_{0}$, $\vert \cdot\vert_{1}$, $\Vert\cdot\Vert$ and $\vert\cdot\vert_{\max}$ denote the $L_{0}$, $L_{1}$, $L_{2}$ (Euclidean) and maximum norms of a vector, respectively. Let ${\mathbf{I}}_{d}$ and ${\mathbf{O}}_{d\times d}$ be a $d\times d$ identity matrix and null matrix, respectively. For a $d\times d$ matrix ${\mathbf{W}}=(w_{ij})_{d\times d}$, we let $\Vert{\mathbf{W}} \Vert=\lambda_{\max}^{1/2}\left( {\mathbf{W}}^{^{\intercal}}{\mathbf{W} }\right) $ be the operator norm, $\Vert{\mathbf{W}}\Vert_{F}=\left[ \mathsf{Tr}\left( {\mathbf{W}}^{^{\intercal}}{\mathbf{W}}\right) \right] ^{1/2}$ the Frobenius norm, $\Vert{\mathbf{W}}\Vert_{1}=\max_{1\leq j\leq d}\sum_{i=1}^{d} |w_{ij}|$, $\Vert{\mathbf{W}}\Vert_{\max}=\max_{1\leq i\leq d}\max_{1\leq j\leq d} |w_{ij}|$, and $\vert{\mathbf{W}}\vert_{1}=\sum _{i=1}^{d}\sum_{j=1}^{d} |w_{ij}|$, where $\lambda_{\max}(\cdot)$ is the maximum eigenvalue of a matrix and $\mathsf{Tr}(\cdot)$ is the trace. Denote the determinant of a square matrix as $\mathsf{det}(\cdot)$. Let $a_{n}\sim b_{n}$, $a_{n}\propto b_{n}$ and $a_{n}\gg b_{n}$ denote that $a_{n} /b_{n}\rightarrow1$, $0<\underline{c}\leq a_{n}/b_{n}\leq\overline{c}<\infty$ and $b_{n}/a_{n}\rightarrow0$, respectively.

Time-varying VAR and network models

\setcounter{equation}{0}

In this section, we first introduce a locally stationary VAR model with time-varying transition and precision matrices, and then define two types of time-varying network structures with Granger causality and partial correlation linkages, respectively. Section (ref) will further generalise them to the factor-adjusted time-varying VAR and network setting.

Time-varying VAR models

Suppose that $(X_{t}:t=1,\mathcal{\ldots},n)$ with $X_{t}=(x_{t,1} ,\mathcal{\ldots},x_{t,d})^{^{\intercal}}$ is a sequence of $d$-dimensional random vectors generated by a time-varying VAR model of order $p$:

where ${\mathbf{A}}_{t,k}={\mathbf{A}}_{k}(t/n)$, $k=1,\mathcal{\ldots},p$, are $d\times d$ time-varying transition matrices with each entry being a smooth deterministic function of scaled times, ${\boldsymbol{\Sigma}} _{t}={\boldsymbol{\Sigma}}(t/n)$ is a $d\times d$ time-varying volatility matrix, and $(\varepsilon_{t})$ is a sequence of independent and identically distributed (i.i.d.) $d$-dimensional random vectors with zero mean and identity covariance matrix. Define ${\boldsymbol{\Omega}}_{t} ={\boldsymbol{\Omega}}(t/n)$ as the inverse of ${\boldsymbol{\Sigma}}_{t}$, the time-varying precision matrix. We consider the ultra large time series setting, i.e., the dimension $d$ is allowed to diverge at an exponential rate of the sample size $n$. The time-varying VAR model ((ref)) is a natural extension of the finite-dimensional time-varying VAR to high-dimensional time series. If ${\boldsymbol{\Sigma}}_{t}$ is replaced by a time-invariant covariance matrix, ((ref)) becomes the same model as that considered by DQC17. Furthermore, when both ${\mathbf{A}}_{t,k}$, $k=1,\mathcal{\ldots},p$, and ${\boldsymbol{\Sigma}}_{t}$ are time-invariant constant matrices, ((ref)) becomes the high-dimensional stable VAR:

which has been extensively studied in the recent literature BM15, HLL15, KC15, BB19, LZ21. Throughout the paper, we assume that the following conditions are satisfied.

The first condition in Assumption (ref)(i) is a natural extension of the stability assumption imposed on the constant transition matrices Lu06, indicating that the time-varying VAR process is locally stationary/stable and leading to the following Wold representation

with the coefficient matrices ${\boldsymbol{\Phi}}_{t,k}$ being absolutely summable (in appropriate matrix norm). For example, when $p=1$, we have ${\boldsymbol{\Phi}}_{t,0}={\mathbf{I}}_{d}$ and ${\boldsymbol{\Phi}} _{t,k}=\Pi_{j=1}^{k}{\mathbf{A}}_{t-j+1,1}$ for $k\geq1$. Assume that, for $k$ sufficiently large,

where $C_{1}$ is a positive constant and $0<\rho<1$. A similar assumption can be found in DQC17. In some special model settings, ((ref)) may be violated, and we refer the interested readers to the discussions in BM15 and LZ21. In fact, the condition ((ref)) may be removed by imposing some high-level conditions (e.g., the sub-Gaussian condition on $x_{t,i}$ proved in Lemma B.1). The smoothness conditions in Assumption (ref)(i)(ii) are common in kernel-based local estimation method and theory. The sub-Gaussian moment condition in Assumption (ref)(iii) is not uncommon in the literature of high-dimensional feature selection and covariance/precision matrix estimation W19, and is weaker than the Gaussian assumption frequently used in the high-dimensional VAR literature BM15, KC15.

Time-varying network structures

Write ${\mathbf{A}}_{t,k}=\left( a_{k,ij|t}\right) _{d\times d}$, ${\boldsymbol{\Omega}}_{t}=\left( \omega_{ij|t}\right) _{d\times d}$, ${\mathbf{A}}_{k}(\tau)=\left( a_{k,ij}(\tau)\right) _{d\times d}$ and ${\boldsymbol{\Omega}}(\tau)=\left( \omega_{ij}(\tau)\right) _{d\times d}$, where $1\leq t\leq n$ and $0\leq\tau\leq1$. We define the network structure via a time-varying graph ${\mathbb{G}}_{t}=({\mathbb{V}},{\mathbb{E}}_{t})$, where ${\mathbb{V}}=\{1,2,\mathcal{\ldots},d\}$ denotes a set of vertices, and ${\mathbb{E}}_{t}=\left\{ (i,j)\in{\mathbb{V}}\times{\mathbb{V}} :\ c_{ij|t}\neq0,\ i\neq j\right\} $ denotes a time-varying set of edges. The choice of $c_{ij|t}$ is determined by the definition of linkage. The construction of ${\mathbb{G}}_{t}$ is similar to that in KSAX10 and ZLW10 for independent network data. Following the stable network analysis in BB19 and BCO22, we next consider two types of time-varying linkages: the directed Granger causality linkage and undirected partial correlation linkage.

The definition of Granger causality is first introduced by G69 to investigate the causal relations in small economic time series systems. In the context of stable VAR (with order $p$), we say that $x_{t,j}$ Granger causes $x_{t,i}$ if there exists $k\in\{1,2,\mathcal{\ldots},p\}$ such that $x_{t-k,j}$ improves predictability of $x_{t,i}$ by reducing the forecasting error. It is a natural idea to use the stable transition matrices ${\mathbf{A}}_{k}=\left( a_{k,ij}\right) _{d\times d}$ in ((ref)) to determine the Granger causality structure, i.e., if there exists at least one $k$ such that $a_{k,ij}\neq0$, then $x_{t,j}$ Granger causes $x_{t,i}$. We may extend the stable Granger causality structure to a more general time-varying version using ((ref)). At a given time point $t$, we say that lags of $x_{t,j}$ Granger cause $x_{t,i}$ if there exists at least one $k$ such that $a_{k,ij|t}\neq0$. Hence, for given $\tau\in(0,1)$, we define the time-varying local graph ${\mathbb{G}}_{\tau}^{G}=\left( {\mathbb{V}},{\mathbb{E}}_{\tau }^{G}\right) $ with

The partial correlation is a commonly-used conditional dependence measure for network time series. We next extend it to the time-varying setting using ${\boldsymbol{\Omega}}_{t}={\boldsymbol{\Omega}}(t/n)$ in ((ref)). Let $\rho_{ij|t}=\mathsf{cor}(e_{t,i}, e_{t,j} | e_{t,k}, k\neq i,j)$ be the time-varying (contemporaneous) partial correlation between the innovations $e_{t,i}$ and $e_{t,j}$, where $e_{t,i}$ is the $i$-th element of $e_{t}$. Following D72, we may show that $\rho_{ij|t}\neq0$ is equivalent to $\omega_{ij|t}\neq0$ for $i\neq j$. Hence, we can construct the set of edges by collecting the index pairs of the non-zero entries in the time-varying precision matrix. For $\tau\in(0, 1)$, define the local graph ${\mathbb{G} }_{\tau}^{P}=\left( {\mathbb{V}}, {\mathbb{E}}_{\tau}^{P}\right) $ with

In practice, the primary interest often lies in the full network structures over the entire time interval. This requires the construction of a uniform version of ${\mathbb{G}}_{\tau}^{G}$ and ${\mathbb{G}}_{\tau}^{P}$. Denote the uniform graphs by ${\mathbb{G}}^{G}=\left( {\mathbb{V}},{\mathbb{E}} ^{G}\right) $ and ${\mathbb{G}}^{P}=\left( {\mathbb{V}},{\mathbb{E}} ^{P}\right) $, with

and

It is easy to verify that ${\mathbb{E}}_{\tau}^{G}\subset{\mathbb{E}}^{G}$ and ${\mathbb{E}}_{\tau}^{P}\subset{\mathbb{E}}^{P}$ for any $\tau\in(0,1)$. Section (ref) below defines the discrete versions of the above uniform networks and provide their estimates.

Methodology

\setcounter{equation}{0}

Let $A_{k,i}^{^{\intercal}}(\cdot)$ and $C_{i}^{^{\intercal}}(\cdot)$ be the $i$-th row of ${\mathbf{A}}_{k}(\cdot)$ and ${\boldsymbol{\Omega}} ^{-1/2}(\cdot)$, respectively,

and $\tau_{t}=t/n$. The time-varying VAR model ((ref)) can be equivalently written as

which is a high-dimensional time-varying coefficient autoregressive model with a scalar response and $pd$ candidate predictors for each $i$. As the dimension of the predictors is allowed to be ultra large, we need to impose an appropriate sparsity restriction on the vector of time-varying parameters ${\boldsymbol{\alpha}}_{i\bullet}(\cdot)$ to limit the number of its significant elements. High-dimensional varying-coefficient models have been systematically studied in the literature and various nonparametric screening and shrinkage methods have been proposed to select the significant covariates, estimate the coefficient functions and identify the model structure under the independent data assumption WLH08, WX09, L12, CHLP14, FMD14, LLW14, LKZ15. In this section, under the high-dimensional locally stationary time series framework, we propose a three-stage procedure to estimate the Granger causality and partial correlation network structures: (i) first obtain preliminary local linear estimates of ${\boldsymbol{\alpha}}_{i\bullet}(\cdot)$ (and its derivatives) using time-varying LASSO, which serves as a first-stage screening of the predictors in ${\mathbf{X}}_{t-1}$; (ii) conduct local linear estimation and feature selection using weighted group LASSO, where the weights are constructed via a local linear approximation to the SCAD penalty using the preliminary estimates of ${\boldsymbol{\alpha}}_{i\bullet}(\cdot)$ from Stage (i); (iii) estimate the error precision matrix ${\boldsymbol{\Omega}}(\cdot)$ via the time-varying CLIME method. The estimated transition and precision matrices are finally used to construct the uniform network structures.

Preliminary time-varying LASSO estimation

For $\tau\in(0,1)$, under the smoothness condition on the transition matrices in Assumption (ref)(i), we have the following local linear approximation to ${\boldsymbol{\alpha}}_{i\bullet}(\tau_{t})$: \[ {\boldsymbol{\alpha}}_{i\bullet}(\tau_{t})\approx{\boldsymbol{\alpha} }_{i\bullet}(\tau)+{\boldsymbol{\alpha}}_{i\bullet}^{\prime}(\tau)(\tau _{t}-\tau),\ \ i=1,\mathcal{\ldots},d, \] when $\tau_{t}$ falls within a small neighbourhood of $\tau$, where ${\boldsymbol{\alpha}}_{i\bullet}^{\prime}(\cdot)$ is a $(pd)$-dimensional vector of the first-order derivatives of the elements in ${\boldsymbol{\alpha }}_{i\bullet}(\cdot)$. Hence, for each $i\in\{1,2,\mathcal{\ldots},d\}$ and a given $\tau\in(0,1)$, we define the following local linear objective function FG96:

where $K_{h}(\cdot)=\frac{1}{h}K(\cdot/h)$ with $K(\cdot)$ being a kernel function and $h$ being a bandwidth or smoothing parameter. The estimates of ${\boldsymbol{\alpha}}_{i\bullet}(\tau)$ and ${\boldsymbol{\alpha}}_{i\bullet }^{\prime}(\tau)$ are obtained by minimising $\mathcal{L}_{i} ({\boldsymbol{\alpha}},{\boldsymbol{\beta}}\ |\ \tau)$ with respect to ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$. However, this local linear estimation is only feasible when the dimension of the predictors is fixed or significantly smaller than the sample size $n$ C07, LCG11. In our high-dimensional setting, as the number of predictors may exceed $n$, it is challenging to obtain satisfactory estimation by directly minimising $\mathcal{L}_{i}({\boldsymbol{\alpha}},{\boldsymbol{\beta}}\ |\ \tau)$. To address this issue, we assume that the number of significant components in ${\boldsymbol{\alpha}}_{i\bullet}(\tau)$ is much smaller than $n$ and then incorporate a LASSO penalty term in the local linear objective function ((ref)).

The LASSO estimation was first introduced by T96 in the context of linear regression and has become one of the most commonly-used tools in high-dimensional variable and feature selection. We next adopt a time-varying version of the LASSO estimation. Define

where $\lambda_{1}$ is a tuning parameter. Let $\widetilde{\boldsymbol{\alpha }}_{i\bullet}(\tau)$ and $\widetilde{\boldsymbol{\alpha}}_{i\bullet}^{\prime }(\tau)$ be the solution to the minimisation of $\mathcal{L}_{i}^{\ast }({\boldsymbol{\alpha}}, {\boldsymbol{\beta}}\ |\ \tau)$ with respect to ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$. We call them the preliminary time-varying LASSO estimates. This LASSO estimation may not accurately identify the true significant predictors, but can remove a large number of irrelevant predictors and hence, serves as a preliminary screening step. Furthermore, the first-stage estimates will be used to construct weights in the weighted group LASSO in the second stage to more precisely estimate the time-varying parameters and accurately select the significant predictors.

Penalised local linear estimation with weighted group LASSO

In order to estimate the uniform Granger causality network, we next introduce a global penalised method to simultaneously estimate the time-varying parameters at $\tau_{t}$, $t=1,\mathcal{\ldots},n$, and identify the non-zero index sets $\mathcal{J}_{i}=\bigcup_{t=1}^{n}\mathcal{J}_{i}(\tau_{t})$ and $\mathcal{J}_{i}^{\prime}=\bigcup_{t=1}^{n}\mathcal{J}_{i}^{\prime}(\tau_{t} )$, where \[ \mathcal{J}_{i}(\tau)=\left\{ 1\leq j\leq pd:\ \alpha_{i,j}(\tau )\neq0\right\} \ \ \mathrm{and}\ \ \mathcal{J}_{i}^{\prime}(\tau)=\left\{ 1\leq j\leq pd:\ \alpha_{i,j}^{\prime}(\tau)\neq0\right\} \] with $\alpha_{i,j}(\cdot)$ and $\alpha_{i,j}^{\prime}(\cdot)$ being the $j$-th element of ${\boldsymbol{\alpha}}_{i\bullet}(\cdot)$ and ${\boldsymbol{\alpha }}_{i\bullet}^{\prime}(\cdot)$, respectively. For each $i$, note that identifying the zero elements in ${\boldsymbol{\alpha}}_{i\bullet}^{\prime }(\tau_{t})$ (uniformly over $t$) is equivalent to identifying the indices $j$, $1\leq j\leq pd$, such that $D_{i,j}=0$, where \[ D_{i,j}^{2}=\sum\limits_{t=1}^{n}\left[ \alpha_{i,j}(\tau_{t})-\frac{1} {n}\sum\limits_{s=1}^{n}\alpha_{i,j}(\tau_{s})\right] ^{2}. \] In practice, $D_{i,j}^{2}$ can be estimated by \[ \widetilde{D}_{i,j}^{2}=\sum\limits_{t=1}^{n}\left[ \widetilde{\alpha} _{i,j}(\tau_{t})-\frac{1}{n}\sum\limits_{s=1}^{n}\widetilde{\alpha}_{i,j} (\tau_{s})\right] ^{2}, \] using the preliminary time-varying LASSO estimates $\widetilde{\alpha} _{i,j}(\tau_{t})$, $t=1,\ldots,n$. Let ${\mathbf{A}}=({\boldsymbol{\alpha} }_{\bullet1},\mathcal{\ldots},{\boldsymbol{\alpha}}_{\bullet n})^{^{\intercal }}$ with ${\boldsymbol{\alpha}}_{\bullet t}=(\alpha_{1|t},\mathcal{\ldots },\alpha_{pd|t})^{^{\intercal}}$, and ${\mathbf{B}}=({\boldsymbol{\beta} }_{\bullet1},\mathcal{\ldots},{\boldsymbol{\beta}}_{\bullet n})^{^{\intercal} }$ with ${\boldsymbol{\beta}}_{\bullet t}=(\beta_{1|t},\mathcal{\ldots} ,\beta_{pd|t})^{^{\intercal}}$. We define a global version of the penalised objective function with weighted group LASSO:

where \[ \widetilde{\boldsymbol{\alpha}}_{i,j}=\left[ \widetilde{\alpha}_{i,j} (\tau_{1}),\mathcal{\ldots},\widetilde{\alpha}_{i,j}(\tau_{n})\right] ^{^{\intercal}},\ \ {\boldsymbol{\alpha}}_{j}=\left( \alpha_{j|1} ,\mathcal{\ldots},\alpha_{j|n}\right) ^{^{\intercal}},\ \ {\boldsymbol{\beta }}_{j}=\left( \beta_{j|1},\mathcal{\ldots},\beta_{j|n}\right) ^{^{\intercal }}, \] while $\lambda_{2}$ is a tuning parameter and $p_{\lambda}^{\prime}(\cdot)$ is the derivative of the SCAD penalty function: \[ p_{\lambda}^{\prime}(z)=\lambda\left[ I(z\leq\lambda)+\frac{(a_{0} \lambda-z)_{+}}{(a_{0}-1)\lambda}I(z>\lambda)\right] , \] with $a_{0}=3.7$ as suggested in FL01 and $I(\cdot)$ being the indicator function. The penalty terms in ((ref)) are motivated by the local linear approximation to the SCAD penalty function ZL08. The terms $p_{\lambda_{2}}^{\prime}\left( \left\Vert \widetilde {\boldsymbol{\alpha}}_{i,j}\right\Vert \right) $ and $p_{\lambda_{2}} ^{\prime}\left( \widetilde{D}_{i,j}\right) $ in ((ref)) serve as the weights for the group LASSO, and their values are determined by the preliminary estimates in Section (ref), i.e., the corresponding weight is heavy when $\left\Vert \widetilde{\boldsymbol{\alpha}}_{i,j}\right\Vert $ or $\widetilde{D}_{i,j}$ is close to zero, whereas it is light or equal to zero when $\left\Vert \widetilde{\boldsymbol{\alpha}}_{i,j}\right\Vert $ or $\widetilde{D}_{i,j}$ is large. An advantage of using $\widetilde{D}_{i,j}$ in the second penalty term over the $L_{2}$-norm of $\widetilde {\boldsymbol{\alpha}}_{j}^{\prime}=\left[ \widetilde{\alpha}_{i,j}^{\prime }(\tau_{1}),\mathcal{\ldots},\widetilde{\alpha}_{i,j}^{\prime}(\tau _{n})\right] ^{^{\intercal}}$ is that the estimates of the time-varying parameters involved in $\widetilde{D}_{i,j}$ often perform more stably than their derivative counterparts.

Let $\widehat{\mathbf{A}}_{i}$ and $\widehat{\mathbf{B}}_{i}$ be the minimiser of $\mathcal{Q}_{i}({\mathbf{A}},{\mathbf{B}})$ with respect to ${\mathbf{A}}$ and ${\mathbf{B}}$, where

The index set $\mathcal{J}_{i}$ is estimated by $\widehat{\mathcal{J}} _{i}=\left\{ j:\ \widehat{\boldsymbol{\alpha}}_{i,j}\neq{\mathbf{0}} _{n}\right\} $, and $\mathcal{J}_{i}^{\prime}$ is estimated by $\widehat {\mathcal{J}}_{i}^{\prime}=\left\{ j:\ \widehat{\boldsymbol{\alpha}} _{i,j}^{\prime}\neq{\mathbf{0}}_{n}\right\} $, where ${\mathbf{0}}_{k}$ is a $k$-dimensional vector of zeros. A similar shrinkage estimation method is used by LKZ15 and CLWZ21 to identify a high-dimensional semi-varying coefficient model structure for independent data. So far as we know, there is no work on such a penalised technique and its relevant theory for high-dimensional locally stationary time series data.

Estimation of the time-varying precision matrix

In this section, we study the estimation of ${\boldsymbol{\Omega}}(\cdot)$ in model ((ref)), which is crucial to uncover the time-varying and uniform network structures of partial correlations. Estimation of large static precision matrices has been extensively studied under the sparsity assumption, and various estimation techniques, such as the penalised likelihood, graphical Danzig selector and CLIME, have been proposed in the literature LF09, Y10, CLL11. XCW20 further introduce a time-varying CLIME method for high-dimensional locally stationary time series which are observable. Note that in this paper, ${\boldsymbol{\Omega}}(\cdot)$ is the time-varying precision matrix for the high-dimensional unobservable error vector $e_{t}$ and hence, its estimation requires substantial modification of the time-varying CLIME methodology and theory.

With $\widehat{\boldsymbol{\alpha}}_{i\bullet}(\cdot)$, $i=1,\mathcal{\ldots },d$, from Section (ref), we can then extract estimates of the time-varying transition matrices, denoted by $\widehat{\mathbf{A}}_{k} (\tau_{t})$, $t=1,\mathcal{\ldots},n$, $k=1,\mathcal{\ldots},p$, and approximate $e_{t}$ by

The approximation accuracy depends on the uniform prediction rates of the time-varying weighted group LASSO estimates. In order to apply the time-varying CLIME, we assume that ${\boldsymbol{\Omega}}(\cdot)$ satisfies a uniform sparsity assumption, a natural extension of the classic sparsity assumption to the locally stationary time series setting. Specifically, we assume $\left\{ {\boldsymbol{\Omega}}(\tau):0\leq\tau\leq1\right\} \in\mathcal{S}(q,\xi_{d})$, where {

} where $0\leq q<1$, \textquotedblleft${\mathbf{W}}\succ0$" denotes that ${\mathbf{W}}$ is positive definite, and $C_{2}$ is a bounded positive constant. Define

where the weight function $\varpi_{n,t}(\cdot)$ is constructed via the local linear smoothing: \[ \varpi_{n,t}(\tau)=K\left( \frac{\tau_{t}-\tau}{b}\right) s_{n,2} (\tau)-K_{1}\left( \frac{\tau_{t}-\tau}{b}\right) s_{n,1}(\tau), \] in which $s_{n,j}(\tau)=\sum_{t=1}^{n}K_{j}\left( \frac{\tau_{t}-\tau} {b}\right) $, $K_{j}(x)=x^{j}K(x)$, and $b$ is a bandwidth. With the uniform sparsity assumption ((ref)), we estimate ${\boldsymbol{\Omega}}(\tau)$ via the time-varying CLIME method:

where $\lambda_{3}$ is a tuning parameter. As the underlying time-varying precision matrix is symmetric, the matrix estimate obtained from ((ref)) needs to be symmetrised to obtain the final estimate, denoted as $\widehat{\boldsymbol{\Omega}}(\tau)=\left[ \widehat{\omega}_{ij} (\tau)\right] _{d\times d}$, where

Estimation of uniform time-varying networks

In practice, when the sample size $n$ is sufficiently large, it is often sensible to approximate the uniform edge sets, ${\mathbb{E}}^{G}$ and ${\mathbb{E}}^{P}$, by the following discrete versions:

and

Hence, we next estimate ${\mathbb{E}}_{n}^{G}$ and ${\mathbb{E}}_{n}^{P}$ instead of ${\mathbb{E}}^{G}$ and ${\mathbb{E}}^{P}$. With the time-varying transition and precision matrix estimates in Sections (ref) and (ref), we can estimate ${\mathbb{E}}_{n}^{G}$ by

where $\widehat{a}_{k,ij}(\tau_{t})$ is the $(i,j)$-entry of $\widehat {\mathbf{A}}_{k}(\tau_{t})$, and estimate ${\mathbb{E}}_{n}^{P}$ by

where $\lambda_{3}$ is the tuning parameter used in the time-varying CLIME.

Main theoretical results

\setcounter{equation}{0}

To ease the notational burden, throughout this section, we focus on the time-varying VAR(1) model:

where ${\mathbf{A}}(\tau)=\left[ \alpha_{ij}(\tau)\right] _{d\times d}$. For a general time-varying VAR($p$) model ((ref)), it can be equivalently re-written as a $(pd)$-dimensional VAR(1) model as follows: \[ {\mathbf{X}}_{t}={\mathbf{A}}_{t}^{\ast}{\mathbf{X}}_{t-1}+{\mathbf{e}}_{t}, \] where ${\mathbf{X}}_{t}$ is defined in ((ref)), ${\mathbf{e}} _{t}=\left( e_{t}^{^{\intercal}},0_{d}^{^{\intercal}},\mathcal{\ldots} ,0_{d}^{^{\intercal}}\right) ^{^{\intercal}}$, and ${\mathbf{A}}_{t}^{\ast}$ is a $(pd)\times(pd)$ time-varying transition matrix: \[ {\mathbf{A}}_{t}^{\ast}=\left(

\right) . \]

Uniform consistency of the time-varying LASSO estimates

Define {

} and \[ \mathcal{B}_{i}(\tau)=\left\{ \left( u_{1}^{^{\intercal}}, u_{2}^{^{\intercal }}\right) ^{^{\intercal}}: \|u_{1}\|^{2}+\|u_{2}\|^{2}=1,\ \sum_{j=1} ^{d}\left( |u_{1,j}|+|u_{2,j}|\right) \leq3 \left( \sum_{j\in\mathcal{J} _{i}(\tau)}|u_{1,j}|+\sum_{j\in\mathcal{J}_{i}^{\prime}(\tau)}|u_{2,j}|\right) \right\} , \] where $\mathcal{J}_{i}(\tau)$ and $\mathcal{J}_{i}^{\prime}(\tau)$ are defined as in Section (ref) but with $p=1$. To derive the uniform consistency property of the preliminary time-varying LASSO estimates defined in Section (ref), we need the following assumptions, some of which may be weakened at the cost of lengthier proofs.

Assumption (ref)(i) is a mild restriction which can be satisfied by some commonly-used kernels such as the uniform kernel and the Epanechnikov kernel. The compact support assumption on the kernel function is not essential and can be replaced by appropriate tail conditions. The bandwidth conditions in Assumption (ref)(ii) are crucial for deriving the uniform convergence properties of the kernel-based quantities. When $s$ is bounded and $d$ diverges at a polynomial rate of $n$, the conditions can be simplified to $nh/\log^{2} n\rightarrow\infty$ and $h^{2}\log n\rightarrow0$. Assumption (ref)(ii) can be seen as a uniform version of the so-called restricted eigenvalue condition widely used in high-dimensional linear regression models BRT09,BM15. Appendix D in the supplement provides sufficient conditions for the high-dimensional locally stationary Gaussian time series to satisfy Assumption (ref)(ii). Furthermore, with the Hanson-Wright inequality for time-varying (non-Gaussian) VAR processes ZW21, we may show that $\max_{1\leq t\leq n}\left\Vert {\boldsymbol{\Psi}}(\tau_{t})-\mathsf{E}[{\boldsymbol{\Psi}} (\tau_{t})]\right\Vert _{\max}=O_{P}\left( \sqrt{\log(n\vee d)/(nh)}\right) $. Then, using Lemma D.1 in Appendix D and assuming $s\sqrt{\log(n\vee d)/(nh)}=o(1)$, a sufficient condition for ((ref)) is \[ \min_{1\leq i\leq d}\min_{1\leq t\leq n}\inf_{u\in\mathcal{B}_{i}(\tau_{t} )}u^{^{\intercal}} \mathsf{E}\left[ {\boldsymbol{\Psi}}(\tau_{t})\right] u\geq\kappa_{0}. \]

\setcounter{theorem}{0}

Theorem (ref) shows that the preliminary time-varying LASSO estimates of the transition matrices are uniformly consistent with the convergence rates relying on $s$ and $\lambda_{1}$. Although the dimension of variates $d$ is allowed to diverge at an exponential rate of $n$, the number of significant elements in ${\boldsymbol{\alpha}}_{i\bullet}(\cdot)$ cannot diverge too fast in order to guarantee the consistency property. Furthermore, the uniform convergence result ((ref)) can be strengthened to

A similar uniform convergence property holds for the first-order derivative function estimates, see (A.1) in the proof of Theorem (ref).

The oracle property of the weighted group LASSO estimates

Denote the complement of $\mathcal{J}_{i}$ and $\mathcal{J}_{i}^{\prime}$ as $\overline{\mathcal{J}}_{i}$ and $\overline{\mathcal{J}}_{i}^{\prime}$, respectively, i.e., $\overline{\mathcal{J}}_{i}=\bigcap_{t=1}^{n}\left\{ j:\ \alpha_{i,j}(\tau_{t})=0\right\} $ and $\overline{\mathcal{J}} _{i}^{\prime}=\bigcap_{t=1}^{n}\left\{ j:\ \alpha_{i,j}^{\prime}(\tau _{t})=0\right\} $. Let ${\mathbf{A}}^{o}=\left( {\boldsymbol{\alpha} }_{\bullet1}^{o},\mathcal{\ldots},{\boldsymbol{\alpha}}_{\bullet n} ^{o}\right) ^{^{\intercal}}$ and ${\mathbf{B}}^{o}=\left( {\boldsymbol{\beta }}_{\bullet1}^{o},\mathcal{\ldots},{\boldsymbol{\beta}}_{\bullet n} ^{o}\right) ^{^{\intercal}}$, where ${\boldsymbol{\alpha}}_{\bullet t} ^{o}=(\alpha_{1|t}^{o},\mathcal{\ldots},\alpha_{d|t}^{o})^{^{\intercal}}$ with $\alpha_{j|t}^{o}=0$ for $j\in\overline{\mathcal{J}}_{i}$ and ${\boldsymbol{\beta}}_{\bullet t}^{o}=(\beta_{1|t}^{o},\mathcal{\ldots} ,\beta_{d|t}^{o})^{^{\intercal}}$ with $\beta_{j|t}^{o}=0$ for $j\in \overline{\mathcal{J}}_{i}^{\prime}$. Define the (infeasible) oracle estimates:

as the values of ${\mathbf{A}}^{o}$ and ${\mathbf{B}}^{o}$ that minimise $\mathcal{Q}_{i}({\mathbf{A}}^{o},{\mathbf{B}}^{o})$. We need to impose the following condition on the tuning parameter $\lambda_{2}$ and the lower bounds for the significant time-varying coefficients in the transition matrix.

When $s$ is a fixed positive integer, $h\propto n^{-1/5}$, $\lambda_{1}\propto n^{-2/5+\eta_{0}}$ with $0<\eta_{0}<1/5$, and $d\sim\exp\left\{ n^{\eta_{1} }\right\} $ with $0<\eta_{1}<\eta_{0}$, it is easy to verify Assumption (ref)(i) by setting $\lambda_{2}\propto n^{1/2-\eta_{2}}$ with $0<\eta_{2}<2/5-[\eta_{0}\vee(2\eta_{1})]$. Assumption (ref)(ii) imposes restrictions on the lower bounds for the time-varying coefficient functions and their deviations from the means. These restrictions are weaker than Assumption 6(ii) in LKZ15 and Assumption 8 in CLWZ21, and they ensure that the significant coefficient functions and derivatives can be detected w.p.a.1.

Since the penalised local linear estimates are identical to the infeasible oracle estimates defined in ((ref)) and ((ref)) w.p.a.1, the sparsity property holds for the global model selection procedures proposed in Section (ref), i.e., the zero elements in the time-varying transition matrix can be estimated exactly as zeros. Following the proof of Theorem (ref), we may verify properties (i)--(iv) for the folded concave penalty function discussed in FXZ14 w.p.a.1. Hence, Theorem (ref) may be regarded as a generalisation of Theorem 1 in FXZ14 and Theorem 3.1 in LKZ15 to high-dimensional locally stationary time series.

With the oracle property in Theorem (ref), it is straightforward to derive the following consistency property of the network estimates for the directed edges of Granger causality linkages.

\setcounter{corollary}{0}

Uniform consistency of the time-varying CLIME estimates

To derive the uniform consistency property of the time-varying CLIME estimates, we need the following conditions on the tuning parameters $b$ and $\lambda_{3}$.

The following theorem gives the uniform convergence rates of the time-varying precision matrix estimate $\widehat{\boldsymbol{\Omega}}(\tau)$ under various matrix norms.

The uniform convergence rates in Theorem (ref) rely on $\nu _{n,d}^{\diamond}$ and $\nu_{n,d}^{\ast}$. The first rate $\nu_{n,d} ^{\diamond}$ is the conventional uniform convergence rate for nonparametric kernel-based quantities, whereas the second rate $\nu_{n,d}^{\ast}$ is from the approximation errors of $\widehat{e}_{t}$ to the latent VAR errors $e_{t} $. Note that the dimension $d$ affects the uniform convergence rates via $\xi_{d}$ and $\log(n\vee d)$, and the uniform consistency property holds in the ultra-high dimensional setting when $d$ diverges at an exponential rate of $n$. Theorem (ref) can be seen as an extension of Theorem 1 in CLL11 to the high-dimensional locally stationary time series setting.

From Theorem (ref), we readily have the following consistency property for the network estimates of the undirected edges of partial correlation linkages.

Factor-adjusted time-varying VAR and networks

\setcounter{equation}{0}

In this section, we let $(Z_{t}:t=1,\mathcal{\ldots},n)$ with $Z_{t} =(z_{t,1},\mathcal{\ldots},z_{t,d})^{^{\intercal}}$ be an observed sequence of $d$-dimensional random vectors. To accommodate strong cross-sectional dependence which is not uncommon for large-scale time series collected in practice, we assume that $Z_{t}$ is generated by an approximate factor model:

where ${\boldsymbol{\Lambda}}=(\Lambda_{1},\mathcal{\ldots},\Lambda _{d})^{^{\intercal}}$ is a $d\times k$ matrix of factor loadings, $F_{t}$ is a $k$-dimensional vector of latent factors and $(X_{t})$ is assumed to satisfy the time-varying VAR model ((ref)). More generally, we may assume the following time-varying factor model structure:

where ${\boldsymbol{\Lambda}}_{t}={\boldsymbol{\Lambda}}(t/n)$ is a time-varying factor loading matrix with each entry being a smooth function of scaled time. The approximate factor model and its time-varying generalisation have been extensively studied in the literature CR83, BN02, SW02, MHvS11, SW17. The primary interest of this section is to estimate the time-varying networks for the idiosyncratic error vector $X_{t}$. Even though the components of $Z_{t}$ may be highly correlated, those of $X_{t}$ are often only weakly correlated. Hence, it is sensible to impose the sparsity assumption on the time-varying transition and precision matrices of the idiosyncratic error process, making it possible to apply the estimation methodology proposed in Section (ref). However, this is non-trivial as neither the common components (${\boldsymbol{\Lambda}}F_{t}$ or ${\boldsymbol{\Lambda}}_{t}F_{t}$) nor the idiosyncratic error components are observable. Motivated by recent work on bridging factor and sparse models for high-dimensional data FMM21, KM22, we next use the principal component analysis (PCA) or its localised version to remove the common components driven by latent factors in the observed time series data.

Let ${\mathbf{Z}}=\left( Z_{1},\mathcal{\ldots},Z_{n}\right) ^{^{\intercal} }$, ${\mathbf{F}}=\left( F_{1},\mathcal{\ldots},F_{n}\right) ^{^{\intercal} }$ and ${\mathbf{X}}=\left( X_{1},\mathcal{\ldots},X_{n}\right) ^{^{\intercal}}$. For the conventional factor model ((ref)), we conduct an eigenanalysis on the $n\times n$ matrix ${\mathbf{Z}}{\mathbf{Z} }^{^{\intercal}}$. The estimate of ${\mathbf{F}}$, denoted as $\widehat {\mathbf{F}}=\left( \widehat{F}_{1},\mathcal{\ldots},\widehat{F}_{n}\right) ^{^{\intercal}}$, is obtained as the $n\times k$ matrix consisting of the eigenvectors (multiplied by $\sqrt{n}$) corresponding to the $k$ largest eigenvalues of ${\mathbf{Z}}{\mathbf{Z}}^{^{\intercal}}$. The factor loading matrix is estimated by $\widehat{\boldsymbol{\Lambda}}=\left( \widehat {\Lambda}_{1},\mathcal{\ldots},\widehat{\Lambda}_{d}\right) ^{^{\intercal} }={\mathbf{Z}}^{^{\intercal}}\widehat{\mathbf{F}}/n$. Consequently, the common component ${\boldsymbol{\Lambda}}F_{t}$is estimated by $\widehat {\boldsymbol{\Lambda}}\widehat{F}_{t}$ and the idiosyncratic error component $X_{t}$ is estimated by

For the time-varying factor model ((ref)), the above PCA estimation procedure needs some amendments. Specifically, let \[ K_{t,h_{\ast}}(\tau)=\frac{K_{h_{\ast}}(\tau_{t}-\tau)}{\sum_{s=1} ^{n}K_{h_{\ast}}(\tau_{s}-\tau)},\ \ 0<\tau<1, \] where $h_{\ast}$ is a bandwidth and $K_{h_{\ast}}(\cdot)$ is defined as in Section (ref), and define the localised data matrix: \[ {\mathbf{Z}}(\tau)=\left[ Z_{1}(\tau),\mathcal{\ldots},Z_{n}(\tau)\right] ^{^{\intercal}}\ \ \mathrm{with}\ \ Z_{t}(\tau)=Z_{t}K_{t,h_{\ast}}^{1/2} (\tau). \] Through an eigenanalysis on the matrix ${\mathbf{Z}}(\tau){\mathbf{Z} }^{^{\intercal}}(\tau)$, we can obtain the local PCA estimates of the factors and factor-loading matrix, denoted by $\widehat{\mathbf{F}}(\tau)=\left[ \widehat{F}_{1}(\tau),\mathcal{\ldots},\widehat{F}_{n}(\tau)\right] ^{^{\intercal}}$ and $\widehat{\boldsymbol{\Lambda}}(\tau)$, respectively. Then, the idiosyncratic error vector $X_{t}$ is approximated by

where we've kept the same notation $\widehat{X}_{t}$ as in ((ref)) to avoid notational burden.

As in Section (ref), we only consider the time-varying VAR(1) model for the idiosyncratic error vector. With the approximation $\widehat{X}_{t}$, we can apply the three-stage estimation procedure proposed in Section (ref). Denote the preliminary time-varying LASSO estimate as $\widetilde{\alpha} _{ij}^{\dagger}(\cdot)$, the second-stage weighted group LASSO estimate as $\widehat{\alpha}_{ij}^{\dagger}(\cdot)$, and the factor-adjusted time-varying precision matrix estimate as $\widehat{\boldsymbol{\Omega}}^{\dagger} (\cdot)=\left[ \widehat\omega_{ij}^{\dagger}(\cdot)\right] _{d\times d}$. Subsequently, we may construct the uniform network estimates $\widehat {\mathbb{E}}_{n}^{G,{\dagger}}$ and $\widehat{\mathbb{E}}_{n}^{P,{\dagger}}$, defined similarly to $\widehat{\mathbb{E}}_{n}^{G}$ and $\widehat{\mathbb{E} }_{n}^{P}$ in ((ref)) and ((ref)), but with $\widehat{\alpha }_{ij}(\cdot)$ and $\widehat\omega_{ij}(\cdot)$ replaced by $\widehat{\alpha }_{ij}^{\dagger}(\cdot)$ and $\widehat\omega_{ij}^{\dagger}(\cdot)$, respectively. To derive the convergence properties of these factor-adjusted estimates, we need the following assumption, which modifies Assumptions (ref)--(ref) to incorporate the approximation error of the idiosyncratic error components.

Assumption (ref)(i) imposes a high-level condition on the approximation of the latent $X_{t}$ in the factor model, i.e., the approximation error $\delta_{X}$ uniformly converges to zero with a rate faster than $s^{-1} [\log(n\vee d)]^{-1/2}$. By Corollary 1 in FLM13, a typical rate for the approximation error from PCA estimation of the conventional factor model ((ref)) is

where $\upsilon>2$ is a positive number related to moment restrictions. From Theorem 3.5 in SW17, we may obtain the typical uniform rate for $\delta_{X}$ under the time-varying factor model ((ref)) when the local PCA estimation is used. In Assumption (ref)(ii)--(iv), we amend Assumptions (ref)(i), (ref)(i) and (ref)(ii) to incorporate the approximation error $\delta_{X}$. However, if we further assume that $h\propto n^{-1/5}$ and $d$ diverges at a polynomial rate of $n$ satisfying $d\gg n^{1+2/\upsilon}$, then the rate in ((ref)) can be simplified to $\delta_{X}=O_{P}\left( (\log d)n^{-1/2}\right) =o_{P}(h^{2})$ and thus $\zeta_{n,d}\propto\zeta_{n,d}^{\dagger}$. Consequently, we may remove Assumption (ref)(ii)--(iv) and $\delta_{X}$ would not be involved in the estimation convergence rates under model ((ref)).

The following two propositions extend the theoretical results in Section (ref) to the factor-adjusted time-varying VAR and networks.

\setcounter{prop}{0}

Monte-Carlo simulation

\setcounter{equation}{0}

In this section, we provide four simulated examples to examine the finite-sample numerical performance of the proposed high-dimensional time-varying VAR and network estimates. Throughout this section, we denote the proposed time-varying weighted group LASSO method as tv-wgLASSO and the time-varying CLIME method as tv-CLIME. We compare the performance of the tv-wgLASSO with the (infeasible) time-varying oracle estimation, denoted as tv-Oracle, which estimates only the true significant coefficient functions (assuming they were known), and the unpenalised full time-varying estimation, denoted as tv-Full, which estimates all the coefficient functions without penalisation. We compare the performance of tv-CLIME with the time-varying graphical LASSO estimation, denoted as tv-GLASSO, which is implemented using the R package “glassoFast" on the VAR residuals. In addition, to investigate the loss of estimation accuracy due to the VAR model error approximation, we also report results from the infeasible tv-CLIME, which directly uses the VAR errors (rather than residuals) in the estimation of the precision matrices.

In the simulation, we use the Epanechnikov kernel $K(t)=0.75(1-t^{2})_{+}$ with bandwidth $h=b=0.75[\log(d)/n]^{1/5}$ as in LKZ15. The bandwidth for the local PCA is set as $h_{\ast}=(2.35/\sqrt{12})[\sqrt{d}/n] ^{1/5}$ as in SW17. We set the sample size $n$ as 200 and 400, and the dimension $d$ as 50 and 100. Although such dimensions are smaller than the sample size, when $n=200$ and $d=100$, the “effective sample size" used in each local linear estimation in ((ref)) is approximately $2nh\approx140$, which is smaller than the combined number of unknown coefficient functions and their derivative, $2d=200$. Consequently, in this case we fail to implement the naive tv-Full estimation. There are three tuning parameters in the proposed estimation procedure: $\lambda_{1}$ in the first stage of preliminary time-varying LASSO estimation, $\lambda_{2}$ in the second stage of time-varying weighted group LASSO, and $\lambda_{3}$ in the third stage of time-varying CLIME. They are selected by the Bayesian information criterion (BIC), the generalised information criterion (GIC), and the extended Bayesian information criterion (EBIC), respectively. Appendix E in the supplement gives definitions of these information criteria.

To evaluate whether the time-varying model structure is accurately estimated, we report the false positive (FP), the false negative (FN), the true positive rate (TPR), the true negative rate (TNR), the positive predictive value (PPV), the negative predictive value (NPV), the F1 score (F1), and the Matthews correlation coefficient (MCC). Definitions of these measures are available in Appendix E of the supplement. To evaluate the performance of the coefficient estimators, we report the average R square (average $R^{2}$) over all the dimensions, the average scaled Frobenius norm of estimation errors of coefficient functions (EE$_{A}$), and the root-mean-squared error of the errors (RMSE$_{e}$). Taking our proposed tv-wgLASSO estimator for time-varying VAR(1) as an example, \[ \mathrm{EE}_{A}=\frac{1}{n\sqrt{d}}\sum_{t=1}^{n}\left\Vert \widehat {\mathbf{A}}_{1}(\tau_{t})-{\mathbf{A}}_{1}(\tau_{t})\right\Vert _{F}\ \ \mathrm{and}\ \ \mathrm{RMSE}_{e}=\sqrt{\frac{1}{nd}\sum_{i=1}^{d} \sum_{t=1}^{n}(\widehat{e}_{t,i}-{e}_{t,i})^{2}}. \] To evaluate the performance of the precision matrix estimators, we report the average scaled Frobenius norm of estimation error ($\mathrm{EE}_{\Omega}$) defined as \[ \mathrm{EE}_{\Omega}=\frac{1}{n\sqrt{d}}\sum_{t=1}^{n}\left\Vert \widehat{\boldsymbol{\Omega}}(\tau_{t})-{\boldsymbol{\Omega}}(\tau _{t})\right\Vert _{F}. \] All the above measures are calculated for each Monte Carlo replication and then averaged over $100$ replications.

\noindentExample 1.\ \ The data is generated from a time-varying VAR(1) model with ${\mathbf{A}}_{1}(\tau)$ being a diagonal matrix for all $\tau\in[0,1]$. Each diagonal entry of ${\mathbf{A}}_{1}(\tau)$ independently takes a value of either $0.64\Phi(5(\tau-1/2))$ or $0.64-0.64\Phi (5(\tau-1/2))$ with an equal probability of 0.5, where $\Phi(\cdot)$ is the standard normal distribution function. We set ${\boldsymbol{\Omega}}(\tau)$ to be a block diagonal matrix: ${\boldsymbol{\Omega}}(\tau)={\mathbf{I}} _{d/2}\otimes{\boldsymbol{\Omega}}_{\ast}(\tau)$, where ${\boldsymbol{\Omega} }_{\ast}(\tau)=\left[ \omega_{ij,\ast}(\tau)\right] _{2\times2}$ with $\omega_{11,\ast}(\tau)=\omega_{22,\ast}(\tau)\equiv1$, and $\omega_{12,\ast }(\tau)=\omega_{21,\ast}(\tau)=1.4\Phi(5(\tau-1/2))-0.7$. The diagonal structure of ${\mathbf{A}}_{1}(\tau)$ implies that no Granger causality exists between variables, whereas the block diagonal structure of ${\boldsymbol{\Omega}}(\tau)$ results in weak cross-sectional dependence between the components of $X_{t}$.

Table (ref) reports the estimation results of the time-varying transition matrices and Granger networks. For the proposed tv-wgLASSO, the FP and FN values are very small compared with $d^{2}$ (the total number of potential directed Granger causality linkages or entries of the transition matrix). This leads to large values of the TPR, TNR, PPV, NPV, F1 and MCC measures, all of which are close to $1$. We can also see that the FP and FN values double when $d$ increases from $50$ to $100$, but decrease substantially when $n$ grows from $200$ to $400$. These results clearly show that tv-wgLASSO can accurately recover the time-varying Granger network as long as the sample size is moderately large. The average $R^{2}$ of tv-wgLASSO is close to that of tv-Oracle, but the naive tv-Full method tends to have large $R^{2}$ due to model over-fitting. Although the EE$_{A}$ values of tv-wgLASSO are larger than those of tv-Oracle when $n=200$, they drop significantly and are even slightly smaller than those of tv-Oracle when $n=400$. A similar pattern can be observed in RMSE$_{e}$, indicating that the proposed tv-wgLASSO is capable of providing good approximations to VAR errors, which are used in the subsequent time-varying precision matrix estimation. Unsurprisingly, the tv-Full method fails to estimate the time-varying transition matrix when $d=100$ and $n=200$.

Table (ref) reports the estimation results of the time-varying precision matrices and partial correlation networks. When $n=200$, both tv-CLIME and tv-GLASSO have zero FP values, whereas tv-CLIME has smaller FN than tv-GLASSO. Hence, the proposed tv-CLIME performs better than tv-GLASSO in terms of the F1 and MCC measures. When $n=400$, both tv-CLIME and tv-GLASSO correctly recover the time-varying partial correlation networks. In terms of the precision matrix estimation accuracy (EE$_{\Omega}$), tv-GLASSO performs slightly better than tv-CLIME. In addition, by comparing the tv-CLIME and the infeasible tv-CLIME, we may conclude that the VAR error approximation has negligible impact on the precision matrix and partial correlation network estimation.

\noindentExample 2.\ \ The data is generated from a time-varying VAR(1) model with ${\mathbf{A}}_{1}(\tau)$ being an upper triangular matrix for all $\tau\in[0,1]$. Each diagonal entry of ${\mathbf{A}}_{1}(\tau)$ takes the value of $0.7\Phi(5(\tau-1/2))$, each super-diagonal entry takes the value of $0.7-0.7\Phi(5(\tau-1/2))$, and the remaining entries take the value of $0$. We set ${\boldsymbol{\Omega}}(\tau)=\left[ \omega_{ij}(\tau)\right] _{d\times d}$ to be a banded symmetric matrix for all $\tau\in[0,1]$ with $\omega_{ii}(\tau)\equiv1$, $\omega_{i,(i+1)}(\tau)=0.7\Phi(5(\tau-1/2))-0.7$, $\omega_{i,(i+2)}(\tau)=0.7-0.7\Phi(5(\tau-1/2))$, and $\omega_{i,j} (\tau)\equiv0$ if $|i-j|>2$.

Table (ref) reports the estimation results of the time-varying transition matrices and Granger networks. Note that the time series variables in this example are more correlated to each other than those in Example 1, which affects the network estimation accuracy. When $d=100$ and $n=200$, the FP and FN values of tv-wgLASSO reach their maximum at 20.73 and 37.55, respectively, whereas the F1 and MCC values are around $0.85$. As in Example 1, the F1 and MCC values increase when $n$ increases from $200$ to $400$, and again the average $R^{2}$ of tv-wgLASSO is close to that of tv-Oracle. However, tv-wgLASSO has much larger EE$_{A}$ and RMSE$_{e}$ than tv-Oracle.

Table (ref) reports the estimation results of the time-varying precision matrices and partial correlation networks. It follows from the EE$_{A}$ and RMSE$_{e}$ results in Table (ref) that the VAR error approximation is poorer than that in Example 1. Consequently the proposed tv-CLIME performs worse than the infeasible tv-CLIME using the true VAR errors directly in the estimation. In particular, FN of the tv-CLIME is much larger than that of the infeasible tv-CLIME when $n=200$. Due to the same reason, the infeasible tv-CLIME also outperforms the tv-GLASSO. In addition, we find that the tv-CLIME is better than the tv-GLASSO in recovering the time-varying precision network when $n=200$, and they perform equally well when $n=400$.

\noindentExample 3.\ \ The data is generated from a VAR(1) model with ${\mathbf{A}}_{1}(\tau)=\left[ a_{ij}(\tau)\right] _{d\times d}$ being a Toeplitz matrix and $a_{ij}(\tau)=(0.4-0.1\tau)^{|i-j|+1}$. We also set ${\boldsymbol{\Omega}}(\tau)=\left[ \omega_{ij}(\tau)\right] _{d\times d}$ to be a Toeplitz matrix with $\omega_{ij}(\tau)= (0.8-0.1\tau)^{|i-j|}$. In this example, both the transition and precision matrices are non-sparse, and we aim to examine how our proposed methods perform when the (exact) sparsity assumption fails.

Table (ref) reports the estimation errors of the various methods considered. In this example, the tv-Oracle is equivalent to tv-Full and both suffer from the curse of dimensionality in the conventional local linear estimation procedure for the time-varying transition matrices (in particular when $d=100$ and $n=200$). Consequently, the EE$_{A}$ and RMSE$_{e}$ of the tv-wgLASSO are much smaller than those of the tv-Oracle. The EE$_{\Omega}$ results of the tv-CLIME are very close to those of the infeasible tv-CLIME, suggesting that the VAR error approximation has little impact on the tv-CLIME performance as discussed in Example 1. In addition, the EE$_{\Omega}$ results of the tv-CLIME and Oracle tv-CLIME are generally close to those of tv-GLASSO. The simulation results show that the proposed tv-wgLASSO and tv-CLIME perform reasonably well when the sparsity assumption on transition and precision matrices is not satisfied.

\noindentExample 4.\ \ The data is generated from a factor-adjusted time-varying VAR model in the form of ((ref)). The idiosyncratic errors of the time-varying factor model are generated from a VAR(1) model in Example 2. The two factors in $F_{t}=(F_{t,1},F_{t,2})^{^{\intercal}}$ are generated from two univariate AR(1) processes: $F_{t,1}=0.6F_{t-1,1}+\sqrt{1-0.6^{2} }u_{t,1}^{F}$ and $F_{t,2}=0.3F_{t-1,2}+\sqrt{1-0.3^{2}}u_{t,2}^{F}$, where $u_{t,1}^{F}$ and $u_{t,2}^{F}$ are independently drawn from a standard normal distribution. The factor-loading matrix is defined as ${\boldsymbol{\Lambda} }_{t}=\left( \Lambda_{t,1},\Lambda_{t,2}\right) $ where $\Lambda_{t,1} \equiv\Lambda_{1}$ is a time-invariant vector drawn from a $d$-dimensional standard multivariate normal distribution and $\Lambda_{t,2}=(\Lambda _{1t,2},\mathcal{\ldots},\Lambda_{dt,2})^{^{\intercal}}$ with $\Lambda _{it,2}=2/\left( 1+\exp\{-2[10(t/n)-5(i/d)-2]\}\right) $ for $i=1,\mathcal{\ldots},d$.

Table (ref) reports the estimation results of the time-varying transition matrices and Granger networks for the idiosyncratic errors, and Table (ref) reports the estimation results of the time-varying precision matrices and partial correlation networks. Comparing with the results in Tables (ref) and (ref), we can observe that the factor-adjusted estimation introduces additional estimation errors, leading to smaller values of F1 and MCC. The impact is more marked when $n=200$ but reduces substantially when $n=400$. As in the previous examples, the F1 and MCC values increase when $n$ increases from $200$ to $400$. Thus we may conclude that, although the factor model estimation errors are passed onto the three-stage estimation procedure, their impact on the estimation of the networks is not significant when the sample size is moderately large ($n=400$).

An empirical application

\setcounter{equation}{0}

In this section, we apply the proposed methods to estimate the Granger causality and partial correlation networks using the FRED-MD macroeconomic dataset. The dataset, available on the Fred-MD website\footnote{https://research.stlouisfed.org/econ/mccracken/fred-databases/} , consists of $127$ U.S. macroeconomic variables observed monthly over the period from January 1959 to July 2022. These macroeconomic variables can be classified into eight groups: consumption, orders and inventories; housing; interest and exchange rates; labour market; money and credit; output and income; prices; and the stock market. More detailed description can be found in MN16.

We follow MN16 and MN20 to remove outliers and fill missing values. Each variable is standardised to have zero mean and unit variance. We consider the two factor modelling methods in Section (ref) to accommodate strong cross-sectional dependence: the approximate factor model ((ref)) with constant factor loadings, and the time-varying factor model ((ref)) with dynamic factor loadings. The information criteria proposed by BN02 and SW17 are used to determine the number of factors in these two models (see Appendix E in the supplement for description of the criteria). Seven factors are selected for the factor model with constant loadings, whereas only four are selected for the time-varying factor model. Since the latter provides a more parsimonious model specification, we hereafter report network estimation results only for this model. The estimated idiosyncratic errors, denoted as $\widehat{x}_{t,i}$, $i=1,\mathcal{\ldots},127$, $t=1,\mathcal{\ldots},763$, are then used for our empirical analysis. MPS22 suggest determining the optimal order of a high-dimensional VAR model via a ratio criterion, comparing the Frobenius norms of the estimated transition matrices over different lags. We extend their criterion to the time-varying VAR model context (see Appendix E in the supplement for detail) and subsequently select the time-varying VAR(1) model for $\widehat{X} _{t}=\left( \widehat{x}_{t,1},\mathcal{\ldots},\widehat{x}_{t,127}\right) ^{^{\intercal}}$.

Figure (ref) plots the estimated Granger networks from the static VAR(1) and the time-varying VAR(1) models. From the estimated time-varying transition matrix, we uncover $190$ directed linkages in the Granger causality network, among which $78$ are self-linkages and $143$ are linkages within the same category. In particular, the self-linkages, which correspond to the significant diagonal entries of the transition matrix, indicate that the macroeconomic variables in the following four categories: consumption, orders and inventories; interest and exchange rates; money and credit; and prices, are more persistent than the others, even though all the variables have been transformed into stationary ones in the preliminary analysis. By contrast, we find 155 directed linkages for the Granger network estimated via static VAR(1) and hence, our time-varying VAR(1) model captures more linkages in the network estimation. Figure (ref) plots the Granger networks estimated without factor adjustment. Compared with the factor-adjusted version, the Granger network via time-varying VAR(1) is more dense with $1118$ directed linkages, among which $104$ are self-linkages and $432$ are within categories. As pointed out by MN16, common factors, which may be interpreted as business cycles, are the main sources of the Granger causalities between macroeconomic variables, leading to a rather dense network structure. On the other hand, the estimated Granger network via static VAR(1) without factor adjustment has only 450 linkages.

We further explore the dynamic smooth structural changes of Gaussian causality linkages. Taking the logarithmic growth rate of S&P PE ratio (S&P PE ratio)\footnote{We show in the parentheses the variable names used in the FRED-MD dataset. The variable transformation is conducted following the guideline in the dataset.} as an example, there are four directed linkages to this variable: acceleration of the logarithmic monetary base (BOGMBASE), the logarithmic return of S&P 500 index (S&P 500), the logarithmic return of S&P 500 industrials index (S&P: indust), and the logarithmic growth rate of the S&P PE ratio which is a self-linkage. We re-estimate the corresponding time-varying coefficients using the nonparametric autoregression model with only the four selected predictors, and draw the 90% confidence bands using the R package “tvReg". Figure (ref) plots the estimated curves of the four coefficient functions. We find that the logarithmic growth rate of S&P PE ratio is generally persistent and positively correlated to BOGMBASE in the most recent two decades. The estimated time-varying coefficient of the S&P 500 industrials index return is significant but close to zero. It is thus unsurprising that the static VAR(1) model with classic LASSO penalty does not detect the Granger causality linkage from this variable. In fact, LASSO tends to select only one variable in a group of highly-correlated predictors. Due to high correlation between the two index returns, only the S&P 500 Index return is selected in the static VAR(1) model. In contrast, the proposed time-varying LASSO selects both of the two index returns at different time periods, and the second-stage weighted group LASSO aggregates the information over time and selects both index returns.

We plot the estimated partial correlation networks in Figure (ref), which are generally sparse. Using the factor-adjusted time-varying CLIME, $234$ undirected linkages are detected in the estimated network, among which $205$ linkages are within the same category. In contrast, the estimated network without factor adjustment contains $236$ linkages with $211$ in the same category. Unlike the Granger network estimation, it seems that whether to make factor adjustment or not has little impact on the partial correlation network estimation.

We next examine the time-varying pattern of partial correlation linkages between S&P PE ratio and four other variables: S&P 500, S&P: indust, S&P div yield (the increment of S&P composite common stock: dividend yield), and BAAFFM (the spread between Moody's seasoned baa corporate bond and effective federal funds rate). We re-estimate the relevant time-varying functions with a 200-month moving window JV15, and draw the 90% confidence bands using R package “SILGGM" in Figure (ref). Note that the partial correlation has a sign opposite to the corresponding entry in the precision matrix. We find that S&P PE ratio is positively (partially) correlated with S&P 500 and S&P: indust, whilst negatively (partially) correlated with S&P div yield. The confidence bands in Figure (ref) suggest that time-invariant partial correlation linkages are inappropriate to describe the network structure of the FRED-MD data.

Conclusion

In this paper we estimate a general time-varying VAR model for high-dimensional locally stationary time series. A three-stage estimation procedure combining time-varying LASSO, weighted group LASSO and time-varying CLIME is developed to estimate both transition and error precision matrices, allowing smooth structural changes over time. The estimated transition and precision matrices are further used to construct dual network structures with directed Granger causality linkages and undirected partial correlation linkages, respectively. Under the sparse structural assumption and other technical conditions, we derive the uniform consistency and oracle properties for the developed estimates. In order to accommodate high correlation among large-scale time series and avoid directly imposing the sparsity assumption, we also extend the methodology and theory to a more general factor-adjusted time-varying VAR and network structures. Both the simulation and empirical studies show that the developed network model and methodology have reliable numerical performance in finite samples.

Supplementary materials

{ The supplement contains proofs of the main asymptotic theorems, some technical lemmas with proofs, verification of Assumption (ref)(ii) and discussions on tuning parameter selection.}